Dynamic central heating temperatures (#1206)

* feat: Add rule to build central heating temperature profiles, adjust build_cop_profiles accordingly

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* style: make new config settings prettier

* docs: add file headers

* docs: update configtables and module docstring

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* docs: update release notes

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* chore: use smooth ambient temperature through rolling window; remove default vals in class

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update doc/configtables/sector.csv

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Amos Schledorn <a.schledorn@tu-berlin.de>
Co-authored-by: Fabian Neumann <fabian.neumann@outlook.de>

config/config.default.yaml

@@ -449,19 +449,26 @@ sector:
       2045: 0.8
       2050: 1.0
     district_heating_loss: 0.15
-    # check these numbers!
+    supply_temperature_approximation:
-    forward_temperature:
+      max_forward_temperature:
         default: 90
         DK: 70
         SE: 70
         NO: 70
-      FI: 70
+      min_forward_temperature:
+        default: 68
+        DK: 54
+        SE: 54
+        NO: 54
       return_temperature:
         default: 50
         DK: 40
         SE: 40
         NO: 40
         FI: 40
+      lower_threshold_ambient_temperature: 0
+      upper_threshold_ambient_temperature: 10
+      rolling_window_ambient_temperature: 72
     heat_source_cooling: 6 #K
     heat_pump_cop_approximation:
       refrigerant: ammonia

doc/configtables/sector.csv

@@ -9,8 +9,13 @@ district_heating,--,,`prepare_sector_network.py <https://github.com/PyPSA/pypsa-
 -- potential,--,float,maximum fraction of urban demand which can be supplied by district heating
 -- progress,--,Dictionary with planning horizons as keys., Increase of today's district heating demand to potential maximum district heating share. Progress = 0 means today's district heating share. Progress = 1 means maximum fraction of urban demand is supplied by district heating
 -- district_heating_loss,--,float,Share increase in district heat demand in urban central due to heat losses
--- forward_temperature,°C,Dictionary with country codes as keys. One key must be 'default'.,Forward temperature in district heating
+-- supply_temperature_approximation,,,
--- return_temperature,°C,Dictionary with country codes as keys. One key must be 'default'.,Return temperature in district heating. Must be lower than forward temperature
+-- -- max_forward_temperature,°C,Dictionary with country codes as keys. One key must be 'default'., Max. forward temperature in district heating (if ambient temperature lower-or-equal `lower_threshold_ambient_temperature`)
+-- -- min_forward_temperature,°C,Dictionary with country codes as keys. One key must be 'default'., Min. forward temperature in district heating (if ambient temperature higher-or-equal `upper_threshold_ambient_temperature`)
+-- -- return_temperature,°C,Dictionary with country codes as keys. One key must be 'default'.,Return temperature in district heating. Must be lower than forward temperature
+-- -- lower_threshold_ambient_temperature,°C,float, Assume `max_forward_temperature` if ambient temperature is below this threshold
+-- -- upper_threshold_ambient_temperature,°C,float, Assume `min_forward_temperature` if ambient temperature is above this threshold
+-- -- rolling_window_ambient_temperature, h, int, Rolling window size for averaging ambient temperature when approximating supply temperature
 -- heat_source_cooling,K,float,Cooling of heat source for heat pumps
 -- heat_pump_cop_approximation,,,
 -- -- refrigerant,--,"{ammonia, isobutane}",Heat pump refrigerant assumed for COP approximation

doc/release_notes.rst

@@ -61,6 +61,8 @@ Upcoming Release
 
 * Update JRC-IDEES-2015 to `JRC-IDEES-2021 <https://publications.jrc.ec.europa.eu/repository/handle/JRC137809>`__. The reference year is changed from 2015 to 2019.
 
+* Made central heating supply temperatures dynamic based on an adaptation of a reference curve from Pieper et al. (2019) (https://www.sciencedirect.com/science/article/pii/S0360544219305857?via%3Dihub).
+
 * Added option to use country-specific district heating forward and return temperatures. Defaults to lower temperatures in Scandinavia.
 
 * Added unsustainable biomass potentials for solid, gaseous, and liquid biomass. The potentials can be phased-out and/or

rules/build_sector.smk

@@ -212,17 +212,72 @@ rule build_temperature_profiles:
         "../scripts/build_temperature_profiles.py"
 
 
+rule build_central_heating_temperature_profiles:
+    params:
+        max_forward_temperature_central_heating=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "max_forward_temperature",
+        ),
+        min_forward_temperature_central_heating=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "min_forward_temperature",
+        ),
+        return_temperature_central_heating=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "return_temperature",
+        ),
+        snapshots=config_provider("snapshots"),
+        lower_threshold_ambient_temperature=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "lower_threshold_ambient_temperature",
+        ),
+        upper_threshold_ambient_temperature=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "upper_threshold_ambient_temperature",
+        ),
+        rolling_window_ambient_temperature=config_provider(
+            "sector",
+            "district_heating",
+            "supply_temperature_approximation",
+            "rolling_window_ambient_temperature",
+        ),
+    input:
+        temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
+        regions_onshore=resources("regions_onshore_elec_s{simpl}_{clusters}.geojson"),
+    output:
+        central_heating_forward_temperature_profiles=resources(
+            "central_heating_forward_temperature_profiles_elec_s{simpl}_{clusters}.nc"
+        ),
+        central_heating_return_temperature_profiles=resources(
+            "central_heating_return_temperature_profiles_elec_s{simpl}_{clusters}.nc"
+        ),
+    resources:
+        mem_mb=20000,
+    log:
+        logs("build_central_heating_temperature_profiles_s{simpl}_{clusters}.log"),
+    benchmark:
+        benchmarks("build_central_heating_temperature_profiles/s{simpl}_{clusters}")
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_central_heating_temperature_profiles/run.py"
+
+
 rule build_cop_profiles:
     params:
         heat_pump_sink_T_decentral_heating=config_provider(
             "sector", "heat_pump_sink_T_individual_heating"
         ),
-        forward_temperature_central_heating=config_provider(
-            "sector", "district_heating", "forward_temperature"
-        ),
-        return_temperature_central_heating=config_provider(
-            "sector", "district_heating", "return_temperature"
-        ),
         heat_source_cooling_central_heating=config_provider(
             "sector", "district_heating", "heat_source_cooling"
         ),
@@ -232,6 +287,12 @@ rule build_cop_profiles:
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
         snapshots=config_provider("snapshots"),
     input:
+        central_heating_forward_temperature_profiles=resources(
+            "central_heating_forward_temperature_profiles_elec_s{simpl}_{clusters}.nc"
+        ),
+        central_heating_return_temperature_profiles=resources(
+            "central_heating_return_temperature_profiles_elec_s{simpl}_{clusters}.nc"
+        ),
         temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"),
         temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
         regions_onshore=resources("regions_onshore_elec_s{simpl}_{clusters}.geojson"),

scripts/build_central_heating_temperature_profiles/central_heating_temperature_approximator.py

@@ -0,0 +1,218 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import pandas as pd
+import xarray as xr
+
+
+class CentralHeatingTemperatureApproximator:
+    """
+    A class to approximate central heating temperatures based on ambient
+    temperature.
+
+    Attributes
+    ----------
+    ambient_temperature : xr.DataArray
+        The ambient temperature data.
+    max_forward_temperature : xr.DataArray
+        The maximum forward temperature.
+    min_forward_temperature : xr.DataArray
+        The minimum forward temperature.
+    fixed_return_temperature : xr.DataArray
+        The fixed return temperature.
+    lower_threshold_ambient_temperature : float
+        Forward temperature is `max_forward_temperature` for ambient temperatures lower-or-equal this threshold.
+    upper_threshold_ambient_temperature : float
+        Forward temperature is `min_forward_temperature` for ambient temperatures higher-or-equal this threshold.
+    """
+
+    def __init__(
+        self,
+        ambient_temperature: xr.DataArray,
+        max_forward_temperature: float,
+        min_forward_temperature: float,
+        fixed_return_temperature: float,
+        lower_threshold_ambient_temperature: float,
+        upper_threshold_ambient_temperature: float,
+        rolling_window_ambient_temperature: int,
+    ) -> None:
+        """
+        Initialize the CentralHeatingTemperatureApproximator.
+
+        Parameters
+        ----------
+        ambient_temperature : xr.DataArray
+            The ambient temperature data.
+        max_forward_temperature : xr.DataArray
+            The maximum forward temperature.
+        min_forward_temperature : xr.DataArray
+            The minimum forward temperature.
+        fixed_return_temperature : xr.DataArray
+            The fixed return temperature.
+        lower_threshold_ambient_temperature : float
+            Forward temperature is `max_forward_temperature` for ambient temperatures lower-or-equal this threshold.
+        upper_threshold_ambient_temperature : float
+            Forward temperature is `min_forward_temperature` for ambient temperatures higher-or-equal this threshold.
+        rolling_window_ambient_temperature : int
+            Rolling window size for averaging ambient temperature.
+        """
+        self._ambient_temperature = ambient_temperature
+        self.max_forward_temperature = max_forward_temperature
+        self.min_forward_temperature = min_forward_temperature
+        self.fixed_return_temperature = fixed_return_temperature
+        self.lower_threshold_ambient_temperature = lower_threshold_ambient_temperature
+        self.upper_threshold_ambient_temperature = upper_threshold_ambient_temperature
+        self.rolling_window_ambient_temperature = rolling_window_ambient_temperature
+
+    def ambient_temperature_rolling_mean(self) -> xr.DataArray:
+        """
+        Property to get ambient temperature.
+
+        Returns
+        -------
+        xr.DataArray
+            Rolling mean of ambient temperature input.
+        """
+        # bfill to avoid NAs in the beginning
+        return (
+            self._ambient_temperature.rolling(
+                time=self.rolling_window_ambient_temperature
+            )
+            .mean(skip_na=True)
+            .bfill(dim="time")
+        )
+
+    @property
+    def forward_temperature(self) -> xr.DataArray:
+        """
+        Property to get dynamic forward temperature.
+
+        Returns
+        -------
+        xr.DataArray
+            Dynamic forward temperatures
+        """
+        return self._approximate_forward_temperature()
+
+    @property
+    def return_temperature(self) -> float:
+        """
+        Property to get return temperature.
+
+        Returns
+        -------
+        float
+            Return temperature.
+        """
+        return self._approximate_return_temperature()
+
+    def _approximate_forward_temperature(self) -> xr.DataArray:
+        """
+        Approximate dynamic forward temperature based on reference curve. Adapted from [Pieper et al. (2019)](https://doi.org/10.1016/j.energy.2019.03.165).
+
+        Returns
+        -------
+        xr.DataArray
+            Dynamic forward temperatures.
+        """
+
+        forward_temperature = xr.where(
+            self.ambient_temperature_rolling_mean()
+            <= self.lower_threshold_ambient_temperature,
+            self.max_forward_temperature,
+            xr.where(
+                self.ambient_temperature_rolling_mean()
+                >= self.upper_threshold_ambient_temperature,
+                self.min_forward_temperature,
+                self.min_forward_temperature
+                + (self.max_forward_temperature - self.min_forward_temperature)
+                * (
+                    self.upper_threshold_ambient_temperature
+                    - self.ambient_temperature_rolling_mean()
+                )
+                / (
+                    self.upper_threshold_ambient_temperature
+                    - self.lower_threshold_ambient_temperature
+                ),
+            ),
+        )
+        return forward_temperature
+
+    def _approximate_return_temperature(self) -> float:
+        """
+        Approximate return temperature.
+
+        Returns
+        -------
+        float
+            Return temperature.
+        """
+        return self.fixed_return_temperature
+
+    @property
+    def forward_temperature(self) -> xr.DataArray:
+        """
+        Property to get dynamic forward temperature.
+
+        Returns
+        -------
+        xr.DataArray
+            Dynamic forward temperatures.
+        """
+        return self._approximate_forward_temperature()
+
+    @property
+    def return_temperature(self) -> float:
+        """
+        Property to get return temperature.
+
+        Returns
+        -------
+        float
+            Return temperature.
+        """
+        return self._approximate_return_temperature()
+
+    def _approximate_forward_temperature(self) -> xr.DataArray:
+        """
+        Approximate dynamic forward temperature.
+
+        Returns
+        -------
+        xr.DataArray
+            Dynamic forward temperatures.
+        """
+        forward_temperature = xr.where(
+            self.ambient_temperature_rolling_mean()
+            <= self.lower_threshold_ambient_temperature,
+            self.max_forward_temperature,
+            xr.where(
+                self.ambient_temperature_rolling_mean()
+                >= self.upper_threshold_ambient_temperature,
+                self.min_forward_temperature,
+                self.min_forward_temperature
+                + (self.max_forward_temperature - self.min_forward_temperature)
+                * (
+                    self.upper_threshold_ambient_temperature
+                    - self.ambient_temperature_rolling_mean()
+                )
+                / (
+                    self.upper_threshold_ambient_temperature
+                    - self.lower_threshold_ambient_temperature
+                ),
+            ),
+        )
+        return forward_temperature
+
+    def _approximate_return_temperature(self) -> float:
+        """
+        Approximate return temperature.
+
+        Returns
+        -------
+        float
+            Return temperature.
+        """
+        return self.fixed_return_temperature

scripts/build_central_heating_temperature_profiles/run.py

@@ -0,0 +1,147 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+"""
+Approximate district heating forward and return temperature profiles based on
+ambient temperature. The method is based on a reference curve from Pieper et
+al. 2019, where for ambient temperatures below 0C, the highest possible forward
+temperature is assumed and vice versa for temperatures above 10C. Between these
+threshold levels, forward temperatures are linearly interpolated.
+
+By default, temperature levels are increased for non-Scandinavian countries.
+The default ratios between min. and max. forward temperatures is based on AGFW-Hauptbericht 2022.
+
+Relevant Settings
+-----------------
+
+.. code:: yaml
+    sector:
+        district_heating:
+            max_forward_temperature:
+            min_forward_temperature:
+            return_temperature:
+Inputs
+------
+- `resources/<run_name>/temp_air_total`: Air temperature
+
+Outputs
+-------
+- `resources/<run_name>/central_heating_temperature_profiles.nc`:
+
+References
+----------
+- Pieper, et al. (2019): "Assessment of a combination of three heat sources for heat pumps to supply district heating" (https://doi.org/10.1016/j.energy.2019.03.165).
+- AGFW (2022): "Hauptbericht 2022" (https://www.agfw.de/zahlen-und-statistiken/agfw-hauptbericht)
+"""
+
+import sys
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import xarray as xr
+from _helpers import set_scenario_config
+from central_heating_temperature_approximator import (
+    CentralHeatingTemperatureApproximator,
+)
+
+
+def get_country_from_node_name(node_name: str) -> str:
+    return node_name[:2]
+
+
+def map_temperature_dict_to_onshore_regions(
+    supply_temperature_by_country: dict,
+    regions_onshore: pd.Index,
+    snapshots: pd.DatetimeIndex,
+) -> xr.DataArray:
+    """
+    Map dictionary of temperatures to onshore regions.
+
+    Parameters:
+    ----------
+    supply_temperature_by_country : dictionary
+        Dictionary with temperatures as values and country keys as keys. One key must be named "default"
+    regions_onshore : pd.Index
+        Names of onshore regions
+    snapshots : pd.DatetimeIndex
+        Time stamps
+
+    Returns:
+    -------
+    xr.DataArray
+        The dictionary values mapped to onshore regions with onshore regions as coordinates.
+    """
+    return xr.DataArray(
+        [
+            [
+                (
+                    supply_temperature_by_country[get_country_from_node_name(node_name)]
+                    if get_country_from_node_name(node_name)
+                    in supply_temperature_by_country.keys()
+                    else supply_temperature_by_country["default"]
+                )
+                for node_name in regions_onshore.values
+            ]
+            # pass both nodes and snapshots as dimensions to preserve correct data structure
+            for _ in snapshots
+        ],
+        dims=["time", "name"],
+        coords={"time": snapshots, "name": regions_onshore},
+    )
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_cop_profiles",
+            simpl="",
+            clusters=48,
+        )
+
+    set_scenario_config(snakemake)
+
+    # map forward and return temperatures specified on country-level to onshore regions
+    regions_onshore = gpd.read_file(snakemake.input.regions_onshore)["name"]
+    snapshots = pd.date_range(freq="h", **snakemake.params.snapshots)
+    max_forward_temperature_central_heating_by_node_and_time: xr.DataArray = (
+        map_temperature_dict_to_onshore_regions(
+            supply_temperature_by_country=snakemake.params.max_forward_temperature_central_heating,
+            regions_onshore=regions_onshore,
+            snapshots=snapshots,
+        )
+    )
+    min_forward_temperature_central_heating_by_node_and_time: xr.DataArray = (
+        map_temperature_dict_to_onshore_regions(
+            supply_temperature_by_country=snakemake.params.min_forward_temperature_central_heating,
+            regions_onshore=regions_onshore,
+            snapshots=snapshots,
+        )
+    )
+    return_temperature_central_heating_by_node_and_time: xr.DataArray = (
+        map_temperature_dict_to_onshore_regions(
+            supply_temperature_by_country=snakemake.params.return_temperature_central_heating,
+            regions_onshore=regions_onshore,
+            snapshots=snapshots,
+        )
+    )
+
+    central_heating_temperature_approximator = CentralHeatingTemperatureApproximator(
+        ambient_temperature=xr.open_dataarray(snakemake.input.temp_air_total),
+        max_forward_temperature=max_forward_temperature_central_heating_by_node_and_time,
+        min_forward_temperature=min_forward_temperature_central_heating_by_node_and_time,
+        fixed_return_temperature=return_temperature_central_heating_by_node_and_time,
+        lower_threshold_ambient_temperature=snakemake.params.lower_threshold_ambient_temperature,
+        upper_threshold_ambient_temperature=snakemake.params.upper_threshold_ambient_temperature,
+        rolling_window_ambient_temperature=snakemake.params.rolling_window_ambient_temperature,
+    )
+
+    central_heating_temperature_approximator.forward_temperature.to_netcdf(
+        snakemake.output.central_heating_forward_temperature_profiles
+    )
+    central_heating_temperature_approximator.return_temperature.to_netcdf(
+        snakemake.output.central_heating_return_temperature_profiles
+    )

scripts/build_cop_profiles/run.py

@@ -53,47 +53,6 @@ from DecentralHeatingCopApproximator import DecentralHeatingCopApproximator
 from scripts.definitions.heat_system_type import HeatSystemType
 
 
-def map_temperature_dict_to_onshore_regions(
-    supply_temperature_by_country: dict,
-    regions_onshore: pd.Index,
-    snapshots: pd.DatetimeIndex,
-) -> xr.DataArray:
-    """
-    Map dictionary of temperatures to onshore regions.
-
-    Parameters:
-    ----------
-    supply_temperature_by_country : dictionary
-        Dictionary with temperatures as values and country keys as keys. One key must be named "default"
-    regions_onshore : pd.Index
-        Names of onshore regions
-    snapshots : pd.DatetimeIndex
-        Time stamps
-
-    Returns:
-    -------
-    xr.DataArray
-        The dictionary values mapped to onshore regions with onshore regions as coordinates.
-    """
-    return xr.DataArray(
-        [
-            [
-                (
-                    supply_temperature_by_country[get_country_from_node_name(node_name)]
-                    if get_country_from_node_name(node_name)
-                    in supply_temperature_by_country.keys()
-                    else supply_temperature_by_country["default"]
-                )
-                for node_name in regions_onshore.values
-            ]
-            # pass both nodes and snapshots as dimensions to preserve correct data structure
-            for _ in snapshots
-        ],
-        dims=["time", "name"],
-        coords={"time": snapshots, "name": regions_onshore},
-    )
-
-
 def get_cop(
     heat_system_type: str,
     heat_source: str,
@@ -154,20 +113,13 @@ if __name__ == "__main__":
     # map forward and return temperatures specified on country-level to onshore regions
     regions_onshore = gpd.read_file(snakemake.input.regions_onshore)["name"]
     snapshots = pd.date_range(freq="h", **snakemake.params.snapshots)
-    forward_temperature_central_heating_by_node_and_time: xr.DataArray = (
+    central_heating_forward_temperature: xr.DataArray = xr.open_dataarray(
-        map_temperature_dict_to_onshore_regions(
+        snakemake.input.central_heating_forward_temperature_profiles
-            supply_temperature_by_country=snakemake.params.forward_temperature_central_heating,
-            regions_onshore=regions_onshore,
-            snapshots=snapshots,
-        )
-    )
-    return_temperature_central_heating_by_node_and_time: xr.DataArray = (
-        map_temperature_dict_to_onshore_regions(
-            supply_temperature_by_country=snakemake.params.return_temperature_central_heating,
-            regions_onshore=regions_onshore,
-            snapshots=snapshots,
     )
+    central_heating_return_temperature: xr.DataArray = xr.open_dataarray(
+        snakemake.input.central_heating_return_temperature_profiles
     )
+
     cop_all_system_types = []
     for heat_system_type, heat_sources in snakemake.params.heat_pump_sources.items():
         cop_this_system_type = []
@@ -179,8 +131,8 @@ if __name__ == "__main__":
                 heat_system_type=heat_system_type,
                 heat_source=heat_source,
                 source_inlet_temperature_celsius=source_inlet_temperature_celsius,
-                forward_temperature_by_node_and_time=forward_temperature_central_heating_by_node_and_time,
+                forward_temperature_by_node_and_time=central_heating_forward_temperature,
-                return_temperature_by_node_and_time=return_temperature_central_heating_by_node_and_time,
+                return_temperature_by_node_and_time=central_heating_return_temperature,
             )
             cop_this_system_type.append(cop_da)
         cop_all_system_types.append(