From 44ea1cd25270ad887afae3c7292a55caa6261c19 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 11 Mar 2026 11:50:01 +0100 Subject: [PATCH 01/17] feat: first draft of layered PTES model (operational for hard-coded layers, dummy heat transfer coef and resistive boosting) --- config/config.default.yaml | 3 + rules/build_sector.smk | 66 ++-- rules/solve_myopic.smk | 1 + rules/solve_overnight.smk | 1 + rules/solve_perfect.smk | 1 + scripts/build_cop_profiles/run.py | 33 +- .../build_heat_source_utilisation_profiles.py | 71 ++-- .../ptes_approximator.py | 182 ++++++++++ .../ptes_temperature_approximator.py | 315 ------------------ scripts/build_ptes_operations/run.py | 105 ++---- scripts/definitions/heat_source.py | 15 +- scripts/definitions/heat_system.py | 15 +- scripts/prepare_sector_network.py | 280 +++++++++------- scripts/solve_network.py | 182 ++++++++++ 14 files changed, 692 insertions(+), 578 deletions(-) create mode 100644 scripts/build_ptes_operations/ptes_approximator.py delete mode 100644 scripts/build_ptes_operations/ptes_temperature_approximator.py diff --git a/config/config.default.yaml b/config/config.default.yaml index bbdb363f72..960a6e00e7 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -665,6 +665,9 @@ sector: bottom_temperature: 35 design_top_temperature: 90 design_bottom_temperature: 35 + layered: + num_layers: 1 + heat_transfer_coefficient: 0.01 # H [W/(m²·K)] ates: enable: false suitable_aquifer_types: diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 66d591f84a..29143389c2 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -529,7 +529,6 @@ def input_heat_source_temperature( replace_names: dict[str, str] = { "air": "air_total", "ground": "soil_total", - "ptes": "ptes_top_profiles", }, ) -> dict[str, str]: """ @@ -549,7 +548,7 @@ def input_heat_source_temperature( for heat sources that require temperature profiles (excludes constant temperature sources). """ - from scripts.definitions.heat_source import HeatSource + from scripts.definitions.heat_source import HeatSource, HeatSourceType heat_sources = set( config_provider("sector", "heat_sources", "urban central")(w) @@ -562,12 +561,11 @@ def input_heat_source_temperature( for heat_source_name in heat_sources: heat_source = HeatSource(heat_source_name) # Skip heat sources with temperatures defined in config (not from file) - if heat_source.temperature_from_config: + if ( + heat_source.temperature_from_config + or heat_source.source_type == HeatSourceType.STORAGE + ): continue - if heat_source_name == "ptes": - file_names[f"temp_{heat_source_name}"] = resources( - f"temp_{replace_names.get(heat_source_name, heat_source_name)}_base_s_{{clusters}}_{{planning_horizons}}.nc" - ) else: file_names[f"temp_{heat_source_name}"] = resources( f"temp_{replace_names.get(heat_source_name, heat_source_name)}_base_s_{{clusters}}.nc" @@ -575,12 +573,11 @@ def input_heat_source_temperature( return file_names -def input_ptes_bottom_temperature(w) -> dict[str, str]: +def input_ptes_operations(w) -> dict[str, str]: """ - Generate conditional input for PTES bottom temperature profiles. + Generate conditional input for PTES operations. - Only includes the input file if PTES is configured as a heat source - for urban central heating. + Only includes the input file if PTES is enabled. Parameters ---------- @@ -593,11 +590,10 @@ def input_ptes_bottom_temperature(w) -> dict[str, str]: Dictionary with "temp_ptes_bottom" key if PTES is a heat source, empty dict otherwise. """ - heat_sources = config_provider("sector", "heat_sources", "urban central")(w) - if "ptes" in heat_sources: + if config_provider("sector", "district_heating", "ptes", "enable")(w): return { - "temp_ptes_bottom": resources( - "temp_ptes_bottom_profiles_base_s_{clusters}_{planning_horizons}.nc" + "ptes_operations": resources( + "ptes_operations_base_s_{clusters}_{planning_horizons}.nc" ) } return {} @@ -678,9 +674,12 @@ rule build_cop_profiles: heat_source_temperatures=config_provider( "sector", "district_heating", "heat_source_temperatures" ), + ptes_enable=config_provider("sector", "district_heating", "ptes", "enable"), + ptes_layered=config_provider("sector", "district_heating", "ptes", "layered"), snapshots=config_provider("snapshots"), input: unpack(input_heat_source_temperature), + unpack(input_ptes_operations), central_heating_forward_temperature_profiles=resources( "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc" ), @@ -738,6 +737,7 @@ rule build_ptes_operations: "ptes", "design_bottom_temperature", ), + layered=config_provider("sector", "district_heating", "ptes", "layered"), input: central_heating_forward_temperature_profiles=resources( "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc" @@ -747,17 +747,8 @@ rule build_ptes_operations: ), regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), output: - ptes_top_temperature_profiles=resources( - "temp_ptes_top_profiles_base_s_{clusters}_{planning_horizons}.nc" - ), - ptes_bottom_temperature_profiles=resources( - "temp_ptes_bottom_profiles_base_s_{clusters}_{planning_horizons}.nc" - ), - ptes_e_max_pu_profiles=resources( - "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc" - ), - ptes_boost_per_discharge_profiles=resources( - "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc" + ptes_operations=resources( + "ptes_operations_base_s_{clusters}_{planning_horizons}.nc" ), resources: mem_mb=2000, @@ -784,7 +775,7 @@ rule build_heat_source_utilisation_profiles: ptes_enable=config_provider("sector", "district_heating", "ptes", "enable"), input: unpack(input_heat_source_temperature), - unpack(input_ptes_bottom_temperature), + unpack(input_ptes_operations), central_heating_forward_temperature_profiles=resources( "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc" ), @@ -1620,6 +1611,7 @@ rule prepare_sector_network: input: unpack(input_profile_offwind), unpack(input_heat_source_power), + unpack(input_ptes_operations), **rules.cluster_gas_network.output, **rules.build_gas_input_locations.output, snapshot_weightings=resources( @@ -1692,26 +1684,6 @@ rule prepare_sector_network: temp_soil_total=resources("temp_soil_total_base_s_{clusters}.nc"), temp_air_total=resources("temp_air_total_base_s_{clusters}.nc"), cop_profiles=resources("cop_profiles_base_s_{clusters}_{planning_horizons}.nc"), - ptes_e_max_pu_profiles=lambda w: ( - resources( - "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc" - ) - if config_provider("sector", "district_heating", "ptes", "enable")(w) - and config_provider( - "sector", "district_heating", "ptes", "temperature_dependent_capacity" - )(w) - else [] - ), - ptes_boost_per_discharge_profiles=lambda w: ( - resources( - "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc" - ) - if config_provider("sector", "district_heating", "ptes", "enable")(w) - and config_provider( - "sector", "district_heating", "ptes", "discharge_resistive_boosting" - )(w) - else [] - ), heat_source_direct_utilisation_profiles=lambda w: ( resources( "heat_source_direct_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc" diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index d90c332b5d..03f3289eae 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -116,6 +116,7 @@ rule solve_sector_network_myopic: ), custom_extra_functionality=input_custom_extra_functionality, input: + unpack(input_ptes_operations), network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc" ), diff --git a/rules/solve_overnight.smk b/rules/solve_overnight.smk index 270609ecf2..c0da468064 100644 --- a/rules/solve_overnight.smk +++ b/rules/solve_overnight.smk @@ -14,6 +14,7 @@ rule solve_sector_network: ), custom_extra_functionality=input_custom_extra_functionality, input: + unpack(input_ptes_operations), network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" ), diff --git a/rules/solve_perfect.smk b/rules/solve_perfect.smk index 59ec549254..4234bbdd0a 100644 --- a/rules/solve_perfect.smk +++ b/rules/solve_perfect.smk @@ -101,6 +101,7 @@ rule solve_sector_network_perfect: ), custom_extra_functionality=input_custom_extra_functionality, input: + unpack(input_ptes_operations), network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_brownfield_all_years.nc" ), diff --git a/scripts/build_cop_profiles/run.py b/scripts/build_cop_profiles/run.py index 6af8b02e8a..b9a1e6f316 100644 --- a/scripts/build_cop_profiles/run.py +++ b/scripts/build_cop_profiles/run.py @@ -65,7 +65,7 @@ from scripts.build_cop_profiles.decentral_heating_cop_approximator import ( DecentralHeatingCopApproximator, ) -from scripts.definitions.heat_source import HeatSource +from scripts.definitions.heat_source import HeatSource, HeatSourceType from scripts.definitions.heat_system_type import HeatSystemType @@ -167,6 +167,13 @@ def get_source_temperature( heat_source = HeatSource(heat_source_name) if heat_source.temperature_from_config: return snakemake_params["heat_source_temperatures"][heat_source_name] + elif heat_source.source_type == HeatSourceType.STORAGE: + # PTES layer temperatures are constants from the ptes_operations dataset + if heat_source_name.startswith("ptes layer"): + layer_idx = int(heat_source_name.split()[-1]) + return float(ptes_ds["layer_temperatures"].sel(layer=layer_idx).item()) + else: + return float(ptes_ds.attrs["top_temperature"]) else: return xr.open_dataarray(snakemake_input[f"temp_{heat_source_name}"]) @@ -268,8 +275,29 @@ def get_sink_inlet_temperature( snakemake.input.central_heating_return_temperature_profiles ) + # Load PTES operations dataset if enabled + if snakemake.params.ptes_enable: + ptes_ds = xr.open_dataset(snakemake.input.ptes_operations) + num_ptes_layers = int(ptes_ds.attrs["num_layers"]) + else: + ptes_ds = None + num_ptes_layers = 0 + + def expand_ptes_layers(heat_source_names): + """Expand 'ptes' into per-layer entries if multi-layer PTES is enabled.""" + expanded = [] + for name in heat_source_names: + if name == "ptes" and num_ptes_layers > 1: + expanded.extend(f"ptes layer {l}" for l in range(num_ptes_layers)) + else: + expanded.append(name) + return expanded + cop_all_system_types = [] for heat_system_type, heat_sources in snakemake.params.heat_sources.items(): + heat_sources = ( + expand_ptes_layers(heat_sources) if heat_sources else heat_sources + ) cop_this_system_type = [] if not heat_sources: cop_all_system_types.append( @@ -309,6 +337,9 @@ def get_sink_inlet_temperature( ) ) + if ptes_ds is not None: + ptes_ds.close() + cop_dataarray = xr.concat( cop_all_system_types, dim=pd.Index(snakemake.params.heat_sources.keys(), name="heat_system"), diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py index 4d18185fcf..4667a3a8ef 100644 --- a/scripts/build_heat_source_utilisation_profiles.py +++ b/scripts/build_heat_source_utilisation_profiles.py @@ -61,11 +61,24 @@ import xarray as xr from scripts._helpers import configure_logging, set_scenario_config -from scripts.definitions.heat_source import HeatSource +from scripts.definitions.heat_source import HeatSource, HeatSourceType logger = logging.getLogger(__name__) +def expand_heat_sources_with_ptes_layers( + heat_source_names: list[str], num_layers: int +) -> list[str]: + """Expand 'ptes' into per-layer entries if multi-layer PTES is enabled.""" + expanded = [] + for name in heat_source_names: + if name == "ptes" and num_layers > 1: + expanded.extend(f"ptes layer {l}" for l in range(num_layers)) + else: + expanded.append(name) + return expanded + + def get_source_temperature( snakemake_params: dict, snakemake_input: dict, heat_source_name: str ) -> float | xr.DataArray: @@ -96,6 +109,13 @@ def get_source_temperature( heat_source = HeatSource(heat_source_name) if heat_source.temperature_from_config: return snakemake_params["heat_source_temperatures"][heat_source_name] + elif heat_source.source_type == HeatSourceType.STORAGE: + # PTES layer temperatures are constants from the ptes_operations dataset + if heat_source_name.startswith("ptes layer"): + layer_idx = int(heat_source_name.split()[-1]) + return float(ptes_ds["layer_temperatures"].sel(layer=layer_idx).item()) + else: + return float(ptes_ds.attrs["top_temperature"]) else: if f"temp_{heat_source_name}" not in snakemake_input.keys(): raise ValueError( @@ -182,15 +202,14 @@ def get_preheater_utilisation_profile( def get_heat_pump_cooling( heat_source_name: str, default_heat_source_cooling: float, - snakemake_input: dict, return_temperature: xr.DataArray = None, ) -> float | xr.DataArray: """ Get the additional heat source cooling (temperature drop) through the heat pump for a heat source. For PTES, this equals the temperature difference between - return flow and bottom layers (return_temperature - bottom_temperature), which can - be time-varying. For other sources, uses the default constant value. + return flow and bottom layer (return_temperature - bottom_temperature). + For other sources, uses the default constant value. Parameters ---------- @@ -198,8 +217,6 @@ def get_heat_pump_cooling( Name of the heat source (e.g., 'ptes', 'geothermal', 'air'). default_heat_source_cooling : float Default heat source cooling in Kelvin, from config. - snakemake_input : dict - Snakemake input files, may contain PTES temperature profiles. return_temperature : xr.DataArray, optional District heating return temperature profiles in °C. Required for PTES. @@ -207,25 +224,16 @@ def get_heat_pump_cooling( ------- float | xr.DataArray Heat source cooling in Kelvin. Returns a float for most sources, - or a DataArray for PTES when temperatures vary with time. - - Raises - ------ - ValueError - If heat source is PTES but bottom temperature profile is not provided. + or a DataArray for PTES. """ - if heat_source_name == "ptes": - if "temp_ptes_bottom" not in snakemake_input.keys(): - raise ValueError( - "PTES heat source requires bottom temperature profile " - "(temp_ptes_bottom) to calculate heat source cooling." - ) + heat_source = HeatSource(heat_source_name) + if heat_source.source_type == HeatSourceType.STORAGE: if return_temperature is None: raise ValueError( "PTES heat source requires return_temperature to calculate heat pump cooling." ) - ptes_bottom_temperature = xr.open_dataarray(snakemake_input["temp_ptes_bottom"]) - return return_temperature - ptes_bottom_temperature + bottom_temperature = float(ptes_ds.attrs["bottom_temperature"]) + return return_temperature - bottom_temperature return default_heat_source_cooling @@ -234,7 +242,7 @@ def get_heat_pump_cooling( from scripts._helpers import mock_snakemake snakemake = mock_snakemake( - "build_cop_profiles", + "build_heat_source_utilisation_profiles", clusters=48, ) configure_logging(snakemake) @@ -244,12 +252,25 @@ def get_heat_pump_cooling( ptes_enable: bool = snakemake.params.ptes_enable # Validate PTES configuration - if ptes_enable and "ptes" not in heat_sources: + has_ptes = any(s == "ptes" or s.startswith("ptes layer") for s in heat_sources) + if ptes_enable and not has_ptes: raise ValueError( "PTES is enabled (district_heating.ptes.enable=true) but 'ptes' " - "is not in heat_sources.urban_central. PTES requires being listed in heat_sources to create the necessary buses and links for heat discharge to the 'urban central heat' bus." + "is not in heat_sources.urban_central. PTES requires being listed " + "in heat_sources to create the necessary buses and links for heat " + "discharge to the 'urban central heat' bus." ) + # Load PTES operations dataset if enabled + if ptes_enable: + ptes_ds = xr.open_dataset(snakemake.input.ptes_operations) + num_ptes_layers = int(ptes_ds.attrs["num_layers"]) + else: + ptes_ds = None + num_ptes_layers = 0 + + heat_sources = expand_heat_sources_with_ptes_layers(heat_sources, num_ptes_layers) + central_heating_forward_temperature: xr.DataArray = xr.open_dataarray( snakemake.input.central_heating_forward_temperature_profiles ) @@ -285,7 +306,6 @@ def get_heat_pump_cooling( heat_source_cooling=get_heat_pump_cooling( heat_source_name=heat_source_key, default_heat_source_cooling=snakemake.params.heat_source_cooling, - snakemake_input=snakemake.input, return_temperature=central_heating_return_temperature, ), ).assign_coords(heat_source=heat_source_key) @@ -293,3 +313,6 @@ def get_heat_pump_cooling( ], dim="heat_source", ).to_netcdf(snakemake.output.heat_source_preheater_utilisation_profiles) + + if ptes_ds is not None: + ptes_ds.close() diff --git a/scripts/build_ptes_operations/ptes_approximator.py b/scripts/build_ptes_operations/ptes_approximator.py new file mode 100644 index 0000000000..6bd20d8d60 --- /dev/null +++ b/scripts/build_ptes_operations/ptes_approximator.py @@ -0,0 +1,182 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT + +import numpy as np +import xarray as xr + + +class PtesApproximator: + def __init__( + self, + forward_temperature: xr.DataArray, + return_temperature: xr.DataArray, + top_temperature: float, + bottom_temperature: float, + num_layers: int, + design_top_temperature: float, + design_bottom_temperature: float, + design_standing_losses: float, + ): + self.forward_temperature = forward_temperature + self.return_temperature = return_temperature + self.top_temperature = top_temperature + self.bottom_temperature = bottom_temperature + self.num_layers = num_layers + self.design_top_temperature = design_top_temperature + self.design_bottom_temperature = design_bottom_temperature + self.design_standing_losses = design_standing_losses + + self._time_dim = [d for d in forward_temperature.dims if d != "name"][0] + + @property + def layer_temperatures(self) -> np.ndarray: + """Calculate the fixed temperature of each layer [°C], hottest first.""" + return np.linspace( + self.top_temperature, self.bottom_temperature, self.num_layers + 1 + )[:-1] + + @property + def volume_weights(self) -> np.ndarray: + """ + Volume weight W_l per layer. + + W_l = (T_top - T_bottom) / (T_l - T_bottom). + A unit of energy in a colder layer occupies more volume. + W_1 = 1 for the hottest layer; W_l ≥ 1 for colder layers. + """ + weights = (self.top_temperature - self.bottom_temperature) / ( + self.layer_temperatures - self.bottom_temperature + ) + return weights + + @property + def charging_availability(self) -> xr.DataArray: + """ + Binary charging availability Φ⁺_{l,t}. + + For each timestep, only the layer whose temperature is closest to the forward temperature can receive charge. Returns a DataArray with dimensions (snapshot, name, layer). + """ + # |T_l - T_fwd_t| for each layer + # Broadcast: T_fwd is (snapshot, name), T_layers is (layer,) + layer_dim = xr.DataArray( + self.layer_temperatures, + dims=["layer"], + coords={"layer": np.arange(len(self.layer_temperatures))}, + ) + abs_diff = np.abs( + layer_dim - self.forward_temperature + ) # (snapshot, name, layer) + + # Find the layer index with minimum distance at each (snapshot, name) + closest_layer = abs_diff.argmin(dim="layer") # (snapshot, name) + + # Build binary availability: 1 where layer == closest_layer + layer_indices = xr.DataArray( + np.arange(len(self.layer_temperatures)), + dims=["layer"], + coords={"layer": np.arange(len(self.layer_temperatures))}, + ) + + binary_availability = (layer_indices == closest_layer).astype(float) + binary_availability = binary_availability.transpose( + self._time_dim, "name", "layer" + ) + + return binary_availability + + @property + def interlayer_heat_transfer_coefficients(self) -> np.ndarray: + """Tbd.""" + return np.full(max(self.num_layers - 1, 0), 0.01) + + @property + def boost_ratios(self) -> xr.DataArray: + """ + Resistive boost ratio α^RH_{l,t} per layer per timestep. + + α^RH_{l,t} = max(0, (T_fwd_t - T_l) / (T_l - T_bottom)) + + Returns DataArray with dimensions (snapshot, name, layer). + """ + layer_dim = xr.DataArray( + self.layer_temperatures, + dims=["layer"], + coords={"layer": np.arange(len(self.layer_temperatures))}, + ) + + numerator = self.forward_temperature - layer_dim # (snapshot, name, layer) + denominator = layer_dim - self.bottom_temperature # (layer,) + + alpha = (numerator / denominator).clip(min=0) + return alpha + + @property + def e_max_pu(self) -> float: + """ + Normalized storage capacity as fraction of design capacity. + + e_max_pu = (T_top - T_bottom) / (T_design_top - T_design_bottom), + clipped to non-negative. + """ + delta_t = self.top_temperature - self.bottom_temperature + max_delta_t = self.design_top_temperature - self.design_bottom_temperature + return max(0, delta_t / max_delta_t) + + @property + def standing_losses(self) -> float: + """Tbd.""" + return np.full(self.num_layers, self.design_standing_losses) + + def to_dataset(self) -> xr.Dataset: + """ + Export all pre-computed parameters as a single xr.Dataset. + + Variables: + - layer_temperatures: (layer,) + - volume_weights: (layer,) + - charging_availability: (snapshot, name, layer) + - interlayer_transfer_coefficients: (layer_pair,) + - boost_ratios: (snapshot, name, layer) + - standing_losses: (layer,) + """ + layer_coord = np.arange(self.num_layers) + pair_coord = np.arange(self.num_layers - 1) + + ds = xr.Dataset( + { + "layer_temperatures": xr.DataArray( + self.layer_temperatures, + dims=["layer"], + coords={"layer": layer_coord}, + ), + "volume_weights": xr.DataArray( + self.volume_weights, + dims=["layer"], + coords={"layer": layer_coord}, + ), + "charging_availability": self.charging_availability, + "interlayer_heat_transfer_coefficients": xr.DataArray( + self.interlayer_heat_transfer_coefficients, + dims=["layer_pair"], + coords={"layer_pair": pair_coord}, + ), + "boost_ratios": self.boost_ratios, + "standing_losses": xr.DataArray( + self.standing_losses, + dims=["layer"], + coords={"layer": layer_coord}, + ), + "e_max_pu": self.e_max_pu, + }, + attrs={ + "num_layers": self.num_layers, + "top_temperature": self.top_temperature, + "bottom_temperature": self.bottom_temperature, + "heat_transfer_coefficient": 0.01, + "design_top_temperature": self.design_top_temperature, + "design_bottom_temperature": self.design_bottom_temperature, + "design_standing_losses": self.design_standing_losses, + }, + ) + return ds diff --git a/scripts/build_ptes_operations/ptes_temperature_approximator.py b/scripts/build_ptes_operations/ptes_temperature_approximator.py deleted file mode 100644 index 61140ef31d..0000000000 --- a/scripts/build_ptes_operations/ptes_temperature_approximator.py +++ /dev/null @@ -1,315 +0,0 @@ -# SPDX-FileCopyrightText: Contributors to PyPSA-Eur -# -# SPDX-License-Identifier: MIT - -import logging - -import xarray as xr - -logger = logging.getLogger(__name__) - - -class PtesTemperatureApproximator: - """ - A unified class to handle pit thermal energy storage (PTES) temperature-related calculations. - - It calculates top temperature profiles, determines when charge or discharge boosting is needed, - and approximates storage capacity based on temperature differences. - - Attributes - ---------- - forward_temperature : xr.DataArray - The forward temperature profile from the district heating network. - return_temperature : xr.DataArray - The return temperature profile from the district heating network. - top_temperature : float | str - Operational temperature specification for top layer in PTES. - bottom_temperature : float | str - Operational temperature specification for bottom layer in PTES. - charge_boosting_required : bool - Whether charge boosting is required/allowed. - discharge_boosting_required : bool - Whether discharge boosting is required/allowed. - temperature_dependent_capacity : bool - Whether storage capacity varies with temperature. If False, assumes constant capacity. - design_top_temperature : float - Maximum design temperature for the top layer of PTES, used for capacity normalization. - design_bottom_temperature : float - Minimum design temperature for the bottom layer of PTES, used for capacity normalization. - """ - - def __init__( - self, - forward_temperature: xr.DataArray, - return_temperature: xr.DataArray, - top_temperature: float | str, - bottom_temperature: float | str, - charge_boosting_required: bool, - discharge_boosting_required: bool, - temperature_dependent_capacity: bool, - design_top_temperature: float, - design_bottom_temperature: float, - ): - """ - Initialize PtesTemperatureApproximator. - - Parameters - ---------- - forward_temperature : xr.DataArray - The forward temperature profile from the district heating network. - return_temperature : xr.DataArray - The return temperature profile from the district heating network. - top_temperature : float | str - Operational temperature of top layer in PTES. Either a float value or 'forward' for dynamic profiles. - bottom_temperature : float | str - Operational temperature of bottom layer in PTES. Either a float value or 'return' for dynamic profiles. - charge_boosting_required : bool - Whether charge boosting is required/allowed. - discharge_boosting_required : bool - Whether discharge boosting is required/allowed. - temperature_dependent_capacity : bool - Whether storage capacity varies with temperature. If False, assumes constant capacity. - design_top_temperature : float - Maximum design temperature for the top layer of PTES, used for capacity normalization - and clipping dynamic top temperature profiles. - design_bottom_temperature : float - Minimum design temperature for the bottom layer of PTES, used for capacity normalization. - """ - self.forward_temperature = forward_temperature - self.return_temperature = return_temperature - self.top_temperature = top_temperature - self.bottom_temperature = bottom_temperature - self.charge_boosting_required = charge_boosting_required - self.discharge_boosting_required = discharge_boosting_required - self.temperature_dependent_capacity = temperature_dependent_capacity - self.design_top_temperature = design_top_temperature - self.design_bottom_temperature = design_bottom_temperature - - if self.charge_boosting_required: - raise NotImplementedError( - "Charge boosting for PTES is currently not supported but might be retintroduced in the future." - ) - - @property - def top_temperature_profile(self) -> xr.DataArray: - """ - PTES top layer temperature profile. - - Returns either the forward temperature (if top_temperature == 'forward'), - clipped to the design_top_temperature, or a constant temperature profile - (if top_temperature is a numeric value). - - Returns - ------- - xr.DataArray - The resulting top temperature profile for PTES. - """ - if self.top_temperature == "forward": - logger.info( - f"PTES top temperature profile: Using dynamic forward temperature from district heating network, clipped by design top temperature to {self.design_top_temperature}°C " - f"Forward temperature range: {float(self.forward_temperature.min().values):.1f}°C to {float(self.forward_temperature.max().values):.1f}°C)" - ) - return self.forward_temperature.where( - self.forward_temperature <= self.design_top_temperature, - self.design_top_temperature, - ) - elif isinstance(self.top_temperature, (int, float)): - logger.info( - f"PTES top temperature profile: Using constant temperature of {self.top_temperature}°C " - f"for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.full_like(self.forward_temperature, self.top_temperature) - else: - raise ValueError( - f"Invalid top_temperature: {self.top_temperature}. " - "Must be 'forward' or a numeric value." - ) - - @property - def bottom_temperature_profile(self) -> xr.DataArray: - """ - PTES bottom layer temperature profile. - - Returns either the return temperature (if bottom_temperature == 'return') - or a constant temperature profile (if bottom_temperature is a numeric value). - - Returns - ------- - xr.DataArray - The resulting bottom temperature profile for PTES. - """ - if self.bottom_temperature == "return": - logger.info( - f"PTES bottom temperature profile: Using dynamic return temperature from district heating network " - f"(shape: {self.return_temperature.shape}, range: {float(self.return_temperature.min().values):.1f}°C to {float(self.return_temperature.max().values):.1f}°C)" - ) - return self.return_temperature - elif isinstance(self.bottom_temperature, (int, float)): - logger.info( - f"PTES bottom temperature profile: Using constant temperature of {self.bottom_temperature}°C " - f"for all {self.return_temperature.size} snapshots and nodes" - ) - return xr.full_like(self.return_temperature, self.bottom_temperature) - else: - raise ValueError( - f"Invalid bottom_temperature: {self.bottom_temperature}. " - "Must be 'return' or a numeric value." - ) - - @property - def e_max_pu(self) -> xr.DataArray: - """ - Calculate e_max_pu for PTES as design_temperature_delta / actual_temperature_delta. - - Returns - ------- - xr.DataArray - Normalized delta T values between 0 and 1, representing the - available storage capacity as a fraction of maximum design capacity. - If temperature_dependent_capacity is False, returns constant capacity of 1.0. - """ - if self.temperature_dependent_capacity: - delta_t = self.top_temperature_profile - self.bottom_temperature_profile - # Get max possible delta_t for normalization - max_delta_t = self.design_top_temperature - self.design_bottom_temperature - normalized_delta_t = delta_t / max_delta_t - result = normalized_delta_t.clip(min=0) # Ensure non-negative values - logger.info( - f"PTES capacity (e_max_pu): Calculating temperature-dependent capacity. " - f"Normalization: max_delta_t={max_delta_t:.2f}K (max_top={self.design_top_temperature:.2f}°C, min_bottom={self.design_bottom_temperature:.2f}°C). " - f"Resulting capacity range: {float(result.min().values):.3f} to {float(result.max().values):.3f} p.u." - ) - return result - else: - logger.info( - f"PTES capacity (e_max_pu): Using constant capacity of 1.0 p.u. (temperature-independent) " - f"for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.ones_like(self.forward_temperature) - - @property - def boost_per_discharge(self) -> xr.DataArray: - """ - Calculate the additional lift required between the store's - current top temperature and the forward temperature with the lift - already achieved inside the store. - - Notes - ----- - The total thermal output required to reach the forward temperature is: - - Q_total = Q_discharge + Q_boost - - The total heat transfer is partitioned into: - - Q_discharge = Ṽ·ρ·cₚ·(T_top − T_bottom) - Q_boost = Ṽ·ρ·cₚ·(T_forward − T_top) - - Solving this to constant Ṽ gives α as the ratio of required boost to available store energy: - - α = Q_boost / Q_discharge - = (T_forward − T_top) / (T_top − T_bottom) - - This expression quantifies the share of PTES output that is covered - by stored energy relative to the additional heating needed to meet - the desired forward temperature. - - Returns - ------- - xr.DataArray - The resulting fraction of PTES charge that must be further heated. - """ - if self.discharge_boosting_required: - result = ( - (self.forward_temperature - self.top_temperature_profile) - / (self.top_temperature_profile - self.bottom_temperature_profile) - ).where(self.forward_temperature > self.top_temperature_profile, 0) - - # Count how many snapshots require boosting - boosting_needed = ( - (self.forward_temperature > self.top_temperature_profile).sum().values - ) - total_snapshots = self.forward_temperature.size - - logger.info( - f"Discharge boosting (boost_per_discharge): Enabled. " - f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations " - f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})" - ) - return result - else: - logger.info( - f"Discharge boosting (boost_per_discharge): Not required. " - f"Returning boost_per_discharge=0 for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.zeros_like(self.forward_temperature) - - @property - def boost_per_charge(self) -> xr.DataArray: - """ - Calculate the additional boost energy ratio required during charging. - - .. note:: - Charge boosting is currently not implemented and will raise - NotImplementedError if charge_boosting_required is True. - - This calculates how much additional energy is needed to raise the PTES - top layer from the forward temperature to the design top temperature, - relative to the energy delivered by charging. - - Notes - ----- - To fill the storage from the return temperature all the way up to its - design top temperature, the total thermal energy required is split into: - - Q_charge = Ṽ·ρ·cₚ·(T_forward − T_bottom) - Q_boost = Ṽ·ρ·cₚ·(T_top − T_forward) - - - Q_charge is the energy delivered by charging to the forward setpoint. - - Q_boost is the extra boost energy needed to reach the design top temperature. - - Defining α as the ratio of boost energy to charge energy: - - α = Q_boost / Q_charge - = (T_top − T_forward) / (T_forward − T_return) - - Wherever the forward temperature meets or exceeds the design top temperature, - α is set to zero since no further boost is needed. - - Returns - ------- - xr.DataArray - The ratio of additional boost energy needed to the energy delivered by charging. - """ - if self.charge_boosting_required: - # Get the max top temperature value - max_top = ( - self.top_temperature - if isinstance(self.top_temperature, (int, float)) - else self.top_temperature_profile - ) - result = ( - ( - (max_top - self.forward_temperature) - / (self.forward_temperature - self.return_temperature) - ) - .where(self.forward_temperature < max_top, 0) - .clip(max=1) - ) - - # Count how many snapshots require boosting - boosting_needed = (self.forward_temperature < max_top).sum().values - total_snapshots = self.forward_temperature.size - - logger.info( - f"Charge boosting (boost_per_charge): Enabled. " - f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations " - f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})" - ) - return result - else: - logger.info( - f"Charge boosting (boost_per_charge): Not required. " - f"Returning boost_per_charge=0 for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.zeros_like(self.forward_temperature) diff --git a/scripts/build_ptes_operations/run.py b/scripts/build_ptes_operations/run.py index 347c5e5efe..d16ddb3565 100644 --- a/scripts/build_ptes_operations/run.py +++ b/scripts/build_ptes_operations/run.py @@ -4,24 +4,12 @@ """ Build PTES (Pit Thermal Energy Storage) operational profiles. -This script calculates temperature and capacity profiles for pit thermal energy -storage systems integrated with district heating networks. It determines: +This script pre-computes temperature-dependent parameters for pit thermal +energy storage systems integrated with district heating networks, using the +unified ``PtesApproximator`` class. -1. **Top temperature profiles**: The operational top layer temperature, either - following the district heating forward temperature (clipped to design limits) - or a constant value. - -2. **Capacity profiles (e_max_pu)**: Normalized storage capacity based on the - temperature difference between top and bottom layers, relative to the design - temperature difference. This captures how storage capacity varies when - operating temperatures deviate from design conditions. - -3. **Discharge boosting profiles**: When the PTES top temperature is below the - required forward temperature, additional heating (boosting) is needed during - discharge. This profile quantifies the ratio of boost energy to stored energy. - -The outputs are used by ``prepare_sector_network.py`` to configure PTES storage -components, charger/discharger links, and optional resistive boosting infrastructure. +Outputs are consumed by ``prepare_sector_network.py`` and by the COP / heat +source utilisation profile builders. Relevant Settings ----------------- @@ -31,31 +19,23 @@ district_heating: ptes: enable: true - temperature_dependent_capacity: false # if true, e_max_pu varies with temperature difference (static but scaled if top/bottom are constant) - charge_boosting_required: false # currently not supported - discharge_resistive_boosting: false # if true, adds resistive boosting link for discharge boosting and disables heat pump boosting - top_temperature: 90 # or "forward" for dynamic - bottom_temperature: 35 # or "return" for dynamic - design_top_temperature: 90 # used to compute design temperature difference for e_max_pu if temperature_dependent_capacity is true - design_bottom_temperature: 35 # used to compute design temperature difference for e_max_pu if temperature_dependent_capacity is true + top_temperature: 90 + bottom_temperature: 35 + design_top_temperature: 90 + design_bottom_temperature: 35 + discharge_resistive_boosting: false + layered: + num_layers: 1 Inputs ------ - ``resources//central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Forward temperature profiles for district heating networks (°C). - ``resources//central_heating_return_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Return temperature profiles for district heating networks (°C). Outputs ------- -- ``resources//temp_ptes_top_profiles_base_s_{clusters}_{planning_horizons}.nc`` - PTES top layer temperature profile (°C), clipped to design_top_temperature. -- ``resources//ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Normalized PTES capacity profiles (0-1 p.u.). -- ``resources//ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Discharge boosting ratio profiles. Values represent the ratio of boost energy - to discharge energy needed when T_forward > T_top. Only used when resistive - boosting is enabled. +- ``resources//ptes_layered_params_base_s_{clusters}_{planning_horizons}.nc`` + Pre-computed parameter dataset for the PTES model. References ---------- @@ -69,9 +49,7 @@ import xarray as xr from scripts._helpers import set_scenario_config -from scripts.build_ptes_operations.ptes_temperature_approximator import ( - PtesTemperatureApproximator, -) +from scripts.build_ptes_operations.ptes_approximator import PtesApproximator logger = logging.getLogger(__name__) @@ -87,43 +65,32 @@ set_scenario_config(snakemake) - logger.info( - "Loading district heating temperature profiles and calculating PTES operational profiles" - ) + logger.info("Building PTES operational profiles") logger.info(f"PTES configuration: {snakemake.params}") - # Discharge boosting profiles are calculated when resistive boosting is enabled. - # For heat pump-based boosting, add "ptes" to central heating heat sources instead and disable resistive boosting. - discharge_boosting_required: bool = snakemake.params.discharge_resistive_boosting - - ptes_temperature_approximator = PtesTemperatureApproximator( - forward_temperature=xr.open_dataarray( - snakemake.input.central_heating_forward_temperature_profiles - ), - return_temperature=xr.open_dataarray( - snakemake.input.central_heating_return_temperature_profiles - ), - top_temperature=snakemake.params.top_temperature, - bottom_temperature=snakemake.params.bottom_temperature, - charge_boosting_required=snakemake.params.charge_boosting_required, - discharge_boosting_required=discharge_boosting_required, - temperature_dependent_capacity=snakemake.params.temperature_dependent_capacity, - design_bottom_temperature=snakemake.params.design_bottom_temperature, - design_top_temperature=snakemake.params.design_top_temperature, + forward_temperature = xr.open_dataarray( + snakemake.input.central_heating_forward_temperature_profiles ) - - ptes_temperature_approximator.top_temperature_profile.to_netcdf( - snakemake.output.ptes_top_temperature_profiles + return_temperature = xr.open_dataarray( + snakemake.input.central_heating_return_temperature_profiles ) - ptes_temperature_approximator.bottom_temperature_profile.to_netcdf( - snakemake.output.ptes_bottom_temperature_profiles - ) + # Resolve top/bottom temperature to constant floats + top_temp = float(snakemake.params.top_temperature) + bottom_temp = float(snakemake.params.bottom_temperature) - ptes_temperature_approximator.e_max_pu.to_netcdf( - snakemake.output.ptes_e_max_pu_profiles - ) + layered_config = snakemake.params.layered - ptes_temperature_approximator.boost_per_discharge.to_netcdf( - snakemake.output.ptes_boost_per_discharge_profiles + approximator = PtesApproximator( + forward_temperature=forward_temperature, + return_temperature=return_temperature, + top_temperature=top_temp, + bottom_temperature=bottom_temp, + num_layers=layered_config["num_layers"], + design_top_temperature=snakemake.params.design_top_temperature, + design_bottom_temperature=snakemake.params.design_bottom_temperature, + design_standing_losses=0.0, # placeholder; actual value set from costs data ) + + # Write single dataset with all pre-computed PTES parameters + approximator.to_dataset().to_netcdf(snakemake.output.ptes_operations) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 6aba369787..a9af0959a9 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -78,6 +78,9 @@ class HeatSource(Enum): AIR = "air" GROUND = "ground" PTES = "ptes" + PTES_LAYER_0 = "ptes layer 0" + PTES_LAYER_1 = "ptes layer 1" + PTES_LAYER_2 = "ptes layer 2" # PTX excess heat sources ELECTROLYSIS_WASTE = "electrolysis_waste" FISCHER_TROPSCH_WASTE = "fischer_tropsch_waste" @@ -111,7 +114,12 @@ def source_type(self) -> HeatSourceType: return HeatSourceType.INEXHAUSTIBLE elif self in [HeatSource.GEOTHERMAL, HeatSource.RIVER_WATER]: return HeatSourceType.SUPPLY_LIMITED - elif self == HeatSource.PTES: + elif self in [ + HeatSource.PTES, + HeatSource.PTES_LAYER_0, + HeatSource.PTES_LAYER_1, + HeatSource.PTES_LAYER_2, + ]: return HeatSourceType.STORAGE else: return HeatSourceType.PROCESS_WASTE @@ -289,7 +297,10 @@ def requires_heat_pump(self, ptes_discharge_resistive_boosting: bool) -> bool: bool False for PTES with resistive boosting, True otherwise. """ - if self == HeatSource.PTES and ptes_discharge_resistive_boosting: + if ( + self.source_type == HeatSourceType.STORAGE + and ptes_discharge_resistive_boosting + ): logging.info( "PTES configured with resistive boosting during discharge; " "heat pump not built for PTES." diff --git a/scripts/definitions/heat_system.py b/scripts/definitions/heat_system.py index b7ad87994b..2f22d2aa13 100644 --- a/scripts/definitions/heat_system.py +++ b/scripts/definitions/heat_system.py @@ -231,15 +231,18 @@ def heat_pump_costs_name(self, heat_source: HeatSource | str) -> str: # Check if this is an excess-heat-sourced heat pump if heat_source in [ HeatSource.PTES, + HeatSource.PTES_LAYER_0, + HeatSource.PTES_LAYER_1, + HeatSource.PTES_LAYER_2, HeatSource.GEOTHERMAL, HeatSource.SEA_WATER, HeatSource.RIVER_WATER, - HeatSource.ELECTROLYSIS_waste, - HeatSource.FISCHER_TROPSCH_waste, - HeatSource.SABATIER_waste, - HeatSource.HABER_BOSCH_waste, - HeatSource.METHANOLISATION_waste, - HeatSource.FUEL_CELL_waste, + HeatSource.ELECTROLYSIS_WASTE, + HeatSource.FISCHER_TROPSCH_WASTE, + HeatSource.SABATIER_WASTE, + HeatSource.HABER_BOSCH_WASTE, + HeatSource.METHANOLISATION_WASTE, + HeatSource.FUEL_CELL_WASTE, ]: return f"{self.central_or_decentral} excess-heat-sourced heat pump" else: diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 1302f10fb3..70af3cdd1e 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -2755,8 +2755,7 @@ def add_heat( heat_source_direct_utilisation_profile_file: str, heat_source_preheater_utilisation_profile_file: str, hourly_heat_demand_total_file: str, - ptes_e_max_pu_file: str, - ptes_boost_per_discharge_profile_file: str, + ptes_operations_file: str, ates_e_nom_max: str, ates_capex_as_fraction_of_geothermal_heat_source: float, ates_recovery_factor: float, @@ -3098,84 +3097,137 @@ def add_heat( ): 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", - ) + ptes_ds = xr.open_dataset(ptes_operations_file) + num_layers = int(ptes_ds.attrs["num_layers"]) energy_to_power_ratio_water_pit = costs.at[ "central water pit storage", "energy to power ratio" ] - n.add( - "Link", - nodes, - suffix=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, - lifetime=costs.at["central water pit storage", "lifetime"], - marginal_cost=costs.at["central water pit charger", "marginal_cost"], - ) + e_max_pu = float(ptes_ds["e_max_pu"]) - n.add( - "Bus", - HeatSource.PTES.resource_bus(nodes, heat_system), - location=nodes, - carrier=f"{heat_system} ptes heat", - unit="MWh_th", - ) + for layer in range(num_layers): + layer_suffix = f" layer {layer}" if num_layers > 1 else "" + n.add( + "Bus", + nodes + f" {heat_system} water pits{layer_suffix}", + location=nodes, + carrier=f"{heat_system} water pits", + unit="MWh_th", + ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} water pits discharger", - bus0=nodes + f" {heat_system} water pits", - bus1=HeatSource.PTES.resource_bus(nodes, heat_system), - carrier=f"{heat_system} water pits discharger", - efficiency=costs.at[ - "central water pit discharger", - "efficiency", - ], - p_nom_extendable=True, - lifetime=costs.at["central water pit storage", "lifetime"], - ) - n.links.loc[ - nodes + f" {heat_system} water pits charger", - "energy to power ratio", - ] = energy_to_power_ratio_water_pit + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits charger{layer_suffix}", + bus0=nodes + f" {heat_system} heat", + bus1=nodes + f" {heat_system} water pits{layer_suffix}", + efficiency=costs.at[ + "central water pit charger", + "efficiency", + ], + carrier=f"{heat_system} water pits charger", + p_nom_extendable=True, + lifetime=costs.at["central water pit storage", "lifetime"], + marginal_cost=costs.at[ + "central water pit charger", "marginal_cost" + ], + p_max_pu=ptes_ds["charging_availability"] + .sel(layer=layer) + .to_pandas(), + ) - if options["district_heating"]["ptes"]["temperature_dependent_capacity"]: - # Load pre-calculated e_max_pu profiles - e_max_pu_data = xr.open_dataarray(ptes_e_max_pu_file) - e_max_pu = ( - e_max_pu_data.sel(name=nodes).to_pandas().reindex(index=n.snapshots) + heat_source = HeatSource(f"ptes{layer_suffix}") + n.add( + "Bus", + heat_source.resource_bus(nodes, heat_system), + location=nodes, + carrier=f"{heat_system} ptes heat", + unit="MWh_th", ) - else: - e_max_pu = 1 - n.add( - "Store", - nodes, - suffix=f" {heat_system} water pits", - bus=nodes + f" {heat_system} water pits", - e_cyclic=True, - e_nom_extendable=True, - e_max_pu=e_max_pu, - carrier=f"{heat_system} water pits", - standing_loss=costs.at["central water pit storage", "standing_losses"] - / 100, # convert %/hour into unit/hour - capital_cost=costs.at["central water pit storage", "capital_cost"], - lifetime=costs.at["central water pit storage", "lifetime"], - ) + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits discharger{layer_suffix}", + bus0=nodes + f" {heat_system} water pits{layer_suffix}", + bus1=heat_source.resource_bus(nodes, heat_system), + carrier=f"{heat_system} water pits discharger", + efficiency=costs.at[ + "central water pit discharger", + "efficiency", + ], + p_nom_extendable=True, + lifetime=costs.at["central water pit storage", "lifetime"], + ) + n.links.loc[ + nodes + f" {heat_system} water pits charger{layer_suffix}", + "energy to power ratio", + ] = energy_to_power_ratio_water_pit + n.add( + "Store", + nodes, + suffix=f" {heat_system} water pits{layer_suffix}", + bus=nodes + f" {heat_system} water pits{layer_suffix}", + # e_cyclic=True, + e_initial=0, + e_nom_extendable=True, + e_max_pu=e_max_pu, + carrier=f"{heat_system} water pits", + standing_loss=costs.at[ + "central water pit storage", "standing_losses" + ] + / 100 + if num_layers == 1 + else 0, + capital_cost=costs.at["central water pit storage", "capital_cost"] + if num_layers == 1 + else 0, + lifetime=costs.at["central water pit storage", "lifetime"], + ) + + # Inter-layer links: bidirectional flow between adjacent layers + for pair_idx in range(num_layers - 1): + upper_layer = f"layer {pair_idx}" + lower_layer = f"layer {pair_idx + 1}" + + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits inter {upper_layer}-{lower_layer}", + bus0=nodes + f" {heat_system} water pits {upper_layer}", + bus1=nodes + f" {heat_system} water pits {lower_layer}", + carrier=f"{heat_system} water pits interlayer", + p_nom_extendable=True, + p_min_pu=-1, + ) + + if num_layers > 1: + # Main bus for the aggregate Store + n.add( + "Bus", + nodes + f" {heat_system} water pits", + location=nodes, + carrier=f"{heat_system} water pits", + unit="MWh_th", + ) + + n.add( + "Store", + nodes, + suffix=f" {heat_system} water pits", + bus=nodes + f" {heat_system} water pits", + e_cyclic=True, + e_nom_extendable=True, + e_max_pu=e_max_pu, + carrier=f"{heat_system} water pits", + standing_loss=costs.at[ + "central water pit storage", "standing_losses" + ] + / 100, + capital_cost=costs.at["central water pit storage", "capital_cost"], + lifetime=costs.at["central water pit storage", "lifetime"], + ) if enable_ates and heat_system == HeatSystem.URBAN_CENTRAL: n.add("Carrier", f"{heat_system} aquifer thermal energy storage") @@ -3238,18 +3290,6 @@ def add_heat( costs_name_heat_pump = heat_system.heat_pump_costs_name(heat_source) - cop_heat_pump = ( - cop.sel( - heat_system=heat_system.system_type.value, - heat_source=heat_source.value, - name=nodes, - ) - .transpose("time", "name") - .to_pandas() - if options["time_dep_hp_cop"] - else costs.loc[[costs_name_heat_pump], ["efficiency"]] - ) - heat_carrier = heat_source.heat_carrier(heat_system) if heat_source.requires_bus: @@ -3376,6 +3416,17 @@ def add_heat( "ptes" ]["discharge_resistive_boosting"] ): + cop_heat_pump = ( + cop.sel( + heat_system=heat_system.system_type.value, + heat_source=heat_source.value, + name=nodes, + ) + .transpose("time", "name") + .to_pandas() + if options["time_dep_hp_cop"] + else costs.loc[[costs_name_heat_pump], ["efficiency"]] + ) n.add( "Link", nodes, @@ -3384,18 +3435,19 @@ def add_heat( bus1=nodes, bus2=heat_source.get_heat_pump_input_bus(nodes, heat_system), carrier=f"{heat_system} {heat_source} heat pump", - efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(), + efficiency=1 / cop_heat_pump.clip(lower=0.001), # .squeeze(), efficiency2=heat_source.get_heat_pump_efficiency2(cop_heat_pump), capital_cost=costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor, - p_min_pu=-(cop_heat_pump > 0).squeeze().astype(float), + p_min_pu=-(cop_heat_pump > 0).astype( + float + ), # .squeeze().astype(float), p_max_pu=0, p_nom_extendable=True, lifetime=costs.at[costs_name_heat_pump, "lifetime"], ) if options["resistive_heaters"]: - ptes_heat_source = HeatSource.PTES key = f"{heat_system.central_or_decentral} resistive heater" if ( @@ -3404,16 +3456,7 @@ def add_heat( and params.sector["district_heating"]["ptes"][ "discharge_resistive_boosting" ] - and ptes_heat_source.value - in params.heat_sources[heat_system.system_type.value] ): - ptes_boost_per_discharge_profiles = ( - xr.open_dataarray(ptes_boost_per_discharge_profile_file) - .sel(name=nodes) - .to_pandas() - .reindex(index=n.snapshots) - ) - n.add( "Bus", nodes, @@ -3422,24 +3465,34 @@ def add_heat( carrier=f"{heat_system} resistive heat", ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} water pits resistive booster", - bus0=nodes + f" {heat_system} heat", - bus1=nodes + f" {heat_system} resistive heat", - bus2=ptes_heat_source.preheater_input_bus(nodes, heat_system), - # eff = 1 - eff2 (energy conservation) - efficiency=ptes_boost_per_discharge_profiles - / (ptes_boost_per_discharge_profiles + 1), - # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating - # (similar to HP balance: p_el x COP = p_source + p_el ) - efficiency2=1 / (ptes_boost_per_discharge_profiles + 1), - p_nom_extendable=True, - p_min_pu=-(ptes_boost_per_discharge_profiles > 0).astype(float), - carrier=f"{heat_system} water pits resistive booster", - p_max_pu=0, - ) + for layer in range(num_layers): + layer_suffix = f" layer {layer}" if num_layers > 1 else "" + ptes_heat_source = HeatSource(f"ptes{layer_suffix}") + layer_boost_per_discharge_profile = ( + ptes_ds["boost_ratios"] + .sel(name=nodes, layer=layer) + .to_pandas() + .reindex(index=n.snapshots) + ) + + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits resistive booster{layer_suffix}", + bus0=nodes + f" {heat_system} heat", + bus1=nodes + f" {heat_system} resistive heat", + bus2=ptes_heat_source.preheater_input_bus(nodes, heat_system), + # eff = 1 - eff2 (energy conservation) + efficiency=layer_boost_per_discharge_profile + / (layer_boost_per_discharge_profile + 1), + # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating + # (similar to HP balance: p_el x COP = p_source + p_el ) + efficiency2=1 / (layer_boost_per_discharge_profile + 1), + p_nom_extendable=True, + p_min_pu=-(layer_boost_per_discharge_profile > 0).astype(float), + carrier=f"{heat_system} water pits resistive booster", + p_max_pu=0, + ) n.add( "Link", @@ -6391,7 +6444,7 @@ def add_import_options( heat_source_direct_utilisation_profile_file=snakemake.input.heat_source_direct_utilisation_profiles, heat_source_preheater_utilisation_profile_file=snakemake.input.heat_source_preheater_utilisation_profiles, hourly_heat_demand_total_file=snakemake.input.hourly_heat_demand_total, - ptes_e_max_pu_file=snakemake.input.ptes_e_max_pu_profiles, + ptes_operations_file=getattr(snakemake.input, "ptes_operations", ""), ates_e_nom_max=snakemake.input.ates_potentials, ates_capex_as_fraction_of_geothermal_heat_source=snakemake.params.sector[ "district_heating" @@ -6403,7 +6456,6 @@ def add_import_options( "recovery_factor" ], enable_ates=snakemake.params.sector["district_heating"]["ates"]["enable"], - ptes_boost_per_discharge_profile_file=snakemake.input.ptes_boost_per_discharge_profiles, district_heat_share_file=snakemake.input.district_heat_share, solar_thermal_total_file=snakemake.input.solar_thermal_total, retro_cost_file=snakemake.input.retro_cost, diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 9bfc69f543..b8d9d0a668 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -864,10 +864,12 @@ def add_TES_energy_to_power_ratio_constraints(n: pypsa.Network) -> None: """ indices_charger_p_nom_extendable = n.links.index[ n.links.index.str.contains("water tanks charger|water pits charger") + & ~n.links.index.str.contains("layer") & n.links.p_nom_extendable ] indices_stores_e_nom_extendable = n.stores.index[ n.stores.index.str.contains("water tanks|water pits") + & ~n.stores.index.str.contains("layer") & n.stores.e_nom_extendable ] @@ -930,12 +932,14 @@ def add_TES_charger_ratio_constraints(n: pypsa.Network) -> None: n.links.index.str.contains( "water tanks charger|water pits charger|aquifer thermal energy storage charger" ) + & ~n.links.index.str.contains("layer") & n.links.p_nom_extendable ] indices_discharger_p_nom_extendable = n.links.index[ n.links.index.str.contains( "water tanks discharger|water pits discharger|aquifer thermal energy storage discharger" ) + & ~n.links.index.str.contains("layer") & n.links.p_nom_extendable ] @@ -970,6 +974,177 @@ def add_TES_charger_ratio_constraints(n: pypsa.Network) -> None: n.model.add_constraints(lhs == 0, name="TES_charger_ratio") +def _has_layered_ptes(n: pypsa.Network) -> bool: + """Return True if the network contains layered PTES stores.""" + return n.stores.index.str.contains("water pits layer").any() + + +def add_layered_ptes_volume_capacity_constraint( + n: pypsa.Network, ptes_ds: xr.Dataset +) -> None: + """ + Enforce volume-weighted layer SOC ≤ aggregate e_nom. + + For each node with layered PTES: + Σ_l W_l · e_{l,t} ≤ ē ∀ t + + where W_l is the volume weight of layer l and ē is the aggregate + store's optimal capacity (e_nom variable). + """ + layer_stores = n.stores.index[ + n.stores.index.str.contains("water pits layer") & n.stores.e_nom_extendable + ] + agg_stores = n.stores.index[ + n.stores.index.str.contains("water pits") + & ~n.stores.index.str.contains("layer") + & n.stores.e_nom_extendable + ] + + if layer_stores.empty or agg_stores.empty: + logger.warning( + "No valid layered or aggregate PTES stores found. " + "Not enforcing volume capacity constraints." + ) + return + + constraints = [] + for agg in agg_stores: + # Find corresponding layer stores (same prefix before "water pits") + prefix = agg.split("water pits")[0] + "water pits layer" + layers = layer_stores[layer_stores.str.startswith(prefix)] + if layers.empty: + continue + + weighted_sum = None + for layer_name in layers: + layer_idx = int(layer_name.split("layer")[-1].strip()) + w = float(ptes_ds["volume_weights"].sel(layer=layer_idx).item()) + layer_e = n.model["Store-e"].loc[:, layer_name] + term = w * layer_e + weighted_sum = term if weighted_sum is None else weighted_sum + term + + agg_e_nom = n.model["Store-e_nom"].loc[agg] + # Σ W_l · e_{l,t} - ē ≤ 0 + constraints.append(weighted_sum - agg_e_nom) + + if not constraints: + return + + merged = linopy.expressions.merge( + constraints, dim="Store-ext" if PYPSA_V1 else "name" + ) + n.model.add_constraints(merged <= 0, name="layered_ptes_volume_capacity") + + +def add_layered_ptes_interlayer_flow_constraint( + n: pypsa.Network, ptes_ds: xr.Dataset +) -> None: + """ + Enforce interlayer heat diffusion proportional to upper-layer SOC. + + For each interlayer link: + p_{l→l+1, t} = -Ψ_{l,l+1} · e_{l, t} ∀ t + + where Ψ is the interlayer transfer coefficient and e_{l,t} is the + state-of-energy of the upper layer store. + """ + interlayer_links = n.links.index[ + n.links.carrier.str.contains("water pits interlayer") + ] + + if interlayer_links.empty: + return + + constraints = [] + for link_name in interlayer_links: + # Extract pair index from link name (e.g., "... inter layer 0-layer 1" → 0) + pair_str = link_name.split("inter ")[-1] + pair_idx = int(pair_str.split("layer ")[1].split("-")[0]) + heat_transfer_coef = float( + ptes_ds["interlayer_heat_transfer_coefficients"] + .sel(layer_pair=pair_idx) + .item() + ) + + upper_bus = n.links.at[link_name, "bus0"] + # Find the store on the upper bus + upper_store = n.stores.index[n.stores.bus == upper_bus] + if upper_store.empty: + continue + upper_store = upper_store[0] + + link_p = n.model["Link-p"].loc[:, link_name] + store_e = n.model["Store-e"].loc[:, upper_store] + # p = Ψ · e + constraints.append(link_p - heat_transfer_coef * store_e) + + if not constraints: + return + + merged = linopy.expressions.merge( + constraints, dim="Link, Store" if PYPSA_V1 else "name" + ) + n.model.add_constraints(merged == 0, name="layered_ptes_interlayer_flow") + + +def add_layered_ptes_aggregate_throughput_constraint( + n: pypsa.Network, ptes_ds: xr.Dataset +) -> None: + """ + Limit total volume-weighted charging + discharging power by aggregate capacity. + + For each node with layered PTES: + R · Σ_l W_l · (p_ch_{l,t} + p_dis_{l,t}) ≤ ē_nom ∀ t + + where R is the energy-to-power ratio, W_l the volume weight, p_ch/p_dis + the layer charger/discharger operational power, and ē_nom the aggregate + store capacity. + """ + layer_chargers = n.links.index[ + n.links.index.str.contains("water pits charger") + & n.links.index.str.contains("layer") + ] + agg_stores = n.stores.index[ + n.stores.index.str.contains("water pits") + & ~n.stores.index.str.contains("layer") + & n.stores.e_nom_extendable + ] + + if layer_chargers.empty or agg_stores.empty: + return + + constraints = [] + for agg in agg_stores: + prefix = agg.split("water pits")[0] + "water pits" + + chargers = layer_chargers[layer_chargers.str.startswith(prefix + " charger")] + dischargers = pd.Index( + [c.replace(" charger ", " discharger ") for c in chargers] + ).intersection(n.links.index) + if chargers.empty or dischargers.empty: + continue + + R = n.links.at[chargers[0], "energy to power ratio"] + + weighted_sum = None + for ch, dis in zip(chargers, dischargers): + layer_idx = int(ch.rsplit("layer ", 1)[-1]) + w = float(ptes_ds["volume_weights"].sel(layer=layer_idx).item()) + term = w * (n.model["Link-p"].loc[:, ch] + n.model["Link-p"].loc[:, dis]) + weighted_sum = term if weighted_sum is None else weighted_sum + term + + agg_e_nom = n.model["Store-e_nom"].loc[agg] + constraints.append(R * weighted_sum - agg_e_nom) + + if not constraints: + return + + merged = linopy.expressions.merge( + constraints, dim="Store-ext" if PYPSA_V1 else "name" + ) + n.model.add_constraints(merged <= 0, name="layered_ptes_aggregate_throughput") + + def add_battery_constraints(n): """ Add constraint ensuring that charger = discharger, i.e. @@ -1216,6 +1391,13 @@ def extra_functionality( add_TES_energy_to_power_ratio_constraints(n) add_TES_charger_ratio_constraints(n) + if _has_layered_ptes(n): + ptes_ds = xr.open_dataset(snakemake.input.ptes_operations) + add_layered_ptes_volume_capacity_constraint(n, ptes_ds) + add_layered_ptes_interlayer_flow_constraint(n, ptes_ds) + # add_layered_ptes_aggregate_throughput_constraint(n, ptes_ds) + ptes_ds.close() + add_battery_constraints(n) add_lossy_bidirectional_link_constraints(n) add_pipe_retrofit_constraint(n) From b6727d073f0309602b75b41111c990e716a7b024 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 12 Mar 2026 17:43:24 +0100 Subject: [PATCH 02/17] refactor: remove expand_ptes_layers function and related validation for heat sources --- scripts/build_cop_profiles/run.py | 13 ---------- .../build_heat_source_utilisation_profiles.py | 25 ------------------- 2 files changed, 38 deletions(-) diff --git a/scripts/build_cop_profiles/run.py b/scripts/build_cop_profiles/run.py index b9a1e6f316..c16581138d 100644 --- a/scripts/build_cop_profiles/run.py +++ b/scripts/build_cop_profiles/run.py @@ -283,21 +283,8 @@ def get_sink_inlet_temperature( ptes_ds = None num_ptes_layers = 0 - def expand_ptes_layers(heat_source_names): - """Expand 'ptes' into per-layer entries if multi-layer PTES is enabled.""" - expanded = [] - for name in heat_source_names: - if name == "ptes" and num_ptes_layers > 1: - expanded.extend(f"ptes layer {l}" for l in range(num_ptes_layers)) - else: - expanded.append(name) - return expanded - cop_all_system_types = [] for heat_system_type, heat_sources in snakemake.params.heat_sources.items(): - heat_sources = ( - expand_ptes_layers(heat_sources) if heat_sources else heat_sources - ) cop_this_system_type = [] if not heat_sources: cop_all_system_types.append( diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py index 4667a3a8ef..3c92c2833c 100644 --- a/scripts/build_heat_source_utilisation_profiles.py +++ b/scripts/build_heat_source_utilisation_profiles.py @@ -66,19 +66,6 @@ logger = logging.getLogger(__name__) -def expand_heat_sources_with_ptes_layers( - heat_source_names: list[str], num_layers: int -) -> list[str]: - """Expand 'ptes' into per-layer entries if multi-layer PTES is enabled.""" - expanded = [] - for name in heat_source_names: - if name == "ptes" and num_layers > 1: - expanded.extend(f"ptes layer {l}" for l in range(num_layers)) - else: - expanded.append(name) - return expanded - - def get_source_temperature( snakemake_params: dict, snakemake_input: dict, heat_source_name: str ) -> float | xr.DataArray: @@ -251,16 +238,6 @@ def get_heat_pump_cooling( heat_sources: list[str] = snakemake.params.heat_sources ptes_enable: bool = snakemake.params.ptes_enable - # Validate PTES configuration - has_ptes = any(s == "ptes" or s.startswith("ptes layer") for s in heat_sources) - if ptes_enable and not has_ptes: - raise ValueError( - "PTES is enabled (district_heating.ptes.enable=true) but 'ptes' " - "is not in heat_sources.urban_central. PTES requires being listed " - "in heat_sources to create the necessary buses and links for heat " - "discharge to the 'urban central heat' bus." - ) - # Load PTES operations dataset if enabled if ptes_enable: ptes_ds = xr.open_dataset(snakemake.input.ptes_operations) @@ -269,8 +246,6 @@ def get_heat_pump_cooling( ptes_ds = None num_ptes_layers = 0 - heat_sources = expand_heat_sources_with_ptes_layers(heat_sources, num_ptes_layers) - central_heating_forward_temperature: xr.DataArray = xr.open_dataarray( snakemake.input.central_heating_forward_temperature_profiles ) From 9bc232e4a19eeceb91d2c2fa5eaec276f0b49dc7 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 12 Mar 2026 17:44:06 +0100 Subject: [PATCH 03/17] feat: add layered PTES heat pump capacity constraint function --- scripts/solve_network.py | 48 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index b8d9d0a668..e8bbfcb533 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1036,6 +1036,54 @@ def add_layered_ptes_volume_capacity_constraint( n.model.add_constraints(merged <= 0, name="layered_ptes_volume_capacity") +def add_layered_ptes_heat_pump_capacity_constraint( + n: pypsa.Network, +) -> None: + """ + Tbd + """ + layer_heat_pumps = n.links.index[ + n.links.index.str.contains("ptes layer") & n.links.p_nom_extendable + ] + agg_heat_pumps = n.links.index[ + n.links.index.str.contains("ptes") + & ~n.links.index.str.contains("layer") + & n.links.p_nom_extendable + ] + + if layer_heat_pumps.empty or agg_heat_pumps.empty: + logger.warning( + "No valid layered or aggregate PTES heat pumps found. " + "Not enforcing volume capacity constraints." + ) + return + + constraints = [] + for agg in agg_heat_pumps: + # Find corresponding layer stores (same prefix before "ptes") + prefix = agg.split("ptes")[0] + "ptes layer" + layers = layer_heat_pumps[layer_heat_pumps.str.startswith(prefix)] + if layers.empty: + continue + + layer_p_sum = None + for layer_name in layers: + layer_p = n.model["Link-p"].loc[:, layer_name] + layer_p_sum = layer_p if layer_p_sum is None else layer_p_sum + layer_p + + agg_p_nom = n.model["Link-p_nom"].loc[agg] + # Σ W_l · e_{l,t} - ē ≤ 0 + constraints.append(layer_p_sum - agg_p_nom) + + if not constraints: + return + + merged = linopy.expressions.merge( + constraints, dim="Link-ext" if PYPSA_V1 else "name" + ) + n.model.add_constraints(merged <= 0, name="layered_ptes_heat_pump_capacity") + + def add_layered_ptes_interlayer_flow_constraint( n: pypsa.Network, ptes_ds: xr.Dataset ) -> None: From 1d2d5d8ba20b6f11c45298b81a81f29d484b66ee Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 12 Mar 2026 17:44:50 +0100 Subject: [PATCH 04/17] introduce HP boosting and cost ptes discharging --- scripts/prepare_sector_network.py | 37 +++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 70af3cdd1e..ed80d9cd63 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3157,6 +3157,9 @@ def add_heat( "central water pit discharger", "efficiency", ], + marginal_cost=costs.at[ + "central water pit discharger", "marginal_cost" + ], p_nom_extendable=True, lifetime=costs.at["central water pit storage", "lifetime"], ) @@ -3177,9 +3180,7 @@ def add_heat( standing_loss=costs.at[ "central water pit storage", "standing_losses" ] - / 100 - if num_layers == 1 - else 0, + / 100, capital_cost=costs.at["central water pit storage", "capital_cost"] if num_layers == 1 else 0, @@ -3427,6 +3428,27 @@ def add_heat( if options["time_dep_hp_cop"] else costs.loc[[costs_name_heat_pump], ["efficiency"]] ) + + p_min_pu = ( + 0 + if heat_source == HeatSource.PTES + and params.sector["district_heating"]["ptes"]["layered"][ + "num_layers" + ] + > 1 + else -(cop_heat_pump > 0).astype(float) + ) + capital_cost = ( + 0 + if heat_source + in [ + HeatSource.PTES_LAYER_0, + HeatSource.PTES_LAYER_1, + HeatSource.PTES_LAYER_2, + ] + else costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor + ) + n.add( "Link", nodes, @@ -3435,13 +3457,10 @@ def add_heat( bus1=nodes, bus2=heat_source.get_heat_pump_input_bus(nodes, heat_system), carrier=f"{heat_system} {heat_source} heat pump", - efficiency=1 / cop_heat_pump.clip(lower=0.001), # .squeeze(), + efficiency=1 / cop_heat_pump.clip(lower=0.001), efficiency2=heat_source.get_heat_pump_efficiency2(cop_heat_pump), - capital_cost=costs.at[costs_name_heat_pump, "capital_cost"] - * overdim_factor, - p_min_pu=-(cop_heat_pump > 0).astype( - float - ), # .squeeze().astype(float), + capital_cost=capital_cost, + p_min_pu=p_min_pu, p_max_pu=0, p_nom_extendable=True, lifetime=costs.at[costs_name_heat_pump, "lifetime"], From 30fc5227bbd88375223444c234edf43b0bb25615 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 12 Mar 2026 17:45:02 +0100 Subject: [PATCH 05/17] fix volume capacity weights --- scripts/build_ptes_operations/ptes_approximator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/build_ptes_operations/ptes_approximator.py b/scripts/build_ptes_operations/ptes_approximator.py index 6bd20d8d60..eee5745486 100644 --- a/scripts/build_ptes_operations/ptes_approximator.py +++ b/scripts/build_ptes_operations/ptes_approximator.py @@ -45,7 +45,7 @@ def volume_weights(self) -> np.ndarray: A unit of energy in a colder layer occupies more volume. W_1 = 1 for the hottest layer; W_l ≥ 1 for colder layers. """ - weights = (self.top_temperature - self.bottom_temperature) / ( + weights = (self.design_top_temperature - self.design_bottom_temperature) / ( self.layer_temperatures - self.bottom_temperature ) return weights From bf30c43f2e644b1324a38fb85d916ea26f5dfb1c Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 26 Mar 2026 10:20:52 +0100 Subject: [PATCH 06/17] use conducive heat transfer for interlayer transfer --- .../ptes_approximator.py | 41 +++++++++++++++++-- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/scripts/build_ptes_operations/ptes_approximator.py b/scripts/build_ptes_operations/ptes_approximator.py index eee5745486..05a7a43e47 100644 --- a/scripts/build_ptes_operations/ptes_approximator.py +++ b/scripts/build_ptes_operations/ptes_approximator.py @@ -17,6 +17,7 @@ def __init__( design_top_temperature: float, design_bottom_temperature: float, design_standing_losses: float, + interlayer_heat_transfer_coefficient: float, ): self.forward_temperature = forward_temperature self.return_temperature = return_temperature @@ -26,6 +27,7 @@ def __init__( self.design_top_temperature = design_top_temperature self.design_bottom_temperature = design_bottom_temperature self.design_standing_losses = design_standing_losses + self.interlayer_heat_transfer_coefficient = interlayer_heat_transfer_coefficient self._time_dim = [d for d in forward_temperature.dims if d != "name"][0] @@ -87,8 +89,41 @@ def charging_availability(self) -> xr.DataArray: @property def interlayer_heat_transfer_coefficients(self) -> np.ndarray: - """Tbd.""" - return np.full(max(self.num_layers - 1, 0), 0.01) + """.""" + return np.full(self.num_layers - 1, self.interlayer_heat_transfer_coefficient) + + def _interlayer_heat_transfer_coefficients( + self, + conductivity=0.6, + density=1000, + heat_capacity=4184, + storage_height=15, + thermocline=1, + seconds_per_hour=3600, + ) -> np.ndarray: + """ + Interlayer heat transfer coefficient relative to layer state-of-charge. Requires "nominal" layer height, which is assumed to be storage_height distributed equally across layers. Returns array of length num_layers - 1, where each entry corresponds to the coefficient between layer l and l+1. + + Args: + ---- + conductivity: Thermal conductivity of the storage medium [W/m/K] + density: Density of the storage medium [kg/m3] + heat_capacity: Specific heat capacity of the storage medium [J/kg/K] + storage_height: Total height of the storage medium [m] + seconds_per_hour: Number of seconds in an hour (for unit conversion) + """ + layer_height = storage_height / ( + self.num_layers + 1 + ) # assuming equal layer heights + + return ( + conductivity + / (density * heat_capacity * layer_height) + * seconds_per_hour + * (self.layer_temperatures[:-1] - self.layer_temperatures[1:]) + / (self.layer_temperatures[:-1] - self.bottom_temperature) + / thermocline + ) @property def boost_ratios(self) -> xr.DataArray: @@ -173,7 +208,7 @@ def to_dataset(self) -> xr.Dataset: "num_layers": self.num_layers, "top_temperature": self.top_temperature, "bottom_temperature": self.bottom_temperature, - "heat_transfer_coefficient": 0.01, + "storage_height": 15, # m, assumed for interlayer coefficient calculation "design_top_temperature": self.design_top_temperature, "design_bottom_temperature": self.design_bottom_temperature, "design_standing_losses": self.design_standing_losses, From 12ce6416499f4a948abe6d67f0521cc5458b4368 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 26 Mar 2026 10:32:09 +0100 Subject: [PATCH 07/17] fix: update marginal cost for central water pit charger --- data/custom_costs.csv | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/data/custom_costs.csv b/data/custom_costs.csv index 312cd0b3da..d7b71e3bab 100644 --- a/data/custom_costs.csv +++ b/data/custom_costs.csv @@ -10,4 +10,5 @@ all,battery,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degene all,battery inverter,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, all,home battery storage,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, all,water tank charger,marginal_cost,0.03,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, -all,central water pit charger,marginal_cost,0.025,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, +all,central water pit charger,marginal_cost,0.04,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, +all,central water pit discharger,marginal_cost,0.05,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation, From d6c377a2ec92bdf91af6e85724c11d0b0c823e8a Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Thu, 26 Mar 2026 10:32:25 +0100 Subject: [PATCH 08/17] feat: add additional PTES layers to heat source enumeration --- scripts/definitions/heat_source.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index a9af0959a9..a9ff972cf6 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -81,6 +81,8 @@ class HeatSource(Enum): PTES_LAYER_0 = "ptes layer 0" PTES_LAYER_1 = "ptes layer 1" PTES_LAYER_2 = "ptes layer 2" + PTES_LAYER_3 = "ptes layer 3" + PTES_LAYER_4 = "ptes layer 4" # PTX excess heat sources ELECTROLYSIS_WASTE = "electrolysis_waste" FISCHER_TROPSCH_WASTE = "fischer_tropsch_waste" @@ -119,6 +121,8 @@ def source_type(self) -> HeatSourceType: HeatSource.PTES_LAYER_0, HeatSource.PTES_LAYER_1, HeatSource.PTES_LAYER_2, + HeatSource.PTES_LAYER_3, + HeatSource.PTES_LAYER_4, ]: return HeatSourceType.STORAGE else: From 1cf12d95ebc0a6e4494a1eff067e525eafa9ef67 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 1 Apr 2026 17:05:24 +0200 Subject: [PATCH 09/17] first draft: simplify boosting --- rules/build_sector.smk | 41 +--- .../build_heat_source_utilisation_profiles.py | 215 ++++++++---------- .../ptes_temperature_approximator.py | 148 +----------- scripts/build_ptes_operations/run.py | 22 +- scripts/definitions/heat_source.py | 133 ++--------- scripts/lib/validation/config/sector.py | 4 +- scripts/prepare_sector_network.py | 96 +++----- 7 files changed, 150 insertions(+), 509 deletions(-) diff --git a/rules/build_sector.smk b/rules/build_sector.smk index b446001b78..2490186fd7 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -725,12 +725,6 @@ rule build_ptes_operations: "bottom_temperature", ), snapshots=config_provider("snapshots"), - charge_boosting_required=config_provider( - "sector", "district_heating", "ptes", "charge_boosting_required" - ), - discharge_resistive_boosting=config_provider( - "sector", "district_heating", "ptes", "discharge_resistive_boosting" - ), temperature_dependent_capacity=config_provider( "sector", "district_heating", "ptes", "temperature_dependent_capacity" ), @@ -764,9 +758,6 @@ rule build_ptes_operations: ptes_e_max_pu_profiles=resources( "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc" ), - ptes_boost_per_discharge_profiles=resources( - "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc" - ), resources: mem_mb=2000, log: @@ -777,9 +768,9 @@ rule build_ptes_operations: scripts("build_ptes_operations/run.py") -rule build_direct_heat_source_utilisation_profiles: +rule build_heat_source_utilisation_profiles: message: - "Building direct heat source utilization profiles for industrial applications for {wildcards.clusters} clusters and {wildcards.planning_horizons} planning horizon" + "Building heat source alpha profiles for {wildcards.clusters} clusters and {wildcards.planning_horizons} planning horizon" params: heat_sources=config_provider("sector", "heat_sources", "urban central"), heat_source_cooling=config_provider( @@ -803,11 +794,8 @@ rule build_direct_heat_source_utilisation_profiles: "central_heating_return_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc" ), output: - heat_source_direct_utilisation_profiles=resources( - "heat_source_direct_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc" - ), - heat_source_preheater_utilisation_profiles=resources( - "heat_source_preheater_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc" + heat_source_boosting_profiles=resources( + "heat_source_boosting_profiles_base_s_{clusters}_{planning_horizons}.nc" ), resources: mem_mb=20000, @@ -1749,26 +1737,9 @@ rule prepare_sector_network: )(w) else [] ), - ptes_boost_per_discharge_profiles=lambda w: ( - resources( - "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc" - ) - if config_provider("sector", "district_heating", "ptes", "enable")(w) - and config_provider( - "sector", "district_heating", "ptes", "discharge_resistive_boosting" - )(w) - else [] - ), - heat_source_direct_utilisation_profiles=lambda w: ( - resources( - "heat_source_direct_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc" - ) - if len(config_provider("sector", "heat_sources", "urban central")(w)) > 0 - else [] - ), - heat_source_preheater_utilisation_profiles=lambda w: ( + heat_source_boosting_profiles=lambda w: ( resources( - "heat_source_preheater_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc" + "heat_source_boosting_profiles_base_s_{clusters}_{planning_horizons}.nc" ) if len(config_provider("sector", "heat_sources", "urban central")(w)) > 0 else [] diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py index c37ab62692..18a4e0399b 100644 --- a/scripts/build_heat_source_utilisation_profiles.py +++ b/scripts/build_heat_source_utilisation_profiles.py @@ -2,27 +2,29 @@ # # SPDX-License-Identifier: MIT """ -Build heat source utilisation profiles for district heating networks. +Build heat source alpha profiles for district heating networks. -This script calculates when and how much heat from various sources (geothermal, -PTES, river water, etc.) can be used, based on the temperature relationship -between the heat source and the district heating network. +This script calculates the boosting ratio (alpha) for each heat source: how +much heat pump output is needed per unit of source heat to reach the forward +temperature of the district heating network. -Two utilisation modes are calculated: +**Alpha profile**: For each heat source and timestep, alpha is: -1. **Direct utilisation**: When the source temperature meets or exceeds the - forward temperature (T_source ≥ T_forward), the heat source can directly - supply the district heating network. Profile value is 1.0 (full utilisation) - or 0.0 (not possible). +- 0 when T_source ≥ T_forward (direct use possible, no HP boost needed) +- (T_forward − T_source) / delta_T otherwise, where delta_T is: + - heat_pump_cooling if T_source < T_return + - T_source − T_return + heat_pump_cooling if T_return ≤ T_source < T_forward -2. **Preheater utilisation**: When the source temperature is between the return - and forward temperatures (T_return < T_source < T_forward), the heat source - can preheat the return flow before a heat pump lifts it to forward temperature. - The profile value represents that share of the heat above the return temperature which is utilised to increase the heat pump's sink inflow temperature. The return flow serves as the source inlet. +heat_pump_cooling is the additional temperature drop achievable by the heat +pump beyond the return temperature (i.e. how far below T_return the HP can +extract from the source). For PTES it equals T_return − T_bottom; for other +sources it is the constant config value ``heat_source_cooling``. -These profiles are used by ``prepare_sector_network.py`` to configure heat -utilisation links that model cascading temperature use: direct supply when -possible, preheating when beneficial, with heat pumps handling the final lift. +The alpha profile is consumed by ``prepare_sector_network.py`` to set the +efficiencies of the heat source utilisation links: + +- bus0 (source) → bus1 (DH heat) at efficiency 1 + alpha +- bus2 (HP output bus) drawn at efficiency −alpha per unit of source Relevant Settings ----------------- @@ -48,12 +50,9 @@ Outputs ------- -- ``resources//heat_source_direct_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Direct utilisation profiles indexed by (time, name, heat_source). - Values: 1.0 when T_source ≥ T_forward, 0.0 otherwise. -- ``resources//heat_source_preheater_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Preheater utilisation profiles indexed by (time, name, heat_source). - Values: heat extraction efficiency when T_return < T_source < T_forward, 0.0 otherwise. +- ``resources//heat_source_boosting_profiles_base_s_{clusters}_{planning_horizons}.nc`` + Boosting ratio profiles indexed by (time, name, heat_source). + Values: 0 when T_source ≥ T_forward; (T_fwd − T_src) / delta_T otherwise. """ import logging @@ -109,81 +108,6 @@ def get_source_temperature( return xr.open_dataarray(snakemake_input[f"temp_{heat_source_name}"]) -def get_direct_utilisation_profile( - source_temperature: float | xr.DataArray, forward_temperature: xr.DataArray -) -> xr.DataArray | float: - """ - Calculate when a heat source can directly supply district heating. - - Direct utilisation is possible when the source temperature meets or exceeds - the required forward temperature of the district heating network. - - Parameters - ---------- - source_temperature : float | xr.DataArray - Heat source temperature in °C. If float, applies uniformly. - If DataArray, indexed by (time, name). - forward_temperature : xr.DataArray - District heating forward temperature profiles in °C, - indexed by (time, name). - - Returns - ------- - xr.DataArray - Binary profile: 1.0 where T_source ≥ T_forward (direct use possible), - 0.0 otherwise. - """ - return xr.where(source_temperature >= forward_temperature, 1.0, 0.0) - - -def get_preheater_utilisation_profile( - source_temperature: float | xr.DataArray, - forward_temperature: xr.DataArray, - return_temperature: xr.DataArray, - heat_source_cooling: float, -) -> xr.DataArray | float: - """ - Calculate preheater utilisation efficiency for intermediate-temperature sources. - - When a heat source temperature is between the return and forward temperatures, - it can preheat the return flow before a heat pump provides the final temperature - lift. This improves overall efficiency by reducing the heat pump's lift. - - The efficiency represents the fraction of heat extracted from the source that - goes into preheating (vs. the additional cooling through the heat pump): - - efficiency = (T_source - T_return) / (T_source - T_return + heat_source_cooling) - - Parameters - ---------- - source_temperature : float | xr.DataArray - Heat source temperature in °C. If float, applies uniformly. - If DataArray, indexed by (time, name). - forward_temperature : xr.DataArray - District heating forward temperature profiles in °C, - indexed by (time, name). - return_temperature : xr.DataArray - District heating return temperature profiles in °C, - indexed by (time, name). - heat_source_cooling : float | xr.DataArray - Additional temperature drop (K) when extracting heat from the source - through the heat pump, beyond the preheating contribution. - - Returns - ------- - xr.DataArray - Preheater efficiency profile: value in (0, 1) where T_return < T_source < T_forward, - 0.0 otherwise (source too cold or hot enough for direct use). - """ - return xr.where( - (source_temperature < forward_temperature) - * (source_temperature > return_temperature), - (source_temperature - return_temperature) - / (source_temperature - return_temperature + heat_source_cooling), - 0.0, - ) - - def get_heat_pump_cooling( heat_source_name: str, default_heat_source_cooling: float, @@ -191,11 +115,14 @@ def get_heat_pump_cooling( return_temperature: xr.DataArray = None, ) -> float | xr.DataArray: """ - Get the additional heat source cooling (temperature drop) through the heat pump for a heat source. + Get the additional heat source cooling (temperature drop) through the heat pump. + + This is the number of Kelvin below the return temperature that the heat pump + can extract from the source, i.e. the HP cools the source from T_source down + to T_return − heat_pump_cooling. - For PTES, this equals the temperature difference between - return flow and bottom layers (return_temperature - bottom_temperature), which can - be time-varying. For other sources, uses the default constant value. + For PTES, equals T_return − T_bottom (the bottom layer is the HP's cold side). + For other sources, uses the constant config value ``heat_source_cooling``. Parameters ---------- @@ -204,20 +131,19 @@ def get_heat_pump_cooling( default_heat_source_cooling : float Default heat source cooling in Kelvin, from config. snakemake_input : dict - Snakemake input files, may contain PTES temperature profiles. + Snakemake input files, may contain PTES bottom temperature profiles. return_temperature : xr.DataArray, optional District heating return temperature profiles in °C. Required for PTES. Returns ------- float | xr.DataArray - Heat source cooling in Kelvin. Returns a float for most sources, - or a DataArray for PTES when temperatures vary with time. + Heat source cooling in Kelvin. Raises ------ ValueError - If heat source is PTES but bottom temperature profile is not provided. + If heat source is PTES but required profiles are not provided. """ if heat_source_name == "ptes": if "temp_ptes_bottom" not in snakemake_input.keys(): @@ -234,12 +160,68 @@ def get_heat_pump_cooling( return default_heat_source_cooling +def get_boosting_profile( + source_temperature: float | xr.DataArray, + forward_temperature: xr.DataArray, + return_temperature: xr.DataArray, + heat_pump_cooling: float | xr.DataArray, +) -> xr.DataArray: + """ + Calculate the boosting ratio: HP heat needed per unit of source heat. + + Alpha represents the ratio of heat pump output required to source heat input + in order to reach the district heating forward temperature: + + alpha = 0 if T_source ≥ T_forward + alpha = (T_fwd − T_src) / delta_T otherwise + + where delta_T is: + + delta_T = heat_pump_cooling if T_source < T_return + delta_T = T_source − T_return + heat_pump_cooling if T_return ≤ T_source < T_forward + + For PTES (where heat_pump_cooling = T_return − T_bottom): + delta_T = T_source − T_bottom (when T_source ≥ T_return, which is the normal case) + + Parameters + ---------- + source_temperature : float | xr.DataArray + Heat source temperature in °C. + forward_temperature : xr.DataArray + District heating forward temperature in °C, indexed by (time, name). + return_temperature : xr.DataArray + District heating return temperature in °C, indexed by (time, name). + heat_pump_cooling : float | xr.DataArray + Additional temperature drop achievable by HP beyond return temperature (K). + + Returns + ------- + xr.DataArray + Alpha profile: 0 where direct use is possible, positive ratio otherwise. + Shape matches forward_temperature. + """ + delta_t = xr.where( + source_temperature < return_temperature, + heat_pump_cooling, + source_temperature - return_temperature + heat_pump_cooling, + ) + # Prevent division by zero (e.g. T_source = T_bottom with resistive boosting). + # In that case the source cannot be economically used: alpha → ∞ is approximated + # by a very large number so the optimiser will not dispatch this configuration. + safe_delta_t = delta_t.where(delta_t > 1e-6, other=1e6) + return xr.where( + source_temperature >= forward_temperature, + 0.0, + (forward_temperature - source_temperature) / safe_delta_t, + ) + + if __name__ == "__main__": if "snakemake" not in globals(): from scripts._helpers import mock_snakemake snakemake = mock_snakemake( - "build_cop_profiles", + "build_heat_source_utilisation_profiles", clusters=48, ) configure_logging(snakemake) @@ -264,30 +246,15 @@ def get_heat_pump_cooling( xr.concat( [ - get_direct_utilisation_profile( + get_boosting_profile( source_temperature=get_source_temperature( snakemake_params=snakemake.params, snakemake_input=snakemake.input, heat_source_name=heat_source_key, ), forward_temperature=central_heating_forward_temperature, - ).assign_coords(heat_source=heat_source_key) - for heat_source_key in heat_sources - ], - dim="heat_source", - ).to_netcdf(snakemake.output.heat_source_direct_utilisation_profiles) - - xr.concat( - [ - get_preheater_utilisation_profile( - source_temperature=get_source_temperature( - heat_source_name=heat_source_key, - snakemake_params=snakemake.params, - snakemake_input=snakemake.input, - ), - forward_temperature=central_heating_forward_temperature, return_temperature=central_heating_return_temperature, - heat_source_cooling=get_heat_pump_cooling( + heat_pump_cooling=get_heat_pump_cooling( heat_source_name=heat_source_key, default_heat_source_cooling=snakemake.params.heat_source_cooling, snakemake_input=snakemake.input, @@ -297,4 +264,4 @@ def get_heat_pump_cooling( for heat_source_key in heat_sources ], dim="heat_source", - ).to_netcdf(snakemake.output.heat_source_preheater_utilisation_profiles) + ).to_netcdf(snakemake.output.heat_source_boosting_profiles) diff --git a/scripts/build_ptes_operations/ptes_temperature_approximator.py b/scripts/build_ptes_operations/ptes_temperature_approximator.py index 61140ef31d..6bd5060c42 100644 --- a/scripts/build_ptes_operations/ptes_temperature_approximator.py +++ b/scripts/build_ptes_operations/ptes_temperature_approximator.py @@ -13,8 +13,8 @@ class PtesTemperatureApproximator: """ A unified class to handle pit thermal energy storage (PTES) temperature-related calculations. - It calculates top temperature profiles, determines when charge or discharge boosting is needed, - and approximates storage capacity based on temperature differences. + It calculates top and bottom temperature profiles and approximates storage capacity based on + temperature differences. Attributes ---------- @@ -26,10 +26,6 @@ class PtesTemperatureApproximator: Operational temperature specification for top layer in PTES. bottom_temperature : float | str Operational temperature specification for bottom layer in PTES. - charge_boosting_required : bool - Whether charge boosting is required/allowed. - discharge_boosting_required : bool - Whether discharge boosting is required/allowed. temperature_dependent_capacity : bool Whether storage capacity varies with temperature. If False, assumes constant capacity. design_top_temperature : float @@ -44,8 +40,6 @@ def __init__( return_temperature: xr.DataArray, top_temperature: float | str, bottom_temperature: float | str, - charge_boosting_required: bool, - discharge_boosting_required: bool, temperature_dependent_capacity: bool, design_top_temperature: float, design_bottom_temperature: float, @@ -63,10 +57,6 @@ def __init__( Operational temperature of top layer in PTES. Either a float value or 'forward' for dynamic profiles. bottom_temperature : float | str Operational temperature of bottom layer in PTES. Either a float value or 'return' for dynamic profiles. - charge_boosting_required : bool - Whether charge boosting is required/allowed. - discharge_boosting_required : bool - Whether discharge boosting is required/allowed. temperature_dependent_capacity : bool Whether storage capacity varies with temperature. If False, assumes constant capacity. design_top_temperature : float @@ -79,17 +69,10 @@ def __init__( self.return_temperature = return_temperature self.top_temperature = top_temperature self.bottom_temperature = bottom_temperature - self.charge_boosting_required = charge_boosting_required - self.discharge_boosting_required = discharge_boosting_required self.temperature_dependent_capacity = temperature_dependent_capacity self.design_top_temperature = design_top_temperature self.design_bottom_temperature = design_bottom_temperature - if self.charge_boosting_required: - raise NotImplementedError( - "Charge boosting for PTES is currently not supported but might be retintroduced in the future." - ) - @property def top_temperature_profile(self) -> xr.DataArray: """ @@ -186,130 +169,3 @@ def e_max_pu(self) -> xr.DataArray: f"for all {self.forward_temperature.size} snapshots and nodes" ) return xr.ones_like(self.forward_temperature) - - @property - def boost_per_discharge(self) -> xr.DataArray: - """ - Calculate the additional lift required between the store's - current top temperature and the forward temperature with the lift - already achieved inside the store. - - Notes - ----- - The total thermal output required to reach the forward temperature is: - - Q_total = Q_discharge + Q_boost - - The total heat transfer is partitioned into: - - Q_discharge = Ṽ·ρ·cₚ·(T_top − T_bottom) - Q_boost = Ṽ·ρ·cₚ·(T_forward − T_top) - - Solving this to constant Ṽ gives α as the ratio of required boost to available store energy: - - α = Q_boost / Q_discharge - = (T_forward − T_top) / (T_top − T_bottom) - - This expression quantifies the share of PTES output that is covered - by stored energy relative to the additional heating needed to meet - the desired forward temperature. - - Returns - ------- - xr.DataArray - The resulting fraction of PTES charge that must be further heated. - """ - if self.discharge_boosting_required: - result = ( - (self.forward_temperature - self.top_temperature_profile) - / (self.top_temperature_profile - self.bottom_temperature_profile) - ).where(self.forward_temperature > self.top_temperature_profile, 0) - - # Count how many snapshots require boosting - boosting_needed = ( - (self.forward_temperature > self.top_temperature_profile).sum().values - ) - total_snapshots = self.forward_temperature.size - - logger.info( - f"Discharge boosting (boost_per_discharge): Enabled. " - f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations " - f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})" - ) - return result - else: - logger.info( - f"Discharge boosting (boost_per_discharge): Not required. " - f"Returning boost_per_discharge=0 for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.zeros_like(self.forward_temperature) - - @property - def boost_per_charge(self) -> xr.DataArray: - """ - Calculate the additional boost energy ratio required during charging. - - .. note:: - Charge boosting is currently not implemented and will raise - NotImplementedError if charge_boosting_required is True. - - This calculates how much additional energy is needed to raise the PTES - top layer from the forward temperature to the design top temperature, - relative to the energy delivered by charging. - - Notes - ----- - To fill the storage from the return temperature all the way up to its - design top temperature, the total thermal energy required is split into: - - Q_charge = Ṽ·ρ·cₚ·(T_forward − T_bottom) - Q_boost = Ṽ·ρ·cₚ·(T_top − T_forward) - - - Q_charge is the energy delivered by charging to the forward setpoint. - - Q_boost is the extra boost energy needed to reach the design top temperature. - - Defining α as the ratio of boost energy to charge energy: - - α = Q_boost / Q_charge - = (T_top − T_forward) / (T_forward − T_return) - - Wherever the forward temperature meets or exceeds the design top temperature, - α is set to zero since no further boost is needed. - - Returns - ------- - xr.DataArray - The ratio of additional boost energy needed to the energy delivered by charging. - """ - if self.charge_boosting_required: - # Get the max top temperature value - max_top = ( - self.top_temperature - if isinstance(self.top_temperature, (int, float)) - else self.top_temperature_profile - ) - result = ( - ( - (max_top - self.forward_temperature) - / (self.forward_temperature - self.return_temperature) - ) - .where(self.forward_temperature < max_top, 0) - .clip(max=1) - ) - - # Count how many snapshots require boosting - boosting_needed = (self.forward_temperature < max_top).sum().values - total_snapshots = self.forward_temperature.size - - logger.info( - f"Charge boosting (boost_per_charge): Enabled. " - f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations " - f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})" - ) - return result - else: - logger.info( - f"Charge boosting (boost_per_charge): Not required. " - f"Returning boost_per_charge=0 for all {self.forward_temperature.size} snapshots and nodes" - ) - return xr.zeros_like(self.forward_temperature) diff --git a/scripts/build_ptes_operations/run.py b/scripts/build_ptes_operations/run.py index 347c5e5efe..c3a8ce6865 100644 --- a/scripts/build_ptes_operations/run.py +++ b/scripts/build_ptes_operations/run.py @@ -16,12 +16,8 @@ temperature difference. This captures how storage capacity varies when operating temperatures deviate from design conditions. -3. **Discharge boosting profiles**: When the PTES top temperature is below the - required forward temperature, additional heating (boosting) is needed during - discharge. This profile quantifies the ratio of boost energy to stored energy. - The outputs are used by ``prepare_sector_network.py`` to configure PTES storage -components, charger/discharger links, and optional resistive boosting infrastructure. +components and charger/discharger links. Relevant Settings ----------------- @@ -32,8 +28,6 @@ ptes: enable: true temperature_dependent_capacity: false # if true, e_max_pu varies with temperature difference (static but scaled if top/bottom are constant) - charge_boosting_required: false # currently not supported - discharge_resistive_boosting: false # if true, adds resistive boosting link for discharge boosting and disables heat pump boosting top_temperature: 90 # or "forward" for dynamic bottom_temperature: 35 # or "return" for dynamic design_top_temperature: 90 # used to compute design temperature difference for e_max_pu if temperature_dependent_capacity is true @@ -52,10 +46,6 @@ PTES top layer temperature profile (°C), clipped to design_top_temperature. - ``resources//ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc`` Normalized PTES capacity profiles (0-1 p.u.). -- ``resources//ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc`` - Discharge boosting ratio profiles. Values represent the ratio of boost energy - to discharge energy needed when T_forward > T_top. Only used when resistive - boosting is enabled. References ---------- @@ -92,10 +82,6 @@ ) logger.info(f"PTES configuration: {snakemake.params}") - # Discharge boosting profiles are calculated when resistive boosting is enabled. - # For heat pump-based boosting, add "ptes" to central heating heat sources instead and disable resistive boosting. - discharge_boosting_required: bool = snakemake.params.discharge_resistive_boosting - ptes_temperature_approximator = PtesTemperatureApproximator( forward_temperature=xr.open_dataarray( snakemake.input.central_heating_forward_temperature_profiles @@ -105,8 +91,6 @@ ), top_temperature=snakemake.params.top_temperature, bottom_temperature=snakemake.params.bottom_temperature, - charge_boosting_required=snakemake.params.charge_boosting_required, - discharge_boosting_required=discharge_boosting_required, temperature_dependent_capacity=snakemake.params.temperature_dependent_capacity, design_bottom_temperature=snakemake.params.design_bottom_temperature, design_top_temperature=snakemake.params.design_top_temperature, @@ -123,7 +107,3 @@ ptes_temperature_approximator.e_max_pu.to_netcdf( snakemake.output.ptes_e_max_pu_profiles ) - - ptes_temperature_approximator.boost_per_discharge.to_netcdf( - snakemake.output.ptes_boost_per_discharge_profiles - ) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 26f1a74399..9a9cafd5e6 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -21,9 +21,9 @@ class HeatSource(Enum): Have spatial/temporal constraints, require resource tracking via buses. May support direct utilisation or preheating depending on temperature. - **Sources with preheater** (PTES): - When source temperature is between return and forward temperatures, - can preheat return flow before heat pump provides final lift. + **Sources with boosting** (PTES, geothermal, river water): + When source temperature is below forward temperature, a heat pump + provides supplemental boost energy routed through the hp_output_bus. Attributes ---------- @@ -215,91 +215,6 @@ def get_lifetime(self, costs, heat_system) -> float: else: return float("inf") - def get_heat_pump_input_bus(self, nodes, heat_system: str) -> str: - """ - Get the secondary input bus for the heat pump link. - - The heat pump link has bus0 (electricity input), bus1 (heat output), - and optionally bus2 (heat source input). This method determines bus2 - based on the heat source type: - - - Inexhaustible sources: No bus2 (empty string) - - Limited sources: bus2 is the heat pump input bus (at return temperature if pre-heated, else at source temperature) - - Parameters - ---------- - nodes : pd.Index or list - The nodes for which to generate the bus name. - heat_system : str - The heat system identifier (e.g., 'urban central'). - - Returns - ------- - str - The bus2 name for the heat pump, or empty string if not applicable. - """ - if self.requires_bus: - return nodes + f" {self.heat_pump_input_carrier(heat_system)}" - else: - return "" - - def get_heat_pump_efficiency2(self, cop_heat_pump) -> float: - """ - Get the efficiency2 (heat source consumption) for the heat pump link. - - The heat pump link uses an inverted bus configuration where bus0 is the - heat output (not input) to attribute capital costs to the heat bus. This - means the link operates with negative flow (p_max_pu=0, p_min_pu < 0). - - - Inexhaustible sources (air, ground, sea_water): Returns 1.0 since no - resource tracking is needed. - - Limited sources (geothermal, river_water, ptes): Returns 1 - 1/COP to - track heat drawn from the source bus. - - For a standard heat pump energy balance: - - Q_output = Q_source + W_electricity - COP = Q_output / W_electricity - - The fraction of output heat from the source is: - - Q_source / Q_output = (COP - 1) / COP - - However, since bus0 is the output in our inverted configuration, PyPSA - interprets efficiency2 as: input_bus2 = p_bus0 * efficiency2 - - With negative p_bus0 (heat flowing out), we need efficiency2 to give the - correct heat drawn from the source. The formula simplifies to: - - efficiency2 = 1 - 1/COP - - This ensures that for each unit of heat output, the link draws the correct - amount from the heat source bus. - - Parameters - ---------- - cop_heat_pump : float or pd.Series - The coefficient of performance of the heat pump. - - Returns - ------- - float or pd.Series - 1 - 1/COP for limited sources (tracks heat drawn from source bus). - - Raises - ------ - NotImplementedError - If called for inexhaustible heat sources. - - See Also - -------- - prepare_sector_network.add_heat : Creates heat pump links with this efficiency. - """ - if self.requires_bus: - return 1 - (1 / cop_heat_pump.clip(lower=0.001)) - else: - return None - def heat_carrier(self, heat_system) -> str: """ Get the carrier name for heat from this source. @@ -317,31 +232,12 @@ def heat_carrier(self, heat_system) -> str: """ return f"{heat_system} {self} heat" - def preheater_input_carrier(self, heat_system) -> str: - """ - Get the carrier name for partially-cooled heat from this source. - - Used in cascading temperature utilisation: heat that has been used - for direct supply but still has usable thermal energy. - - Parameters - ---------- - heat_system : HeatSystem or str - The heat system (e.g., 'urban central'). - - Returns - ------- - str - Carrier name with '-pre-heater input' suffix in format '{heat_system} {source} pre-heater input'. - """ - return f"{self.heat_carrier(heat_system)} pre-heater input" - - def heat_pump_input_carrier(self, heat_system) -> str: + def hp_output_carrier(self, heat_system) -> str: """ - Get the carrier name for fully-cooled heat from this source. + Get the carrier name for heat produced by the boosting heat pump. - Represents heat pump input, for - final temperature lift by heat pump. + This bus receives HP output and is consumed by the utilisation link + at rate ``boosting_ratio`` per unit of source heat. Parameters ---------- @@ -351,13 +247,13 @@ def heat_pump_input_carrier(self, heat_system) -> str: Returns ------- str - Carrier name with '-heat-pump input' suffix in format '{heat_system} {source} heat-pump input'. + Carrier name in format '{heat_system} {source} heat hp output'. """ - return f"{self.heat_carrier(heat_system)} heat-pump input" + return f"{self.heat_carrier(heat_system)} hp output" - def preheater_input_bus(self, nodes, heat_system) -> str: + def hp_output_bus(self, nodes, heat_system) -> str: """ - Get bus name for partially-cooled heat at the given nodes. + Get the bus name for heat pump output at the given nodes. Parameters ---------- @@ -369,9 +265,12 @@ def preheater_input_bus(self, nodes, heat_system) -> str: Returns ------- str - Bus name combining nodes with medium-temperature carrier in format 'nodes + {heat_system} {source} pre-heater input'. + Bus name in format 'nodes + {heat_system} {source} heat hp output'. """ - return nodes + f" {self.preheater_input_carrier(heat_system)}" + if self.requires_bus: + return nodes + f" {self.hp_output_carrier(heat_system)}" + else: + return nodes + f" {heat_system} heat" def resource_bus(self, nodes, heat_system) -> str: """ diff --git a/scripts/lib/validation/config/sector.py b/scripts/lib/validation/config/sector.py index ef50600597..c833aacda3 100644 --- a/scripts/lib/validation/config/sector.py +++ b/scripts/lib/validation/config/sector.py @@ -47,8 +47,8 @@ class _PtesConfig(BaseModel): description="If True, enables boosting by resistive heaters instead of heat pumps. " "`prepare_sector_network` then adds the links ` urban central water pits resistive booster` " "and ` urban central water pits resistive heater stand-alone` and reroutes heat generation " - "from resistive heaters accordingly. The required boosting energy is computed in `build_ptes_operations` " - "as `ptes_boost_per_discharge_profiles_base_s_.nc`.", + "from resistive heaters accordingly. The required boosting alpha is computed from the unified " + "heat source utilisation profiles in `build_heat_source_utilisation_profiles`.", ) top_temperature: float | Literal["forward"] = Field( 90, diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 6b5edf6bcb..1c1a1f373f 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -2753,11 +2753,9 @@ def add_heat( n: pypsa.Network, costs: pd.DataFrame, cop_profiles_file: str, - heat_source_direct_utilisation_profile_file: str, - heat_source_preheater_utilisation_profile_file: str, + heat_source_boosting_profile_file: str, hourly_heat_demand_total_file: str, ptes_e_max_pu_file: str, - ptes_boost_per_discharge_profile_file: str, ates_e_nom_max: str, ates_capex_as_fraction_of_geothermal_heat_source: float, ates_recovery_factor: float, @@ -2788,8 +2786,9 @@ def add_heat( DataFrame containing cost information for different technologies cop_profiles_file : str Path to NetCDF file containing coefficient of performance (COP) values for heat pumps - direct_heat_source_utilisation_profile_file : str - Path to NetCDF file containing direct heat source utilisation profiles + heat_source_boosting_profile_file : str + Path to NetCDF file containing heat source boosting ratio profiles (indexed by + time, name, heat_source). Values are the HP heat required per unit of source heat. hourly_heat_demand_total_file : str Path to CSV file containing hourly heat demand data ptes_supplemental_heating_required_file: str @@ -2853,15 +2852,9 @@ def add_heat( cop = xr.open_dataarray(cop_profiles_file) - heat_source_direct_utilisation_profile = ( - xr.open_dataarray(heat_source_direct_utilisation_profile_file) - if len(heat_source_direct_utilisation_profile_file) > 0 - else None - ) - - heat_source_preheater_utilisation_profile = ( - xr.open_dataarray(heat_source_preheater_utilisation_profile_file) - if len(heat_source_preheater_utilisation_profile_file) > 0 + heat_source_boosting_profile = ( + xr.open_dataarray(heat_source_boosting_profile_file) + if len(heat_source_boosting_profile_file) > 0 else None ) @@ -3254,7 +3247,7 @@ def add_heat( heat_carrier = heat_source.heat_carrier(heat_system) if heat_source.requires_bus: - # add heat source carrier and bus + # add heat source carrier and buses n.add("Carrier", heat_carrier) n.add( "Bus", @@ -3263,60 +3256,37 @@ def add_heat( carrier=heat_carrier, ) - preheater_utilisation_profile = ( - heat_source_preheater_utilisation_profile.sel( - heat_source=heat_source.value, name=nodes - ) - .transpose("time", "name") - .to_pandas() - ) - n.add( - "Bus", - heat_source.preheater_input_bus(nodes, heat_system), - location=nodes, - carrier=heat_source.preheater_input_carrier(heat_system), + "Carrier", + heat_source.hp_output_carrier(heat_system), ) - n.add( "Bus", - heat_source.get_heat_pump_input_bus(nodes, heat_system), + heat_source.hp_output_bus(nodes, heat_system), location=nodes, - carrier=heat_source.heat_pump_input_carrier(heat_system), + carrier=heat_source.hp_output_carrier(heat_system), ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} {heat_source} heat preheater", - bus0=heat_source.preheater_input_bus(nodes, heat_system), - bus1=nodes + f" {heat_system} heat", - bus2=heat_source.get_heat_pump_input_bus(nodes, heat_system), - efficiency=preheater_utilisation_profile, - efficiency2=1 - preheater_utilisation_profile, - carrier=f"{heat_system} {heat_source} heat preheater", - p_nom_extendable=True, - ) - - direct_utilisation_profile = ( - heat_source_direct_utilisation_profile.sel( + boosting_profile = ( + heat_source_boosting_profile.sel( heat_source=heat_source.value, name=nodes ) .transpose("time", "name") .to_pandas() ) - # add link for direct usage of heat source when source temperature exceeds forward temperature - + # Utilisation link: source heat (bus0) + HP boost (bus2) → DH (bus1). + # Per unit of source heat, eff1 = 1 + boosting_ratio goes to DH and + # eff2 = -boosting_ratio is drawn from the hp_output_bus. n.add( "Link", nodes, suffix=f" {heat_system} {heat_source} heat utilisation", bus0=heat_source.resource_bus(nodes, heat_system), bus1=nodes + f" {heat_system} heat", - bus2=heat_source.preheater_input_bus(nodes, heat_system), - efficiency=direct_utilisation_profile, - efficiency2=1 - direct_utilisation_profile, + bus2=heat_source.hp_output_bus(nodes, heat_system), + efficiency=1 + boosting_profile, + efficiency2=-boosting_profile, carrier=f"{heat_system} {heat_source} heat utilisation", p_nom_extendable=True, ) @@ -3358,16 +3328,16 @@ def add_heat( "ptes" ]["discharge_resistive_boosting"] ): + # Inexhaustible sources (air, ground): simple HP → DH heat bus. + # For limited sources the HP is added inside the requires_bus block above. n.add( "Link", nodes, suffix=f" {heat_system} {heat_source} heat pump", bus0=nodes + f" {heat_system} heat", bus1=nodes, - bus2=heat_source.get_heat_pump_input_bus(nodes, heat_system), carrier=f"{heat_system} {heat_source} heat pump", efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(), - efficiency2=heat_source.get_heat_pump_efficiency2(cop_heat_pump), capital_cost=costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor, p_min_pu=-(cop_heat_pump > 0).squeeze().astype(float), @@ -3387,11 +3357,12 @@ def add_heat( "discharge_resistive_boosting" ] ): - ptes_boost_per_discharge_profiles = ( - xr.open_dataarray(ptes_boost_per_discharge_profile_file) - .sel(name=nodes) + boosting_profile = ( + heat_source_boosting_profile.sel( + heat_source=HeatSource.PTES.value, name=nodes + ) + .transpose("time", "name") .to_pandas() - .reindex(index=n.snapshots) ) n.add( @@ -3408,15 +3379,14 @@ def add_heat( suffix=f" {heat_system} water pits resistive booster", bus0=nodes + f" {heat_system} heat", bus1=nodes + f" {heat_system} resistive heat", - bus2=ptes_heat_source.preheater_input_bus(nodes, heat_system), + bus2=ptes_heat_source.hp_output_bus(nodes, heat_system), # eff = 1 - eff2 (energy conservation) - efficiency=ptes_boost_per_discharge_profiles - / (ptes_boost_per_discharge_profiles + 1), + efficiency=boosting_profile / (boosting_profile + 1), # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating # (similar to HP balance: p_el x COP = p_source + p_el ) - efficiency2=1 / (ptes_boost_per_discharge_profiles + 1), + efficiency2=1 / (boosting_profile + 1), p_nom_extendable=True, - p_min_pu=-(ptes_boost_per_discharge_profiles > 0).astype(float), + p_min_pu=-(boosting_profile > 0).astype(float), carrier=f"{heat_system} water pits resistive booster", p_max_pu=0, ) @@ -6416,8 +6386,7 @@ def add_import_options( n=n, costs=costs, cop_profiles_file=snakemake.input.cop_profiles, - heat_source_direct_utilisation_profile_file=snakemake.input.heat_source_direct_utilisation_profiles, - heat_source_preheater_utilisation_profile_file=snakemake.input.heat_source_preheater_utilisation_profiles, + heat_source_boosting_profile_file=snakemake.input.heat_source_boosting_profiles, hourly_heat_demand_total_file=snakemake.input.hourly_heat_demand_total, ptes_e_max_pu_file=snakemake.input.ptes_e_max_pu_profiles, ates_e_nom_max=snakemake.input.ates_potentials, @@ -6431,7 +6400,6 @@ def add_import_options( "recovery_factor" ], enable_ates=snakemake.params.sector["district_heating"]["ates"]["enable"], - ptes_boost_per_discharge_profile_file=snakemake.input.ptes_boost_per_discharge_profiles, district_heat_share_file=snakemake.input.district_heat_share, solar_thermal_total_file=snakemake.input.solar_thermal_total, retro_cost_file=snakemake.input.retro_cost, From c253e0d71846243076906f860fb482c381f1f2c9 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 7 Apr 2026 15:03:21 +0200 Subject: [PATCH 10/17] feat: add interlayer heat transfer coefficient to PTES operations --- rules/build_sector.smk | 6 ++++++ scripts/build_ptes_operations/run.py | 1 + 2 files changed, 7 insertions(+) diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 29143389c2..eca8197c46 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -738,6 +738,12 @@ rule build_ptes_operations: "design_bottom_temperature", ), layered=config_provider("sector", "district_heating", "ptes", "layered"), + interlayer_heat_transfer_coefficient=config_provider( + "sector", + "district_heating", + "ptes", + "interlayer_heat_transfer_coefficient", + ), input: central_heating_forward_temperature_profiles=resources( "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc" diff --git a/scripts/build_ptes_operations/run.py b/scripts/build_ptes_operations/run.py index d16ddb3565..683738e530 100644 --- a/scripts/build_ptes_operations/run.py +++ b/scripts/build_ptes_operations/run.py @@ -90,6 +90,7 @@ design_top_temperature=snakemake.params.design_top_temperature, design_bottom_temperature=snakemake.params.design_bottom_temperature, design_standing_losses=0.0, # placeholder; actual value set from costs data + interlayer_heat_transfer_coefficient=snakemake.params.interlayer_heat_transfer_coefficient, ) # Write single dataset with all pre-computed PTES parameters From ad5f1de471a82c63b374ac6e6e6f5074a8b01889 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 7 Apr 2026 18:11:00 +0200 Subject: [PATCH 11/17] update boosting logic --- .../build_heat_source_utilisation_profiles.py | 78 +++---------------- scripts/prepare_sector_network.py | 20 +++-- 2 files changed, 21 insertions(+), 77 deletions(-) diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py index 18a4e0399b..9eb05f7fa1 100644 --- a/scripts/build_heat_source_utilisation_profiles.py +++ b/scripts/build_heat_source_utilisation_profiles.py @@ -108,63 +108,10 @@ def get_source_temperature( return xr.open_dataarray(snakemake_input[f"temp_{heat_source_name}"]) -def get_heat_pump_cooling( - heat_source_name: str, - default_heat_source_cooling: float, - snakemake_input: dict, - return_temperature: xr.DataArray = None, -) -> float | xr.DataArray: - """ - Get the additional heat source cooling (temperature drop) through the heat pump. - - This is the number of Kelvin below the return temperature that the heat pump - can extract from the source, i.e. the HP cools the source from T_source down - to T_return − heat_pump_cooling. - - For PTES, equals T_return − T_bottom (the bottom layer is the HP's cold side). - For other sources, uses the constant config value ``heat_source_cooling``. - - Parameters - ---------- - heat_source_name : str - Name of the heat source (e.g., 'ptes', 'geothermal', 'air'). - default_heat_source_cooling : float - Default heat source cooling in Kelvin, from config. - snakemake_input : dict - Snakemake input files, may contain PTES bottom temperature profiles. - return_temperature : xr.DataArray, optional - District heating return temperature profiles in °C. Required for PTES. - - Returns - ------- - float | xr.DataArray - Heat source cooling in Kelvin. - - Raises - ------ - ValueError - If heat source is PTES but required profiles are not provided. - """ - if heat_source_name == "ptes": - if "temp_ptes_bottom" not in snakemake_input.keys(): - raise ValueError( - "PTES heat source requires bottom temperature profile " - "(temp_ptes_bottom) to calculate heat source cooling." - ) - if return_temperature is None: - raise ValueError( - "PTES heat source requires return_temperature to calculate heat pump cooling." - ) - ptes_bottom_temperature = xr.open_dataarray(snakemake_input["temp_ptes_bottom"]) - return return_temperature - ptes_bottom_temperature - return default_heat_source_cooling - - def get_boosting_profile( source_temperature: float | xr.DataArray, forward_temperature: xr.DataArray, return_temperature: xr.DataArray, - heat_pump_cooling: float | xr.DataArray, ) -> xr.DataArray: """ Calculate the boosting ratio: HP heat needed per unit of source heat. @@ -200,19 +147,18 @@ def get_boosting_profile( Alpha profile: 0 where direct use is possible, positive ratio otherwise. Shape matches forward_temperature. """ - delta_t = xr.where( - source_temperature < return_temperature, - heat_pump_cooling, - source_temperature - return_temperature + heat_pump_cooling, - ) - # Prevent division by zero (e.g. T_source = T_bottom with resistive boosting). - # In that case the source cannot be economically used: alpha → ∞ is approximated - # by a very large number so the optimiser will not dispatch this configuration. - safe_delta_t = delta_t.where(delta_t > 1e-6, other=1e6) return xr.where( source_temperature >= forward_temperature, + # no boosting needed if source_temp > forward_temp 0.0, - (forward_temperature - source_temperature) / safe_delta_t, + xr.where( + source_temperature < return_temperature, + # source does not pre-heat return flow if below return temp + 1, + # if source is between return and forward temp, it pre-heats return flow (part of heat is utilised directly, part is boosted by HP) + (forward_temperature - source_temperature) + / (source_temperature - return_temperature), + ).clip(min=0, max=1), ) @@ -254,12 +200,6 @@ def get_boosting_profile( ), forward_temperature=central_heating_forward_temperature, return_temperature=central_heating_return_temperature, - heat_pump_cooling=get_heat_pump_cooling( - heat_source_name=heat_source_key, - default_heat_source_cooling=snakemake.params.heat_source_cooling, - snakemake_input=snakemake.input, - return_temperature=central_heating_return_temperature, - ), ).assign_coords(heat_source=heat_source_key) for heat_source_key in heat_sources ], diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 1c1a1f373f..415e15eb4c 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3275,20 +3275,24 @@ def add_heat( .to_pandas() ) - # Utilisation link: source heat (bus0) + HP boost (bus2) → DH (bus1). - # Per unit of source heat, eff1 = 1 + boosting_ratio goes to DH and - # eff2 = -boosting_ratio is drawn from the hp_output_bus. + # Utilisation link (reverse operation, p < 0): + # Consumes from resource bus (bus1) and hp_output bus (bus2), + # produces onto DH heat bus (bus0). + # Per unit of source heat: 1 source + b HP → (1 + b) DH heat. + # Energy conserved: eff + eff2 = 1/(1+b) + b/(1+b) = 1. n.add( "Link", nodes, suffix=f" {heat_system} {heat_source} heat utilisation", - bus0=heat_source.resource_bus(nodes, heat_system), - bus1=nodes + f" {heat_system} heat", + bus0=nodes + f" {heat_system} heat", + bus1=heat_source.resource_bus(nodes, heat_system), bus2=heat_source.hp_output_bus(nodes, heat_system), - efficiency=1 + boosting_profile, - efficiency2=-boosting_profile, + efficiency=1 / (1 + boosting_profile), + efficiency2=boosting_profile / (1 + boosting_profile), carrier=f"{heat_system} {heat_source} heat utilisation", p_nom_extendable=True, + p_min_pu=-1, + p_max_pu=0, ) if heat_source.requires_generator: @@ -3334,7 +3338,7 @@ def add_heat( "Link", nodes, suffix=f" {heat_system} {heat_source} heat pump", - bus0=nodes + f" {heat_system} heat", + bus0=heat_source.hp_output_bus(nodes, heat_system), bus1=nodes, carrier=f"{heat_system} {heat_source} heat pump", efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(), From c4168244b0e2367909356ea29c8d2cfa8cf8b8fe Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 8 Apr 2026 17:43:50 +0200 Subject: [PATCH 12/17] fix boosting logic: 1 utilisation link distributes to UCH and HP input --- scripts/definitions/heat_source.py | 8 ++++---- scripts/prepare_sector_network.py | 24 ++++++++++++------------ 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 9a9cafd5e6..4f831c0095 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -232,7 +232,7 @@ def heat_carrier(self, heat_system) -> str: """ return f"{heat_system} {self} heat" - def hp_output_carrier(self, heat_system) -> str: + def hp_input_carrier(self, heat_system) -> str: """ Get the carrier name for heat produced by the boosting heat pump. @@ -249,9 +249,9 @@ def hp_output_carrier(self, heat_system) -> str: str Carrier name in format '{heat_system} {source} heat hp output'. """ - return f"{self.heat_carrier(heat_system)} hp output" + return f"{self.heat_carrier(heat_system)} heat pump input" - def hp_output_bus(self, nodes, heat_system) -> str: + def hp_input_bus(self, nodes, heat_system) -> str: """ Get the bus name for heat pump output at the given nodes. @@ -268,7 +268,7 @@ def hp_output_bus(self, nodes, heat_system) -> str: Bus name in format 'nodes + {heat_system} {source} heat hp output'. """ if self.requires_bus: - return nodes + f" {self.hp_output_carrier(heat_system)}" + return nodes + f" {self.hp_input_carrier(heat_system)}" else: return nodes + f" {heat_system} heat" diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 415e15eb4c..7c78fe85e0 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3258,13 +3258,13 @@ def add_heat( n.add( "Carrier", - heat_source.hp_output_carrier(heat_system), + heat_source.hp_input_carrier(heat_system), ) n.add( "Bus", - heat_source.hp_output_bus(nodes, heat_system), + heat_source.hp_input_bus(nodes, heat_system), location=nodes, - carrier=heat_source.hp_output_carrier(heat_system), + carrier=heat_source.hp_input_carrier(heat_system), ) boosting_profile = ( @@ -3284,15 +3284,13 @@ def add_heat( "Link", nodes, suffix=f" {heat_system} {heat_source} heat utilisation", - bus0=nodes + f" {heat_system} heat", - bus1=heat_source.resource_bus(nodes, heat_system), - bus2=heat_source.hp_output_bus(nodes, heat_system), - efficiency=1 / (1 + boosting_profile), - efficiency2=boosting_profile / (1 + boosting_profile), + bus0=heat_source.resource_bus(nodes, heat_system), + bus1=nodes + f" {heat_system} heat", + bus2=heat_source.hp_input_bus(nodes, heat_system), + efficiency=1 - boosting_profile, + efficiency2=boosting_profile, carrier=f"{heat_system} {heat_source} heat utilisation", p_nom_extendable=True, - p_min_pu=-1, - p_max_pu=0, ) if heat_source.requires_generator: @@ -3338,10 +3336,12 @@ def add_heat( "Link", nodes, suffix=f" {heat_system} {heat_source} heat pump", - bus0=heat_source.hp_output_bus(nodes, heat_system), + bus0=nodes + f" {heat_system} heat", bus1=nodes, + bus2=heat_source.hp_input_bus(nodes, heat_system), carrier=f"{heat_system} {heat_source} heat pump", efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(), + efficiency2=1 - 1 / cop_heat_pump.clip(lower=0.001).squeeze(), capital_cost=costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor, p_min_pu=-(cop_heat_pump > 0).squeeze().astype(float), @@ -3383,7 +3383,7 @@ def add_heat( suffix=f" {heat_system} water pits resistive booster", bus0=nodes + f" {heat_system} heat", bus1=nodes + f" {heat_system} resistive heat", - bus2=ptes_heat_source.hp_output_bus(nodes, heat_system), + bus2=ptes_heat_source.hp_input_bus(nodes, heat_system), # eff = 1 - eff2 (energy conservation) efficiency=boosting_profile / (boosting_profile + 1), # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating From d52a69f05d97b7dcbed02f6a366abd3ca2cdac65 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 8 Apr 2026 17:59:02 +0200 Subject: [PATCH 13/17] fix faulty HP input for inexhaustible heat sources --- scripts/definitions/heat_source.py | 2 +- scripts/prepare_sector_network.py | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 4f831c0095..2815778ed6 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -270,7 +270,7 @@ def hp_input_bus(self, nodes, heat_system) -> str: if self.requires_bus: return nodes + f" {self.hp_input_carrier(heat_system)}" else: - return nodes + f" {heat_system} heat" + return "" def resource_bus(self, nodes, heat_system) -> str: """ diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 7c78fe85e0..d71cda7c6d 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3275,11 +3275,7 @@ def add_heat( .to_pandas() ) - # Utilisation link (reverse operation, p < 0): - # Consumes from resource bus (bus1) and hp_output bus (bus2), - # produces onto DH heat bus (bus0). - # Per unit of source heat: 1 source + b HP → (1 + b) DH heat. - # Energy conserved: eff + eff2 = 1/(1+b) + b/(1+b) = 1. + # Utilisation link distributes part of heat above return temperature to forward flow and rest to heat pump for boosting n.add( "Link", nodes, From fd178c269675e9448e944e8e6932ade139e93147 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 8 Apr 2026 18:11:33 +0200 Subject: [PATCH 14/17] clean up docstrings --- .../build_heat_source_utilisation_profiles.py | 51 ++++++++----------- scripts/definitions/heat_source.py | 29 ++++++----- scripts/prepare_sector_network.py | 13 ++--- 3 files changed, 40 insertions(+), 53 deletions(-) diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py index 9eb05f7fa1..59986fa1cc 100644 --- a/scripts/build_heat_source_utilisation_profiles.py +++ b/scripts/build_heat_source_utilisation_profiles.py @@ -2,29 +2,25 @@ # # SPDX-License-Identifier: MIT """ -Build heat source alpha profiles for district heating networks. +Build heat source boosting ratio profiles for district heating networks. -This script calculates the boosting ratio (alpha) for each heat source: how +This script calculates the boosting ratio (b) for each heat source: how much heat pump output is needed per unit of source heat to reach the forward temperature of the district heating network. -**Alpha profile**: For each heat source and timestep, alpha is: +**Boosting ratio profile**: For each heat source and timestep, b is: -- 0 when T_source ≥ T_forward (direct use possible, no HP boost needed) -- (T_forward − T_source) / delta_T otherwise, where delta_T is: - - heat_pump_cooling if T_source < T_return - - T_source − T_return + heat_pump_cooling if T_return ≤ T_source < T_forward +- 0 when T_source ≥ T_forward (direct use, no HP boost needed) +- 1 when T_source < T_return (source cannot preheat return flow) +- (T_forward − T_source) / (T_source − T_return) otherwise, clipped to [0, 1] -heat_pump_cooling is the additional temperature drop achievable by the heat -pump beyond the return temperature (i.e. how far below T_return the HP can -extract from the source). For PTES it equals T_return − T_bottom; for other -sources it is the constant config value ``heat_source_cooling``. +The boosting ratio is consumed by ``prepare_sector_network.py`` to set the +efficiencies of the heat source utilisation link (forward, p ≥ 0): -The alpha profile is consumed by ``prepare_sector_network.py`` to set the -efficiencies of the heat source utilisation links: +- bus0 (source) → bus1 (DH heat) at efficiency (1 − b) +- bus0 (source) → bus2 (HP input bus) at efficiency2 = b -- bus0 (source) → bus1 (DH heat) at efficiency 1 + alpha -- bus2 (HP output bus) drawn at efficiency −alpha per unit of source +Energy is conserved: (1 − b) + b = 1. Relevant Settings ----------------- @@ -36,7 +32,6 @@ - air - geothermal district_heating: - heat_source_cooling: 6 # K geothermal: constant_temperature_celsius: 65 @@ -52,7 +47,7 @@ ------- - ``resources//heat_source_boosting_profiles_base_s_{clusters}_{planning_horizons}.nc`` Boosting ratio profiles indexed by (time, name, heat_source). - Values: 0 when T_source ≥ T_forward; (T_fwd − T_src) / delta_T otherwise. + Values in [0, 1]: 0 = direct use, 1 = full HP boosting required. """ import logging @@ -116,19 +111,15 @@ def get_boosting_profile( """ Calculate the boosting ratio: HP heat needed per unit of source heat. - Alpha represents the ratio of heat pump output required to source heat input - in order to reach the district heating forward temperature: + The boosting ratio b represents the fraction of source heat that must be + routed through the heat pump (rather than used directly) to reach the + district heating forward temperature: - alpha = 0 if T_source ≥ T_forward - alpha = (T_fwd − T_src) / delta_T otherwise + b = 0 if T_source ≥ T_forward + b = 1 if T_source < T_return + b = (T_forward − T_source) / (T_source − T_return) otherwise - where delta_T is: - - delta_T = heat_pump_cooling if T_source < T_return - delta_T = T_source − T_return + heat_pump_cooling if T_return ≤ T_source < T_forward - - For PTES (where heat_pump_cooling = T_return − T_bottom): - delta_T = T_source − T_bottom (when T_source ≥ T_return, which is the normal case) + The result is clipped to [0, 1]. Parameters ---------- @@ -138,13 +129,11 @@ def get_boosting_profile( District heating forward temperature in °C, indexed by (time, name). return_temperature : xr.DataArray District heating return temperature in °C, indexed by (time, name). - heat_pump_cooling : float | xr.DataArray - Additional temperature drop achievable by HP beyond return temperature (K). Returns ------- xr.DataArray - Alpha profile: 0 where direct use is possible, positive ratio otherwise. + Boosting ratio profile in [0, 1]: 0 = direct use, 1 = full HP boosting. Shape matches forward_temperature. """ return xr.where( diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 2815778ed6..fe7b3306d4 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -10,7 +10,7 @@ class HeatSource(Enum): """ - Enumeration representing different heat sources for heat pumps and direct utilisation. + Enumeration representing different heat sources for heat pumps and utilisation. Heat sources are categorized by their characteristics: @@ -19,11 +19,9 @@ class HeatSource(Enum): **Limited sources requiring a bus** (GEOTHERMAL, RIVER_WATER, PTES): Have spatial/temporal constraints, require resource tracking via buses. - May support direct utilisation or preheating depending on temperature. - - **Sources with boosting** (PTES, geothermal, river water): - When source temperature is below forward temperature, a heat pump - provides supplemental boost energy routed through the hp_output_bus. + A utilisation link splits source heat into direct DH contribution + (fraction 1 − b) and HP input (fraction b), where b is the + boosting ratio from ``build_heat_source_utilisation_profiles``. Attributes ---------- @@ -43,7 +41,7 @@ class HeatSource(Enum): See Also -------- HeatSystem : Defines heat system types (urban central, urban decentral, rural). - build_heat_source_utilisation_profiles : Calculates utilisation profiles for heat sources. + build_heat_source_utilisation_profiles : Calculates boosting ratio profiles for heat sources. build_cop_profiles : Calculates COP profiles for heat pumps using these sources. """ @@ -234,10 +232,10 @@ def heat_carrier(self, heat_system) -> str: def hp_input_carrier(self, heat_system) -> str: """ - Get the carrier name for heat produced by the boosting heat pump. + Get the carrier name for the boosting heat pump input bus. - This bus receives HP output and is consumed by the utilisation link - at rate ``boosting_ratio`` per unit of source heat. + For limited sources, the utilisation link routes a fraction (the boosting + ratio) of source heat to this bus, where it is consumed by the heat pump. Parameters ---------- @@ -247,13 +245,18 @@ def hp_input_carrier(self, heat_system) -> str: Returns ------- str - Carrier name in format '{heat_system} {source} heat hp output'. + Carrier name in format '{heat_system} {source} heat heat pump input'. """ return f"{self.heat_carrier(heat_system)} heat pump input" def hp_input_bus(self, nodes, heat_system) -> str: """ - Get the bus name for heat pump output at the given nodes. + Get the bus name for the heat pump input at the given nodes. + + For limited sources (requires_bus=True), returns the dedicated HP input + bus where the utilisation link deposits the boosting fraction of source + heat. For inexhaustible sources, returns an empty string (no bus2 + connection needed since the HP draws from ambient implicitly). Parameters ---------- @@ -265,7 +268,7 @@ def hp_input_bus(self, nodes, heat_system) -> str: Returns ------- str - Bus name in format 'nodes + {heat_system} {source} heat hp output'. + Bus name for limited sources, empty string for inexhaustible sources. """ if self.requires_bus: return nodes + f" {self.hp_input_carrier(heat_system)}" diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index d71cda7c6d..8e9eb098a7 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -2788,11 +2788,12 @@ def add_heat( Path to NetCDF file containing coefficient of performance (COP) values for heat pumps heat_source_boosting_profile_file : str Path to NetCDF file containing heat source boosting ratio profiles (indexed by - time, name, heat_source). Values are the HP heat required per unit of source heat. + time, name, heat_source). Values in [0, 1] represent the fraction of source + heat that must be routed through the heat pump for temperature boosting. hourly_heat_demand_total_file : str Path to CSV file containing hourly heat demand data - ptes_supplemental_heating_required_file: str - Path to CSV file indicating when supplemental heating for thermal energy storage (TES) is needed + ptes_e_max_pu_file : str + Path to NetCDF file containing PTES energy capacity profiles district_heat_share_file : str Path to CSV file containing district heating share information solar_thermal_total_file : str @@ -2808,8 +2809,6 @@ def add_heat( params : dict Dictionary containing parameters including: - heat_sources - - heat_utilisation_potentials - - direct_utilisation_heat_sources pop_weighted_energy_totals : pd.DataFrame Population-weighted energy totals by region heating_efficiencies : pd.DataFrame @@ -3380,15 +3379,11 @@ def add_heat( bus0=nodes + f" {heat_system} heat", bus1=nodes + f" {heat_system} resistive heat", bus2=ptes_heat_source.hp_input_bus(nodes, heat_system), - # eff = 1 - eff2 (energy conservation) efficiency=boosting_profile / (boosting_profile + 1), - # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating - # (similar to HP balance: p_el x COP = p_source + p_el ) efficiency2=1 / (boosting_profile + 1), p_nom_extendable=True, p_min_pu=-(boosting_profile > 0).astype(float), carrier=f"{heat_system} water pits resistive booster", - p_max_pu=0, ) n.add( From 02b66f12a64c712431d1cb46f9106a5244cd3487 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Mon, 13 Apr 2026 15:50:52 +0200 Subject: [PATCH 15/17] fix util link ratios --- scripts/definitions/heat_source.py | 101 ++++++++++++++++++++++++----- scripts/prepare_sector_network.py | 76 +++++++++++++--------- 2 files changed, 130 insertions(+), 47 deletions(-) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index fe7b3306d4..058ca2a23d 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -106,6 +106,25 @@ def requires_bus(self) -> bool: else: return False + @property + def supports_preheating(self) -> bool: + """ + Check if the heat source supports preheating district heating return flow. + + Preheating sources have temperatures between return and forward + (T_ret <= T_src < T_fwd). The source directly heats return flow, and a + heat pump boosts only the remaining fraction to forward temperature. + + Non-preheating limited sources (e.g., river_water) act as HP evaporator + input: all source heat enters the HP cold side. + + Returns + ------- + bool + True for PTES and GEOTHERMAL, False otherwise. + """ + return self in [HeatSource.PTES, HeatSource.GEOTHERMAL] + @property def requires_generator(self) -> bool: """ @@ -230,12 +249,13 @@ def heat_carrier(self, heat_system) -> str: """ return f"{heat_system} {self} heat" - def hp_input_carrier(self, heat_system) -> str: + def intermediate_carrier(self, heat_system) -> str: """ - Get the carrier name for the boosting heat pump input bus. + Get the carrier name for the intermediate bus between utilisation link and HP. - For limited sources, the utilisation link routes a fraction (the boosting - ratio) of source heat to this bus, where it is consumed by the heat pump. + For preheating sources, the HP produces onto this bus and the utilisation + link consumes from it. For evaporator sources, the utilisation link + produces onto this bus and the HP draws from it as cold-side input. Parameters ---------- @@ -245,18 +265,19 @@ def hp_input_carrier(self, heat_system) -> str: Returns ------- str - Carrier name in format '{heat_system} {source} heat heat pump input'. + Carrier name in format '{heat_system} {source} heat heat pump output'. """ - return f"{self.heat_carrier(heat_system)} heat pump input" + return f"{self.heat_carrier(heat_system)} heat pump output" - def hp_input_bus(self, nodes, heat_system) -> str: + def intermediate_bus(self, nodes, heat_system) -> str: """ - Get the bus name for the heat pump input at the given nodes. + Get the intermediate bus connecting the utilisation link and heat pump. - For limited sources (requires_bus=True), returns the dedicated HP input - bus where the utilisation link deposits the boosting fraction of source - heat. For inexhaustible sources, returns an empty string (no bus2 - connection needed since the HP draws from ambient implicitly). + For limited sources (requires_bus=True), returns the dedicated + intermediate bus. For preheating sources, the HP produces onto this bus + and the utilisation link consumes from it. For evaporator sources, the + utilisation link produces onto this bus and the HP draws from it as + cold-side input. For inexhaustible sources, returns an empty string. Parameters ---------- @@ -268,13 +289,63 @@ def hp_input_bus(self, nodes, heat_system) -> str: Returns ------- str - Bus name for limited sources, empty string for inexhaustible sources. + Bus name for evaporator sources, empty string otherwise. """ - if self.requires_bus: - return nodes + f" {self.hp_input_carrier(heat_system)}" + if self.requires_bus and not self.supports_preheating: + return nodes + f" {self.intermediate_carrier(heat_system)}" else: return "" + def hp_output_bus(self, nodes, heat_system) -> str: + """ + Get the bus where the heat pump outputs its heat (bus0). + + For preheating sources, the HP outputs to the intermediate bus. + For all other sources, the HP outputs directly to the DH heat bus. + This always represents bus0 of the (reverse-operating) HP link, + preserving correct investment sizing. + + Parameters + ---------- + nodes : pd.Index or str + Node identifier(s). + heat_system : HeatSystem or str + The heat system (e.g., 'urban central'). + + Returns + ------- + str + The HP output bus name. + """ + if self.supports_preheating: + return nodes + f" {self.intermediate_carrier(heat_system)}" + else: + return nodes + f" {heat_system} heat" + + def hp_eff2(self, cop): + """ + Get efficiency2 for the heat pump link. + + For preheating sources, eff2=0 (bus2 is unused). For evaporator + sources, eff2 = 1 - 1/COP represents the fraction of HP output + drawn from the cold side. For inexhaustible sources the same + formula applies but bus2="" so the value is irrelevant. + + Parameters + ---------- + cop : float or pd.DataFrame + The coefficient of performance of the heat pump. + + Returns + ------- + float or pd.DataFrame + 0 for preheating sources, 1 - 1/COP otherwise. + """ + if self.supports_preheating: + return 0 + else: + return 1 - 1 / cop + def resource_bus(self, nodes, heat_system) -> str: """ Get the primary resource bus for heat from this source, e.g. `DE0 0 urban central geothermal heat` or `DE0 0 urban central ptes heat`. diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 8e9eb098a7..d4e01714a4 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3232,15 +3232,19 @@ def add_heat( costs_name_heat_pump = heat_system.heat_pump_costs_name(heat_source) cop_heat_pump = ( - cop.sel( - heat_system=heat_system.system_type.value, - heat_source=heat_source.value, - name=nodes, + ( + cop.sel( + heat_system=heat_system.system_type.value, + heat_source=heat_source.value, + name=nodes, + ) + .transpose("time", "name") + .to_pandas() + if options["time_dep_hp_cop"] + else costs.loc[[costs_name_heat_pump], ["efficiency"]] ) - .transpose("time", "name") - .to_pandas() - if options["time_dep_hp_cop"] - else costs.loc[[costs_name_heat_pump], ["efficiency"]] + .clip(lower=0.001) + .squeeze() ) heat_carrier = heat_source.heat_carrier(heat_system) @@ -3257,33 +3261,41 @@ def add_heat( n.add( "Carrier", - heat_source.hp_input_carrier(heat_system), + heat_source.intermediate_carrier(heat_system), ) n.add( "Bus", - heat_source.hp_input_bus(nodes, heat_system), + nodes + f" {heat_source.intermediate_carrier(heat_system)}", location=nodes, - carrier=heat_source.hp_input_carrier(heat_system), + carrier=heat_source.intermediate_carrier(heat_system), ) - boosting_profile = ( - heat_source_boosting_profile.sel( - heat_source=heat_source.value, name=nodes + if heat_source.supports_preheating: + boosting_profile = ( + heat_source_boosting_profile.sel( + heat_source=heat_source.value, name=nodes + ) + .transpose("time", "name") + .to_pandas() ) - .transpose("time", "name") - .to_pandas() - ) - - # Utilisation link distributes part of heat above return temperature to forward flow and rest to heat pump for boosting + # Utilisation link: resource_bus → DH heat + intermediate_bus + # Preheating: source directly heats DH (eff=1+b), HP boosts remainder (eff2=-b) + # Evaporator: all source → HP cold side (eff=0, eff2=1) n.add( "Link", nodes, suffix=f" {heat_system} {heat_source} heat utilisation", bus0=heat_source.resource_bus(nodes, heat_system), bus1=nodes + f" {heat_system} heat", - bus2=heat_source.hp_input_bus(nodes, heat_system), - efficiency=1 - boosting_profile, - efficiency2=boosting_profile, + bus2=nodes + f" {heat_source.intermediate_carrier(heat_system)}", + efficiency=1 + boosting_profile / cop_heat_pump + if heat_source.supports_preheating + else 0, + efficiency2=-boosting_profile + * cop_heat_pump + / (boosting_profile + cop_heat_pump) + if heat_source.supports_preheating + else 1, carrier=f"{heat_system} {heat_source} heat utilisation", p_nom_extendable=True, ) @@ -3325,18 +3337,16 @@ def add_heat( "ptes" ]["discharge_resistive_boosting"] ): - # Inexhaustible sources (air, ground): simple HP → DH heat bus. - # For limited sources the HP is added inside the requires_bus block above. n.add( "Link", nodes, suffix=f" {heat_system} {heat_source} heat pump", - bus0=nodes + f" {heat_system} heat", + bus0=heat_source.hp_output_bus(nodes, heat_system), bus1=nodes, - bus2=heat_source.hp_input_bus(nodes, heat_system), + bus2=heat_source.intermediate_bus(nodes, heat_system), carrier=f"{heat_system} {heat_source} heat pump", - efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(), - efficiency2=1 - 1 / cop_heat_pump.clip(lower=0.001).squeeze(), + efficiency=1 / cop_heat_pump, + efficiency2=heat_source.hp_eff2(cop_heat_pump), capital_cost=costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor, p_min_pu=-(cop_heat_pump > 0).squeeze().astype(float), @@ -3378,11 +3388,13 @@ def add_heat( suffix=f" {heat_system} water pits resistive booster", bus0=nodes + f" {heat_system} heat", bus1=nodes + f" {heat_system} resistive heat", - bus2=ptes_heat_source.hp_input_bus(nodes, heat_system), - efficiency=boosting_profile / (boosting_profile + 1), - efficiency2=1 / (boosting_profile + 1), + bus2=nodes + + f" {ptes_heat_source.intermediate_carrier(heat_system)}", + efficiency=0.5, + efficiency2=0.5, p_nom_extendable=True, - p_min_pu=-(boosting_profile > 0).astype(float), + p_max_pu=0, + p_min_pu=-1, carrier=f"{heat_system} water pits resistive booster", ) From e0f2051a529ef89d2868028e8ef1921554a82619 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 28 Apr 2026 16:34:36 +0200 Subject: [PATCH 16/17] fix efficiencies --- .../ptes_approximator.py | 88 +++-- scripts/prepare_sector_network.py | 327 +++++++++--------- scripts/solve_network.py | 135 +++----- 3 files changed, 258 insertions(+), 292 deletions(-) diff --git a/scripts/build_ptes_operations/ptes_approximator.py b/scripts/build_ptes_operations/ptes_approximator.py index 09e199c39f..f7f96bb8ef 100644 --- a/scripts/build_ptes_operations/ptes_approximator.py +++ b/scripts/build_ptes_operations/ptes_approximator.py @@ -105,11 +105,17 @@ def heat_pump_efficiency(self) -> xr.DataArray: """ return xr.where( - self.layer_temperatures_da > self.return_temperature, - self.rho - * self.c_p - * (self.forward_temperature - self.layer_temperatures_da), - self.rho * self.c_p * (self.forward_temperature - self.return_temperature), + self.layer_temperatures_da > self.forward_temperature, + 0.0, + xr.where( + self.layer_temperatures_da > self.return_temperature, + self.rho + * self.c_p + * (self.forward_temperature - self.layer_temperatures_da), + self.rho + * self.c_p + * (self.forward_temperature - self.return_temperature), + ), ) @property @@ -123,21 +129,17 @@ def layer_needs_boosting(self) -> xr.DataArray: @property def charger_efficiency(self) -> xr.DataArray: - """Heat-to-volume conversion efficiency for PTES charging [m3/MWh].""" - return ( - 1 - / ( - self.rho - * self.c_p - * (self.layer_temperatures_da - self.bottom_temperature) - ) - * self.charging_availability + """Heat-to-volume conversion efficiency for PTES charging [MWh to m3].""" + return xr.where( + self.charging_availability, self.charging_availability / self.m3_to_mwh, 0 ) @property def m3_to_mwh(self) -> np.ndarray: """Per-layer conversion from volume state to energy equivalent [MWh/m3].""" - return self.rho * self.c_p * (self.layer_temperatures - self.bottom_temperature) + return ( + self.rho * self.c_p * (self.layer_temperatures_da - self.bottom_temperature) + ) @property def e_max_pu(self) -> float: @@ -161,6 +163,38 @@ def return_temp_layer(self) -> xr.DataArray: .astype(int) ) + @property + def preheater_return_layer(self) -> xr.DataArray: + """ + Reinjection layer index for the preheater's no-boost return volume. + + For each source layer l and node: + a) T_l > T_ret and l is not the closest-to-T_ret layer: target = closest-to-T_ret + b) T_l > T_ret and l is the closest-to-T_ret layer: target = min(l + 1, num_layers - 1) + (next colder modelled layer; prevents the self-loop that produces phantom heat) + c) T_l <= T_ret: no physical feedback; placeholder = num_layers - 1. + The preheater's efficiency3 (= 1 - layer_needs_boosting) masks this flow + to 0 because any layer colder than T_ret always needs boosting + (T_fwd > T_ret >= T_l). + + Returns dimensions (layer, name). + """ + ret_layer = self.return_temp_layer + next_colder = xr.where( + ret_layer < self.num_layers - 1, ret_layer + 1, self.num_layers - 1 + ) + layer_idx = self.layer_temperatures_da.layer + + return xr.where( + layer_idx < ret_layer, + ret_layer, + xr.where( + layer_idx == ret_layer, + next_colder, + self.num_layers - 1, + ), + ).astype(int) + @property def hp_return_layer(self) -> xr.DataArray: """ @@ -198,15 +232,10 @@ def hp_return_layer(self) -> xr.DataArray: self.layer_temperatures_da - target_temp.sel(layer=l) ).argmin(dim="layer") closest_layer = xr.where( - closest_layer <= l and closest_layer < self.num_layers - 1, + (closest_layer <= l) & (closest_layer < self.num_layers - 1), closest_layer + 1, closest_layer, ) - closest_layer = xr.where( - closest_layer == self.num_layers - 1, # bottom layer - None, - closest_layer.astype(int), - ) ret_val.append(closest_layer) return xr.concat( @@ -214,11 +243,13 @@ def hp_return_layer(self) -> xr.DataArray: ) @property - def interlayer_heat_transfer_coefficients(self) -> np.ndarray: + def interlayer_heat_transfer_coefficients(self) -> xr.DataArray: """.""" - return np.full(self.num_layers - 1, self.interlayer_heat_transfer_coefficient) + return xr.full_like( + self.layer_temperatures_da, self._interlayer_heat_transfer_coefficient() + ) - def _interlayer_heat_transfer_coefficients( + def _interlayer_heat_transfer_coefficient( self, conductivity=0.6, density=1000, @@ -238,6 +269,7 @@ def _interlayer_heat_transfer_coefficients( storage_height: Total height of the storage medium [m] seconds_per_hour: Number of seconds in an hour (for unit conversion) """ + return 0.001 layer_height = storage_height / ( self.num_layers + 1 ) # assuming equal layer heights @@ -274,7 +306,6 @@ def to_dataset(self) -> xr.Dataset: - e_max_pu: scalar """ layer_coord = np.arange(self.num_layers) - pair_coord = np.arange(self.num_layers - 1) ds = xr.Dataset( { @@ -284,6 +315,7 @@ def to_dataset(self) -> xr.Dataset: coords={"layer": layer_coord}, ), "return_temperature_layer": self.return_temp_layer, + "preheater_return_layer": self.preheater_return_layer, "hp_return_layer": self.hp_return_layer, "m3_to_mwh": xr.DataArray( self.m3_to_mwh, @@ -295,11 +327,7 @@ def to_dataset(self) -> xr.Dataset: "layer_needs_boosting": self.layer_needs_boosting, "charger_efficiency": self.charger_efficiency, "charging_availability": self.charging_availability, - "interlayer_heat_transfer_coefficients": xr.DataArray( - self.interlayer_heat_transfer_coefficients, - dims=["layer_pair"], - coords={"layer_pair": pair_coord}, - ), + "interlayer_heat_transfer_coefficients": self.interlayer_heat_transfer_coefficients, "standing_losses": xr.DataArray( self.standing_losses, dims=["layer"], diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 94b70c3600..9b0ce06bcb 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -3126,171 +3126,185 @@ def add_heat( unit="m3", ) - # if not bottom layer - if layer < num_layers - 1: - n.add( - "Bus", - nodes + f" {heat_system} ptes preheater{layer_suffix}", - location=nodes, - carrier=f"{heat_system} ptes preheater", - unit="m3", - ) + n.add( + "Bus", + nodes + f" {heat_system} ptes preheater{layer_suffix}", + location=nodes, + carrier=f"{heat_system} ptes preheater", + unit="m3", + ) - n.add( - "Bus", - nodes + f" {heat_system} ptes hp input{layer_suffix}", - location=nodes, - carrier=f"{heat_system} ptes hp input", - unit="m3", - ) + n.add( + "Bus", + nodes + f" {heat_system} ptes hp input{layer_suffix}", + location=nodes, + carrier=f"{heat_system} ptes hp input", + unit="m3", + ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} water pits charger{layer_suffix}", - bus0=nodes + f" {heat_system} heat", - bus1=nodes + f" {heat_system} water pits{layer_suffix}", - efficiency=charger_efficiency, - carrier=f"{heat_system} water pits charger", - p_nom_extendable=True, - lifetime=costs.at["central water pit storage", "lifetime"], - marginal_cost=costs.at[ - "central water pit charger", "marginal_cost" - ], - p_max_pu=ptes_ds["charging_availability"] - .sel(layer=layer) - .to_pandas(), - ) + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits charger{layer_suffix}", + bus0=nodes + f" {heat_system} heat", + bus1=nodes + f" {heat_system} water pits{layer_suffix}", + efficiency=charger_efficiency, + carrier=f"{heat_system} water pits charger{layer_suffix}", + p_nom_extendable=True, + lifetime=costs.at["central water pit storage", "lifetime"], + marginal_cost=costs.at[ + "central water pit charger", "marginal_cost" + ], + p_max_pu=ptes_ds["charging_availability"] + .sel(layer=layer) + .to_pandas(), + ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} water pits discharger{layer_suffix}", - bus0=nodes + f" {heat_system} water pits{layer_suffix}", - bus1=nodes + f" {heat_system} ptes preheater{layer_suffix}", - carrier=f"{heat_system} water pits discharger", - efficiency=1, - marginal_cost=costs.at[ - "central water pit discharger", "marginal_cost" - ], - p_nom_extendable=True, - lifetime=costs.at["central water pit storage", "lifetime"], - ) + n.add( + "Link", + nodes, + suffix=f" {heat_system} water pits discharger{layer_suffix}", + bus0=nodes + f" {heat_system} water pits{layer_suffix}", + bus1=nodes + f" {heat_system} ptes preheater{layer_suffix}", + carrier=f"{heat_system} water pits discharger{layer_suffix}", + efficiency=1, + marginal_cost=costs.at[ + "central water pit discharger", "marginal_cost" + ], + p_nom_extendable=True, + p_max_pu=0 if layer == num_layers - 1 else 1, + lifetime=costs.at["central water pit storage", "lifetime"], + ) - needs_boosting = ptes_ds["layer_needs_boosting"].sel(layer=layer) + needs_boosting = ptes_ds["layer_needs_boosting"].sel(layer=layer) - return_temperature_layer = pd.Index( - [ - f"{name} {heat_system} water pits layer {ptes_ds['return_temperature_layer'].sel(name=name).values}" - for name in nodes - ] - ) + preheater_return_layer = ( + nodes + + f" {heat_system} water pits layer " + + ptes_ds["preheater_return_layer"] + .sel(layer=layer) + .to_pandas() + .astype(int) + .astype(str) + ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} ptes preheater{layer_suffix}", - bus0=nodes + f" {heat_system} ptes preheater{layer_suffix}", - bus1=nodes + f" {heat_system} heat", - bus2=nodes + f" {heat_system} ptes hp input{layer_suffix}", - bus3=return_temperature_layer, - efficiency=preheater_efficiency.to_pandas(), - efficiency2=needs_boosting.to_pandas(), - efficiency3=(1 - needs_boosting).to_pandas(), - p_nom_extendable=True, - carrier=f"{heat_system} ptes preheater{layer_suffix}", - ) + n.add( + "Link", + nodes, + suffix=f" {heat_system} ptes preheater{layer_suffix}", + bus0=nodes + f" {heat_system} ptes preheater{layer_suffix}", + bus1=nodes + f" {heat_system} heat", + bus2=nodes + f" {heat_system} ptes hp input{layer_suffix}", + bus3=preheater_return_layer, + efficiency=preheater_efficiency.to_pandas(), + efficiency2=needs_boosting.to_pandas(), + efficiency3=(1 - needs_boosting).to_pandas(), + p_nom_extendable=True, + carrier=f"{heat_system} ptes preheater{layer_suffix}", + ) - heat_source = HeatSource(f"ptes{layer_suffix}") - cop_heat_pump = ( - cop.sel( - heat_system=heat_system.system_type.value, - heat_source=heat_source.value, - name=nodes, - ) - .transpose("time", "name") - .to_pandas() - .clip(lower=0.001) + heat_source = HeatSource(f"ptes{layer_suffix}") + cop_heat_pump = ( + cop.sel( + heat_system=heat_system.system_type.value, + heat_source=heat_source.value, + name=nodes, ) + .transpose("time", "name") + .to_pandas() + .clip(lower=0.001) + ) - heat_pump_efficiency = ( - ptes_ds["heat_pump_efficiency"] - .sel(layer=layer) - .transpose("time", "name") - .to_pandas() - .clip(lower=0.001) - ) + heat_pump_efficiency = ( + ptes_ds["heat_pump_efficiency"] + .sel(layer=layer) + .transpose("time", "name") + .to_pandas() + .clip(lower=0.001) + ) - hp_return_layer = pd.Index( - [ - f"{name} {heat_system} water pits layer {ptes_ds['hp_return_layer'].sel(name=name, layer=layer).values.astype(int)}" - for name in nodes - ] - ) + hp_return_layer = ( + nodes + + f" {heat_system} water pits layer " + + ptes_ds["hp_return_layer"] + .sel(layer=layer) + .to_pandas() + .astype(int) + .astype(str) + ) + n.add( + "Link", + nodes, + suffix=f" {heat_system} {heat_source} heat pump", + bus0=nodes + f" {heat_system} heat", + bus1=nodes, + bus2=nodes + f" {heat_system} ptes hp input{layer_suffix}", + bus3=hp_return_layer, + efficiency=1 / cop_heat_pump, + efficiency2=1 / heat_pump_efficiency, + efficiency3=-1 / heat_pump_efficiency, + capital_cost=0, + p_min_pu=-needs_boosting.to_pandas(), + p_max_pu=0, + p_nom_extendable=True, + carrier=f"{heat_system} {heat_source} heat pump", + ) - n.add( - "Link", - nodes, - suffix=f" {heat_system} {heat_source} heat pump", - bus0=nodes + f" {heat_system} heat", - bus1=nodes, - bus2=nodes + f" {heat_system} ptes hp input{layer_suffix}", - bus3=hp_return_layer, - efficiency=1 / cop_heat_pump, - efficiency2=1 / heat_pump_efficiency, - efficiency3=1, - capital_cost=0, - p_min_pu=-needs_boosting.to_pandas(), - p_max_pu=0, - p_nom_extendable=True, - carrier=f"{heat_system} {heat_source} heat pump", - ) + n.links.loc[ + nodes + f" {heat_system} water pits charger{layer_suffix}", + "energy to power ratio", + ] = energy_to_power_ratio_water_pit + n.add( + "Store", + nodes, + suffix=f" {heat_system} water pits{layer_suffix}", + bus=nodes + f" {heat_system} water pits{layer_suffix}", + e_initial=0, + e_nom_extendable=True, + carrier=f"{heat_system} water pits{layer_suffix}", + standing_loss=costs.at[ + "central water pit storage", "standing_losses" + ] + / 100, + ) - n.links.loc[ - nodes + f" {heat_system} water pits charger{layer_suffix}", - "energy to power ratio", - ] = energy_to_power_ratio_water_pit + # Inter-layer links: bidirectional flow between adjacent layers + if layer < num_layers - 1: + lower_layer = f" layer {layer + 1}" n.add( - "Store", + "Link", nodes, - suffix=f" {heat_system} water pits{layer_suffix}", - bus=nodes + f" {heat_system} water pits{layer_suffix}", - # e_cyclic=True, - e_initial=0, - e_nom_extendable=True, - carrier=f"{heat_system} water pits", - standing_loss=costs.at[ - "central water pit storage", "standing_losses" - ] - / 100, - ) - else: - n.add( - "Generator", - nodes + f" {heat_system} ptes bottom layer sink", - bus=nodes + f" {heat_system} water pits{layer_suffix}", - location=nodes, - carrier=f"{heat_system} water pits", - unit="m3", - p_max_pu=0, - p_min_pu=-1, + suffix=f" {heat_system} water pits inter{layer_suffix}-{lower_layer}", + bus0=nodes + f" {heat_system} water pits{layer_suffix}", + bus1=nodes + f" {heat_system} water pits{lower_layer}", + carrier=f"{heat_system} water pits interlayer", p_nom_extendable=True, ) - # Inter-layer links: bidirectional flow between adjacent layers - for pair_idx in range(num_layers - 1): - upper_layer = f"layer {pair_idx}" - lower_layer = f"layer {pair_idx + 1}" - n.add( - "Link", - nodes, - suffix=f" {heat_system} water pits inter {upper_layer}-{lower_layer}", - bus0=nodes + f" {heat_system} water pits {upper_layer}", - bus1=nodes + f" {heat_system} water pits {lower_layer}", - carrier=f"{heat_system} water pits interlayer", - p_nom_extendable=True, + "Generator", + nodes + f" {heat_system} ptes{layer_suffix} vent", + bus=nodes + f" {heat_system} water pits{layer_suffix}", + location=nodes, + carrier=f"{heat_system} water pits", + unit="m3", + p_max_pu=0, p_min_pu=-1, + marginal_cost=-1, + p_nom_extendable=True, + ) + if layer == num_layers - 1: + n.add( + "Generator", + nodes + f" {heat_system} ptes bottom layer water", + bus=nodes + f" {heat_system} water pits{layer_suffix}", + location=nodes, + carrier=f"{heat_system} water pits", + unit="m3", + p_max_pu=1, + p_min_pu=0, + marginal_cost=1, + p_nom_extendable=True, ) if num_layers > 1: @@ -3534,18 +3548,6 @@ def add_heat( "ptes" ]["discharge_resistive_boosting"] ): - cop_heat_pump = ( - cop.sel( - heat_system=heat_system.system_type.value, - heat_source=heat_source.value, - name=nodes, - ) - .transpose("time", "name") - .to_pandas() - if options["time_dep_hp_cop"] - else costs.loc[[costs_name_heat_pump], ["efficiency"]] - ) - p_min_pu = ( 0 if heat_source == HeatSource.PTES @@ -3556,14 +3558,7 @@ def add_heat( else -(cop_heat_pump > 0).astype(float) ) capital_cost = ( - 0 - if heat_source - in [ - HeatSource.PTES_LAYER_0, - HeatSource.PTES_LAYER_1, - HeatSource.PTES_LAYER_2, - ] - else costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor + costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor ) n.add( diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 21720892b7..3f575d1a3f 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1019,21 +1019,22 @@ def add_layered_ptes_volume_capacity_constraint( continue weighted_sum = None + m3_to_mwh_top_layer = ptes_ds["m3_to_mwh"].sel(layer=0).item() for layer_name in layers: - layer_idx = int(layer_name.split("layer")[-1].strip()) - c = float(ptes_ds["m3_to_mwh"].sel(layer=layer_idx).item()) - layer_e = n.model["Store-e"].loc[:, layer_name] - term = c * layer_e - weighted_sum = term if weighted_sum is None else weighted_sum + term + # layer_idx = int(layer_name.split("layer")[-1].strip()) + layer_soc_m3 = n.model["Store-e"].loc[:, layer_name] + weighted_sum = ( + layer_soc_m3 if weighted_sum is None else weighted_sum + layer_soc_m3 + ) - agg_e_nom = n.model["Store-e_nom"].loc[agg] - # Σ W_l · e_{l,t} - ē ≤ 0 - constraints.append(weighted_sum - agg_e_nom) + agg_e_nom_mwh = n.model["Store-e_nom"].loc[agg] + agg_e_nom_m3 = agg_e_nom_mwh / m3_to_mwh_top_layer + constraints.append(weighted_sum - agg_e_nom_m3) merged = linopy.expressions.merge( constraints, dim="Store-ext" if PYPSA_V1 else "name" ) - n.model.add_constraints(merged <= 0, name="layered_ptes_volume_capacity") + n.model.add_constraints(merged == 0, name="layered_ptes_volume_capacity") def add_layered_ptes_heat_pump_capacity_constraint( @@ -1068,66 +1069,6 @@ def add_layered_ptes_heat_pump_capacity_constraint( n.model.add_constraints(merged <= 0, name="layered_ptes_heat_pump_capacity") -def add_layered_ptes_chargers_capacity_constraint( - n: pypsa.Network, -) -> None: - """Tbd""" - layer_chargers = n.links.index[ - n.links.index.str.contains("water pits charger layer") - & n.links.p_nom_extendable - ] - dummy_chargers = n.links.index[ - n.links.index.str.contains("water pits charger") - & ~n.links.index.str.contains("layer") - & n.links.p_nom_extendable - ] - - exprs = [] - for dummy in dummy_chargers: - node_prefix = dummy.replace(" water pits charger", "") - matching = layer_chargers[ - layer_chargers.str.startswith(node_prefix + " water pits charger layer") - ] - layer_sum = sum(n.model["Link-p_nom"].loc[hp] for hp in matching) - dummy_p_nom = n.model["Link-p_nom"].loc[dummy] - exprs.append(layer_sum - dummy_p_nom) - - dim = "Link-ext" if PYPSA_V1 else "name" - merged = linopy.expressions.merge(exprs, dim=dim, cls=type(exprs[0])) - n.model.add_constraints(merged <= 0, name="layered_ptes_chargers_capacity") - - -def add_layered_ptes_dischargers_capacity_constraint( - n: pypsa.Network, -) -> None: - """Tbd""" - layer_dischargers = n.links.index[ - n.links.index.str.contains("water pits discharger layer") - & n.links.p_nom_extendable - ] - dummy_dischargers = n.links.index[ - n.links.index.str.contains("water pits discharger") - & ~n.links.index.str.contains("layer") - & n.links.p_nom_extendable - ] - - exprs = [] - for dummy in dummy_dischargers: - node_prefix = dummy.replace(" water pits discharger", "") - matching = layer_dischargers[ - layer_dischargers.str.startswith( - node_prefix + " water pits discharger layer" - ) - ] - layer_sum = sum(n.model["Link-p_nom"].loc[hp] for hp in matching) - dummy_p_nom = n.model["Link-p_nom"].loc[dummy] - exprs.append(layer_sum - dummy_p_nom) - - dim = "Link-ext" if PYPSA_V1 else "name" - merged = linopy.expressions.merge(exprs, dim=dim, cls=type(exprs[0])) - n.model.add_constraints(merged <= 0, name="layered_ptes_dischargers_capacity") - - def add_layered_ptes_interlayer_flow_constraint( n: pypsa.Network, ptes_ds: xr.Dataset ) -> None: @@ -1140,28 +1081,34 @@ def add_layered_ptes_interlayer_flow_constraint( where Ψ is the interlayer transfer coefficient and e_{l,t} is the state-of-energy of the upper layer store. """ - interlayer_links = n.links.index[ - n.links.carrier.str.contains("water pits interlayer") - ] + interlayer_links = n.links.index[n.links.carrier.str.contains("water pits inter")] + if interlayer_links.empty: + raise RuntimeError( + "No interlayer links found for layered PTES interlayer flow constraints." + ) + elif len(interlayer_links) > len(ptes_ds.layer - 1) * len(ptes_ds.name): + raise RuntimeError( + f"More interlayer links ({len(interlayer_links)}) found than layers in dataset ({len(ptes_ds.layer)})." + ) constraints = [] for link_name in interlayer_links: - # Extract pair index from link name (e.g., "... inter layer 0-layer 1" → 0) - pair_str = link_name.split("inter ")[-1] - pair_idx = int(pair_str.split("layer ")[1].split("-")[0]) + upper_layer_bus = n.links.at[link_name, "bus0"] + # lower_layer_bus = n.links.at[link_name, "bus1"] + + upper_store = n.stores.index[n.stores.bus == upper_layer_bus] + if len(upper_store) == 1: + upper_store = upper_store[0] + else: + raise RuntimeError( + f"None or multiple stores found for upper layer bus {upper_layer_bus} in interlayer link {link_name}." + ) + + layer = int(upper_layer_bus.split("layer")[-1].strip()) heat_transfer_coef = float( - ptes_ds["interlayer_heat_transfer_coefficients"] - .sel(layer_pair=pair_idx) - .item() + ptes_ds["interlayer_heat_transfer_coefficients"].sel(layer=layer).item() ) - upper_bus = n.links.at[link_name, "bus0"] - # Find the store on the upper bus - upper_store = n.stores.index[n.stores.bus == upper_bus] - if upper_store.empty: - continue - upper_store = upper_store[0] - link_p = n.model["Link-p"].loc[:, link_name] store_e = n.model["Store-e"].loc[:, upper_store] constraints.append(link_p - heat_transfer_coef * store_e) @@ -1195,9 +1142,6 @@ def add_layered_ptes_aggregate_throughput_constraint( & n.stores.e_nom_extendable ] - if layer_chargers.empty or agg_stores.empty: - return - constraints = [] for agg in agg_stores: prefix = agg.split("water pits")[0] + "water pits" @@ -1209,20 +1153,19 @@ def add_layered_ptes_aggregate_throughput_constraint( if chargers.empty or dischargers.empty: continue - R = n.links.at[chargers[0], "energy to power ratio"] + e2p_ratio = n.links.at[chargers[0], "energy to power ratio"] weighted_sum = None for ch, dis in zip(chargers, dischargers): layer_idx = int(ch.rsplit("layer ", 1)[-1]) - c = float(ptes_ds["m3_to_mwh"].sel(layer=layer_idx).item()) - term = c * (n.model["Link-p"].loc[:, ch] + n.model["Link-p"].loc[:, dis]) + m3_to_mwh = float(ptes_ds["m3_to_mwh"].sel(layer=layer_idx).item()) + term = ( + n.model["Link-p"].loc[:, ch] + m3_to_mwh * n.model["Link-p"].loc[:, dis] + ) weighted_sum = term if weighted_sum is None else weighted_sum + term agg_e_nom = n.model["Store-e_nom"].loc[agg] - constraints.append(R * weighted_sum - agg_e_nom) - - if not constraints: - return + constraints.append(weighted_sum - e2p_ratio * agg_e_nom) merged = linopy.expressions.merge( constraints, dim="Store-ext" if PYPSA_V1 else "name" @@ -1481,7 +1424,7 @@ def extra_functionality( add_layered_ptes_volume_capacity_constraint(n, ptes_ds) add_layered_ptes_interlayer_flow_constraint(n, ptes_ds) add_layered_ptes_heat_pump_capacity_constraint(n) - # add_layered_ptes_aggregate_throughput_constraint(n, ptes_ds) + add_layered_ptes_aggregate_throughput_constraint(n, ptes_ds) ptes_ds.close() add_battery_constraints(n) From ddad176108703a89254530282c9eb2d37b681120 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 28 Apr 2026 16:34:47 +0200 Subject: [PATCH 17/17] fix incorrekt HP sink inlet temperature --- scripts/build_cop_profiles/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/build_cop_profiles/run.py b/scripts/build_cop_profiles/run.py index 3029691246..d4712c2661 100644 --- a/scripts/build_cop_profiles/run.py +++ b/scripts/build_cop_profiles/run.py @@ -215,7 +215,7 @@ def get_source_inlet_temperature( # When source temperature <= return temperature, no preheating: # heat pump draws directly from the source. return xr.where( - source_temperature > central_heating_return_temperature, + source_temperature < central_heating_return_temperature, source_temperature, central_heating_return_temperature, )