"""
Energy System Module
=====================
Multi-technology energy system modeling with optimization and visualization.
:author: Dipl.-Ing. (FH) Jonas Pfeiffer
"""
import copy
import json
import logging
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib import cm
from scipy.optimize import minimize as scipy_minimize
from districtheatingsim.heat_generators import (
TECH_CLASS_BY_TYPE,
TECH_CLASS_REGISTRY,
ThermalStorageAdapter,
)
from districtheatingsim.heat_generators.json_encoder import CustomJSONEncoder
from districtheatingsim.heat_generators.results import TechnologyResult
from districtheatingsim.utilities.schema import add_meta, check_version
logging.basicConfig(level=logging.INFO)
[docs]
class EnergySystem:
"""
Multi-technology district heating system integration.
:param time_steps: Simulation time steps
:type time_steps: numpy.ndarray
:param load_profile: Hourly thermal load [kW]
:type load_profile: numpy.ndarray
:param VLT_L: Supply temperature profile [°C]
:type VLT_L: numpy.ndarray
:param RLT_L: Return temperature profile [°C]
:type RLT_L: numpy.ndarray
:param TRY_data: Test Reference Year meteorological data
:type TRY_data: object
:param COP_data: Heat pump performance data
:type COP_data: object
:param economic_parameters: Economic parameters dict
:type economic_parameters: dict
.. note::
Supports multi-technology dispatch, storage integration and optimization.
"""
[docs]
def __init__(
self,
time_steps: np.ndarray,
load_profile: np.ndarray,
VLT_L: np.ndarray,
RLT_L: np.ndarray,
TRY_data: object,
COP_data: object,
economic_parameters: dict,
):
"""
Initialize energy system.
:param time_steps: Time steps for simulation
:type time_steps: numpy.ndarray
:param load_profile: Hourly thermal load [kW]
:type load_profile: numpy.ndarray
:param VLT_L: Supply temperature [°C]
:type VLT_L: numpy.ndarray
:param RLT_L: Return temperature [°C]
:type RLT_L: numpy.ndarray
:param TRY_data: Test Reference Year data
:type TRY_data: object
:param COP_data: Heat pump performance data
:type COP_data: object
:param economic_parameters: Economic parameters
:type economic_parameters: dict
"""
self.time_steps = time_steps
self.load_profile = load_profile
self.VLT_L = VLT_L
self.RLT_L = RLT_L
self.TRY_data = TRY_data
self.COP_data = COP_data
self.economic_parameters = economic_parameters
self.technologies = [] # List to store generator objects
self.storage = None
# One TechnologyResult per result row — the single source of truth from
# which the legacy German parallel lists in self.results are projected
# (see _add_tech_result / _project_results). Not serialized.
self.tech_results = []
self.results = {}
# Hours per time step, inferred from the first interval. A single-step
# profile carries no interval to infer from — fall back to 1 h (hourly
# resolution) instead of crashing on the empty np.diff.
if len(self.time_steps) >= 2:
self.duration = (np.diff(self.time_steps[:2]) / np.timedelta64(1, "h"))[0]
else:
self.duration = 1.0
[docs]
def add_technology(self, tech) -> None:
"""
Add a heat generation technology to the energy system.
:param tech: Technology object to add.
:type tech: BaseHeatGenerator
.. note:: Technologies operate based on priority and control strategies.
"""
self.technologies.append(tech)
[docs]
def add_storage(self, storage) -> None:
"""
Add a seasonal thermal energy storage system to the energy system.
:param storage: Seasonal Thermal Energy Storage object.
:type storage: ThermalStorageAdapter
.. note:: Enables temporal decoupling of generation and demand for improved efficiency.
"""
self.storage = storage
[docs]
def initialize_results(self) -> None:
"""
Initialize the results dictionary for energy system calculations.
.. note:: Sets up structure for energy balance, economic, environmental, and performance results.
"""
if not hasattr(self, "results") or not isinstance(self.results, dict):
self.results = {}
self.results.update(
{
"time_steps": self.time_steps,
"Last_L": self.load_profile,
"VLT_L": self.VLT_L,
"RLT_L": self.RLT_L,
"Jahreswärmebedarf": (np.sum(self.load_profile) / 1000) * self.duration,
"Restlast_L": self.load_profile.copy(),
"Restwärmebedarf": (np.sum(self.load_profile) / 1000) * self.duration,
"WGK_Gesamt": 0,
"Strombedarf": 0,
"Strommenge": 0,
"el_Leistungsbedarf_L": np.zeros_like(self.load_profile),
"el_Leistung_L": np.zeros_like(self.load_profile),
"el_Leistung_ges_L": np.zeros_like(self.load_profile),
"specific_emissions_Gesamt": 0,
"primärenergiefaktor_Gesamt": 0,
}
)
# Records are the source of truth; the German lists below are projected
# from them by _project_results() at the end of calculate_mix().
self.tech_results = []
# Ensure lists are initialized or cleared
for key in [
"Wärmeleistung_L",
"colors",
"Wärmemengen",
"Anteile",
"WGK",
"specific_emissions_L",
"primärenergie_L",
"techs",
]:
if key not in self.results:
self.results[key] = []
else:
self.results[key].clear()
[docs]
def set_optimization_variables(self, variables: list, variables_order: list) -> None:
"""
Set optimization variables for technologies.
:param variables: Optimization variable values
:type variables: list
:param variables_order: Variable names
:type variables_order: list
"""
for tech in self.technologies:
if len(variables) > 0:
idx = tech.name.split("_")[-1]
tech.set_parameters(variables, variables_order, idx)
def _add_tech_result(self, record: TechnologyResult) -> None:
"""Append one result row. The single append point for result rows, so the
projected German lists (see _project_results) can never diverge."""
self.tech_results.append(record)
def _project_results(self) -> None:
"""Rebuild the legacy German parallel lists from self.tech_results.
The GUI (results table, pie chart, Sankey) and serialization still read
these keys; they are a pure projection of the records, filled in one pass.
"""
rows = self.tech_results
self.results["techs"] = [r.name for r in rows]
self.results["Wärmeleistung_L"] = [r.heat_output_kW for r in rows]
self.results["Wärmemengen"] = [r.heat_amount_MWh for r in rows]
self.results["Anteile"] = [r.share for r in rows]
self.results["WGK"] = [r.heat_generation_cost for r in rows]
self.results["specific_emissions_L"] = [r.specific_co2 for r in rows]
self.results["primärenergie_L"] = [r.primary_energy for r in rows]
self.results["colors"] = [r.color for r in rows]
[docs]
def aggregate_results(self, tech_results: dict) -> None:
"""
Aggregate technology results into system-level metrics.
:param tech_results: Technology results dictionary
:type tech_results: dict
"""
self._add_tech_result(
TechnologyResult(
name=tech_results.get("tech_name", "unknown"),
heat_output_kW=tech_results.get("Wärmeleistung_L", np.zeros_like(self.load_profile)),
heat_amount_MWh=tech_results.get("Wärmemenge", 0),
share=tech_results.get("Wärmemenge", 0) / self.results["Jahreswärmebedarf"],
heat_generation_cost=tech_results.get("WGK", 0),
specific_co2=tech_results.get("spec_co2_total", 0),
primary_energy=tech_results.get("primärenergie", 0),
color=tech_results.get("color", "gray"),
)
)
if tech_results.get("Wärmemenge", 0) > 1e-6:
self.results["Restlast_L"] -= tech_results.get("Wärmeleistung_L", np.zeros_like(self.load_profile))
self.results["Restwärmebedarf"] -= tech_results.get("Wärmemenge", 0)
self.results["WGK_Gesamt"] += (tech_results["Wärmemenge"] * tech_results["WGK"]) / self.results[
"Jahreswärmebedarf"
]
self.results["specific_emissions_Gesamt"] += (
tech_results["Wärmemenge"] * tech_results["spec_co2_total"]
) / self.results["Jahreswärmebedarf"]
self.results["primärenergiefaktor_Gesamt"] += (
tech_results["primärenergie"] / self.results["Jahreswärmebedarf"]
)
if tech_results.get("Strommenge"):
self.results["Strommenge"] += tech_results["Strommenge"]
self.results["el_Leistung_L"] += tech_results["el_Leistung_L"]
self.results["el_Leistung_ges_L"] += tech_results["el_Leistung_L"]
if tech_results.get("Strombedarf"):
self.results["Strombedarf"] += tech_results["Strombedarf"]
self.results["el_Leistungsbedarf_L"] += tech_results["el_Leistung_L"]
self.results["el_Leistung_ges_L"] -= tech_results["el_Leistung_L"]
if "Wärmeleistung_Speicher_L" in tech_results.keys():
self.results["Restlast_L"] -= tech_results["Wärmeleistung_Speicher_L"]
self._add_tech_result(
TechnologyResult(
name=f"{tech_results['tech_name']}_Speicher",
heat_output_kW=tech_results["Wärmeleistung_Speicher_L"],
heat_amount_MWh=0,
share=0,
heat_generation_cost=0,
specific_co2=0,
primary_energy=0,
color="gray",
)
)
[docs]
def calculate_mix(self, variables: list | None = None, variables_order: list | None = None) -> dict:
"""
Calculate energy generation mix with technology dispatch and storage.
:param variables: Optimization variables, defaults to []
:type variables: list
:param variables_order: Variable order, defaults to []
:type variables_order: list
:return: System results dictionary
:rtype: dict
"""
if variables is None:
variables = []
if variables_order is None:
variables_order = []
self.initialize_results()
# Every share / WGK / emissions term divides by Jahreswärmebedarf; a zero
# (empty or all-zero load profile) would silently turn the whole result set
# to NaN/inf. Fail loud instead.
if self.results["Jahreswärmebedarf"] <= 0:
raise ValueError(
"Jahreswärmebedarf ist 0 (leeres oder Null-Lastprofil) — die "
"Wärmegestehungskosten und Anteile sind nicht definiert."
)
# Initialize optimization variables
self.set_optimization_variables(variables, variables_order)
# Separate STES from generators before iterating (modifying a list during
# iteration skips elements and is undefined behaviour).
stes_list = [t for t in self.technologies if isinstance(t, ThermalStorageAdapter)]
self.technologies = [t for t in self.technologies if not isinstance(t, ThermalStorageAdapter)]
if stes_list:
self.storage = stes_list[-1] # Use last-added storage if multiple were added
for tech in self.technologies:
# Initialize each technology
tech.init_operation(8760)
if self.storage:
self.storage_state = np.zeros(len(self.time_steps))
# Initialize results for each time step
time_steps = len(self.time_steps)
for t in range(time_steps):
Q_in_total = 0 # Total heat input
T_Q_in_flow = self.VLT_L[t] # Supply temperature
T_Q_out_return = self.RLT_L[t] # Return temperature
Q_out_total = self.load_profile[t] # Heat demand
remaining_load = Q_out_total
# Get storage state and temperatures
upper_storage_temperature, lower_storage_temperature = (
self.storage.current_storage_temperatures(t - 1) if t > 0 else (0, 0)
)
# Get storage state and available energy
current_storage_state, available_energy, max_energy = (
self.storage.current_storage_state(t - 1, T_Q_out_return, T_Q_in_flow) if t > 0 else (0, 0, 0)
)
# Calculate storage losses
Q_loss = self.storage.Q_loss[t - 1] if t > 0 else 0
# Control generators based on priority
for _i, tech in enumerate(self.technologies):
tech.active = tech.strategy.decide_operation(
tech.active, upper_storage_temperature, lower_storage_temperature, remaining_load
)
if tech.active:
# Create kwargs dictionary with technology-specific data
kwargs = {
"remaining_load": remaining_load,
"VLT_L": self.VLT_L[t],
"COP_data": self.COP_data,
"time_steps": self.time_steps,
"duration": self.duration,
"TRY_data": self.TRY_data,
"RLT_L": self.RLT_L[t],
"upper_storage_temperature": upper_storage_temperature,
"lower_storage_temperature": lower_storage_temperature,
"current_storage_state": current_storage_state,
"available_energy": available_energy,
"max_energy": max_energy,
"Q_loss": Q_loss,
}
Q_in, _ = tech.generate(t, **kwargs)
remaining_load -= Q_in
Q_in_total += Q_in
tech.calculated = True # Mark technology as calculated
# Update storage
self.storage.simulate_stratified_temperature_mass_flows(
t, Q_in_total, Q_out_total, T_Q_in_flow, T_Q_out_return
)
# Calculate storage results
self.storage.calculate_efficiency(self.load_profile)
self.results["storage_class"] = self.storage
for tech in self.technologies:
# Perform technology-specific calculation
tech_results = tech.calculate(
economic_parameters=self.economic_parameters,
duration=self.duration,
load_profile=self.results["Restlast_L"],
VLT_L=self.VLT_L,
RLT_L=self.RLT_L,
TRY_data=self.TRY_data,
COP_data=self.COP_data,
time_steps=self.time_steps,
)
if tech_results["Wärmemenge"] > 1e-6:
self.aggregate_results(tech_results)
else:
# Add technology as inactive with zero contribution
self.aggregate_results({"tech_name": tech.name})
# Credit network storage discharge against Restlast_L.
# Positive _Q_net_storage_flow = storage discharging = demand covered by storage.
if self.storage:
storage_discharge = np.maximum(self.storage._Q_net_storage_flow, 0.0)
storage_discharge_mwh = np.sum(storage_discharge) / 1000.0 * self.duration
if storage_discharge_mwh > 1e-3:
# Storage cost: capital annuity + maintenance (VDI 2067), spread over the
# discharged energy. No fuel cost – the loss energy is paid on the generator side.
self.storage.calculate_costs(storage_discharge_mwh, self.economic_parameters)
self.results["Restlast_L"] -= storage_discharge
self._add_tech_result(
TechnologyResult(
name=f"{self.storage.name} (Entladung)",
heat_output_kW=storage_discharge,
heat_amount_MWh=storage_discharge_mwh,
share=storage_discharge_mwh / self.results["Jahreswärmebedarf"],
heat_generation_cost=self.storage.WGK,
specific_co2=0.0,
primary_energy=0.0,
color="steelblue",
)
)
self.results["WGK_Gesamt"] += self.storage.A_N / self.results["Jahreswärmebedarf"]
self.results["Restwärmebedarf"] -= storage_discharge_mwh
# Calculate unmet demand after processing all technologies.
# Use np.maximum(..., 0) so that over-production (negative residual, absorbed by
# seasonal storage) does not produce a negative unmet-demand entry.
if np.any(self.results["Restlast_L"] > 1e-6):
unmet_demand = np.sum(np.maximum(self.results["Restlast_L"], 0)) / 1000 * self.duration
self._add_tech_result(
TechnologyResult(
name="Ungedeckter Bedarf",
heat_output_kW=np.maximum(self.results["Restlast_L"], 0),
heat_amount_MWh=unmet_demand,
share=unmet_demand / self.results["Jahreswärmebedarf"],
heat_generation_cost=0,
specific_co2=0,
primary_energy=0,
color="black",
)
)
# Project the records into the legacy German parallel lists consumed by the
# GUI and serialization. Single pass → all eight lists stay in lockstep.
self._project_results()
self.getInitialPlotData()
return self.results
[docs]
def optimize_mix(self, weights: dict, num_restarts: int = 5, unmet_demand_penalty: float = 1e6, seed=None):
"""
Optimize energy mix for multi-objective performance.
:param weights: Optimization weights (WGK_Gesamt, specific_emissions_Gesamt, primärenergiefaktor_Gesamt)
:type weights: dict
:param num_restarts: Number of random restarts, defaults to 5
:type num_restarts: int
:param unmet_demand_penalty: Penalty weight on the uncovered-demand fraction added to the
objective, defaults to 1e6 (see EnergySystemOptimizer for the rationale)
:type unmet_demand_penalty: float
:param seed: Seed for the random-restart draws; ``None`` (default) is non-deterministic,
an int makes ``optimize_mix`` reproducible.
:type seed: int or None
:return: Optimized energy system
:rtype: EnergySystem
"""
optimizer = EnergySystemOptimizer(self, weights, num_restarts, unmet_demand_penalty, seed)
self.optimized_energy_system = optimizer.optimize()
return self.optimized_energy_system
[docs]
def getInitialPlotData(self) -> tuple:
"""
Extract and prepare data for visualization.
:return: (extracted_data, initial_vars)
:rtype: tuple
"""
# Extract data using the explicit get_plot_data() interface instead of dir()
# to avoid picking up inherited attributes and internal helpers.
self.extracted_data = {}
for tech_class in self.technologies:
for var_name, var_value in tech_class.get_plot_data().items():
if len(var_value) == len(self.time_steps):
unique_var_name = f"{tech_class.name}_{var_name}"
self.extracted_data[unique_var_name] = var_value
# Add storage data
if self.storage:
Q_net_storage_flow = self.storage.Q_net_storage_flow
# Separate storage charging (negative values) and discharging (positive values)
Q_net_positive = np.maximum(Q_net_storage_flow, 0) # Storage discharging
Q_net_negative = np.minimum(Q_net_storage_flow, 0) # Storage charging
# Add storage data to extracted data structure
self.extracted_data["Speicherbeladung_kW"] = Q_net_negative
self.extracted_data["Speicherentladung_kW"] = Q_net_positive
# Additional storage state variables for plotting
self.extracted_data["Speicher_SOC_%"] = self.storage._soc * 100.0
self.extracted_data["Speicher_T_oben_°C"] = self.storage._T_supply
self.extracted_data["Speicher_T_mitte_°C"] = self.storage._T_middle
self.extracted_data["Speicher_T_unten_°C"] = self.storage._T_return
self.extracted_data["Speicher_Verluste_kW"] = self.storage.Q_loss
if "Ungedeckter Bedarf" in self.results["techs"]:
# Find index of "Ungedeckter Bedarf" in technology list
if isinstance(self.results["techs"], list):
unmet_demand_index = self.results["techs"].index("Ungedeckter Bedarf")
elif isinstance(self.results["techs"], np.ndarray):
unmet_demand_index = np.where(self.results["techs"] == "Ungedeckter Bedarf")[0][0]
else:
# Skip unmet demand if data type is unknown
unmet_demand_index = None
# Add unmet demand to extracted data structure if index was found
if unmet_demand_index is not None:
self.extracted_data["Ungedeckter_Bedarf_kW"] = self.results["Wärmeleistung_L"][unmet_demand_index]
# Initial selection
self.initial_vars = [var_name for var_name in self.extracted_data.keys() if "_Wärmeleistung" in var_name]
self.initial_vars.append("Last_L")
if self.storage:
self.initial_vars.append("Speicherbeladung_kW")
self.initial_vars.append("Speicherentladung_kW")
return self.extracted_data, self.initial_vars
[docs]
def plot_stack_plot(self, figure=None, selected_vars=None, second_y_axis=False) -> None:
"""
Create stack plot visualization of energy system operation.
:param figure: Figure object, defaults to None
:type figure: matplotlib.figure.Figure, optional
:param selected_vars: Selected variables, defaults to None
:type selected_vars: list, optional
:param second_y_axis: Use second y-axis, defaults to False
:type second_y_axis: bool, optional
"""
if figure is None:
figure = plt.figure()
if selected_vars is None:
selected_vars = self.initial_vars
# X-Achse: Jahresstunden als int
n_steps = len(self.time_steps)
x = np.arange(n_steps)
import matplotlib.gridspec as gridspec
figure.clear()
# Breitere Legenden-Spalten für lange Namen
gs = gridspec.GridSpec(1, 3, width_ratios=[0.22, 0.56, 0.22], figure=figure)
ax_legend_left = figure.add_subplot(gs[0, 0])
ax_main = figure.add_subplot(gs[0, 1])
ax_legend_right = figure.add_subplot(gs[0, 2])
ax_legend_left.axis("off")
ax_legend_right.axis("off")
ax_main.set_prop_cycle(color=cm.tab10.colors)
# Stackplot- und Linienplot-Logik auf ax_main
stackplot_vars = []
if "Speicherbeladung_kW" in selected_vars:
stackplot_vars.append("Speicherbeladung_kW")
if "Speicherentladung_kW" in selected_vars:
stackplot_vars.append("Speicherentladung_kW")
stackplot_vars += [var for var in selected_vars if var not in stackplot_vars and "_Wärmeleistung" in var]
if "Ungedeckter_Bedarf_kW" in selected_vars:
stackplot_vars.append("Ungedeckter_Bedarf_kW")
if "Speicherentladung_kW" in stackplot_vars:
stackplot_vars.remove("Speicherentladung_kW")
stackplot_vars.append("Speicherentladung_kW")
line_vars = [var for var in selected_vars if var not in stackplot_vars and var != "Last_L"]
def _clean_label(var: str) -> str:
"""Strip internal suffixes for display in legend."""
label = var.replace("_Wärmeleistung", "")
label = label.replace("_kW", " (kW)").replace("_%", " (%)")
label = label.replace("_°C", " (°C)")
return label
stackplot_data = []
stackplot_labels = []
for var in stackplot_vars:
if var == "Speicherbeladung_kW" and var in self.extracted_data:
ax_main.fill_between(
x,
0,
self.extracted_data[var],
label=_clean_label(var),
step="mid",
color="gray",
alpha=1.0,
)
elif var in self.extracted_data:
stackplot_data.append(self.extracted_data[var])
stackplot_labels.append(_clean_label(var))
if stackplot_data:
ax_main.stackplot(x, stackplot_data, labels=stackplot_labels, step="mid", edgecolor="none")
ax2 = ax_main.twinx() if second_y_axis else None
lines_ax1 = []
labels_ax1 = []
lines_ax2 = []
labels_ax2 = []
import itertools
color_cycle = itertools.cycle(cm.Dark2.colors)
for var_name in line_vars:
if var_name in self.extracted_data:
display_label = _clean_label(var_name)
if ax2:
(line,) = ax2.plot(x, self.extracted_data[var_name], label=display_label, color=next(color_cycle))
lines_ax2.append(line)
labels_ax2.append(display_label)
else:
(line,) = ax_main.plot(x, self.extracted_data[var_name], label=display_label)
lines_ax1.append(line)
labels_ax1.append(display_label)
if "Last_L" in selected_vars:
(line,) = ax_main.plot(x, self.results["Last_L"], color="black", label="Wärmebedarf", linewidth=0.25)
lines_ax1.append(line)
labels_ax1.append("Wärmebedarf")
# Achsenbeschriftung und Grid
ax_main.set_title("Jahresganglinie", fontsize=16)
ax_main.set_xlabel("", fontsize=14)
ax_main.set_ylabel("Wärmeleistung [kW]", fontsize=14)
ax_main.grid(True, alpha=0.3)
if ax2:
ax2.set_ylabel("Temperatur (°C)", fontsize=14)
ax2.tick_params(axis="y", labelsize=14)
# Use month labels when simulation covers a full year (≥8000 h)
month_names = ["Jan", "Feb", "Mär", "Apr", "Mai", "Jun", "Jul", "Aug", "Sep", "Okt", "Nov", "Dez"]
month_starts = [0, 744, 1416, 2160, 2880, 3624, 4344, 5088, 5832, 6552, 7296, 8016]
if n_steps >= 8000:
valid = [(s, m) for s, m in zip(month_starts, month_names, strict=False) if s < n_steps]
ax_main.set_xticks([s for s, _ in valid])
ax_main.set_xticklabels([m for _, m in valid], fontsize=12)
else:
step = max(1, n_steps // 10)
ax_main.set_xticks(np.arange(0, n_steps + step, step))
ax_main.set_xticklabels([str(i) for i in np.arange(0, n_steps + step, step)])
# Legenden in eigenen Achsen
def get_ncol(n):
return 1 if n <= 18 else 2
if lines_ax1 or stackplot_labels:
ncol_left = get_ncol(len(lines_ax1) + len(stackplot_labels))
ax_legend_left.legend(
ax_main.get_legend_handles_labels()[0],
ax_main.get_legend_handles_labels()[1],
loc="best",
fontsize=12,
frameon=False,
ncol=ncol_left,
)
if lines_ax2:
ncol_right = get_ncol(len(lines_ax2))
ax_legend_right.legend(lines_ax2, labels_ax2, loc="best", fontsize=12, frameon=False, ncol=ncol_right)
# Weniger Rand, damit die Daten direkt an den Achsen anliegen
figure.subplots_adjust(left=0.08, right=0.92, wspace=0.18)
# X-Achse: min/max exakt an Daten
ax_main.set_xlim(x[0], x[-1])
# Y-Achse: min/max exakt an Daten
y_data_ax1 = []
for arr in stackplot_data:
y_data_ax1.append(np.asarray(arr))
for var_name in line_vars:
if var_name in self.extracted_data:
y_data_ax1.append(np.asarray(self.extracted_data[var_name]))
if "Last_L" in selected_vars:
y_data_ax1.append(np.asarray(self.results["Last_L"]))
if y_data_ax1:
y_min = min(arr.min() for arr in y_data_ax1)
y_max = max(arr.max() for arr in y_data_ax1)
ax_main.set_ylim(y_min, y_max)
if ax2:
y_data_ax2 = []
for var_name in line_vars:
if var_name in self.extracted_data:
y_data_ax2.append(np.asarray(self.extracted_data[var_name]))
if y_data_ax2:
y2_min = min(arr.min() for arr in y_data_ax2)
y2_max = max(arr.max() for arr in y_data_ax2)
ax2.set_ylim(y2_min, y2_max)
[docs]
def plot_pie_chart(self, figure=None) -> None:
"""
Create pie chart visualization of technology contributions.
:param figure: Figure object, defaults to None
:type figure: matplotlib.figure.Figure, optional
"""
if figure is None:
figure = plt.figure()
figure.clear()
ax = figure.add_subplot(111)
labels_all = self.results["techs"]
anteile_all = np.asarray(self.results["Anteile"], dtype=float)
colors_all = self.results["colors"]
waermemengen_all = np.asarray(self.results["Wärmemengen"], dtype=float)
# Filter out zero or negative shares
mask = anteile_all > 0
labels = [label for label, m in zip(labels_all, mask, strict=False) if m]
anteile = anteile_all[mask]
colors = [c for c, m in zip(colors_all, mask, strict=False) if m]
waermemengen = waermemengen_all[mask]
# Donut-style pie chart
wedges, _ = ax.pie(
anteile,
labels=None,
colors=colors,
startangle=90,
wedgeprops=dict(width=0.55, edgecolor="white", linewidth=1.5),
pctdistance=0.75,
)
ax.set_title("Anteile Wärmeerzeugung", fontsize=13, pad=12)
ax.axis("equal")
# Legend: name, MWh, percentage
legend_labels = [
f"{label}: {mwh:.0f} MWh ({100 * share:.1f} %)"
for label, mwh, share in zip(labels, waermemengen, anteile, strict=False)
]
ax.legend(wedges, legend_labels, loc="center left", bbox_to_anchor=(0.85, 0, 0.5, 1), fontsize=9, frameon=False)
[docs]
def copy(self):
"""
Create deep copy of EnergySystem instance.
:return: Deep copy of energy system
:rtype: EnergySystem
"""
# Create a new EnergySystem instance with copied basic attributes
copied_system = EnergySystem(
time_steps=self.time_steps.copy(),
load_profile=self.load_profile.copy(),
VLT_L=self.VLT_L.copy(),
RLT_L=self.RLT_L.copy(),
TRY_data=copy.deepcopy(self.TRY_data),
COP_data=copy.deepcopy(self.COP_data),
economic_parameters=copy.deepcopy(self.economic_parameters),
)
# Deep-copy the technologies
copied_system.technologies = [copy.deepcopy(tech) for tech in self.technologies]
# Deep-copy the storage, if it exists
if self.storage:
copied_system.storage = copy.deepcopy(self.storage)
# Deep-copy the results dictionary
copied_system.results = copy.deepcopy(self.results)
# Copy any additional attributes that may have been added dynamically
for attr_name, attr_value in self.__dict__.items():
if attr_name not in copied_system.__dict__:
copied_system.__dict__[attr_name] = copy.deepcopy(attr_value)
return copied_system
[docs]
def to_dict(self) -> dict:
"""
Convert EnergySystem to dictionary for serialization and storage.
Returns
-------
dict
Dictionary representation of the complete energy system.
"""
return add_meta(
{
"time_steps": self.time_steps.astype(str).tolist(), # Convert datetime64 to string
"load_profile": self.load_profile.tolist(),
"VLT_L": self.VLT_L.tolist(),
"RLT_L": self.RLT_L.tolist(),
"TRY_data": [data.tolist() for data in self.TRY_data],
"COP_data": self.COP_data.tolist(),
"economic_parameters": self.economic_parameters,
"technologies": [tech.to_dict() for tech in self.technologies],
"storage": self.storage.to_dict() if self.storage else None,
"results": {
key: (value.to_dict(orient="split") if isinstance(value, pd.DataFrame) else value)
for key, value in self.results.items()
},
},
"energy_system",
)
[docs]
@classmethod
def from_dict(cls, data: dict):
"""
Recreate EnergySystem instance from dictionary representation.
Parameters
----------
data : dict
Dictionary representation of the EnergySystem.
Returns
-------
EnergySystem
Fully initialized EnergySystem object.
"""
# Schema version: tolerates the _meta block, the legacy top-level "version"
# field, and pre-versioning files (0); warns if newer than this app.
check_version(data, "energy_system")
# Restore basic attributes
time_steps = np.array(data["time_steps"], dtype="datetime64")
load_profile = np.array(data["load_profile"])
VLT_L = np.array(data["VLT_L"])
RLT_L = np.array(data["RLT_L"])
TRY_data = [np.array(item) for item in data["TRY_data"]]
COP_data = np.array(data["COP_data"])
economic_parameters = data["economic_parameters"]
# Create the EnergySystem object
obj = cls(
time_steps=time_steps,
load_profile=load_profile,
VLT_L=VLT_L,
RLT_L=RLT_L,
TRY_data=TRY_data,
COP_data=COP_data,
economic_parameters=economic_parameters,
)
# Restore technologies — prefer tech_type (class name) stored since v1.0.2,
# fall back to prefix-matching for JSON files written by older versions.
obj.technologies = []
for tech_data in data.get("technologies", []):
tech_class = TECH_CLASS_BY_TYPE.get(tech_data.get("tech_type"))
if tech_class is None:
for prefix, cls in TECH_CLASS_REGISTRY.items():
if tech_data["name"].startswith(prefix):
tech_class = cls
break
if tech_class is not None:
obj.technologies.append(tech_class.from_dict(tech_data))
else:
logging.warning(
"Could not restore technology '%s': unknown type '%s'",
tech_data.get("name"),
tech_data.get("tech_type"),
)
# Restore storage
if data.get("storage"):
obj.storage = ThermalStorageAdapter.from_dict(data["storage"])
if obj.storage is None:
logging.warning(
"Thermal storage could not be loaded (outdated format). Please re-configure the storage in the GUI."
)
# Records are not serialized (the German lists in results carry the data);
# the GUI reads the lists after a load, so an empty record list is fine.
obj.tech_results = []
# Restore results (if available)
obj.results = {}
if "results" in data:
for key, value in data["results"].items():
if isinstance(value, dict) and "columns" in value and "data" in value:
obj.results[key] = pd.DataFrame(**value)
elif isinstance(value, list):
if all(isinstance(v, list) for v in value):
obj.results[key] = [np.array(v) for v in value]
else:
obj.results[key] = np.array(value)
else:
obj.results[key] = value
return obj
[docs]
def save_to_csv(self, file_path: str) -> None:
"""
Save energy system results to CSV file.
:param file_path: Path for CSV output
:type file_path: str
"""
if not self.results:
raise ValueError("No results available to save.")
# Initialize the DataFrame with the timestamps
df = pd.DataFrame({"time_steps": self.results["time_steps"]})
# Add the load data
df["Last_L"] = self.results["Last_L"]
# Add the heat generation data for each technology
for tech_results, techs in zip(self.results["Wärmeleistung_L"], self.results["techs"], strict=False):
df[techs] = tech_results
# Add the electrical power data
df["el_Leistungsbedarf_L"] = self.results["el_Leistungsbedarf_L"]
df["el_Leistung_L"] = self.results["el_Leistung_L"]
df["el_Leistung_ges_L"] = self.results["el_Leistung_ges_L"]
# Save the DataFrame as a CSV file
df.to_csv(file_path, index=False, sep=";", encoding="utf-8-sig")
[docs]
def save_to_json(self, file_path: str) -> None:
"""
Save complete EnergySystem object to JSON file for persistence.
Parameters
----------
file_path : str
Path for JSON file output.
"""
with open(file_path, "w") as json_file:
json.dump(self.to_dict(), json_file, indent=4, cls=CustomJSONEncoder)
[docs]
@classmethod
def load_from_json(cls, file_path: str):
"""
Load complete EnergySystem object from JSON file.
Parameters
----------
file_path : str
Path to JSON file for loading.
Returns
-------
EnergySystem
Loaded EnergySystem object with complete configuration.
"""
try:
with open(file_path) as json_file:
data_loaded = json.load(json_file)
return cls.from_dict(data_loaded)
except Exception as e:
raise ValueError(f"Error loading JSON file: {e}") from e
[docs]
class EnergySystemOptimizer:
"""
Multi-objective optimizer for energy system configuration.
:param initial_energy_system: Initial system configuration
:type initial_energy_system: EnergySystem
:param weights: Optimization weights dict with 'WGK_Gesamt', 'specific_emissions_Gesamt', 'primärenergiefaktor_Gesamt'
:type weights: dict
:param num_restarts: Number of random restart runs, defaults to 5
:type num_restarts: int, optional
.. note::
Uses SLSQP with random restarts for multi-objective optimization.
"""
[docs]
def __init__(
self,
initial_energy_system: "EnergySystem",
weights: dict[str, float],
num_restarts: int = 5,
unmet_demand_penalty: float = 1e6,
seed=None,
):
"""
Initialize multi-objective optimizer.
:param initial_energy_system: Initial system configuration
:type initial_energy_system: EnergySystem
:param weights: Optimization weights
:type weights: dict
:param num_restarts: Number of random restarts, defaults to 5
:type num_restarts: int
:param unmet_demand_penalty: Penalty weight applied to the uncovered-demand fraction
(Restwärmebedarf / Jahreswärmebedarf) and added to the objective, defaults to 1e6.
Without it the objective (WGK + emissions + primary-energy, all divided by the full
annual demand) is minimised by *shrinking* generators: less generation → the uncovered
load lands in the cost-free "Ungedeckter Bedarf" row → every term falls toward 0, so the
optimum is an empty/non-covering system (verified: a CHP collapses to 0 kW / 0 % coverage).
A large penalty makes covering demand strictly dominate the cost saving from undersizing.
:type unmet_demand_penalty: float
:raises ValueError: If required weights missing or negative
"""
self.initial_energy_system = initial_energy_system
self.weights = weights
self.num_restarts = num_restarts
self.unmet_demand_penalty = unmet_demand_penalty
# Local RNG so the random restarts are controllable (seed=None stays non-deterministic);
# avoids seeding the global np.random and lets optimize_mix be golden-mastered.
self.rng = np.random.default_rng(seed)
# Validate optimization weights
required_weights = ["WGK_Gesamt", "specific_emissions_Gesamt", "primärenergiefaktor_Gesamt"]
for weight_key in required_weights:
if weight_key not in weights:
raise ValueError(f"Required weight '{weight_key}' missing from weights dictionary")
if weights[weight_key] < 0:
raise ValueError(f"Weight '{weight_key}' must be non-negative")
[docs]
def optimize(self) -> "EnergySystem":
"""
Perform multi-objective optimization with random restarts.
:return: Optimized energy system
:rtype: EnergySystem
:raises ValueError: If no optimization parameters available
:raises RuntimeError: If optimization fails in all restarts
"""
best_solution = None
best_objective_value = float("inf")
# Validate that technologies have optimization parameters
has_optimization_params = False
for tech in self.initial_energy_system.technologies:
idx = tech.name.split("_")[-1] if "_" in tech.name else "0"
tech_values, tech_variables, tech_bounds = tech.add_optimization_parameters(idx)
if tech_values and tech_variables and tech_bounds:
has_optimization_params = True
break
if not has_optimization_params:
raise ValueError(
"No optimization parameters available. Energy system optimization requires "
"technologies with configurable parameters (e.g., capacity, storage volume)."
)
for restart in range(self.num_restarts):
logging.info("Starting optimization run %d/%d", restart + 1, self.num_restarts)
# Create fresh copy for this optimization run
self.energy_system_copy = self.initial_energy_system.copy()
# Extract optimization parameters from all technologies
initial_values = []
bounds = []
variables_mapping = {}
for tech in self.energy_system_copy.technologies:
idx = tech.name.split("_")[-1] if "_" in tech.name else "0"
tech_values, tech_variables, tech_bounds = tech.add_optimization_parameters(idx)
# Skip technologies without optimization parameters
if not tech_values or not tech_variables or not tech_bounds:
continue
initial_values.extend(tech_values)
bounds.extend(tech_bounds)
# Map variables to technology for solution interpretation
for var in tech_variables:
variables_mapping[var] = tech.name
variables_order = list(variables_mapping.keys())
if not initial_values:
logging.warning("No optimization parameters found. Skipping optimization.")
return self.initial_energy_system
# Generate random initial values within parameter bounds
random_initial_values = [
self.rng.uniform(low=bound[0], high=bound[1]) if bound[1] > bound[0] else bound[0] for bound in bounds
]
logging.debug("Initial values for restart %d: %s", restart + 1, random_initial_values)
def objective_function(variables, variables_order=variables_order):
"""
Multi-objective function for energy system optimization.
Parameters
----------
variables : array_like
Technology parameter values for evaluation.
Returns
-------
float
Weighted sum of optimization criteria.
"""
try:
# Create fresh copy for objective evaluation
fresh_energy_system = self.energy_system_copy.copy()
# Calculate energy system performance with given parameters
results = fresh_energy_system.calculate_mix(variables, variables_order)
# Penalise uncovered demand so the optimiser cannot lower the objective by
# undersizing generators (the uncovered load would otherwise land in the
# cost-free "Ungedeckter Bedarf" row and drag every term toward 0). The penalty
# is proportional to the uncovered *fraction*, so it stays well-defined even when
# full coverage is physically impossible (it then minimises the gap, then cost).
jahresbedarf = results["Jahreswärmebedarf"]
unmet_fraction = max(results["Restwärmebedarf"], 0.0) / jahresbedarf if jahresbedarf > 0 else 0.0
# Calculate weighted multi-objective value
weighted_sum = (
self.weights["WGK_Gesamt"] * results["WGK_Gesamt"]
+ self.weights["specific_emissions_Gesamt"] * results["specific_emissions_Gesamt"]
+ self.weights["primärenergiefaktor_Gesamt"] * results["primärenergiefaktor_Gesamt"]
+ self.unmet_demand_penalty * unmet_fraction
)
return weighted_sum
except Exception as e:
logging.debug("Error in objective function evaluation: %s", e)
return float("inf") # Return large value for infeasible solutions
# Perform optimization with SLSQP algorithm
try:
result = scipy_minimize(
objective_function,
random_initial_values,
method="SLSQP",
bounds=bounds,
options={"maxiter": 1000, "ftol": 1e-6},
)
# Check if current solution is better than previous best
if result.success and result.fun < best_objective_value:
best_objective_value = result.fun
best_solution = result
logging.info("New best solution found in restart %d: %.4f", restart + 1, result.fun)
except Exception as e:
logging.warning("Optimization failed in restart %d: %s", restart + 1, e)
continue
# Apply best solution if found
if best_solution is not None:
logging.info("Optimization completed. Best objective value: %.4f", best_objective_value)
# Apply optimal parameters to energy system
for tech in self.energy_system_copy.technologies:
idx = tech.name.split("_")[-1] if "_" in tech.name else "0"
tech.set_parameters(best_solution.x, variables_order, idx)
# Store optimization results
self.best_solution = best_solution
self.best_objective_value = best_objective_value
return self.energy_system_copy
else:
raise RuntimeError(
"Optimization failed to find valid solution in all restart attempts. "
"Consider adjusting parameter bounds, weights, or increasing restart attempts."
)
[docs]
def get_optimization_summary(self) -> dict[str, float | int | bool]:
"""
Generate optimization summary report.
:return: Summary dict with success, best_objective_value, num_restarts, etc.
:rtype: dict
"""
if hasattr(self, "best_solution") and self.best_solution is not None:
return {
"success": True,
"best_objective_value": self.best_objective_value,
"num_restarts": self.num_restarts,
"optimization_message": f"Optimization successful with {self.num_restarts} restarts",
"solution_variables": self.best_solution.x.tolist(),
"function_evaluations": getattr(self.best_solution, "nfev", 0),
"iterations": getattr(self.best_solution, "nit", 0),
}
else:
return {
"success": False,
"best_objective_value": float("inf"),
"num_restarts": self.num_restarts,
"optimization_message": "Optimization failed to find valid solution",
"solution_variables": [],
"function_evaluations": 0,
"iterations": 0,
}