Add operational reserve margin constraint analogous to GenX

Co-authored-by: FabianHofmann <hofmann@fias.uni-frankfurt.de>

config.default.yaml

@@ -44,6 +44,11 @@ electricity:
   co2base: 1.487e+9
   agg_p_nom_limits: data/agg_p_nom_minmax.csv
 
+  operational_reserve: # like https://genxproject.github.io/GenX/dev/core/#Reserves
+    epsilon_load: 0.02 # share of total load
+    epsilon_vres: 0.02 # share of total renewable supply
+    contingency: 4000 # fixed capacity in MW
+
   extendable_carriers:
     Generator: []
     StorageUnit: [] # battery, H2

doc/release_notes.rst

@@ -22,6 +22,9 @@ Energy Security Release (April 2022)
     * old: ``estimate_renewable_capacities_from_capacity_stats``
     * new: ``estimate_renewable_capacities``
 
+* Add operational reserve margin constraint analogous to `GenX implementation <https://genxproject.github.io/GenX/dev/core/#Reserves>`_.
+  Can be activated with config setting ``electricity: operational_reserve:``.
+
 * Add function to add global constraint on use of gas in :mod:`prepare_network`. This can be activated with `electricity: gaslimit:` given in MWh. 
 
 * The powerplants that have been shut down before 2021 are filtered out. 

scripts/solve_network.py

@@ -84,8 +84,9 @@ import pandas as pd
 import re
 
 import pypsa
-from pypsa.linopf import (get_var, define_constraints, linexpr, join_exprs,
-                          network_lopf, ilopf)
+from pypsa.linopf import (get_var, define_constraints, define_variables,
+                          linexpr, join_exprs, network_lopf, ilopf)
+from pypsa.descriptors import get_switchable_as_dense as get_as_dense
 
 from pathlib import Path
 from vresutils.benchmark import memory_logger
@@ -211,6 +212,73 @@ def add_SAFE_constraints(n, config):
     define_constraints(n, lhs, '>=', rhs, 'Safe', 'mintotalcap')
 
 
+def add_operational_reserve_margin_constraint(n, config):
+    
+    reserve_config = config["electricity"]["operational_reserve"]
+    EPSILON_LOAD = reserve_config["epsilon_load"]
+    EPSILON_VRES = reserve_config["epsilon_vres"]
+    CONTINGENCY = reserve_config["contingency"]
+
+    # Reserve Variables 
+    reserve = get_var(n, 'Generator', 'r')
+
+    # Share of extendable renewable capacities
+    ext_i = n.generators.query('p_nom_extendable').index
+    vres_i = n.generators_t.p_max_pu.columns
+    capacity_factor = n.generators_t.p_max_pu[vres_i.intersection(ext_i)]
+    renewable_capacity_variables = get_var(n, 'Generator', 'p_nom')[vres_i.intersection(ext_i)]
+
+    # Left-hand-side
+    lhs = (
+        linexpr((1, reserve)).sum(1) + 
+        linexpr((-EPSILON_VRES * capacity_factor, renewable_capacity_variables)).sum(1)
+    )
+
+    
+    # Total demand at t
+    demand =  n.loads_t.p.sum(1)
+    
+    # VRES potential of non extendable generators
+    capacity_factor = n.generators_t.p_max_pu[vres_i.difference(ext_i)]
+    renewable_capacity = n.generators.p_nom[vres_i.difference(ext_i)]
+    potential = (capacity_factor * renewable_capacity).sum(1)
+    
+    # Right-hand-side
+    rhs = EPSILON_LOAD * demand + EPSILON_VRES * potential + CONTINGENCY
+        
+    define_constraints(n, lhs, '>=', rhs, "Reserve margin")
+
+
+def update_capacity_constraint(n):
+    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 = get_var(n, 'Generator', 'p')
+    reserve = get_var(n, 'Generator', 'r')
+    
+    capacity_variable = get_var(n, 'Generator', 'p_nom')
+    capacity_fixed = n.generators.p_nom[fix_i]
+    
+    p_max_pu = get_as_dense(n, 'Generator', 'p_max_pu')
+    
+    lhs = linexpr((1, dispatch), (1, reserve))
+    lhs += linexpr((-p_max_pu[ext_i], capacity_variable)).reindex(columns=gen_i, fill_value='')
+    
+    rhs = (p_max_pu[fix_i] * capacity_fixed).reindex(columns=gen_i, fill_value=0)
+    
+    define_constraints(n, lhs, '<=', rhs, 'Generators', 'updated_capacity_constraint')
+
+
+def add_operational_reserve_margin(n, sns, config):
+
+    define_variables(n, 0, np.inf, 'Generator', 'r', axes=[sns, n.generators.index])
+
+    add_operational_reserve_margin_constraint(n, config)
+    
+    update_capacity_constraint(n)
+
+
 def add_battery_constraints(n):
     nodes = n.buses.index[n.buses.carrier == "battery"]
     if nodes.empty or ('Link', 'p_nom') not in n.variables.index:
@@ -236,6 +304,8 @@ def extra_functionality(n, snapshots):
         add_SAFE_constraints(n, config)
     if 'CCL' in opts and n.generators.p_nom_extendable.any():
         add_CCL_constraints(n, config)
+    if config["electricity"].get("operational_reserve"):
+        add_operational_reserve_margin(n, snapshots, config)
     for o in opts:
         if "EQ" in o:
             add_EQ_constraints(n, o)