diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py index 09add0e7d..687e66f64 100644 --- a/examples/municipality/src/run_municipality.py +++ b/examples/municipality/src/run_municipality.py @@ -1,8 +1,8 @@ from pathlib import Path from mesido.esdl.esdl_parser import ESDLFileParser -from mesido.workflows import EndScenarioSizingStagedHIGHS, run_end_scenario_sizing - +from mesido.esdl.profile_parser import ProfileReaderFromFile +from mesido.workflows import EndScenarioSizingHeadLossStaged, run_end_scenario_sizing if __name__ == "__main__": import time @@ -11,10 +11,12 @@ base_folder = Path(__file__).resolve().parent.parent solution = run_end_scenario_sizing( - EndScenarioSizingStagedHIGHS, + EndScenarioSizingHeadLossStaged, base_folder=base_folder, esdl_file_name="GROW_withATES_Prod_install.esdl", esdl_parser=ESDLFileParser, + profile_reader=ProfileReaderFromFile, + input_timeseries_file="demand_profiles.csv", ) print("Execution time: " + time.strftime("%M:%S", time.gmtime(time.time() - start_time))) diff --git a/src/mesido/asset_sizing_mixin.py b/src/mesido/asset_sizing_mixin.py index 1d9d8848a..d8ae0db03 100644 --- a/src/mesido/asset_sizing_mixin.py +++ b/src/mesido/asset_sizing_mixin.py @@ -1367,7 +1367,7 @@ def __pipe_topology_constraints(self, ensemble_member): ) ) - # These are the constraints to order the milp loss of the pipe classes. + # These are the constraints to order the heat loss of the pipe classes. if not self.energy_system_options()["neglect_pipe_heat_losses"]: for ( pipe, diff --git a/src/mesido/head_loss_class.py b/src/mesido/head_loss_class.py index ae2b64365..fdcf0882d 100644 --- a/src/mesido/head_loss_class.py +++ b/src/mesido/head_loss_class.py @@ -719,6 +719,11 @@ def _hn_pipe_head_loss( ): n_linear_lines = network_settings["n_linearization_lines"] n_timesteps = len(optimization_problem.times()) + # TODO: check if we want to add a simplification for the linearisation of headloss when + # the pipe is relatively short compared to other pipes, thus with relatively small + # effects. + # if length<100: + # n_linear_lines = 1 a, b = darcy_weisbach.get_linear_pipe_dh_vs_q_fit( diameter, diff --git a/src/mesido/workflows/__init__.py b/src/mesido/workflows/__init__.py index 9790c989f..86344cabe 100644 --- a/src/mesido/workflows/__init__.py +++ b/src/mesido/workflows/__init__.py @@ -3,6 +3,7 @@ EndScenarioSizingDiscountedHIGHS, EndScenarioSizingDiscountedStagedHIGHS, EndScenarioSizingHIGHS, + EndScenarioSizingHeadLossStaged, EndScenarioSizingStaged, EndScenarioSizingStagedHIGHS, SolverGurobi, @@ -23,6 +24,7 @@ "EndScenarioSizingDiscountedHIGHS", "EndScenarioSizingDiscountedStagedHIGHS", "EndScenarioSizingHIGHS", + "EndScenarioSizingHeadLossStaged", "EndScenarioSizingStaged", "EndScenarioSizingStagedHIGHS", "SolverGurobi", diff --git a/src/mesido/workflows/grow_workflow.py b/src/mesido/workflows/grow_workflow.py index 2842bff27..0fa8b20b0 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -271,7 +271,7 @@ def solver_options(self): if options["solver"] == "highs": highs_options = options["highs"] if self.__priority == 1: - highs_options["time_limit"] = 100 + highs_options["time_limit"] = 1000 else: highs_options["time_limit"] = 100000 return options @@ -422,6 +422,13 @@ def __init__(self, *args, **kwargs): HeadLossOption.LINEARIZED_N_LINES_WEAK_INEQUALITY ) self.heat_network_settings["minimize_head_losses"] = True + # TODO: check with we want to set a default minimum and maximum pressure for this workflow + # self.heat_network_settings["pipe_minimum_pressure"] = 4 + # self.heat_network_settings["pipe_maximum_pressure"] = 25 + + # TODO: if headloss and thus also hydraulic power is included we might want to use some default + # costs for electricity profile as otherwise the minimisation of hydraulic power as a separate + # goal might be cmputationally intensive. class EndScenarioSizingHeadLossDiscounted(EndScenarioSizingHeadLoss, EndScenarioSizingDiscounted): @@ -448,9 +455,11 @@ def __init__( *args, **kwargs, ): + # _stage needs to be set before the super to ensure that the correct stage is passed on + # to the other classes. + self._stage = stage super().__init__(*args, **kwargs) - self._stage = stage self.__boolean_bounds = boolean_bounds if self._stage == 1: @@ -458,7 +467,7 @@ def __init__( self.heat_network_settings["head_loss_option"] = HeadLossOption.NO_HEADLOSS self.heat_network_settings["minimize_head_losses"] = False - if self._stage == 2 and priorities_output: + if self._stage != 1 and priorities_output: self._priorities_output = priorities_output def energy_system_options(self): @@ -472,7 +481,7 @@ def energy_system_options(self): def bounds(self): bounds = super().bounds() - if self._stage == 2: + if self._stage != 1: bounds.update(self.__boolean_bounds) return bounds @@ -558,6 +567,88 @@ def run_end_scenario_sizing_no_heat_losses( return solution +def create_staged_bounds(solution, results, parameters, bounds, new_stage): + """ + The additional bounds required for the staged approach are added here. + + new_stage: passes the new stage number that will be started after these bounds are added. This + parameter can be used to select additional/different bounds for different stages + + """ + boolean_bounds = {} + # We give bounds for stage 2 by allowing one DN sizes larger than what was found in the + # stage 1 optimization. + pc_map = solution.get_pipe_class_map() + for pipe_classes in pc_map.values(): + v_prev = 0.0 + first_pipe_class = True + for var_name in pipe_classes.values(): + v = results[var_name][0] + if new_stage == 2: + if first_pipe_class and abs(v) == 1.0: + boolean_bounds[var_name] = (abs(v), abs(v)) + elif abs(v) == 1.0: + boolean_bounds[var_name] = (0.0, abs(v)) + elif v_prev == 1.0: + boolean_bounds[var_name] = (0.0, 1.0) + else: + boolean_bounds[var_name] = (abs(v), abs(v)) + else: + boolean_bounds[var_name] = (abs(v), abs(v)) + v_prev = v + first_pipe_class = False + + for asset in [ + *solution.energy_system_components.get("heat_source", []), + *solution.energy_system_components.get("heat_buffer", []), + ]: + var_name = f"{asset}_aggregation_count" + lb = results[var_name][0] + ub = solution.bounds()[var_name][1] + if round(lb) >= 1: + boolean_bounds[var_name] = (lb, ub) + + t = solution.times() + from rtctools.optimization.timeseries import Timeseries + + for p in solution.energy_system_components.get("heat_pipe", []): + if p in solution.hot_pipes and parameters[f"{p}.area"] > 0.0: + lb = [] + ub = [] + bounds_pipe = bounds[f"{p}__flow_direct_var"] + for i in range(len(t)): + r = results[f"{p}__flow_direct_var"][i] + # bound to roughly represent 40km of heat losses in pipes as other new producers + # might also be required and change the direction + lb.append( + r + if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 + else ( + bounds_pipe[0] + if not isinstance(bounds_pipe[0], Timeseries) + else bounds_pipe[0].values[i] + ) + ) + ub.append( + r + if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 + else ( + bounds_pipe[1] + if not isinstance(bounds_pipe[1], Timeseries) + else bounds_pipe[1].values[i] + ) + ) + + boolean_bounds[f"{p}__flow_direct_var"] = (Timeseries(t, lb), Timeseries(t, ub)) + try: + r = results[f"{p}__is_disconnected"] + boolean_bounds[f"{p}__is_disconnected"] = (Timeseries(t, r), Timeseries(t, r)) + except KeyError: + pass + + return boolean_bounds + + def run_end_scenario_sizing( end_scenario_problem_class, staged_pipe_optimization=True, @@ -598,63 +689,7 @@ def run_end_scenario_sizing( parameters = solution.parameters(0) bounds = solution.bounds() - # We give bounds for stage 2 by allowing one DN sizes larger than what was found in the - # stage 1 optimization. - pc_map = solution.get_pipe_class_map() - for pipe_classes in pc_map.values(): - v_prev = 0.0 - first_pipe_class = True - for var_name in pipe_classes.values(): - v = results[var_name][0] - if first_pipe_class and abs(v) == 1.0: - boolean_bounds[var_name] = (abs(v), abs(v)) - elif abs(v) == 1.0: - boolean_bounds[var_name] = (0.0, abs(v)) - elif v_prev == 1.0: - boolean_bounds[var_name] = (0.0, 1.0) - else: - boolean_bounds[var_name] = (abs(v), abs(v)) - v_prev = v - first_pipe_class = False - - for asset in [ - *solution.energy_system_components.get("heat_source", []), - *solution.energy_system_components.get("heat_buffer", []), - ]: - var_name = f"{asset}_aggregation_count" - lb = results[var_name][0] - ub = solution.bounds()[var_name][1] - if round(lb) >= 1: - boolean_bounds[var_name] = (lb, ub) - - t = solution.times() - from rtctools.optimization.timeseries import Timeseries - - for p in solution.energy_system_components.get("heat_pipe", []): - if p in solution.hot_pipes and parameters[f"{p}.area"] > 0.0: - lb = [] - ub = [] - bounds_pipe = bounds[f"{p}__flow_direct_var"] - for i in range(len(t)): - r = results[f"{p}__flow_direct_var"][i] - # bound to roughly represent 4km of milp losses in pipes - lb.append( - r - if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-2 - else bounds_pipe[0] - ) - ub.append( - r - if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-2 - else bounds_pipe[1] - ) - - boolean_bounds[f"{p}__flow_direct_var"] = (Timeseries(t, lb), Timeseries(t, ub)) - try: - r = results[f"{p}__is_disconnected"] - boolean_bounds[f"{p}__is_disconnected"] = (Timeseries(t, r), Timeseries(t, r)) - except KeyError: - pass + boolean_bounds = create_staged_bounds(solution, results, parameters, bounds, new_stage=2) priorities_output = solution._priorities_output solution = run_optimization_problem( @@ -665,6 +700,21 @@ def run_end_scenario_sizing( **kwargs, ) + # if issubclass(end_scenario_problem_class, EndScenarioSizingHeadLoss): + # results = solution.extract_results() + # parameters = solution.parameters(0) + # bounds = solution.bounds() + # + # boolean_bounds, priorities_output = create_staged_bounds(solution, results, parameters, + # bounds, new_stage=3) + # solution = run_optimization_problem( + # end_scenario_problem_class, + # stage=3, + # boolean_bounds=boolean_bounds, + # priorities_output=priorities_output, + # **kwargs, + # ) + print("Execution time: " + time.strftime("%M:%S", time.gmtime(time.time() - start_time))) return solution