"""
Pandapipes Network Initialization Module
=========================================
This module provides comprehensive network initialization capabilities for district heating
systems using GeoJSON-based geographic data.
:author: Dipl.-Ing. (FH) Jonas Pfeiffer
It handles the complete workflow from GeoJSON
data processing to pandapipes network creation, including heat demand integration, temperature
calculation, and multi-producer system configuration.
The module supports various network configurations including traditional hot water networks,
cold networks with decentralized heat pumps, and hybrid systems with multiple heat generators.
It automatically processes building heat demands, calculates temperature requirements, and
creates appropriate network topologies with proper controller configurations.
"""
import json
import logging
import warnings
from typing import Any
import geopandas as gpd
import numpy as np
import pandapipes as pp
import pandas as pd
from pandapipes.control.run_control import run_control
from districtheatingsim.constants import CP_WATER_KJ_KGK, KELVIN_OFFSET
from districtheatingsim.net_generation.network_geojson_schema import NetworkGeoJSONSchema
from districtheatingsim.net_simulation_pandapipes.pipe_std_types import resolve_pipe_u_w_per_m2k
from districtheatingsim.net_simulation_pandapipes.result_validation import (
validate_net_results,
validate_pressure_plausibility,
)
from districtheatingsim.net_simulation_pandapipes.utilities import (
COP_WP,
correct_flow_directions,
create_controllers,
init_diameter_types,
)
[docs]
def initialize_geojson(NetworkGenerationData) -> Any:
"""
Initialize district heating network from unified GeoJSON and heat demand data.
:param NetworkGenerationData: Configuration with network_geojson_path, heat_demand_json_path, temperatures, pipe specs, producer config
:type NetworkGenerationData: object
:return: Updated NetworkGenerationData with initialized net, time series, mass flows, building data
:rtype: Any
:raises FileNotFoundError: If GeoJSON or JSON files not found
:raises ValueError: If temperature constraints violated (return > supply)
:raises KeyError: If required JSON fields missing
.. note::
Loads unified GeoJSON (Wärmenetz.geojson), processes heat demands, validates temperatures.
Handles cold networks (COP calculation), applies 2% min load. Calculates mass flows
(main: total demand/ΔT, secondary: percentage-based). Creates complete pandapipes network.
"""
# Load unified network GeoJSON data
# Read unified GeoJSON file
network_gdf = gpd.read_file(NetworkGenerationData.network_geojson_path, driver="GeoJSON")
# Separate features by type
gdf_dict = {
"flow_line": network_gdf[network_gdf["feature_type"] == NetworkGeoJSONSchema.FEATURE_TYPE_FLOW].copy(),
"return_line": network_gdf[network_gdf["feature_type"] == NetworkGeoJSONSchema.FEATURE_TYPE_RETURN].copy(),
"heat_consumer": network_gdf[network_gdf["feature_type"] == NetworkGeoJSONSchema.FEATURE_TYPE_BUILDING].copy(),
"heat_producer": network_gdf[network_gdf["feature_type"] == NetworkGeoJSONSchema.FEATURE_TYPE_GENERATOR].copy(),
}
print(f"Loaded unified network GeoJSON with {len(network_gdf)} features")
print(f" Flow lines: {len(gdf_dict['flow_line'])}")
print(f" Return lines: {len(gdf_dict['return_line'])}")
print(f" Heat consumers: {len(gdf_dict['heat_consumer'])}")
print(f" Heat producers: {len(gdf_dict['heat_producer'])}")
print(f"Max supply temperature heat generator: {NetworkGenerationData.max_supply_temperature_heat_generator} °C")
# Load and process heat demand data
with open(NetworkGenerationData.heat_demand_json_path, encoding="utf-8") as f:
loaded_data = json.load(f)
results = {k: v for k, v in loaded_data.items() if isinstance(v, dict) and "wärme" in v}
heat_demand_df = pd.DataFrame.from_dict({k: v for k, v in loaded_data.items() if k.isdigit()}, orient="index")
# Extract building temperature data
supply_temperature_buildings = heat_demand_df["VLT_max"].values.astype(float)
return_temperature_buildings = heat_demand_df["RLT_max"].values.astype(float)
# Extract time series data
yearly_time_steps = np.array(heat_demand_df["zeitschritte"].values[0]).astype(np.datetime64)
total_building_heat_demand_W = np.array([results[str(i)]["wärme"] for i in range(len(results))]) * 1000
total_building_heating_demand_W = np.array([results[str(i)]["heizwärme"] for i in range(len(results))]) * 1000
total_building_hot_water_demand_W = (
np.array([results[str(i)]["warmwasserwärme"] for i in range(len(results))]) * 1000
)
supply_temperature_building_curve = np.array([results[str(i)]["vorlauftemperatur"] for i in range(len(results))])
return_temperature_building_curve = np.array([results[str(i)]["rücklauftemperatur"] for i in range(len(results))])
maximum_building_heat_load_W = np.array(results["0"]["max_last"]) * 1000
print(f"Max heat demand buildings (W): {maximum_building_heat_load_W}")
# Calculate return temperature for heat consumers
if NetworkGenerationData.fixed_return_temperature_heat_consumer is None:
return_temperature_heat_consumer = return_temperature_buildings + NetworkGenerationData.dT_RL
print(f"Return temperature heat consumers: {return_temperature_heat_consumer} °C")
else:
return_temperature_heat_consumer = np.full_like(
return_temperature_buildings, NetworkGenerationData.fixed_return_temperature_heat_consumer
)
print(f"Return temperature heat consumers: {return_temperature_heat_consumer} °C")
# Validate temperature constraints
if np.any(return_temperature_heat_consumer >= NetworkGenerationData.max_supply_temperature_heat_generator):
raise ValueError(
"Return temperature must not be higher than the supply temperature at the injection point. Please check your inputs."
)
# Calculate minimum supply temperature for heat consumers
if NetworkGenerationData.min_supply_temperature_building is None:
min_supply_temperature_heat_consumer = np.zeros_like(
supply_temperature_buildings, NetworkGenerationData.min_supply_temperature_building
)
print(f"Minimum supply temperature heat consumers: {min_supply_temperature_heat_consumer} °C")
else:
min_supply_temperature_heat_consumer = np.full_like(
supply_temperature_buildings,
NetworkGenerationData.min_supply_temperature_building + NetworkGenerationData.dT_RL,
)
print(f"Minimum supply temperature heat consumers: {min_supply_temperature_heat_consumer} °C")
# Validate minimum supply temperature constraints
if np.any(min_supply_temperature_heat_consumer >= NetworkGenerationData.max_supply_temperature_heat_generator):
raise ValueError(
"Supply temperature at the heat consumer cannot be higher than the supply temperature at the injection point. Please check your inputs."
)
# Initialize heat and power arrays
waerme_hast_ges_W = []
max_waerme_hast_ges_W = []
strombedarf_hast_ges_W = []
max_el_leistung_hast_ges_W = []
# Process heat demands based on network configuration
if NetworkGenerationData.netconfiguration == "kaltes Netz":
# Cold network: Calculate heat pump performance
COP_file_values = np.genfromtxt(NetworkGenerationData.COP_filename, delimiter=";")
COP, _ = COP_WP(supply_temperature_buildings, return_temperature_heat_consumer, COP_file_values)
print(f"COP dezentrale Wärmepumpen Gebäude: {COP}")
# Calculate heat pump electricity consumption and network heat demand
for waerme_gebaeude, leistung_gebaeude, cop in zip(
total_building_heat_demand_W, maximum_building_heat_load_W, COP, strict=False
):
strombedarf_wp = waerme_gebaeude / cop
waerme_hast = waerme_gebaeude - strombedarf_wp
waerme_hast_ges_W.append(waerme_hast)
strombedarf_hast_ges_W.append(strombedarf_wp)
el_leistung_wp = leistung_gebaeude / cop
waerme_leistung_hast = leistung_gebaeude - el_leistung_wp
max_waerme_hast_ges_W.append(waerme_leistung_hast)
max_el_leistung_hast_ges_W.append(el_leistung_wp)
waerme_hast_ges_W = np.array(waerme_hast_ges_W)
max_waerme_hast_ges_W = np.array(max_waerme_hast_ges_W)
strombedarf_hast_ges_W = np.array(strombedarf_hast_ges_W)
max_el_leistung_hast_ges_W = np.array(max_el_leistung_hast_ges_W)
else:
# Traditional network: Direct heat transfer
waerme_hast_ges_W = total_building_heat_demand_W
max_waerme_hast_ges_W = maximum_building_heat_load_W
strombedarf_hast_ges_W = np.zeros_like(total_building_heat_demand_W)
max_el_leistung_hast_ges_W = np.zeros_like(maximum_building_heat_load_W)
# Prepare data dictionaries for network creation
consumer_dict = {
"qext_w": max_waerme_hast_ges_W,
"min_supply_temperature_heat_consumer": min_supply_temperature_heat_consumer,
"return_temperature_heat_consumer": return_temperature_heat_consumer,
}
pipe_dict = {
"pipetype": NetworkGenerationData.pipetype,
"v_max_pipe": NetworkGenerationData.max_velocity_pipe,
"material_filter": NetworkGenerationData.material_filter_pipe,
"pipe_creation_mode": "type",
"k_mm": NetworkGenerationData.k_mm_pipe,
}
# Calculate mass flows for secondary producers
if NetworkGenerationData.secondary_producers:
cp = CP_WATER_KJ_KGK # kJ/kgK - specific heat capacity of water
print(f"Specific heat capacity of water: {cp} kJ/kgK")
print(f"maximum_building_heat_load_W: {maximum_building_heat_load_W}")
sum_maximum_building_heat_load_W = np.sum(maximum_building_heat_load_W)
print(f"sum_maximum_building_heat_load_W: {sum_maximum_building_heat_load_W}")
print(
f"Max supply temperature heat generator: {NetworkGenerationData.max_supply_temperature_heat_generator} °C"
)
print(f"Return temperature heat consumer: {np.average(return_temperature_heat_consumer)} °C")
mass_flow = (sum_maximum_building_heat_load_W / 1000) / (
cp
* (
NetworkGenerationData.max_supply_temperature_heat_generator
- np.average(return_temperature_heat_consumer)
)
)
print(f"Mass flow of main producer: {mass_flow} kg/s")
for secondary_producer in NetworkGenerationData.secondary_producers:
secondary_producer.mass_flow = secondary_producer.load_percentage / 100 * mass_flow
print(f"Mass flow of secondary producer {secondary_producer.index}: {secondary_producer.mass_flow} kg/s")
producer_dict = {
"supply_temperature": NetworkGenerationData.max_supply_temperature_heat_generator,
"flow_pressure_pump": NetworkGenerationData.flow_pressure_pump,
"lift_pressure_pump": NetworkGenerationData.lift_pressure_pump,
"main_producer_location_index": NetworkGenerationData.main_producer_location_index,
"secondary_producers": NetworkGenerationData.secondary_producers,
}
# Create the pandapipes network
net = create_network(gdf_dict, consumer_dict, pipe_dict, producer_dict)
# Store processed data in NetworkGenerationData object
NetworkGenerationData.supply_temperature_buildings = supply_temperature_buildings
NetworkGenerationData.return_temperature_buildings = return_temperature_buildings
NetworkGenerationData.supply_temperature_building_curve = supply_temperature_building_curve
NetworkGenerationData.return_temperature_building_curve = return_temperature_building_curve
NetworkGenerationData.yearly_time_steps = yearly_time_steps
NetworkGenerationData.waerme_gebaeude_ges_W = total_building_heat_demand_W
NetworkGenerationData.heizwaerme_gebaeude_ges_W = total_building_heating_demand_W
NetworkGenerationData.ww_waerme_gebaeude_ges_W = total_building_hot_water_demand_W
NetworkGenerationData.max_waerme_gebaeude_ges_W = maximum_building_heat_load_W
NetworkGenerationData.return_temperature_heat_consumer = return_temperature_heat_consumer
NetworkGenerationData.min_supply_temperature_heat_consumer = min_supply_temperature_heat_consumer
NetworkGenerationData.waerme_hast_ges_W = waerme_hast_ges_W
NetworkGenerationData.strombedarf_hast_ges_W = strombedarf_hast_ges_W
NetworkGenerationData.max_waerme_hast_ges_W = max_waerme_hast_ges_W
NetworkGenerationData.max_el_leistung_hast_ges_W = max_el_leistung_hast_ges_W
# Apply minimum load constraints (2% of maximum)
max_heat = np.max(NetworkGenerationData.waerme_hast_ges_W)
max_power = np.max(NetworkGenerationData.strombedarf_hast_ges_W)
NetworkGenerationData.waerme_hast_ges_W = np.where(
NetworkGenerationData.waerme_hast_ges_W < 0.02 * max_heat,
0.02 * max_heat,
NetworkGenerationData.waerme_hast_ges_W,
)
NetworkGenerationData.strombedarf_hast_ges_W = np.where(
NetworkGenerationData.waerme_hast_ges_W < 0.02 * max_heat,
0.02 * max_power,
NetworkGenerationData.strombedarf_hast_ges_W,
)
# Convert to kW units
NetworkGenerationData.waerme_hast_ges_kW = np.where(
NetworkGenerationData.waerme_hast_ges_W == 0, 0, NetworkGenerationData.waerme_hast_ges_W / 1000
)
NetworkGenerationData.strombedarf_hast_ges_kW = np.where(
NetworkGenerationData.strombedarf_hast_ges_W == 0, 0, NetworkGenerationData.strombedarf_hast_ges_W / 1000
)
# Calculate total system demands
NetworkGenerationData.waerme_ges_kW = np.sum(NetworkGenerationData.waerme_hast_ges_kW, axis=0)
NetworkGenerationData.strombedarf_ges_kW = np.sum(NetworkGenerationData.strombedarf_hast_ges_kW, axis=0)
NetworkGenerationData.net = net
return NetworkGenerationData
[docs]
def get_line_coords_and_lengths(gdf: gpd.GeoDataFrame) -> tuple[list[list[tuple]], list[float]]:
"""
Extract 2-D coordinates and lengths from LineString geometries.
Z-coordinates (elevation) are stripped here so that the returned
coordinate tuples are always ``(x, y)``. Elevation data is extracted
separately via :func:`build_elevation_lookup_from_gdf` and passed to
the junction-creation step as ``height_m``.
:param gdf: GeoDataFrame with LineString geometries (flow/return lines)
:type gdf: gpd.GeoDataFrame
:return: (all_line_coords, all_line_lengths) - 2-D coordinate sequences and lengths
:rtype: Tuple[List[List[Tuple]], List[float]]
.. note::
Only processes LineString geometries, skips others with warning. Uses
GeoPandas length property for geodetic calculation.
"""
all_line_coords, all_line_lengths = [], []
gdf["length"] = gdf.geometry.length
for _index, row in gdf.iterrows():
line = row["geometry"]
if line.geom_type == "LineString":
# Strip Z — keep only (x, y) for junction dict keys
coords = [(c[0], c[1]) for c in line.coords]
length = row["length"]
all_line_coords.append(coords)
all_line_lengths.append(length)
else:
print(f"Geometrie ist kein LineString: {line.geom_type}")
return all_line_coords, all_line_lengths
[docs]
def build_elevation_lookup_from_gdf(gdf: gpd.GeoDataFrame) -> dict[tuple, float]:
"""
Build a ``{(x, y): z_m}`` elevation lookup from a GeoDataFrame with 3-D geometries.
Only vertices that actually carry a Z-coordinate contribute to the lookup.
Vertices without Z are silently skipped (they will fall back to ``height_m=0``).
:param gdf: GeoDataFrame with Point or LineString geometries (may be 2-D or 3-D)
:type gdf: gpd.GeoDataFrame
:return: Mapping from 2-D coordinate tuple to elevation [m above sea level]
:rtype: Dict[Tuple[float, float], float]
"""
lookup: dict[tuple, float] = {}
for geom in gdf.geometry:
if geom is None:
continue
if geom.geom_type == "Point":
if geom.has_z:
lookup[(geom.x, geom.y)] = geom.z
elif geom.geom_type == "LineString":
for coord in geom.coords:
if len(coord) > 2:
lookup[(coord[0], coord[1])] = coord[2]
return lookup
[docs]
def get_all_point_coords_from_line_cords(all_line_coords: list[list[tuple]]) -> list[tuple]:
"""
Extract unique point coordinates for network junction creation.
:param all_line_coords: List of 2-D coordinate sequences [(x1,y1), (x2,y2), ...]
:type all_line_coords: List[List[Tuple]]
:return: Unique point coordinates (x, y) for junction locations
:rtype: List[Tuple]
.. note::
Removes duplicates using set operations. Order not guaranteed.
Essential for proper network topology without duplicate junctions.
"""
point_coords = [koordinate for paar in all_line_coords for koordinate in paar]
unique_point_coords = list(set(point_coords))
return unique_point_coords
[docs]
def create_network(
gdf_dict: dict[str, gpd.GeoDataFrame],
consumer_dict: dict[str, Any],
pipe_dict: dict[str, Any],
producer_dict: dict[str, Any],
) -> pp.pandapipesNet:
"""
Create complete pandapipes network with junctions, pipes, consumers, and producers.
:param gdf_dict: GeoDataFrames with keys flow_line, return_line, heat_consumer, heat_producer
:type gdf_dict: Dict[str, gpd.GeoDataFrame]
:param consumer_dict: Heat consumer config (qext_w, min_supply_temperature_heat_consumer, return_temperature_heat_consumer)
:type consumer_dict: Dict[str, Any]
:param pipe_dict: Pipe config (pipetype, v_max_pipe, material_filter, pipe_creation_mode, k_mm)
:type pipe_dict: Dict[str, Any]
:param producer_dict: Producer config (supply_temperature, pressures, main_producer_location_index, secondary_producers)
:type producer_dict: Dict[str, Any]
:return: Complete pandapipes network with optimized diameters and controllers
:rtype: pp.pandapipesNet
.. note::
Steps: 1) junctions from coords, 2) pipes (supply/return), 3) heat consumers,
4) producers (main=circ_pump_pressure, secondary=circ_pump_mass), 5) pipeflow,
6) controllers, diameter optimization. Corrects flow directions automatically.
"""
# Extract data from dictionaries
gdf_flow_line, gdf_return_line, gdf_heat_exchanger, gdf_heat_producer = (
gdf_dict["flow_line"],
gdf_dict["return_line"],
gdf_dict["heat_consumer"],
gdf_dict["heat_producer"],
)
qext_w, min_supply_temperature_heat_consumer, return_temperature_heat_consumer = (
consumer_dict["qext_w"],
consumer_dict["min_supply_temperature_heat_consumer"],
consumer_dict["return_temperature_heat_consumer"],
)
supply_temperature, flow_pressure_pump, lift_pressure_pump, main_producer_location_index, secondary_producers = (
producer_dict["supply_temperature"],
producer_dict["flow_pressure_pump"],
producer_dict["lift_pressure_pump"],
producer_dict["main_producer_location_index"],
producer_dict["secondary_producers"],
)
pipetype, v_max_pipe, material_filter, pipe_creation_mode, k_mm = (
pipe_dict["pipetype"],
pipe_dict["v_max_pipe"],
pipe_dict["material_filter"],
pipe_dict["pipe_creation_mode"],
pipe_dict["k_mm"],
)
# Create empty network and get pipe properties
net = pp.create_empty_network(fluid="water")
pipe_std_types = pp.std_types.available_std_types(net, "pipe")
properties = pipe_std_types.loc[pipetype]
diameter_mm = properties["inner_diameter_mm"]
u_w_per_m2k_pipe = resolve_pipe_u_w_per_m2k(properties)
# Convert temperatures to Kelvin
supply_temperature_k = supply_temperature + KELVIN_OFFSET
return_temperature_heat_consumer_k = return_temperature_heat_consumer + KELVIN_OFFSET
def create_junctions_from_coords(
net_i: pp.pandapipesNet, all_coords: list[tuple], elevation_lookup: dict[tuple, float]
) -> dict[tuple, int]:
"""
Create junctions in the network from coordinate points.
Each junction receives a ``height_m`` value looked up from
*elevation_lookup*. This enables pandapipes to include the
hydrostatic pressure contribution ``ρ·g·Δh`` in the flow equations.
Junctions without an entry in the lookup default to ``height_m=0``.
Parameters
----------
net_i : pp.pandapipesNet
The pandapipes network object.
all_coords : List[Tuple]
List of 2-D coordinate tuples for junction locations.
elevation_lookup : Dict[Tuple[float, float], float]
Mapping from ``(x, y)`` to elevation [m above sea level].
Returns
-------
Dict[Tuple, int]
Dictionary mapping coordinates to junction IDs.
"""
junction_dict = {}
for i, coords in enumerate(all_coords, start=0):
height = elevation_lookup.get(coords, 0.0)
junction_id = pp.create_junction(
net_i, pn_bar=1.05, tfluid_k=supply_temperature_k, height_m=height, name=f"Junction {i}", geodata=coords
)
junction_dict[coords] = junction_id
return junction_dict
def create_pipes(
net_i: pp.pandapipesNet,
all_line_coords: list[list[tuple]],
all_line_lengths: list[float],
junction_dict: dict[tuple, int],
pipe_mode: str,
pipe_type_or_diameter: str | float,
line_type: str,
) -> None:
"""
Create pipes in the network from line geometries.
Parameters
----------
net_i : pp.pandapipesNet
The pandapipes network object.
all_line_coords : List[List[Tuple]]
List of line coordinate sequences.
all_line_lengths : List[float]
List of corresponding line lengths.
junction_dict : Dict[Tuple, int]
Dictionary mapping coordinates to junction IDs.
pipe_mode : str
Pipe creation mode ("type" or "diameter").
pipe_type_or_diameter : Union[str, float]
Pipe type name or diameter value.
line_type : str
Description of line type for naming.
"""
for coords, length_m, i in zip(all_line_coords, all_line_lengths, range(len(all_line_coords)), strict=False):
if pipe_mode == "diameter":
diameter_mm = pipe_type_or_diameter
pp.create_pipe_from_parameters(
net_i,
from_junction=junction_dict[coords[0]],
to_junction=junction_dict[coords[1]],
length_km=length_m / 1000,
diameter_m=diameter_mm / 1000,
k_mm=k_mm,
u_w_per_m2k=u_w_per_m2k_pipe,
name=f"{line_type} {i}",
geodata=coords,
sections=5,
text_k=283,
)
elif pipe_mode == "type":
pp.create_pipe(
net_i,
from_junction=junction_dict[coords[0]],
to_junction=junction_dict[coords[1]],
std_type=pipe_type_or_diameter,
length_km=length_m / 1000,
k_mm=k_mm,
name=f"{line_type} {i}",
geodata=coords,
sections=5,
text_k=283,
)
def create_heat_consumers(
net_i: pp.pandapipesNet, all_coords: list[list[tuple]], junction_dict: dict[tuple, int], name_prefix: str
) -> None:
"""Create heat consumers in the network."""
for i, (coords, q, t) in enumerate(zip(all_coords, qext_w, return_temperature_heat_consumer_k, strict=False)):
pp.create_heat_consumer(
net_i,
from_junction=junction_dict[coords[0]],
to_junction=junction_dict[coords[1]],
loss_coefficient=0,
qext_w=q,
treturn_k=t,
name=f"{name_prefix} {i}",
)
def _resolve_pump_junctions(coords, jd_vl, jd_rl):
"""Return (return_junction_idx, flow_junction_idx) for a generator connection line.
A generator-connection line has one endpoint on the Vorlauf (supply) network
and one on the Rücklauf (return) network. We look up each coordinate in both
dicts so the assignment is correct regardless of the line's digitisation direction.
:param coords: [coord0, coord1] of the generator-connection line.
:param jd_vl: Junction dict for the Vorlauf (flow/supply) network.
:param jd_rl: Junction dict for the Rücklauf (return) network.
:returns: (return_junction_idx, flow_junction_idx)
:raises KeyError: if neither coordinate can be found in the expected dicts.
"""
if coords[0] in jd_vl:
return jd_rl[coords[1]], jd_vl[coords[0]] # return=RL end, flow=VL end
elif coords[1] in jd_vl:
return jd_rl[coords[0]], jd_vl[coords[1]] # return=RL end, flow=VL end
else:
# Fallback – merged dict, original (possibly wrong) order
logging.warning(
"Could not determine VL/RL side for generator connection; "
"using original coord order (may cause pump-direction warning)"
)
merged = {**jd_vl, **jd_rl}
return merged[coords[1]], merged[coords[0]]
def create_circulation_pump_pressure(
net_i: pp.pandapipesNet,
all_coords: list[list[tuple]],
jd_vl: dict[tuple, int],
jd_rl: dict[tuple, int],
name_prefix: str,
) -> None:
"""Create pressure-controlled circulation pumps with correct VL/RL junction assignment."""
for i, coords in enumerate(all_coords, start=0):
return_junc, flow_junc = _resolve_pump_junctions(coords, jd_vl, jd_rl)
pp.create_circ_pump_const_pressure(
net_i,
return_junc,
flow_junc,
p_flow_bar=flow_pressure_pump,
plift_bar=lift_pressure_pump,
t_flow_k=supply_temperature_k,
type="auto",
name=f"{name_prefix} {i}",
)
def create_circulation_pump_mass_flow(
net_i: pp.pandapipesNet,
all_coords: list[list[tuple]],
jd_vl: dict[tuple, int],
jd_rl: dict[tuple, int],
name_prefix: str,
mass_flows: list[float],
elevation_lookup: dict[tuple, float],
) -> None:
"""Create mass-flow-controlled circulation pumps with correct VL/RL junction assignment.
The intermediate junction inserted between pump and flow-control element
receives an elevation equal to the average of its two endpoint elevations.
"""
for i, (coords, mass_flow) in enumerate(zip(all_coords, mass_flows, strict=False), start=0):
return_junc, flow_junc = _resolve_pump_junctions(coords, jd_vl, jd_rl)
mid_coord = ((coords[0][0] + coords[1][0]) / 2, (coords[0][1] + coords[1][1]) / 2)
# Interpolate elevation for the synthetic mid-point junction
z0 = elevation_lookup.get(coords[0], 0.0)
z1 = elevation_lookup.get(coords[1], 0.0)
mid_height = (z0 + z1) / 2.0
mid_junction_idx = pp.create_junction(
net_i,
pn_bar=1.05,
tfluid_k=supply_temperature_k,
height_m=mid_height,
name=f"Junction {name_prefix}",
geodata=mid_coord,
)
pp.create_circ_pump_const_mass_flow(
net_i,
return_junc,
mid_junction_idx,
p_flow_bar=flow_pressure_pump,
mdot_flow_kg_per_s=mass_flow,
t_flow_k=supply_temperature_k,
type="auto",
name=f"{name_prefix} {i}",
in_service=True,
)
pp.create_flow_control(net_i, mid_junction_idx, flow_junc, controlled_mdot_kg_per_s=mass_flow)
# Build elevation lookup from 3-D GeoJSON geometries (Z-coord = height above NN).
# All four GDFs are merged so that every junction—whether on a network line,
# building connection, or generator connection—receives a correct height_m.
elevation_lookup: dict[tuple, float] = {}
for gdf_src in (gdf_flow_line, gdf_return_line, gdf_heat_exchanger, gdf_heat_producer):
elevation_lookup.update(build_elevation_lookup_from_gdf(gdf_src))
if elevation_lookup:
z_vals = list(elevation_lookup.values())
dh = max(z_vals) - min(z_vals)
logging.info(
"Elevation data found: min=%.1f m, max=%.1f m, Δh=%.1f m → hydrostatic pressure offset ≈ %.2f bar",
min(z_vals),
max(z_vals),
dh,
1000 * 9.81 * dh / 1e5,
)
if dh > 0.5 * lift_pressure_pump * 1e5 / (1000 * 9.81):
logging.warning(
"Hydrostatic head (%.1f m) exceeds 50%% of pump lift (%.2f bar). Verify pressure zone design.",
dh,
lift_pressure_pump,
)
else:
logging.info("No elevation data in GeoJSON — all junctions set to height_m=0.")
# Create network topology
flow_line_2d_coords = get_line_coords_and_lengths(gdf_flow_line)[0]
return_line_2d_coords = get_line_coords_and_lengths(gdf_return_line)[0]
junction_dict_vl = create_junctions_from_coords(
net, get_all_point_coords_from_line_cords(flow_line_2d_coords), elevation_lookup
)
junction_dict_rl = create_junctions_from_coords(
net, get_all_point_coords_from_line_cords(return_line_2d_coords), elevation_lookup
)
# Create pipes
create_pipes(
net,
*get_line_coords_and_lengths(gdf_flow_line),
junction_dict_vl,
pipe_creation_mode,
diameter_mm if pipe_creation_mode == "diameter" else pipetype,
"flow line",
)
create_pipes(
net,
*get_line_coords_and_lengths(gdf_return_line),
junction_dict_rl,
pipe_creation_mode,
diameter_mm if pipe_creation_mode == "diameter" else pipetype,
"return line",
)
# Create heat consumers
create_heat_consumers(
net,
get_line_coords_and_lengths(gdf_heat_exchanger)[0],
{**junction_dict_vl, **junction_dict_rl},
"heat consumer",
)
# Create heat producers
all_heat_producer_coords, all_heat_producer_lengths = get_line_coords_and_lengths(gdf_heat_producer)
if all_heat_producer_coords:
# Main producer (pressure controlled)
create_circulation_pump_pressure(
net,
[all_heat_producer_coords[main_producer_location_index]],
junction_dict_vl,
junction_dict_rl,
"heat source",
)
# Secondary producers (mass flow controlled)
if secondary_producers:
mass_flows = [producer.mass_flow for producer in secondary_producers]
secondary_coords = [all_heat_producer_coords[producer.index] for producer in secondary_producers]
create_circulation_pump_mass_flow(
net,
secondary_coords,
junction_dict_vl,
junction_dict_rl,
"heat source slave",
mass_flows,
elevation_lookup,
)
print(f"secondary_producers: {secondary_producers}")
# Initial flow simulation – catch the pump-direction UserWarning so that
# correct_flow_directions() below can fix the topology.
with warnings.catch_warnings():
warnings.simplefilter("ignore", UserWarning)
try:
pp.pipeflow(net, mode="bidirectional", iter=100)
except UserWarning as e:
logging.warning(f"Initial pipeflow UserWarning (will be corrected): {e}")
# Network optimization
net = create_controllers(
net,
qext_w,
supply_temperature,
min_supply_temperature_heat_consumer,
return_temperature_heat_consumer,
secondary_producers,
)
try:
run_control(net, mode="bidirectional", iter=100)
except UserWarning as e:
logging.warning(f"run_control UserWarning (will be corrected by correct_flow_directions): {e}")
net = correct_flow_directions(net)
net = init_diameter_types(net, v_max_pipe=v_max_pipe, material_filter=material_filter, k=k_mm)
# Fail loudly at build time if the design state did not converge, instead of
# letting NaN propagate into the time series (BACKLOG C2).
validate_net_results(net, context="network generation")
# Soft check: warn (don't raise) on physically-impossible negative pressures —
# a sign the pump head is too low for the network losses (BACKLOG C14).
validate_pressure_plausibility(net, context="network generation")
return net