Introduce PTES and TES constraints

config/config.default.yaml

@@ -1223,6 +1223,12 @@ plotting:
     residential urban decentral water tanks discharger: '#baac9e'
     services rural water tanks discharger: '#bbc2b8'
     services urban decentral water tanks discharger: '#bdd8d3'
+    water pits: "#cc826a"
+    water pits charger: "#b36a5e"
+    water pits discharger: "#b37468"
+    urban central water pits: "#d96f4c"
+    urban central water pits charger: "#a85d47"
+    urban central water pits discharger: "#b36452"
     # heat demand
     Heat load: '#cc1f1f'
     heat: '#cc1f1f'

doc/release_notes.rst

@@ -14,6 +14,7 @@ Upcoming Release
 * Bugfix: Geothermal heat potentials are now restricted to those in close proximity to future district heating areas as projected by Manz et al. 2024. Includes a refactoring change: Building of generic technical potentials from heat utilisation potentials was changed to specific computation of geothermal heat potentials.
 
 - ...
+- Implemented an energy-to-power ratio constraint for TES, linking storage capacity to the corresponding charger capacity. Additionally, chargers and dischargers are now sized identically through a unified constraint.
 
 - Bug fix: Added setting ``run: use_shadow_directory:`` (default: ``true``) which sets the ``shadow`` parameter of the snakemake workflow. Configuring to ``true`` sets snakemake ``shadow`` parameter to ``shalloow``, ``false`` to `Ǹone``. Should be set to ``false`` for those cases, where snakemake has an issue with finding missing input/output files in solving rules.
 

scripts/prepare_sector_network.py

@@ -2607,12 +2607,20 @@ def add_heat(
                 unit="MWh_th",
             )
 
+            etpr_water_tanks = costs.at[
+                heat_system.central_or_decentral + " water tank storage",
+                "energy to power ratio",
+            ]
+
             n.add(
                 "Link",
                 nodes + f" {heat_system} water tanks charger",
                 bus0=nodes + f" {heat_system} heat",
                 bus1=nodes + f" {heat_system} water tanks",
-                efficiency=costs.at["water tank charger", "efficiency"],
+                efficiency=costs.at[
+                    heat_system.central_or_decentral + " water tank charger",
+                    "efficiency",
+                ],
                 carrier=f"{heat_system} water tanks charger",
                 p_nom_extendable=True,
             )
@@ -2623,10 +2631,17 @@ def add_heat(
                 bus0=nodes + f" {heat_system} water tanks",
                 bus1=nodes + f" {heat_system} heat",
                 carrier=f"{heat_system} water tanks discharger",
-                efficiency=costs.at["water tank discharger", "efficiency"],
+                efficiency=costs.at[
+                    heat_system.central_or_decentral + " water tank discharger",
+                    "efficiency",
+                ],
                 p_nom_extendable=True,
             )
 
+            n.links.loc[nodes + f" {heat_system} water tanks charger", "etpr"] = (
+                etpr_water_tanks
+            )
+
             tes_time_constant_days = options["tes_tau"][
                 heat_system.central_or_decentral
             ]
@@ -2647,6 +2662,63 @@ def add_heat(
                 ],
             )
 
+            if heat_system == HeatSystem.URBAN_CENTRAL:
+                n.add("Carrier", f"{heat_system} water pits")
+
+                n.add(
+                    "Bus",
+                    nodes + f" {heat_system} water pits",
+                    location=nodes,
+                    carrier=f"{heat_system} water pits",
+                    unit="MWh_th",
+                )
+
+                etpr_water_pit = costs.at[
+                    "central water pit storage", "energy to power ratio"
+                ]
+
+                n.add(
+                    "Link",
+                    nodes + f" {heat_system} water pits charger",
+                    bus0=nodes + f" {heat_system} heat",
+                    bus1=nodes + f" {heat_system} water pits",
+                    efficiency=costs.at[
+                        "central water pit charger",
+                        "efficiency",
+                    ],
+                    carrier=f"{heat_system} water pits charger",
+                    p_nom_extendable=True,
+                )
+
+                n.add(
+                    "Link",
+                    nodes + f" {heat_system} water pits discharger",
+                    bus0=nodes + f" {heat_system} water pits",
+                    bus1=nodes + f" {heat_system} heat",
+                    carrier=f"{heat_system} water pits discharger",
+                    efficiency=costs.at[
+                        "central water pit discharger",
+                        "efficiency",
+                    ],
+                    p_nom_extendable=True,
+                )
+
+                n.links.loc[nodes + f" {heat_system} water pits charger", "etpr"] = (
+                    etpr_water_pit
+                )
+
+                n.add(
+                    "Store",
+                    nodes + f" {heat_system} water pits",
+                    bus=nodes + f" {heat_system} water pits",
+                    e_cyclic=True,
+                    e_nom_extendable=True,
+                    carrier=f"{heat_system} water pits",
+                    standing_loss=1 - np.exp(-1 / 24 / tes_time_constant_days),
+                    capital_cost=costs.at["central water pit storage", "fixed"],
+                    lifetime=costs.at["central water pit storage", "lifetime"],
+                )
+
         if options["resistive_heaters"]:
             key = f"{heat_system.central_or_decentral} resistive heater"
 

scripts/solve_network.py

@@ -33,6 +33,7 @@
 import sys
 from typing import Any
 
+import linopy
 import numpy as np
 import pandas as pd
 import pypsa
@@ -794,6 +795,74 @@ def add_operational_reserve_margin(n, sns, config):
     n.model.add_constraints(lhs <= rhs, name="Generator-p-reserve-upper")
 
 
+def add_TES_etpr_constraints(n):
+    """
+    Add a constraint for each TES storage unit enforcing:
+        Store-e_nom - etpr * Link-p_nom == 0
+
+    This constraint is only applied if both:
+      - n.config["sector"]["heating"] is True
+      - n.config["sector"]["tes"] is True
+    """
+    sector_config = n.config.get("sector", {})
+    if not (sector_config.get("heating", False) and sector_config.get("tes", False)):
+        return
+
+    tes_charger_bool = n.links.index.str.contains(
+        "water tanks charger|water pits charger"
+    )
+    tes_store_bool = n.stores.index.str.contains("water tanks|water pits")
+
+    chargers_ext = n.links[tes_charger_bool].query("p_nom_extendable").index
+    tes_ext = n.stores[tes_store_bool].query("e_nom_extendable").index
+    etpr_values = n.links.loc[chargers_ext, "etpr"].values
+
+    linear_expr_list = []
+    for discharger, tes, etpr_value in zip(chargers_ext, tes_ext, etpr_values):
+        char_var = n.model["Link-p_nom"].loc[discharger]
+        store_var = n.model["Store-e_nom"].loc[tes]
+        linear_expr = store_var - etpr_value * char_var
+        linear_expr_list.append(linear_expr)
+
+    # Merge the individual expressions
+    merged_expr = linopy.expressions.merge(
+        linear_expr_list, dim="Store-ext, Link-ext", cls=type(linear_expr_list[0])
+    )
+
+    n.model.add_constraints(merged_expr == 0, name="TES_etpr")
+
+
+def add_TES_charger_ratio_constraints(n):
+    """
+    Add constraints ensuring that for each TES unit, the charger and discharger are sized the same.
+
+        Link-p_nom(charger) - efficiency * Link-p_nom(discharger) == 0
+
+    This constraint is only applied if both:
+      - n.config["sector"]["heating"] is True
+      - n.config["sector"]["tes"] is True
+    """
+    sector_config = n.config.get("sector", {})
+    if not (sector_config.get("heating", False) and sector_config.get("tes", False)):
+        return
+
+    discharger_bool = n.links.index.str.contains(
+        "water tanks discharger|water pits discharger"
+    )
+    charger_bool = n.links.index.str.contains("water tanks charger|water pits charger")
+
+    dischargers_ext = n.links[discharger_bool].query("p_nom_extendable").index
+    chargers_ext = n.links[charger_bool].query("p_nom_extendable").index
+
+    eff = n.links.efficiency[dischargers_ext].values
+    lhs = (
+        n.model["Link-p_nom"].loc[chargers_ext]
+        - n.model["Link-p_nom"].loc[dischargers_ext] * eff
+    )
+
+    n.model.add_constraints(lhs == 0, name="TES_charger_ratio")
+
+
 def add_battery_constraints(n):
     """
     Add constraint ensuring that charger = discharger, i.e.
@@ -997,6 +1066,8 @@ def extra_functionality(
     ):
         add_solar_potential_constraints(n, config)
 
+    add_TES_charger_ratio_constraints(n)
+    add_TES_etpr_constraints(n)
     add_battery_constraints(n)
     add_lossy_bidirectional_link_constraints(n)
     add_pipe_retrofit_constraint(n)