Merge remote-tracking branch 'origin/Clustering_fix'

dispaset/GAMS/UCM_h.gms

@@ -1162,7 +1162,7 @@ OutputHeatSlack(au,z)=HeatSlack.L(au,z);
 OutputStorageInput(s,z)=StorageInput.L(s,z);
 OutputStorageInput(th,z)=StorageInput.L(th,z);
 OutputStorageLevel(s,z)=StorageLevel.L(s,z)/StorageCapacity(s);
-OutputStorageLevel(th,z)=StorageLevel.L(th,z);
+OutputStorageLevel(th,z)=StorageLevel.L(th,z)/StorageCapacity(th);
 OutputSystemCost(z)=SystemCost.L(z);
 OutputSpillage(s,z)  = Spillage.L(s,z) ;
 OutputShedLoad(n,z) = ShedLoad.L(n,z);

dispaset/postprocessing/plot.py

@@ -103,7 +103,7 @@ def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None,
             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)
-        axes[1].set_ylabel('Level [TWh]')
+        axes[1].set_ylabel('Level [GWh]')
         axes[1].yaxis.label.set_fontsize(12)
         line_SOC = mlines.Line2D([], [], color='black', alpha=alpha, label='Reservoir', linestyle=':')
 
@@ -146,20 +146,20 @@ 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
+    load_change = pd.Series(0, index=demand.index)
+    load_changed = False
     if isinstance(shedload, pd.Series):
         if not shedload.index.equals(demand.index):
             logging.critical('The shedload time series must have the same index as the demand')
             sys.exit(1)
         load_change += -shedload
-        load_changed=True
+        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
+        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')
@@ -174,7 +174,7 @@ def plot_dispatch(demand, plotdata, level=None, curtailment=None, shedload=None,
     elif level is None:
         plt.legend(handles=[line_demand] + [line_shedload] + patches[::-1], loc=4, bbox_to_anchor=(1.2, 0.5))
         axes[0].fill_between(demand.index, demand, reduced_demand, facecolor="none", hatch="X", edgecolor="k",
-                         linestyle='dashed')
+                             linestyle='dashed')
     else:
         plt.legend(title='Dispatch for ' + demand.name[1], handles=[line_demand] + [line_shedload] + [line_SOC] +
                                                                    patches[::-1], loc=4, bbox_to_anchor=(1.2, 0.5))
@@ -299,7 +299,7 @@ def plot_energy_zone_fuel(inputs, results, PPindicators):
     ax.barh(demand, left=ax.get_xticks() - 0.4, width=[0.8] * len(demand), height=ax.get_ylim()[1] * 0.005, linewidth=2,
             color='k')
     plt.show()
-    return ax
+    return GenPerZone
 
 
 def plot_zone_capacities(inputs, plot=True):
@@ -356,13 +356,14 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
 
     aggregation = False
     if 'OutputStorageLevel' in results:
-        lev = filter_by_zone(results['OutputStorageLevel'], inputs, z) / 1E6  # TWh
+        lev = filter_by_zone(results['OutputStorageLevel'], inputs, z)
+        lev = lev * inputs['units']['StorageCapacity'].loc[lev.columns] / 1e3  # GWh of storage
         level = filter_by_storage(lev, inputs, StorageSubset='s')
-        levels = pd.DataFrame(index=results['OutputStorageLevel'].index,columns=inputs['sets']['t'])
-        for t in ['HDAM','HPHS','BEVS','BATS']:
-            temp = filter_by_tech(level,inputs,t)
+        levels = pd.DataFrame(index=results['OutputStorageLevel'].index, columns=inputs['sets']['t'])
+        for t in ['HDAM', 'HPHS', 'BEVS', 'BATS']:
+            temp = filter_by_tech(level, inputs, t)
             levels[t] = temp.sum(axis=1)
-        levels.dropna(axis=1,inplace=True)
+        levels.dropna(axis=1, inplace=True)
         for col in levels.columns:
             if levels[col].max() == 0 and levels[col].min() == 0:
                 del levels[col]
@@ -383,7 +384,6 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
     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 + demand_flex, columns=[('DA', z)])
     demand = demand[('DA', z)]
@@ -396,7 +396,7 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
     else:
         shed_load = pd.Series(0, index=demand.index) / 1000  # GW
     if 'OutputDemandModulation' in results and z in results['OutputDemandModulation']:
-        shifted_load = -results['OutputDemandModulation'][z] / 1000 # GW
+        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
@@ -411,7 +411,8 @@ def plot_zone(inputs, results, z='', rng=None, rug_plot=True):
     else:
         curtailment = None
 
-    plot_dispatch(demand, plotdata, level, curtailment=curtailment, shedload = shed_load, shiftedload=shifted_load, rng=rng, alpha=0.5)
+    plot_dispatch(demand, plotdata, level, curtailment=curtailment, shedload=shed_load, shiftedload=shifted_load,
+                  rng=rng, alpha=0.5)
 
     # Generation plot:
     if rug_plot:
@@ -445,4 +446,4 @@ def storage_levels(inputs, results):
                                title='Difference between the calculated storage Levels and the (imposed) minimum level')
         ax.set_ylabel('%')
 
-    return True
+    return True

dispaset/postprocessing/postprocessing.py

@@ -18,8 +18,6 @@
 from .data_handler import ds_to_df
 
 
-
-
 def get_load_data(inputs, z):
     """ 
     Get the load curve, the residual load curve, and the net residual load curve of a specific zone
@@ -404,11 +402,11 @@ def get_indicators_powerplant(inputs, results):
 
 
 def CostExPost(inputs,results):
-    '''
+    """
     Ex post computation of the operational costs with plotting. This allows breaking down
-    the cost into its different components and check that it matches with the objective 
+    the cost into its different components and check that it matches with the objective
     function from the optimization.
-    
+
     The cost objective function is the following:
              SystemCost(i)
              =E=
@@ -425,9 +423,9 @@ def CostExPost(inputs,results):
              +0.7*Config("ValueOfLostLoad","val")*sum(u,LL_RampUp(u,i)+LL_RampDown(u,i))
              +Config("CostOfSpillage","val")*sum(s,spillage(s,i));
 
-             
+
     :returns: tuple with the cost components and their cumulative sums in two dataframes.
-    '''
+    """
     import datetime
 
     dfin = inputs['param_df']
@@ -555,7 +553,6 @@ def get_units_operation_cost(inputs, results):
 
     Main Author: @AbdullahAlawad
     """
-
     datain = ds_to_df(inputs)
 
     #DataFrame with startup times for each unit (1 for startup)
@@ -610,4 +607,4 @@ def get_units_operation_cost(inputs, results):
 
     UnitOperationCost = FiexedCost+StartUpCost+ShutDownCost+RampUpCost+RampDownCost+VariableCost
 
-    return UnitOperationCost
+    return UnitOperationCost

dispaset/preprocessing/build.py

@@ -164,12 +164,24 @@ def build_single_run(config, profiles=None):
     :profiles:      if turned on reservoir levels are overwritten by newly calculated ones 
                     from the mid term scheduling simulations
     """
-    if profiles is None:
-        ReservoirLevels = UnitBasedTable(plants_sto,'ReservoirLevels',config,fallbacks=['Unit','Technology','Zone'],default=0)
-    else:
-        ReservoirLevels = UnitBasedTable(plants_sto,'ReservoirLevels',config,fallbacks=['Unit','Technology','Zone'],default=0)
-        MidTermSchedulingProfiles = profiles
-        ReservoirLevels.update(MidTermSchedulingProfiles)
+    ReservoirLevels = UnitBasedTable(plants_sto,'ReservoirLevels',config,fallbacks=['Unit','Technology','Zone'],default=0)
+    if profiles is not None:
+        if config['HydroScheduling'] == 'Zonal':
+            if config['SimulationType'] == 'Standard':
+                logging.critical("SimulationType: 'Standard' and HydroScheduling: 'Zonal' not supported! Please"
+                                 " choose different input options")
+                sys.exit(1)
+            else:
+                # Remove the unit number (e.g. [1] - xxxxx)
+                profiles = profiles.rename(columns={col: col.split(' - ')[1] for col in profiles.columns})
+                logging.info('Temporary MTS profile names trimmed back to original unit names')
+
+            if all(profiles.columns.isin(ReservoirLevels.columns)):
+                ReservoirLevels.update(profiles)
+                logging.info('New reservoir levels from Mid-term scheduling are now imposed instead of historic values')
+            else:
+                logging.critical('MTS and finalTS column names do not match!')
+                sys.exit(1)
 
     ReservoirScaledInflows = UnitBasedTable(plants_sto,'ReservoirScaledInflows',config,fallbacks=['Unit','Technology','Zone'],default=0)
     HeatDemand = UnitBasedTable(plants_heat,'HeatDemand',config,fallbacks=['Unit'],default=0)
@@ -318,11 +330,30 @@ def build_single_run(config, profiles=None):
                'PriceTransmission':PriceTransmission}
 
     # Merge the following time series with weighted averages
-    for key in ['ScaledInflows','ReservoirLevels','Outages','AvailabilityFactors','CostHeatSlack']:
+    for key in ['ScaledInflows','Outages','AvailabilityFactors','CostHeatSlack']:
         finalTS[key] = merge_series(plants, finalTS[key], mapping, tablename=key)
     # Merge the following time series by summing
     for key in ['HeatDemand']:
         finalTS[key] = merge_series(plants, finalTS[key], mapping, tablename=key, method='Sum')
+    # Merge the following time series by weighted average based on storage capacity
+    for key in ['ReservoirLevels']:
+        finalTS[key] = merge_series(plants, finalTS[key], mapping, tablename=key, method='StorageWeightedAverage')
+        # Update reservoir levels with newly computed ones from the mid-term scheduling
+        if profiles is not None:
+            if config['HydroScheduling'] != 'Zonal':
+                if all(profiles.columns.isin(finalTS[key].columns)):
+                    finalTS[key].update(profiles)
+                    logging.info('New reservoir levels from Mid-term scheduling are now imposed instead of historic values')
+                else:
+                    logging.critical('MTS and finalTS column names do not match!')
+                    sys.exit(1)
+            else:
+                if config['SimulationType'] == 'Standard':
+                    logging.critical("SimulationType: 'Standard' and HydroScheduling: 'Zonal' not supported! Please"
+                                     " choose different input options")
+                    sys.exit(1)
+                else:
+                    continue
 
 # Check that all times series data is available with the specified data time step:
     for key in FuelPrices:
@@ -464,7 +495,7 @@ def build_single_run(config, profiles=None):
 
     # Storage profile and initial state:
     for i, s in enumerate(sets['s']):
-        if profiles is not None:
+        if profiles is None:
             if (config['InitialFinalReservoirLevel'] == 0) or (config['InitialFinalReservoirLevel'] == ""):
                 if s in finalTS['ReservoirLevels']:
                     # get the time
@@ -479,9 +510,13 @@ def build_single_run(config, profiles=None):
                     parameters['StorageProfile']['val'][i, :] = 0.5
             else:
                 if (config['default']['ReservoirLevelInitial'] > 1) or (config['default']['ReservoirLevelFinal'] > 1):
-                    logging.warning(s + ': The initial or final reservoir levels are higher than its capacity!' )
-                parameters['StorageInitial']['val'][i] = config['default']['ReservoirLevelInitial'] * Plants_sto['StorageCapacity'][s]
-                parameters['StorageProfile']['val'][i, :] = config['default']['ReservoirLevelFinal']
+                    logging.critical(s + ': The initial or final reservoir levels are higher than its capacity!' )
+                    sys.exit(1)
+                else:
+                    parameters['StorageInitial']['val'][i] = config['default']['ReservoirLevelInitial'] * Plants_sto['StorageCapacity'][s]
+                    parameters['StorageProfile']['val'][i, :] = config['default']['ReservoirLevelFinal']
+                    logging.info(s + ': New initial reservoir level is set to: ' + str(config['default']['ReservoirLevelInitial']) +
+                                    ', new final reservoir level is set to: ' + str(config['default']['ReservoirLevelFinal']))
         else:
             if s in finalTS['ReservoirLevels']:
                 # get the time

dispaset/preprocessing/data_handler.py

@@ -239,6 +239,12 @@ def merge_series(plants, data, mapping, method='WeightedAverage', tablename=''):
                         value = value + subunits[idx] * np.maximum(1e-9, plants['PowerCapacity'][idx]*plants['Nunits'][idx])
                     P_j = np.sum(np.maximum(1e-9, plants['PowerCapacity'][oldindexes]*plants['Nunits'][oldindexes]))
                     merged[newunit] = value / P_j
+                elif method == 'StorageWeightedAverage':
+                    for idx in oldindexes:
+                       name = plants['Unit'][idx]
+                       value = value + subunits[idx] * np.maximum(1e-9, plants['STOCapacity'][idx]*plants['Nunits'][idx])
+                    P_j = np.sum(np.maximum(1e-9, plants['STOCapacity'][oldindexes]*plants['Nunits'][oldindexes]))
+                    merged[newunit] = value / P_j
                 elif method == 'Sum':
                     merged[newunit] = subunits.sum(axis=1)
                 else:
@@ -309,7 +315,7 @@ def load_time_series(config,path,header='infer'):
                          '. However, its length does not allow guessing its timestamps. Please use a 8760 elements time series')
             sys.exit(1)
 
-    if data.index.is_all_dates:   
+    if data.index.is_all_dates:
         data.index = data.index.tz_localize(None)   # removing locational data
         # Checking if the required index entries are in the data:
         common = data.index.intersection(config['idx'])
@@ -431,10 +437,10 @@ def load_config_excel(ConfigFile,AbsPath=True):
             config['default'][p] = sheet.cell_value(default[p], 5)
 
         #True/Falst values:
-        config['zones'] = read_truefalse(sheet, 225, 1, 246, 3)
-        config['zones'] = config['zones'] + read_truefalse(sheet, 225, 4, 246, 6)
-        config['mts_zones'] = read_truefalse(sheet, 225, 1, 246, 3, 2)
-        config['mts_zones'] = config['mts_zones'] + read_truefalse(sheet, 225, 4, 246, 6, 2)
+        config['zones'] = read_truefalse(sheet, 225, 1, 247, 3)
+        config['zones'] = config['zones'] + read_truefalse(sheet, 225, 4, 247, 6)
+        config['mts_zones'] = read_truefalse(sheet, 225, 1, 247, 3, 2)
+        config['mts_zones'] = config['mts_zones'] + read_truefalse(sheet, 225, 4, 247, 6, 2)
         config['ReserveParticipation'] = read_truefalse(sheet, 305, 1, 319, 3)
 
         # Set default values (for backward compatibility):