From 86ea2aed00b8790a9cc883fcfd1b7a9f44024b8a Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 10:24:25 +0100
Subject: [PATCH 01/20] refactoring to tighten the coupling between FlowSystem
 and Elements

---
 flixopt/components.py  | 50 +++++++++++++++---------------
 flixopt/effects.py     | 26 +++++++++-------
 flixopt/elements.py    | 38 +++++++++++++----------
 flixopt/flow_system.py | 31 ++++++++++++++++++-
 flixopt/interface.py   | 70 +++++++++++++++++++++---------------------
 flixopt/structure.py   | 55 +++++++++++++++++++++++++++++++--
 6 files changed, 180 insertions(+), 90 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index e4209c8ac..54b5e4ea1 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -212,23 +212,25 @@ def _plausibility_checks(self) -> None:
                         f'({flow.label_full}).'
                     )
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        super().transform_data(flow_system, prefix)
+        super().transform_data(prefix)
         if self.conversion_factors:
-            self.conversion_factors = self._transform_conversion_factors(flow_system)
+            self.conversion_factors = self._transform_conversion_factors()
         if self.piecewise_conversion:
             self.piecewise_conversion.has_time_dim = True
-            self.piecewise_conversion.transform_data(flow_system, f'{prefix}|PiecewiseConversion')
+            self.piecewise_conversion.transform_data(f'{prefix}|PiecewiseConversion')
 
-    def _transform_conversion_factors(self, flow_system: FlowSystem) -> list[dict[str, xr.DataArray]]:
+    def _transform_conversion_factors(self) -> list[dict[str, xr.DataArray]]:
         """Converts all conversion factors to internal datatypes"""
         list_of_conversion_factors = []
         for idx, conversion_factor in enumerate(self.conversion_factors):
             transformed_dict = {}
             for flow, values in conversion_factor.items():
                 # TODO: Might be better to use the label of the component instead of the flow
-                ts = flow_system.fit_to_model_coords(f'{self.flows[flow].label_full}|conversion_factor{idx}', values)
+                ts = self.flow_system.fit_to_model_coords(
+                    f'{self.flows[flow].label_full}|conversion_factor{idx}', values
+                )
                 if ts is None:
                     raise PlausibilityError(f'{self.label_full}: conversion factor for flow "{flow}" must not be None')
                 transformed_dict[flow] = ts
@@ -434,46 +436,46 @@ def create_model(self, model: FlowSystemModel) -> StorageModel:
         self.submodel = StorageModel(model, self)
         return self.submodel
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        super().transform_data(flow_system, prefix)
-        self.relative_minimum_charge_state = flow_system.fit_to_model_coords(
+        super().transform_data(prefix)
+        self.relative_minimum_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|relative_minimum_charge_state',
             self.relative_minimum_charge_state,
         )
-        self.relative_maximum_charge_state = flow_system.fit_to_model_coords(
+        self.relative_maximum_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|relative_maximum_charge_state',
             self.relative_maximum_charge_state,
         )
-        self.eta_charge = flow_system.fit_to_model_coords(f'{prefix}|eta_charge', self.eta_charge)
-        self.eta_discharge = flow_system.fit_to_model_coords(f'{prefix}|eta_discharge', self.eta_discharge)
-        self.relative_loss_per_hour = flow_system.fit_to_model_coords(
+        self.eta_charge = self.flow_system.fit_to_model_coords(f'{prefix}|eta_charge', self.eta_charge)
+        self.eta_discharge = self.flow_system.fit_to_model_coords(f'{prefix}|eta_discharge', self.eta_discharge)
+        self.relative_loss_per_hour = self.flow_system.fit_to_model_coords(
             f'{prefix}|relative_loss_per_hour', self.relative_loss_per_hour
         )
         if not isinstance(self.initial_charge_state, str):
-            self.initial_charge_state = flow_system.fit_to_model_coords(
+            self.initial_charge_state = self.flow_system.fit_to_model_coords(
                 f'{prefix}|initial_charge_state', self.initial_charge_state, dims=['period', 'scenario']
             )
-        self.minimal_final_charge_state = flow_system.fit_to_model_coords(
+        self.minimal_final_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|minimal_final_charge_state', self.minimal_final_charge_state, dims=['period', 'scenario']
         )
-        self.maximal_final_charge_state = flow_system.fit_to_model_coords(
+        self.maximal_final_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|maximal_final_charge_state', self.maximal_final_charge_state, dims=['period', 'scenario']
         )
-        self.relative_minimum_final_charge_state = flow_system.fit_to_model_coords(
+        self.relative_minimum_final_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|relative_minimum_final_charge_state',
             self.relative_minimum_final_charge_state,
             dims=['period', 'scenario'],
         )
-        self.relative_maximum_final_charge_state = flow_system.fit_to_model_coords(
+        self.relative_maximum_final_charge_state = self.flow_system.fit_to_model_coords(
             f'{prefix}|relative_maximum_final_charge_state',
             self.relative_maximum_final_charge_state,
             dims=['period', 'scenario'],
         )
         if isinstance(self.capacity_in_flow_hours, InvestParameters):
-            self.capacity_in_flow_hours.transform_data(flow_system, f'{prefix}|InvestParameters')
+            self.capacity_in_flow_hours.transform_data(f'{prefix}|InvestParameters')
         else:
-            self.capacity_in_flow_hours = flow_system.fit_to_model_coords(
+            self.capacity_in_flow_hours = self.flow_system.fit_to_model_coords(
                 f'{prefix}|capacity_in_flow_hours', self.capacity_in_flow_hours, dims=['period', 'scenario']
             )
 
@@ -720,11 +722,11 @@ def create_model(self, model) -> TransmissionModel:
         self.submodel = TransmissionModel(model, self)
         return self.submodel
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        super().transform_data(flow_system, prefix)
-        self.relative_losses = flow_system.fit_to_model_coords(f'{prefix}|relative_losses', self.relative_losses)
-        self.absolute_losses = flow_system.fit_to_model_coords(f'{prefix}|absolute_losses', self.absolute_losses)
+        super().transform_data(prefix)
+        self.relative_losses = self.flow_system.fit_to_model_coords(f'{prefix}|relative_losses', self.relative_losses)
+        self.absolute_losses = self.flow_system.fit_to_model_coords(f'{prefix}|absolute_losses', self.absolute_losses)
 
 
 class TransmissionModel(ComponentModel):
diff --git a/flixopt/effects.py b/flixopt/effects.py
index ddf8eadeb..4b0a6b1d1 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -340,43 +340,47 @@ def maximum_operation_per_hour(self, value):
         )
         self.maximum_per_hour = value
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.minimum_per_hour = flow_system.fit_to_model_coords(f'{prefix}|minimum_per_hour', self.minimum_per_hour)
+        self.minimum_per_hour = self.flow_system.fit_to_model_coords(
+            f'{prefix}|minimum_per_hour', self.minimum_per_hour
+        )
 
-        self.maximum_per_hour = flow_system.fit_to_model_coords(f'{prefix}|maximum_per_hour', self.maximum_per_hour)
+        self.maximum_per_hour = self.flow_system.fit_to_model_coords(
+            f'{prefix}|maximum_per_hour', self.maximum_per_hour
+        )
 
-        self.share_from_temporal = flow_system.fit_effects_to_model_coords(
+        self.share_from_temporal = self.flow_system.fit_effects_to_model_coords(
             label_prefix=None,
             effect_values=self.share_from_temporal,
             label_suffix=f'(temporal)->{prefix}(temporal)',
             dims=['time', 'period', 'scenario'],
         )
-        self.share_from_periodic = flow_system.fit_effects_to_model_coords(
+        self.share_from_periodic = self.flow_system.fit_effects_to_model_coords(
             label_prefix=None,
             effect_values=self.share_from_periodic,
             label_suffix=f'(periodic)->{prefix}(periodic)',
             dims=['period', 'scenario'],
         )
 
-        self.minimum_temporal = flow_system.fit_to_model_coords(
+        self.minimum_temporal = self.flow_system.fit_to_model_coords(
             f'{prefix}|minimum_temporal', self.minimum_temporal, dims=['period', 'scenario']
         )
-        self.maximum_temporal = flow_system.fit_to_model_coords(
+        self.maximum_temporal = self.flow_system.fit_to_model_coords(
             f'{prefix}|maximum_temporal', self.maximum_temporal, dims=['period', 'scenario']
         )
-        self.minimum_periodic = flow_system.fit_to_model_coords(
+        self.minimum_periodic = self.flow_system.fit_to_model_coords(
             f'{prefix}|minimum_periodic', self.minimum_periodic, dims=['period', 'scenario']
         )
-        self.maximum_periodic = flow_system.fit_to_model_coords(
+        self.maximum_periodic = self.flow_system.fit_to_model_coords(
             f'{prefix}|maximum_periodic', self.maximum_periodic, dims=['period', 'scenario']
         )
-        self.minimum_total = flow_system.fit_to_model_coords(
+        self.minimum_total = self.flow_system.fit_to_model_coords(
             f'{prefix}|minimum_total',
             self.minimum_total,
             dims=['period', 'scenario'],
         )
-        self.maximum_total = flow_system.fit_to_model_coords(
+        self.maximum_total = self.flow_system.fit_to_model_coords(
             f'{prefix}|maximum_total', self.maximum_total, dims=['period', 'scenario']
         )
 
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 337f34fce..aaa3e6cb4 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -100,13 +100,13 @@ def create_model(self, model: FlowSystemModel) -> ComponentModel:
         self.submodel = ComponentModel(model, self)
         return self.submodel
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         if self.on_off_parameters is not None:
-            self.on_off_parameters.transform_data(flow_system, prefix)
+            self.on_off_parameters.transform_data(prefix)
 
         for flow in self.inputs + self.outputs:
-            flow.transform_data(flow_system)  # Flow doesnt need the name_prefix
+            flow.transform_data()  # Flow doesnt need the name_prefix
 
     def _check_unique_flow_labels(self):
         all_flow_labels = [flow.label for flow in self.inputs + self.outputs]
@@ -241,9 +241,9 @@ def create_model(self, model: FlowSystemModel) -> BusModel:
         self.submodel = BusModel(model, self)
         return self.submodel
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.excess_penalty_per_flow_hour = flow_system.fit_to_model_coords(
+        self.excess_penalty_per_flow_hour = self.flow_system.fit_to_model_coords(
             f'{prefix}|excess_penalty_per_flow_hour', self.excess_penalty_per_flow_hour
         )
 
@@ -468,35 +468,39 @@ def create_model(self, model: FlowSystemModel) -> FlowModel:
         self.submodel = FlowModel(model, self)
         return self.submodel
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.relative_minimum = flow_system.fit_to_model_coords(f'{prefix}|relative_minimum', self.relative_minimum)
-        self.relative_maximum = flow_system.fit_to_model_coords(f'{prefix}|relative_maximum', self.relative_maximum)
-        self.fixed_relative_profile = flow_system.fit_to_model_coords(
+        self.relative_minimum = self.flow_system.fit_to_model_coords(
+            f'{prefix}|relative_minimum', self.relative_minimum
+        )
+        self.relative_maximum = self.flow_system.fit_to_model_coords(
+            f'{prefix}|relative_maximum', self.relative_maximum
+        )
+        self.fixed_relative_profile = self.flow_system.fit_to_model_coords(
             f'{prefix}|fixed_relative_profile', self.fixed_relative_profile
         )
-        self.effects_per_flow_hour = flow_system.fit_effects_to_model_coords(
+        self.effects_per_flow_hour = self.flow_system.fit_effects_to_model_coords(
             prefix, self.effects_per_flow_hour, 'per_flow_hour'
         )
-        self.flow_hours_total_max = flow_system.fit_to_model_coords(
+        self.flow_hours_total_max = self.flow_system.fit_to_model_coords(
             f'{prefix}|flow_hours_total_max', self.flow_hours_total_max, dims=['period', 'scenario']
         )
-        self.flow_hours_total_min = flow_system.fit_to_model_coords(
+        self.flow_hours_total_min = self.flow_system.fit_to_model_coords(
             f'{prefix}|flow_hours_total_min', self.flow_hours_total_min, dims=['period', 'scenario']
         )
-        self.load_factor_max = flow_system.fit_to_model_coords(
+        self.load_factor_max = self.flow_system.fit_to_model_coords(
             f'{prefix}|load_factor_max', self.load_factor_max, dims=['period', 'scenario']
         )
-        self.load_factor_min = flow_system.fit_to_model_coords(
+        self.load_factor_min = self.flow_system.fit_to_model_coords(
             f'{prefix}|load_factor_min', self.load_factor_min, dims=['period', 'scenario']
         )
 
         if self.on_off_parameters is not None:
-            self.on_off_parameters.transform_data(flow_system, prefix)
+            self.on_off_parameters.transform_data(prefix)
         if isinstance(self.size, InvestParameters):
-            self.size.transform_data(flow_system, prefix)
+            self.size.transform_data(prefix)
         else:
-            self.size = flow_system.fit_to_model_coords(f'{prefix}|size', self.size, dims=['period', 'scenario'])
+            self.size = self.flow_system.fit_to_model_coords(f'{prefix}|size', self.size, dims=['period', 'scenario'])
 
     def _plausibility_checks(self) -> None:
         # TODO: Incorporate into Variable? (Lower_bound can not be greater than upper bound
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 1fc280226..af50a758d 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -604,7 +604,7 @@ def connect_and_transform(self):
 
         self._connect_network()
         for element in chain(self.components.values(), self.effects.values(), self.buses.values()):
-            element.transform_data(self)
+            element.transform_data()
         self._connected_and_transformed = True
 
     def add_elements(self, *elements: Element) -> None:
@@ -644,6 +644,8 @@ def create_model(self, normalize_weights: bool = True) -> FlowSystemModel:
             raise RuntimeError(
                 'FlowSystem is not connected_and_transformed. Call FlowSystem.connect_and_transform() first.'
             )
+        # Validate cross-element references before creating model
+        self._validate_system_integrity()
         self.model = FlowSystemModel(self, normalize_weights)
         return self.model
 
@@ -777,13 +779,39 @@ def _check_if_element_is_unique(self, element: Element) -> None:
         if element.label_full in self:
             raise ValueError(f'Label of Element {element.label_full} already used in another element!')
 
+    def _validate_system_integrity(self) -> None:
+        """
+        Validate cross-element references to ensure system consistency.
+
+        This performs system-level validation that requires knowledge of multiple elements:
+        - Validates that all Flow.bus references point to existing buses
+        - Can be extended for other cross-element validations
+
+        Should be called after connect_and_transform and before create_model.
+
+        Raises:
+            ValueError: If any cross-element reference is invalid
+        """
+        # Validate bus references in flows
+        for flow in self.flows.values():
+            if flow.bus not in self.buses:
+                available_buses = list(self.buses.keys())
+                raise ValueError(
+                    f'Flow "{flow.label_full}" references bus "{flow.bus}" which does not exist in FlowSystem. '
+                    f'Available buses: {available_buses}. '
+                    f'Did you forget to add the bus using flow_system.add_elements(Bus("{flow.bus}"))?'
+                )
+
     def _add_effects(self, *args: Effect) -> None:
+        for effect in args:
+            effect._set_flow_system(self)  # Link element to FlowSystem
         self.effects.add_effects(*args)
 
     def _add_components(self, *components: Component) -> None:
         for new_component in list(components):
             logger.info(f'Registered new Component: {new_component.label_full}')
             self._check_if_element_is_unique(new_component)  # check if already exists:
+            new_component._set_flow_system(self)  # Link element to FlowSystem
             self.components.add(new_component)  # Add to existing components
             self._flows_cache = None  # Invalidate flows cache
 
@@ -791,6 +819,7 @@ def _add_buses(self, *buses: Bus):
         for new_bus in list(buses):
             logger.info(f'Registered new Bus: {new_bus.label_full}')
             self._check_if_element_is_unique(new_bus)  # check if already exists:
+            new_bus._set_flow_system(self)  # Link element to FlowSystem
             self.buses.add(new_bus)  # Add to existing buses
             self._flows_cache = None  # Invalidate flows cache
 
diff --git a/flixopt/interface.py b/flixopt/interface.py
index 21cbc82b9..e1724715d 100644
--- a/flixopt/interface.py
+++ b/flixopt/interface.py
@@ -78,10 +78,10 @@ def __init__(self, start: TemporalDataUser, end: TemporalDataUser):
         self.end = end
         self.has_time_dim = False
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         dims = None if self.has_time_dim else ['period', 'scenario']
-        self.start = flow_system.fit_to_model_coords(f'{name_prefix}|start', self.start, dims=dims)
-        self.end = flow_system.fit_to_model_coords(f'{name_prefix}|end', self.end, dims=dims)
+        self.start = self._flow_system.fit_to_model_coords(f'{name_prefix}|start', self.start, dims=dims)
+        self.end = self._flow_system.fit_to_model_coords(f'{name_prefix}|end', self.end, dims=dims)
 
 
 @register_class_for_io
@@ -224,9 +224,9 @@ def __getitem__(self, index) -> Piece:
     def __iter__(self) -> Iterator[Piece]:
         return iter(self.pieces)  # Enables iteration like for piece in piecewise: ...
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         for i, piece in enumerate(self.pieces):
-            piece.transform_data(flow_system, f'{name_prefix}|Piece{i}')
+            piece.transform_data(f'{name_prefix}|Piece{i}')
 
 
 @register_class_for_io
@@ -450,9 +450,9 @@ def items(self):
         """
         return self.piecewises.items()
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         for name, piecewise in self.piecewises.items():
-            piecewise.transform_data(flow_system, f'{name_prefix}|{name}')
+            piecewise.transform_data(f'{name_prefix}|{name}')
 
 
 @register_class_for_io
@@ -662,10 +662,10 @@ def has_time_dim(self, value):
         for piecewise in self.piecewise_shares.values():
             piecewise.has_time_dim = value
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
-        self.piecewise_origin.transform_data(flow_system, f'{name_prefix}|PiecewiseEffects|origin')
+    def transform_data(self, name_prefix: str = '') -> None:
+        self.piecewise_origin.transform_data(f'{name_prefix}|PiecewiseEffects|origin')
         for effect, piecewise in self.piecewise_shares.items():
-            piecewise.transform_data(flow_system, f'{name_prefix}|PiecewiseEffects|{effect}')
+            piecewise.transform_data(f'{name_prefix}|PiecewiseEffects|{effect}')
 
 
 @register_class_for_io
@@ -928,20 +928,20 @@ def __init__(
         self.maximum_size = maximum_size if maximum_size is not None else CONFIG.Modeling.big  # default maximum
         self.linked_periods = linked_periods
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
-        self.effects_of_investment = flow_system.fit_effects_to_model_coords(
+    def transform_data(self, name_prefix: str = '') -> None:
+        self.effects_of_investment = self._flow_system.fit_effects_to_model_coords(
             label_prefix=name_prefix,
             effect_values=self.effects_of_investment,
             label_suffix='effects_of_investment',
             dims=['period', 'scenario'],
         )
-        self.effects_of_retirement = flow_system.fit_effects_to_model_coords(
+        self.effects_of_retirement = self._flow_system.fit_effects_to_model_coords(
             label_prefix=name_prefix,
             effect_values=self.effects_of_retirement,
             label_suffix='effects_of_retirement',
             dims=['period', 'scenario'],
         )
-        self.effects_of_investment_per_size = flow_system.fit_effects_to_model_coords(
+        self.effects_of_investment_per_size = self._flow_system.fit_effects_to_model_coords(
             label_prefix=name_prefix,
             effect_values=self.effects_of_investment_per_size,
             label_suffix='effects_of_investment_per_size',
@@ -950,12 +950,12 @@ def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None
 
         if self.piecewise_effects_of_investment is not None:
             self.piecewise_effects_of_investment.has_time_dim = False
-            self.piecewise_effects_of_investment.transform_data(flow_system, f'{name_prefix}|PiecewiseEffects')
+            self.piecewise_effects_of_investment.transform_data(f'{name_prefix}|PiecewiseEffects')
 
-        self.minimum_size = flow_system.fit_to_model_coords(
+        self.minimum_size = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|minimum_size', self.minimum_size, dims=['period', 'scenario']
         )
-        self.maximum_size = flow_system.fit_to_model_coords(
+        self.maximum_size = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|maximum_size', self.maximum_size, dims=['period', 'scenario']
         )
         # Convert tuple (first_period, last_period) to DataArray if needed
@@ -964,28 +964,28 @@ def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None
                 raise TypeError(
                     f'If you provide a tuple to "linked_periods", it needs to be len=2. Got {len(self.linked_periods)=}'
                 )
-            if flow_system.periods is None:
+            if self._flow_system.periods is None:
                 raise ValueError(
                     f'Cannot use linked_periods={self.linked_periods} when FlowSystem has no periods defined. '
                     f'Please define periods in FlowSystem or use linked_periods=None.'
                 )
             logger.debug(f'Computing linked_periods from {self.linked_periods}')
             start, end = self.linked_periods
-            if start not in flow_system.periods.values:
+            if start not in self._flow_system.periods.values:
                 logger.warning(
-                    f'Start of linked periods ({start} not found in periods directly: {flow_system.periods.values}'
+                    f'Start of linked periods ({start} not found in periods directly: {self._flow_system.periods.values}'
                 )
-            if end not in flow_system.periods.values:
+            if end not in self._flow_system.periods.values:
                 logger.warning(
-                    f'End of linked periods ({end} not found in periods directly: {flow_system.periods.values}'
+                    f'End of linked periods ({end} not found in periods directly: {self._flow_system.periods.values}'
                 )
-            self.linked_periods = self.compute_linked_periods(start, end, flow_system.periods)
+            self.linked_periods = self.compute_linked_periods(start, end, self._flow_system.periods)
             logger.debug(f'Computed {self.linked_periods=}')
 
-        self.linked_periods = flow_system.fit_to_model_coords(
+        self.linked_periods = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|linked_periods', self.linked_periods, dims=['period', 'scenario']
         )
-        self.fixed_size = flow_system.fit_to_model_coords(
+        self.fixed_size = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|fixed_size', self.fixed_size, dims=['period', 'scenario']
         )
 
@@ -1294,32 +1294,32 @@ def __init__(
         self.switch_on_total_max: Scalar = switch_on_total_max
         self.force_switch_on: bool = force_switch_on
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
-        self.effects_per_switch_on = flow_system.fit_effects_to_model_coords(
+    def transform_data(self, name_prefix: str = '') -> None:
+        self.effects_per_switch_on = self._flow_system.fit_effects_to_model_coords(
             name_prefix, self.effects_per_switch_on, 'per_switch_on'
         )
-        self.effects_per_running_hour = flow_system.fit_effects_to_model_coords(
+        self.effects_per_running_hour = self._flow_system.fit_effects_to_model_coords(
             name_prefix, self.effects_per_running_hour, 'per_running_hour'
         )
-        self.consecutive_on_hours_min = flow_system.fit_to_model_coords(
+        self.consecutive_on_hours_min = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|consecutive_on_hours_min', self.consecutive_on_hours_min
         )
-        self.consecutive_on_hours_max = flow_system.fit_to_model_coords(
+        self.consecutive_on_hours_max = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|consecutive_on_hours_max', self.consecutive_on_hours_max
         )
-        self.consecutive_off_hours_min = flow_system.fit_to_model_coords(
+        self.consecutive_off_hours_min = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|consecutive_off_hours_min', self.consecutive_off_hours_min
         )
-        self.consecutive_off_hours_max = flow_system.fit_to_model_coords(
+        self.consecutive_off_hours_max = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|consecutive_off_hours_max', self.consecutive_off_hours_max
         )
-        self.on_hours_total_max = flow_system.fit_to_model_coords(
+        self.on_hours_total_max = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|on_hours_total_max', self.on_hours_total_max, dims=['period', 'scenario']
         )
-        self.on_hours_total_min = flow_system.fit_to_model_coords(
+        self.on_hours_total_min = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|on_hours_total_min', self.on_hours_total_min, dims=['period', 'scenario']
         )
-        self.switch_on_total_max = flow_system.fit_to_model_coords(
+        self.switch_on_total_max = self._flow_system.fit_to_model_coords(
             f'{name_prefix}|switch_on_total_max', self.switch_on_total_max, dims=['period', 'scenario']
         )
 
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 2bce6aa52..09ad5291e 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -272,18 +272,51 @@ class Interface:
         transform_data(flow_system): Transform data to match FlowSystem dimensions
     """
 
-    def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None:
+    def transform_data(self, name_prefix: str = '') -> None:
         """Transform the data of the interface to match the FlowSystem's dimensions.
 
         Args:
-            flow_system: The FlowSystem containing timing and dimensional information
             name_prefix: The prefix to use for the names of the variables. Defaults to '', which results in no prefix.
 
         Raises:
             NotImplementedError: Must be implemented by subclasses
+
+        Note:
+            The FlowSystem reference is available via self._flow_system (for Interface objects)
+            or self.flow_system property (for Element objects). Elements must be registered
+            to a FlowSystem before calling this method.
         """
         raise NotImplementedError('Every Interface subclass needs a transform_data() method')
 
+    def _set_flow_system(self, flow_system: FlowSystem) -> None:
+        """Store flow_system reference and propagate to nested Interface objects.
+
+        This method is called automatically during element registration to enable
+        elements to access FlowSystem properties without passing the reference
+        through every method call.
+
+        Args:
+            flow_system: The FlowSystem that this interface belongs to
+        """
+        # Always set _flow_system (creates attribute if it doesn't exist)
+        self._flow_system = flow_system
+
+        # Recursively set for nested Interface objects
+        for attr_name, attr_value in self.__dict__.items():
+            if attr_name.startswith('_'):
+                continue  # Skip private attributes
+
+            if isinstance(attr_value, Interface):
+                attr_value._set_flow_system(flow_system)
+            elif isinstance(attr_value, list):
+                for item in attr_value:
+                    if isinstance(item, Interface):
+                        item._set_flow_system(flow_system)
+            elif isinstance(attr_value, dict):
+                for item in attr_value.values():
+                    if isinstance(item, Interface):
+                        item._set_flow_system(flow_system)
+
     def _create_reference_structure(self) -> tuple[dict, dict[str, xr.DataArray]]:
         """
         Convert all DataArrays to references and extract them.
@@ -861,6 +894,24 @@ def __init__(self, label: str, meta_data: dict | None = None):
         self.label = Element._valid_label(label)
         self.meta_data = meta_data if meta_data is not None else {}
         self.submodel = None
+        self._flow_system: FlowSystem | None = None
+
+    @property
+    def flow_system(self) -> FlowSystem:
+        """Access the FlowSystem this element is registered to.
+
+        Returns:
+            The FlowSystem instance this element belongs to.
+
+        Raises:
+            RuntimeError: If element has not been registered to a FlowSystem yet.
+        """
+        if self._flow_system is None:
+            raise RuntimeError(
+                f'Element "{self.label_full}" is not registered to a FlowSystem. '
+                f'Call flow_system.add_elements() to register this element first.'
+            )
+        return self._flow_system
 
     def _plausibility_checks(self) -> None:
         """This function is used to do some basic plausibility checks for each Element during initialization.

From 76a65a7d398b23820a04b5db23593ef9782ca2f1 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 10:53:20 +0100
Subject: [PATCH 02/20]  implemented the helper methods in the Interface base
 class

---
 flixopt/components.py | 38 +++++++++++--------------
 flixopt/effects.py    | 37 ++++++++++--------------
 flixopt/elements.py   | 28 +++++++------------
 flixopt/interface.py  | 64 ++++++++++++++++++++----------------------
 flixopt/structure.py  | 65 ++++++++++++++++++++++++++++++++-----------
 5 files changed, 119 insertions(+), 113 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 54b5e4ea1..2046e2916 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -228,9 +228,7 @@ def _transform_conversion_factors(self) -> list[dict[str, xr.DataArray]]:
             transformed_dict = {}
             for flow, values in conversion_factor.items():
                 # TODO: Might be better to use the label of the component instead of the flow
-                ts = self.flow_system.fit_to_model_coords(
-                    f'{self.flows[flow].label_full}|conversion_factor{idx}', values
-                )
+                ts = self._fit_coords(f'{self.flows[flow].label_full}|conversion_factor{idx}', values)
                 if ts is None:
                     raise PlausibilityError(f'{self.label_full}: conversion factor for flow "{flow}" must not be None')
                 transformed_dict[flow] = ts
@@ -439,35 +437,31 @@ def create_model(self, model: FlowSystemModel) -> StorageModel:
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         super().transform_data(prefix)
-        self.relative_minimum_charge_state = self.flow_system.fit_to_model_coords(
-            f'{prefix}|relative_minimum_charge_state',
-            self.relative_minimum_charge_state,
-        )
-        self.relative_maximum_charge_state = self.flow_system.fit_to_model_coords(
-            f'{prefix}|relative_maximum_charge_state',
-            self.relative_maximum_charge_state,
+        self.relative_minimum_charge_state = self._fit_coords(
+            f'{prefix}|relative_minimum_charge_state', self.relative_minimum_charge_state
         )
-        self.eta_charge = self.flow_system.fit_to_model_coords(f'{prefix}|eta_charge', self.eta_charge)
-        self.eta_discharge = self.flow_system.fit_to_model_coords(f'{prefix}|eta_discharge', self.eta_discharge)
-        self.relative_loss_per_hour = self.flow_system.fit_to_model_coords(
-            f'{prefix}|relative_loss_per_hour', self.relative_loss_per_hour
+        self.relative_maximum_charge_state = self._fit_coords(
+            f'{prefix}|relative_maximum_charge_state', self.relative_maximum_charge_state
         )
+        self.eta_charge = self._fit_coords(f'{prefix}|eta_charge', self.eta_charge)
+        self.eta_discharge = self._fit_coords(f'{prefix}|eta_discharge', self.eta_discharge)
+        self.relative_loss_per_hour = self._fit_coords(f'{prefix}|relative_loss_per_hour', self.relative_loss_per_hour)
         if not isinstance(self.initial_charge_state, str):
-            self.initial_charge_state = self.flow_system.fit_to_model_coords(
+            self.initial_charge_state = self._fit_coords(
                 f'{prefix}|initial_charge_state', self.initial_charge_state, dims=['period', 'scenario']
             )
-        self.minimal_final_charge_state = self.flow_system.fit_to_model_coords(
+        self.minimal_final_charge_state = self._fit_coords(
             f'{prefix}|minimal_final_charge_state', self.minimal_final_charge_state, dims=['period', 'scenario']
         )
-        self.maximal_final_charge_state = self.flow_system.fit_to_model_coords(
+        self.maximal_final_charge_state = self._fit_coords(
             f'{prefix}|maximal_final_charge_state', self.maximal_final_charge_state, dims=['period', 'scenario']
         )
-        self.relative_minimum_final_charge_state = self.flow_system.fit_to_model_coords(
+        self.relative_minimum_final_charge_state = self._fit_coords(
             f'{prefix}|relative_minimum_final_charge_state',
             self.relative_minimum_final_charge_state,
             dims=['period', 'scenario'],
         )
-        self.relative_maximum_final_charge_state = self.flow_system.fit_to_model_coords(
+        self.relative_maximum_final_charge_state = self._fit_coords(
             f'{prefix}|relative_maximum_final_charge_state',
             self.relative_maximum_final_charge_state,
             dims=['period', 'scenario'],
@@ -475,7 +469,7 @@ def transform_data(self, name_prefix: str = '') -> None:
         if isinstance(self.capacity_in_flow_hours, InvestParameters):
             self.capacity_in_flow_hours.transform_data(f'{prefix}|InvestParameters')
         else:
-            self.capacity_in_flow_hours = self.flow_system.fit_to_model_coords(
+            self.capacity_in_flow_hours = self._fit_coords(
                 f'{prefix}|capacity_in_flow_hours', self.capacity_in_flow_hours, dims=['period', 'scenario']
             )
 
@@ -725,8 +719,8 @@ def create_model(self, model) -> TransmissionModel:
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         super().transform_data(prefix)
-        self.relative_losses = self.flow_system.fit_to_model_coords(f'{prefix}|relative_losses', self.relative_losses)
-        self.absolute_losses = self.flow_system.fit_to_model_coords(f'{prefix}|absolute_losses', self.absolute_losses)
+        self.relative_losses = self._fit_coords(f'{prefix}|relative_losses', self.relative_losses)
+        self.absolute_losses = self._fit_coords(f'{prefix}|absolute_losses', self.absolute_losses)
 
 
 class TransmissionModel(ComponentModel):
diff --git a/flixopt/effects.py b/flixopt/effects.py
index 4b0a6b1d1..b55792469 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -342,45 +342,38 @@ def maximum_operation_per_hour(self, value):
 
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.minimum_per_hour = self.flow_system.fit_to_model_coords(
-            f'{prefix}|minimum_per_hour', self.minimum_per_hour
-        )
-
-        self.maximum_per_hour = self.flow_system.fit_to_model_coords(
-            f'{prefix}|maximum_per_hour', self.maximum_per_hour
-        )
+        self.minimum_per_hour = self._fit_coords(f'{prefix}|minimum_per_hour', self.minimum_per_hour)
+        self.maximum_per_hour = self._fit_coords(f'{prefix}|maximum_per_hour', self.maximum_per_hour)
 
-        self.share_from_temporal = self.flow_system.fit_effects_to_model_coords(
-            label_prefix=None,
+        self.share_from_temporal = self._fit_effect_coords(
+            prefix=None,
             effect_values=self.share_from_temporal,
-            label_suffix=f'(temporal)->{prefix}(temporal)',
+            suffix=f'(temporal)->{prefix}(temporal)',
             dims=['time', 'period', 'scenario'],
         )
-        self.share_from_periodic = self.flow_system.fit_effects_to_model_coords(
-            label_prefix=None,
+        self.share_from_periodic = self._fit_effect_coords(
+            prefix=None,
             effect_values=self.share_from_periodic,
-            label_suffix=f'(periodic)->{prefix}(periodic)',
+            suffix=f'(periodic)->{prefix}(periodic)',
             dims=['period', 'scenario'],
         )
 
-        self.minimum_temporal = self.flow_system.fit_to_model_coords(
+        self.minimum_temporal = self._fit_coords(
             f'{prefix}|minimum_temporal', self.minimum_temporal, dims=['period', 'scenario']
         )
-        self.maximum_temporal = self.flow_system.fit_to_model_coords(
+        self.maximum_temporal = self._fit_coords(
             f'{prefix}|maximum_temporal', self.maximum_temporal, dims=['period', 'scenario']
         )
-        self.minimum_periodic = self.flow_system.fit_to_model_coords(
+        self.minimum_periodic = self._fit_coords(
             f'{prefix}|minimum_periodic', self.minimum_periodic, dims=['period', 'scenario']
         )
-        self.maximum_periodic = self.flow_system.fit_to_model_coords(
+        self.maximum_periodic = self._fit_coords(
             f'{prefix}|maximum_periodic', self.maximum_periodic, dims=['period', 'scenario']
         )
-        self.minimum_total = self.flow_system.fit_to_model_coords(
-            f'{prefix}|minimum_total',
-            self.minimum_total,
-            dims=['period', 'scenario'],
+        self.minimum_total = self._fit_coords(
+            f'{prefix}|minimum_total', self.minimum_total, dims=['period', 'scenario']
         )
-        self.maximum_total = self.flow_system.fit_to_model_coords(
+        self.maximum_total = self._fit_coords(
             f'{prefix}|maximum_total', self.maximum_total, dims=['period', 'scenario']
         )
 
diff --git a/flixopt/elements.py b/flixopt/elements.py
index aaa3e6cb4..64e0113bf 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -243,7 +243,7 @@ def create_model(self, model: FlowSystemModel) -> BusModel:
 
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.excess_penalty_per_flow_hour = self.flow_system.fit_to_model_coords(
+        self.excess_penalty_per_flow_hour = self._fit_coords(
             f'{prefix}|excess_penalty_per_flow_hour', self.excess_penalty_per_flow_hour
         )
 
@@ -470,28 +470,20 @@ def create_model(self, model: FlowSystemModel) -> FlowModel:
 
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
-        self.relative_minimum = self.flow_system.fit_to_model_coords(
-            f'{prefix}|relative_minimum', self.relative_minimum
-        )
-        self.relative_maximum = self.flow_system.fit_to_model_coords(
-            f'{prefix}|relative_maximum', self.relative_maximum
-        )
-        self.fixed_relative_profile = self.flow_system.fit_to_model_coords(
-            f'{prefix}|fixed_relative_profile', self.fixed_relative_profile
-        )
-        self.effects_per_flow_hour = self.flow_system.fit_effects_to_model_coords(
-            prefix, self.effects_per_flow_hour, 'per_flow_hour'
-        )
-        self.flow_hours_total_max = self.flow_system.fit_to_model_coords(
+        self.relative_minimum = self._fit_coords(f'{prefix}|relative_minimum', self.relative_minimum)
+        self.relative_maximum = self._fit_coords(f'{prefix}|relative_maximum', self.relative_maximum)
+        self.fixed_relative_profile = self._fit_coords(f'{prefix}|fixed_relative_profile', self.fixed_relative_profile)
+        self.effects_per_flow_hour = self._fit_effect_coords(prefix, self.effects_per_flow_hour, 'per_flow_hour')
+        self.flow_hours_total_max = self._fit_coords(
             f'{prefix}|flow_hours_total_max', self.flow_hours_total_max, dims=['period', 'scenario']
         )
-        self.flow_hours_total_min = self.flow_system.fit_to_model_coords(
+        self.flow_hours_total_min = self._fit_coords(
             f'{prefix}|flow_hours_total_min', self.flow_hours_total_min, dims=['period', 'scenario']
         )
-        self.load_factor_max = self.flow_system.fit_to_model_coords(
+        self.load_factor_max = self._fit_coords(
             f'{prefix}|load_factor_max', self.load_factor_max, dims=['period', 'scenario']
         )
-        self.load_factor_min = self.flow_system.fit_to_model_coords(
+        self.load_factor_min = self._fit_coords(
             f'{prefix}|load_factor_min', self.load_factor_min, dims=['period', 'scenario']
         )
 
@@ -500,7 +492,7 @@ def transform_data(self, name_prefix: str = '') -> None:
         if isinstance(self.size, InvestParameters):
             self.size.transform_data(prefix)
         else:
-            self.size = self.flow_system.fit_to_model_coords(f'{prefix}|size', self.size, dims=['period', 'scenario'])
+            self.size = self._fit_coords(f'{prefix}|size', self.size, dims=['period', 'scenario'])
 
     def _plausibility_checks(self) -> None:
         # TODO: Incorporate into Variable? (Lower_bound can not be greater than upper bound
diff --git a/flixopt/interface.py b/flixopt/interface.py
index e1724715d..3f21c027b 100644
--- a/flixopt/interface.py
+++ b/flixopt/interface.py
@@ -80,8 +80,8 @@ def __init__(self, start: TemporalDataUser, end: TemporalDataUser):
 
     def transform_data(self, name_prefix: str = '') -> None:
         dims = None if self.has_time_dim else ['period', 'scenario']
-        self.start = self._flow_system.fit_to_model_coords(f'{name_prefix}|start', self.start, dims=dims)
-        self.end = self._flow_system.fit_to_model_coords(f'{name_prefix}|end', self.end, dims=dims)
+        self.start = self._fit_coords(f'{name_prefix}|start', self.start, dims=dims)
+        self.end = self._fit_coords(f'{name_prefix}|end', self.end, dims=dims)
 
 
 @register_class_for_io
@@ -929,22 +929,22 @@ def __init__(
         self.linked_periods = linked_periods
 
     def transform_data(self, name_prefix: str = '') -> None:
-        self.effects_of_investment = self._flow_system.fit_effects_to_model_coords(
-            label_prefix=name_prefix,
+        self.effects_of_investment = self._fit_effect_coords(
+            prefix=name_prefix,
             effect_values=self.effects_of_investment,
-            label_suffix='effects_of_investment',
+            suffix='effects_of_investment',
             dims=['period', 'scenario'],
         )
-        self.effects_of_retirement = self._flow_system.fit_effects_to_model_coords(
-            label_prefix=name_prefix,
+        self.effects_of_retirement = self._fit_effect_coords(
+            prefix=name_prefix,
             effect_values=self.effects_of_retirement,
-            label_suffix='effects_of_retirement',
+            suffix='effects_of_retirement',
             dims=['period', 'scenario'],
         )
-        self.effects_of_investment_per_size = self._flow_system.fit_effects_to_model_coords(
-            label_prefix=name_prefix,
+        self.effects_of_investment_per_size = self._fit_effect_coords(
+            prefix=name_prefix,
             effect_values=self.effects_of_investment_per_size,
-            label_suffix='effects_of_investment_per_size',
+            suffix='effects_of_investment_per_size',
             dims=['period', 'scenario'],
         )
 
@@ -952,10 +952,10 @@ def transform_data(self, name_prefix: str = '') -> None:
             self.piecewise_effects_of_investment.has_time_dim = False
             self.piecewise_effects_of_investment.transform_data(f'{name_prefix}|PiecewiseEffects')
 
-        self.minimum_size = self._flow_system.fit_to_model_coords(
+        self.minimum_size = self._fit_coords(
             f'{name_prefix}|minimum_size', self.minimum_size, dims=['period', 'scenario']
         )
-        self.maximum_size = self._flow_system.fit_to_model_coords(
+        self.maximum_size = self._fit_coords(
             f'{name_prefix}|maximum_size', self.maximum_size, dims=['period', 'scenario']
         )
         # Convert tuple (first_period, last_period) to DataArray if needed
@@ -964,30 +964,28 @@ def transform_data(self, name_prefix: str = '') -> None:
                 raise TypeError(
                     f'If you provide a tuple to "linked_periods", it needs to be len=2. Got {len(self.linked_periods)=}'
                 )
-            if self._flow_system.periods is None:
+            if self.flow_system.periods is None:
                 raise ValueError(
                     f'Cannot use linked_periods={self.linked_periods} when FlowSystem has no periods defined. '
                     f'Please define periods in FlowSystem or use linked_periods=None.'
                 )
             logger.debug(f'Computing linked_periods from {self.linked_periods}')
             start, end = self.linked_periods
-            if start not in self._flow_system.periods.values:
+            if start not in self.flow_system.periods.values:
                 logger.warning(
-                    f'Start of linked periods ({start} not found in periods directly: {self._flow_system.periods.values}'
+                    f'Start of linked periods ({start} not found in periods directly: {self.flow_system.periods.values}'
                 )
-            if end not in self._flow_system.periods.values:
+            if end not in self.flow_system.periods.values:
                 logger.warning(
-                    f'End of linked periods ({end} not found in periods directly: {self._flow_system.periods.values}'
+                    f'End of linked periods ({end} not found in periods directly: {self.flow_system.periods.values}'
                 )
-            self.linked_periods = self.compute_linked_periods(start, end, self._flow_system.periods)
+            self.linked_periods = self.compute_linked_periods(start, end, self.flow_system.periods)
             logger.debug(f'Computed {self.linked_periods=}')
 
-        self.linked_periods = self._flow_system.fit_to_model_coords(
+        self.linked_periods = self._fit_coords(
             f'{name_prefix}|linked_periods', self.linked_periods, dims=['period', 'scenario']
         )
-        self.fixed_size = self._flow_system.fit_to_model_coords(
-            f'{name_prefix}|fixed_size', self.fixed_size, dims=['period', 'scenario']
-        )
+        self.fixed_size = self._fit_coords(f'{name_prefix}|fixed_size', self.fixed_size, dims=['period', 'scenario'])
 
     @property
     def optional(self) -> bool:
@@ -1295,31 +1293,29 @@ def __init__(
         self.force_switch_on: bool = force_switch_on
 
     def transform_data(self, name_prefix: str = '') -> None:
-        self.effects_per_switch_on = self._flow_system.fit_effects_to_model_coords(
-            name_prefix, self.effects_per_switch_on, 'per_switch_on'
-        )
-        self.effects_per_running_hour = self._flow_system.fit_effects_to_model_coords(
+        self.effects_per_switch_on = self._fit_effect_coords(name_prefix, self.effects_per_switch_on, 'per_switch_on')
+        self.effects_per_running_hour = self._fit_effect_coords(
             name_prefix, self.effects_per_running_hour, 'per_running_hour'
         )
-        self.consecutive_on_hours_min = self._flow_system.fit_to_model_coords(
+        self.consecutive_on_hours_min = self._fit_coords(
             f'{name_prefix}|consecutive_on_hours_min', self.consecutive_on_hours_min
         )
-        self.consecutive_on_hours_max = self._flow_system.fit_to_model_coords(
+        self.consecutive_on_hours_max = self._fit_coords(
             f'{name_prefix}|consecutive_on_hours_max', self.consecutive_on_hours_max
         )
-        self.consecutive_off_hours_min = self._flow_system.fit_to_model_coords(
+        self.consecutive_off_hours_min = self._fit_coords(
             f'{name_prefix}|consecutive_off_hours_min', self.consecutive_off_hours_min
         )
-        self.consecutive_off_hours_max = self._flow_system.fit_to_model_coords(
+        self.consecutive_off_hours_max = self._fit_coords(
             f'{name_prefix}|consecutive_off_hours_max', self.consecutive_off_hours_max
         )
-        self.on_hours_total_max = self._flow_system.fit_to_model_coords(
+        self.on_hours_total_max = self._fit_coords(
             f'{name_prefix}|on_hours_total_max', self.on_hours_total_max, dims=['period', 'scenario']
         )
-        self.on_hours_total_min = self._flow_system.fit_to_model_coords(
+        self.on_hours_total_min = self._fit_coords(
             f'{name_prefix}|on_hours_total_min', self.on_hours_total_min, dims=['period', 'scenario']
         )
-        self.switch_on_total_max = self._flow_system.fit_to_model_coords(
+        self.switch_on_total_max = self._fit_coords(
             f'{name_prefix}|switch_on_total_max', self.switch_on_total_max, dims=['period', 'scenario']
         )
 
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 09ad5291e..b6e9e91a5 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -317,6 +317,54 @@ def _set_flow_system(self, flow_system: FlowSystem) -> None:
                     if isinstance(item, Interface):
                         item._set_flow_system(flow_system)
 
+    @property
+    def flow_system(self) -> FlowSystem:
+        """Access the FlowSystem this interface is linked to.
+
+        Returns:
+            The FlowSystem instance this interface belongs to.
+
+        Raises:
+            RuntimeError: If interface has not been linked to a FlowSystem yet.
+
+        Note:
+            For Elements, this is set during add_elements().
+            For parameter classes, this is set recursively when the parent Element is registered.
+        """
+        if not hasattr(self, '_flow_system') or self._flow_system is None:
+            raise RuntimeError(
+                f'{self.__class__.__name__} is not linked to a FlowSystem. '
+                f'Ensure the parent element is registered via flow_system.add_elements() first.'
+            )
+        return self._flow_system
+
+    def _fit_coords(self, name: str, data, dims=None):
+        """Convenience wrapper for FlowSystem.fit_to_model_coords().
+
+        Args:
+            name: The name for the data variable
+            data: The data to transform
+            dims: Optional dimension names
+
+        Returns:
+            Transformed data aligned to FlowSystem coordinates
+        """
+        return self.flow_system.fit_to_model_coords(name, data, dims=dims)
+
+    def _fit_effect_coords(self, prefix: str, effect_values, suffix: str = None, dims=None):
+        """Convenience wrapper for FlowSystem.fit_effects_to_model_coords().
+
+        Args:
+            prefix: Label prefix for effect names
+            effect_values: The effect values to transform
+            suffix: Optional label suffix
+            dims: Optional dimension names
+
+        Returns:
+            Transformed effect values aligned to FlowSystem coordinates
+        """
+        return self.flow_system.fit_effects_to_model_coords(prefix, effect_values, suffix, dims=dims)
+
     def _create_reference_structure(self) -> tuple[dict, dict[str, xr.DataArray]]:
         """
         Convert all DataArrays to references and extract them.
@@ -896,23 +944,6 @@ def __init__(self, label: str, meta_data: dict | None = None):
         self.submodel = None
         self._flow_system: FlowSystem | None = None
 
-    @property
-    def flow_system(self) -> FlowSystem:
-        """Access the FlowSystem this element is registered to.
-
-        Returns:
-            The FlowSystem instance this element belongs to.
-
-        Raises:
-            RuntimeError: If element has not been registered to a FlowSystem yet.
-        """
-        if self._flow_system is None:
-            raise RuntimeError(
-                f'Element "{self.label_full}" is not registered to a FlowSystem. '
-                f'Call flow_system.add_elements() to register this element first.'
-            )
-        return self._flow_system
-
     def _plausibility_checks(self) -> None:
         """This function is used to do some basic plausibility checks for each Element during initialization.
         This is run after all data is transformed to the correct format/type"""

From 11fcd56e9603d0e56e5f1810fb7133d81871e97b Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 11:21:00 +0100
Subject: [PATCH 03/20] made calling .do_modeling() optional

---
 flixopt/calculation.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/flixopt/calculation.py b/flixopt/calculation.py
index 5de2c8870..072e9204d 100644
--- a/flixopt/calculation.py
+++ b/flixopt/calculation.py
@@ -198,6 +198,7 @@ def do_modeling(self) -> FullCalculation:
         self.model.do_modeling()
 
         self.durations['modeling'] = round(timeit.default_timer() - t_start, 2)
+        self._modeled = True
         return self
 
     def fix_sizes(self, ds: xr.Dataset, decimal_rounding: int | None = 5) -> FullCalculation:
@@ -230,6 +231,11 @@ def fix_sizes(self, ds: xr.Dataset, decimal_rounding: int | None = 5) -> FullCal
     def solve(
         self, solver: _Solver, log_file: pathlib.Path | None = None, log_main_results: bool | None = None
     ) -> FullCalculation:
+        # Auto-call do_modeling() if not already done
+        if not self._modeled:
+            logger.info('Model not yet created. Calling do_modeling() automatically.')
+            self.do_modeling()
+
         t_start = timeit.default_timer()
 
         self.model.solve(
@@ -328,6 +334,7 @@ def do_modeling(self) -> AggregatedCalculation:
         )
         self.aggregation_model.do_modeling()
         self.durations['modeling'] = round(timeit.default_timer() - t_start, 2)
+        self._modeled = True
         return self
 
     def _perform_aggregation(self):

From f509ec0066c66177486bd959850743a5871b0291 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 12:07:01 +0100
Subject: [PATCH 04/20] eliminated all circular dependencies in the Submodel
 architecture by implementing a two-phase modeling pattern:

  Phase 1 - Variable Creation (_create_variables()):
  - Called during __init__
  - Creates all variables and submodels
  - No constraints allowed (to avoid circular dependencies)

  Phase 2 - Constraint Creation (_create_constraints()):
  - Called from FlowSystemModel.do_modeling() after all models exist
  - Safely accesses variables from child/sibling models
  - Creates all constraints and relationships
---
 flixopt/components.py |  84 +++++++++++++++-----------
 flixopt/effects.py    |  24 +++++++-
 flixopt/elements.py   | 103 ++++++++++++++++++++++----------
 flixopt/features.py   | 133 ++++++++++++++++++++++++++++--------------
 flixopt/structure.py  |  38 +++++++++++-
 5 files changed, 270 insertions(+), 112 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 2046e2916..c215b1cd2 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -734,9 +734,9 @@ def __init__(self, model: FlowSystemModel, element: Transmission):
 
         super().__init__(model, element)
 
-    def _do_modeling(self):
-        """Initiates all FlowModels"""
-        super()._do_modeling()
+    def _create_constraints(self):
+        """Phase 2: Create constraints (FlowModels already created by parent)"""
+        super()._create_constraints()
 
         # first direction
         self.create_transmission_equation('dir1', self.element.in1, self.element.out1)
@@ -775,26 +775,12 @@ def __init__(self, model: FlowSystemModel, element: LinearConverter):
         self.piecewise_conversion: PiecewiseConversion | None = None
         super().__init__(model, element)
 
-    def _do_modeling(self):
-        super()._do_modeling()
-        # conversion_factors:
-        if self.element.conversion_factors:
-            all_input_flows = set(self.element.inputs)
-            all_output_flows = set(self.element.outputs)
-
-            # für alle linearen Gleichungen:
-            for i, conv_factors in enumerate(self.element.conversion_factors):
-                used_flows = set([self.element.flows[flow_label] for flow_label in conv_factors])
-                used_inputs: set[Flow] = all_input_flows & used_flows
-                used_outputs: set[Flow] = all_output_flows & used_flows
-
-                self.add_constraints(
-                    sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_inputs])
-                    == sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_outputs]),
-                    short_name=f'conversion_{i}',
-                )
+    def _create_variables(self):
+        """Phase 1: Create variables (parent creates FlowModels)"""
+        super()._create_variables()
 
-        else:
+        # Create PiecewiseModel if needed (after FlowModels are created)
+        if not self.element.conversion_factors:
             # TODO: Improve Inclusion of OnOffParameters. Instead of creating a Binary in every flow, the binary could only be part of the Piece itself
             piecewise_conversion = {
                 self.element.flows[flow].submodel.flow_rate.name: piecewise
@@ -813,6 +799,27 @@ def _do_modeling(self):
                 short_name='PiecewiseConversion',
             )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Create conversion factor constraints if specified
+        if self.element.conversion_factors:
+            all_input_flows = set(self.element.inputs)
+            all_output_flows = set(self.element.outputs)
+
+            # für alle linearen Gleichungen:
+            for i, conv_factors in enumerate(self.element.conversion_factors):
+                used_flows = set([self.element.flows[flow_label] for flow_label in conv_factors])
+                used_inputs: set[Flow] = all_input_flows & used_flows
+                used_outputs: set[Flow] = all_output_flows & used_flows
+
+                self.add_constraints(
+                    sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_inputs])
+                    == sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_outputs]),
+                    short_name=f'conversion_{i}',
+                )
+
 
 class StorageModel(ComponentModel):
     """Submodel of Storage"""
@@ -822,8 +829,9 @@ class StorageModel(ComponentModel):
     def __init__(self, model: FlowSystemModel, element: Storage):
         super().__init__(model, element)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables (parent creates FlowModels)"""
+        super()._create_variables()
 
         lb, ub = self._absolute_charge_state_bounds
         self.add_variables(
@@ -835,6 +843,22 @@ def _do_modeling(self):
 
         self.add_variables(coords=self._model.get_coords(), short_name='netto_discharge')
 
+        # Create InvestmentModel for capacity if needed
+        if isinstance(self.element.capacity_in_flow_hours, InvestParameters):
+            self.add_submodels(
+                InvestmentModel(
+                    model=self._model,
+                    label_of_element=self.label_of_element,
+                    label_of_model=self.label_of_element,
+                    parameters=self.element.capacity_in_flow_hours,
+                ),
+                short_name='investment',
+            )
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints (can now access flow.submodel.flow_rate)"""
+        super()._create_constraints()
+
         # netto_discharge:
         # eq: nettoFlow(t) - discharging(t) + charging(t) = 0
         self.add_constraints(
@@ -859,17 +883,8 @@ def _do_modeling(self):
             short_name='charge_state',
         )
 
+        # Bounding constraints for investment
         if isinstance(self.element.capacity_in_flow_hours, InvestParameters):
-            self.add_submodels(
-                InvestmentModel(
-                    model=self._model,
-                    label_of_element=self.label_of_element,
-                    label_of_model=self.label_of_element,
-                    parameters=self.element.capacity_in_flow_hours,
-                ),
-                short_name='investment',
-            )
-
             BoundingPatterns.scaled_bounds(
                 self,
                 variable=self.charge_state,
@@ -880,6 +895,7 @@ def _do_modeling(self):
         # Initial charge state
         self._initial_and_final_charge_state()
 
+        # Balanced sizes
         if self.element.balanced:
             self.add_constraints(
                 self.element.charging.submodel._investment.size * 1
diff --git a/flixopt/effects.py b/flixopt/effects.py
index b55792469..312e8244a 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -393,7 +393,10 @@ class EffectModel(ElementModel):
     def __init__(self, model: FlowSystemModel, element: Effect):
         super().__init__(model, element)
 
-    def _do_modeling(self):
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
         self.total: linopy.Variable | None = None
         self.periodic: ShareAllocationModel = self.add_submodels(
             ShareAllocationModel(
@@ -428,6 +431,10 @@ def _do_modeling(self):
             name=self.label_full,
         )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
         self.add_constraints(
             self.total == self.temporal.total + self.periodic.total, name=self.label_full, short_name='total'
         )
@@ -672,17 +679,28 @@ def add_share_to_penalty(self, name: str, expression: linopy.LinearExpression) -
             raise TypeError(f'Penalty shares must be scalar expressions! ({expression.ndim=})')
         self.penalty.add_share(name, expression, dims=())
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
+        # Create EffectModel for each effect
         for effect in self.effects.values():
             effect.create_model(self._model)
+
+        # Create penalty allocation model
         self.penalty = self.add_submodels(
             ShareAllocationModel(self._model, dims=(), label_of_element='Penalty'),
             short_name='penalty',
         )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Add cross-effect shares
         self._add_share_between_effects()
 
+        # Set objective
         self._model.add_objective(
             (self.effects.objective_effect.submodel.total * self._model.weights).sum() + self.penalty.total.sum()
         )
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 64e0113bf..0b821c4b0 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -553,8 +553,10 @@ class FlowModel(ElementModel):
     def __init__(self, model: FlowSystemModel, element: Flow):
         super().__init__(model, element)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
         # Main flow rate variable
         self.add_variables(
             lower=self.absolute_flow_rate_bounds[0],
@@ -563,6 +565,28 @@ def _do_modeling(self):
             short_name='flow_rate',
         )
 
+        # Create investment and/or on_off submodels (creates their variables)
+        if self.with_investment:
+            self._create_investment_model()
+        if self.with_on_off:
+            on = self.add_variables(binary=True, short_name='on', coords=self._model.get_coords())
+            self.add_submodels(
+                OnOffModel(
+                    model=self._model,
+                    label_of_element=self.label_of_element,
+                    parameters=self.element.on_off_parameters,
+                    on_variable=on,
+                    previous_states=self.previous_states,
+                    label_of_model=self.label_of_element,
+                ),
+                short_name='on_off',
+            )
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Create flow rate bounding constraints
         self._constraint_flow_rate()
 
         # Total flow hours tracking
@@ -610,13 +634,13 @@ def _create_investment_model(self):
         )
 
     def _constraint_flow_rate(self):
+        """Create bounding constraints for flow_rate (models already created in _create_variables)"""
         if not self.with_investment and not self.with_on_off:
             # Most basic case. Already covered by direct variable bounds
             pass
 
         elif self.with_on_off and not self.with_investment:
             # OnOff, but no Investment
-            self._create_on_off_model()
             bounds = self.relative_flow_rate_bounds
             BoundingPatterns.bounds_with_state(
                 self,
@@ -627,7 +651,6 @@ def _constraint_flow_rate(self):
 
         elif self.with_investment and not self.with_on_off:
             # Investment, but no OnOff
-            self._create_investment_model()
             BoundingPatterns.scaled_bounds(
                 self,
                 variable=self.flow_rate,
@@ -637,9 +660,6 @@ def _constraint_flow_rate(self):
 
         elif self.with_investment and self.with_on_off:
             # Investment and OnOff
-            self._create_investment_model()
-            self._create_on_off_model()
-
             BoundingPatterns.scaled_bounds_with_state(
                 model=self,
                 variable=self.flow_rate,
@@ -785,25 +805,32 @@ def __init__(self, model: FlowSystemModel, element: Bus):
         self.excess_output: linopy.Variable | None = None
         super().__init__(model, element)
 
-    def _do_modeling(self) -> None:
-        super()._do_modeling()
-        # inputs == outputs
+    def _create_variables(self):
+        """Phase 1: Create variables"""
+        super()._create_variables()
+
+        # Create excess variables if needed
+        if self.element.with_excess:
+            self.excess_input = self.add_variables(lower=0, coords=self._model.get_coords(), short_name='excess_input')
+            self.excess_output = self.add_variables(
+                lower=0, coords=self._model.get_coords(), short_name='excess_output'
+            )
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints (can now access flow.submodel.flow_rate)"""
+        super()._create_constraints()
+
+        # Register flow variables and create balance constraint
         for flow in self.element.inputs + self.element.outputs:
             self.register_variable(flow.submodel.flow_rate, flow.label_full)
+
         inputs = sum([flow.submodel.flow_rate for flow in self.element.inputs])
         outputs = sum([flow.submodel.flow_rate for flow in self.element.outputs])
         eq_bus_balance = self.add_constraints(inputs == outputs, short_name='balance')
 
-        # Fehlerplus/-minus:
+        # Add excess to balance and penalty if needed
         if self.element.with_excess:
             excess_penalty = np.multiply(self._model.hours_per_step, self.element.excess_penalty_per_flow_hour)
-
-            self.excess_input = self.add_variables(lower=0, coords=self._model.get_coords(), short_name='excess_input')
-
-            self.excess_output = self.add_variables(
-                lower=0, coords=self._model.get_coords(), short_name='excess_output'
-            )
-
             eq_bus_balance.lhs -= -self.excess_input + self.excess_output
 
             self._model.effects.add_share_to_penalty(self.label_of_element, (self.excess_input * excess_penalty).sum())
@@ -831,10 +858,13 @@ def __init__(self, model: FlowSystemModel, element: Component):
         self.on_off: OnOffModel | None = None
         super().__init__(model, element)
 
-    def _do_modeling(self):
-        """Initiates all FlowModels"""
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
         all_flows = self.element.inputs + self.element.outputs
+
+        # Set on_off_parameters on flows if needed
         if self.element.on_off_parameters:
             for flow in all_flows:
                 if flow.on_off_parameters is None:
@@ -845,20 +875,13 @@ def _do_modeling(self):
                 if flow.on_off_parameters is None:
                     flow.on_off_parameters = OnOffParameters()
 
+        # Create FlowModels (which creates their variables)
         for flow in all_flows:
             self.add_submodels(flow.create_model(self._model), short_name=flow.label)
 
+        # Create component on variable and OnOffModel if needed
         if self.element.on_off_parameters:
             on = self.add_variables(binary=True, short_name='on', coords=self._model.get_coords())
-            if len(all_flows) == 1:
-                self.add_constraints(on == all_flows[0].submodel.on_off.on, short_name='on')
-            else:
-                flow_ons = [flow.submodel.on_off.on for flow in all_flows]
-                # TODO: Is the EPSILON even necessary?
-                self.add_constraints(on <= sum(flow_ons) + CONFIG.Modeling.epsilon, short_name='on|ub')
-                self.add_constraints(
-                    on >= sum(flow_ons) / (len(flow_ons) + CONFIG.Modeling.epsilon), short_name='on|lb'
-                )
 
             self.on_off = self.add_submodels(
                 OnOffModel(
@@ -872,6 +895,26 @@ def _do_modeling(self):
                 short_name='on_off',
             )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints (can now access flow.submodel.on_off.on)"""
+        super()._create_constraints()
+
+        all_flows = self.element.inputs + self.element.outputs
+
+        # Link component on to flow on states
+        if self.element.on_off_parameters:
+            on = self['on']
+            if len(all_flows) == 1:
+                self.add_constraints(on == all_flows[0].submodel.on_off.on, short_name='on')
+            else:
+                flow_ons = [flow.submodel.on_off.on for flow in all_flows]
+                # TODO: Is the EPSILON even necessary?
+                self.add_constraints(on <= sum(flow_ons) + CONFIG.Modeling.epsilon, short_name='on|ub')
+                self.add_constraints(
+                    on >= sum(flow_ons) / (len(flow_ons) + CONFIG.Modeling.epsilon), short_name='on|lb'
+                )
+
+        # Prevent simultaneous flows
         if self.element.prevent_simultaneous_flows:
             # Simultanious Useage --> Only One FLow is On at a time, but needs a Binary for every flow
             ModelingPrimitives.mutual_exclusivity_constraint(
diff --git a/flixopt/features.py b/flixopt/features.py
index 0d1fc7784..0d61489ba 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -47,12 +47,10 @@ def __init__(
         self.parameters = parameters
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _do_modeling(self):
-        super()._do_modeling()
-        self._create_variables_and_constraints()
-        self._add_effects()
+    def _create_variables(self):
+        """Phase 1: Create variables only"""
+        super()._create_variables()
 
-    def _create_variables_and_constraints(self):
         size_min, size_max = (self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size)
         if self.parameters.linked_periods is not None:
             # Mask size bounds: linked_periods is a binary DataArray that zeros out non-linked periods
@@ -72,6 +70,13 @@ def _create_variables_and_constraints(self):
                 coords=self._model.get_coords(['period', 'scenario']),
                 short_name='invested',
             )
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Create bounding constraints if not mandatory
+        if not self.parameters.mandatory:
             BoundingPatterns.bounds_with_state(
                 self,
                 variable=self.size,
@@ -79,6 +84,7 @@ def _create_variables_and_constraints(self):
                 bounds=(self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size),
             )
 
+        # Create linked periods constraints
         if self.parameters.linked_periods is not None:
             masked_size = self.size.where(self.parameters.linked_periods, drop=True)
             self.add_constraints(
@@ -86,6 +92,9 @@ def _create_variables_and_constraints(self):
                 short_name='linked_periods',
             )
 
+        # Add effects
+        self._add_effects()
+
     def _add_effects(self):
         """Add investment effects"""
         if self.parameters.effects_of_investment:
@@ -173,14 +182,34 @@ def __init__(
         self.parameters = parameters
         super().__init__(model, label_of_element, label_of_model=label_of_model)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables only"""
+        super()._create_variables()
 
         if self.parameters.use_off:
-            off = self.add_variables(binary=True, short_name='off', coords=self._model.get_coords())
-            self.add_constraints(self.on + off == 1, short_name='complementary')
+            self.add_variables(binary=True, short_name='off', coords=self._model.get_coords())
 
-        # 3. Total duration tracking using existing pattern
+        # Switch tracking variables
+        if self.parameters.use_switch_on:
+            self.add_variables(binary=True, short_name='switch|on', coords=self.get_coords())
+            self.add_variables(binary=True, short_name='switch|off', coords=self.get_coords())
+
+            if self.parameters.switch_on_total_max is not None:
+                self.add_variables(
+                    lower=0,
+                    upper=self.parameters.switch_on_total_max,
+                    coords=self._model.get_coords(('period', 'scenario')),
+                    short_name='switch|count',
+                )
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        if self.parameters.use_off:
+            self.add_constraints(self.on + self.off == 1, short_name='complementary')
+
+        # Total duration tracking
         ModelingPrimitives.expression_tracking_variable(
             self,
             tracked_expression=(self.on * self._model.hours_per_step).sum('time'),
@@ -192,11 +221,8 @@ def _do_modeling(self):
             coords=['period', 'scenario'],
         )
 
-        # 4. Switch tracking using existing pattern
+        # Switch tracking constraints
         if self.parameters.use_switch_on:
-            self.add_variables(binary=True, short_name='switch|on', coords=self.get_coords())
-            self.add_variables(binary=True, short_name='switch|off', coords=self.get_coords())
-
             BoundingPatterns.state_transition_bounds(
                 self,
                 state_variable=self.on,
@@ -208,15 +234,9 @@ def _do_modeling(self):
             )
 
             if self.parameters.switch_on_total_max is not None:
-                count = self.add_variables(
-                    lower=0,
-                    upper=self.parameters.switch_on_total_max,
-                    coords=self._model.get_coords(('period', 'scenario')),
-                    short_name='switch|count',
-                )
-                self.add_constraints(count == self.switch_on.sum('time'), short_name='switch|count')
+                self.add_constraints(self.switch_on_nr == self.switch_on.sum('time'), short_name='switch|count')
 
-        # 5. Consecutive on duration using existing pattern
+        # Consecutive on duration tracking
         if self.parameters.use_consecutive_on_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -229,7 +249,7 @@ def _do_modeling(self):
                 previous_duration=self._get_previous_on_duration(),
             )
 
-        # 6. Consecutive off duration using existing pattern
+        # Consecutive off duration tracking
         if self.parameters.use_consecutive_off_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -241,7 +261,6 @@ def _do_modeling(self):
                 duration_dim='time',
                 previous_duration=self._get_previous_off_duration(),
             )
-            # TODO:
 
         self._add_effects()
 
@@ -337,8 +356,9 @@ def __init__(
 
         super().__init__(model, label_of_element, label_of_model)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables only"""
+        super()._create_variables()
         self.inside_piece = self.add_variables(
             binary=True,
             short_name='inside_piece',
@@ -358,6 +378,9 @@ def _do_modeling(self):
             coords=self._model.get_coords(dims=self.dims),
         )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
         # eq:  lambda0(t) + lambda1(t) = inside_piece(t)
         self.add_constraints(self.inside_piece == self.lambda0 + self.lambda1, short_name='inside_piece')
 
@@ -393,13 +416,16 @@ def __init__(
         self.zero_point: linopy.Variable | None = None
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
         # Validate all piecewise variables have the same number of segments
         segment_counts = [len(pw) for pw in self._piecewise_variables.values()]
         if not all(count == segment_counts[0] for count in segment_counts):
             raise ValueError(f'All piecewises must have the same number of pieces, got {segment_counts}')
 
+        # Create PieceModel submodels (which creates their variables)
         for i in range(len(list(self._piecewise_variables.values())[0])):
             new_piece = self.add_submodels(
                 PieceModel(
@@ -412,6 +438,20 @@ def _do_modeling(self):
             )
             self.pieces.append(new_piece)
 
+        # Create zero_point variable if needed
+        if self._zero_point is True:
+            self.zero_point = self.add_variables(
+                coords=self._model.get_coords(self.dims),
+                binary=True,
+                short_name='zero_point',
+            )
+        elif isinstance(self._zero_point, linopy.Variable):
+            self.zero_point = self._zero_point
+
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
         for var_name in self._piecewise_variables:
             variable = self._model.variables[var_name]
             self.add_constraints(
@@ -430,15 +470,7 @@ def _do_modeling(self):
 
             # a) eq: Segment1.onSeg(t) + Segment2.onSeg(t) + ... = 1                Aufenthalt nur in Segmenten erlaubt
             # b) eq: -On(t) + Segment1.onSeg(t) + Segment2.onSeg(t) + ... = 0       zusätzlich kann alles auch Null sein
-            if isinstance(self._zero_point, linopy.Variable):
-                self.zero_point = self._zero_point
-                rhs = self.zero_point
-            elif self._zero_point is True:
-                self.zero_point = self.add_variables(
-                    coords=self._model.get_coords(self.dims),
-                    binary=True,
-                    short_name='zero_point',
-                )
+            if isinstance(self._zero_point, linopy.Variable) or self._zero_point is True:
                 rhs = self.zero_point
             else:
                 rhs = 1
@@ -476,7 +508,10 @@ def __init__(
 
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _do_modeling(self):
+    def _create_variables(self):
+        """Phase 1: Create variables and submodels"""
+        super()._create_variables()
+
         self.shares = {
             effect: self.add_variables(coords=self._model.get_coords(['period', 'scenario']), short_name=effect)
             for effect in self._piecewise_shares
@@ -502,7 +537,11 @@ def _do_modeling(self):
             short_name='PiecewiseEffects',
         )
 
-        # Shares
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Add shares to effects
         self._model.effects.add_share_to_effects(
             name=self.label_of_element,
             expressions={effect: variable * 1 for effect, variable in self.shares.items()},
@@ -542,8 +581,10 @@ def __init__(
 
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _do_modeling(self):
-        super()._do_modeling()
+    def _create_variables(self):
+        """Phase 1: Create variables"""
+        super()._create_variables()
+
         self.total = self.add_variables(
             lower=self._total_min,
             upper=self._total_max,
@@ -551,8 +592,6 @@ def _do_modeling(self):
             name=self.label_full,
             short_name='total',
         )
-        # eq: sum = sum(share_i) # skalar
-        self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
 
         if 'time' in self._dims:
             self.total_per_timestep = self.add_variables(
@@ -562,6 +601,14 @@ def _do_modeling(self):
                 short_name='per_timestep',
             )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # eq: sum = sum(share_i) # skalar
+        self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
+
+        if 'time' in self._dims:
             self._eq_total_per_timestep = self.add_constraints(self.total_per_timestep == 0, short_name='per_timestep')
 
             # Add it to the total
diff --git a/flixopt/structure.py b/flixopt/structure.py
index b6e9e91a5..99a14048d 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -100,12 +100,20 @@ def __init__(self, flow_system: FlowSystem, normalize_weights: bool):
         self.submodels: Submodels = Submodels({})
 
     def do_modeling(self):
+        # Phase 1: Create all models and their variables
         self.effects = self.flow_system.effects.create_model(self)
         for component in self.flow_system.components.values():
             component.create_model(self)
         for bus in self.flow_system.buses.values():
             bus.create_model(self)
 
+        # Phase 2: Create all constraints (all variables exist now)
+        self.effects._create_constraints()
+        for component in self.flow_system.components.values():
+            component.submodel._create_constraints()
+        for bus in self.flow_system.buses.values():
+            bus.submodel._create_constraints()
+
         # Add scenario equality constraints after all elements are modeled
         self._add_scenario_equality_constraints()
 
@@ -1353,7 +1361,7 @@ def __init__(self, model: FlowSystemModel, label_of_element: str, label_of_model
         self.submodels: Submodels = Submodels({})
 
         logger.debug(f'Creating {self.__class__.__name__}  "{self.label_full}"')
-        self._do_modeling()
+        self._create_variables()  # Phase 1: Create variables only
 
     def add_variables(self, short_name: str = None, **kwargs) -> linopy.Variable:
         """Create and register a variable in one step"""
@@ -1503,8 +1511,34 @@ def __repr__(self) -> str:
     def hours_per_step(self):
         return self._model.hours_per_step
 
+    def _create_variables(self):
+        """
+        Phase 1: Create variables only.
+
+        Override in subclasses to create variables and submodels.
+        Do NOT create constraints here - use _create_constraints() for that.
+        """
+        pass
+
+    def _create_constraints(self):
+        """
+        Phase 2: Create constraints.
+
+        Override in subclasses to create constraints.
+        At this point, ALL models and their variables exist, so you can safely
+        reference variables from child models, sibling models, etc.
+        """
+        pass
+
     def _do_modeling(self):
-        """Called at the end of initialization. Override in subclasses to create variables and constraints."""
+        """
+        DEPRECATED: Use _create_variables() and _create_constraints() instead.
+
+        This method is kept for backward compatibility but will be removed in a future version.
+        If you override this in a subclass, split your logic:
+        - Variable creation → _create_variables()
+        - Constraint creation → _create_constraints()
+        """
         pass
 
 

From 5b66a9722f40f813d76b6b094a43a2b536558685 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 12:36:30 +0100
Subject: [PATCH 05/20] Added explicit calls to _create_constraints() for
 nested submodel

---
 flixopt/effects.py  | 9 +++++++++
 flixopt/elements.py | 4 ++++
 flixopt/features.py | 9 +++++++++
 3 files changed, 22 insertions(+)

diff --git a/flixopt/effects.py b/flixopt/effects.py
index 312e8244a..f57eb34bd 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -435,6 +435,10 @@ def _create_constraints(self):
         """Phase 2: Create constraints"""
         super()._create_constraints()
 
+        # Create constraints for share allocation submodels
+        self.periodic._create_constraints()
+        self.temporal._create_constraints()
+
         self.add_constraints(
             self.total == self.temporal.total + self.periodic.total, name=self.label_full, short_name='total'
         )
@@ -697,6 +701,11 @@ def _create_constraints(self):
         """Phase 2: Create constraints"""
         super()._create_constraints()
 
+        # Create constraints for all effect submodels and penalty
+        for effect in self.effects.values():
+            effect.submodel._create_constraints()
+        self.penalty._create_constraints()
+
         # Add cross-effect shares
         self._add_share_between_effects()
 
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 0b821c4b0..38d717ad8 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -901,6 +901,10 @@ def _create_constraints(self):
 
         all_flows = self.element.inputs + self.element.outputs
 
+        # Create constraints for all flow submodels
+        for flow in all_flows:
+            flow.submodel._create_constraints()
+
         # Link component on to flow on states
         if self.element.on_off_parameters:
             on = self['on']
diff --git a/flixopt/features.py b/flixopt/features.py
index 0d61489ba..20cbd7ce4 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -139,6 +139,8 @@ def _add_effects(self):
                 ),
                 short_name='segments',
             )
+            # Create constraints for piecewise effects model
+            self.piecewise_effects._create_constraints()
 
     @property
     def size(self) -> linopy.Variable:
@@ -452,6 +454,10 @@ def _create_constraints(self):
         """Phase 2: Create constraints"""
         super()._create_constraints()
 
+        # Create constraints for all piece submodels
+        for piece in self.pieces:
+            piece._create_constraints()
+
         for var_name in self._piecewise_variables:
             variable = self._model.variables[var_name]
             self.add_constraints(
@@ -541,6 +547,9 @@ def _create_constraints(self):
         """Phase 2: Create constraints"""
         super()._create_constraints()
 
+        # Create constraints for piecewise model
+        self.piecewise_model._create_constraints()
+
         # Add shares to effects
         self._model.effects.add_share_to_effects(
             name=self.label_of_element,

From bb7ee686ffc0fec550d37a91c9e8d0ff259cd058 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 13:04:12 +0100
Subject: [PATCH 06/20] Fix

---
 flixopt/elements.py | 22 ++++++++++++--------
 flixopt/features.py | 50 ++++++++++++++++++++++-----------------------
 2 files changed, 39 insertions(+), 33 deletions(-)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index 38d717ad8..ba593e681 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -582,14 +582,7 @@ def _create_variables(self):
                 short_name='on_off',
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create flow rate bounding constraints
-        self._constraint_flow_rate()
-
-        # Total flow hours tracking
+        # Total flow hours tracking (creates variable + constraint)
         ModelingPrimitives.expression_tracking_variable(
             model=self,
             name=f'{self.label_full}|total_flow_hours',
@@ -602,6 +595,19 @@ def _create_constraints(self):
             short_name='total_flow_hours',
         )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        # Create constraints for nested submodels
+        if self.with_investment:
+            self._investment._create_constraints()
+        if self.with_on_off:
+            self.on_off._create_constraints()
+
+        # Create flow rate bounding constraints
+        self._constraint_flow_rate()
+
         # Load factor constraints
         self._create_bounds_for_load_factor()
 
diff --git a/flixopt/features.py b/flixopt/features.py
index 20cbd7ce4..2bf107fbd 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -204,14 +204,7 @@ def _create_variables(self):
                     short_name='switch|count',
                 )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        if self.parameters.use_off:
-            self.add_constraints(self.on + self.off == 1, short_name='complementary')
-
-        # Total duration tracking
+        # Total duration tracking (creates variable + constraint)
         ModelingPrimitives.expression_tracking_variable(
             self,
             tracked_expression=(self.on * self._model.hours_per_step).sum('time'),
@@ -223,22 +216,7 @@ def _create_constraints(self):
             coords=['period', 'scenario'],
         )
 
-        # Switch tracking constraints
-        if self.parameters.use_switch_on:
-            BoundingPatterns.state_transition_bounds(
-                self,
-                state_variable=self.on,
-                switch_on=self.switch_on,
-                switch_off=self.switch_off,
-                name=f'{self.label_of_model}|switch',
-                previous_state=self._previous_states.isel(time=-1) if self._previous_states is not None else 0,
-                coord='time',
-            )
-
-            if self.parameters.switch_on_total_max is not None:
-                self.add_constraints(self.switch_on_nr == self.switch_on.sum('time'), short_name='switch|count')
-
-        # Consecutive on duration tracking
+        # Consecutive on duration tracking (creates variable + constraints)
         if self.parameters.use_consecutive_on_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -251,7 +229,7 @@ def _create_constraints(self):
                 previous_duration=self._get_previous_on_duration(),
             )
 
-        # Consecutive off duration tracking
+        # Consecutive off duration tracking (creates variable + constraints)
         if self.parameters.use_consecutive_off_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -264,6 +242,28 @@ def _create_constraints(self):
                 previous_duration=self._get_previous_off_duration(),
             )
 
+    def _create_constraints(self):
+        """Phase 2: Create constraints"""
+        super()._create_constraints()
+
+        if self.parameters.use_off:
+            self.add_constraints(self.on + self.off == 1, short_name='complementary')
+
+        # Switch tracking constraints
+        if self.parameters.use_switch_on:
+            BoundingPatterns.state_transition_bounds(
+                self,
+                state_variable=self.on,
+                switch_on=self.switch_on,
+                switch_off=self.switch_off,
+                name=f'{self.label_of_model}|switch',
+                previous_state=self._previous_states.isel(time=-1) if self._previous_states is not None else 0,
+                coord='time',
+            )
+
+            if self.parameters.switch_on_total_max is not None:
+                self.add_constraints(self.switch_on_nr == self.switch_on.sum('time'), short_name='switch|count')
+
         self._add_effects()
 
     def _add_effects(self):

From 01efda51a97d1c897c7de0cb6f05c33c099cfcc9 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 13:23:41 +0100
Subject: [PATCH 07/20] Revert changes

---
 flixopt/components.py | 52 ++++++++++++---------------
 flixopt/effects.py    | 29 ++++-----------
 flixopt/elements.py   | 46 ++++++------------------
 flixopt/features.py   | 84 ++++++++++++++++---------------------------
 flixopt/structure.py  | 38 +++-----------------
 5 files changed, 75 insertions(+), 174 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index c215b1cd2..557bd1d86 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -734,9 +734,9 @@ def __init__(self, model: FlowSystemModel, element: Transmission):
 
         super().__init__(model, element)
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints (FlowModels already created by parent)"""
-        super()._create_constraints()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # first direction
         self.create_transmission_equation('dir1', self.element.in1, self.element.out1)
@@ -775,9 +775,9 @@ def __init__(self, model: FlowSystemModel, element: LinearConverter):
         self.piecewise_conversion: PiecewiseConversion | None = None
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables (parent creates FlowModels)"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # Create PiecewiseModel if needed (after FlowModels are created)
         if not self.element.conversion_factors:
@@ -799,10 +799,6 @@ def _create_variables(self):
                 short_name='PiecewiseConversion',
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
         # Create conversion factor constraints if specified
         if self.element.conversion_factors:
             all_input_flows = set(self.element.inputs)
@@ -829,10 +825,11 @@ class StorageModel(ComponentModel):
     def __init__(self, model: FlowSystemModel, element: Storage):
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables (parent creates FlowModels)"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
+        # Create variables
         lb, ub = self._absolute_charge_state_bounds
         self.add_variables(
             lower=lb,
@@ -843,22 +840,7 @@ def _create_variables(self):
 
         self.add_variables(coords=self._model.get_coords(), short_name='netto_discharge')
 
-        # Create InvestmentModel for capacity if needed
-        if isinstance(self.element.capacity_in_flow_hours, InvestParameters):
-            self.add_submodels(
-                InvestmentModel(
-                    model=self._model,
-                    label_of_element=self.label_of_element,
-                    label_of_model=self.label_of_element,
-                    parameters=self.element.capacity_in_flow_hours,
-                ),
-                short_name='investment',
-            )
-
-    def _create_constraints(self):
-        """Phase 2: Create constraints (can now access flow.submodel.flow_rate)"""
-        super()._create_constraints()
-
+        # Create constraints (can now access flow.submodel.flow_rate)
         # netto_discharge:
         # eq: nettoFlow(t) - discharging(t) + charging(t) = 0
         self.add_constraints(
@@ -883,8 +865,18 @@ def _create_constraints(self):
             short_name='charge_state',
         )
 
-        # Bounding constraints for investment
+        # Create InvestmentModel and bounding constraints for investment
         if isinstance(self.element.capacity_in_flow_hours, InvestParameters):
+            self.add_submodels(
+                InvestmentModel(
+                    model=self._model,
+                    label_of_element=self.label_of_element,
+                    label_of_model=self.label_of_element,
+                    parameters=self.element.capacity_in_flow_hours,
+                ),
+                short_name='investment',
+            )
+
             BoundingPatterns.scaled_bounds(
                 self,
                 variable=self.charge_state,
diff --git a/flixopt/effects.py b/flixopt/effects.py
index f57eb34bd..074809a5a 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -393,9 +393,9 @@ class EffectModel(ElementModel):
     def __init__(self, model: FlowSystemModel, element: Effect):
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         self.total: linopy.Variable | None = None
         self.periodic: ShareAllocationModel = self.add_submodels(
@@ -431,14 +431,6 @@ def _create_variables(self):
             name=self.label_full,
         )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create constraints for share allocation submodels
-        self.periodic._create_constraints()
-        self.temporal._create_constraints()
-
         self.add_constraints(
             self.total == self.temporal.total + self.periodic.total, name=self.label_full, short_name='total'
         )
@@ -683,9 +675,9 @@ def add_share_to_penalty(self, name: str, expression: linopy.LinearExpression) -
             raise TypeError(f'Penalty shares must be scalar expressions! ({expression.ndim=})')
         self.penalty.add_share(name, expression, dims=())
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # Create EffectModel for each effect
         for effect in self.effects.values():
@@ -697,15 +689,6 @@ def _create_variables(self):
             short_name='penalty',
         )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create constraints for all effect submodels and penalty
-        for effect in self.effects.values():
-            effect.submodel._create_constraints()
-        self.penalty._create_constraints()
-
         # Add cross-effect shares
         self._add_share_between_effects()
 
diff --git a/flixopt/elements.py b/flixopt/elements.py
index ba593e681..1d91ef454 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -553,9 +553,9 @@ class FlowModel(ElementModel):
     def __init__(self, model: FlowSystemModel, element: Flow):
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # Main flow rate variable
         self.add_variables(
@@ -565,7 +565,7 @@ def _create_variables(self):
             short_name='flow_rate',
         )
 
-        # Create investment and/or on_off submodels (creates their variables)
+        # Create investment and/or on_off submodels (creates their variables and constraints)
         if self.with_investment:
             self._create_investment_model()
         if self.with_on_off:
@@ -595,16 +595,6 @@ def _create_variables(self):
             short_name='total_flow_hours',
         )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create constraints for nested submodels
-        if self.with_investment:
-            self._investment._create_constraints()
-        if self.with_on_off:
-            self.on_off._create_constraints()
-
         # Create flow rate bounding constraints
         self._constraint_flow_rate()
 
@@ -811,9 +801,9 @@ def __init__(self, model: FlowSystemModel, element: Bus):
         self.excess_output: linopy.Variable | None = None
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # Create excess variables if needed
         if self.element.with_excess:
@@ -822,10 +812,6 @@ def _create_variables(self):
                 lower=0, coords=self._model.get_coords(), short_name='excess_output'
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints (can now access flow.submodel.flow_rate)"""
-        super()._create_constraints()
-
         # Register flow variables and create balance constraint
         for flow in self.element.inputs + self.element.outputs:
             self.register_variable(flow.submodel.flow_rate, flow.label_full)
@@ -864,9 +850,9 @@ def __init__(self, model: FlowSystemModel, element: Component):
         self.on_off: OnOffModel | None = None
         super().__init__(model, element)
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         all_flows = self.element.inputs + self.element.outputs
 
@@ -881,7 +867,7 @@ def _create_variables(self):
                 if flow.on_off_parameters is None:
                     flow.on_off_parameters = OnOffParameters()
 
-        # Create FlowModels (which creates their variables)
+        # Create FlowModels (which creates their variables and constraints)
         for flow in all_flows:
             self.add_submodels(flow.create_model(self._model), short_name=flow.label)
 
@@ -901,16 +887,6 @@ def _create_variables(self):
                 short_name='on_off',
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints (can now access flow.submodel.on_off.on)"""
-        super()._create_constraints()
-
-        all_flows = self.element.inputs + self.element.outputs
-
-        # Create constraints for all flow submodels
-        for flow in all_flows:
-            flow.submodel._create_constraints()
-
         # Link component on to flow on states
         if self.element.on_off_parameters:
             on = self['on']
diff --git a/flixopt/features.py b/flixopt/features.py
index 2bf107fbd..877d0efcc 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -47,10 +47,11 @@ def __init__(
         self.parameters = parameters
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables only"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
+        # Create variables
         size_min, size_max = (self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size)
         if self.parameters.linked_periods is not None:
             # Mask size bounds: linked_periods is a binary DataArray that zeros out non-linked periods
@@ -71,11 +72,7 @@ def _create_variables(self):
                 short_name='invested',
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create bounding constraints if not mandatory
+        # Create constraints
         if not self.parameters.mandatory:
             BoundingPatterns.bounds_with_state(
                 self,
@@ -84,7 +81,6 @@ def _create_constraints(self):
                 bounds=(self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size),
             )
 
-        # Create linked periods constraints
         if self.parameters.linked_periods is not None:
             masked_size = self.size.where(self.parameters.linked_periods, drop=True)
             self.add_constraints(
@@ -139,8 +135,6 @@ def _add_effects(self):
                 ),
                 short_name='segments',
             )
-            # Create constraints for piecewise effects model
-            self.piecewise_effects._create_constraints()
 
     @property
     def size(self) -> linopy.Variable:
@@ -184,10 +178,11 @@ def __init__(
         self.parameters = parameters
         super().__init__(model, label_of_element, label_of_model=label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables only"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
+        # Create variables
         if self.parameters.use_off:
             self.add_variables(binary=True, short_name='off', coords=self._model.get_coords())
 
@@ -242,10 +237,7 @@ def _create_variables(self):
                 previous_duration=self._get_previous_off_duration(),
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
+        # Create constraints
         if self.parameters.use_off:
             self.add_constraints(self.on + self.off == 1, short_name='complementary')
 
@@ -358,9 +350,11 @@ def __init__(
 
         super().__init__(model, label_of_element, label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables only"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
+
+        # Create variables
         self.inside_piece = self.add_variables(
             binary=True,
             short_name='inside_piece',
@@ -380,9 +374,7 @@ def _create_variables(self):
             coords=self._model.get_coords(dims=self.dims),
         )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
+        # Create constraints
         # eq:  lambda0(t) + lambda1(t) = inside_piece(t)
         self.add_constraints(self.inside_piece == self.lambda0 + self.lambda1, short_name='inside_piece')
 
@@ -418,16 +410,16 @@ def __init__(
         self.zero_point: linopy.Variable | None = None
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
         # Validate all piecewise variables have the same number of segments
         segment_counts = [len(pw) for pw in self._piecewise_variables.values()]
         if not all(count == segment_counts[0] for count in segment_counts):
             raise ValueError(f'All piecewises must have the same number of pieces, got {segment_counts}')
 
-        # Create PieceModel submodels (which creates their variables)
+        # Create PieceModel submodels (which creates their variables and constraints)
         for i in range(len(list(self._piecewise_variables.values())[0])):
             new_piece = self.add_submodels(
                 PieceModel(
@@ -450,14 +442,7 @@ def _create_variables(self):
         elif isinstance(self._zero_point, linopy.Variable):
             self.zero_point = self._zero_point
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create constraints for all piece submodels
-        for piece in self.pieces:
-            piece._create_constraints()
-
+        # Create constraints
         for var_name in self._piecewise_variables:
             variable = self._model.variables[var_name]
             self.add_constraints(
@@ -514,10 +499,11 @@ def __init__(
 
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables and submodels"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
+        # Create variables
         self.shares = {
             effect: self.add_variables(coords=self._model.get_coords(['period', 'scenario']), short_name=effect)
             for effect in self._piecewise_shares
@@ -531,6 +517,7 @@ def _create_variables(self):
             },
         }
 
+        # Create piecewise model (which creates its variables and constraints)
         self.piecewise_model = self.add_submodels(
             PiecewiseModel(
                 model=self._model,
@@ -543,13 +530,6 @@ def _create_variables(self):
             short_name='PiecewiseEffects',
         )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
-        # Create constraints for piecewise model
-        self.piecewise_model._create_constraints()
-
         # Add shares to effects
         self._model.effects.add_share_to_effects(
             name=self.label_of_element,
@@ -590,10 +570,11 @@ def __init__(
 
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
-    def _create_variables(self):
-        """Phase 1: Create variables"""
-        super()._create_variables()
+    def _do_modeling(self):
+        """Create variables, constraints, and nested submodels"""
+        super()._do_modeling()
 
+        # Create variables
         self.total = self.add_variables(
             lower=self._total_min,
             upper=self._total_max,
@@ -610,10 +591,7 @@ def _create_variables(self):
                 short_name='per_timestep',
             )
 
-    def _create_constraints(self):
-        """Phase 2: Create constraints"""
-        super()._create_constraints()
-
+        # Create constraints
         # eq: sum = sum(share_i) # skalar
         self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
 
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 99a14048d..794608ea0 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -100,20 +100,13 @@ def __init__(self, flow_system: FlowSystem, normalize_weights: bool):
         self.submodels: Submodels = Submodels({})
 
     def do_modeling(self):
-        # Phase 1: Create all models and their variables
+        # Create all element models
         self.effects = self.flow_system.effects.create_model(self)
         for component in self.flow_system.components.values():
             component.create_model(self)
         for bus in self.flow_system.buses.values():
             bus.create_model(self)
 
-        # Phase 2: Create all constraints (all variables exist now)
-        self.effects._create_constraints()
-        for component in self.flow_system.components.values():
-            component.submodel._create_constraints()
-        for bus in self.flow_system.buses.values():
-            bus.submodel._create_constraints()
-
         # Add scenario equality constraints after all elements are modeled
         self._add_scenario_equality_constraints()
 
@@ -1361,7 +1354,7 @@ def __init__(self, model: FlowSystemModel, label_of_element: str, label_of_model
         self.submodels: Submodels = Submodels({})
 
         logger.debug(f'Creating {self.__class__.__name__}  "{self.label_full}"')
-        self._create_variables()  # Phase 1: Create variables only
+        self._do_modeling()
 
     def add_variables(self, short_name: str = None, **kwargs) -> linopy.Variable:
         """Create and register a variable in one step"""
@@ -1511,33 +1504,12 @@ def __repr__(self) -> str:
     def hours_per_step(self):
         return self._model.hours_per_step
 
-    def _create_variables(self):
-        """
-        Phase 1: Create variables only.
-
-        Override in subclasses to create variables and submodels.
-        Do NOT create constraints here - use _create_constraints() for that.
-        """
-        pass
-
-    def _create_constraints(self):
-        """
-        Phase 2: Create constraints.
-
-        Override in subclasses to create constraints.
-        At this point, ALL models and their variables exist, so you can safely
-        reference variables from child models, sibling models, etc.
-        """
-        pass
-
     def _do_modeling(self):
         """
-        DEPRECATED: Use _create_variables() and _create_constraints() instead.
+        Override in subclasses to create variables, constraints, and submodels.
 
-        This method is kept for backward compatibility but will be removed in a future version.
-        If you override this in a subclass, split your logic:
-        - Variable creation → _create_variables()
-        - Constraint creation → _create_constraints()
+        This method is called during __init__. Create all nested submodels first
+        (so their variables exist), then create constraints that reference those variables.
         """
         pass
 

From 87f08e6daa7c53240208dcb09406ff8a4e1cff78 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 13:39:37 +0100
Subject: [PATCH 08/20] Revert changes

---
 flixopt/components.py | 37 ++++++++++++++++++-------------------
 1 file changed, 18 insertions(+), 19 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 557bd1d86..6cd04b7e3 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -779,8 +779,24 @@ def _do_modeling(self):
         """Create variables, constraints, and nested submodels"""
         super()._do_modeling()
 
-        # Create PiecewiseModel if needed (after FlowModels are created)
-        if not self.element.conversion_factors:
+        # Create conversion factor constraints if specified
+        if self.element.conversion_factors:
+            all_input_flows = set(self.element.inputs)
+            all_output_flows = set(self.element.outputs)
+
+            # für alle linearen Gleichungen:
+            for i, conv_factors in enumerate(self.element.conversion_factors):
+                used_flows = set([self.element.flows[flow_label] for flow_label in conv_factors])
+                used_inputs: set[Flow] = all_input_flows & used_flows
+                used_outputs: set[Flow] = all_output_flows & used_flows
+
+                self.add_constraints(
+                    sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_inputs])
+                    == sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_outputs]),
+                    short_name=f'conversion_{i}',
+                )
+
+        else:
             # TODO: Improve Inclusion of OnOffParameters. Instead of creating a Binary in every flow, the binary could only be part of the Piece itself
             piecewise_conversion = {
                 self.element.flows[flow].submodel.flow_rate.name: piecewise
@@ -799,23 +815,6 @@ def _do_modeling(self):
                 short_name='PiecewiseConversion',
             )
 
-        # Create conversion factor constraints if specified
-        if self.element.conversion_factors:
-            all_input_flows = set(self.element.inputs)
-            all_output_flows = set(self.element.outputs)
-
-            # für alle linearen Gleichungen:
-            for i, conv_factors in enumerate(self.element.conversion_factors):
-                used_flows = set([self.element.flows[flow_label] for flow_label in conv_factors])
-                used_inputs: set[Flow] = all_input_flows & used_flows
-                used_outputs: set[Flow] = all_output_flows & used_flows
-
-                self.add_constraints(
-                    sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_inputs])
-                    == sum([flow.submodel.flow_rate * conv_factors[flow.label] for flow in used_outputs]),
-                    short_name=f'conversion_{i}',
-                )
-
 
 class StorageModel(ComponentModel):
     """Submodel of Storage"""

From e2d09d51eb7a8801fdf1d5a4794e43a0c8257898 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 13:42:48 +0100
Subject: [PATCH 09/20] Revert changes

---
 flixopt/elements.py | 66 +++++++++++++++------------------------------
 1 file changed, 22 insertions(+), 44 deletions(-)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index 1d91ef454..3e3bf5d2f 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -565,23 +565,6 @@ def _do_modeling(self):
             short_name='flow_rate',
         )
 
-        # Create investment and/or on_off submodels (creates their variables and constraints)
-        if self.with_investment:
-            self._create_investment_model()
-        if self.with_on_off:
-            on = self.add_variables(binary=True, short_name='on', coords=self._model.get_coords())
-            self.add_submodels(
-                OnOffModel(
-                    model=self._model,
-                    label_of_element=self.label_of_element,
-                    parameters=self.element.on_off_parameters,
-                    on_variable=on,
-                    previous_states=self.previous_states,
-                    label_of_model=self.label_of_element,
-                ),
-                short_name='on_off',
-            )
-
         # Total flow hours tracking (creates variable + constraint)
         ModelingPrimitives.expression_tracking_variable(
             model=self,
@@ -595,9 +578,6 @@ def _do_modeling(self):
             short_name='total_flow_hours',
         )
 
-        # Create flow rate bounding constraints
-        self._constraint_flow_rate()
-
         # Load factor constraints
         self._create_bounds_for_load_factor()
 
@@ -637,6 +617,7 @@ def _constraint_flow_rate(self):
 
         elif self.with_on_off and not self.with_investment:
             # OnOff, but no Investment
+            self._create_on_off_model()
             bounds = self.relative_flow_rate_bounds
             BoundingPatterns.bounds_with_state(
                 self,
@@ -647,6 +628,7 @@ def _constraint_flow_rate(self):
 
         elif self.with_investment and not self.with_on_off:
             # Investment, but no OnOff
+            self._create_investment_model()
             BoundingPatterns.scaled_bounds(
                 self,
                 variable=self.flow_rate,
@@ -656,6 +638,9 @@ def _constraint_flow_rate(self):
 
         elif self.with_investment and self.with_on_off:
             # Investment and OnOff
+            self._create_investment_model()
+            self._create_on_off_model()
+
             BoundingPatterns.scaled_bounds_with_state(
                 model=self,
                 variable=self.flow_rate,
@@ -804,18 +789,9 @@ def __init__(self, model: FlowSystemModel, element: Bus):
     def _do_modeling(self):
         """Create variables, constraints, and nested submodels"""
         super()._do_modeling()
-
-        # Create excess variables if needed
-        if self.element.with_excess:
-            self.excess_input = self.add_variables(lower=0, coords=self._model.get_coords(), short_name='excess_input')
-            self.excess_output = self.add_variables(
-                lower=0, coords=self._model.get_coords(), short_name='excess_output'
-            )
-
-        # Register flow variables and create balance constraint
+        # inputs == outputs
         for flow in self.element.inputs + self.element.outputs:
             self.register_variable(flow.submodel.flow_rate, flow.label_full)
-
         inputs = sum([flow.submodel.flow_rate for flow in self.element.inputs])
         outputs = sum([flow.submodel.flow_rate for flow in self.element.outputs])
         eq_bus_balance = self.add_constraints(inputs == outputs, short_name='balance')
@@ -823,6 +799,13 @@ def _do_modeling(self):
         # Add excess to balance and penalty if needed
         if self.element.with_excess:
             excess_penalty = np.multiply(self._model.hours_per_step, self.element.excess_penalty_per_flow_hour)
+
+            self.excess_input = self.add_variables(lower=0, coords=self._model.get_coords(), short_name='excess_input')
+
+            self.excess_output = self.add_variables(
+                lower=0, coords=self._model.get_coords(), short_name='excess_output'
+            )
+
             eq_bus_balance.lhs -= -self.excess_input + self.excess_output
 
             self._model.effects.add_share_to_penalty(self.label_of_element, (self.excess_input * excess_penalty).sum())
@@ -874,6 +857,15 @@ def _do_modeling(self):
         # Create component on variable and OnOffModel if needed
         if self.element.on_off_parameters:
             on = self.add_variables(binary=True, short_name='on', coords=self._model.get_coords())
+            if len(all_flows) == 1:
+                self.add_constraints(on == all_flows[0].submodel.on_off.on, short_name='on')
+            else:
+                flow_ons = [flow.submodel.on_off.on for flow in all_flows]
+                # TODO: Is the EPSILON even necessary?
+                self.add_constraints(on <= sum(flow_ons) + CONFIG.Modeling.epsilon, short_name='on|ub')
+                self.add_constraints(
+                    on >= sum(flow_ons) / (len(flow_ons) + CONFIG.Modeling.epsilon), short_name='on|lb'
+                )
 
             self.on_off = self.add_submodels(
                 OnOffModel(
@@ -887,20 +879,6 @@ def _do_modeling(self):
                 short_name='on_off',
             )
 
-        # Link component on to flow on states
-        if self.element.on_off_parameters:
-            on = self['on']
-            if len(all_flows) == 1:
-                self.add_constraints(on == all_flows[0].submodel.on_off.on, short_name='on')
-            else:
-                flow_ons = [flow.submodel.on_off.on for flow in all_flows]
-                # TODO: Is the EPSILON even necessary?
-                self.add_constraints(on <= sum(flow_ons) + CONFIG.Modeling.epsilon, short_name='on|ub')
-                self.add_constraints(
-                    on >= sum(flow_ons) / (len(flow_ons) + CONFIG.Modeling.epsilon), short_name='on|lb'
-                )
-
-        # Prevent simultaneous flows
         if self.element.prevent_simultaneous_flows:
             # Simultanious Useage --> Only One FLow is On at a time, but needs a Binary for every flow
             ModelingPrimitives.mutual_exclusivity_constraint(

From b1d754410dac8356afa921d30d2e8180e16d2a51 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 13:45:26 +0100
Subject: [PATCH 10/20] Revert changes

---
 flixopt/features.py | 106 ++++++++++++++++++--------------------------
 1 file changed, 44 insertions(+), 62 deletions(-)

diff --git a/flixopt/features.py b/flixopt/features.py
index 877d0efcc..d01c0c8ab 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -48,10 +48,11 @@ def __init__(
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
     def _do_modeling(self):
-        """Create variables, constraints, and nested submodels"""
         super()._do_modeling()
+        self._create_variables_and_constraints()
+        self._add_effects()
 
-        # Create variables
+    def _create_variables_and_constraints(self):
         size_min, size_max = (self.parameters.minimum_or_fixed_size, self.parameters.maximum_or_fixed_size)
         if self.parameters.linked_periods is not None:
             # Mask size bounds: linked_periods is a binary DataArray that zeros out non-linked periods
@@ -71,9 +72,6 @@ def _do_modeling(self):
                 coords=self._model.get_coords(['period', 'scenario']),
                 short_name='invested',
             )
-
-        # Create constraints
-        if not self.parameters.mandatory:
             BoundingPatterns.bounds_with_state(
                 self,
                 variable=self.size,
@@ -88,9 +86,6 @@ def _do_modeling(self):
                 short_name='linked_periods',
             )
 
-        # Add effects
-        self._add_effects()
-
     def _add_effects(self):
         """Add investment effects"""
         if self.parameters.effects_of_investment:
@@ -182,24 +177,11 @@ def _do_modeling(self):
         """Create variables, constraints, and nested submodels"""
         super()._do_modeling()
 
-        # Create variables
         if self.parameters.use_off:
-            self.add_variables(binary=True, short_name='off', coords=self._model.get_coords())
+            off = self.add_variables(binary=True, short_name='off', coords=self._model.get_coords())
+            self.add_constraints(self.on + off == 1, short_name='complementary')
 
-        # Switch tracking variables
-        if self.parameters.use_switch_on:
-            self.add_variables(binary=True, short_name='switch|on', coords=self.get_coords())
-            self.add_variables(binary=True, short_name='switch|off', coords=self.get_coords())
-
-            if self.parameters.switch_on_total_max is not None:
-                self.add_variables(
-                    lower=0,
-                    upper=self.parameters.switch_on_total_max,
-                    coords=self._model.get_coords(('period', 'scenario')),
-                    short_name='switch|count',
-                )
-
-        # Total duration tracking (creates variable + constraint)
+        # 3. Total duration tracking using existing pattern
         ModelingPrimitives.expression_tracking_variable(
             self,
             tracked_expression=(self.on * self._model.hours_per_step).sum('time'),
@@ -211,7 +193,31 @@ def _do_modeling(self):
             coords=['period', 'scenario'],
         )
 
-        # Consecutive on duration tracking (creates variable + constraints)
+        # 4. Switch tracking using existing pattern
+        if self.parameters.use_switch_on:
+            self.add_variables(binary=True, short_name='switch|on', coords=self.get_coords())
+            self.add_variables(binary=True, short_name='switch|off', coords=self.get_coords())
+
+            BoundingPatterns.state_transition_bounds(
+                self,
+                state_variable=self.on,
+                switch_on=self.switch_on,
+                switch_off=self.switch_off,
+                name=f'{self.label_of_model}|switch',
+                previous_state=self._previous_states.isel(time=-1) if self._previous_states is not None else 0,
+                coord='time',
+            )
+
+            if self.parameters.switch_on_total_max is not None:
+                count = self.add_variables(
+                    lower=0,
+                    upper=self.parameters.switch_on_total_max,
+                    coords=self._model.get_coords(('period', 'scenario')),
+                    short_name='switch|count',
+                )
+                self.add_constraints(count == self.switch_on.sum('time'), short_name='switch|count')
+
+        # 5. Consecutive on duration using existing pattern
         if self.parameters.use_consecutive_on_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -224,7 +230,7 @@ def _do_modeling(self):
                 previous_duration=self._get_previous_on_duration(),
             )
 
-        # Consecutive off duration tracking (creates variable + constraints)
+        # 6. Consecutive off duration using existing pattern
         if self.parameters.use_consecutive_off_hours:
             ModelingPrimitives.consecutive_duration_tracking(
                 self,
@@ -236,25 +242,7 @@ def _do_modeling(self):
                 duration_dim='time',
                 previous_duration=self._get_previous_off_duration(),
             )
-
-        # Create constraints
-        if self.parameters.use_off:
-            self.add_constraints(self.on + self.off == 1, short_name='complementary')
-
-        # Switch tracking constraints
-        if self.parameters.use_switch_on:
-            BoundingPatterns.state_transition_bounds(
-                self,
-                state_variable=self.on,
-                switch_on=self.switch_on,
-                switch_off=self.switch_off,
-                name=f'{self.label_of_model}|switch',
-                previous_state=self._previous_states.isel(time=-1) if self._previous_states is not None else 0,
-                coord='time',
-            )
-
-            if self.parameters.switch_on_total_max is not None:
-                self.add_constraints(self.switch_on_nr == self.switch_on.sum('time'), short_name='switch|count')
+            # TODO:
 
         self._add_effects()
 
@@ -432,17 +420,6 @@ def _do_modeling(self):
             )
             self.pieces.append(new_piece)
 
-        # Create zero_point variable if needed
-        if self._zero_point is True:
-            self.zero_point = self.add_variables(
-                coords=self._model.get_coords(self.dims),
-                binary=True,
-                short_name='zero_point',
-            )
-        elif isinstance(self._zero_point, linopy.Variable):
-            self.zero_point = self._zero_point
-
-        # Create constraints
         for var_name in self._piecewise_variables:
             variable = self._model.variables[var_name]
             self.add_constraints(
@@ -461,7 +438,15 @@ def _do_modeling(self):
 
             # a) eq: Segment1.onSeg(t) + Segment2.onSeg(t) + ... = 1                Aufenthalt nur in Segmenten erlaubt
             # b) eq: -On(t) + Segment1.onSeg(t) + Segment2.onSeg(t) + ... = 0       zusätzlich kann alles auch Null sein
-            if isinstance(self._zero_point, linopy.Variable) or self._zero_point is True:
+            if isinstance(self._zero_point, linopy.Variable):
+                self.zero_point = self._zero_point
+                rhs = self.zero_point
+            elif self._zero_point is True:
+                self.zero_point = self.add_variables(
+                    coords=self._model.get_coords(self.dims),
+                    binary=True,
+                    short_name='zero_point',
+                )
                 rhs = self.zero_point
             else:
                 rhs = 1
@@ -582,6 +567,8 @@ def _do_modeling(self):
             name=self.label_full,
             short_name='total',
         )
+        # eq: sum = sum(share_i) # skalar
+        self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
 
         if 'time' in self._dims:
             self.total_per_timestep = self.add_variables(
@@ -591,11 +578,6 @@ def _do_modeling(self):
                 short_name='per_timestep',
             )
 
-        # Create constraints
-        # eq: sum = sum(share_i) # skalar
-        self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
-
-        if 'time' in self._dims:
             self._eq_total_per_timestep = self.add_constraints(self.total_per_timestep == 0, short_name='per_timestep')
 
             # Add it to the total

From b63250d1bdc9d6477a888e536fb41d92e8bc3bf9 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 16:46:46 +0100
Subject: [PATCH 11/20] Revert changes

---
 flixopt/elements.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index 3e3bf5d2f..b72128cc3 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -565,6 +565,8 @@ def _do_modeling(self):
             short_name='flow_rate',
         )
 
+        self._constraint_flow_rate()
+
         # Total flow hours tracking (creates variable + constraint)
         ModelingPrimitives.expression_tracking_variable(
             model=self,

From 22f364bd8f513e0c2e98d10b39f6bcd9377a58f4 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 14 Nov 2025 17:07:00 +0100
Subject: [PATCH 12/20] Add more validation to add_elements()

---
 flixopt/flow_system.py | 42 ++++++++++++++++++++++++++++++++++--------
 1 file changed, 34 insertions(+), 8 deletions(-)

diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index af50a758d..03dfc7a54 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -621,17 +621,29 @@ def add_elements(self, *elements: Element) -> None:
                 stacklevel=2,
             )
             self._connected_and_transformed = False
+
         for new_element in list(elements):
+            # Validate element type first
+            if not isinstance(new_element, (Component, Effect, Bus)):
+                raise TypeError(
+                    f'Tried to add incompatible object to FlowSystem: {type(new_element)=}: {new_element=} '
+                )
+
+            # Common validations for all element types (before any state changes)
+            self._check_if_element_already_assigned(new_element)
+            self._check_if_element_is_unique(new_element)
+
+            # Dispatch to type-specific handlers
             if isinstance(new_element, Component):
                 self._add_components(new_element)
             elif isinstance(new_element, Effect):
                 self._add_effects(new_element)
             elif isinstance(new_element, Bus):
                 self._add_buses(new_element)
-            else:
-                raise TypeError(
-                    f'Tried to add incompatible object to FlowSystem: {type(new_element)=}: {new_element=} '
-                )
+
+            # Log registration
+            element_type = type(new_element).__name__
+            logger.info(f'Registered new {element_type}: {new_element.label_full}')
 
     def create_model(self, normalize_weights: bool = True) -> FlowSystemModel:
         """
@@ -779,6 +791,24 @@ def _check_if_element_is_unique(self, element: Element) -> None:
         if element.label_full in self:
             raise ValueError(f'Label of Element {element.label_full} already used in another element!')
 
+    def _check_if_element_already_assigned(self, element: Element) -> None:
+        """
+        Check if element already belongs to another FlowSystem.
+
+        Args:
+            element: Element to check
+
+        Raises:
+            ValueError: If element is already assigned to a different FlowSystem
+        """
+        if element._flow_system is not None and element._flow_system is not self:
+            raise ValueError(
+                f'Element "{element.label_full}" is already assigned to another FlowSystem. '
+                f'Each element can only belong to one FlowSystem at a time. '
+                f'To use this element in multiple systems, create a copy: '
+                f'flow_system.add_elements(element.copy())'
+            )
+
     def _validate_system_integrity(self) -> None:
         """
         Validate cross-element references to ensure system consistency.
@@ -809,16 +839,12 @@ def _add_effects(self, *args: Effect) -> None:
 
     def _add_components(self, *components: Component) -> None:
         for new_component in list(components):
-            logger.info(f'Registered new Component: {new_component.label_full}')
-            self._check_if_element_is_unique(new_component)  # check if already exists:
             new_component._set_flow_system(self)  # Link element to FlowSystem
             self.components.add(new_component)  # Add to existing components
             self._flows_cache = None  # Invalidate flows cache
 
     def _add_buses(self, *buses: Bus):
         for new_bus in list(buses):
-            logger.info(f'Registered new Bus: {new_bus.label_full}')
-            self._check_if_element_is_unique(new_bus)  # check if already exists:
             new_bus._set_flow_system(self)  # Link element to FlowSystem
             self.buses.add(new_bus)  # Add to existing buses
             self._flows_cache = None  # Invalidate flows cache

From 155428ad67dc11a4193a03598897fc346400c781 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 16 Nov 2025 17:51:22 +0100
Subject: [PATCH 13/20]  Fixed inconsistent argument passing for
 _fit_effect_coords

---
 flixopt/interface.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/flixopt/interface.py b/flixopt/interface.py
index 3f21c027b..e03738422 100644
--- a/flixopt/interface.py
+++ b/flixopt/interface.py
@@ -1293,9 +1293,15 @@ def __init__(
         self.force_switch_on: bool = force_switch_on
 
     def transform_data(self, name_prefix: str = '') -> None:
-        self.effects_per_switch_on = self._fit_effect_coords(name_prefix, self.effects_per_switch_on, 'per_switch_on')
+        self.effects_per_switch_on = self._fit_effect_coords(
+            prefix=name_prefix,
+            effect_values=self.effects_per_switch_on,
+            suffix='per_switch_on',
+        )
         self.effects_per_running_hour = self._fit_effect_coords(
-            name_prefix, self.effects_per_running_hour, 'per_running_hour'
+            prefix=name_prefix,
+            effect_values=self.effects_per_running_hour,
+            suffix='per_running_hour',
         )
         self.consecutive_on_hours_min = self._fit_coords(
             f'{name_prefix}|consecutive_on_hours_min', self.consecutive_on_hours_min

From 705dd3f1261d8040e84578b0edffb3f463d126a5 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 16 Nov 2025 17:52:11 +0100
Subject: [PATCH 14/20]  Refactored _set_flow_system to be a method on each
 Interface subclass

---
 flixopt/components.py |  6 ++++++
 flixopt/elements.py   | 24 +++++++++++++++++++++++-
 flixopt/interface.py  | 25 +++++++++++++++++++++++++
 flixopt/structure.py  | 20 +++-----------------
 4 files changed, 57 insertions(+), 18 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 6cd04b7e3..2b091d7fa 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -434,6 +434,12 @@ def create_model(self, model: FlowSystemModel) -> StorageModel:
         self.submodel = StorageModel(model, self)
         return self.submodel
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to parent Component and capacity_in_flow_hours if it's InvestParameters."""
+        super()._set_flow_system(flow_system)
+        if isinstance(self.capacity_in_flow_hours, InvestParameters):
+            self.capacity_in_flow_hours._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         super().transform_data(prefix)
diff --git a/flixopt/elements.py b/flixopt/elements.py
index b72128cc3..b57cda292 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -17,7 +17,7 @@
 from .features import InvestmentModel, OnOffModel
 from .interface import InvestParameters, OnOffParameters
 from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilitiesAbstract
-from .structure import Element, ElementModel, FlowSystemModel, register_class_for_io
+from .structure import Element, ElementModel, FlowSystemModel, Interface, register_class_for_io
 
 if TYPE_CHECKING:
     import linopy
@@ -100,6 +100,14 @@ def create_model(self, model: FlowSystemModel) -> ComponentModel:
         self.submodel = ComponentModel(model, self)
         return self.submodel
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested Interface objects and flows."""
+        super()._set_flow_system(flow_system)
+        if self.on_off_parameters is not None:
+            self.on_off_parameters._set_flow_system(flow_system)
+        for flow in self.inputs + self.outputs:
+            flow._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         if self.on_off_parameters is not None:
@@ -241,6 +249,12 @@ def create_model(self, model: FlowSystemModel) -> BusModel:
         self.submodel = BusModel(model, self)
         return self.submodel
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested flows."""
+        super()._set_flow_system(flow_system)
+        for flow in self.inputs + self.outputs:
+            flow._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         self.excess_penalty_per_flow_hour = self._fit_coords(
@@ -468,6 +482,14 @@ def create_model(self, model: FlowSystemModel) -> FlowModel:
         self.submodel = FlowModel(model, self)
         return self.submodel
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested Interface objects."""
+        super()._set_flow_system(flow_system)
+        if self.on_off_parameters is not None:
+            self.on_off_parameters._set_flow_system(flow_system)
+        if isinstance(self.size, Interface):
+            self.size._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         prefix = '|'.join(filter(None, [name_prefix, self.label_full]))
         self.relative_minimum = self._fit_coords(f'{prefix}|relative_minimum', self.relative_minimum)
diff --git a/flixopt/interface.py b/flixopt/interface.py
index e03738422..847aea96f 100644
--- a/flixopt/interface.py
+++ b/flixopt/interface.py
@@ -224,6 +224,12 @@ def __getitem__(self, index) -> Piece:
     def __iter__(self) -> Iterator[Piece]:
         return iter(self.pieces)  # Enables iteration like for piece in piecewise: ...
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested Piece objects."""
+        super()._set_flow_system(flow_system)
+        for piece in self.pieces:
+            piece._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         for i, piece in enumerate(self.pieces):
             piece.transform_data(f'{name_prefix}|Piece{i}')
@@ -450,6 +456,12 @@ def items(self):
         """
         return self.piecewises.items()
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested Piecewise objects."""
+        super()._set_flow_system(flow_system)
+        for piecewise in self.piecewises.values():
+            piecewise._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         for name, piecewise in self.piecewises.items():
             piecewise.transform_data(f'{name_prefix}|{name}')
@@ -662,6 +674,13 @@ def has_time_dim(self, value):
         for piecewise in self.piecewise_shares.values():
             piecewise.has_time_dim = value
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested Piecewise objects."""
+        super()._set_flow_system(flow_system)
+        self.piecewise_origin._set_flow_system(flow_system)
+        for piecewise in self.piecewise_shares.values():
+            piecewise._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         self.piecewise_origin.transform_data(f'{name_prefix}|PiecewiseEffects|origin')
         for effect, piecewise in self.piecewise_shares.items():
@@ -928,6 +947,12 @@ def __init__(
         self.maximum_size = maximum_size if maximum_size is not None else CONFIG.Modeling.big  # default maximum
         self.linked_periods = linked_periods
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to nested PiecewiseEffects object if present."""
+        super()._set_flow_system(flow_system)
+        if self.piecewise_effects_of_investment is not None:
+            self.piecewise_effects_of_investment._set_flow_system(flow_system)
+
     def transform_data(self, name_prefix: str = '') -> None:
         self.effects_of_investment = self._fit_effect_coords(
             prefix=name_prefix,
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 794608ea0..452094d9a 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -296,28 +296,14 @@ def _set_flow_system(self, flow_system: FlowSystem) -> None:
         elements to access FlowSystem properties without passing the reference
         through every method call.
 
+        Subclasses with nested Interface objects should override this method
+        to explicitly propagate the reference to their nested interfaces.
+
         Args:
             flow_system: The FlowSystem that this interface belongs to
         """
-        # Always set _flow_system (creates attribute if it doesn't exist)
         self._flow_system = flow_system
 
-        # Recursively set for nested Interface objects
-        for attr_name, attr_value in self.__dict__.items():
-            if attr_name.startswith('_'):
-                continue  # Skip private attributes
-
-            if isinstance(attr_value, Interface):
-                attr_value._set_flow_system(flow_system)
-            elif isinstance(attr_value, list):
-                for item in attr_value:
-                    if isinstance(item, Interface):
-                        item._set_flow_system(flow_system)
-            elif isinstance(attr_value, dict):
-                for item in attr_value.values():
-                    if isinstance(item, Interface):
-                        item._set_flow_system(flow_system)
-
     @property
     def flow_system(self) -> FlowSystem:
         """Access the FlowSystem this interface is linked to.

From ded18c841113b9cbf3ddc0ac585302229f105fef Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Sun, 16 Nov 2025 18:17:43 +0100
Subject: [PATCH 15/20] Add missing _set_flow_system

---
 flixopt/components.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/flixopt/components.py b/flixopt/components.py
index 2b091d7fa..1f7f4ab60 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -182,6 +182,12 @@ def create_model(self, model: FlowSystemModel) -> LinearConverterModel:
         self.submodel = LinearConverterModel(model, self)
         return self.submodel
 
+    def _set_flow_system(self, flow_system) -> None:
+        """Propagate flow_system reference to parent Component and piecewise_conversion."""
+        super()._set_flow_system(flow_system)
+        if self.piecewise_conversion is not None:
+            self.piecewise_conversion._set_flow_system(flow_system)
+
     def _plausibility_checks(self) -> None:
         super()._plausibility_checks()
         if not self.conversion_factors and not self.piecewise_conversion:

From 7c2a36b76204cc1c7c2ed6ef6366789f8d0092b2 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 17 Nov 2025 18:15:42 +0100
Subject: [PATCH 16/20] Add missing type hints

---
 flixopt/structure.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/flixopt/structure.py b/flixopt/structure.py
index 296dc8d0d..c768004be 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -24,7 +24,7 @@
 from loguru import logger
 
 from . import io as fx_io
-from .core import TimeSeriesData, get_dataarray_stats
+from .core import FlowSystemDimensions, TimeSeriesData, get_dataarray_stats
 
 if TYPE_CHECKING:  # for type checking and preventing circular imports
     import pathlib
@@ -32,6 +32,7 @@
 
     from .effects import EffectCollectionModel
     from .flow_system import FlowSystem
+    from .types import Effect_TPS, Numeric_TPS, NumericOrBool
 
 
 CLASS_REGISTRY = {}
@@ -320,7 +321,9 @@ def flow_system(self) -> FlowSystem:
             )
         return self._flow_system
 
-    def _fit_coords(self, name: str, data, dims=None):
+    def _fit_coords(
+        self, name: str, data: NumericOrBool | None, dims: Collection[FlowSystemDimensions] | None = None
+    ) -> xr.DataArray | None:
         """Convenience wrapper for FlowSystem.fit_to_model_coords().
 
         Args:
@@ -333,7 +336,13 @@ def _fit_coords(self, name: str, data, dims=None):
         """
         return self.flow_system.fit_to_model_coords(name, data, dims=dims)
 
-    def _fit_effect_coords(self, prefix: str, effect_values, suffix: str = None, dims=None):
+    def _fit_effect_coords(
+        self,
+        prefix: str | None,
+        effect_values: Effect_TPS | Numeric_TPS | None,
+        suffix: str | None = None,
+        dims: Collection[FlowSystemDimensions] | None = None,
+    ) -> Effect_TPS | None:
         """Convenience wrapper for FlowSystem.fit_effects_to_model_coords().
 
         Args:

From 72c525efdb7227999f4868ff1e286d79f46c0dba Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 17 Nov 2025 18:16:24 +0100
Subject: [PATCH 17/20] Change timing of validate_system_integrity() and
 improve cache invalidation

---
 flixopt/flow_system.py | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index b7725ed87..12970a23b 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -594,6 +594,10 @@ def connect_and_transform(self):
         self._connect_network()
         for element in chain(self.components.values(), self.effects.values(), self.buses.values()):
             element.transform_data()
+
+        # Validate cross-element references immediately after transformation
+        self._validate_system_integrity()
+
         self._connected_and_transformed = True
 
     def add_elements(self, *elements: Element) -> None:
@@ -645,8 +649,7 @@ def create_model(self, normalize_weights: bool = True) -> FlowSystemModel:
             raise RuntimeError(
                 'FlowSystem is not connected_and_transformed. Call FlowSystem.connect_and_transform() first.'
             )
-        # Validate cross-element references before creating model
-        self._validate_system_integrity()
+        # System integrity was already validated in connect_and_transform()
         self.model = FlowSystemModel(self, normalize_weights)
         return self.model
 
@@ -830,13 +833,17 @@ def _add_components(self, *components: Component) -> None:
         for new_component in list(components):
             new_component._set_flow_system(self)  # Link element to FlowSystem
             self.components.add(new_component)  # Add to existing components
-            self._flows_cache = None  # Invalidate flows cache
+        # Invalidate cache once after all additions
+        if components:
+            self._flows_cache = None
 
     def _add_buses(self, *buses: Bus):
         for new_bus in list(buses):
             new_bus._set_flow_system(self)  # Link element to FlowSystem
             self.buses.add(new_bus)  # Add to existing buses
-            self._flows_cache = None  # Invalidate flows cache
+        # Invalidate cache once after all additions
+        if buses:
+            self._flows_cache = None
 
     def _connect_network(self):
         """Connects the network of components and buses. Can be rerun without changes if no elements were added"""
@@ -867,9 +874,12 @@ def _connect_network(self):
                     bus.outputs.append(flow)
                 elif not flow.is_input_in_component and flow not in bus.inputs:
                     bus.inputs.append(flow)
+
+        # Count flows manually to avoid triggering cache rebuild
+        flow_count = sum(len(c.inputs) + len(c.outputs) for c in self.components.values())
         logger.debug(
             f'Connected {len(self.buses)} Buses and {len(self.components)} '
-            f'via {len(self.flows)} Flows inside the FlowSystem.'
+            f'via {flow_count} Flows inside the FlowSystem.'
         )
 
     def __repr__(self) -> str:

From afb7de364c450fd85b4a29342af890b71cfa27f2 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 17 Nov 2025 19:05:36 +0100
Subject: [PATCH 18/20]   1. calculation.py - Modeling State Consistency

  Updated the modeled property to use the _modeled flag instead of checking self.model is not None, ensuring consistent state tracking throughout
  the Calculation class.

  2. features.py - PieceModel Type Hints

  Improved the dims parameter type from FlowSystemDimensions | None to Collection[FlowSystemDimensions] | None in both PieceModel and
  PiecewiseModel to accurately reflect actual usage with tuples like ('period', 'scenario').

  3. features.py - ShareAllocationModel Error Message

  Updated the validation error message from "Both max_per_hour and min_per_hour cannot be used..." to "max_per_hour and min_per_hour require
  'time' dimension in dims" to match the actual condition being checked.

  4. features.py - PiecewiseModel Zero-Point Documentation

  Added clarifying comments explaining that the zero_point binary variable acts as a gate: when enabled (=1), at most one segment is active; when
  disabled (=0), all segments remain inactive.

  5. effects.py - Penalty Temporal Coupling

  Added a comment in EffectCollectionModel._do_modeling noting that penalty shares are added after the objective is set, explaining the temporal
  coupling in the design.

  6. structure.py - Interface Documentation

  Updated the class-level docstring to reflect the correct transform_data(name_prefix='') signature instead of the outdated
  transform_data(flow_system).

  7. components.py - Model Class Docstrings

  Replaced generic docstrings with specific descriptions:
  - TransmissionModel: "Create transmission efficiency equations and optional absolute loss constraints for both flow directions"
  - LinearConverterModel: "Create linear conversion equations or piecewise conversion constraints between input and output flows"
  - StorageModel: "Create charge state variables, energy balance equations, and optional investment submodels"

  All changes improve code clarity, consistency, and maintainability without altering functionality.
---
 flixopt/calculation.py |  6 +-----
 flixopt/components.py  |  6 +++---
 flixopt/effects.py     |  3 +++
 flixopt/features.py    | 12 +++++++++---
 flixopt/structure.py   |  2 +-
 5 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/flixopt/calculation.py b/flixopt/calculation.py
index 072e9204d..6fb4108db 100644
--- a/flixopt/calculation.py
+++ b/flixopt/calculation.py
@@ -100,8 +100,6 @@ def __init__(
             raise NotADirectoryError(f'Path {self.folder} exists and is not a directory.')
         self.folder.mkdir(parents=False, exist_ok=True)
 
-        self._modeled = False
-
     @property
     def main_results(self) -> dict[str, Scalar | dict]:
         from flixopt.features import InvestmentModel
@@ -198,7 +196,6 @@ def do_modeling(self) -> FullCalculation:
         self.model.do_modeling()
 
         self.durations['modeling'] = round(timeit.default_timer() - t_start, 2)
-        self._modeled = True
         return self
 
     def fix_sizes(self, ds: xr.Dataset, decimal_rounding: int | None = 5) -> FullCalculation:
@@ -232,7 +229,7 @@ def solve(
         self, solver: _Solver, log_file: pathlib.Path | None = None, log_main_results: bool | None = None
     ) -> FullCalculation:
         # Auto-call do_modeling() if not already done
-        if not self._modeled:
+        if not self.modeled:
             logger.info('Model not yet created. Calling do_modeling() automatically.')
             self.do_modeling()
 
@@ -334,7 +331,6 @@ def do_modeling(self) -> AggregatedCalculation:
         )
         self.aggregation_model.do_modeling()
         self.durations['modeling'] = round(timeit.default_timer() - t_start, 2)
-        self._modeled = True
         return self
 
     def _perform_aggregation(self):
diff --git a/flixopt/components.py b/flixopt/components.py
index 1f7f4ab60..2af5d427d 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -747,7 +747,7 @@ def __init__(self, model: FlowSystemModel, element: Transmission):
         super().__init__(model, element)
 
     def _do_modeling(self):
-        """Create variables, constraints, and nested submodels"""
+        """Create transmission efficiency equations and optional absolute loss constraints for both flow directions"""
         super()._do_modeling()
 
         # first direction
@@ -788,7 +788,7 @@ def __init__(self, model: FlowSystemModel, element: LinearConverter):
         super().__init__(model, element)
 
     def _do_modeling(self):
-        """Create variables, constraints, and nested submodels"""
+        """Create linear conversion equations or piecewise conversion constraints between input and output flows"""
         super()._do_modeling()
 
         # Create conversion factor constraints if specified
@@ -837,7 +837,7 @@ def __init__(self, model: FlowSystemModel, element: Storage):
         super().__init__(model, element)
 
     def _do_modeling(self):
-        """Create variables, constraints, and nested submodels"""
+        """Create charge state variables, energy balance equations, and optional investment submodels"""
         super()._do_modeling()
 
         # Create variables
diff --git a/flixopt/effects.py b/flixopt/effects.py
index 074809a5a..306f726c9 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -693,6 +693,9 @@ def _do_modeling(self):
         self._add_share_between_effects()
 
         # Set objective
+        # Note: penalty.total is used here, but penalty shares from buses/components
+        # are added later via add_share_to_penalty(). The ShareAllocationModel supports
+        # this pattern - shares can be added after the objective is defined.
         self._model.add_objective(
             (self.effects.objective_effect.submodel.total * self._model.weights).sum() + self.penalty.total.sum()
         )
diff --git a/flixopt/features.py b/flixopt/features.py
index d01c0c8ab..57ead98be 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -15,6 +15,8 @@
 from .structure import FlowSystemModel, Submodel
 
 if TYPE_CHECKING:
+    from collections.abc import Collection
+
     from .core import FlowSystemDimensions, Scalar, TemporalData
     from .interface import InvestParameters, OnOffParameters, Piecewise
 
@@ -329,7 +331,7 @@ def __init__(
         model: FlowSystemModel,
         label_of_element: str,
         label_of_model: str,
-        dims: FlowSystemDimensions | None,
+        dims: Collection[FlowSystemDimensions] | None,
     ):
         self.inside_piece: linopy.Variable | None = None
         self.lambda0: linopy.Variable | None = None
@@ -375,7 +377,7 @@ def __init__(
         label_of_model: str,
         piecewise_variables: dict[str, Piecewise],
         zero_point: bool | linopy.Variable | None,
-        dims: FlowSystemDimensions | None,
+        dims: Collection[FlowSystemDimensions] | None,
     ):
         """
         Modeling a Piecewise relation between miultiple variables.
@@ -451,6 +453,10 @@ def _do_modeling(self):
             else:
                 rhs = 1
 
+            # This constraint ensures at most one segment is active at a time.
+            # When zero_point is a binary variable, it acts as a gate:
+            # - zero_point=1: at most one segment can be active (normal piecewise operation)
+            # - zero_point=0: all segments must be inactive (effectively disables the piecewise)
             self.add_constraints(
                 sum([piece.inside_piece for piece in self.pieces]) <= rhs,
                 name=f'{self.label_full}|{variable.name}|single_segment',
@@ -536,7 +542,7 @@ def __init__(
         min_per_hour: TemporalData | None = None,
     ):
         if 'time' not in dims and (max_per_hour is not None or min_per_hour is not None):
-            raise ValueError('Both max_per_hour and min_per_hour cannot be used when has_time_dim is False')
+            raise ValueError("max_per_hour and min_per_hour require 'time' dimension in dims")
 
         self._dims = dims
         self.total_per_timestep: linopy.Variable | None = None
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 452094d9a..9d5f2fef6 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -270,7 +270,7 @@ class Interface:
         - Recursive handling of complex nested structures
 
     Subclasses must implement:
-        transform_data(flow_system): Transform data to match FlowSystem dimensions
+        transform_data(name_prefix=''): Transform data to match FlowSystem dimensions
     """
 
     def transform_data(self, name_prefix: str = '') -> None:

From ae616a62c41f904d49aee77d8720a704a01612fa Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Tue, 18 Nov 2025 18:51:48 +0100
Subject: [PATCH 19/20] Update CHANGELOG.md

---
 CHANGELOG.md | 34 +++++++++++++++++++---------------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7953efb5d..210fd9fb7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -51,31 +51,35 @@ If upgrading from v2.x, see the [v3.0.0 release notes](https://github.com/flixOp
 
 ## [Unreleased] - ????-??-??
 
-**Summary**:
+**Summary**: Internal architecture improvements to simplify FlowSystem-Element coupling and eliminate circular dependencies.
 
 If upgrading from v2.x, see the [v3.0.0 release notes](https://github.com/flixOpt/flixOpt/releases/tag/v3.0.0) and [Migration Guide](https://flixopt.github.io/flixopt/latest/user-guide/migration-guide-v3/).
 
 ### ✨ Added
-
-### 💥 Breaking Changes
+- **Auto-modeling**: `Calculation.solve()` now automatically calls `do_modeling()` if not already done, making the explicit `do_modeling()` call optional for simpler workflows
+- **System validation**: Added `_validate_system_integrity()` to validate cross-element references (e.g., Flow.bus) immediately after transformation, providing clearer error messages
+- **Element registration validation**: Added checks to prevent elements from being assigned to multiple FlowSystems simultaneously
+- **Helper methods in Interface base class**: Added `_fit_coords()` and `_fit_effect_coords()` convenience wrappers for cleaner data transformation code
+- **FlowSystem property in Interface**: Added `flow_system` property to access the linked FlowSystem with clear error messages if not yet linked
 
 ### ♻️ Changed
-
-### 🗑️ Deprecated
-
-### 🔥 Removed
+- **Refactored FlowSystem-Element coupling**:
+    - Introduced `_set_flow_system()` method in Interface base class to propagate FlowSystem reference to nested Interface objects
+    - Each Interface subclass now explicitly propagates the reference to its nested interfaces (e.g., Component → OnOffParameters, Flow → InvestParameters)
+    - Elements can now access FlowSystem via `self.flow_system` property instead of passing it through every method call
+- **Simplified transform_data() signature**: Removed `flow_system` parameter from `transform_data()` methods - FlowSystem reference is now accessed via `self.flow_system` property
+- **Two-phase modeling pattern**: Clarified the separation between variable creation (in `__init__`) and constraint creation (in `_do_modeling()`) to eliminate circular dependencies in Submodel architecture
+- **Improved cache invalidation**: Cache invalidation in `add_elements()` now happens once after all additions rather than per element
+- **Better logging**: Enhanced element registration logging to show element type and full label
 
 ### 🐛 Fixed
-
-### 🔒 Security
-
-### 📦 Dependencies
-
-### 📝 Docs
+- Fixed inconsistent argument passing in `_fit_effect_coords()` - changed parameters from `label_prefix`/`label_suffix` to `prefix`/`suffix` for consistency
 
 ### 👷 Development
-
-### 🚧 Known Issues
+- **Eliminated circular dependencies**: Implemented clean two-phase modeling pattern where Phase 1 creates all variables and Phase 2 creates constraints that reference those variables
+- Added comprehensive docstrings to `_do_modeling()` methods explaining "Create variables, constraints, and nested submodels"
+- Added missing type hints throughout the codebase
+- Improved code organization by making FlowSystem reference propagation explicit and traceable
 
 ---
 

From 354f068596430271dd4ab2effc76116ba38c4578 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Tue, 18 Nov 2025 18:55:39 +0100
Subject: [PATCH 20/20] Update CHANGELOG.md

---
 CHANGELOG.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 210fd9fb7..b5c73395b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -68,16 +68,16 @@ If upgrading from v2.x, see the [v3.0.0 release notes](https://github.com/flixOp
     - Each Interface subclass now explicitly propagates the reference to its nested interfaces (e.g., Component → OnOffParameters, Flow → InvestParameters)
     - Elements can now access FlowSystem via `self.flow_system` property instead of passing it through every method call
 - **Simplified transform_data() signature**: Removed `flow_system` parameter from `transform_data()` methods - FlowSystem reference is now accessed via `self.flow_system` property
-- **Two-phase modeling pattern**: Clarified the separation between variable creation (in `__init__`) and constraint creation (in `_do_modeling()`) to eliminate circular dependencies in Submodel architecture
+- **Two-phase modeling pattern within _do_modeling()**: Clarified the pattern where `_do_modeling()` creates nested submodels first (so their variables exist), then creates constraints that reference those variables - eliminates circular dependencies in Submodel architecture
 - **Improved cache invalidation**: Cache invalidation in `add_elements()` now happens once after all additions rather than per element
-- **Better logging**: Enhanced element registration logging to show element type and full label
+- **Better logging**: Centralized element registration logging to show element type and full label
 
 ### 🐛 Fixed
-- Fixed inconsistent argument passing in `_fit_effect_coords()` - changed parameters from `label_prefix`/`label_suffix` to `prefix`/`suffix` for consistency
+- Fixed inconsistent argument passing in `_fit_effect_coords()` - standardized all calls to use named arguments (`prefix=`, `effect_values=`, `suffix=`) instead of mix of positional and named arguments
 
 ### 👷 Development
-- **Eliminated circular dependencies**: Implemented clean two-phase modeling pattern where Phase 1 creates all variables and Phase 2 creates constraints that reference those variables
-- Added comprehensive docstrings to `_do_modeling()` methods explaining "Create variables, constraints, and nested submodels"
+- **Eliminated circular dependencies**: Implemented two-phase modeling pattern within `_do_modeling()` where nested submodels are created first (creating their variables), then constraints are created that can safely reference those submodel variables
+- Added comprehensive docstrings to `_do_modeling()` methods explaining the pattern: "Create variables, constraints, and nested submodels"
 - Added missing type hints throughout the codebase
 - Improved code organization by making FlowSystem reference propagation explicit and traceable