Merge pull request #63 from pielube/dev_calamatti

Operational period for electrolyzer and fuel cell / New Italian incentive system for REC
pielube · May 2, 2024 · 7d50adc · 7d50adc
2 parents 99f55d7 + acd61d4
commit 7d50adc
Show file tree

Hide file tree

Showing 18 changed files with 36,250 additions and 56 deletions.
diff --git a/core/economics.py b/core/economics.py
@@ -98,6 +98,8 @@ def NPV(file_studycase,
 
     results = {}                        # dictionary initialise economic results of each locations
 
+    peakP = sum(location_name.get('PV', {}).get('peakP',0) for location_name in studycase.values()) + sum(location_name.get('wind', {}).get('Npower',0)for location_name in studycase.values())/1000  #calculate total nominal power of REC
+
     for location_name in tc:            # for reach location
 
         results[location_name] = {}     # dictionary initialise economic results
@@ -254,14 +256,43 @@ def NPV(file_studycase,
 
         # REC incentives redistribution
         if 'REC' in economic_data:
-            csc     = balances[location_name]['electricity']['collective self consumption']
-
-            inc_pro = - csc * economic_data['REC']['incentives redistribution']['producers']/100 * economic_data['REC']['collective self consumption incentives']
+            csc = balances[location_name]['electricity']['collective self consumption']
+            if type(economic_data['REC']['collective self consumption incentives']) == str:     # if the pun price series is given (Italian legislation)
+                pun = pd.read_csv(path+'/energy_price/'+economic_data['REC']['collective self consumption incentives'])['0'].to_numpy()
+                if len(pun) < (c.HOURS_YEAR+2):                   # it means that the serie must be repeated for the simulation_length selected
+                    pun = np.tile(pun,int(c.simulation_years))
+                # check su potenza nominale impianto
+                if peakP > 600:
+                    csc_inc = [min(0.100, (0.060 + max(0,0.180-i))) for i in pun]
+                elif peakP > 200 and peakP <= 600:
+                    csc_inc = [min(0.110, (0.070 + max(0,0.180-i))) for i in pun]
+                else:
+                    csc_inc = [min(0.120, (0.080 + max(0,0.180-i))) for i in pun]
+                for i in range(len(csc_inc)):
+                    csc_inc[i] += 0.008    #ARERA
+
+                #PV incentives adjustment for the central and northern regions
+                if 'region' in economic_data['REC']:
+                    region = economic_data['REC']['region']
+                    centro = ['Lazio', 'Marche', 'Toscana', 'Umbria', 'Abruzzo']
+                    nord = ['Emilia-Romagna', 'Friuli-Venezia Giulia', 'Liguria', 'Lombardia', 'Piemonte',
+                            'Trentino-Alto Adige/Südtirol', "Valle d'Aosta/Vallée d'Aoste", "Veneto"]
+                    if region in centro:
+                        for i in range(len(csc_inc)):
+                            csc_inc[i] += 0.004
+                    elif region in nord:
+                        for i in range(len(csc_inc)):
+                            csc_inc[i] += 0.010
+
+            else:
+                csc_inc = economic_data['REC']['collective self consumption incentives']
+
+            inc_pro = - csc * economic_data['REC']['incentives redistribution']['producers']/100 * csc_inc
             inc_pro = np.tile(inc_pro,years_factor)
             inc_pro = np.reshape(inc_pro,(-1,8760))    
             results[location_name]['CF_studycase']['CSC'] += inc_pro.sum(axis=1,where=inc_pro>0) 
 
-            inc_con = csc * economic_data['REC']['incentives redistribution']['consumers']/100 * economic_data['REC']['collective self consumption incentives']
+            inc_con = csc * economic_data['REC']['incentives redistribution']['consumers']/100 * csc_inc
             inc_con= np.tile(inc_con,years_factor)
             inc_con = np.reshape(inc_con,(-1,8760))
             results[location_name]['CF_studycase']['CSC'] += inc_con.sum(axis=1,where=inc_con>0)   

diff --git a/core/rec.py b/core/rec.py
@@ -72,6 +72,7 @@ def __init__(self,structure,general,file_structure,file_general,path):
         else:
             check_pv = False
         self.weather = self.weather_generation(general,path,check,file_general) # check if metereological data have to been downloaded from PVgis or has already been done in a previous simulation
+        self.weather = pd.concat([self.weather] * c.simulation_years, ignore_index = True)
         if check == False:
             with open(f"previous_simulation/{file_general}.pkl", 'wb') as f: pickle.dump(general, f)
         if check_pv == False:

diff --git a/input_test_1/energy_market.json b/input_test_1/energy_market.json
@@ -1,5 +1,5 @@
 {								 					 
-"REC": {"collective self consumption incentives": 0.118,
+"REC": {"collective self consumption incentives": "pun_2023.csv",
 		"incentives redistribution": {	"producers": 20,
 										"consumers": 80,
 										"manager": 0},

diff --git a/input_test_1/energy_price/PUN_curve.py b/input_test_1/energy_price/PUN_curve.py
@@ -0,0 +1,268 @@
+# -*- coding: utf-8 -*-
+"""
+Create a daily PUN curve from 2023 daily data
+
+@author: Matti
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+
+#%% Jan ############
+jan = np.zeros((24,7))
+
+jan[:,0] = (195,190,180,180,185,195,200,200,205,200,200,200,200,200,200,200,200,260,240,240,205,200,195,190)
+jan[:,1] = (175,150,120,115,115,150,195,195,200,200,195,175,175,175,175,190,195,200,200,200,195,190,175,150)
+jan[:,2] = (120,120,100,120,125,140,170,190,220,190,175,170,170,170,175,195,200,205,210,205,200,190,175,140)
+jan[:,3] = (170,120,105,100,105,150,185,200,220,250,190,180,170,170,180,200,220,250,240,245,200,190,180,165)
+jan[:,4] = (175,175,170,155,155,160,175,175,175,175,165,150,150,135,150,175,175,205,220,250,220,180,175,170)
+jan[:,5] = (185,185,185,185,185,185,190,190,195,200,190,185,165,165,185,185,190,250,260,270,220,190,185,185)
+jan[:,6] = (200,200,200,200,200,200,200,200,200,200,185,165,165,165,185,200,200,205,215,240,240,205,200,200)
+
+jan_mean = np.mean(jan, axis = 1)
+jan_tot = np.sum(jan_mean)
+index_jan = np.array([hour / jan_tot for hour in jan_mean])
+weight_jan = np.array([index / index_jan.mean() for index in index_jan])
+
+#%% Feb ############
+feb = np.zeros((24,7))
+feb[:,0] = (145,145,145,145,145,145,160,200,230,195,170,150,150,150,150,160,170,190,220,205,190,170,160,150)
+feb[:,1] = (150,145,140,130,135,140,165,200,205,180,160,150,135,135,140,150,165,180,205,205,170,160,155,145)
+feb[:,2] = (145,140,135,130,135,150,170,185,215,165,155,150,150,150,150,155,165,200,225,215,185,165,160,155)
+feb[:,3] = (150,150,145,145,145,150,165,195,200,165,160,155,150,150,150,160,165,170,200,190,170,165,155,150)
+feb[:,4] = (150,140,135,135,135,145,165,200,200,170,165,160,150,150,155,160,165,195,220,215,185,165,160,150)
+feb[:,5] = (150,150,145,140,140,140,150,150,150,145,140,115,100,85,90,110,145,150,170,185,155,150,145,135)
+feb[:,6] = (140,135,130,125,125,115,135,140,120,115,110,100,95,85,100,110,145,150,190,200,190,180,160,150)
+
+feb_mean = np.mean(feb, axis = 1)
+feb_tot = np.sum(feb_mean)
+index_feb = np.array([hour / feb_tot for hour in feb_mean])
+weight_feb = np.array([index / index_feb.mean() for index in index_feb])
+
+#%% Mar ############
+mar = np.zeros((24,7))
+mar[:,0] = (130,110,110,135,115,140,150,150,150,150,145,140,115,100,120,140,150,155,220,235,185,150,150,140)
+mar[:,1] = (135,115,90,85,95,130,150,150,170,145,125,120,105,105,120,130,145,150,155,195,180,150,135,110)
+mar[:,2] = (130,115,115,110,110,140,150,165,180,155,135,110,105,100,105,110,140,155,205,210,175,150,145,140)
+mar[:,3] = (130,125,100,100,100,130,150,160,150,110,100,100,85,90,95,105,110,145,195,215,200,145,140,135)
+mar[:,4] = (125,115,110,100,110,125,130,135,140,115,95,90,80,80,105,130,130,155,195,205,185,145,135,130)
+mar[:,5] = (135,135,130,130,135,135,140,140,135,130,100,95,85,80,90,95,135,140,195,215,200,175,140,135)
+mar[:,6] = (135,130,125,120,125,125,125,110,100,100,100,100,100,95,95,100,125,130,145,180,175,145,135,125)
+
+mar_mean = np.mean(mar, axis = 1)
+mar_tot = np.sum(mar_mean)
+index_mar = np.array([hour / mar_tot for hour in mar_mean])
+weight_mar = np.array([index / index_mar.mean() for index in index_mar])
+
+#%% Apr ############
+
+apr = np.zeros((24,7))
+apr[:,0] = (130,125,125,120,125,130,145,180,200,160,140,135,130,125,130,130,135,140,145,195,205,180,145,140)
+apr[:,1] = (140,135,130,130,135,140,145,190,200,140,130,125,115,115,115,125,130,135,145,205,225,200,145,140)
+apr[:,2] = (135,130,125,125,125,135,145,200,205,145,140,135,115,110,115,125,130,135,145,195,205,205,145,140)
+apr[:,3] = (135,125,115,115,120,130,145,190,205,190,145,140,125,125,125,130,135,140,155,200,215,205,145,140)
+apr[:,4] = (135,130,125,120,120,130,145,200,210,190,140,135,110,105,110,115,125,135,145,190,200,200,150,145)
+apr[:,5] = (165,170,145,135,130,130,165,165,125,120,110,85,70,65,55,65,80,115,125,160,185,175,135,130)
+apr[:,6] = (125,125,120,120,120,120,125,130,125,120,125,120,120,100,95,100,100,120,130,155,195,155,130,125)
+
+apr_mean = np.mean(apr, axis = 1)
+apr_tot = np.sum(apr_mean)
+index_apr = np.array([hour / apr_tot for hour in apr_mean])
+weight_apr = np.array([index / index_apr.mean() for index in index_apr])
+
+#%% May ############
+
+may = np.zeros((24,7))
+may[:,0] = (110,105,100,100,100,110,125,160,150,130,115,110,100,100,95,100,100,105,115,115,115,110,105,100)
+may[:,1] = (90,75,75,70,75,90,110,125,145,125,115,110,100,95,100,100,105,105,110,120,130,120,110,100)
+may[:,2] = (100,100,95,95,90,95,110,120,150,140,120,110,95,90,100,100,110,120,140,145,150,135,115,115)
+may[:,3] = (110,105,105,105,105,110,115,120,130,115,110,100,95,95,95,105,110,110,120,140,140,120,110,110)
+may[:,4] = (115,100,100,95,95,100,115,125,145,145,125,115,105,100,105,110,115,110,120,140,145,130,110,105)
+may[:,5] = (100,100,95,90,85,90,100,95,105,110,105,100,100,95,100,100,105,105,110,120,125,110,110,105)
+may[:,6] = (105,100,95,90,90,95,90,100,105,100,95,90,80,65,70,75,60,95,100,110,130,125,115,110)
+
+may_mean = np.mean(may, axis = 1)
+may_tot = np.sum(may_mean)
+index_may = np.array([hour / may_tot for hour in may_mean])
+weight_may = np.array([index / index_may.mean() for index in index_may])
+
+#%% Jun ############
+
+jun = np.zeros((24,7))
+jun[:,0] = (95,90,85,95,95,100,105,115,125,105,95,90,85,85,90,90,90,95,110,125,130,115,105,100)
+jun[:,1] = (100,100,100,100,100,100,105,125,125,125,105,105,100,95,100,105,105,105,105,125,125,115,105,100)
+jun[:,2] = (100,90,85,85,90,95,105,120,125,120,115,110,105,100,105,105,105,105,110,140,140,125,115,110)
+jun[:,3] = (110,105,100,100,100,110,115,125,145,115,105,100,95,90,95,95,100,105,120,145,150,145,130,115)
+jun[:,4] = (125,120,115,115,115,120,130,150,145,125,105,100,80,80,80,80,95,115,135,155,175,160,150,140)
+jun[:,5] = (130,125,120,115,110,105,105,105,100,95,85,75,75,60,55,60,85,100,110,140,145,145,140,135)
+jun[:,6] = (125,115,105,105,105,105,100,105,100,100,95,95,80,50,55,50,75,95,105,140,150,150,140,120)
+
+jun_mean = np.mean(jun, axis = 1)
+jun_tot = np.sum(jun_mean)
+index_jun = np.array([hour / jun_tot for hour in jun_mean])
+weight_jun = np.array([index / index_jun.mean() for index in index_jun])
+
+#%% Jul ############
+
+jul = np.zeros((24,7))
+jul[:,0] = (105,100,100,100,95,100,110,110,110,105,100,105,100,100,105,115,125,130,140,160,170,160,140,120)
+jul[:,1] = (105,100,100,95,95,95,100,110,110,110,105,105,100,100,115,120,130,145,155,180,190,165,140,120)
+jul[:,2] = (105,100,100,100,100,100,105,120,120,120,110,110,105,110,125,140,150,155,195,205,195,170,130,115)
+jul[:,3] = (105,105,100,95,95,100,105,115,115,120,110,120,105,105,120,120,125,140,140,165,185,155,135,115)
+jul[:,4] = (110,105,100,100,100,100,110,110,115,110,110,110,100,100,105,110,115,120,135,145,145,135,120,110)
+jul[:,5] = (105,105,105,100,100,100,100,100,105,105,100,95,90,90,90,90,95,105,125,130,140,135,125,110)
+jul[:,6] = (115,110,105,100,100,100,100,100,100,95,90,85,85,75,85,90,95,100,110,120,155,160,125,120)
+
+jul_mean = np.mean(jul, axis = 1)
+jul_tot = np.sum(jul_mean)
+index_jul = np.array([hour / jul_tot for hour in jul_mean])
+weight_jul = np.array([index / index_jul.mean() for index in index_jul]) 
+
+#%% Aug ############
+
+aug = np.zeros((24,7))
+aug[:,0] = (100,100,100,100,95,100,105,110,105,100,100,100,100,95,100,100,100,100,110,145,150,130,110,105)
+aug[:,1] = (110,105,105,100,105,105,105,110,105,105,100,95,90,80,85,95,105,105,110,125,140,130,125,105)
+aug[:,2] = (105,100,100,95,95,100,110,120,130,110,100,95,90,90,90,95,100,110,120,150,155,140,120,110)
+aug[:,3] = (115,110,110,105,105,105,115,125,120,115,100,95,95,90,90,95,100,110,115,140,145,140,125,110)
+aug[:,4] = (115,110,105,110,110,110,115,115,115,115,110,100,90,85,85,90,100,115,130,155,155,130,115,115)
+aug[:,5] = (110,105,105,100,100,100,100,100,100,90,85,70,70,55,55,55,100,100,120,140,145,140,125,115)
+aug[:,6] = (115,110,105,105,100,100,100,100,100,90,80,50,45,45,45,50,85,100,110,130,140,140,130,110)
+
+aug_mean = np.mean(aug,axis = 1)
+aug_tot = np.sum(aug_mean)
+index_aug = np.array([hour / aug_tot for hour in aug_mean])
+weight_aug = np.array([index / index_aug.mean() for index in index_aug])
+
+#%% Sep ############
+
+sep = np.zeros((24,7))
+sep[:,0] = (110,110,100,105,100,110,120,130,145,140,125,120,115,115,120,120,120,125,140,160,145,120,115,110)
+sep[:,1] = (120,115,110,110,110,110,120,130,160,150,130,120,110,105,110,115,120,130,130,175,145,125,115,110)
+sep[:,2] = (115,110,110,110,110,115,125,125,150,150,130,125,115,115,120,120,130,140,150,195,155,125,115,110)
+sep[:,3] = (110,110,105,100,100,105,120,120,135,130,125,120,110,100,110,115,115,120,145,180,150,115,110,110)
+sep[:,4] = (115,105,105,105,105,110,125,150,160,150,130,125,115,115,115,120,120,125,130,155,135,120,115,110)
+sep[:,5] = (115,115,115,110,105,115,120,115,120,120,115,105,95,85,90,100,110,115,130,155,145,130,115,115)
+sep[:,6] = (105,100,100,100,100,100,100,100,95,90,85,90,95,80,50,45,55,95,110,120,135,120,110,100)
+
+sep_mean = np.mean(sep, axis = 1)
+sep_tot = np.sum(sep_mean)
+index_sep = np.array([hour / sep_tot for hour in sep_mean])
+weight_sep = np.array([index / index_sep.mean() for index in index_sep])
+
+#%% Oct ############
+
+ott = np.zeros((24,7))
+ott[:,0] = (160,155,155,155,150,150,170,175,185,190,195,175,155,160,165,190,190,195,210,240,205,175,160,155)
+ott[:,1] = (150,145,140,140,140,140,165,185,190,170,160,150,150,150,150,150,160,175,195,205,175,160,150,150)
+ott[:,2] = (140,140,130,135,140,140,150,165,165,155,150,145,140,140,145,150,150,170,190,205,175,140,135,130)
+ott[:,3] = (140,140,135,130,130,135,150,190,200,180,160,150,150,145,150,150,150,155,175,190,175,150,140,130)
+ott[:,4] = (130,125,120,120,120,130,135,150,170,160,160,150,130,130,135,135,140,135,150,165,140,135,130,130)
+ott[:,5] = (135,130,125,130,130,130,135,150,160,175,150,135,130,115,115,115,130,135,150,175,150,135,130,125)
+ott[:,6] = (150,150,150,145,145,145,150,155,150,150,145,135,110,90,80,105,140,145,150,165,155,150,145,140)
+
+ott_mean = np.mean(ott, axis = 1)
+ott_tot = np.sum(ott_mean)
+index_ott = np.array([hour / ott_tot for hour in ott_mean])
+weight_oct = np.array([index / index_ott.mean() for index in index_ott])
+
+#%% Nov ############
+
+nov = np.zeros((24,7))
+nov[:,0] = (125,125,125,120,120,115,120,130,145,130,120,115,110,110,125,130,135,150,150,145,130,125,120,115)
+nov[:,1] = (110,105,105,105,110,105,125,135,150,140,125,125,120,125,130,145,150,150,145,140,125,120,115,105)
+nov[:,2] = (105,100,90,85,95,105,115,125,130,125,115,110,95,105,120,125,150,155,160,155,145,125,120,120)
+nov[:,3] = (110,105,105,105,110,125,135,155,170,145,135,135,130,135,140,155,185,180,175,170,140,135,130,110)
+nov[:,4] = (105,100,95,90,95,105,125,140,140,135,115,110,100,105,125,130,155,165,165,160,140,130,125,120)
+nov[:,5] = (125,125,125,125,125,125,125,130,125,125,120,105,115,110,125,125,130,150,160,165,150,130,125,125)
+nov[:,6] = (125,125,120,120,125,125,125,125,125,120,115,120,115,115,110,120,115,130,130,125,125,125,120,120)
+
+nov_mean = np.mean(nov, axis = 1)
+nov_tot = np.sum(nov_mean)
+index_nov = np.array([hour / nov_tot for hour in nov_mean])
+weight_nov = np.array([index / index_nov.mean() for index in index_nov])
+
+# nov_mean2 = np.mean(nov, axis = 0)
+# index_nov = np.zeros((24,7))
+# for i in range(24):
+#     for j in range(7):
+#         index_nov[i,j] = nov[i,j] / nov_mean2[j] 
+# weight_nov2 = np.mean(index_nov, axis = 1)
+
+#%% Dec ############
+
+dec = np.zeros((24,7))
+dec[:,0] = (115,115,110,110,95,100,120,125,145,140,125,120,120,120,125,125,135,130,135,125,120,115,110,95)
+dec[:,1] = (90,90,85,85,85,95,115,125,155,140,130,125,120,120,125,140,155,150,140,135,115,115,115,100)
+dec[:,2] = (90,85,80,80,85,90,110,115,125,120,115,115,110,115,120,125,135,140,130,125,115,115,110,100)
+dec[:,3] = (100,95,90,95,95,105,115,125,155,150,150,125,115,115,120,125,140,140,135,125,115,115,110,100)
+dec[:,4] = (95,90,90,90,90,90,105,120,140,120,105,105,100,105,110,110,120,125,115,115,110,110,95,85)
+dec[:,5] = (90,85,70,65,55,55,80,95,100,95,90,85,80,75,85,95,105,115,115,115,110,105,100,100)
+dec[:,6] = (110,100,90,85,80,95,90,95,95,100,95,95,100,90,95,100,110,115,125,135,130,115,115,110)
+
+dec_mean = np.mean(dec, axis = 1)
+dec_tot = np.sum(dec_mean)
+index_dec = np.array([hour / dec_tot for hour in dec_mean])
+weight_dec = np.array([index / index_dec.mean() for index in index_dec])
+
+#%% Weight Tot #####
+
+weight = np.array((weight_jan, weight_feb, weight_mar, weight_apr, weight_may, weight_jun,
+                  weight_jul, weight_aug, weight_sep, weight_oct, weight_nov, weight_dec))
+
+dm=[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]  # duration of months [days]
+
+weight_jan_t = np.tile(weight_jan, dm[0])
+weight_feb_t = np.tile(weight_feb, dm[1])
+weight_mar_t = np.tile(weight_mar, dm[2])
+weight_apr_t = np.tile(weight_apr, dm[3])
+weight_may_t = np.tile(weight_may, dm[4])
+weight_jun_t = np.tile(weight_jun, dm[5])
+weight_jul_t = np.tile(weight_jul, dm[6])
+weight_aug_t = np.tile(weight_aug, dm[7])
+weight_sep_t = np.tile(weight_sep, dm[8])
+weight_oct_t = np.tile(weight_oct, dm[9])
+weight_nov_t = np.tile(weight_nov, dm[10])
+weight_dec_t = np.tile(weight_dec, dm[11])
+
+# Obtaining hourly array
+
+weight_hour = np.concatenate((weight_jan_t, weight_feb_t, weight_mar_t, weight_apr_t,
+                              weight_may_t, weight_jun_t, weight_jul_t, weight_aug_t,
+                              weight_sep_t, weight_oct_t, weight_nov_t, weight_dec_t))
+
+weight_frame = pd.DataFrame(weight_hour)
+weight_frame.to_csv('weight.csv', index = False, header = False)
+
+#%% Test ###########
+
+
+x = np.arange(1,25)
+# plt.figure(dpi = 1000, figsize = (10,5))
+# plt.subplot(1,2,1)
+# plt.plot(x, dec_mean)
+# for day in range(dec.shape[1]):
+#     plt.plot(x, dec[:,day], color = 'tab:grey', alpha = 0.3)
+# plt.xlabel('Hour [-]')
+# plt.ylabel('PUN [€/MWh]')
+# plt.xticks([1,7,13,20])
+# plt.grid()
+# plt.subplot(1,2,2)
+# plt.plot(x, weight_dec)
+# plt.xlabel('Hour [-]')
+# plt.ylabel('PUN Weight [-]')
+# plt.grid()
+# plt.xticks([1,7,13,20])
+# plt.tight_layout() 
+# plt.show()
+
+# month = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+# plt.figure(dpi = 1000, figsize = (9,6))
+# for i in range(weight.shape[0]):
+#     plt.plot(x, weight[i,:], label = month[i])
+# plt.xlabel('Hour [-]')
+# plt.ylabel('PUN Weight [-]')
+# plt.grid(axis = 'y', alpha = 0.5)
+# plt.xticks([1,7,13,20])
+# plt.legend(ncol = 3)
+# plt.show()