diff --git a/examples/municipality/src/run_municipality.py b/examples/municipality/src/run_municipality.py
index 153ab70f1..27691cb1b 100644
--- a/examples/municipality/src/run_municipality.py
+++ b/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
 
diff --git a/examples/pipe_diameter_sizing/src/example.py b/examples/pipe_diameter_sizing/src/example.py
index 95a0ebf81..de7418062 100644
--- a/examples/pipe_diameter_sizing/src/example.py
+++ b/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
 
 
diff --git a/src/mesido/component_type_mixin.py b/src/mesido/component_type_mixin.py
index e48936491..f68d814fe 100644
--- a/src/mesido/component_type_mixin.py
+++ b/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 = {}
 
diff --git a/src/mesido/electricity_physics_mixin.py b/src/mesido/electricity_physics_mixin.py
index 14c16347b..a0b619e5e 100644
--- a/src/mesido/electricity_physics_mixin.py
+++ b/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.
diff --git a/src/mesido/esdl/esdl_heat_model.py b/src/mesido/esdl/esdl_heat_model.py
index bd60d9dc7..7102dc8b7 100644
--- a/src/mesido/esdl/esdl_heat_model.py
+++ b/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
diff --git a/src/mesido/esdl/esdl_mixin.py b/src/mesido/esdl/esdl_mixin.py
index 10a21743e..e827c9bc7 100644
--- a/src/mesido/esdl/esdl_mixin.py
+++ b/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)
 
diff --git a/src/mesido/gas_physics_mixin.py b/src/mesido/gas_physics_mixin.py
index b9efb5f51..ebb8129a1 100644
--- a/src/mesido/gas_physics_mixin.py
+++ b/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,10 +329,7 @@ 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
 
@@ -351,10 +337,7 @@ 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)
 
diff --git a/src/mesido/heat_physics_mixin.py b/src/mesido/heat_physics_mixin.py
index eb4c7cfba..db91e0960 100644
--- a/src/mesido/heat_physics_mixin.py
+++ b/src/mesido/heat_physics_mixin.py
@@ -118,30 +118,15 @@ def __init__(self, *args, **kwargs):
         self._hn_pipe_to_head_loss_map = {}
 
         # Boolean path-variable for the direction of the flow, inport to outport is positive flow.
-        self.__heat_flow_direct_var = {}
         self.__heat_flow_direct_bounds = {}
         self._heat_pipe_to_flow_direct_map = {}
 
         # Boolean path-variable to determine whether flow is going through a pipe.
-        self.__heat_pipe_disconnect_var = {}
-        self.__heat_pipe_disconnect_var_bounds = {}
         self._heat_pipe_disconnect_map = {}
 
-        # Boolean path-variable for the status of the check valve
-        self.__check_valve_status_var = {}
-        self.__check_valve_status_var_bounds = {}
-        self.__check_valve_status_map = {}
-
         # Boolean path-variable for the status of the control valve
-        self.__control_valve_direction_var = {}
-        self.__control_valve_direction_var_bounds = {}
         self.__control_valve_direction_map = {}
 
-        # Boolean path-variable to disable the hex for certain moments in time
-        self.__disabled_hex_map = {}
-        self.__disabled_hex_var = {}
-        self.__disabled_hex_var_bounds = {}
-
         # To avoid artificial energy generation at t0
         self.__buffer_t0_bounds = {}
 
@@ -252,7 +237,6 @@ def _get_min_bound(bound):
 
         # Set structure used instead of list. Purpose: to make lookup faster when there are many
         # pipes in the network.
-        set_self_hot_pipes = set(self.hot_pipes)
 
         for pipe_name in self.energy_system_components.get("heat_pipe", []):
             commodity = self.energy_system_components_commodity.get(pipe_name)
@@ -295,14 +279,9 @@ def _get_min_bound(bound):
                         pipe_linear_line_segment_var_name
                     ] = initialized_vars[10][pipe_linear_line_segment_var_name]
 
-            neighbour = self.has_related_pipe(pipe_name)
-            if neighbour and pipe_name not in set_self_hot_pipes:
-                flow_dir_var = f"{self.cold_to_hot_pipe(pipe_name)}__flow_direct_var"
-            else:
-                flow_dir_var = f"{pipe_name}__flow_direct_var"
+            flow_dir_var = f"{pipe_name}.__flow_direct_var"
 
             self._heat_pipe_to_flow_direct_map[pipe_name] = flow_dir_var
-            self.__heat_flow_direct_var[flow_dir_var] = ca.MX.sym(flow_dir_var)
 
             # Fix the directions that are already implied by the bounds on heat
             # Nonnegative heat implies that flow direction Boolean is equal to one.
@@ -321,46 +300,19 @@ def _get_min_bound(bound):
                 heat_out_lb <= 0.0 and heat_out_ub <= 0.0
             ):
                 self.__heat_flow_direct_bounds[flow_dir_var] = (0.0, 0.0)
-            else:
-                self.__heat_flow_direct_bounds[flow_dir_var] = (0.0, 1.0)
+            # else:
+            #     self.__heat_flow_direct_bounds[flow_dir_var] = (0.0, 1.0)
 
             if parameters[f"{pipe_name}.disconnectable"]:
-                neighbour = self.has_related_pipe(pipe_name)
-                if neighbour and pipe_name not in set_self_hot_pipes:
-                    disconnected_var = f"{self.cold_to_hot_pipe(pipe_name)}__is_disconnected"
-                else:
-                    disconnected_var = f"{pipe_name}__is_disconnected"
-
+                disconnected_var = f"{pipe_name}.__is_disconnected"
                 self._heat_pipe_disconnect_map[pipe_name] = disconnected_var
-                self.__heat_pipe_disconnect_var[disconnected_var] = ca.MX.sym(disconnected_var)
-                self.__heat_pipe_disconnect_var_bounds[disconnected_var] = (0.0, 1.0)
 
             if heat_in_ub <= 0.0 and heat_out_lb >= 0.0:
                 raise Exception(f"Heat flow rate in/out of pipe '{pipe_name}' cannot be zero.")
 
-        # Integers for disabling the HEX temperature constraints
-        for hex in [
-            *self.energy_system_components.get("heat_exchanger", []),
-            *self.energy_system_components.get("heat_pump", []),
-        ]:
-            disabeld_hex_var = f"{hex}__disabled"
-            self.__disabled_hex_map[hex] = disabeld_hex_var
-            self.__disabled_hex_var[disabeld_hex_var] = ca.MX.sym(disabeld_hex_var)
-            self.__disabled_hex_var_bounds[disabeld_hex_var] = (0, 1.0)
-
-        for v in self.energy_system_components.get("check_valve", []):
-            status_var = f"{v}__status_var"
-
-            self.__check_valve_status_map[v] = status_var
-            self.__check_valve_status_var[status_var] = ca.MX.sym(status_var)
-            self.__check_valve_status_var_bounds[status_var] = (0.0, 1.0)
-
         for v in self.energy_system_components.get("control_valve", []):
-            flow_dir_var = f"{v}__flow_direct_var"
-
+            flow_dir_var = f"{v}.__flow_direct_var"
             self.__control_valve_direction_map[v] = flow_dir_var
-            self.__control_valve_direction_var[flow_dir_var] = ca.MX.sym(flow_dir_var)
-            self.__control_valve_direction_var_bounds[flow_dir_var] = (0.0, 1.0)
 
         for ates, (
             (hot_pipe, _hot_pipe_orientation),
@@ -674,17 +626,12 @@ def path_variables(self):
         """
         variables = super().path_variables.copy()
         variables.extend(self.__pipe_head_loss_var.values())
-        variables.extend(self.__heat_flow_direct_var.values())
-        variables.extend(self.__heat_pipe_disconnect_var.values())
-        variables.extend(self.__check_valve_status_var.values())
-        variables.extend(self.__control_valve_direction_var.values())
         variables.extend(self.__demand_insulation_class_var.values())
         variables.extend(self.__pipe_linear_line_segment_var.values())
         variables.extend(self.__temperature_regime_var.values())
         variables.extend(self.__carrier_selected_var.values())
         variables.extend(self.__ates_temperature_disc_var.values())
         variables.extend(self.__ates_temperature_selected_var.values())
-        variables.extend(self.__disabled_hex_var.values())
         variables.extend(self.__pipe_heat_loss_path_var.values())
         variables.extend(self.__ates_temperature_ordering_var.values())
         variables.extend(self.__ates_temperature_disc_ordering_var.values())
@@ -696,15 +643,10 @@ def variable_is_discrete(self, variable):
         All variables that only can take integer values should be added to this function.
         """
         if (
-            variable in self.__heat_flow_direct_var
-            or variable in self.__heat_pipe_disconnect_var
-            or variable in self.__check_valve_status_var
-            or variable in self.__control_valve_direction_var
-            or variable in self.__demand_insulation_class_var
+            variable in self.__demand_insulation_class_var
             or variable in self.__pipe_linear_line_segment_var
             or variable in self.__carrier_selected_var
             or variable in self.__ates_temperature_selected_var
-            or variable in self.__disabled_hex_var
             or variable in self.__ates_temperature_ordering_var
             or variable in self.__ates_temperature_disc_ordering_var
             or variable in self.__carrier_temperature_disc_ordering_var
@@ -735,9 +677,6 @@ def bounds(self):
         """
         bounds = super().bounds()
         bounds.update(self.__heat_flow_direct_bounds)
-        bounds.update(self.__heat_pipe_disconnect_var_bounds)
-        bounds.update(self.__check_valve_status_var_bounds)
-        bounds.update(self.__control_valve_direction_var_bounds)
         bounds.update(self.__buffer_t0_bounds)
         bounds.update(self.__demand_insulation_class_var_bounds)
         bounds.update(self.__pipe_linear_line_segment_var_bounds)
@@ -746,7 +685,6 @@ def bounds(self):
         bounds.update(self.__carrier_selected_var_bounds)
         bounds.update(self.__ates_temperature_disc_var_bounds)
         bounds.update(self.__ates_temperature_selected_var_bounds)
-        bounds.update(self.__disabled_hex_var_bounds)
 
         bounds.update(self.__pipe_head_loss_bounds)
         bounds.update(self.__pipe_head_loss_zero_bounds)
@@ -1263,6 +1201,8 @@ def __flow_direction_path_constraints(self, ensemble_member):
         minimum_velocity = self.heat_network_settings["minimum_velocity"]
         maximum_velocity = self.heat_network_settings["maximum_velocity"]
 
+        set_self_hot_pipes = set(self.hot_pipes)
+
         # Also ensure that the discharge has the same sign as the heat.
         for p in self.energy_system_components.get("heat_pipe", []):
             flow_dir_var = self._heat_pipe_to_flow_direct_map[p]
@@ -1275,6 +1215,16 @@ def __flow_direction_path_constraints(self, ensemble_member):
             else:
                 is_disconnected = self.state(is_disconnected_var)
 
+            if self.has_related_pipe(p) and p not in set_self_hot_pipes:
+                hot_pipe = self.cold_to_hot_pipe(p)
+
+                hot_pipe_flow_var = self.state(f"{hot_pipe}.__flow_direct_var")
+                constraints.append((flow_dir - hot_pipe_flow_var, 0.0, 0.0))
+
+                if is_disconnected_var is not None:
+                    disconnected_var_hot = self.state(f"{hot_pipe}.__is_disconnected")
+                    constraints.append((is_disconnected - disconnected_var_hot, 0.0, 0.0))
+
             q_pipe = self.state(f"{p}.Q")
             heat_in = self.state(f"{p}.HeatIn.Heat")
             heat_out = self.state(f"{p}.HeatOut.Heat")
@@ -1334,7 +1284,7 @@ def __flow_direction_path_constraints(self, ensemble_member):
                     np.inf,
                 )
             )
-            big_m = 2.0 * np.max(
+            big_m = np.max(
                 np.abs(
                     (
                         *self.bounds()[f"{p}.HeatIn.Heat"],
@@ -1655,6 +1605,9 @@ def __pipe_heat_to_discharge_path_constraints(self, ensemble_member):
             carrier = parameters[f"{p}.carrier_id"]
             temperatures = self.temperature_regimes(carrier)
 
+            # TODO: flowdir can be 1 or 0 when Q==0.0, so heat needs to be explicitely set to be
+            # negative or possitive based on flowdir
+
             for heat in [scaled_heat_in, scaled_heat_out]:
                 if self.energy_system_options()["neglect_pipe_heat_losses"]:
                     temp = parameters[f"{p}.temperature"]
@@ -2892,7 +2845,7 @@ def __constraints_temperature_heat_to_discharge_bypass_secondary(
             # Getting var for disabled constraints
             small_m = 0  # 0W
             tol = 1e-5 * big_m  # W
-            is_disabled = self.state(self.__disabled_hex_map[heat_exchanger])
+            is_disabled = self.state(f"{heat_exchanger}.__disabled")
 
             # Constraints to set the disabled integer, note we only set it for the primary
             # side as the secondary side implicetly follows from the energy balance constraints.
@@ -3078,7 +3031,7 @@ def __check_valve_head_discharge_path_constraints(self, ensemble_member):
         maximum_velocity = self.heat_network_settings["maximum_velocity"]
 
         for v in self.energy_system_components.get("check_valve", []):
-            status_var = self.__check_valve_status_map[v]
+            status_var = f"{v}.__status_var"
             status = self.state(status_var)
 
             q = self.state(f"{v}.Q")
diff --git a/src/mesido/pycml/__init__.py b/src/mesido/pycml/__init__.py
index c1192f04e..dcf60bf83 100644
--- a/src/mesido/pycml/__init__.py
+++ b/src/mesido/pycml/__init__.py
@@ -3,6 +3,7 @@
     Connector,
     ConstantInput,
     ControlInput,
+    DiscreteVariable,
     FlattenedModel,
     Model,
     SymbolicParameter,
@@ -14,6 +15,7 @@
     "Connector",
     "ConstantInput",
     "ControlInput",
+    "DiscreteVariable",
     "FlattenedModel",
     "Model",
     "SymbolicParameter",
diff --git a/src/mesido/pycml/component_library/milp/electricity/electricity_storage.py b/src/mesido/pycml/component_library/milp/electricity/electricity_storage.py
index 5f3a22fa7..250fc08c6 100644
--- a/src/mesido/pycml/component_library/milp/electricity/electricity_storage.py
+++ b/src/mesido/pycml/component_library/milp/electricity/electricity_storage.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
 
 from numpy import nan
@@ -32,6 +32,7 @@ def __init__(self, name, **modifiers):
         self.min_voltage = nan
         self.charge_efficiency = 1.0
         self.discharge_efficiency = 1.0
+        self.discharge_var = False
 
         self.add_variable(ElectricityPort, "ElectricityIn")
 
@@ -50,6 +51,16 @@ def __init__(self, name, **modifiers):
             nominal=self.ElectricityIn.Power.nominal,
             max=self.ElectricityIn.Power.max,
         )
+        self.add_variable(DiscreteVariable, "__is_charging", min=0.0, max=1.0)
+
+        if self.discharge_var:
+            self.add_variable(
+                Variable,
+                "__effective_power_discharging",
+                nominal=self.ElectricityIn.Power.nominal,
+                min=0.0,
+                max=-self.ElectricityIn.Power.min,
+            )
 
         self.add_equation(
             (
diff --git a/src/mesido/pycml/component_library/milp/gas/gas_pipe.py b/src/mesido/pycml/component_library/milp/gas/gas_pipe.py
index c9f4ddece..c487c6374 100644
--- a/src/mesido/pycml/component_library/milp/gas/gas_pipe.py
+++ b/src/mesido/pycml/component_library/milp/gas/gas_pipe.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
 
 from numpy import nan, pi
@@ -44,6 +44,8 @@ def __init__(self, name, **modifiers):
         self.add_variable(Variable, "dH")
         self.add_variable(Variable, "Q", nominal=self.Q_nominal)
 
+        self.add_variable(DiscreteVariable, "__gas_flow_direct_var", min=0.0, max=1.0)
+
         # Flow should be preserved
         self.add_equation(((self.GasIn.Q - self.GasOut.Q) / self.Q_nominal))
         self.add_equation(((self.Q - self.GasOut.Q) / self.Q_nominal))
diff --git a/src/mesido/pycml/component_library/milp/gas/gas_tank_storage.py b/src/mesido/pycml/component_library/milp/gas/gas_tank_storage.py
index 7e791bef4..74a8b46ec 100644
--- a/src/mesido/pycml/component_library/milp/gas/gas_tank_storage.py
+++ b/src/mesido/pycml/component_library/milp/gas/gas_tank_storage.py
@@ -26,6 +26,7 @@ def __init__(self, name, **modifiers):
         super().__init__(name, **modifiers)
 
         self.component_type = "gas_tank_storage"
+        self.discharge_var = False
 
         self.min_head = 30.0
         self.density = 2.5e3  # H2 density [g/m3] at 30bar
@@ -47,6 +48,14 @@ def __init__(self, name, **modifiers):
             max=self.density_max_storage * self.volume,
             nominal=self._nominal_stored_gas,
         )
+        if self.discharge_var:
+            self.add_variable(
+                Variable,
+                "__Q_discharge",
+                nominal=self.GasIn.Q.nominal,
+                min=0.0,
+                max=-self.GasIn.Q.min,
+            )
 
         self.add_equation(
             ((self.GasIn.mass_flow - self.Gas_tank_flow) / (self.Q_nominal * self.density))
diff --git a/src/mesido/pycml/component_library/milp/heat/check_valve.py b/src/mesido/pycml/component_library/milp/heat/check_valve.py
index e34acef6d..959e62748 100644
--- a/src/mesido/pycml/component_library/milp/heat/check_valve.py
+++ b/src/mesido/pycml/component_library/milp/heat/check_valve.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
 
 from ._non_storage_component import _NonStorageComponent
@@ -34,6 +34,8 @@ def __init__(self, name, **modifiers):
 
         self.add_variable(Variable, "dH", min=0.0)
 
+        self.add_variable(DiscreteVariable, "__status_var", min=0.0, max=1.0)
+
         self.add_equation(self.dH - (self.HeatOut.H - self.HeatIn.H))
 
         self.add_equation((self.HeatIn.Heat - self.HeatOut.Heat) / self.Heat_nominal)
diff --git a/src/mesido/pycml/component_library/milp/heat/control_valve.py b/src/mesido/pycml/component_library/milp/heat/control_valve.py
index 324d86a4a..ee07d48b1 100644
--- a/src/mesido/pycml/component_library/milp/heat/control_valve.py
+++ b/src/mesido/pycml/component_library/milp/heat/control_valve.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
 
 from ._non_storage_component import _NonStorageComponent
@@ -26,6 +26,8 @@ def __init__(self, name, **modifiers):
 
         self.add_variable(Variable, "dH")
 
+        self.add_variable(DiscreteVariable, "__flow_direct_var", min=0.0, max=1.0)
+
         self.add_equation(self.dH - (self.HeatOut.H - self.HeatIn.H))
         self.add_equation(
             (self.HeatOut.Hydraulic_power - self.HeatIn.Hydraulic_power)
diff --git a/src/mesido/pycml/component_library/milp/heat/heat_exchanger.py b/src/mesido/pycml/component_library/milp/heat/heat_exchanger.py
index a8c45231e..13e8c7a10 100644
--- a/src/mesido/pycml/component_library/milp/heat/heat_exchanger.py
+++ b/src/mesido/pycml/component_library/milp/heat/heat_exchanger.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.component_library.milp._internal import BaseAsset
 from mesido.pycml.component_library.milp.heat.heat_four_port import HeatFourPort
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
@@ -58,6 +58,7 @@ def __init__(self, name, **modifiers):
         self.add_variable(Variable, "Heat_flow", nominal=self.nominal)
         self.add_variable(Variable, "dH_prim")
         self.add_variable(Variable, "dH_sec")
+        self.add_variable(DiscreteVariable, "__disabled", min=0.0, max=1.0)
 
         # Hydraulically decoupled so Heads remain the same
         self.add_equation(self.dH_prim - (self.Primary.HeatOut.H - self.Primary.HeatIn.H))
diff --git a/src/mesido/pycml/component_library/milp/heat/heat_pipe.py b/src/mesido/pycml/component_library/milp/heat/heat_pipe.py
index b1d345e25..1adde8dc8 100644
--- a/src/mesido/pycml/component_library/milp/heat/heat_pipe.py
+++ b/src/mesido/pycml/component_library/milp/heat/heat_pipe.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
 
 from numpy import nan, pi
@@ -61,6 +61,11 @@ def __init__(self, name, **modifiers):
 
         self.add_variable(Variable, "dH")
 
+        self.add_variable(DiscreteVariable, "__flow_direct_var", min=0.0, max=1.0)
+
+        if self.disconnectable:
+            self.add_variable(DiscreteVariable, "__is_disconnected", min=0.0, max=1.0)
+
         # rho * ff * length * area / 2 / diameter * velocity**3
         ff = 0.02  # Order of magnitude expected with 0.05-2.5m/s in 20mm-1200mm diameter pipe
         velo = self.Q_nominal / self.area
diff --git a/src/mesido/pycml/component_library/milp/heat/heat_pump.py b/src/mesido/pycml/component_library/milp/heat/heat_pump.py
index 23d4e7ed2..369b7581e 100644
--- a/src/mesido/pycml/component_library/milp/heat/heat_pump.py
+++ b/src/mesido/pycml/component_library/milp/heat/heat_pump.py
@@ -1,4 +1,4 @@
-from mesido.pycml import Variable
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.component_library.milp._internal import BaseAsset
 from mesido.pycml.component_library.milp.heat.heat_four_port import HeatFourPort
 from mesido.pycml.pycml_mixin import add_variables_documentation_automatically
@@ -54,6 +54,7 @@ def __init__(self, name, **modifiers):
         self.add_variable(Variable, "Power_elec", min=0.0)
         self.add_variable(Variable, "dH_prim", max=0.0)
         self.add_variable(Variable, "dH_sec", min=0.0)
+        self.add_variable(DiscreteVariable, "__disabled", min=0.0, max=1.0)
 
         # Hydraulically decoupled so Heads remain the same
         # #TODO: can't these two equations be moved to the non_storagecomponent?
diff --git a/src/mesido/pycml/component_library/milp/multicommodity/electrolyzer.py b/src/mesido/pycml/component_library/milp/multicommodity/electrolyzer.py
index 08eade878..4507aea12 100644
--- a/src/mesido/pycml/component_library/milp/multicommodity/electrolyzer.py
+++ b/src/mesido/pycml/component_library/milp/multicommodity/electrolyzer.py
@@ -1,4 +1,5 @@
-from mesido.pycml import Variable
+from mesido.electricity_physics_mixin import ElectrolyzerOption
+from mesido.pycml import DiscreteVariable, Variable
 from mesido.pycml.component_library.milp._internal import BaseAsset
 from mesido.pycml.component_library.milp._internal.electricity_component import (
     ElectricityComponent,
@@ -42,6 +43,11 @@ def __init__(self, name, **modifiers):
         self.b_eff_coefficient = nan
         self.c_eff_coefficient = nan
 
+        self.include_asset_is_switched_on = False
+        self.electrolyzer_efficiency_option = (
+            ElectrolyzerOption.LINEARIZED_THREE_LINES_WEAK_INEQUALITY
+        )
+
         self.minimum_load = nan
 
         self.density = 2.5  # H2 density [kg/m3] at 30bar
@@ -67,3 +73,12 @@ def __init__(self, name, **modifiers):
         self.add_equation(
             (self.GasOut.mass_flow - self.Gas_mass_flow_out) / self.nominal_gass_mass_out
         )
+
+        if self.include_asset_is_switched_on:
+            self.add_variable(DiscreteVariable, "__asset_is_switched_on", min=0.0, max=1.0)
+        if (
+            self.electrolyzer_efficiency_option
+            == ElectrolyzerOption.LINEARIZED_THREE_LINES_EQUALITY
+        ):
+            for n in range(3):
+                self.add_variable(DiscreteVariable, f"__line_{n}_active", min=0.0, max=1.0)
diff --git a/src/mesido/pycml/model_base.py b/src/mesido/pycml/model_base.py
index 4d011664f..cea91fccc 100644
--- a/src/mesido/pycml/model_base.py
+++ b/src/mesido/pycml/model_base.py
@@ -108,6 +108,13 @@ def has_derivative(self):
         return hasattr(self, "_derivative")
 
 
+class DiscreteVariable(BaseVariable):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.discrete = True
+        self.python_type = int
+
+
 class ControlInput(BaseVariable):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -201,7 +208,7 @@ def add_variable(self, type_, var_name, *dimensions, **kwargs):
         else:
             var = self._variables[var_name] = type_(f"{self.__prefix}{var_name}", **kwargs)
 
-        if isinstance(var, (Variable, ControlInput, ConstantInput)) and (
+        if isinstance(var, (DiscreteVariable, Variable, ControlInput, ConstantInput)) and (
             isinstance(var.value, (ca.MX, BaseVariable)) or not np.isnan(var.value)
         ):
             # For states and algebraic states, we move the "value" part to an equation
diff --git a/src/mesido/pycml/pycml_mixin.py b/src/mesido/pycml/pycml_mixin.py
index 1579dbe9f..7caffc024 100644
--- a/src/mesido/pycml/pycml_mixin.py
+++ b/src/mesido/pycml/pycml_mixin.py
@@ -95,9 +95,24 @@ def get_names_for_class(current_class_: Type) -> List[str]:
                             dynamic_names_of_port = DYNAMIC_NAME_CACHE[port_class_name]
                             for dynamic_name_of_port in dynamic_names_of_port:
                                 dynamic_names.append(f"{port_name}.{dynamic_name_of_port}")
-                        elif node.args[0].id == "Variable":
+                        elif node.args[0].id == "Variable" or node.args[0].id == "DiscreteVariable":
                             # This is a variable for this component, save its name.
-                            dynamic_names.append(node.args[1].value)
+                            # TODO: later we might want to specify the restrictions on the
+                            # variables, discrete/continuous, min/max, or the if statement
+                            # related to it.
+                            if isinstance(node.args[1], ast.Constant):
+                                dynamic_names.append(node.args[1].value)
+                            elif isinstance(node.args[1], ast.JoinedStr):
+                                str_full = ""
+                                for i in range(len(node.args[1].values)):
+                                    str_part = node.args[1].values[i]
+                                    if isinstance(str_part, ast.Constant):
+                                        str_full += str_part.value
+                                    elif isinstance(str_part, ast.FormattedValue):
+                                        str_full += f"{{{str_part.value.id}}}"
+                                    else:
+                                        raise RuntimeError(f"Unknown case:\n{ast.dump(node)}")
+                                dynamic_names.append(str_full)
                         else:
                             raise RuntimeError(f"Unknown case:\n{ast.dump(node)}")
                     else:
diff --git a/src/mesido/workflows/grow_workflow.py b/src/mesido/workflows/grow_workflow.py
index 2044bf01a..537bb4d61 100644
--- a/src/mesido/workflows/grow_workflow.py
+++ b/src/mesido/workflows/grow_workflow.py
@@ -680,9 +680,9 @@ def run_end_scenario_sizing(
             if p in solution.hot_pipes and parameters[f"{p}.area"] > 0.0:
                 lb = []
                 ub = []
-                bounds_pipe = bounds[f"{p}__flow_direct_var"]
+                bounds_pipe = bounds[f"{p}.__flow_direct_var"]
                 for i in range(len(t)):
-                    r = results[f"{p}__flow_direct_var"][i]
+                    r = results[f"{p}.__flow_direct_var"][i]
                     # bound to roughly represent 4km of milp losses in pipes
                     lb.append(
                         r
@@ -695,7 +695,7 @@ def run_end_scenario_sizing(
                         else bounds_pipe[1]
                     )
 
-                boolean_bounds[f"{p}__flow_direct_var"] = (Timeseries(t, lb), Timeseries(t, ub))
+                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))
diff --git a/src/mesido/workflows/multicommodity_simulator_workflow.py b/src/mesido/workflows/multicommodity_simulator_workflow.py
index a713bd937..e90d69e7a 100644
--- a/src/mesido/workflows/multicommodity_simulator_workflow.py
+++ b/src/mesido/workflows/multicommodity_simulator_workflow.py
@@ -299,10 +299,10 @@ def __create_asset_list_controls(self, asset_types_to_include, assets_without_co
             "electricity_source": "Electricity_source",
             "gas_demand": "Gas_demand_mass_flow",
             "gas_source": "Gas_source_mass_flow",
-            "gas_tank_storage": {"charge": "Gas_tank_flow", "discharge": "__Q_discharge"},
+            "gas_tank_storage": {"charge": "Gas_tank_flow", "discharge": ".__Q_discharge"},
             "electricity_storage": {
                 "charge": "Effective_power_charging",
-                "discharge": "__effective_power_discharging",
+                "discharge": ".__effective_power_discharging",
             },
             "electrolyzer": "Power_consumed",
         }
diff --git a/tests/models/ates_temperature/src/run_ates_temperature.py b/tests/models/ates_temperature/src/run_ates_temperature.py
index 69024ff52..724556978 100644
--- a/tests/models/ates_temperature/src/run_ates_temperature.py
+++ b/tests/models/ates_temperature/src/run_ates_temperature.py
@@ -181,7 +181,7 @@ def path_constraints(self, ensemble_member):
             *self.energy_system_components.get("heat_pump"),
             *self.energy_system_components.get("heat_exchanger"),
         ]:
-            disabled_var = self.state(f"{asset}__disabled")
+            disabled_var = self.state(f"{asset}.__disabled")
             sum_disabled_vars += disabled_var
 
         constraints.append((sum_disabled_vars, 1.0, 2.0))
diff --git a/tests/models/heatpump/src/run_heat_pump.py b/tests/models/heatpump/src/run_heat_pump.py
index f32243ce5..f4334c7e3 100644
--- a/tests/models/heatpump/src/run_heat_pump.py
+++ b/tests/models/heatpump/src/run_heat_pump.py
@@ -124,7 +124,7 @@ def solver_options(self):
         options = super().solver_options()
         options["solver"] = "highs"
         highs_options = options["highs"] = {}
-        highs_options["mip_rel_gap"] = 0.001
+        highs_options["mip_rel_gap"] = 0.0001
         return options
 
     def temperature_carriers(self):
diff --git a/tests/test_electricity_storage.py b/tests/test_electricity_storage.py
index fa0b592eb..4bfd16e41 100644
--- a/tests/test_electricity_storage.py
+++ b/tests/test_electricity_storage.py
@@ -51,9 +51,9 @@ def test_source_sink(self):
         storage_name = solution.energy_system_components.get("electricity_storage")[0]
         charge_eff = parameters[f"{storage_name}.charge_efficiency"]
         discharge_eff = parameters[f"{storage_name}.discharge_efficiency"]
-        is_charging = results[f"{storage_name}__is_charging"]
+        is_charging = results[f"{storage_name}.__is_charging"]
         eff_power_change_bat = results[f"{storage_name}.Effective_power_charging"]
-        eff_power_change_discharge_bat = results[f"{storage_name}__effective_power_discharging"]
+        eff_power_change_discharge_bat = results[f"{storage_name}.__effective_power_discharging"]
         power_bat_network = results[f"{storage_name}.ElectricityIn.Power"]
         stored_el = results[f"{storage_name}.Stored_electricity"]
 
diff --git a/tests/test_electrolyzer.py b/tests/test_electrolyzer.py
index eedc3019c..37af3918f 100644
--- a/tests/test_electrolyzer.py
+++ b/tests/test_electrolyzer.py
@@ -283,7 +283,7 @@ def solver_options(self):
         )
         # Check that the electrolyzer is switched off
         np.testing.assert_allclose(
-            results["Electrolyzer_fc66__asset_is_switched_on"][-1],
+            results["Electrolyzer_fc66.__asset_is_switched_on"][-1],
             0,
         )
         # Check that the input power is greater than 0
@@ -298,7 +298,7 @@ def solver_options(self):
         )
         # Check that the electrolyzer is switched off
         np.testing.assert_allclose(
-            results["Electrolyzer_fc66__asset_is_switched_on"][:-1],
+            results["Electrolyzer_fc66.__asset_is_switched_on"][:-1],
             np.ones(2),
         )
         # Check electrolyzer input power
@@ -395,13 +395,13 @@ def test_electrolyzer_equality_constraint(self):
         results = solution.extract_results()
 
         # Check that there is only one activated line per timestep
-        for timestep in range(len(results["Electrolyzer_fc66__line_0_active"])):
+        for timestep in range(len(results["Electrolyzer_fc66.__line_0_active"])):
             np.testing.assert_allclose(
                 (
-                    results["Electrolyzer_fc66__line_0_active"][timestep]
-                    + results["Electrolyzer_fc66__line_1_active"][timestep]
-                    + results["Electrolyzer_fc66__line_2_active"][timestep]
-                    + (1 - results["Electrolyzer_fc66__asset_is_switched_on"][timestep])
+                    results["Electrolyzer_fc66.__line_0_active"][timestep]
+                    + results["Electrolyzer_fc66.__line_1_active"][timestep]
+                    + results["Electrolyzer_fc66.__line_2_active"][timestep]
+                    + (1 - results["Electrolyzer_fc66.__asset_is_switched_on"][timestep])
                 ),
                 1.0,
             )
@@ -411,8 +411,8 @@ def test_electrolyzer_equality_constraint(self):
         for idx in range(3):
             np.testing.assert_allclose(
                 (
-                    results[f"Electrolyzer_fc66__line_{idx}_active"][idx]
-                    + (1 - results["Electrolyzer_fc66__asset_is_switched_on"][idx])
+                    results[f"Electrolyzer_fc66.__line_{idx}_active"][idx]
+                    + (1 - results["Electrolyzer_fc66.__asset_is_switched_on"][idx])
                 ),
                 1.0,
             )
@@ -477,10 +477,10 @@ def test_electrolyzer_equality_constraint_inactive_line(self):
         # such that no line should be active
         np.testing.assert_allclose(
             (
-                results["Electrolyzer_fc66__line_0_active"][-1]
-                + results["Electrolyzer_fc66__line_1_active"][-1]
-                + results["Electrolyzer_fc66__line_2_active"][-1]
-                + (1 - results["Electrolyzer_fc66__asset_is_switched_on"][-1])
+                results["Electrolyzer_fc66.__line_0_active"][-1]
+                + results["Electrolyzer_fc66.__line_1_active"][-1]
+                + results["Electrolyzer_fc66.__line_2_active"][-1]
+                + (1 - results["Electrolyzer_fc66.__asset_is_switched_on"][-1])
             ),
             1.0,
         )
diff --git a/tests/test_head_loss.py b/tests/test_head_loss.py
index 67ba3b39d..dc5e5e1b5 100644
--- a/tests/test_head_loss.py
+++ b/tests/test_head_loss.py
@@ -439,7 +439,7 @@ def energy_system_options(self):
                     results["Pipe4.HeatIn.Q"]
                     / (solution.parameters(0)["Pipe4.diameter"] ** 2 / 4.0 * np.pi),
                     0.0,
-                    atol=1e-07,
+                    atol=1e-03,
                 )
             elif (
                 solution.heat_network_settings["head_loss_option"]
diff --git a/tests/test_heat.py b/tests/test_heat.py
index b1acaa7cf..0d9c78f6a 100644
--- a/tests/test_heat.py
+++ b/tests/test_heat.py
@@ -321,8 +321,8 @@ def test_disconnected_network_pipe(self):
 
         self.assertLess(q_disconnected[1], q_connected[1])
         self.assertAlmostEqual(q_disconnected[1], 0.0, 5)
-        np.testing.assert_allclose(results_connected["Pipe1__is_disconnected"], 0.0)
-        np.testing.assert_allclose(results_disconnected["Pipe1__is_disconnected"][1], 1.0)
+        np.testing.assert_allclose(results_connected["Pipe1.__is_disconnected"], 0.0)
+        np.testing.assert_allclose(results_disconnected["Pipe1.__is_disconnected"][1], 1.0)
 
         np.testing.assert_allclose(q_connected[2:], q_disconnected[2:])
 
@@ -370,4 +370,4 @@ def test_disconnected_pipe_darcy_weisbach(self):
         # Without any constraints on the maximum or minimum head/pressure
         # (loss) in the system, we expect equal results.
         np.testing.assert_allclose(q_linear, q_dw)
-        np.testing.assert_allclose(results_dw["Pipe1__is_disconnected"][1], 1.0)
+        np.testing.assert_allclose(results_dw["Pipe1.__is_disconnected"][1], 1.0)
diff --git a/tests/test_multicommodity.py b/tests/test_multicommodity.py
index 010537eed..297690a0e 100644
--- a/tests/test_multicommodity.py
+++ b/tests/test_multicommodity.py
@@ -194,11 +194,11 @@ def test_heat_pump_elec_min_elec(self):
         heatpump_power = results["GenericConversion_3d3f.Power_elec"]
         heatpump_heat_prim = results["GenericConversion_3d3f.Primary_heat"]
         heatpump_heat_sec = results["GenericConversion_3d3f.Secondary_heat"]
-        heatpump_disabled = results["GenericConversion_3d3f__disabled"]
+        heatpump_disabled = results["GenericConversion_3d3f.__disabled"]
         # heatdemand_sec = results["HeatingDemand_18aa.Heat_demand"]
         heatdemand_prim = results["HeatingDemand_3322.Heat_demand"]
         elec_prod_power = results["ElectricityProducer_ac2e.ElectricityOut.Power"]
-        # pipe_sec_out_hp_disconnected = results["Pipe_408e__is_disconnected"]
+        # pipe_sec_out_hp_disconnected = results["Pipe_408e.__is_disconnected"]
 
         # check that heatpump is not used:
         np.testing.assert_allclose(heatpump_power, np.zeros(len(heatpump_power)), atol=tol)
@@ -253,10 +253,10 @@ def solver_options(self):
 
         heatpump_power = results["GenericConversion_3d3f.Power_elec"]
         heatpump_heat_sec = results["GenericConversion_3d3f.Secondary_heat"]
-        heatpump_disabled = results["GenericConversion_3d3f__disabled"]
+        heatpump_disabled = results["GenericConversion_3d3f.__disabled"]
         heatdemand_sec = results["HeatingDemand_18aa.Heat_demand"]
         var_opex_hp = results["GenericConversion_3d3f__variable_operational_cost"]
-        # pipe_sec_out_hp_disconnected = results["Pipe_408e__is_disconnected"]
+        # pipe_sec_out_hp_disconnected = results["Pipe_408e.__is_disconnected"]
 
         # check that heatpump is not used when electricity price is high:
         price_profile = solution.get_timeseries("Electr.price_profile").values
diff --git a/tests/test_multicommodity_simulator.py b/tests/test_multicommodity_simulator.py
index f9b5437cc..bf176b737 100644
--- a/tests/test_multicommodity_simulator.py
+++ b/tests/test_multicommodity_simulator.py
@@ -625,8 +625,8 @@ def test_multi_commodity_simulator_sequential_staged(self):
             if len(key.split("__")) > 1 and key.split("__")[1] == "gas_flow_direct_var":
                 # For the scenario when Q is -0 and 0 then gas_flow_direct_var 0 or 1 for the same
                 # volumetric flow rate of zero. So only check gas_flow_direct_var when Q != zero
-                zero_staged = results_staged[f"{key.split('__')[0]}.GasIn.Q"] != 0
-                zero_unstaged = results_unstaged[f"{key.split('__')[0]}.GasIn.Q"] != 0
+                zero_staged = results_staged[f"{key.split('__')[0]}GasIn.Q"] != 0
+                zero_unstaged = results_unstaged[f"{key.split('__')[0]}GasIn.Q"] != 0
                 np.testing.assert_allclose(value[zero_staged], value_staged[zero_unstaged])
             else:
                 np.testing.assert_allclose(value, value_staged)
diff --git a/tests/test_multiple_in_and_out_port_components.py b/tests/test_multiple_in_and_out_port_components.py
index f59a3e700..000d4e634 100644
--- a/tests/test_multiple_in_and_out_port_components.py
+++ b/tests/test_multiple_in_and_out_port_components.py
@@ -75,7 +75,7 @@ def energy_system_options(self):
 
         prim_heat = results["HeatExchange_39ed.Primary_heat"]
         sec_heat = results["HeatExchange_39ed.Secondary_heat"]
-        disabled = results["HeatExchange_39ed__disabled"]
+        disabled = results["HeatExchange_39ed.__disabled"]
 
         # We check the energy converted betweeen the commodities
         eff = parameters["HeatExchange_39ed.efficiency"]
@@ -163,8 +163,8 @@ def energy_system_options(self):
         hex_active = "HeatExchange_e410_copy"
         hex_bypass = "HeatExchange_e410"
 
-        np.testing.assert_allclose(results[f"{hex_active}__disabled"][:-1], 0)
-        np.testing.assert_allclose(results[f"{hex_bypass}__disabled"][:-1], 1)
+        np.testing.assert_allclose(results[f"{hex_active}.__disabled"][:-1], 0)
+        np.testing.assert_allclose(results[f"{hex_bypass}.__disabled"][:-1], 1)
 
         np.testing.assert_array_less(0.001, results[f"{hex_active}.Primary.Q"][:-1])
         np.testing.assert_array_less(0.001, results[f"{hex_bypass}.Primary.Q"][:-1])
@@ -285,8 +285,8 @@ def times(self, variable=None) -> np.ndarray:
         hex_active = "HeatExchange_e410_copy"
         hex_bypass = "HeatExchange_e410"
 
-        np.testing.assert_allclose(results[f"{hex_active}__disabled"][:-1], 0)
-        np.testing.assert_allclose(results[f"{hex_bypass}__disabled"][:-1], 1)
+        np.testing.assert_allclose(results[f"{hex_active}.__disabled"][:-1], 0)
+        np.testing.assert_allclose(results[f"{hex_bypass}.__disabled"][:-1], 1)
 
         np.testing.assert_array_less(0.001, results[f"{hex_active}.Primary.Q"][:-1])
         np.testing.assert_array_less(0.001, results[f"{hex_bypass}.Primary.Q"][:-1])
diff --git a/tests/test_pipe_diameter_sizing.py b/tests/test_pipe_diameter_sizing.py
index ed791f706..d7488fd93 100644
--- a/tests/test_pipe_diameter_sizing.py
+++ b/tests/test_pipe_diameter_sizing.py
@@ -117,6 +117,7 @@ def test_half_network_gone(self):
                     2.0e-4,
                     parameters[f"{pipe}.temperature"],
                 )
+
                 c_v = parameters[f"{pipe}.length"] * ff / (2 * 9.81) / pc.inner_diameter
                 dh_max = c_v * pc.maximum_velocity**2
                 dh_manual = dh_max * results[f"{pipe}.Q"][1:] / pc.area / pc.maximum_velocity
@@ -127,7 +128,9 @@ def test_half_network_gone(self):
         for pipe in diameters.keys():
             if pipe in pipes_removed:
                 hydraulic_power_sum += sum(abs(results[f"{pipe}.Hydraulic_power"]))
-        self.assertEqual(hydraulic_power_sum, 0.0, "Hydraulic power exists for a removed pipe")
+        np.testing.assert_allclose(
+            hydraulic_power_sum, 0.0, err_msg="Hydraulic power exists for a removed pipe", atol=1e-9
+        )
 
         # Hydraulic power = delta pressure * Q = f(Q^3), where delta pressure = f(Q^2)
         # The linear approximation of the 3rd order function should overestimate the hydraulic
diff --git a/tests/test_temperature_ates_hp.py b/tests/test_temperature_ates_hp.py
index ef323f449..9fcd0a212 100644
--- a/tests/test_temperature_ates_hp.py
+++ b/tests/test_temperature_ates_hp.py
@@ -61,7 +61,7 @@ def test_ates_temperature(self):
         energy_conservation_test(solution, results)
         heat_to_discharge_test(solution, results)
 
-        ates_charging = results["Pipe1__flow_direct_var"]  # =1 if charging
+        ates_charging = results["Pipe1.__flow_direct_var"]  # =1 if charging
         ates_temperature = results["ATES_cb47.Temperature_ates"]
         ates_temperature_disc = results["ATES_cb47__temperature_ates_disc"]
         carrier_temperature = results["41770304791669983859190_temperature"]
@@ -74,8 +74,8 @@ def test_ates_temperature(self):
         heat_ates = results["ATES_cb47.Heat_ates"]
         heat_loss_ates = results["ATES_cb47.Heat_loss"]
         ates_stored_heat = results["ATES_cb47.Stored_heat"]
-        hex_disabled = results["HeatExchange_32ba__disabled"]
-        hp_disabled = results["HeatPump_7f2c__disabled"]
+        hex_disabled = results["HeatExchange_32ba.__disabled"]
+        hp_disabled = results["HeatPump_7f2c.__disabled"]
 
         # geo_source = results["GeothermalSource_4e5b.Heat_source"]
         objective = solution.objective_value
diff --git a/tests/test_varying_temperature.py b/tests/test_varying_temperature.py
index 2db1d2488..6328b1c34 100644
--- a/tests/test_varying_temperature.py
+++ b/tests/test_varying_temperature.py
@@ -49,6 +49,8 @@ def test_1a_temperature_variation(self):
         test = TestCase()
         test.assertTrue(heat_problem.solver_stats["success"], msg="Optimisation did not succeed")
         # Check that the highest supply temperature is selected
+        heat_to_discharge_test(heat_problem, results)
+        demand_matching_test(heat_problem, results)
         np.testing.assert_allclose(results[f"{3625334968694477359}_temperature"], 85.0)
         np.testing.assert_allclose(results[f"{3625334968694477359}_85.0"], 1.0)
 
@@ -317,7 +319,7 @@ def test_hex_temperature_variation_disablehex(self):
         # Check that the problem has an infeasible temperature for the hex
         np.testing.assert_allclose(results[f"{33638164429859421}_temperature"], 69.0)
         # Verify that the hex is disabled
-        np.testing.assert_allclose(results["HeatExchange_39ed__disabled"], 1.0)
+        np.testing.assert_allclose(results["HeatExchange_39ed.__disabled"], 1.0)
         np.testing.assert_allclose(results["HeatExchange_39ed.Primary_heat"], 0.0)
 
         demand_matching_test(heat_problem, results)
diff --git a/tests/test_warmingup_unit_cases.py b/tests/test_warmingup_unit_cases.py
index ad1ac9876..52b542b76 100644
--- a/tests/test_warmingup_unit_cases.py
+++ b/tests/test_warmingup_unit_cases.py
@@ -156,9 +156,9 @@ def test_3a(self):
         heat_to_discharge_test(heat_problem, results)
 
         # We only check the flow directions for the time-steps that there is flow in the pipe.
-        inds = np.round(1 - results["Pipe_e53a__is_disconnected"]).astype(bool)
+        inds = np.round(1 - results["Pipe_e53a.__is_disconnected"]).astype(bool)
         np.testing.assert_allclose(
-            np.round(results["Pipe_e53a__flow_direct_var"][inds]) * 2.0 - 1.0,
+            np.round(results["Pipe_e53a.__flow_direct_var"][inds]) * 2.0 - 1.0,
             np.sign(results["Pipe_e53a.Q"][inds]),
         )
 
diff --git a/tests/utils_tests.py b/tests/utils_tests.py
index a77320fcc..385102ed4 100644
--- a/tests/utils_tests.py
+++ b/tests/utils_tests.py
@@ -113,7 +113,7 @@ def heat_to_discharge_test(solution, results):
      discharge and heatflow can be negative.
     """
     test = TestCase()
-    tol = 1.0e-2
+    tol_heat = 50  # W
     for d in [
         *solution.energy_system_components.get("heat_demand", []),
         *solution.energy_system_components.get("airco", []),
@@ -136,7 +136,7 @@ def heat_to_discharge_test(solution, results):
         np.testing.assert_allclose(
             results[f"{d}.Heat_flow"],
             results[f"{d}.HeatIn.Heat"] - results[f"{d}.HeatOut.Heat"],
-            atol=tol,
+            atol=tol_heat,
         )
         np.testing.assert_allclose(
             results[f"{d}.HeatOut.Heat"], results[f"{d}.Q"] * rho * cp * return_t
@@ -149,7 +149,7 @@ def heat_to_discharge_test(solution, results):
         np.testing.assert_allclose(
             results[f"{d}.Cold_demand"],
             results[f"{d}.HeatOut.Heat"] - results[f"{d}.HeatIn.Heat"],
-            atol=tol,
+            atol=tol_heat,
         )
         np.testing.assert_allclose(
             results[f"{d}.HeatOut.Heat"], results[f"{d}.Q"] * rho * cp * supply_t
@@ -167,12 +167,13 @@ def heat_to_discharge_test(solution, results):
         # supply_t = solution.parameters(0)[f"{d}.T_supply"]
         # return_t = solution.parameters(0)[f"{d}.T_return"]
         supply_t, return_t, dt = _get_component_temperatures(solution, results, d)
+        sup_t_for_error = supply_t if isinstance(supply_t, float) else min(supply_t)
 
         print(d, max(abs(results[f"{d}.HeatOut.Heat"] - results[f"{d}.Q"] * rho * cp * supply_t)))
         np.testing.assert_allclose(
             results[f"{d}.HeatOut.Heat"],
             results[f"{d}.Q"] * rho * cp * supply_t,
-            atol=5.0,
+            atol=1.0e-6 * rho * cp * sup_t_for_error,
             rtol=1.0e-4,
         )
 
@@ -194,13 +195,13 @@ def heat_to_discharge_test(solution, results):
             np.testing.assert_allclose(
                 results[f"{d}.Heat_ates"],
                 results[f"{d}.HeatIn.Heat"] - results[f"{d}.HeatOut.Heat"],
-                atol=tol,
+                atol=tol_heat,
             )
         except KeyError:
             np.testing.assert_allclose(
                 results[f"{d}.Heat_buffer"],
                 results[f"{d}.HeatIn.Heat"] - results[f"{d}.HeatOut.Heat"],
-                atol=tol,
+                atol=tol_heat,
             )
         supply_temp, return_temp, dt = _get_component_temperatures(solution, results, d)
         indices = results[f"{d}.HeatIn.Q"] >= 0
@@ -215,13 +216,13 @@ def heat_to_discharge_test(solution, results):
         test.assertTrue(
             expr=all(
                 results[f"{d}.HeatIn.Heat"][indices]
-                <= (results[f"{d}.HeatIn.Q"][indices] * rho * cp * supply_t + tol)
+                <= (results[f"{d}.HeatIn.Q"][indices] * rho * cp * supply_t + tol_heat)
             )
         )
         np.testing.assert_allclose(
             results[f"{d}.HeatOut.Heat"][indices],
             results[f"{d}.HeatIn.Q"][indices] * rho * cp * return_t,
-            atol=tol,
+            atol=tol_heat,
         )
         indices = results[f"{d}.HeatIn.Q"] <= 0
         if isinstance(supply_t, float):
@@ -235,13 +236,13 @@ def heat_to_discharge_test(solution, results):
         np.testing.assert_allclose(
             results[f"{d}.HeatIn.Heat"][indices],
             results[f"{d}.HeatIn.Q"][indices] * rho * cp * supply_t,
-            atol=tol,
+            atol=tol_heat,
         )
 
         test.assertTrue(
             expr=all(
                 results[f"{d}.HeatOut.Heat"][indices]
-                >= (results[f"{d}.HeatIn.Q"][indices] * rho * cp * return_t - tol)
+                >= (results[f"{d}.HeatIn.Q"][indices] * rho * cp * return_t - tol_heat)
             )
         )
 
@@ -269,13 +270,13 @@ def heat_to_discharge_test(solution, results):
             heat = results[f"{d}.{p}_heat"]
 
             if p == "Primary":
-                np.testing.assert_allclose(heat_out, discharge * rho * cp * return_t, atol=tol)
-                test.assertTrue(expr=all(heat_in <= discharge * rho * cp * supply_t + tol))
-                test.assertTrue(expr=all(heat <= discharge * rho * cp * dt + tol))
+                np.testing.assert_allclose(heat_out, discharge * rho * cp * return_t, atol=tol_heat)
+                test.assertTrue(expr=all(heat_in <= discharge * rho * cp * supply_t + tol_heat))
+                test.assertTrue(expr=all(heat <= discharge * rho * cp * dt + tol_heat))
             elif p == "Secondary":
-                test.assertTrue(expr=all(heat >= discharge * rho * cp * dt - tol))
-                np.testing.assert_allclose(heat_out, discharge * rho * cp * supply_t, atol=tol)
-                test.assertTrue(expr=all(heat_in <= discharge * rho * cp * return_t + tol))
+                test.assertTrue(expr=all(heat >= discharge * rho * cp * dt - tol_heat))
+                np.testing.assert_allclose(heat_out, discharge * rho * cp * supply_t, atol=tol_heat)
+                test.assertTrue(expr=all(heat_in <= discharge * rho * cp * return_t + tol_heat))
 
     for p in solution.energy_system_components.get("heat_pipe", []):
         cp = solution.parameters(0)[f"{p}.cp"]
@@ -301,7 +302,7 @@ def heat_to_discharge_test(solution, results):
         test.assertTrue(
             expr=all(
                 results[f"{p}.HeatOut.Heat"][indices]
-                <= results[f"{p}.Q"][indices] * rho * cp * temperature + tol
+                <= results[f"{p}.Q"][indices] * rho * cp * temperature + tol_heat
             )
         )
         indices = results[f"{p}.Q"] < 0
@@ -318,13 +319,13 @@ def heat_to_discharge_test(solution, results):
         test.assertTrue(
             expr=all(
                 results[f"{p}.HeatIn.Heat"][indices]
-                >= results[f"{p}.Q"][indices] * rho * cp * temperature - tol
+                >= results[f"{p}.Q"][indices] * rho * cp * temperature - tol_heat
             )
         )
         test.assertTrue(
             expr=all(
                 results[f"{p}.HeatOut.Heat"][indices]
-                >= results[f"{p}.Q"][indices] * rho * cp * temperature - tol
+                >= results[f"{p}.Q"][indices] * rho * cp * temperature - tol_heat
             )
         )
         indices = results[f"{p}.Q"] == 0
@@ -335,13 +336,13 @@ def heat_to_discharge_test(solution, results):
         np.testing.assert_allclose(
             results[f"{p}.HeatIn.Heat"][indices],
             results[f"{p}.Q"][indices] * rho * cp * temperature,
-            atol=tol,
+            atol=tol_heat,
             err_msg=f"{p} has mismatch in milp to discharge",
         )
         np.testing.assert_allclose(
             results[f"{p}.HeatOut.Heat"][indices],
             results[f"{p}.Q"][indices] * rho * cp * temperature,
-            atol=tol,
+            atol=tol_heat,
             err_msg=f"{p} has mismatch in milp to discharge",
         )
 
@@ -437,13 +438,13 @@ def energy_conservation_test(solution, results):
 
     for p in solution.energy_system_components.get("heat_pipe", []):
         if (
-            f"{p}__is_disconnected" in results.keys()
-            or f"{solution.cold_to_hot_pipe(p)}__is_disconnected" in results.keys()
+            f"{p}.__is_disconnected" in results.keys()
+            or f"{solution.cold_to_hot_pipe(p)}.__is_disconnected" in results.keys()
         ):
             if p in solution.hot_pipes:
-                p_discon = results[f"{p}__is_disconnected"].copy()
+                p_discon = results[f"{p}.__is_disconnected"].copy()
             else:
-                p_discon = results[f"{solution.cold_to_hot_pipe(p)}__is_disconnected"].copy()
+                p_discon = results[f"{solution.cold_to_hot_pipe(p)}.__is_disconnected"].copy()
 
             p_discon[p_discon < 0.5] = 0  # fix for discrete value sometimes being 0.003 or so.
             np.testing.assert_allclose(