Merge remote-tracking branch 'origin/FlexibleDemand'

ConfigFiles/configEU.yml

@@ -60,5 +60,4 @@ default: {CostHeatSlack: 50.0, CostLoadShedding: 400.0, LoadShedding: 0.05, Pric
   PriceOfLignite: 8.0, PriceOfNuclear: 3.0, PriceOfPeat: 8.0}
 
 #---------------- Multipliers
-modifiers: {Demand: 1.0, Solar: 1.0, Storage: 1.0, Wind: 1.0}
-
+modifiers: {Demand: 1.0, Solar: 1.0, Storage: 1.0, Wind: 1.0}

ConfigFiles/configEU_flex.yml

@@ -0,0 +1,64 @@
+#---------------- Simulation optinos (Input/output)
+SimulationDirectory:  'Simulations/simulationEU'
+WriteExcel: 0.0
+WriteGDX: 1.0
+WritePickle: 1.0
+
+GAMS_folder:  ' '
+cplex_path:  ''
+
+#----------------- Time range and horizon
+StartDate: !!python/tuple [2016, 1, 1, 0, 0, 0]
+StopDate: !!python/tuple [2016, 12, 31, 0, 0, 0]
+HorizonLength: 2
+LookAhead: 1
+
+#------------------ Model options
+Clustering: 1.0
+SimulationType:  'Integer clustering'  #LP, LP clustered
+ReserveCalculation:  'Generic'
+AllowCurtailment: 1.0
+
+#------------------- Timeseries
+Demand:  'Database/TotalLoadValue/##/1h/2016.csv'
+HeatDemand:  'Database/Heat_demand/2016.csv'
+LoadShedding:  ''
+
+Interconnections:  'Database/CrossBorderFlows/1h/2016.csv'
+NTC:  'Database/DayAheadNTC/1h/2016.csv'
+PowerPlantData:  'Database/PowerPlants/##/clustered.csv'
+Outages:  'Database/OutageFactors/##/2016.csv'
+
+RenewablesAF:  'Database/AvailabilityFactors/##/1h/2016.csv'
+
+ReservoirLevels:  'Database/HydroData/ScaledLevels/##/1h/2016.csv'
+ReservoirScaledInflows:  'Database/HydroData/ScaledInflows/##/1h/2016_profile_from_2012.csv'
+ShareOfFlexibleDemand: 'Database/TotalLoadValue/ShareFlexibleDemandEU.csv'
+
+#---------------- Prices
+PriceOfBiomass:  ''
+PriceOfBlackCoal:  ''
+PriceOfCO2:  ''
+PriceOfFuelOil:  ''
+PriceOfGas:  ''
+PriceOfLignite:  ''
+PriceOfNuclear:  ''
+PriceOfPeat:  ''
+CostHeatSlack:  ''
+CostLoadShedding:  ''
+
+#---------------- Countries
+countries: ['AT', 'BE', 'BG', 'CH', 'CZ', 'DE', 'DK', 'EE', 'EL', 'ES',
+   'FI', 'FR', 'HR','HU', 'IE', 'IT', 'LT', 'LV', 'NL', 'NO', 'PL',
+   'PT', 'RO', 'SE', 'SI', 'SK', 'UK']
+
+#---------------- Reserves
+ReserveParticipation: ['COMC', 'GTUR', 'HDAM', 'HPHS', 'ICEN', 'CAES', 'BATS', 'BEVS']
+
+#---------------- Default prices
+default: {CostHeatSlack: 50.0, CostLoadShedding: 400.0, LoadShedding: 0.05, PriceOfBiomass: 13.0,
+  PriceOfBlackCoal: 11.0, PriceOfCO2: 7.0, PriceOfFuelOil: 35.0, PriceOfGas: 20.0,
+  PriceOfLignite: 8.0, PriceOfNuclear: 3.0, PriceOfPeat: 8.0}
+
+#---------------- Multipliers
+modifiers: {Demand: 1.0, Solar: 1.0, Storage: 1.0, Wind: 1.0}

dispaset/GAMS/UCM_h.gms

@@ -49,6 +49,9 @@ $setglobal LPFormulation 0
 * (1 to retrieve 0 to not)
 $setglobal RetrieveStatus 0
 
+* Activate the flexible demand equations
+$setglobal ActivateFlexibleDemand 1
+
 *===============================================================================
 *Definition of   sets and parameters
 *===============================================================================
@@ -152,6 +155,9 @@ CostLoadShedding(n,h)            [EUR\MW]  Value of lost load
 LoadMaximum(u,h)                 [%]     Maximum load given AF and OF
 PowerMustRun(u,h)                [MW\u]    Minimum power output
 StorageFinalMin(s)               [MWh]   Minimum storage level at the end of the optimization horizon
+MaxFlexDemand(n)                 [MW]    Maximum value of the flexible demand parameter
+MaxOverSupply(n,h)               [MWh]   Maximum flexible demand accumultation
+AccumulatedOverSupply_inital(n)  [MWh]   Initial value of the flexible demand accumulation
 ;
 
 * Scalar variables necessary to the loop:
@@ -313,6 +319,7 @@ $If %LPFormulation% == 1 POSITIVE VARIABLES Committed (u,h) ; Committed.UP(u,h)
 $If not %LPFormulation% == 1 INTEGER VARIABLES Committed (u,h), StartUp(u,h), ShutDown(u,h) ; Committed.UP(u,h) = Nunits(u) ; StartUp.UP(u,h) = Nunits(u) ; ShutDown.UP(u,h) = Nunits(u) ;
 
 POSITIVE VARIABLES
+AccumulatedOverSupply(n,h) [MWh]   Accumulated oversupply due to the flexible demand
 CostStartUpH(u,h)          [EUR]   Cost of starting up
 CostShutDownH(u,h)         [EUR]   cost of shutting down
 CostRampUpH(u,h)           [EUR]   Ramping cost
@@ -345,6 +352,7 @@ WaterSlack(s)              [MWh]   Unsatisfied water level constraint
 
 free variable
 SystemCostD                ![EUR]   Total system cost for one optimization period
+DemandModulation(n,h)      [MW] Difference between the flexible demand and the baseline
 ;
 
 *===============================================================================
@@ -373,6 +381,11 @@ PowerMustRun(u,h)$(sum(tr,Technology(u,tr))>=1 and smin(n,Location(u,n)*(1-Curta
 * Part of the reserve that can be provided by offline quickstart units:
 K_QuickStart(n) = Config("QuickStartShare","val");
 
+* Flexible Demand
+MaxFlexDemand(n) = smax(h,Demand("Flex",n,h));
+MaxOverSupply(n,h) = Config("DemandFlexibility","val") * Demand("Flex",n,h);
+AccumulatedOverSupply_inital(n) = 0;
+
 $If %Verbose% == 1 Display RampStartUpMaximum, RampShutDownMaximum, CommittedInitial;
 
 $offorder
@@ -425,6 +438,11 @@ EQ_Flow_limits_upper
 EQ_Force_Commitment
 EQ_Force_DeCommitment
 EQ_LoadShedding
+EQ_Flexible_Demand
+EQ_Flexible_Demand_Max
+EQ_Flexible_Demand_Modulation_Min
+EQ_Flexible_Demand_Modulation_Max
+EQ_No_Flexible_Demand
 $If %RetrieveStatus% == 1 EQ_CommittedCalc
 ;
 
@@ -553,14 +571,52 @@ EQ_Demand_balance_DA(n,i)..
          sum(u,Power(u,i)*Location(u,n))
           +sum(l,Flow(l,i)*LineNode(l,n))
          =E=
-         Demand("DA",n,i)
+         Demand("DA",n,i) + Demand("Flex",n,i)
+         +DemandModulation(n,i)
          +sum(s,StorageInput(s,i)*Location(s,n))
          -ShedLoad(n,i)
          +sum(p2h,PowerConsumption(p2h,i)*Location(p2h,n))
          -LL_MaxPower(n,i)
          +LL_MinPower(n,i)
 ;
 
+* Energy balance at the level of the flexible demand, considered as a storage capacity
+EQ_Flexible_Demand(n,i)..
+         DemandModulation(n,i)
+         =e=
+         AccumulatedOverSupply(n,i)
+         - AccumulatedOverSupply_inital(n)$(ord(i) = 1)
+         - AccumulatedOverSupply(n,i-1)$(ord(i) > 1)
+;
+
+* The accumulated oversupply is limited by size of flexible demand (i.e. storage) capacity
+EQ_Flexible_Demand_max(n,i)..
+         AccumulatedOverSupply(n,i)
+         =l=
+         MaxOverSupply(n,i)
+;
+
+* The maximum downards demand modulation  at each time step is limited by the value of the flexible demand
+EQ_Flexible_Demand_Modulation_Min(n,i)..
+         DemandModulation(n,i)
+         =g=
+         -Demand("Flex",n,i)
+;
+
+* The upwards demand modulation is arbitrarily limited by the maximum recorded value of the flexible demand
+EQ_Flexible_Demand_Modulation_max(n,i)..
+         DemandModulation(n,i)
+         =l=
+         MaxFlexDemand(n) - Demand("flex",n,i)
+;
+
+* If flexible demand is not considered, set the accumulated over supply to zero:
+EQ_No_Flexible_Demand(n,i)..
+         AccumulatedOverSupply(n,i)
+         =e=
+         0
+;
+
 *Hourly demand balance in the upwards spinning reserve market for each node
 EQ_Demand_balance_2U(n,i)..
          sum((u,t),Reserve_2U(u,i)*Technology(u,t)*Reserve(t)*Location(u,n))
@@ -842,6 +898,11 @@ EQ_Flow_limits_upper,
 EQ_Force_Commitment,
 EQ_Force_DeCommitment,
 EQ_LoadShedding,
+$If %ActivateFlexibleDemand% == 1 EQ_Flexible_Demand,
+$If %ActivateFlexibleDemand% == 1 EQ_Flexible_Demand_Max,
+$if not %ActivateFlexibleDemand% == 1 EQ_No_Flexible_Demand,
+EQ_Flexible_Demand_Modulation_Min,
+EQ_Flexible_Demand_Modulation_Max,
 $If %RetrieveStatus% == 1 EQ_CommittedCalc
 /
 ;
@@ -912,6 +973,8 @@ if(UCM_SIMPLE.Modelstat <> 1 and UCM_SIMPLE.Modelstat <> 8 and not failed, Commi
          StorageInitial(s) =   sum(i$(ord(i)=LastKeptHour-FirstHour+1),StorageLevel.L(s,i));
          StorageInitial(chp) =   sum(i$(ord(i)=LastKeptHour-FirstHour+1),StorageLevel.L(chp,i));
 
+$If %ActivateFlexibleDemand% == 1 AccumulatedOverSupply_inital(n) = sum(i$(ord(i)=LastKeptHour-FirstHour+1),AccumulatedOverSupply.L(n,i));
+
 
 *Loop variables to display after solving:
 $If %Verbose% == 1 Display LastKeptHour,PowerInitial,CostStartUpH.L,CostShutDownH.L,CostRampUpH.L;
@@ -937,8 +1000,9 @@ OutputSystemCost(h)
 OutputSpillage(s,h)
 OutputShedLoad(n,h)
 OutputCurtailedPower(n,h)
+$If %ActivateFlexibleDemand% == 1 OutputDemandModulation(n,h)
 ShadowPrice(n,h)
-HeatShadowPrice(au,h)
+HeatShadowPrice(au,z)
 LostLoad_MaxPower(n,h)
 LostLoad_MinPower(n,h)
 LostLoad_2D(n,h)
@@ -967,6 +1031,7 @@ OutputSystemCost(z)=SystemCost.L(z);
 OutputSpillage(s,z)  = Spillage.L(s,z) ;
 OutputShedLoad(n,z) = ShedLoad.L(n,z);
 OutputCurtailedPower(n,z)=CurtailedPower.L(n,z);
+$If %ActivateFlexibleDemand% == 1 OutputDemandModulation(n,z)=DemandModulation.L(n,z);
 LostLoad_MaxPower(n,z)  = LL_MaxPower.L(n,z);
 LostLoad_MinPower(n,z)  = LL_MinPower.L(n,z);
 LostLoad_2D(n,z) = LL_2D.L(n,z);
@@ -993,6 +1058,7 @@ OutputSystemCost,
 OutputSpillage,
 OutputShedLoad,
 OutputCurtailedPower,
+$If %ActivateFlexibleDemand% == 1 OutputDemandModulation,
 OutputGenMargin,
 LostLoad_MaxPower,
 LostLoad_MinPower,

dispaset/postprocessing/data_handler.py

@@ -21,6 +21,7 @@
             'OutputSpillage':('u','h'),
             'OutputShedLoad':('n','h'),
             'OutputCurtailedPower':('n','h'),
+            'OutputDemandModulation':('n','h'),
             'LostLoad_MaxPower':('n','h'),
             'LostLoad_MinPower':('n','h'),
             'LostLoad_2D':('n','h'),
@@ -132,7 +133,7 @@ def get_sim_results(path='.', cache=None, temp_path=None, return_xarray=False, r
 
     keys_sparse = ['OutputPower','OutputPowerConsumption', 'OutputSystemCost', 'OutputCommitted', 'OutputCurtailedPower', 'OutputFlow',
                    'OutputShedLoad', 'OutputSpillage', 'OutputStorageLevel', 'OutputStorageInput', 'OutputHeat',
-                   'OutputHeatSlack']
+                   'OutputHeatSlack','OutputDemandModulation']
 
     # Setting the proper index to the result dataframes:
     from itertools import chain

dispaset/postprocessing/plot.py

@@ -11,7 +11,7 @@
 from .postprocessing import get_imports, get_plot_data, filter_by_zone, filter_by_tech, filter_by_storage
 
 
-def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None, rng=None,
+def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None, shiftedload=None, rng=None,
                   alpha=None, figsize=(13, 6)):
     """
     Function that plots the dispatch data and the reservoir level as a cumulative sum
@@ -87,19 +87,19 @@ def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None,
             sumplot_lev = level[cols_lvl[0:]].cumsum(axis=1)
             sumplot_lev['zero'] = 0
             sumplot_lev = sumplot_lev[['zero'] + sumplot_lev.columns[:-1].tolist()]
-            print('Its dataframe')
             for j in range(len(sumplot_lev.columns) - 1):
                 col3 = sumplot_lev.columns[j]
                 col4 = sumplot_lev.columns[j + 1]
                 rez_color = commons['colors'][col4]
                 rez_hatch = commons['hatches'][col4]
+                axes[1].plot(pdrng, sumplot_lev.loc[pdrng, col4], color='k', alpha=alpha, linestyle=':')
                 axes[1].fill_between(pdrng, sumplot_lev.loc[pdrng, col3], sumplot_lev.loc[pdrng, col4],
-                                     facecolor=rez_color, alpha=alpha, hatch=rez_hatch)
+                                     facecolor=rez_color, alpha=0.3)
                 labels.append(col4)
-                patches.append(mpatches.Patch(facecolor=rez_color, alpha=alpha, hatch=rez_hatch, label=col4))
+                patches.append(mpatches.Patch(facecolor=rez_color, alpha=0.3, label=col4))
                 colorlist.append(rez_color)
         elif isinstance(level, pd.Series):
-            # Create right axis:
+            # Create lower axis:
             axes[1].plot(pdrng, level[pdrng], color='k', alpha=alpha, linestyle=':')
             axes[1].fill_between(pdrng, 0, level[pdrng],
                                  facecolor=commons['colors']['WAT'], alpha=.3)
@@ -146,20 +146,30 @@ def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None,
     axes[0].set_ylabel('Power [GW]')
     axes[0].yaxis.label.set_fontsize(12)
 
+    load_change = pd.Series(0,index=demand.index)
+    load_changed=False
     if isinstance(shedload, pd.Series):
         if not shedload.index.equals(demand.index):
-            logging.error('The shedload time series must have the same index as the demand')
+            logging.critical('The shedload time series must have the same index as the demand')
             sys.exit(1)
-        reduced_demand = demand - shedload
-        axes[0].plot(pdrng, reduced_demand[pdrng], color='k', alpha=alpha, linestyle='dashed')
-        line_shedload = mlines.Line2D([], [], color='black', alpha=alpha, label='Shed Load', linestyle='dashed')
+        load_change += -shedload
+        load_changed=True
+    if isinstance(shiftedload, pd.Series):
+        if not shiftedload.index.equals(demand.index):
+            logging.critical('The shiftedload time series must have the same index as the demand')
+            sys.exit(1)
+        load_change += -shiftedload
+        load_changed=True
+    reduced_demand = demand + load_change
+    axes[0].plot(pdrng, reduced_demand[pdrng], color='k', alpha=alpha, linestyle='dashed')
+    line_shedload = mlines.Line2D([], [], color='black', alpha=alpha, label='New load', linestyle='dashed')
 
     line_demand = mlines.Line2D([], [], color='black', label='Load')
     # plt.legend(handles=[line_demand] + patches[::-1], loc=4)
 
-    if shedload is None and level is None:
+    if not load_changed and level is None:
         plt.legend(handles=[line_demand] + patches[::-1], loc=4, bbox_to_anchor=(1.2, 0.5))
-    if shedload is None:
+    elif not load_changed:
         plt.legend(handles=[line_demand] + [line_SOC] + patches[::-1], loc=4, bbox_to_anchor=(1.2, 0.5))
     elif level is None:
         plt.legend(handles=[line_demand] + [line_shedload] + patches[::-1], loc=4, bbox_to_anchor=(1.2, 0.5))
@@ -368,9 +378,14 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
         demand_p2h = demand_p2h.sum(axis=1)
     else:
         demand_p2h = pd.Series(0, index=results['OutputPower'].index)
+    if ('Flex', z) in inputs['param_df']['Demand']:
+        demand_flex = inputs['param_df']['Demand'][('Flex', z)] / 1000
+    else:
+        demand_flex = pd.Series(0, index=results['OutputPower'].index)
 
+
     demand_da = inputs['param_df']['Demand'][('DA', z)] / 1000  # GW
-    demand = pd.DataFrame(demand_da + demand_p2h, columns=[('DA', z)])
+    demand = pd.DataFrame(demand_da + demand_p2h + demand_flex, columns=[('DA', z)])
     demand = demand[('DA', z)]
 
     sum_generation = plotdata.sum(axis=1)
@@ -380,22 +395,23 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
         shed_load = pd.Series(shed_load, index=demand.index).fillna(0)
     else:
         shed_load = pd.Series(0, index=demand.index) / 1000  # GW
-    diff = (sum_generation - demand + shed_load).abs()
+    if 'OutputDemandModulation' in results and z in results['OutputDemandModulation']:
+        shifted_load = -results['OutputDemandModulation'][z] / 1000 # GW
+        shifted_load = pd.Series(shifted_load, index=demand.index).fillna(0)
+    else:
+        shifted_load = pd.Series(0, index=demand.index) / 1000  # GW
+    diff = (sum_generation - demand + shifted_load + shed_load).abs()
+
     if diff.max() > 0.01 * demand.max():
         logging.critical('There is up to ' + str(
             diff.max() / demand.max() * 100) + '% difference in the instantaneous energy balance of zone ' + z)
 
-    if 'OutputCurtailedPower' in results and z in results['OutputCurtailedPower'] \
-            and 'OutputShedLoad' in results and z in results['OutputShedLoad']:
-        curtailment = results['OutputCurtailedPower'][z] / 1000  # GW
-        plot_dispatch(demand, plotdata, level, curtailment=curtailment, shedload = shed_load, rng=rng, alpha=0.8)
-    elif 'OutputCurtailedPower' in results and z in results['OutputCurtailedPower']:
+    if 'OutputCurtailedPower' in results and z in results['OutputCurtailedPower']:
         curtailment = results['OutputCurtailedPower'][z] / 1000  # GW
-        plot_dispatch(demand, plotdata, level, curtailment=curtailment, rng=rng, alpha=0.8)
-    elif 'OutputShedLoad' in results and z in results['OutputShedLoad']:
-        plot_dispatch(demand, plotdata, level, shedload = shed_load, rng=rng, alpha=0.8)
     else:
-        plot_dispatch(demand, plotdata, level, rng=rng, alpha=0.8)
+        curtailment = None
+
+    plot_dispatch(demand, plotdata, level, curtailment=curtailment, shedload = shed_load, shiftedload=shifted_load, rng=rng, alpha=0.5)
 
     # Generation plot:
     if rug_plot: