Speedup staged headloss

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)))

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,

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,

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",

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,17 +455,19 @@ 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:
             self.heat_network_settings["minimum_velocity"] = 0.0
             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