Add constraint to set soc cyclic over the year for storage_units in r…

…olling horizon

etrago/execute/market_optimization.py

@@ -105,7 +105,7 @@ def market_optimization(self):
             "overlap"
         ],
         solver_name=self.args["solver"],
-        extra_functionality=Constraints(self.args, False).functionality,
+        extra_functionality=Constraints(self.args, False, apply_on="market_model").functionality,
     )
 
     # Reset formulation to previous setting of args
@@ -120,7 +120,8 @@ def market_optimization(self):
 
 
 def optimize_with_rolling_horizon(
-    n, pre_market, snapshots=None, horizon=2, overlap=0, **kwargs
+    n, pre_market, snapshots, horizon, overlap, solver_name,
+    extra_functionality, **kwargs
 ):
     """
     Optimizes the network in a rolling horizon fashion.
@@ -200,6 +201,19 @@ def optimize_with_rolling_horizon(
             n.storage_units.state_of_charge_initial = (
                 n.storage_units_t.state_of_charge.loc[snapshots[start - 1]]
             )
+            print(i)
+            # Make sure that state of charge of batteries and pumped hydro
+            # plants are cyclic over the year by using the state_of_charges
+            # from the pre_market_model
+            if i==0:
+                n.storage_units.state_of_charge_initial = (
+                    pre_market.storage_units_t.state_of_charge.iloc[-1]
+                    )
+
+            elif i == len(starting_points)-1:
+                extra_functionality=Constraints(
+                    {"extra_functionality": []},
+                    False, apply_on="last_market_model").functionality,
 
         status, condition = n.optimize(sns, **kwargs)
 

etrago/tools/constraints.py

@@ -3181,10 +3181,57 @@ def split_dispatch_disaggregation_constraints_nmp(self, n, sns):
     # TODO: implementieren
 
 
+
+def fixed_storage_unit_soc_at_the_end(n, sns):
+    """
+    Defines energy balance constraints for storage units. In principal the
+    constraints states:
+
+    previous_soc + p_store - p_dispatch + inflow - spill == soc
+    """
+    from xarray import DataArray
+
+    sns = n.snapshots[-1]
+    m = n.model
+    c = "StorageUnit"
+    assets = n.df(c)
+    if assets.empty:
+        return
+
+    # elapsed hours
+    eh = n.snapshot_weightings.stores[sns]
+    # efficiencies
+    eff_stand = (1 - n.storage_units.standing_loss).pow(eh)
+    eff_dispatch = n.storage_units.efficiency_dispatch
+    eff_store = n.storage_units.efficiency_store
+
+    # SOC first hour of the year
+    post_soc = m[f"{c}-state_of_charge"].loc[n.snapshots[0]]
+
+    # SOC last hour of the year
+    soc = m[f"{c}-state_of_charge"].loc[sns]
+
+    lhs = [
+        (-1, soc),
+        (-1 / eff_dispatch * eh, m[f"{c}-p_dispatch"].loc[sns]),
+        (eff_store * eh, m[f"{c}-p_store"].loc[sns]),
+    ]
+
+    if f"{c}-spill" in m.variables:
+        lhs += [(-eh, m[f"{c}-spill"])]
+
+    # We add inflow and initial soc for noncyclic assets to rhs
+    rhs = DataArray((-n.storage_units.inflow).mul(eh))
+
+    lhs += [(eff_stand, post_soc)]
+
+    m.add_constraints(lhs, "=", rhs, name=f"{c}-energy_balance_end")
+
 class Constraints:
-    def __init__(self, args, conduct_dispatch_disaggregation):
+    def __init__(self, args, conduct_dispatch_disaggregation, apply_on="grid_model"):
         self.args = args
         self.conduct_dispatch_disaggregation = conduct_dispatch_disaggregation
+        self.apply_on = apply_on
 
     def functionality(self, network, snapshots):
         """Add constraints to pypsa-model using extra-functionality.
@@ -3211,6 +3258,9 @@ def functionality(self, network, snapshots):
                     len(snapshots) > 1500
                 ):
                     add_ch4_constraints_linopy(self, network, snapshots)
+
+                if self.apply_on=="last_market_model":
+                    fixed_storage_unit_soc_at_the_end(network, snapshots)
                 add_chp_constraints_linopy(network, snapshots)
             else:
                 add_chp_constraints_nmp(network)