From 03372d4865e7921f70ed0d759196967b7ded042d Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Mon, 5 Jan 2026 08:26:50 +0100 Subject: [PATCH 01/19] first work on ensemble optimization --- src/mesido/esdl/esdl_heat_model.py | 1 + src/mesido/esdl/esdl_mixin.py | 84 ++- src/mesido/esdl/profile_parser.py | 5 +- src/mesido/heat_physics_mixin.py | 13 + .../milp/heat/heat_source.py | 1 + .../case_2a_ensemble/input/ensemble.csv | 3 + .../input/forecast1/timeseries_import.csv | 4 + .../input/forecast2/timeseries_import.csv | 4 + .../input/timeseries_import.csv | 4 + .../input/timeseries_import.xml | 52 ++ .../unit_cases/case_2a_ensemble/model/2a.esdl | 521 ++++++++++++++++++ .../case_2a_ensemble/output/.gitignore | 2 + .../unit_cases/case_2a_ensemble/src/run_2a.py | 207 +++++++ tests/test_warmingup_unit_cases.py | 33 ++ 14 files changed, 932 insertions(+), 2 deletions(-) create mode 100644 tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv create mode 100644 tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv create mode 100644 tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv create mode 100644 tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv create mode 100644 tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml create mode 100644 tests/models/unit_cases/case_2a_ensemble/model/2a.esdl create mode 100644 tests/models/unit_cases/case_2a_ensemble/output/.gitignore create mode 100644 tests/models/unit_cases/case_2a_ensemble/src/run_2a.py diff --git a/src/mesido/esdl/esdl_heat_model.py b/src/mesido/esdl/esdl_heat_model.py index 88c36417c..89f3cedae 100644 --- a/src/mesido/esdl/esdl_heat_model.py +++ b/src/mesido/esdl/esdl_heat_model.py @@ -1280,6 +1280,7 @@ def convert_heat_source(self, asset: Asset) -> Tuple[Type[HeatSource], MODIFIERS modifiers = dict( Heat_source=dict(min=0.0, max=max_supply, nominal=max_supply / 2.0), + Max_heat=max_supply, **self._generic_modifiers(asset), **self._generic_heat_modifiers(0.0, max_supply, q_nominal), **self._supply_return_temperature_modifiers(asset), diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index 208a9f378..f21a35efa 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -2,12 +2,14 @@ import copy import dataclasses import logging +import os import sys from datetime import timedelta from pathlib import Path from typing import Any, Dict, List, Optional, Tuple import esdl.esdl_handler +from rtctools.data import csv from mesido.component_type_mixin import ( ModelicaComponentTypeMixin, @@ -691,6 +693,75 @@ def read(self) -> None: None """ super().read() + # TODO: update for ensembles, eg ensemble_size is now still 1, if io gets updated, + # then the size also updates automatically + if self.csv_ensemble_mode: + ensemble_size=2 + self.__ensemble = np.genfromtxt( + os.path.join(self._input_folder, "ensemble.csv"), + delimiter=",", + deletechars="", + dtype=None, + names=True, + encoding=None, + ) + # + # for ensemble_member_index, ensemble_member_name in enumerate(self.__ensemble["name"]): + # _timeseries = csv.load( + # os.path.join( + # self._input_folder, + # ensemble_member_name, + # "timeseries_import.csv", + # ), + # delimiter=",", + # with_time=True, + # ) + # self.__timeseries_times = _timeseries[_timeseries.dtype.names[0]] + # + # self.io.reference_datetime = self.__timeseries_times[0] + # + # for key in _timeseries.dtype.names[1:]: + # self.io.set_timeseries( + # key, + # self.__timeseries_times, + # np.asarray(_timeseries[key], dtype=np.float64), + # ensemble_member_index, + # ) + # for ensemble_member_index, ensemble_member_name in enumerate(self.__ensemble["name"]): + # try: + # _parameters = csv.load( + # os.path.join( + # self._input_folder, + # ensemble_member_name, + # self.csv_parameters_basename + ".csv", + # ), + # delimiter=self.csv_delimiter, + # ) + # for key in _parameters.dtype.names: + # self.io.set_parameter(key, float(_parameters[key]), ensemble_member_index) + # except IOError: + # pass + # logger.debug("CSVMixin: Read parameters.") + # + # for ensemble_member_name in self.__ensemble["name"]: + # try: + # _initial_state = csv.load( + # os.path.join( + # self._input_folder, + # ensemble_member_name, + # self.csv_initial_state_basename + ".csv", + # ), + # delimiter=self.csv_delimiter, + # ) + # check_initial_state_array(_initial_state) + # _initial_state = { + # key: float(_initial_state[key]) for key in _initial_state.dtype.names + # } + # except IOError: + # _initial_state = {} + # self.__initial_state.append(AliasDict(self.alias_relation, _initial_state)) + # logger.debug("CSVMixin: Read initial state.") + energy_system_components = self.energy_system_components esdl_carriers = self.esdl_carriers self.hot_cold_pipe_relations() @@ -701,8 +772,13 @@ def read(self) -> None: esdl_asset_id_to_name_map=self.esdl_asset_id_to_name_map, esdl_assets=self.esdl_assets, carrier_properties=esdl_carriers, - ensemble_size=self.ensemble_size, + ensemble_size=ensemble_size, ) + for ensemble_member_index in range(ensemble_size): + self.io.set_parameter("GeothermalSource_fafd.Max_heat", 10e5, ensemble_member_index) + if ensemble_member_index == 1: + self.io.set_parameter('GeothermalSource_fafd.Max_heat', 2e5, + ensemble_member_index) def write(self) -> None: """ @@ -826,3 +902,9 @@ def filter_asset_measures( filtered_assets[asset_id] = asset_type return filtered_assets + + def ensemble_member_probability(self, ensemble_member): + if self.csv_ensemble_mode: + return self.__ensemble["probability"][ensemble_member] + else: + return 1.0 diff --git a/src/mesido/esdl/profile_parser.py b/src/mesido/esdl/profile_parser.py index 74a156d56..5fbbe9d2e 100644 --- a/src/mesido/esdl/profile_parser.py +++ b/src/mesido/esdl/profile_parser.py @@ -112,7 +112,8 @@ def read_profiles( "power" ] if profile is not None: - values = profile + correction = 1-ensemble_member * 0.1 + values = profile *correction else: if "heat_demand" not in component_type: # We don't set a default profile for source targets @@ -639,6 +640,8 @@ def _load_csv( carrier_properties: Dict[str, Dict], ensemble_size: int, ) -> None: + #TODO: this data should be loaded for the different ensemble members if they each have + # their own profile. data = pd.read_csv(self._file_path) if len(data.filter(like="Unnamed").columns) > 0: diff --git a/src/mesido/heat_physics_mixin.py b/src/mesido/heat_physics_mixin.py index 43441cfd8..a12e2f655 100644 --- a/src/mesido/heat_physics_mixin.py +++ b/src/mesido/heat_physics_mixin.py @@ -1506,6 +1506,12 @@ def __source_heat_to_discharge_path_constraints(self, ensemble_member): discharge = self.state(f"{s}.Q") heat_out = self.state(f"{s}.HeatOut.Heat") + #somehow the ensemblemember dependent path constraint is only allowed in the bounds of + # the problem, e.g. the upper or lower bound. For non path constraints, it is allowed + # in the equation. + heat_flow = self.state(f"{s}.Heat_flow") + max_heat = parameters[f"{s}.Max_heat"] + constraints.append((heat_flow/heat_nominal, -np.inf, max_heat/heat_nominal)) constraint_nominal = (heat_nominal * cp * rho * dt * q_nominal) ** 0.5 @@ -3755,6 +3761,13 @@ def constraints(self, ensemble_member): """ constraints = super().constraints(ensemble_member) + # parameters = self.parameters(ensemble_member) + # for s in self.energy_system_components.get("heat_source", []): + # heat_nominal = parameters[f"{s}.Heat_nominal"] + # heat_flow = self.__state_vector_scaled(f"{s}.Heat_flow", ensemble_member) + # max_heat = parameters[f"{s}.Max_heat"] + # constraints.append(((heat_flow-max_heat)/heat_nominal, -np.inf, 0.0)) + if self.heat_network_settings["head_loss_option"] != HeadLossOption.NO_HEADLOSS: # constraints.extend(self._hn_pipe_head_loss_constraints(ensemble_member)) constraints.extend( diff --git a/src/mesido/pycml/component_library/milp/heat/heat_source.py b/src/mesido/pycml/component_library/milp/heat/heat_source.py index 634c35f5f..c1ecdb8b3 100644 --- a/src/mesido/pycml/component_library/milp/heat/heat_source.py +++ b/src/mesido/pycml/component_library/milp/heat/heat_source.py @@ -44,6 +44,7 @@ def __init__(self, name, **modifiers): self.price = nan # TODO: delete not needed anymore self.co2_coeff = 1.0 self.pump_efficiency = 0.5 + self.Max_heat = nan # Assumption: heat in/out and added is nonnegative # Heat in the return (i.e. cold) line is zero diff --git a/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv b/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv new file mode 100644 index 000000000..11521b7ac --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv @@ -0,0 +1,3 @@ +Unnamed: 0,name,probability +0,forecast1,0.7 +1,forecast2,0.3 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv new file mode 100644 index 000000000..1e9ee86c7 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv @@ -0,0 +1,4 @@ +DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 +2013-05-19 22:00:00,350000,350000,350000 +2013-05-19 23:00:00,350000,350000,350000 +2013-05-20 00:00:00,350000,350000,350000 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv new file mode 100644 index 000000000..1e9ee86c7 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv @@ -0,0 +1,4 @@ +DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 +2013-05-19 22:00:00,350000,350000,350000 +2013-05-19 23:00:00,350000,350000,350000 +2013-05-20 00:00:00,350000,350000,350000 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv new file mode 100644 index 000000000..1e9ee86c7 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv @@ -0,0 +1,4 @@ +DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 +2013-05-19 22:00:00,350000,350000,350000 +2013-05-19 23:00:00,350000,350000,350000 +2013-05-20 00:00:00,350000,350000,350000 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml new file mode 100644 index 000000000..b1defa019 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml @@ -0,0 +1,52 @@ + + + 0.0 + +
+ instantaneous + + 74846440-61d9-4f57-9bce-77e3664ba832 + b145eba8-612d-42b7-bea7-7e22a2440509 + + + + -999.0 + W +
+ + + + + +
+ instantaneous + + c6c8d8f7-b6ae-4e8b-be9a-72282a5c3f1e + 874dc021-09f4-46ec-bb5a-c3276916fbfd + + + + -999.0 + W +
+ + + + + +
+ instantaneous + + 6f99424a-004e-4fb4-8edb-17c374d5be6b + 30713528-878b-46fd-82e2-9da8eb95075e + + + + -999.0 + W +
+ + + + + diff --git a/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl b/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl new file mode 100644 index 000000000..ac27eb46a --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl @@ -0,0 +1,521 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/models/unit_cases/case_2a_ensemble/output/.gitignore b/tests/models/unit_cases/case_2a_ensemble/output/.gitignore new file mode 100644 index 000000000..c4697e99f --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/output/.gitignore @@ -0,0 +1,2 @@ +diag.xml +timeseries_export.xml diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py new file mode 100644 index 000000000..13e709a59 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -0,0 +1,207 @@ +from rtctools.optimization.control_tree_mixin import ControlTreeMixin +from rtctools.optimization.csv_mixin import CSVMixin + +from mesido.asset_sizing_mixin import AssetSizingMixin +from mesido.esdl.esdl_mixin import ESDLMixin +from mesido.esdl.esdl_parser import ESDLFileParser +from mesido.esdl.profile_parser import ProfileReaderFromFile +from mesido.physics_mixin import PhysicsMixin +from mesido.qth_not_maintained.qth_mixin import QTHMixin + +import numpy as np + +from rtctools.optimization.collocated_integrated_optimization_problem import ( + CollocatedIntegratedOptimizationProblem, +) +from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin +from rtctools.optimization.homotopy_mixin import HomotopyMixin +from rtctools.optimization.linearized_order_goal_programming_mixin import ( + LinearizedOrderGoalProgrammingMixin, +) +from rtctools.optimization.single_pass_goal_programming_mixin import SinglePassGoalProgrammingMixin +from rtctools.util import run_optimization_problem + + +class TargetDemandGoal(Goal): + priority = 1 + + order = 2 + + def __init__(self, state, target): + self.state = state + + self.targets = target + self.target_min = target[0] + self.target_max = target[0] + self.function_range = (0.0, 2.0 * max(target[0].values)) + self.function_nominal = np.median(target[0].values) + + def function(self, optimization_problem, ensemble_member): + """Note that for path goals, the ensemble member index is not passed to the call + to :func:`OptimizationProblem.state`. This call returns a time-independent symbol + that is also independent of the active ensemble member. Path goals are + applied to all times and all ensemble members simultaneously.""" + self.target_min = self.targets[ensemble_member] + self.target_max = self.targets[ensemble_member] + return optimization_problem.state(self.state) + + +class MinimizeSourcesHeatGoal(Goal): + priority = 2 + + order = 2 + + def __init__(self, source): + self.target_max = 0.0 + self.function_range = (0.0, 2e6) + self.source = source + self.function_nominal = 1e5 + + def function(self, optimization_problem, ensemble_member): + return optimization_problem.state(f"{self.source}.Heat_source") + +class MinimizeSourcesSizeGoal(Goal): + priority = 2 + + order = 2 + + def __init__(self, source): + self.target_max = 0.0 + self.function_range = (0.0, 2e6) + self.source = source + self.function_nominal = 1e5 + + def function(self, optimization_problem, ensemble_member): + return optimization_problem.extra_variable(f"{self.source}__max_size", ensemble_member) + + +class MinimizeSourcesQTHGoal(Goal): + priority = 2 + + order = 2 + + def __init__(self, source): + self.source = source + self.function_nominal = 1e5 + + def function(self, optimization_problem, ensemble_member): + return optimization_problem.state(f"{self.source}.Heat_source") + + +class _GoalsAndOptions: + def path_goals(self): + goals = super().path_goals().copy() + + for demand in self.energy_system_components["heat_demand"]: + target = [self.get_timeseries(f"{demand}.target_heat_demand", + ensemble_member=ensemble_member) for ensemble_member in + range(self.ensemble_size)] + state = f"{demand}.Heat_demand" + + goals.append(TargetDemandGoal(state, target)) + + return goals + + +class HeatProblem( + _GoalsAndOptions, + PhysicsMixin, + LinearizedOrderGoalProgrammingMixin, + GoalProgrammingMixin, + ESDLMixin, + CollocatedIntegratedOptimizationProblem, +): + # def path_goals(self): + # goals = super().path_goals().copy() + # + # for s in self.energy_system_components["heat_source"]: + # goals.append(MinimizeSourcesHeatGoal(s)) + # + # return goals + + def energy_system_options(self): + options = super().energy_system_options() + options["heat_loss_disconnected_pipe"] = True + self.heat_network_settings["minimum_velocity"] = 0.0001 + + return options + + def solver_options(self): + options = super().solver_options() + options["solver"] = "highs" + return options + + +class HeatProblemEnsemble(AssetSizingMixin, + HeatProblem, + # CSVMixin, + ControlTreeMixin): + csv_ensemble_mode = True + + def pre(self): + super().pre() + + # Empty dict for intermediate ensemble results + self.intermediate_results = {ensemble_member: [] for ensemble_member in range( + self.ensemble_size)} + + def goals(self): + goals = super().goals().copy() + for s in self.energy_system_components["heat_source"]: + goals.append(MinimizeSourcesSizeGoal(s)) + return goals + + def bounds(self): + bounds = super().bounds() + print('a') + return bounds + + def parameters(self, ensemble_member): + parameters = super().parameters(ensemble_member) + if ensemble_member == 1: + parameters['GeothermalSource_fafd.Max_heat']=1e5 + return parameters + + def __fixed_max_size(self): + constraints=[] + for heat_source in self.energy_system_components["heat_source"]: + max_size_prev = self.extra_variable(f"{heat_source}__max_size", ensemble_member=0) + for e_m in range(self.ensemble_size-1): + max_size = self.extra_variable(f"{heat_source}__max_size", ensemble_member=e_m+1) + constraints.append((max_size - max_size_prev, 0.0, 0.0)) + max_size_prev = max_size + return constraints + + def constraints(self, ensemble_member): + constraints = super().constraints(ensemble_member).copy() + constraints.extend(self.__fixed_max_size()) + return constraints + + + +class QTHProblem( + _GoalsAndOptions, + QTHMixin, + HomotopyMixin, + SinglePassGoalProgrammingMixin, + ESDLMixin, + CollocatedIntegratedOptimizationProblem, +): + def path_goals(self): + goals = super().path_goals().copy() + + for s in self.energy_system_components["heat_source"]: + goals.append(MinimizeSourcesQTHGoal(s)) + + return goals + + +if __name__ == "__main__": + run_optimization_problem( + HeatProblem, + esdl_file_name="2a.esdl", + esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="timeseries_import.xml", + ) + # run_heat_network_optimization(HeatProblem, QTHProblem) diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 3a53f1665..0de02841c 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -193,6 +193,39 @@ def test_3a(self): atol=1.0e-6, ) + def test_2a_ensemble(self): + """ + This is the most basic check where we have a simple network and check for the basic physics. + This simple network includes two source, pipes, nodes, and 3 demands. + + Checks; + - Demand matching + - Energy conservation + - Heat to discharge + + """ + import models.unit_cases.case_2a_ensemble.src.run_2a as run_2a + from models.unit_cases.case_2a_ensemble.src.run_2a import HeatProblemEnsemble + + base_folder = Path(run_2a.__file__).resolve().parent.parent + + # Just a "problem is not infeasible" + heat_problem = run_esdl_mesido_optimization( + HeatProblemEnsemble, + base_folder=base_folder, + esdl_file_name="2a.esdl", + esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="timeseries_import.csv", + ) + + results = {} + for ensemble_member in range(heat_problem.ensemble_size): + results[ensemble_member] = heat_problem.extract_results(ensemble_member=ensemble_member) + + demand_matching_test(heat_problem, heat_problem.extract_results()) + energy_conservation_test(heat_problem, heat_problem.extract_results()) + heat_to_discharge_test(heat_problem, heat_problem.extract_results()) if __name__ == "__main__": From e726c2569449132e4153b168f0b8ec4f4a9a1b53 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 30 Jan 2026 20:01:36 +0100 Subject: [PATCH 02/19] clean up parsing profiles ensembles --- src/mesido/esdl/esdl_mixin.py | 61 ++------------------------ src/mesido/esdl/profile_parser.py | 72 ++++++++++++++++++------------- 2 files changed, 46 insertions(+), 87 deletions(-) diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index f21a35efa..0d66faf1a 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -695,8 +695,9 @@ def read(self) -> None: super().read() # TODO: update for ensembles, eg ensemble_size is now still 1, if io gets updated, # then the size also updates automatically + ensemble_size=1 + self.__ensemble=None if self.csv_ensemble_mode: - ensemble_size=2 self.__ensemble = np.genfromtxt( os.path.join(self._input_folder, "ensemble.csv"), delimiter=",", @@ -705,62 +706,7 @@ def read(self) -> None: names=True, encoding=None, ) - # - # for ensemble_member_index, ensemble_member_name in enumerate(self.__ensemble["name"]): - # _timeseries = csv.load( - # os.path.join( - # self._input_folder, - # ensemble_member_name, - # "timeseries_import.csv", - # ), - # delimiter=",", - # with_time=True, - # ) - # self.__timeseries_times = _timeseries[_timeseries.dtype.names[0]] - # - # self.io.reference_datetime = self.__timeseries_times[0] - # - # for key in _timeseries.dtype.names[1:]: - # self.io.set_timeseries( - # key, - # self.__timeseries_times, - # np.asarray(_timeseries[key], dtype=np.float64), - # ensemble_member_index, - # ) - # for ensemble_member_index, ensemble_member_name in enumerate(self.__ensemble["name"]): - # try: - # _parameters = csv.load( - # os.path.join( - # self._input_folder, - # ensemble_member_name, - # self.csv_parameters_basename + ".csv", - # ), - # delimiter=self.csv_delimiter, - # ) - # for key in _parameters.dtype.names: - # self.io.set_parameter(key, float(_parameters[key]), ensemble_member_index) - # except IOError: - # pass - # logger.debug("CSVMixin: Read parameters.") - # - # for ensemble_member_name in self.__ensemble["name"]: - # try: - # _initial_state = csv.load( - # os.path.join( - # self._input_folder, - # ensemble_member_name, - # self.csv_initial_state_basename + ".csv", - # ), - # delimiter=self.csv_delimiter, - # ) - # check_initial_state_array(_initial_state) - # _initial_state = { - # key: float(_initial_state[key]) for key in _initial_state.dtype.names - # } - # except IOError: - # _initial_state = {} - # self.__initial_state.append(AliasDict(self.alias_relation, _initial_state)) - # logger.debug("CSVMixin: Read initial state.") + ensemble_size = len(self.__ensemble) energy_system_components = self.energy_system_components esdl_carriers = self.esdl_carriers @@ -773,6 +719,7 @@ def read(self) -> None: esdl_assets=self.esdl_assets, carrier_properties=esdl_carriers, ensemble_size=ensemble_size, + ensemble=self.__ensemble, ) for ensemble_member_index in range(ensemble_size): self.io.set_parameter("GeothermalSource_fafd.Max_heat", 10e5, ensemble_member_index) diff --git a/src/mesido/esdl/profile_parser.py b/src/mesido/esdl/profile_parser.py index 5fbbe9d2e..89fef25af 100644 --- a/src/mesido/esdl/profile_parser.py +++ b/src/mesido/esdl/profile_parser.py @@ -58,6 +58,7 @@ def read_profiles( esdl_assets: Dict[str, Asset], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble, ) -> None: """ This function takes a datastore and a dictionary of milp network components and loads a @@ -91,6 +92,7 @@ def read_profiles( esdl_asset_id_to_name_map=esdl_asset_id_to_name_map, carrier_properties=carrier_properties, ensemble_size=ensemble_size, + ensemble=ensemble, ) try: @@ -112,8 +114,7 @@ def read_profiles( "power" ] if profile is not None: - correction = 1-ensemble_member * 0.1 - values = profile *correction + values = profile else: if "heat_demand" not in component_type: # We don't set a default profile for source targets @@ -614,6 +615,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble, ) -> None: if self._file_path.suffix == ".xml": logger.warning( @@ -628,6 +630,7 @@ def _load_profiles_from_source( energy_system_components=energy_system_components, carrier_properties=carrier_properties, ensemble_size=ensemble_size, + ensemble=ensemble, ) else: raise _ProfileParserException( @@ -639,27 +642,27 @@ def _load_csv( energy_system_components: Dict[str, Set[str]], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble, ) -> None: - #TODO: this data should be loaded for the different ensemble members if they each have - # their own profile. - data = pd.read_csv(self._file_path) - if len(data.filter(like="Unnamed").columns) > 0: - raise Exception( - f"An unnamed column has been found in profile source file: {self._file_path}" - ) + data_dict = {} + if ensemble_size > 1: + input_folder = self._file_path.parent + file_name = self._file_path.name + for i, ensemble_name, _ in ensemble: + data_dict[i] = pd.read_csv(Path(input_folder /ensemble_name/ file_name)) + else: + data_dict[0] = pd.read_csv(self._file_path) - try: - timeseries_import_times = [ - datetime.datetime.strptime(entry.replace("Z", ""), "%Y-%m-%d %H:%M:%S").replace( - tzinfo=datetime.timezone.utc + for data in data_dict.values(): + if len(data.filter(like="Unnamed").columns) > 0: + raise Exception( + f"An unnamed column has been found in profile source file: {self._file_path}" ) - for entry in data["DateTime"].to_numpy() - ] - except ValueError: + try: timeseries_import_times = [ - datetime.datetime.strptime(entry.replace("Z", ""), "%Y-%m-%dT%H:%M:%S").replace( + datetime.datetime.strptime(entry.replace("Z", ""), "%Y-%m-%d %H:%M:%S").replace( tzinfo=datetime.timezone.utc ) for entry in data["DateTime"].to_numpy() @@ -667,47 +670,56 @@ def _load_csv( except ValueError: try: timeseries_import_times = [ - datetime.datetime.strptime( - entry.replace("Z", ""), "%d-%m-%Y %H:%M" - ).replace(tzinfo=datetime.timezone.utc) + datetime.datetime.strptime(entry.replace("Z", ""), "%Y-%m-%dT%H:%M:%S").replace( + tzinfo=datetime.timezone.utc + ) for entry in data["DateTime"].to_numpy() ] except ValueError: - raise _ProfileParserException("Date time string is not in a supported format") + try: + timeseries_import_times = [ + datetime.datetime.strptime( + entry.replace("Z", ""), "%d-%m-%Y %H:%M" + ).replace(tzinfo=datetime.timezone.utc) + for entry in data["DateTime"].to_numpy() + ] + except ValueError: + raise _ProfileParserException("Date time string is not in a supported format") - logger.warning("Timezone specification not supported yet: default UTC has been used") + logger.warning("Timezone specification not supported yet: default UTC has been used") - self._reference_datetimes = timeseries_import_times + self._reference_datetimes = timeseries_import_times - for ensemble_member in range(ensemble_size): + for e_m in range(ensemble_size): + data_em = data_dict[e_m] for component_type, var_name in self.component_type_to_var_name_map.items(): for component_name in energy_system_components.get(component_type, []): try: column_name = f"{component_name.replace(' ', '')}" - values = data[column_name].to_numpy() + values = data_em[column_name].to_numpy() if np.isnan(values).any(): raise Exception( f"Column name: {column_name}, NaN exists in the profile source" f" file {self._file_path}." - f" Detials: {data[data[column_name].isnull()]}" + f" Detials: {data_em[data_em[column_name].isnull()]}" ) except KeyError: pass else: - self._profiles[ensemble_member][component_name + var_name] = values + self._profiles[e_m][component_name + var_name] = values for properties in carrier_properties.values(): carrier_name = properties.get("name") try: - values = data[carrier_name].to_numpy() + values = data_em[carrier_name].to_numpy() if np.isnan(values).any(): raise Exception( f"Carrier name: {carrier_name}, NaN exists in the profile source file" - f" {self._file_path}. Details: {data[data[carrier_name].isnull()]}" + f" {self._file_path}. Details: {data_em[data_em[carrier_name].isnull()]}" ) except KeyError: pass else: - self._profiles[ensemble_member][ + self._profiles[e_m][ carrier_name + self.carrier_profile_var_name ] = values From a0b80434c2e568b55726f0ef6a60600483e7c853 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 30 Jan 2026 20:03:03 +0100 Subject: [PATCH 03/19] working test case, both for profiles on demand and for a parameter on the max size --- .../input/forecast2/timeseries_import.csv | 6 ++-- .../unit_cases/case_2a_ensemble/model/2a.esdl | 2 +- .../unit_cases/case_2a_ensemble/src/run_2a.py | 29 +++++++++------- tests/test_warmingup_unit_cases.py | 33 +++++++++++++++++-- 4 files changed, 51 insertions(+), 19 deletions(-) diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv index 1e9ee86c7..74e7a9e13 100644 --- a/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv +++ b/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv @@ -1,4 +1,4 @@ DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 -2013-05-19 22:00:00,350000,350000,350000 -2013-05-19 23:00:00,350000,350000,350000 -2013-05-20 00:00:00,350000,350000,350000 +2013-05-19 22:00:00,350000,350000,300000 +2013-05-19 23:00:00,350000,350000,300000 +2013-05-20 00:00:00,350000,350000,300000 diff --git a/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl b/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl index ac27eb46a..1410c2520 100644 --- a/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl +++ b/tests/models/unit_cases/case_2a_ensemble/model/2a.esdl @@ -252,7 +252,7 @@ - + diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index 13e709a59..d2e4a45a6 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -67,7 +67,7 @@ class MinimizeSourcesSizeGoal(Goal): def __init__(self, source): self.target_max = 0.0 - self.function_range = (0.0, 2e6) + self.function_range = (0.0, 10e6) self.source = source self.function_nominal = 1e5 @@ -151,17 +151,6 @@ def goals(self): goals.append(MinimizeSourcesSizeGoal(s)) return goals - def bounds(self): - bounds = super().bounds() - print('a') - return bounds - - def parameters(self, ensemble_member): - parameters = super().parameters(ensemble_member) - if ensemble_member == 1: - parameters['GeothermalSource_fafd.Max_heat']=1e5 - return parameters - def __fixed_max_size(self): constraints=[] for heat_source in self.energy_system_components["heat_source"]: @@ -172,11 +161,27 @@ def __fixed_max_size(self): max_size_prev = max_size return constraints + def __update_target_demand_constraint(self, ensemble_member): + constraints = [] + for demand in self.energy_system_components["heat_demand"]: + var = self.state_vector(f"{demand}.Heat_demand") + nom = self.variable_nominal(f"{demand}.Heat_demand") + target = np.asarray(self.get_timeseries(f"{demand}.target_heat_demand", + ensemble_member).values) + + constraints.append((var*nom/nom, 0.0, target/nom)) + return constraints + def constraints(self, ensemble_member): constraints = super().constraints(ensemble_member).copy() constraints.extend(self.__fixed_max_size()) + constraints.extend(self.__update_target_demand_constraint(ensemble_member)) return constraints + def path_constraints(self, ensemble_member): + constraints = super().path_constraints(ensemble_member).copy() + + return constraints class QTHProblem( diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 0de02841c..75e5d676c 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -209,7 +209,13 @@ def test_2a_ensemble(self): base_folder = Path(run_2a.__file__).resolve().parent.parent - # Just a "problem is not infeasible" + class HeatProblemEnsembleLimit(HeatProblemEnsemble): + def parameters(self, ensemble_member): + parameters = super().parameters(ensemble_member) + if ensemble_member == 1: + parameters['GeothermalSource_fafd.Max_heat']=1e4 + return parameters + heat_problem = run_esdl_mesido_optimization( HeatProblemEnsemble, base_folder=base_folder, @@ -218,15 +224,36 @@ def test_2a_ensemble(self): profile_reader=ProfileReaderFromFile, input_timeseries_file="timeseries_import.csv", ) + heat_problem_limit = run_esdl_mesido_optimization( + HeatProblemEnsembleLimit, + base_folder=base_folder, + esdl_file_name="2a.esdl", + esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="timeseries_import.csv", + ) - results = {} + results, results_limit = {}, {} for ensemble_member in range(heat_problem.ensemble_size): results[ensemble_member] = heat_problem.extract_results(ensemble_member=ensemble_member) + for ensemble_member in range(heat_problem_limit.ensemble_size): + results_limit[ensemble_member] = heat_problem_limit.extract_results(ensemble_member=ensemble_member) - demand_matching_test(heat_problem, heat_problem.extract_results()) + # demand_matching_test(heat_problem, heat_problem.extract_results()) energy_conservation_test(heat_problem, heat_problem.extract_results()) heat_to_discharge_test(heat_problem, heat_problem.extract_results()) + for result in [results, results_limit]: + for ensemble_member in range(heat_problem.ensemble_size): + print("### Ensemble member", ensemble_member) + for demand in heat_problem.energy_system_components.get("heat_demand",[]): + print(demand, result[ensemble_member][f"{demand}.Heat_demand"]) + for prod in heat_problem.energy_system_components.get("heat_source", []): + print(prod, "Heat_source", result[ensemble_member][f"{prod}.Heat_source"]) + print(prod, "Max size", result[ensemble_member][f"{prod}__max_size"]) + # for results_limit, demand in ensemble 2 is now not matched, because the summed heat + # source capacities are insufficient. + if __name__ == "__main__": a = TestWarmingUpUnitCases() From a5eb184c5f786217fe77cc338f8e6bb4264ede52 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Mon, 2 Feb 2026 08:02:50 +0100 Subject: [PATCH 04/19] start on assert statements for demand matching to include ensembles --- tests/test_warmingup_unit_cases.py | 30 +++++++++++------------------- tests/utils_tests.py | 2 +- 2 files changed, 12 insertions(+), 20 deletions(-) diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 75e5d676c..440f1faaa 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -224,33 +224,25 @@ def parameters(self, ensemble_member): profile_reader=ProfileReaderFromFile, input_timeseries_file="timeseries_import.csv", ) - heat_problem_limit = run_esdl_mesido_optimization( - HeatProblemEnsembleLimit, - base_folder=base_folder, - esdl_file_name="2a.esdl", - esdl_parser=ESDLFileParser, - profile_reader=ProfileReaderFromFile, - input_timeseries_file="timeseries_import.csv", - ) results, results_limit = {}, {} - for ensemble_member in range(heat_problem.ensemble_size): - results[ensemble_member] = heat_problem.extract_results(ensemble_member=ensemble_member) - for ensemble_member in range(heat_problem_limit.ensemble_size): - results_limit[ensemble_member] = heat_problem_limit.extract_results(ensemble_member=ensemble_member) + for e_m in range(heat_problem.ensemble_size): + results[e_m] = heat_problem.extract_results(ensemble_member=e_m) - # demand_matching_test(heat_problem, heat_problem.extract_results()) + # for e_m in range(heat_problem.ensemble_size): + # demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) energy_conservation_test(heat_problem, heat_problem.extract_results()) heat_to_discharge_test(heat_problem, heat_problem.extract_results()) - for result in [results, results_limit]: - for ensemble_member in range(heat_problem.ensemble_size): - print("### Ensemble member", ensemble_member) + for result in [results]: + for e_m in range(heat_problem.ensemble_size): + print("### Ensemble member", e_m) for demand in heat_problem.energy_system_components.get("heat_demand",[]): - print(demand, result[ensemble_member][f"{demand}.Heat_demand"]) + print(demand, result[e_m][f"{demand}.Heat_demand"]) + print(demand, heat_problem.get_timeseries(f"{demand}.target_heat_demand", e_m)) for prod in heat_problem.energy_system_components.get("heat_source", []): - print(prod, "Heat_source", result[ensemble_member][f"{prod}.Heat_source"]) - print(prod, "Max size", result[ensemble_member][f"{prod}__max_size"]) + print(prod, "Heat_source", result[e_m][f"{prod}.Heat_source"]) + print(prod, "Max size", result[e_m][f"{prod}__max_size"]) # for results_limit, demand in ensemble 2 is now not matched, because the summed heat # source capacities are insufficient. diff --git a/tests/utils_tests.py b/tests/utils_tests.py index b42118929..109ec434f 100644 --- a/tests/utils_tests.py +++ b/tests/utils_tests.py @@ -21,7 +21,7 @@ def feasibility_test(solution): ) -def demand_matching_test(solution, results, atol=1.0e-3, rtol=1.0e-6): +def demand_matching_test(solution, results, ensemble_member=0, atol=1.0e-3, rtol=1.0e-6): """ "Test function to check whether the milp demand of each consumer is matched""" for d in solution.energy_system_components.get("heat_demand", []): if len(solution.times()) > 0: From a3628ba89a665b808216fe770040b0b037817f3d Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Mon, 2 Feb 2026 09:25:02 +0100 Subject: [PATCH 05/19] updates to test case with working demand matching timeseries for ensembles --- .../unit_cases/case_2a_ensemble/src/run_2a.py | 31 ++++++++++++++++--- tests/test_warmingup_unit_cases.py | 4 +-- tests/utils_tests.py | 9 +++--- 3 files changed, 33 insertions(+), 11 deletions(-) diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index d2e4a45a6..45b7497eb 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -41,8 +41,28 @@ def function(self, optimization_problem, ensemble_member): to :func:`OptimizationProblem.state`. This call returns a time-independent symbol that is also independent of the active ensemble member. Path goals are applied to all times and all ensemble members simultaneously.""" - self.target_min = self.targets[ensemble_member] - self.target_max = self.targets[ensemble_member] + + return optimization_problem.state(self.state) + +class TargetDemandPathGoal(Goal): + priority = 1 + + order = 2 + + def __init__(self, state, target): + self.state = state + + self.target_min = target + self.target_max = target + self.function_range = (0.0, 2.0 * max(target.values)) + self.function_nominal = np.median(target.values) + + def function(self, optimization_problem, ensemble_member): + """Note that for path goals, the ensemble member index is not passed to the call + to :func:`OptimizationProblem.state`. This call returns a time-independent symbol + that is also independent of the active ensemble member. Path goals are + applied to all times and all ensemble members simultaneously.""" + return optimization_problem.state(self.state) @@ -53,7 +73,7 @@ class MinimizeSourcesHeatGoal(Goal): def __init__(self, source): self.target_max = 0.0 - self.function_range = (0.0, 2e6) + self.function_range = (0.0, 10e6) self.source = source self.function_nominal = 1e5 @@ -97,8 +117,9 @@ def path_goals(self): ensemble_member=ensemble_member) for ensemble_member in range(self.ensemble_size)] state = f"{demand}.Heat_demand" + target = self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member=0) - goals.append(TargetDemandGoal(state, target)) + goals.append(TargetDemandPathGoal(state, target)) return goals @@ -164,7 +185,7 @@ def __fixed_max_size(self): def __update_target_demand_constraint(self, ensemble_member): constraints = [] for demand in self.energy_system_components["heat_demand"]: - var = self.state_vector(f"{demand}.Heat_demand") + var = self.state_vector(f"{demand}.Heat_demand", ensemble_member=ensemble_member) nom = self.variable_nominal(f"{demand}.Heat_demand") target = np.asarray(self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member).values) diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 440f1faaa..3d5bdb22c 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -229,8 +229,8 @@ def parameters(self, ensemble_member): for e_m in range(heat_problem.ensemble_size): results[e_m] = heat_problem.extract_results(ensemble_member=e_m) - # for e_m in range(heat_problem.ensemble_size): - # demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) + for e_m in range(heat_problem.ensemble_size): + demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) energy_conservation_test(heat_problem, heat_problem.extract_results()) heat_to_discharge_test(heat_problem, heat_problem.extract_results()) diff --git a/tests/utils_tests.py b/tests/utils_tests.py index 109ec434f..24af4d922 100644 --- a/tests/utils_tests.py +++ b/tests/utils_tests.py @@ -28,14 +28,15 @@ def demand_matching_test(solution, results, ensemble_member=0, atol=1.0e-3, rtol len_times = len(solution.times()) else: len_times = len(solution.get_timeseries(f"{d}.target_heat_demand").values) - target = solution.get_timeseries(f"{d}.target_heat_demand").values[0:len_times] + target = solution.get_timeseries(f"{d}.target_heat_demand", ensemble_member).values[ + 0:len_times] np.testing.assert_allclose(target, results[f"{d}.Heat_demand"], atol=atol, rtol=rtol) for d in solution.energy_system_components.get("cold_demand", []): if len(solution.times()) > 0: len_times = len(solution.times()) else: len_times = len(solution.get_timeseries(f"{d}.target_cold_demand").values) - target = solution.get_timeseries(f"{d}.target_cold_demand").values[0:len_times] + target = solution.get_timeseries(f"{d}.target_cold_demand", ensemble_member).values[0:len_times] np.testing.assert_allclose(target, results[f"{d}.Cold_demand"], atol=atol, rtol=rtol) for d in solution.energy_system_components.get("gas_demand", []): timeseries_name = f"{d}.target_gas_demand" @@ -44,7 +45,7 @@ def demand_matching_test(solution, results, ensemble_member=0, atol=1.0e-3, rtol len_times = len(solution.times()) else: len_times = len(solution.get_timeseries(timeseries_name).values) - target = solution.get_timeseries(timeseries_name).values[0:len_times] + target = solution.get_timeseries(timeseries_name, ensemble_member).values[0:len_times] np.testing.assert_allclose( target, results[f"{d}.Gas_demand_mass_flow"], atol=atol, rtol=rtol ) @@ -55,7 +56,7 @@ def demand_matching_test(solution, results, ensemble_member=0, atol=1.0e-3, rtol len_times = len(solution.times()) else: len_times = len(solution.get_timeseries(timeseries_name).values) - target = solution.get_timeseries(timeseries_name).values[0:len_times] + target = solution.get_timeseries(timeseries_name, ensemble_member).values[0:len_times] np.testing.assert_allclose( target, results[f"{d}.Electricity_demand"], atol=atol, rtol=rtol ) From d6ec6d36dd920bbde90b231c04b471731fdc931d Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Mon, 2 Feb 2026 17:06:18 +0100 Subject: [PATCH 06/19] wip for generalizing both ensemble profile methods --- .../unit_cases/case_2a_ensemble/src/run_2a.py | 66 +++++++++++++------ tests/test_warmingup_unit_cases.py | 11 +--- 2 files changed, 48 insertions(+), 29 deletions(-) diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index 45b7497eb..ae56c7a1e 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -25,24 +25,30 @@ class TargetDemandGoal(Goal): priority = 1 - order = 2 + order = 1 - def __init__(self, state, target): + def __init__(self, state, target, index): self.state = state + self.index = index self.targets = target - self.target_min = target[0] - self.target_max = target[0] - self.function_range = (0.0, 2.0 * max(target[0].values)) - self.function_nominal = np.median(target[0].values) + # self.target_min = target[0].values[index] + # self.target_max = target[0].values[index] + # self.function_range = (0.0, 2.0 * max(target[0].values)) + self.function_nominal = np.median(target[0].values[index]) def function(self, optimization_problem, ensemble_member): """Note that for path goals, the ensemble member index is not passed to the call to :func:`OptimizationProblem.state`. This call returns a time-independent symbol that is also independent of the active ensemble member. Path goals are applied to all times and all ensemble members simultaneously.""" + nom = optimization_problem.variable_nominal(self.state) + vector_state = nom * optimization_problem.state_vector(self.state, ensemble_member)[ + self.index] + # self.target_min = self.targets[ensemble_member].values[self.index] + # self.target_max = self.targets[ensemble_member].values[self.index] - return optimization_problem.state(self.state) + return self.targets[ensemble_member].values[self.index] - vector_state class TargetDemandPathGoal(Goal): priority = 1 @@ -58,10 +64,10 @@ def __init__(self, state, target): self.function_nominal = np.median(target.values) def function(self, optimization_problem, ensemble_member): - """Note that for path goals, the ensemble member index is not passed to the call - to :func:`OptimizationProblem.state`. This call returns a time-independent symbol - that is also independent of the active ensemble member. Path goals are - applied to all times and all ensemble members simultaneously.""" + """Path goals are applied to all timeseries and ensemble members simultaneously, + therefore the goals ensemble independent and only one target for all the ensembles can be + applied.""" + return optimization_problem.state(self.state) @@ -109,17 +115,31 @@ def function(self, optimization_problem, ensemble_member): class _GoalsAndOptions: - def path_goals(self): - goals = super().path_goals().copy() + # def path_goals(self): + # goals = super().path_goals().copy() + # + # for demand in self.energy_system_components["heat_demand"]: + # target = [self.get_timeseries(f"{demand}.target_heat_demand", + # ensemble_member=ensemble_member) for ensemble_member in + # range(self.ensemble_size)] + # state = f"{demand}.Heat_demand" + # target = self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member=0) + # + # goals.append(TargetDemandPathGoal(state, target)) + # + # return goals - for demand in self.energy_system_components["heat_demand"]: - target = [self.get_timeseries(f"{demand}.target_heat_demand", - ensemble_member=ensemble_member) for ensemble_member in - range(self.ensemble_size)] - state = f"{demand}.Heat_demand" - target = self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member=0) + def goals(self): + goals = super().goals().copy() - goals.append(TargetDemandPathGoal(state, target)) + for i in range(len(self.times())): + for demand in self.energy_system_components["heat_demand"]: + target = [self.get_timeseries(f"{demand}.target_heat_demand", + ensemble_member=ensemble_member) for ensemble_member in + range(self.ensemble_size)] + state = f"{demand}.Heat_demand" + + goals.append(TargetDemandGoal(state, target, i)) return goals @@ -170,6 +190,7 @@ def goals(self): goals = super().goals().copy() for s in self.energy_system_components["heat_source"]: goals.append(MinimizeSourcesSizeGoal(s)) + return goals def __fixed_max_size(self): @@ -183,6 +204,11 @@ def __fixed_max_size(self): return constraints def __update_target_demand_constraint(self, ensemble_member): + """ + This constraint method adds upper limits to the demand heat production based on the + target of a specific ensemble. This requires the first ensemble to always have the + timeseries with the largest values. + """ constraints = [] for demand in self.energy_system_components["heat_demand"]: var = self.state_vector(f"{demand}.Heat_demand", ensemble_member=ensemble_member) diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 3d5bdb22c..91943f116 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -209,13 +209,6 @@ def test_2a_ensemble(self): base_folder = Path(run_2a.__file__).resolve().parent.parent - class HeatProblemEnsembleLimit(HeatProblemEnsemble): - def parameters(self, ensemble_member): - parameters = super().parameters(ensemble_member) - if ensemble_member == 1: - parameters['GeothermalSource_fafd.Max_heat']=1e4 - return parameters - heat_problem = run_esdl_mesido_optimization( HeatProblemEnsemble, base_folder=base_folder, @@ -229,8 +222,8 @@ def parameters(self, ensemble_member): for e_m in range(heat_problem.ensemble_size): results[e_m] = heat_problem.extract_results(ensemble_member=e_m) - for e_m in range(heat_problem.ensemble_size): - demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) + # for e_m in range(heat_problem.ensemble_size): + # demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) energy_conservation_test(heat_problem, heat_problem.extract_results()) heat_to_discharge_test(heat_problem, heat_problem.extract_results()) From 0882f55eec8bc557204e214b42569c110b199f32 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 11 Feb 2026 16:30:37 +0100 Subject: [PATCH 07/19] working test case for ensemble timeseries parsing --- .../input/timeseries_import.csv | 4 -- .../input/timeseries_import.xml | 52 ------------------- tests/test_profile_parsing.py | 45 ++++++++++++++++ 3 files changed, 45 insertions(+), 56 deletions(-) delete mode 100644 tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv delete mode 100644 tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml diff --git a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv deleted file mode 100644 index 1e9ee86c7..000000000 --- a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.csv +++ /dev/null @@ -1,4 +0,0 @@ -DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 -2013-05-19 22:00:00,350000,350000,350000 -2013-05-19 23:00:00,350000,350000,350000 -2013-05-20 00:00:00,350000,350000,350000 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml b/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml deleted file mode 100644 index b1defa019..000000000 --- a/tests/models/unit_cases/case_2a_ensemble/input/timeseries_import.xml +++ /dev/null @@ -1,52 +0,0 @@ - - - 0.0 - -
- instantaneous - - 74846440-61d9-4f57-9bce-77e3664ba832 - b145eba8-612d-42b7-bea7-7e22a2440509 - - - - -999.0 - W -
- - - - - -
- instantaneous - - c6c8d8f7-b6ae-4e8b-be9a-72282a5c3f1e - 874dc021-09f4-46ec-bb5a-c3276916fbfd - - - - -999.0 - W -
- - - - - -
- instantaneous - - 6f99424a-004e-4fb4-8edb-17c374d5be6b - 30713528-878b-46fd-82e2-9da8eb95075e - - - - -999.0 - W -
- - - - - diff --git a/tests/test_profile_parsing.py b/tests/test_profile_parsing.py index fd7439185..1c3734f36 100644 --- a/tests/test_profile_parsing.py +++ b/tests/test_profile_parsing.py @@ -315,6 +315,51 @@ def test_loading_profile_from_esdl(self): atol=1e-2, ) + def test_loading_profiles_ensemble_members(self): + """ + This is the most basic check where we have a simple network and check for the basic physics. + This simple network includes two source, pipes, nodes, and 3 demands. + + Checks; + - Demand matching + - Energy conservation + - Heat to discharge + + """ + import models.unit_cases.case_2a_ensemble.src.run_2a as run_2a + from models.unit_cases.case_2a_ensemble.src.run_2a import HeatProblemEnsemble + + base_folder = Path(run_2a.__file__).resolve().parent.parent + model_folder = base_folder / "model" + input_folder = base_folder / "input" + + problem = HeatProblemEnsemble( + base_folder=base_folder, + model_folder=model_folder, + input_folder=input_folder, + esdl_file_name="2a.esdl", + esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="timeseries_import.csv", + ) + + problem.pre() + + #check that the ensemble size is set at 2, which is based on the ensemble.csv + np.testing.assert_equal(problem.ensemble_size, 2) + + # check that the timeseries are loaded for all ensemble sizes and that the timeseries are + # equal to a heating demand of 350000 except for HeatingDemand_6f99 at the second + # ensemble, where it is equal to 300000. + timeseries_names = problem.io.get_timeseries_names() + for t_name in timeseries_names: + for e_m in range(problem.ensemble_size): + t_series = problem.get_timeseries(t_name, e_m) + if e_m == 1 and t_name == 'HeatingDemand_6f99.target_heat_demand': + np.testing.assert_allclose(t_series.values, [300000]*3) + else: + np.testing.assert_allclose(t_series.values, [350000]*3) + if __name__ == "__main__": # unittest.main() From 4ec39922293a4e2b746ae1822ef69609e955b15f Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 11 Feb 2026 16:34:46 +0100 Subject: [PATCH 08/19] style --- src/mesido/esdl/esdl_mixin.py | 8 ++- src/mesido/esdl/profile_parser.py | 19 +++---- src/mesido/heat_physics_mixin.py | 4 +- .../unit_cases/case_2a_ensemble/src/run_2a.py | 49 +++++++++++-------- tests/test_profile_parsing.py | 8 +-- tests/test_warmingup_unit_cases.py | 45 ----------------- tests/utils_tests.py | 7 ++- 7 files changed, 52 insertions(+), 88 deletions(-) diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index 0d66faf1a..556a85944 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -9,7 +9,6 @@ from typing import Any, Dict, List, Optional, Tuple import esdl.esdl_handler -from rtctools.data import csv from mesido.component_type_mixin import ( ModelicaComponentTypeMixin, @@ -695,8 +694,8 @@ def read(self) -> None: super().read() # TODO: update for ensembles, eg ensemble_size is now still 1, if io gets updated, # then the size also updates automatically - ensemble_size=1 - self.__ensemble=None + ensemble_size = 1 + self.__ensemble = None if self.csv_ensemble_mode: self.__ensemble = np.genfromtxt( os.path.join(self._input_folder, "ensemble.csv"), @@ -724,8 +723,7 @@ def read(self) -> None: for ensemble_member_index in range(ensemble_size): self.io.set_parameter("GeothermalSource_fafd.Max_heat", 10e5, ensemble_member_index) if ensemble_member_index == 1: - self.io.set_parameter('GeothermalSource_fafd.Max_heat', 2e5, - ensemble_member_index) + self.io.set_parameter("GeothermalSource_fafd.Max_heat", 2e5, ensemble_member_index) def write(self) -> None: """ diff --git a/src/mesido/esdl/profile_parser.py b/src/mesido/esdl/profile_parser.py index 89fef25af..8e82177c6 100644 --- a/src/mesido/esdl/profile_parser.py +++ b/src/mesido/esdl/profile_parser.py @@ -650,7 +650,7 @@ def _load_csv( input_folder = self._file_path.parent file_name = self._file_path.name for i, ensemble_name, _ in ensemble: - data_dict[i] = pd.read_csv(Path(input_folder /ensemble_name/ file_name)) + data_dict[i] = pd.read_csv(Path(input_folder / ensemble_name / file_name)) else: data_dict[0] = pd.read_csv(self._file_path) @@ -670,9 +670,9 @@ def _load_csv( except ValueError: try: timeseries_import_times = [ - datetime.datetime.strptime(entry.replace("Z", ""), "%Y-%m-%dT%H:%M:%S").replace( - tzinfo=datetime.timezone.utc - ) + datetime.datetime.strptime( + entry.replace("Z", ""), "%Y-%m-%dT%H:%M:%S" + ).replace(tzinfo=datetime.timezone.utc) for entry in data["DateTime"].to_numpy() ] except ValueError: @@ -684,7 +684,9 @@ def _load_csv( for entry in data["DateTime"].to_numpy() ] except ValueError: - raise _ProfileParserException("Date time string is not in a supported format") + raise _ProfileParserException( + "Date time string is not in a supported format" + ) logger.warning("Timezone specification not supported yet: default UTC has been used") @@ -714,14 +716,13 @@ def _load_csv( if np.isnan(values).any(): raise Exception( f"Carrier name: {carrier_name}, NaN exists in the profile source file" - f" {self._file_path}. Details: {data_em[data_em[carrier_name].isnull()]}" + f" {self._file_path}. Details: " + f"{data_em[data_em[carrier_name].isnull()]}" ) except KeyError: pass else: - self._profiles[e_m][ - carrier_name + self.carrier_profile_var_name - ] = values + self._profiles[e_m][carrier_name + self.carrier_profile_var_name] = values def _load_xml(self, energy_system_components, esdl_asset_id_to_name_map): timeseries_import_basename = self._file_path.stem diff --git a/src/mesido/heat_physics_mixin.py b/src/mesido/heat_physics_mixin.py index a12e2f655..ac94fa650 100644 --- a/src/mesido/heat_physics_mixin.py +++ b/src/mesido/heat_physics_mixin.py @@ -1506,12 +1506,12 @@ def __source_heat_to_discharge_path_constraints(self, ensemble_member): discharge = self.state(f"{s}.Q") heat_out = self.state(f"{s}.HeatOut.Heat") - #somehow the ensemblemember dependent path constraint is only allowed in the bounds of + # somehow the ensemblemember dependent path constraint is only allowed in the bounds of # the problem, e.g. the upper or lower bound. For non path constraints, it is allowed # in the equation. heat_flow = self.state(f"{s}.Heat_flow") max_heat = parameters[f"{s}.Max_heat"] - constraints.append((heat_flow/heat_nominal, -np.inf, max_heat/heat_nominal)) + constraints.append((heat_flow / heat_nominal, -np.inf, max_heat / heat_nominal)) constraint_nominal = (heat_nominal * cp * rho * dt * q_nominal) ** 0.5 diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index ae56c7a1e..a982179e9 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -1,6 +1,3 @@ -from rtctools.optimization.control_tree_mixin import ControlTreeMixin -from rtctools.optimization.csv_mixin import CSVMixin - from mesido.asset_sizing_mixin import AssetSizingMixin from mesido.esdl.esdl_mixin import ESDLMixin from mesido.esdl.esdl_parser import ESDLFileParser @@ -13,6 +10,7 @@ from rtctools.optimization.collocated_integrated_optimization_problem import ( CollocatedIntegratedOptimizationProblem, ) +from rtctools.optimization.control_tree_mixin import ControlTreeMixin from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin from rtctools.optimization.homotopy_mixin import HomotopyMixin from rtctools.optimization.linearized_order_goal_programming_mixin import ( @@ -43,13 +41,15 @@ def function(self, optimization_problem, ensemble_member): that is also independent of the active ensemble member. Path goals are applied to all times and all ensemble members simultaneously.""" nom = optimization_problem.variable_nominal(self.state) - vector_state = nom * optimization_problem.state_vector(self.state, ensemble_member)[ - self.index] + vector_state = ( + nom * optimization_problem.state_vector(self.state, ensemble_member)[self.index] + ) # self.target_min = self.targets[ensemble_member].values[self.index] # self.target_max = self.targets[ensemble_member].values[self.index] return self.targets[ensemble_member].values[self.index] - vector_state + class TargetDemandPathGoal(Goal): priority = 1 @@ -68,7 +68,6 @@ def function(self, optimization_problem, ensemble_member): therefore the goals ensemble independent and only one target for all the ensembles can be applied.""" - return optimization_problem.state(self.state) @@ -86,6 +85,7 @@ def __init__(self, source): def function(self, optimization_problem, ensemble_member): return optimization_problem.state(f"{self.source}.Heat_source") + class MinimizeSourcesSizeGoal(Goal): priority = 2 @@ -134,9 +134,12 @@ def goals(self): for i in range(len(self.times())): for demand in self.energy_system_components["heat_demand"]: - target = [self.get_timeseries(f"{demand}.target_heat_demand", - ensemble_member=ensemble_member) for ensemble_member in - range(self.ensemble_size)] + target = [ + self.get_timeseries( + f"{demand}.target_heat_demand", ensemble_member=ensemble_member + ) + for ensemble_member in range(self.ensemble_size) + ] state = f"{demand}.Heat_demand" goals.append(TargetDemandGoal(state, target, i)) @@ -173,18 +176,21 @@ def solver_options(self): return options -class HeatProblemEnsemble(AssetSizingMixin, - HeatProblem, - # CSVMixin, - ControlTreeMixin): +class HeatProblemEnsemble( + AssetSizingMixin, + HeatProblem, + # CSVMixin, + ControlTreeMixin, +): csv_ensemble_mode = True def pre(self): super().pre() # Empty dict for intermediate ensemble results - self.intermediate_results = {ensemble_member: [] for ensemble_member in range( - self.ensemble_size)} + self.intermediate_results = { + ensemble_member: [] for ensemble_member in range(self.ensemble_size) + } def goals(self): goals = super().goals().copy() @@ -194,11 +200,11 @@ def goals(self): return goals def __fixed_max_size(self): - constraints=[] + constraints = [] for heat_source in self.energy_system_components["heat_source"]: max_size_prev = self.extra_variable(f"{heat_source}__max_size", ensemble_member=0) - for e_m in range(self.ensemble_size-1): - max_size = self.extra_variable(f"{heat_source}__max_size", ensemble_member=e_m+1) + for e_m in range(self.ensemble_size - 1): + max_size = self.extra_variable(f"{heat_source}__max_size", ensemble_member=e_m + 1) constraints.append((max_size - max_size_prev, 0.0, 0.0)) max_size_prev = max_size return constraints @@ -213,10 +219,11 @@ def __update_target_demand_constraint(self, ensemble_member): for demand in self.energy_system_components["heat_demand"]: var = self.state_vector(f"{demand}.Heat_demand", ensemble_member=ensemble_member) nom = self.variable_nominal(f"{demand}.Heat_demand") - target = np.asarray(self.get_timeseries(f"{demand}.target_heat_demand", - ensemble_member).values) + target = np.asarray( + self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member).values + ) - constraints.append((var*nom/nom, 0.0, target/nom)) + constraints.append((var * nom / nom, 0.0, target / nom)) return constraints def constraints(self, ensemble_member): diff --git a/tests/test_profile_parsing.py b/tests/test_profile_parsing.py index 1c3734f36..513d319d4 100644 --- a/tests/test_profile_parsing.py +++ b/tests/test_profile_parsing.py @@ -345,7 +345,7 @@ def test_loading_profiles_ensemble_members(self): problem.pre() - #check that the ensemble size is set at 2, which is based on the ensemble.csv + # check that the ensemble size is set at 2, which is based on the ensemble.csv np.testing.assert_equal(problem.ensemble_size, 2) # check that the timeseries are loaded for all ensemble sizes and that the timeseries are @@ -355,10 +355,10 @@ def test_loading_profiles_ensemble_members(self): for t_name in timeseries_names: for e_m in range(problem.ensemble_size): t_series = problem.get_timeseries(t_name, e_m) - if e_m == 1 and t_name == 'HeatingDemand_6f99.target_heat_demand': - np.testing.assert_allclose(t_series.values, [300000]*3) + if e_m == 1 and t_name == "HeatingDemand_6f99.target_heat_demand": + np.testing.assert_allclose(t_series.values, [300000] * 3) else: - np.testing.assert_allclose(t_series.values, [350000]*3) + np.testing.assert_allclose(t_series.values, [350000] * 3) if __name__ == "__main__": diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py index 91943f116..3a53f1665 100644 --- a/tests/test_warmingup_unit_cases.py +++ b/tests/test_warmingup_unit_cases.py @@ -193,51 +193,6 @@ def test_3a(self): atol=1.0e-6, ) - def test_2a_ensemble(self): - """ - This is the most basic check where we have a simple network and check for the basic physics. - This simple network includes two source, pipes, nodes, and 3 demands. - - Checks; - - Demand matching - - Energy conservation - - Heat to discharge - - """ - import models.unit_cases.case_2a_ensemble.src.run_2a as run_2a - from models.unit_cases.case_2a_ensemble.src.run_2a import HeatProblemEnsemble - - base_folder = Path(run_2a.__file__).resolve().parent.parent - - heat_problem = run_esdl_mesido_optimization( - HeatProblemEnsemble, - base_folder=base_folder, - esdl_file_name="2a.esdl", - esdl_parser=ESDLFileParser, - profile_reader=ProfileReaderFromFile, - input_timeseries_file="timeseries_import.csv", - ) - - results, results_limit = {}, {} - for e_m in range(heat_problem.ensemble_size): - results[e_m] = heat_problem.extract_results(ensemble_member=e_m) - - # for e_m in range(heat_problem.ensemble_size): - # demand_matching_test(heat_problem, results[e_m], ensemble_member=e_m) - energy_conservation_test(heat_problem, heat_problem.extract_results()) - heat_to_discharge_test(heat_problem, heat_problem.extract_results()) - - for result in [results]: - for e_m in range(heat_problem.ensemble_size): - print("### Ensemble member", e_m) - for demand in heat_problem.energy_system_components.get("heat_demand",[]): - print(demand, result[e_m][f"{demand}.Heat_demand"]) - print(demand, heat_problem.get_timeseries(f"{demand}.target_heat_demand", e_m)) - for prod in heat_problem.energy_system_components.get("heat_source", []): - print(prod, "Heat_source", result[e_m][f"{prod}.Heat_source"]) - print(prod, "Max size", result[e_m][f"{prod}__max_size"]) - # for results_limit, demand in ensemble 2 is now not matched, because the summed heat - # source capacities are insufficient. if __name__ == "__main__": diff --git a/tests/utils_tests.py b/tests/utils_tests.py index 24af4d922..3946b4b8b 100644 --- a/tests/utils_tests.py +++ b/tests/utils_tests.py @@ -29,14 +29,17 @@ def demand_matching_test(solution, results, ensemble_member=0, atol=1.0e-3, rtol else: len_times = len(solution.get_timeseries(f"{d}.target_heat_demand").values) target = solution.get_timeseries(f"{d}.target_heat_demand", ensemble_member).values[ - 0:len_times] + 0:len_times + ] np.testing.assert_allclose(target, results[f"{d}.Heat_demand"], atol=atol, rtol=rtol) for d in solution.energy_system_components.get("cold_demand", []): if len(solution.times()) > 0: len_times = len(solution.times()) else: len_times = len(solution.get_timeseries(f"{d}.target_cold_demand").values) - target = solution.get_timeseries(f"{d}.target_cold_demand", ensemble_member).values[0:len_times] + target = solution.get_timeseries(f"{d}.target_cold_demand", ensemble_member).values[ + 0:len_times + ] np.testing.assert_allclose(target, results[f"{d}.Cold_demand"], atol=atol, rtol=rtol) for d in solution.energy_system_components.get("gas_demand", []): timeseries_name = f"{d}.target_gas_demand" From fac0fed1867a227aed38cc065320729b0f6bd078 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 11 Feb 2026 16:47:50 +0100 Subject: [PATCH 09/19] clean up --- src/mesido/esdl/esdl_heat_model.py | 1 - src/mesido/heat_physics_mixin.py | 13 ------------- .../component_library/milp/heat/heat_source.py | 1 - tests/test_profile_parsing.py | 4 ++++ 4 files changed, 4 insertions(+), 15 deletions(-) diff --git a/src/mesido/esdl/esdl_heat_model.py b/src/mesido/esdl/esdl_heat_model.py index 89f3cedae..88c36417c 100644 --- a/src/mesido/esdl/esdl_heat_model.py +++ b/src/mesido/esdl/esdl_heat_model.py @@ -1280,7 +1280,6 @@ def convert_heat_source(self, asset: Asset) -> Tuple[Type[HeatSource], MODIFIERS modifiers = dict( Heat_source=dict(min=0.0, max=max_supply, nominal=max_supply / 2.0), - Max_heat=max_supply, **self._generic_modifiers(asset), **self._generic_heat_modifiers(0.0, max_supply, q_nominal), **self._supply_return_temperature_modifiers(asset), diff --git a/src/mesido/heat_physics_mixin.py b/src/mesido/heat_physics_mixin.py index ac94fa650..43441cfd8 100644 --- a/src/mesido/heat_physics_mixin.py +++ b/src/mesido/heat_physics_mixin.py @@ -1506,12 +1506,6 @@ def __source_heat_to_discharge_path_constraints(self, ensemble_member): discharge = self.state(f"{s}.Q") heat_out = self.state(f"{s}.HeatOut.Heat") - # somehow the ensemblemember dependent path constraint is only allowed in the bounds of - # the problem, e.g. the upper or lower bound. For non path constraints, it is allowed - # in the equation. - heat_flow = self.state(f"{s}.Heat_flow") - max_heat = parameters[f"{s}.Max_heat"] - constraints.append((heat_flow / heat_nominal, -np.inf, max_heat / heat_nominal)) constraint_nominal = (heat_nominal * cp * rho * dt * q_nominal) ** 0.5 @@ -3761,13 +3755,6 @@ def constraints(self, ensemble_member): """ constraints = super().constraints(ensemble_member) - # parameters = self.parameters(ensemble_member) - # for s in self.energy_system_components.get("heat_source", []): - # heat_nominal = parameters[f"{s}.Heat_nominal"] - # heat_flow = self.__state_vector_scaled(f"{s}.Heat_flow", ensemble_member) - # max_heat = parameters[f"{s}.Max_heat"] - # constraints.append(((heat_flow-max_heat)/heat_nominal, -np.inf, 0.0)) - if self.heat_network_settings["head_loss_option"] != HeadLossOption.NO_HEADLOSS: # constraints.extend(self._hn_pipe_head_loss_constraints(ensemble_member)) constraints.extend( diff --git a/src/mesido/pycml/component_library/milp/heat/heat_source.py b/src/mesido/pycml/component_library/milp/heat/heat_source.py index c1ecdb8b3..634c35f5f 100644 --- a/src/mesido/pycml/component_library/milp/heat/heat_source.py +++ b/src/mesido/pycml/component_library/milp/heat/heat_source.py @@ -44,7 +44,6 @@ def __init__(self, name, **modifiers): self.price = nan # TODO: delete not needed anymore self.co2_coeff = 1.0 self.pump_efficiency = 0.5 - self.Max_heat = nan # Assumption: heat in/out and added is nonnegative # Heat in the return (i.e. cold) line is zero diff --git a/tests/test_profile_parsing.py b/tests/test_profile_parsing.py index 513d319d4..59ce9386b 100644 --- a/tests/test_profile_parsing.py +++ b/tests/test_profile_parsing.py @@ -347,6 +347,10 @@ def test_loading_profiles_ensemble_members(self): # check that the ensemble size is set at 2, which is based on the ensemble.csv np.testing.assert_equal(problem.ensemble_size, 2) + prob_0 = problem.ensemble_member_probability(0) + prob_1 = problem.ensemble_member_probability(1) + np.testing.assert_allclose(prob_0, 0.7) + np.testing.assert_allclose(prob_1, 0.3) # check that the timeseries are loaded for all ensemble sizes and that the timeseries are # equal to a heating demand of 350000 except for HeatingDemand_6f99 at the second From 3b4b547e76d99cdf8c9805e0f34b5668fba803d5 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 11 Feb 2026 17:00:56 +0100 Subject: [PATCH 10/19] added description --- .../unit_cases/case_2a_ensemble/src/run_2a.py | 64 ------------------- tests/test_profile_parsing.py | 14 ++-- 2 files changed, 6 insertions(+), 72 deletions(-) diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index a982179e9..0b75b08e4 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -3,7 +3,6 @@ from mesido.esdl.esdl_parser import ESDLFileParser from mesido.esdl.profile_parser import ProfileReaderFromFile from mesido.physics_mixin import PhysicsMixin -from mesido.qth_not_maintained.qth_mixin import QTHMixin import numpy as np @@ -12,11 +11,9 @@ ) from rtctools.optimization.control_tree_mixin import ControlTreeMixin from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin -from rtctools.optimization.homotopy_mixin import HomotopyMixin from rtctools.optimization.linearized_order_goal_programming_mixin import ( LinearizedOrderGoalProgrammingMixin, ) -from rtctools.optimization.single_pass_goal_programming_mixin import SinglePassGoalProgrammingMixin from rtctools.util import run_optimization_problem @@ -71,21 +68,6 @@ def function(self, optimization_problem, ensemble_member): return optimization_problem.state(self.state) -class MinimizeSourcesHeatGoal(Goal): - priority = 2 - - order = 2 - - def __init__(self, source): - self.target_max = 0.0 - self.function_range = (0.0, 10e6) - self.source = source - self.function_nominal = 1e5 - - def function(self, optimization_problem, ensemble_member): - return optimization_problem.state(f"{self.source}.Heat_source") - - class MinimizeSourcesSizeGoal(Goal): priority = 2 @@ -101,19 +83,6 @@ def function(self, optimization_problem, ensemble_member): return optimization_problem.extra_variable(f"{self.source}__max_size", ensemble_member) -class MinimizeSourcesQTHGoal(Goal): - priority = 2 - - order = 2 - - def __init__(self, source): - self.source = source - self.function_nominal = 1e5 - - def function(self, optimization_problem, ensemble_member): - return optimization_problem.state(f"{self.source}.Heat_source") - - class _GoalsAndOptions: # def path_goals(self): # goals = super().path_goals().copy() @@ -155,13 +124,6 @@ class HeatProblem( ESDLMixin, CollocatedIntegratedOptimizationProblem, ): - # def path_goals(self): - # goals = super().path_goals().copy() - # - # for s in self.energy_system_components["heat_source"]: - # goals.append(MinimizeSourcesHeatGoal(s)) - # - # return goals def energy_system_options(self): options = super().energy_system_options() @@ -179,19 +141,10 @@ def solver_options(self): class HeatProblemEnsemble( AssetSizingMixin, HeatProblem, - # CSVMixin, ControlTreeMixin, ): csv_ensemble_mode = True - def pre(self): - super().pre() - - # Empty dict for intermediate ensemble results - self.intermediate_results = { - ensemble_member: [] for ensemble_member in range(self.ensemble_size) - } - def goals(self): goals = super().goals().copy() for s in self.energy_system_components["heat_source"]: @@ -238,23 +191,6 @@ def path_constraints(self, ensemble_member): return constraints -class QTHProblem( - _GoalsAndOptions, - QTHMixin, - HomotopyMixin, - SinglePassGoalProgrammingMixin, - ESDLMixin, - CollocatedIntegratedOptimizationProblem, -): - def path_goals(self): - goals = super().path_goals().copy() - - for s in self.energy_system_components["heat_source"]: - goals.append(MinimizeSourcesQTHGoal(s)) - - return goals - - if __name__ == "__main__": run_optimization_problem( HeatProblem, diff --git a/tests/test_profile_parsing.py b/tests/test_profile_parsing.py index 59ce9386b..0184657c6 100644 --- a/tests/test_profile_parsing.py +++ b/tests/test_profile_parsing.py @@ -317,14 +317,12 @@ def test_loading_profile_from_esdl(self): def test_loading_profiles_ensemble_members(self): """ - This is the most basic check where we have a simple network and check for the basic physics. - This simple network includes two source, pipes, nodes, and 3 demands. - - Checks; - - Demand matching - - Energy conservation - - Heat to discharge - + This test constructs multiple ensemble members based on an "ensemble_member" CSV file + that describes the probability of the ensemble member and the name and number. + The profiles related to each ensemble member are read from the respective CSV files and + saved in for each member. + The test checks if the profiles read match the profiles from the CVS files and if the + ensemble_member_size is set accordingly. """ import models.unit_cases.case_2a_ensemble.src.run_2a as run_2a from models.unit_cases.case_2a_ensemble.src.run_2a import HeatProblemEnsemble From 22f1ed72d8ce0a50ee0ef0b29bacc7ef14867162 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 13 Feb 2026 12:29:12 +0100 Subject: [PATCH 11/19] clean up --- src/mesido/esdl/esdl_mixin.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index 556a85944..9fc5f3b9d 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -692,8 +692,6 @@ def read(self) -> None: None """ super().read() - # TODO: update for ensembles, eg ensemble_size is now still 1, if io gets updated, - # then the size also updates automatically ensemble_size = 1 self.__ensemble = None if self.csv_ensemble_mode: @@ -720,10 +718,6 @@ def read(self) -> None: ensemble_size=ensemble_size, ensemble=self.__ensemble, ) - for ensemble_member_index in range(ensemble_size): - self.io.set_parameter("GeothermalSource_fafd.Max_heat", 10e5, ensemble_member_index) - if ensemble_member_index == 1: - self.io.set_parameter("GeothermalSource_fafd.Max_heat", 2e5, ensemble_member_index) def write(self) -> None: """ From 8bbede830f4534228a0cbf59d473e682e4bf2d12 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 18 Feb 2026 11:30:09 +0100 Subject: [PATCH 12/19] added csv_ensemble_mode as a by default False boolean --- src/mesido/esdl/esdl_mixin.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index 9fc5f3b9d..8c413b2ee 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -72,6 +72,8 @@ class ESDLMixin( for example demand profiles. """ + csv_ensemble_mode: bool = False + esdl_run_info_path: Path = None esdl_pi_validate_timeseries: bool = False From 64d1e82544762d4053dda0ccbe68c881ae06f310 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 18 Feb 2026 12:29:23 +0100 Subject: [PATCH 13/19] adding ensemble parsing to all load_profiles_from_source --- src/mesido/esdl/profile_parser.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mesido/esdl/profile_parser.py b/src/mesido/esdl/profile_parser.py index 8e82177c6..8b9ac5120 100644 --- a/src/mesido/esdl/profile_parser.py +++ b/src/mesido/esdl/profile_parser.py @@ -183,6 +183,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble: None|np.ndarray, ) -> None: """ This function must be implemented by the child. It must load the available @@ -242,6 +243,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble: None | np.ndarray, ) -> None: profiles: Dict[str, np.ndarray] = dict() logger.info("Reading profiles from InfluxDB") @@ -615,7 +617,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, - ensemble, + ensemble: None|np.ndarray, ) -> None: if self._file_path.suffix == ".xml": logger.warning( From 2f80b347fa2888f80f36089f8b0f69819d1a581d Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 18 Feb 2026 12:30:00 +0100 Subject: [PATCH 14/19] style --- src/mesido/esdl/profile_parser.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesido/esdl/profile_parser.py b/src/mesido/esdl/profile_parser.py index 8b9ac5120..116131b7a 100644 --- a/src/mesido/esdl/profile_parser.py +++ b/src/mesido/esdl/profile_parser.py @@ -183,7 +183,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, - ensemble: None|np.ndarray, + ensemble: None | np.ndarray, ) -> None: """ This function must be implemented by the child. It must load the available @@ -617,7 +617,7 @@ def _load_profiles_from_source( esdl_asset_id_to_name_map: Dict[str, str], carrier_properties: Dict[str, Dict], ensemble_size: int, - ensemble: None|np.ndarray, + ensemble: None | np.ndarray, ) -> None: if self._file_path.suffix == ".xml": logger.warning( From bfaa0518ca2069c92f15b6b7d2698232becc6b4a Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 20 Feb 2026 09:33:01 +0100 Subject: [PATCH 15/19] update pandapipes version --- tox.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index 2df141d2d..b72ff5e2f 100644 --- a/tox.ini +++ b/tox.ini @@ -15,7 +15,7 @@ deps = pytest-ordering pytest-timeout numpy - pandapipes==0.10.0 + pandapipes==0.11.0 pandapower==2.14.6 deepdiff==7.0.1 extras = all From 4a73141efce59925a31a1c7fee186a764cd0dafd Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 20 Feb 2026 10:03:08 +0100 Subject: [PATCH 16/19] test updating pandapower and pandapipes requirements --- tox.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tox.ini b/tox.ini index b72ff5e2f..543bf52ba 100644 --- a/tox.ini +++ b/tox.ini @@ -15,8 +15,8 @@ deps = pytest-ordering pytest-timeout numpy - pandapipes==0.11.0 - pandapower==2.14.6 + pandapipes>=0.11.0 + pandapower>=2.14.6 deepdiff==7.0.1 extras = all From 1fe8979dfd38ccfcd2aaf418f6dbad7437e3c4e2 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 25 Feb 2026 09:45:38 +0100 Subject: [PATCH 17/19] code clean up after review --- .../case_2a_ensemble/input/ensemble.csv | 4 ++-- .../timeseries_import.csv | 0 .../timeseries_import.csv | 0 .../unit_cases/case_2a_ensemble/src/run_2a.py | 19 ------------------- 4 files changed, 2 insertions(+), 21 deletions(-) rename tests/models/unit_cases/case_2a_ensemble/input/{forecast1 => forecast_high_prob}/timeseries_import.csv (100%) rename tests/models/unit_cases/case_2a_ensemble/input/{forecast2 => forecast_low_prob}/timeseries_import.csv (100%) diff --git a/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv b/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv index 11521b7ac..c26cb83ee 100644 --- a/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv +++ b/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv @@ -1,3 +1,3 @@ Unnamed: 0,name,probability -0,forecast1,0.7 -1,forecast2,0.3 +0,forecast_high_prob,0.7 +1,forecast_low_prob,0.3 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast_high_prob/timeseries_import.csv similarity index 100% rename from tests/models/unit_cases/case_2a_ensemble/input/forecast1/timeseries_import.csv rename to tests/models/unit_cases/case_2a_ensemble/input/forecast_high_prob/timeseries_import.csv diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast_low_prob/timeseries_import.csv similarity index 100% rename from tests/models/unit_cases/case_2a_ensemble/input/forecast2/timeseries_import.csv rename to tests/models/unit_cases/case_2a_ensemble/input/forecast_low_prob/timeseries_import.csv diff --git a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py index 0b75b08e4..ed69bd82d 100644 --- a/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -27,9 +27,6 @@ def __init__(self, state, target, index): self.index = index self.targets = target - # self.target_min = target[0].values[index] - # self.target_max = target[0].values[index] - # self.function_range = (0.0, 2.0 * max(target[0].values)) self.function_nominal = np.median(target[0].values[index]) def function(self, optimization_problem, ensemble_member): @@ -41,8 +38,6 @@ def function(self, optimization_problem, ensemble_member): vector_state = ( nom * optimization_problem.state_vector(self.state, ensemble_member)[self.index] ) - # self.target_min = self.targets[ensemble_member].values[self.index] - # self.target_max = self.targets[ensemble_member].values[self.index] return self.targets[ensemble_member].values[self.index] - vector_state @@ -84,19 +79,6 @@ def function(self, optimization_problem, ensemble_member): class _GoalsAndOptions: - # def path_goals(self): - # goals = super().path_goals().copy() - # - # for demand in self.energy_system_components["heat_demand"]: - # target = [self.get_timeseries(f"{demand}.target_heat_demand", - # ensemble_member=ensemble_member) for ensemble_member in - # range(self.ensemble_size)] - # state = f"{demand}.Heat_demand" - # target = self.get_timeseries(f"{demand}.target_heat_demand", ensemble_member=0) - # - # goals.append(TargetDemandPathGoal(state, target)) - # - # return goals def goals(self): goals = super().goals().copy() @@ -199,4 +181,3 @@ def path_constraints(self, ensemble_member): profile_reader=ProfileReaderFromFile, input_timeseries_file="timeseries_import.xml", ) - # run_heat_network_optimization(HeatProblem, QTHProblem) From 64b29b1ba179942e7145bceae441dab7312c3c11 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 25 Feb 2026 09:52:01 +0100 Subject: [PATCH 18/19] revert pandapipes package --- tox.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tox.ini b/tox.ini index 543bf52ba..2df141d2d 100644 --- a/tox.ini +++ b/tox.ini @@ -15,8 +15,8 @@ deps = pytest-ordering pytest-timeout numpy - pandapipes>=0.11.0 - pandapower>=2.14.6 + pandapipes==0.10.0 + pandapower==2.14.6 deepdiff==7.0.1 extras = all From fadb859d88288ecd16e3bf9b356edc849d5d5a87 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 25 Feb 2026 10:12:16 +0100 Subject: [PATCH 19/19] changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e2ec2f9b..bae68ff7f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,7 @@ # [Unreleased-main] - 2026-02-25 ## Added -- xxx +- Parsing of ensemble profiles when using input profiles from csv. ## Changed - Speed-up timeseries check in from InfluxDB