Discrete variables

examples/municipality/src/run_municipality.py

@@ -3,7 +3,6 @@
 from mesido.esdl.esdl_parser import ESDLFileParser
 from mesido.workflows import EndScenarioSizingStaged, run_end_scenario_sizing
 
-
 if __name__ == "__main__":
     import time
 

examples/pipe_diameter_sizing/src/example.py

@@ -122,6 +122,7 @@ def solver_options(self):
         self._qpsol = CachingQPSol()
         options["casadi_solver"] = self._qpsol
         options["solver"] = "highs"
+
         return options
 
 

src/mesido/component_type_mixin.py

@@ -102,6 +102,9 @@ def pre(self):
                         other_pipe = pipes_map[other_pipe_port]
                         if f"{other_pipe}.Q" not in alias_relation.canonical_variables:
                             alias_relation.add(f"{p}.Q", f"{sign_prefix}{other_pipe}.Q")
+                if self.has_related_pipe(p):
+                    cold_pipe = self.hot_to_cold_pipe(p)
+                    alias_relation.add(f"{p}.__flow_direct_var", f"{cold_pipe}.__flow_direct_var")
 
         node_to_node_logical_link_map = {}
 

src/mesido/electricity_physics_mixin.py

@@ -56,16 +56,12 @@ def __init__(self, *args, **kwargs):
 
         # Variable for when in time an asset switched on due to meeting a requirement
         self.__asset_is_switched_on_map = {}
-        self.__asset_is_switched_on_var = {}
-        self.__asset_is_switched_on_bounds = {}
 
         self.__electricity_producer_upper_bounds = {}
 
         self._electricity_cable_topo_cable_class_map = {}
 
         # Boolean path-variable for the charging of storage assets
-        self.__storage_charging_var = {}
-        self.__storage_charging_bounds = {}
         self.__storage_charging_map = {}
 
         self.__set_point_var = {}
@@ -74,11 +70,8 @@ def __init__(self, *args, **kwargs):
 
         # Boolean path-variable for the equality constraint of the electrolyzer
         self.__electrolyzer_is_active_linear_segment_map = {}
-        self.__electrolyzer_is_active_linear_segment_var = {}
-        self.__electrolyzer_is_active_linear_segment_bounds = {}
-        self.__electricity_storage_discharge_var = {}
+
         self.__electricity_storage_discharge_bounds = {}
-        self.__electricity_storage_discharge_nominals = {}
         self.__electricity_storage_discharge_map = {}
 
         # Map for setting node nominals in case of logical links.
@@ -123,10 +116,7 @@ def pre(self):
             for asset in [
                 *self.energy_system_components.get("electrolyzer", []),
             ]:
-                var_name = f"{asset}__asset_is_switched_on"
-                self.__asset_is_switched_on_map[asset] = var_name
-                self.__asset_is_switched_on_var[var_name] = ca.MX.sym(var_name)
-                self.__asset_is_switched_on_bounds[var_name] = (0.0, 1.0)
+                self.__asset_is_switched_on_map[asset] = f"{asset}.__asset_is_switched_on"
 
         if options["electrolyzer_efficiency"] == ElectrolyzerOption.LINEARIZED_THREE_LINES_EQUALITY:
             for asset in [
@@ -135,32 +125,23 @@ def pre(self):
                 self.__electrolyzer_is_active_linear_segment_map[asset] = {}
                 n_lines = 3
                 for n_line in range(n_lines):
-                    var_name = f"{asset}__line_{n_line}_active"
-
+                    var_name = f"{asset}.__line_{n_line}_active"
                     self.__electrolyzer_is_active_linear_segment_map[asset][
                         f"line_{n_line}"
                     ] = var_name
-                    self.__electrolyzer_is_active_linear_segment_var[var_name] = ca.MX.sym(var_name)
-                    self.__electrolyzer_is_active_linear_segment_bounds[var_name] = (0.0, 1.0)
 
         for asset in [*self.energy_system_components.get("electricity_storage", [])]:
-            var_name = f"{asset}__is_charging"
+            var_name = f"{asset}.__is_charging"
             self.__storage_charging_map[asset] = var_name
-            self.__storage_charging_var[var_name] = ca.MX.sym(var_name)
-            self.__storage_charging_bounds[var_name] = (0.0, 1.0)
 
             if options["electricity_storage_discharge_variables"]:
                 bound_storage = -self.bounds()[f"{asset}.Effective_power_charging"][0]
                 if isinstance(bound_storage, Timeseries):
                     bound_storage = copy.deepcopy(bound_storage)
                     bound_storage.values[bound_storage.values < 0] = 0.0
-                var_name = f"{asset}__effective_power_discharging"
+                var_name = f"{asset}.__effective_power_discharging"
                 self.__electricity_storage_discharge_map[asset] = var_name
-                self.__electricity_storage_discharge_var[var_name] = ca.MX.sym(var_name)
                 self.__electricity_storage_discharge_bounds[var_name] = (0, bound_storage)
-                self.__electricity_storage_discharge_nominals[var_name] = self.variable_nominal(
-                    f"{asset}.Effective_power_charging"
-                )
 
         for asset in [*self.energy_system_components.get("electricity_source", [])]:
             if isinstance(self.bounds()[f"{asset}.Electricity_source"][1], Timeseries):
@@ -223,11 +204,7 @@ def path_variables(self):
         """
         variables = super().path_variables.copy()
 
-        variables.extend(self.__asset_is_switched_on_var.values())
-        variables.extend(self.__storage_charging_var.values())
         variables.extend(self.__set_point_var.values())
-        variables.extend(self.__electrolyzer_is_active_linear_segment_var.values())
-        variables.extend(self.__electricity_storage_discharge_var.values())
 
         return variables
 
@@ -236,22 +213,14 @@ def variable_is_discrete(self, variable):
         All variables that only can take integer values should be added to this function.
         """
 
-        if variable in self.__electrolyzer_is_active_linear_segment_var:
-            return True
-        if variable in self.__asset_is_switched_on_var:
-            return True
-        if variable in self.__storage_charging_var:
-            return True
-        else:
-            return super().variable_is_discrete(variable)
+        return super().variable_is_discrete(variable)
 
     def variable_nominal(self, variable):
         """
         In this function we add all the nominals for the variables defined/added in the HeatMixin.
         """
-        if variable in self.__electricity_storage_discharge_nominals:
-            return self.__electricity_storage_discharge_nominals[variable]
-        elif variable in self.__bus_variable_nominal:
+
+        if variable in self.__bus_variable_nominal:
             return self.__bus_variable_nominal[variable]
         else:
             return super().variable_nominal(variable)
@@ -263,9 +232,6 @@ def bounds(self):
         """
         bounds = super().bounds()
 
-        bounds.update(self.__electrolyzer_is_active_linear_segment_bounds)
-        bounds.update(self.__asset_is_switched_on_bounds)
-        bounds.update(self.__storage_charging_bounds)
         bounds.update(self.__electricity_producer_upper_bounds)
         bounds.update(self.__set_point_bounds)
         bounds.update(self.__electricity_storage_discharge_bounds)
@@ -568,7 +534,7 @@ def __electricity_storage_path_constraints(self, ensemble_member):
 
             # is_charging is 1 if charging and powerin>0
             big_m = 2 * max(np.abs(self.bounds()[f"{asset}.ElectricityIn.Power"]))
-            is_charging = self.state(f"{asset}__is_charging")
+            is_charging = self.state(f"{asset}.__is_charging")
             constraints.append(((power_in + (1 - is_charging) * big_m) / power_nom, 0.0, np.inf))
             constraints.append(((power_in - is_charging * big_m) / power_nom, -np.inf, 0.0))
 
@@ -647,7 +613,7 @@ def __electricity_storage_discharge_var_path_constraints(
             for storage in self.energy_system_components.get("electricity_storage", []):
                 storage_eff_power_charge_var = self.state(f"{storage}.Effective_power_charging")
                 discharge_var_name = self.__electricity_storage_discharge_map[storage]
-                storage_discharge_var = self.__electricity_storage_discharge_var[discharge_var_name]
+                storage_discharge_var = self.state(discharge_var_name)
                 nominal = self.variable_nominal(discharge_var_name)
 
                 # P_effective_charge represents both charging and discharing based on the sign.

src/mesido/esdl/esdl_heat_model.py

@@ -76,6 +76,7 @@ def __init__(
         cp=4200.0,
         min_fraction_tank_volume=0.05,
         v_max_gas=15.0,
+        energy_system_options=None,
         **kwargs,
     ):
         super().__init__(*args, **kwargs)
@@ -86,6 +87,7 @@ def __init__(
         self.cp = cp
         self.v_max_gas = v_max_gas
         self.min_fraction_tank_volume = min_fraction_tank_volume
+        self.energy_system_options = dict() if not energy_system_options else energy_system_options
         if "primary_port_name_convention" in kwargs.keys():
             self.primary_port_name_convention = kwargs["primary_port_name_convention"]
         if "secondary_port_name_convention" in kwargs.keys():
@@ -1721,6 +1723,9 @@ def convert_electricity_storage(
             min_voltage=v_min,
             max_capacity=max_capacity,
             Stored_electricity=dict(min=0.0, max=max_capacity),
+            discharge_var=self.energy_system_options.get(
+                "electricity_storage_discharge_variables", False
+            ),
             ElectricityIn=dict(
                 V=dict(min=v_min, nominal=v_min),
                 I=dict(min=-i_max, max=i_max, nominal=i_nom),
@@ -2155,6 +2160,8 @@ def equations(x):
             Q_nominal=q_nominal,
             density=density,
             efficiency=eff_max,
+            include_asset_is_switched_on=self.energy_system_options["include_asset_is_switched_on"],
+            electrolyzer_efficiency_option=self.energy_system_options["electrolyzer_efficiency"],
             GasOut=dict(
                 Q=dict(
                     min=0.0,
@@ -2222,6 +2229,7 @@ def convert_gas_tank_storage(self, asset: Asset) -> Tuple[Type[GasTankStorage],
             Q_nominal=q_nominal,
             density=density,
             volume=asset.attributes["workingVolume"],
+            discharge_var=self.energy_system_options.get("gas_storage_discharge_variables", False),
             # Gas_tank_flow=dict(min=-hydrogen_specific_energy*asset.attributes["maxDischargeRate"],
             # max=hydrogen_specific_energy*asset.attributes["maxChargeRate"]),
             # TODO: Fix -> Gas network is currenlty non-limiting, mass flow is decoupled from the

src/mesido/esdl/esdl_mixin.py

@@ -118,9 +118,15 @@ def __init__(self, *args, **kwargs) -> None:
         # Although we work with the names, the FEWS import data uses the component IDs
         self.__timeseries_id_map = {a.id: a.name for a in assets.values()}
         name_to_id_map = {a.name: a.id for a in assets.values()}
-
         if isinstance(self, PhysicsMixin):
-            self.__model = ESDLHeatModel(assets, name_to_id_map, **self.esdl_heat_model_options())
+            self.__model = ESDLHeatModel(
+                assets,
+                name_to_id_map,
+                **dict(
+                    energy_system_options=self.energy_system_options(),
+                    **self.esdl_heat_model_options(),
+                ),
+            )
         else:
             assert isinstance(self, QTHMixin)
 

src/mesido/gas_physics_mixin.py

@@ -1,8 +1,6 @@
 import copy
 import logging
 
-import casadi as ca
-
 from mesido.base_component_type_mixin import BaseComponentTypeMixin
 from mesido.head_loss_class import HeadLossClass, HeadLossOption
 from mesido.network_common import NetworkSettings
@@ -113,7 +111,6 @@ def __init__(self, *args, **kwargs):
         self._gn_pipe_to_head_loss_map = {}
 
         # Boolean path-variable for the direction of the flow, inport to outport is positive flow.
-        self.__gas_flow_direct_var = {}
         self.__gas_flow_direct_bounds = {}
         self._gas_pipe_to_flow_direct_map = {}
 
@@ -124,7 +121,6 @@ def __init__(self, *args, **kwargs):
         # self._gas_pipe_disconnect_map = {}
 
         # Boolean variables for the linear line segment options per pipe.
-        # TDOD: change name to _gas_pipe_...
         self.__gas_pipe_linear_line_segment_var = {}  # value 0/1: line segment - not active/active
         self.__gas_pipe_linear_line_segment_var_bounds = {}
         self._gas_pipe_linear_line_segment_map = {}
@@ -133,13 +129,9 @@ def __init__(self, *args, **kwargs):
 
         self._gas_pipe_topo_pipe_class_map = {}
 
-        # self.__gas_pipe_disconnect_var = {}
-        # self.__gas_pipe_disconnect_var_bounds = {}
         self._gas_pipe_disconnect_map = {}
 
-        self.__gas_storage_discharge_var = {}
         self.__gas_storage_discharge_bounds = {}
-        self.__gas_storage_discharge_nominals = {}
         self.__gas_storage_discharge_map = {}
 
         # Map for setting port variable nominals in the case they were not set during the model
@@ -228,10 +220,10 @@ def _get_min_bound(bound):
                     ] = initialized_vars[10][pipe_linear_line_segment_var_name]
 
             # Integer variables
-            flow_dir_var = f"{pipe_name}__gas_flow_direct_var"
+            flow_dir_var = f"{pipe_name}.__gas_flow_direct_var"
 
             self._gas_pipe_to_flow_direct_map[pipe_name] = flow_dir_var
-            self.__gas_flow_direct_var[flow_dir_var] = ca.MX.sym(flow_dir_var)
+            # self.__gas_flow_direct_var[flow_dir_var] = ca.MX.sym(flow_dir_var)
 
             # Fix the directions that are already implied by the bounds on milp
             # Nonnegative milp implies that flow direction Boolean is equal to one.
@@ -260,17 +252,14 @@ def _get_min_bound(bound):
 
         if options["gas_storage_discharge_variables"]:
             for storage in self.energy_system_components.get("gas_tank_storage", []):
+                # updating bounds
                 bound_storage_q = -self.bounds()[f"{storage}.GasIn.Q"][0]
                 if isinstance(bound_storage_q, Timeseries):
                     bound_storage_q = copy.deepcopy(bound_storage_q)
                     bound_storage_q.values[bound_storage_q.values < 0] = 0.0
-                var_name = f"{storage}__Q_discharge"
+                var_name = f"{storage}.__Q_discharge"
                 self.__gas_storage_discharge_map[storage] = var_name
-                self.__gas_storage_discharge_var[var_name] = ca.MX.sym(var_name)
                 self.__gas_storage_discharge_bounds[var_name] = (0, bound_storage_q)
-                self.__gas_storage_discharge_nominals[var_name] = self.variable_nominal(
-                    f"{storage}.GasIn.Q"
-                )
 
         # Setting the node nominals using the connected assets.
         for node, connected_assets in self.energy_system_topology.gas_nodes.items():
@@ -340,21 +329,15 @@ def path_variables(self):
         """
         variables = super().path_variables.copy()
         variables.extend(self.__gas_pipe_head_loss_var.values())
-        variables.extend(self.__gas_flow_direct_var.values())
-        # variables.extend(self.__gas_pipe_disconnect_var.values())  # still to be implemented
         variables.extend(self.__gas_pipe_linear_line_segment_var.values())
-        variables.extend(self.__gas_storage_discharge_var.values())
 
         return variables
 
     def variable_is_discrete(self, variable):
         """
         All variables that only can take integer values should be added to this function.
         """
-        if (
-            variable in self.__gas_flow_direct_var
-            or variable in self.__gas_pipe_linear_line_segment_var
-        ):
+        if variable in self.__gas_pipe_linear_line_segment_var:
             return True
         else:
             return super().variable_is_discrete(variable)
@@ -366,8 +349,6 @@ def variable_nominal(self, variable):
 
         if variable in self.__gas_pipe_head_loss_nominals:
             return self.__gas_pipe_head_loss_nominals[variable]
-        elif variable in self.__gas_storage_discharge_nominals:
-            return self.__gas_storage_discharge_nominals[variable]
         elif variable in self.__gas_node_variable_nominal:
             return self.__gas_node_variable_nominal[variable]
         else:
@@ -606,7 +587,7 @@ def __gas_storage_discharge_path_constraints(self, ensemble_member):
         if options["gas_storage_discharge_variables"]:
             for storage in self.energy_system_components.get("gas_tank_storage", []):
                 storage_charge_var = self.state(f"{storage}.GasIn.Q")
-                storage_discharge_var_name = f"{storage}__Q_discharge"
+                storage_discharge_var_name = f"{storage}.__Q_discharge"
                 storage_discharge_var = self.state(storage_discharge_var_name)
                 nominal = self.variable_nominal(storage_discharge_var_name)