From 27dc9cedfdbe4aa716552d212eb2749f77ae0a00 Mon Sep 17 00:00:00 2001
From: lisazeyen <lisa.zeyen@web.de>
Date: Thu, 20 Apr 2023 09:35:53 +0200
Subject: [PATCH] fix bugs in reserve constraint

---
 scripts/solve_network.py | 30 +++++++++++++++++++-----------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index cfb95bfe..c08c7b09 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -43,6 +43,7 @@ from vresutils.benchmark import memory_logger
 
 logger = logging.getLogger(__name__)
 pypsa.pf.logger.setLevel(logging.WARNING)
+from pypsa.descriptors import get_switchable_as_dense as get_as_dense
 
 
 def add_land_use_constraint(n, config):
@@ -414,7 +415,7 @@ def add_operational_reserve_margin(n, sns, config):
         0, np.inf, coords=[sns, n.generators.index], name="Generator-r"
     )
     reserve = n.model["Generator-r"]
-    lhs = reserve.sum("Generator")
+    summed_reserve = reserve.sum("Generator")
 
     # Share of extendable renewable capacities
     ext_i = n.generators.query("p_nom_extendable").index
@@ -426,10 +427,10 @@ def add_operational_reserve_margin(n, sns, config):
             .loc[vres_i.intersection(ext_i)]
             .rename({"Generator-ext": "Generator"})
         )
-        lhs = lhs + (p_nom_vres * (-EPSILON_VRES * capacity_factor)).sum()
+        lhs = summed_reserve + (p_nom_vres * (-EPSILON_VRES * capacity_factor)).sum("Generator")
 
     # Total demand per t
-    demand = n.loads_t.p_set.sum(axis=1)
+    demand = get_as_dense(n, "Load", "p_set").sum(axis=1)
 
     # VRES potential of non extendable generators
     capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)]
@@ -441,17 +442,24 @@ def add_operational_reserve_margin(n, sns, config):
 
     n.model.add_constraints(lhs >= rhs, name="reserve_margin")
 
+    # additional contraint that capacity is not exceeded
+    gen_i = n.generators.index
+    ext_i = n.generators.query('p_nom_extendable').index
+    fix_i = n.generators.query('not p_nom_extendable').index
+
+    dispatch = n.model["Generator-p"]
     reserve = n.model["Generator-r"]
 
-    lhs = n.model.constraints["Generator-fix-p-upper"].lhs
-    lhs = lhs + reserve.loc[:, lhs.coords["Generator-fix"]].drop("Generator")
-    rhs = n.model.constraints["Generator-fix-p-upper"].rhs
-    n.model.add_constraints(lhs <= rhs, name="Generator-fix-p-upper-reserve")
+    capacity_variable = n.model["Generator-p_nom"].rename({"Generator-ext": "Generator"})
+    capacity_fixed = n.generators.p_nom[fix_i]
 
-    lhs = n.model.constraints["Generator-ext-p-upper"].lhs
-    lhs = lhs + reserve.loc[:, lhs.coords["Generator-ext"]].drop("Generator")
-    rhs = n.model.constraints["Generator-ext-p-upper"].rhs
-    n.model.add_constraints(lhs >= rhs, name="Generator-ext-p-upper-reserve")
+    p_max_pu = get_as_dense(n, 'Generator', 'p_max_pu')
+
+    lhs = dispatch + reserve - capacity_variable * p_max_pu[ext_i]
+
+    rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i, fill_value=0)
+
+    n.model.add_constraints(lhs <= rhs, name="updated_capacity_constraint")
 
 
 def add_battery_constraints(n):