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 diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py index d8d4ced2b..52cc77ad2 100644 --- a/src/mesido/esdl/esdl_mixin.py +++ b/src/mesido/esdl/esdl_mixin.py @@ -2,6 +2,7 @@ import copy import dataclasses import logging +import os import sys from datetime import timedelta from pathlib import Path @@ -71,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 @@ -691,6 +694,19 @@ def read(self) -> None: None """ super().read() + ensemble_size = 1 + self.__ensemble = None + if self.csv_ensemble_mode: + self.__ensemble = np.genfromtxt( + os.path.join(self._input_folder, "ensemble.csv"), + delimiter=",", + deletechars="", + dtype=None, + names=True, + encoding=None, + ) + ensemble_size = len(self.__ensemble) + energy_system_components = self.energy_system_components esdl_carriers = self.esdl_carriers self.hot_cold_pipe_relations() @@ -701,7 +717,8 @@ 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, + ensemble=self.__ensemble, ) def write(self) -> None: @@ -826,3 +843,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 23d939b52..0c7385089 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: @@ -181,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 @@ -240,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") @@ -616,6 +620,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: if self._file_path.suffix == ".xml": logger.warning( @@ -630,6 +635,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( @@ -641,25 +647,27 @@ def _load_csv( energy_system_components: Dict[str, Set[str]], carrier_properties: Dict[str, Dict], ensemble_size: int, + ensemble, ) -> None: - 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() @@ -668,48 +676,58 @@ def _load_csv( try: timeseries_import_times = [ datetime.datetime.strptime( - entry.replace("Z", ""), "%d-%m-%Y %H:%M" + 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: " + f"{data_em[data_em[carrier_name].isnull()]}" ) except KeyError: pass else: - self._profiles[ensemble_member][ - 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/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..c26cb83ee --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/ensemble.csv @@ -0,0 +1,3 @@ +Unnamed: 0,name,probability +0,forecast_high_prob,0.7 +1,forecast_low_prob,0.3 diff --git a/tests/models/unit_cases/case_2a_ensemble/input/forecast_high_prob/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast_high_prob/timeseries_import.csv new file mode 100644 index 000000000..1e9ee86c7 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/forecast_high_prob/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/forecast_low_prob/timeseries_import.csv b/tests/models/unit_cases/case_2a_ensemble/input/forecast_low_prob/timeseries_import.csv new file mode 100644 index 000000000..74e7a9e13 --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/input/forecast_low_prob/timeseries_import.csv @@ -0,0 +1,4 @@ +DateTime,HeatingDemand_7484,HeatingDemand_c6c8,HeatingDemand_6f99 +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 new file mode 100644 index 000000000..1410c2520 --- /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..ed69bd82d --- /dev/null +++ b/tests/models/unit_cases/case_2a_ensemble/src/run_2a.py @@ -0,0 +1,183 @@ +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 + +import numpy as np + +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.linearized_order_goal_programming_mixin import ( + LinearizedOrderGoalProgrammingMixin, +) +from rtctools.util import run_optimization_problem + + +class TargetDemandGoal(Goal): + priority = 1 + + order = 1 + + def __init__(self, state, target, index): + self.state = state + self.index = index + + self.targets = target + 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] + ) + + return self.targets[ensemble_member].values[self.index] - vector_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): + """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) + + +class MinimizeSourcesSizeGoal(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.extra_variable(f"{self.source}__max_size", ensemble_member) + + +class _GoalsAndOptions: + + def goals(self): + goals = super().goals().copy() + + 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 + + +class HeatProblem( + _GoalsAndOptions, + PhysicsMixin, + LinearizedOrderGoalProgrammingMixin, + GoalProgrammingMixin, + ESDLMixin, + CollocatedIntegratedOptimizationProblem, +): + + 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, + ControlTreeMixin, +): + csv_ensemble_mode = True + + 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): + 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 __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) + 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 + + +if __name__ == "__main__": + run_optimization_problem( + HeatProblem, + esdl_file_name="2a.esdl", + esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="timeseries_import.xml", + ) diff --git a/tests/test_profile_parsing.py b/tests/test_profile_parsing.py index fd7439185..0184657c6 100644 --- a/tests/test_profile_parsing.py +++ b/tests/test_profile_parsing.py @@ -315,6 +315,53 @@ def test_loading_profile_from_esdl(self): atol=1e-2, ) + def test_loading_profiles_ensemble_members(self): + """ + 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 + + 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) + 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 + # 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() diff --git a/tests/utils_tests.py b/tests/utils_tests.py index b42118929..3946b4b8b 100644 --- a/tests/utils_tests.py +++ b/tests/utils_tests.py @@ -21,21 +21,25 @@ 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: 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 +48,7 @@ def demand_matching_test(solution, results, atol=1.0e-3, rtol=1.0e-6): 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 +59,7 @@ def demand_matching_test(solution, results, atol=1.0e-3, rtol=1.0e-6): 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 )