Allow a storage unit to be energy neutral over its EnergyType period

eESSInventory links a unit's inventory to its own past and bounds it, but
nothing requires it to end a period where it started. For a store that is the
right default. For demand-side management it is not: shifted consumption is
meant to be recovered within a window, not displaced across the horizon, and a
unit ending each day away from where it started quietly relaxes the shifting it
represents.

eMinimumEnergy and eMaximumEnergy already bound the energy a unit produces over
an EnergyType period, using the same window. They cannot express neutrality,
because they compare output against an exogenous profile whereas neutrality ties
output to charge, and the level the two must agree on is decided by the
optimisation rather than given.

This reuses that period rather than introducing a second way to say daily or
weekly: EnergyNeutrality opts a unit in, EnergyType sets the block. Both columns
are optional and default to off, so cases written before them load unchanged.

Author: erikfilias

openTEPES/openTEPES_DataConfiguration.py

@@ -518,6 +518,9 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
     par['pStorageTimeStep']  = par['pStorageType' ].map(idxCycle   ).fillna(1)                                                                                          .astype('int')
     par['pOutflowsTimeStep'] = par['pOutflowsType'].map(idxOutflows).fillna(1).where(par['pEnergyOutflows'   ].sum()                                   > 0.0, other = 1).astype('int')
     par['pEnergyTimeStep']   = par['pEnergyType'  ].map(idxEnergy  ).fillna(1).where(par['pVariableMinEnergy'].sum() + par['pVariableMaxEnergy'].sum() > 0.0, other = 1).astype('int')
+    # Same period vocabulary, but not gated on the min/max energy profiles: neutrality needs a block
+    # length whether or not the unit also carries an energy bound.
+    par['pNeutralityTimeStep'] = par['pEnergyType'].map(idxEnergy).fillna(1).astype('int')
 
     par['pStorageTimeStep']  = pd.concat([par['pStorageTimeStep'], par['pOutflowsTimeStep'], par['pEnergyTimeStep']], axis=1).min(axis=1)
     # cycle time step can't exceed the stage duration
@@ -618,6 +621,8 @@ def Create_ESS_RES_Sets(mTEPES) -> None:
     par['pStorageTimeStep']     = par['pStorageTimeStep'].loc      [mTEPES.es   ]
     par['pOutflowsTimeStep']    = par['pOutflowsTimeStep'].loc     [mTEPES.es   ]
     par['pStorageType']         = par['pStorageType'].loc          [mTEPES.es   ]
+    par['pIndEnergyNeutrality'] = par['pIndEnergyNeutrality'].loc  [mTEPES.es   ]
+    par['pNeutralityTimeStep']  = par['pNeutralityTimeStep'].loc   [mTEPES.es   ]
 
     # separate positive and negative demands to avoid converting negative values to 0
     par['pDemandElecPos']  = par['pDemandElec'].where(par['pDemandElec'] >= 0.0, 0.0)
@@ -1113,6 +1118,8 @@ def filter_rows(df, set):
     mTEPES.pEnergyTimeStep       = Param(mTEPES.gg,    initialize=par['pEnergyTimeStep'].to_dict()           , within=PositiveIntegers,    doc='Unit energy cycle'                                   )
     mTEPES.pIniInventory         = Param(mTEPES.psnes, initialize=par['pIniInventory'].to_dict()             , within=NonNegativeReals,    doc='ESS Initial storage',                    mutable=True)
     mTEPES.pStorageType          = Param(mTEPES.es,    initialize=par['pStorageType'].to_dict()              , within=Any             ,    doc='ESS Storage type'                                    )
+    mTEPES.pIndEnergyNeutrality = Param(mTEPES.es, initialize=par['pIndEnergyNeutrality'].to_dict(), within=Binary, doc='Energy neutral over its EnergyType period')
+    mTEPES.pNeutralityTimeStep  = Param(mTEPES.es, initialize=par['pNeutralityTimeStep'].to_dict() , within=PositiveIntegers, doc='Energy neutrality period [load levels]')
     mTEPES.pGenLoInvest          = Param(mTEPES.eb,    initialize=par['pGenLoInvest'].to_dict()              , within=NonNegativeReals,    doc='Lower bound of the investment decision', mutable=True)
     mTEPES.pGenUpInvest          = Param(mTEPES.eb,    initialize=par['pGenUpInvest'].to_dict()              , within=NonNegativeReals,    doc='Upper bound of the investment decision', mutable=True)
     mTEPES.pGenLoRetire          = Param(mTEPES.gd,    initialize=par['pGenLoRetire'].to_dict()              , within=NonNegativeReals,    doc='Lower bound of the retirement decision', mutable=True)

openTEPES/openTEPES_InputData.py

@@ -381,6 +381,10 @@ def ProcessParameter(pDataFrame: pd.DataFrame, pTimeStep: int) -> pd.DataFrame:
     par['pEfficiency']                 = dfs['dfGeneration']  ['Efficiency'                ]                                                             #               ESS round-trip efficiency      [p.u.]
     par['pStorageType']                = dfs['dfGeneration']  ['StorageType'               ]                                                             #               ESS storage  type
     par['pOutflowsType']               = dfs['dfGeneration']  ['OutflowsType'              ]                                                             #               ESS outflows type
+    # Optional opt-in to energy neutrality. The period is EnergyType, read just below.
+    par['pIndEnergyNeutrality']        = (dfs['dfGeneration']['EnergyNeutrality']
+                                          if 'EnergyNeutrality' in dfs['dfGeneration'].columns
+                                          else pd.Series(0, index=dfs['dfGeneration'].index)).fillna(0).astype('int')
     par['pEnergyType']                 = dfs['dfGeneration']  ['EnergyType'                ]                                                             #               unit  energy type
     par['pRMaxReactivePower']          = dfs['dfGeneration']  ['MaximumReactivePower'      ] * 1e-3                                                      # rated maximum reactive power                 [Gvar]
     par['pGenLoInvest']                = dfs['dfGeneration']  ['InvestmentLo'              ]                                                             # Lower bound of the investment decision       [p.u.]

openTEPES/openTEPES_ModelFormulationElectricity.py

@@ -311,6 +311,25 @@ def eESSInventory(OptModel,n,es):
             return Constraint.Skip
     setattr(OptModel, f'eESSInventory_{p}_{sc}_{st}', Constraint(mTEPES.nesc, rule=eESSInventory, doc='ESS inventory balance [GWh]'))
 
+    # eMinimumEnergy and eMaximumEnergy bound the energy a unit produces over an EnergyType period
+    # against an exogenous profile. This is the same period and the same window, but it ties output
+    # to charge rather than to a profile: over the block the two must net to zero. That cannot be
+    # written as a min/max pair, since the level they must agree on is decided by the optimisation.
+    ne = [es for es in mTEPES.es if mTEPES.pIndEnergyNeutrality[es]]
+    if ne:
+        def eEnergyNeutrality(OptModel,n,es):
+            win = mTEPES.pNeutralityTimeStep[es]
+            i   = mTEPES.n.ord(n)
+            if i % win != 0 or (p,es) not in mTEPES.pes or (p,sc,st,n) not in mTEPES.s2n:
+                return Constraint.Skip
+            return sum(mTEPES.pDuration[p,sc,n2]() * (OptModel.vTotalOutput[p,sc,n2,es]
+                       - OptModel.vESSTotalCharge[p,sc,n2,es]) for n2 in n2list[i-win:i]) == 0.0
+        setattr(OptModel, f'eEnergyNeutrality_{p}_{sc}_{st}',
+                Constraint(mTEPES.n, ne, rule=eEnergyNeutrality, doc='energy neutrality over the EnergyType period [GWh]'))
+
+        if pIndLogConsole:
+            print('eEnergyNeutrality         ... ', len(getattr(OptModel, f'eEnergyNeutrality_{p}_{sc}_{st}')), ' rows')
+
     if pIndLogConsole:
         print('eESSInventory             ... ', len(getattr(OptModel, f'eESSInventory_{p}_{sc}_{st}')), ' rows')