"""
Heat demand profile generation from CSV building data.
Integrates VDI 4655 and BDEW calculation methods for batch processing
of building portfolios with temperature curves for district heating design.
:author: Dipl.-Ing. (FH) Jonas Pfeiffer
"""
import numpy as np
import pandas as pd
from pyslpheat import bdew_calculate, vdi4655_calculate
from districtheatingsim.utilities.csv_schemas import validate_csv_columns
def _easter_sunday(year: int) -> pd.Timestamp:
"""
Compute Easter Sunday for a given year using the Anonymous Gregorian algorithm.
:param year: Four-digit year
:type year: int
:return: Easter Sunday date
:rtype: pd.Timestamp
"""
a = year % 19
b, c = divmod(year, 100)
d, e = divmod(b, 4)
f = (b + 8) // 25
g = (b - f + 1) // 3
h = (19 * a + b - d - g + 15) % 30
i, k = divmod(c, 4)
ll = (32 + 2 * e + 2 * i - h - k) % 7
m = (a + 11 * h + 22 * ll) // 451
month = (h + ll - 7 * m + 114) // 31
day = ((h + ll - 7 * m + 114) % 31) + 1
return pd.Timestamp(year=year, month=month, day=day)
def _german_national_holidays(year: int) -> np.ndarray:
"""
Return German national public holidays for *year* as a ``datetime64[D]`` array.
Includes only holidays that apply across all Bundesländer.
:param year: Four-digit year
:type year: int
:return: Array of holiday dates
:rtype: numpy.ndarray
"""
easter = _easter_sunday(year)
dates = [
pd.Timestamp(year, 1, 1), # Neujahr
easter - pd.Timedelta(days=2), # Karfreitag
easter + pd.Timedelta(days=1), # Ostermontag
pd.Timestamp(year, 5, 1), # Tag der Arbeit
easter + pd.Timedelta(days=39), # Christi Himmelfahrt
easter + pd.Timedelta(days=50), # Pfingstmontag
pd.Timestamp(year, 10, 3), # Tag der Deutschen Einheit
pd.Timestamp(year, 12, 25), # 1. Weihnachtstag
pd.Timestamp(year, 12, 26), # 2. Weihnachtstag
]
return np.array([d.date() for d in dates], dtype="datetime64[D]")
[docs]
def generate_profiles_from_csv(
data: pd.DataFrame, TRY: str, calc_method: str, year: int = 2023
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Generate heat demand profiles from CSV building data.
:param data: Building data with required columns (Wärmebedarf, Gebäudetyp, Subtyp,
WW_Anteil, Normaußentemperatur, VLT_max, RLT_max, Steigung_Heizkurve)
and optional BDEW columns (Heizgrenztemperatur, Heizexponent, P_max).
:type data: pd.DataFrame
:param TRY: Path to Test Reference Year weather data file
:type TRY: str
:param calc_method: Calculation method ('Datensatz', 'VDI4655', or 'BDEW')
:type calc_method: str
:param year: Year for profile calculation (affects weekday/holiday pattern and TRY mapping),
defaults to 2023
:type year: int
:return: Tuple of (time_steps, total_heat_W, heating_heat_W, warmwater_heat_W, max_heat_W,
supply_temp, return_temp, air_temp)
:rtype: Tuple[np.ndarray, ...]
:raises KeyError: If required CSV columns are missing
:raises ValueError: If data types are invalid
:raises FileNotFoundError: If TRY file not found
.. note::
'Datensatz' mode auto-selects VDI4655 for residential (EFH/MFH), BDEW for commercial
buildings. Optional BDEW columns (Heizgrenztemperatur, Heizexponent, P_max) are read
per building if present; missing values fall back to pyslpheat defaults.
"""
# Fail up front with one clear message naming every missing required column,
# instead of an opaque KeyError deep in the calculation (BACKLOG D4 step 4).
validate_csv_columns(data, "building")
holidays = _german_national_holidays(year)
climate_zone = "9" # Climate zone 9: Germany (VDI 4655)
number_people_household = 2 # Number of people per household (VDI 4655)
# Extract and validate CSV data
try:
YEU_total_heat_kWh = data["Wärmebedarf"].values.astype(float)
data["Gebäudetyp"].values.astype(str)
data["Subtyp"].values.astype(str)
data["WW_Anteil"].values.astype(float)
data["Normaußentemperatur"].values.astype(float)
except KeyError as e:
raise KeyError(f"Missing column in CSV: {e}. Please check CSV file completeness.") from e
except ValueError as e:
raise ValueError(f"Invalid data types in CSV: {e}. Please ensure data is correctly formatted.") from e
# Initialize result containers
total_heat_W = []
heating_heat_W = []
warmwater_heat_W = []
max_heat_requirement_W = []
yearly_time_steps = None
# Mapping of building types to calculation methods
building_type_to_method = {
"EFH": "VDI4655", # Single family house
"MFH": "VDI4655", # Multi-family house
"HEF": "BDEW", # Commercial single family
"HMF": "BDEW", # Commercial multi-family
"GKO": "BDEW", # Office building
"GHA": "BDEW", # Retail building
"GMK": "BDEW", # School building
"GBD": "BDEW", # Hotel building
"GBH": "BDEW", # Restaurant building
"GWA": "BDEW", # Hospital building
"GGA": "BDEW", # Sports facility
"GBA": "BDEW", # Cultural building
"GGB": "BDEW", # Public building
"GPD": "BDEW", # Production building
"GMF": "BDEW", # Mixed-use building
"GHD": "BDEW", # Service building
}
# Process each building in the dataset
for idx, YEU in enumerate(YEU_total_heat_kWh):
current_building_type = str(data.at[idx, "Gebäudetyp"])
current_subtype = str(data.at[idx, "Subtyp"])
current_ww_demand = float(data.at[idx, "WW_Anteil"])
# A building with no (or invalid/NaN) annual demand contributes an all-zero profile.
# pyslpheat rejects annual_heat_kWh <= 0, so we still call it with a 1 kWh placeholder to
# obtain the correctly-shaped time steps + temperatures, then zero the demand arrays below.
zero_demand = not (float(YEU) > 0)
# Determine calculation method
if calc_method == "Datensatz":
try:
current_calc_method = building_type_to_method.get(current_building_type, "VDI4655")
except KeyError:
print(f"Building type '{current_building_type}' not found in mapping, using VDI4655")
current_calc_method = "VDI4655"
else:
current_calc_method = calc_method
# Execute appropriate calculation method
if current_calc_method == "VDI4655":
# Split total demand into heating and hot water components.
if zero_demand:
heating, hot_water = 1.0, 1.0 # placeholder shape; zeroed after the call
else:
heating = YEU * (1 - current_ww_demand)
hot_water = YEU * current_ww_demand
# Calculate VDI 4655 profiles via pyslpheat
df_vdi = vdi4655_calculate(
annual_heating_kWh=heating,
annual_dhw_kWh=hot_water,
annual_electricity_kWh=1, # placeholder; electricity not used downstream
building_type=current_building_type,
number_people_household=number_people_household,
year=year,
climate_zone=climate_zone,
TRY=TRY,
holidays=holidays,
)
yearly_time_steps = df_vdi.index.values
# kWh per 15 min → kW (×4)
hourly_heat_demand_total_kW = df_vdi["Q_total_kWh"].values * 4
hourly_heat_demand_heating_kW = df_vdi["Q_heat_kWh"].values * 4
hourly_heat_demand_warmwater_kW = df_vdi["Q_dhw_kWh"].values * 4
# temperature is hourly in TRY; VDI DataFrame repeats each value 4 times
hourly_air_temperatures = df_vdi["temperature_C"].values[::4]
elif current_calc_method == "BDEW":
# Read optional per-building BDEW parameters (None → pyslpheat uses its defaults)
heating_limit_temp = (
float(data.at[idx, "Heizgrenztemperatur"])
if "Heizgrenztemperatur" in data.columns and pd.notna(data.at[idx, "Heizgrenztemperatur"])
else None
)
heating_exponent = (
float(data.at[idx, "Heizexponent"])
if "Heizexponent" in data.columns and pd.notna(data.at[idx, "Heizexponent"])
else 1.0
)
peak_design_kw = (
float(data.at[idx, "P_max"])
if "P_max" in data.columns
and pd.notna(data.at[idx, "P_max"])
and str(data.at[idx, "P_max"]).strip() not in ("", "None")
else None
)
# Calculate BDEW profiles via pyslpheat (1 kWh placeholder for a zero-demand building).
df_bdew = bdew_calculate(
annual_heat_kWh=(1.0 if zero_demand else YEU),
profile_type=current_building_type,
subtype=current_subtype,
TRY_file_path=TRY,
year=year,
dhw_share=current_ww_demand,
heating_limit_temp=heating_limit_temp,
heating_exponent=heating_exponent,
peak_design_kW=peak_design_kw,
)
yearly_time_steps = df_bdew.index.values
hourly_heat_demand_total_kW = df_bdew["Q_total_kWh"].values
hourly_heat_demand_heating_kW = df_bdew["Q_heat_kWh"].values
hourly_heat_demand_warmwater_kW = df_bdew["Q_dhw_kWh"].values
hourly_air_temperatures = df_bdew["temperature_C"].values
# A zero-demand building keeps the time steps + temperatures but an all-zero profile.
if zero_demand:
hourly_heat_demand_total_kW = np.zeros_like(hourly_heat_demand_total_kW)
hourly_heat_demand_heating_kW = np.zeros_like(hourly_heat_demand_heating_kW)
hourly_heat_demand_warmwater_kW = np.zeros_like(hourly_heat_demand_warmwater_kW)
# Ensure non-negative demand values (clip physical impossible negative values)
hourly_heat_demand_total_kW = np.clip(hourly_heat_demand_total_kW, 0, None)
hourly_heat_demand_heating_kW = np.clip(hourly_heat_demand_heating_kW, 0, None)
hourly_heat_demand_warmwater_kW = np.clip(hourly_heat_demand_warmwater_kW, 0, None)
# Convert to Watts and store results
total_heat_W.append(hourly_heat_demand_total_kW * 1000)
heating_heat_W.append(hourly_heat_demand_heating_kW * 1000)
warmwater_heat_W.append(hourly_heat_demand_warmwater_kW * 1000)
max_heat_requirement_W.append(np.max(hourly_heat_demand_total_kW * 1000))
# Convert lists to numpy arrays for efficient processing
total_heat_W = np.array(total_heat_W)
heating_heat_W = np.array(heating_heat_W)
warmwater_heat_W = np.array(warmwater_heat_W)
max_heat_requirement_W = np.array(max_heat_requirement_W)
# Calculate supply and return temperature curves
supply_temperature_curve, return_temperature_curve = calculate_temperature_curves(data, hourly_air_temperatures)
return (
yearly_time_steps,
total_heat_W,
heating_heat_W,
warmwater_heat_W,
max_heat_requirement_W,
supply_temperature_curve,
return_temperature_curve,
hourly_air_temperatures,
)
[docs]
def calculate_temperature_curves(
data: pd.DataFrame, hourly_air_temperatures: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
"""
Calculate supply and return temperature curves for district heating systems.
:param data: Building data (VLT_max, RLT_max, Steigung_Heizkurve, Normaußentemperatur)
:type data: pd.DataFrame
:param hourly_air_temperatures: Hourly outdoor temperature [°C]
:type hourly_air_temperatures: np.ndarray
:return: Tuple of (supply_temperature_curve, return_temperature_curve)
:rtype: Tuple[np.ndarray, np.ndarray]
.. note::
Weather-compensated curves: T_supply = T_max + slope × (T_outdoor - T_design)
"""
# Extract heating system parameters from building data
supply_temperature_buildings = data["VLT_max"].values.astype(float)
return_temperature_buildings = data["RLT_max"].values.astype(float)
slope = -data["Steigung_Heizkurve"].values.astype(float) # Negative for decreasing curve
min_air_temperatures = data["Normaußentemperatur"].values.astype(float)
# Initialize temperature curve containers
supply_temperature_curve = []
return_temperature_curve = []
# Calculate system temperature difference (constant for each building)
dT = np.expand_dims(supply_temperature_buildings - return_temperature_buildings, axis=1)
# Generate supply temperature curves for each building
for st, s, min_air_temperature in zip(supply_temperature_buildings, slope, min_air_temperatures, strict=False):
# Apply heating curve equation
st_curve = np.where(
hourly_air_temperatures <= min_air_temperature,
st, # Maximum temperature at/below design conditions
st + (s * (hourly_air_temperatures - min_air_temperature)), # Modulated temperature above design
)
supply_temperature_curve.append(st_curve)
# Convert to numpy arrays for efficient operations
supply_temperature_curve = np.array(supply_temperature_curve)
# Calculate return temperature curves (constant spread from supply)
return_temperature_curve = supply_temperature_curve - dT
return supply_temperature_curve, return_temperature_curve