From de1bf22b6796f7f8934a4a1497a483fa24a3665c Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 8 May 2024 15:51:01 +0200 Subject: [PATCH 1/5] added third stage for headloss minimisation --- examples/municipality/src/run_municipality.py | 8 +- src/mesido/asset_sizing_mixin.py | 2 +- src/mesido/head_loss_class.py | 2 + src/mesido/workflows/grow_workflow.py | 231 +++++++++++++----- 4 files changed, 176 insertions(+), 67 deletions(-) diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py index 09add0e7d..d7a42a648 100644 --- a/examples/municipality/src/run_municipality.py +++ b/examples/municipality/src/run_municipality.py @@ -1,8 +1,10 @@ from pathlib import Path from mesido.esdl.esdl_parser import ESDLFileParser +from mesido.esdl.profile_parser import ProfileReaderFromFile from mesido.workflows import EndScenarioSizingStagedHIGHS, run_end_scenario_sizing - +from mesido.workflows.grow_workflow import EndScenarioSizingHeadLossStagedGurobi, \ + EndScenarioSizingStagedGurobi, EndScenarioSizingHeadLossStaged if __name__ == "__main__": import time @@ -11,10 +13,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..613c77804 100644 --- a/src/mesido/head_loss_class.py +++ b/src/mesido/head_loss_class.py @@ -719,6 +719,8 @@ def _hn_pipe_head_loss( ): n_linear_lines = network_settings["n_linearization_lines"] n_timesteps = len(optimization_problem.times()) + 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/grow_workflow.py b/src/mesido/workflows/grow_workflow.py index 2842bff27..801bf54d0 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -157,7 +157,7 @@ def __init__(self, *args, **kwargs): self.__heat_demand_bounds = dict() self.__heat_demand_nominal = dict() - self._save_json = False + self._save_json = True def parameters(self, ensemble_member): parameters = super().parameters(ensemble_member) @@ -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 @@ -421,7 +421,11 @@ def __init__(self, *args, **kwargs): self.heat_network_settings["head_loss_option"] = ( HeadLossOption.LINEARIZED_N_LINES_WEAK_INEQUALITY ) - self.heat_network_settings["minimize_head_losses"] = True + if self._stage == 3: + self.heat_network_settings["minimize_head_losses"] = True + self.heat_network_settings["pipe_minimum_pressure"] = 4 + self.heat_network_settings["pipe_maximum_pressure"] = 25 + self.heat_network_settings["n_linearization_lines"] = 1 class EndScenarioSizingHeadLossDiscounted(EndScenarioSizingHeadLoss, EndScenarioSizingDiscounted): @@ -448,9 +452,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 +464,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 +478,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 +564,74 @@ def run_end_scenario_sizing_no_heat_losses( return solution +def create_staged_bounds(solution, results, parameters, bounds, new_stage): + 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 + + print(f"number of boolean bounds for stage 2 before flow directions: {len(boolean_bounds)}") + + 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 heat losses in pipes + 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 + print(f"number of boolean bounds for stage 2 after flow directions: {len(boolean_bounds)}") + priorities_output = solution._priorities_output + return boolean_bounds, priorities_output + def run_end_scenario_sizing( end_scenario_problem_class, staged_pipe_optimization=True, @@ -598,64 +672,77 @@ 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 - priorities_output = solution._priorities_output + boolean_bounds, priorities_output = create_staged_bounds(solution, results, parameters, bounds, new_stage=2) + + # # 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 + # v_prev_2 = 0.0 + # v_prev_3 = 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) + # elif v_prev_2 == 1.0: + # boolean_bounds[var_name] = (0.0, 1.0) + # # elif v_prev_3 == 1.0: + # # boolean_bounds[var_name] = (0.0, 1.0) + # else: + # boolean_bounds[var_name] = (abs(v), abs(v)) + # v_prev = v + # v_prev_2 = v_prev + # v_prev_3 = v_prev_2 + # 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 + # + # print(f"number of boolean bounds for stage 2 before flow directions: {len(boolean_bounds)}") + # + # 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 heat losses in pipes + # lb.append( + # r + # if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 + # else bounds_pipe[0] + # ) + # ub.append( + # r + # if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 + # 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 + # print(f"number of boolean bounds for stage 2 after flow directions: {len(boolean_bounds)}") + # priorities_output = solution._priorities_output solution = run_optimization_problem( end_scenario_problem_class, @@ -665,6 +752,22 @@ 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 From ed2ddeb5c83e2a2b145d68163118d982bee6b95e Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Fri, 10 May 2024 18:08:05 +0200 Subject: [PATCH 2/5] wip --- examples/municipality/src/run_municipality.py | 2 +- src/mesido/financial_mixin.py | 2 +- src/mesido/workflows/grow_workflow.py | 32 +++++++++---------- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py index d7a42a648..f737eed26 100644 --- a/examples/municipality/src/run_municipality.py +++ b/examples/municipality/src/run_municipality.py @@ -13,7 +13,7 @@ base_folder = Path(__file__).resolve().parent.parent solution = run_end_scenario_sizing( - EndScenarioSizingHeadLossStaged, + EndScenarioSizingHeadLossStagedGurobi, base_folder=base_folder, esdl_file_name="GROW_withATES_Prod_install.esdl", esdl_parser=ESDLFileParser, diff --git a/src/mesido/financial_mixin.py b/src/mesido/financial_mixin.py index 5e01a9dc7..60a7ed39a 100644 --- a/src/mesido/financial_mixin.py +++ b/src/mesido/financial_mixin.py @@ -825,7 +825,7 @@ def __variable_operational_cost_constraints(self, ensemble_member): f"{list(self.get_electricity_carriers().values())[0]['name']}.price_profile" ) else: - price_profile = Timeseries(self.times(), np.zeros(len(self.times()))) + price_profile = Timeseries(self.times(), 0.2e-3*np.ones(len(self.times()))) timesteps = np.diff(self.times()) / 3600.0 diff --git a/src/mesido/workflows/grow_workflow.py b/src/mesido/workflows/grow_workflow.py index 801bf54d0..1884f09c4 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -421,8 +421,8 @@ def __init__(self, *args, **kwargs): self.heat_network_settings["head_loss_option"] = ( HeadLossOption.LINEARIZED_N_LINES_WEAK_INEQUALITY ) - if self._stage == 3: - self.heat_network_settings["minimize_head_losses"] = True + # if self._stage == 3: + self.heat_network_settings["minimize_head_losses"] = True self.heat_network_settings["pipe_minimum_pressure"] = 4 self.heat_network_settings["pipe_maximum_pressure"] = 25 self.heat_network_settings["n_linearization_lines"] = 1 @@ -752,20 +752,20 @@ 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, - ) + # 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))) From 9ed415cab1f809ae23c462b9fba4a45bf0273ccc Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 22 May 2024 15:13:19 +0200 Subject: [PATCH 3/5] clean up remove commented code --- examples/municipality/src/run_municipality.py | 2 +- src/mesido/head_loss_class.py | 7 +- src/mesido/workflows/grow_workflow.py | 76 +------------------ 3 files changed, 9 insertions(+), 76 deletions(-) diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py index f737eed26..d7a42a648 100644 --- a/examples/municipality/src/run_municipality.py +++ b/examples/municipality/src/run_municipality.py @@ -13,7 +13,7 @@ base_folder = Path(__file__).resolve().parent.parent solution = run_end_scenario_sizing( - EndScenarioSizingHeadLossStagedGurobi, + EndScenarioSizingHeadLossStaged, base_folder=base_folder, esdl_file_name="GROW_withATES_Prod_install.esdl", esdl_parser=ESDLFileParser, diff --git a/src/mesido/head_loss_class.py b/src/mesido/head_loss_class.py index 613c77804..fc51819e6 100644 --- a/src/mesido/head_loss_class.py +++ b/src/mesido/head_loss_class.py @@ -719,8 +719,11 @@ def _hn_pipe_head_loss( ): n_linear_lines = network_settings["n_linearization_lines"] n_timesteps = len(optimization_problem.times()) - if length<100: - n_linear_lines = 1 + #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/grow_workflow.py b/src/mesido/workflows/grow_workflow.py index 1884f09c4..e6c46d8c4 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -601,8 +601,6 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): t = solution.times() from rtctools.optimization.timeseries import Timeseries - print(f"number of boolean bounds for stage 2 before flow directions: {len(boolean_bounds)}") - for p in solution.energy_system_components.get("heat_pipe", []): if p in solution.hot_pipes and parameters[f"{p}.area"] > 0.0: lb = [] @@ -610,7 +608,8 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): 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 heat losses in pipes + # 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 @@ -628,7 +627,7 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): boolean_bounds[f"{p}__is_disconnected"] = (Timeseries(t, r), Timeseries(t, r)) except KeyError: pass - print(f"number of boolean bounds for stage 2 after flow directions: {len(boolean_bounds)}") + priorities_output = solution._priorities_output return boolean_bounds, priorities_output @@ -674,75 +673,6 @@ def run_end_scenario_sizing( boolean_bounds, priorities_output = create_staged_bounds(solution, results, parameters, bounds, new_stage=2) - # # 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 - # v_prev_2 = 0.0 - # v_prev_3 = 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) - # elif v_prev_2 == 1.0: - # boolean_bounds[var_name] = (0.0, 1.0) - # # elif v_prev_3 == 1.0: - # # boolean_bounds[var_name] = (0.0, 1.0) - # else: - # boolean_bounds[var_name] = (abs(v), abs(v)) - # v_prev = v - # v_prev_2 = v_prev - # v_prev_3 = v_prev_2 - # 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 - # - # print(f"number of boolean bounds for stage 2 before flow directions: {len(boolean_bounds)}") - # - # 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 heat losses in pipes - # lb.append( - # r - # if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 - # else bounds_pipe[0] - # ) - # ub.append( - # r - # if abs(results[f"{p}.Q"][i] / parameters[f"{p}.area"]) > 2.5e-1 - # 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 - # print(f"number of boolean bounds for stage 2 after flow directions: {len(boolean_bounds)}") - # priorities_output = solution._priorities_output solution = run_optimization_problem( end_scenario_problem_class, From 1e2d7076c17a2e9a851e53b7edf2a2a931360b84 Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 22 May 2024 15:26:12 +0200 Subject: [PATCH 4/5] docstring and clean up --- src/mesido/workflows/grow_workflow.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/mesido/workflows/grow_workflow.py b/src/mesido/workflows/grow_workflow.py index e6c46d8c4..35dce0ef0 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -157,7 +157,7 @@ def __init__(self, *args, **kwargs): self.__heat_demand_bounds = dict() self.__heat_demand_nominal = dict() - self._save_json = True + self._save_json = False def parameters(self, ensemble_member): parameters = super().parameters(ensemble_member) @@ -421,11 +421,10 @@ def __init__(self, *args, **kwargs): self.heat_network_settings["head_loss_option"] = ( HeadLossOption.LINEARIZED_N_LINES_WEAK_INEQUALITY ) - # if self._stage == 3: self.heat_network_settings["minimize_head_losses"] = True - self.heat_network_settings["pipe_minimum_pressure"] = 4 - self.heat_network_settings["pipe_maximum_pressure"] = 25 - self.heat_network_settings["n_linearization_lines"] = 1 + #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 class EndScenarioSizingHeadLossDiscounted(EndScenarioSizingHeadLoss, EndScenarioSizingDiscounted): @@ -565,6 +564,13 @@ def run_end_scenario_sizing_no_heat_losses( 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. @@ -628,8 +634,7 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): except KeyError: pass - priorities_output = solution._priorities_output - return boolean_bounds, priorities_output + return boolean_bounds def run_end_scenario_sizing( end_scenario_problem_class, @@ -671,8 +676,8 @@ def run_end_scenario_sizing( parameters = solution.parameters(0) bounds = solution.bounds() - boolean_bounds, priorities_output = create_staged_bounds(solution, results, parameters, bounds, new_stage=2) - + boolean_bounds = create_staged_bounds(solution, results, parameters, bounds, new_stage=2) + priorities_output = solution._priorities_output solution = run_optimization_problem( end_scenario_problem_class, From 44a555c304e5ea09b7049501b8c2427fd699482d Mon Sep 17 00:00:00 2001 From: Femke Janssen Date: Wed, 22 May 2024 15:36:26 +0200 Subject: [PATCH 5/5] style --- examples/municipality/src/run_municipality.py | 4 +--- src/mesido/financial_mixin.py | 2 +- src/mesido/head_loss_class.py | 2 +- src/mesido/workflows/__init__.py | 2 ++ src/mesido/workflows/grow_workflow.py | 20 +++++++++++++++---- 5 files changed, 21 insertions(+), 9 deletions(-) diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py index d7a42a648..687e66f64 100644 --- a/examples/municipality/src/run_municipality.py +++ b/examples/municipality/src/run_municipality.py @@ -2,9 +2,7 @@ from mesido.esdl.esdl_parser import ESDLFileParser from mesido.esdl.profile_parser import ProfileReaderFromFile -from mesido.workflows import EndScenarioSizingStagedHIGHS, run_end_scenario_sizing -from mesido.workflows.grow_workflow import EndScenarioSizingHeadLossStagedGurobi, \ - EndScenarioSizingStagedGurobi, EndScenarioSizingHeadLossStaged +from mesido.workflows import EndScenarioSizingHeadLossStaged, run_end_scenario_sizing if __name__ == "__main__": import time diff --git a/src/mesido/financial_mixin.py b/src/mesido/financial_mixin.py index 60a7ed39a..5e01a9dc7 100644 --- a/src/mesido/financial_mixin.py +++ b/src/mesido/financial_mixin.py @@ -825,7 +825,7 @@ def __variable_operational_cost_constraints(self, ensemble_member): f"{list(self.get_electricity_carriers().values())[0]['name']}.price_profile" ) else: - price_profile = Timeseries(self.times(), 0.2e-3*np.ones(len(self.times()))) + price_profile = Timeseries(self.times(), np.zeros(len(self.times()))) timesteps = np.diff(self.times()) / 3600.0 diff --git a/src/mesido/head_loss_class.py b/src/mesido/head_loss_class.py index fc51819e6..fdcf0882d 100644 --- a/src/mesido/head_loss_class.py +++ b/src/mesido/head_loss_class.py @@ -719,7 +719,7 @@ 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 + # 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: 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 35dce0ef0..0fa8b20b0 100644 --- a/src/mesido/workflows/grow_workflow.py +++ b/src/mesido/workflows/grow_workflow.py @@ -422,10 +422,14 @@ 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 + # 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): pass @@ -619,12 +623,20 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): 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] + 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] + 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)) @@ -636,6 +648,7 @@ def create_staged_bounds(solution, results, parameters, bounds, new_stage): return boolean_bounds + def run_end_scenario_sizing( end_scenario_problem_class, staged_pipe_optimization=True, @@ -702,7 +715,6 @@ def run_end_scenario_sizing( # **kwargs, # ) - print("Execution time: " + time.strftime("%M:%S", time.gmtime(time.time() - start_time))) return solution