diff --git a/CHANGELOG.md b/CHANGELOG.md index 7953efb5d..b5c73395b 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 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**: Centralized element registration logging to show element type and full label ### 🐛 Fixed - -### 🔒 Security - -### 📦 Dependencies - -### 📝 Docs +- 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 - -### 🚧 Known Issues +- **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 --- diff --git a/flixopt/calculation.py b/flixopt/calculation.py index 64c589e3a..fcde018b7 100644 --- a/flixopt/calculation.py +++ b/flixopt/calculation.py @@ -98,8 +98,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, int | float | dict]: from flixopt.features import InvestmentModel @@ -228,6 +226,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( diff --git a/flixopt/components.py b/flixopt/components.py index c51b4b7d2..73e0c7972 100644 --- a/flixopt/components.py +++ b/flixopt/components.py @@ -181,6 +181,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: @@ -211,23 +217,23 @@ 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._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 @@ -433,46 +439,48 @@ 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 _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(flow_system, prefix) - self.relative_minimum_charge_state = flow_system.fit_to_model_coords( - f'{prefix}|relative_minimum_charge_state', - self.relative_minimum_charge_state, + super().transform_data(prefix) + self.relative_minimum_charge_state = self._fit_coords( + f'{prefix}|relative_minimum_charge_state', self.relative_minimum_charge_state ) - self.relative_maximum_charge_state = 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( - 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 = 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 = 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 = 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 = 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 = 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'], ) 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._fit_coords( f'{prefix}|capacity_in_flow_hours', self.capacity_in_flow_hours, dims=['period', 'scenario'] ) @@ -719,11 +727,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._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): @@ -738,7 +746,7 @@ def __init__(self, model: FlowSystemModel, element: Transmission): super().__init__(model, element) def _do_modeling(self): - """Initiates all FlowModels""" + """Create transmission efficiency equations and optional absolute loss constraints for both flow directions""" super()._do_modeling() # first direction @@ -779,8 +787,10 @@ def __init__(self, model: FlowSystemModel, element: LinearConverter): super().__init__(model, element) def _do_modeling(self): + """Create linear conversion equations or piecewise conversion constraints between input and output flows""" super()._do_modeling() - # 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) @@ -826,8 +836,10 @@ def __init__(self, model: FlowSystemModel, element: Storage): super().__init__(model, element) def _do_modeling(self): + """Create charge state variables, energy balance equations, and optional investment submodels""" super()._do_modeling() + # Create variables lb, ub = self._absolute_charge_state_bounds self.add_variables( lower=lb, @@ -838,6 +850,7 @@ def _do_modeling(self): self.add_variables(coords=self._model.get_coords(), short_name='netto_discharge') + # Create constraints (can now access flow.submodel.flow_rate) # netto_discharge: # eq: nettoFlow(t) - discharging(t) + charging(t) = 0 self.add_constraints( @@ -862,6 +875,7 @@ def _do_modeling(self): short_name='charge_state', ) + # Create InvestmentModel and bounding constraints for investment if isinstance(self.element.capacity_in_flow_hours, InvestParameters): self.add_submodels( InvestmentModel( @@ -883,6 +897,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 ebfc2c906..e011b185d 100644 --- a/flixopt/effects.py +++ b/flixopt/effects.py @@ -339,43 +339,40 @@ 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._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.maximum_per_hour = 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( - 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 = 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 = 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 = 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 = 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 = 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 = 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 = flow_system.fit_to_model_coords( + self.maximum_total = self._fit_coords( f'{prefix}|maximum_total', self.maximum_total, dims=['period', 'scenario'] ) @@ -396,6 +393,9 @@ def __init__(self, model: FlowSystemModel, element: Effect): super().__init__(model, element) 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( ShareAllocationModel( @@ -661,16 +661,26 @@ def add_share_to_penalty(self, name: str, expression: linopy.LinearExpression) - self.penalty.add_share(name, expression, dims=()) def _do_modeling(self): + """Create variables, constraints, and nested submodels""" super()._do_modeling() + + # 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', ) + # Add cross-effect shares 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/elements.py b/flixopt/elements.py index f47002b3a..8f88cbfbb 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 @@ -109,13 +109,21 @@ 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 _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: - 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] @@ -250,9 +258,15 @@ 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 _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 = 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 ) @@ -477,35 +491,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 _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 = 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( - f'{prefix}|fixed_relative_profile', self.fixed_relative_profile - ) - self.effects_per_flow_hour = 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.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 = 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 = 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 = 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'] ) 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._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 @@ -567,7 +585,9 @@ def __init__(self, model: FlowSystemModel, element: Flow): super().__init__(model, element) def _do_modeling(self): + """Create variables, constraints, and nested submodels""" super()._do_modeling() + # Main flow rate variable self.add_variables( lower=self.absolute_flow_rate_bounds[0], @@ -578,7 +598,7 @@ def _do_modeling(self): 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', @@ -623,6 +643,7 @@ 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 @@ -798,7 +819,8 @@ def __init__(self, model: FlowSystemModel, element: Bus): self.excess_output: linopy.Variable | None = None super().__init__(model, element) - def _do_modeling(self) -> None: + def _do_modeling(self): + """Create variables, constraints, and nested submodels""" super()._do_modeling() # inputs == outputs for flow in self.element.inputs + self.element.outputs: @@ -807,7 +829,7 @@ def _do_modeling(self) -> None: 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) @@ -845,9 +867,12 @@ def __init__(self, model: FlowSystemModel, element: Component): super().__init__(model, element) def _do_modeling(self): - """Initiates all FlowModels""" + """Create variables, constraints, and nested submodels""" super()._do_modeling() + 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: @@ -858,9 +883,11 @@ def _do_modeling(self): if flow.on_off_parameters is None: flow.on_off_parameters = OnOffParameters() + # 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) + # 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: diff --git a/flixopt/features.py b/flixopt/features.py index fd9796ba1..0ec935db6 100644 --- a/flixopt/features.py +++ b/flixopt/features.py @@ -14,6 +14,8 @@ from .structure import FlowSystemModel, Submodel if TYPE_CHECKING: + from collections.abc import Collection + from .core import FlowSystemDimensions from .interface import InvestParameters, OnOffParameters, Piecewise from .types import Numeric_PS, Numeric_TPS @@ -172,6 +174,7 @@ def __init__( super().__init__(model, label_of_element, label_of_model=label_of_model) def _do_modeling(self): + """Create variables, constraints, and nested submodels""" super()._do_modeling() if self.parameters.use_off: @@ -326,7 +329,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 @@ -336,7 +339,10 @@ def __init__( super().__init__(model, label_of_element, label_of_model) 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', @@ -356,6 +362,7 @@ def _do_modeling(self): coords=self._model.get_coords(dims=self.dims), ) + # Create constraints # eq: lambda0(t) + lambda1(t) = inside_piece(t) self.add_constraints(self.inside_piece == self.lambda0 + self.lambda1, short_name='inside_piece') @@ -368,7 +375,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. @@ -392,12 +399,15 @@ 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() + # 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 and constraints) for i in range(len(list(self._piecewise_variables.values())[0])): new_piece = self.add_submodels( PieceModel( @@ -441,6 +451,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', @@ -475,6 +489,10 @@ 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() + + # Create variables self.shares = { effect: self.add_variables(coords=self._model.get_coords(['period', 'scenario']), short_name=effect) for effect in self._piecewise_shares @@ -488,6 +506,7 @@ def _do_modeling(self): }, } + # Create piecewise model (which creates its variables and constraints) self.piecewise_model = self.add_submodels( PiecewiseModel( model=self._model, @@ -500,7 +519,7 @@ def _do_modeling(self): short_name='PiecewiseEffects', ) - # Shares + # 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()}, @@ -521,7 +540,7 @@ def __init__( min_per_hour: Numeric_TPS | 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 @@ -541,7 +560,10 @@ 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() + + # Create variables self.total = self.add_variables( lower=self._total_min if self._total_min is not None else -np.inf, upper=self._total_max if self._total_max is not None else np.inf, diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py index cf112a608..12970a23b 100644 --- a/flixopt/flow_system.py +++ b/flixopt/flow_system.py @@ -593,7 +593,11 @@ 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() + + # Validate cross-element references immediately after transformation + self._validate_system_integrity() + self._connected_and_transformed = True def add_elements(self, *elements: Element) -> None: @@ -610,17 +614,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: """ @@ -633,6 +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.' ) + # System integrity was already validated in connect_and_transform() self.model = FlowSystemModel(self, normalize_weights) return self.model @@ -766,22 +783,67 @@ 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. + + 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 + # Invalidate cache once after all additions + if components: + self._flows_cache = None 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 + # 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""" @@ -812,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: diff --git a/flixopt/interface.py b/flixopt/interface.py index f67f501ba..89e94fa55 100644 --- a/flixopt/interface.py +++ b/flixopt/interface.py @@ -74,10 +74,10 @@ def __init__(self, start: Numeric_TPS, end: Numeric_TPS): 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._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 @@ -220,9 +220,15 @@ 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 _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(flow_system, f'{name_prefix}|Piece{i}') + piece.transform_data(f'{name_prefix}|Piece{i}') @register_class_for_io @@ -446,9 +452,15 @@ def items(self): """ return self.piecewises.items() - def transform_data(self, flow_system: FlowSystem, name_prefix: str = '') -> None: + 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(flow_system, f'{name_prefix}|{name}') + piecewise.transform_data(f'{name_prefix}|{name}') @register_class_for_io @@ -658,10 +670,17 @@ 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 _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(): - piecewise.transform_data(flow_system, f'{name_prefix}|PiecewiseEffects|{effect}') + piecewise.transform_data(f'{name_prefix}|PiecewiseEffects|{effect}') @register_class_for_io @@ -920,34 +939,40 @@ 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( - label_prefix=name_prefix, + 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, effect_values=self.effects_of_investment, - label_suffix='effects_of_investment', + suffix='effects_of_investment', dims=['period', 'scenario'], ) - self.effects_of_retirement = 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 = 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'], ) 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._fit_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._fit_coords( f'{name_prefix}|maximum_size', self.maximum_size, dims=['period', 'scenario'] ) # Convert tuple (first_period, last_period) to DataArray if needed @@ -956,30 +981,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._fit_coords( f'{name_prefix}|linked_periods', self.linked_periods, dims=['period', 'scenario'] ) - self.fixed_size = 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: @@ -1282,32 +1305,36 @@ def __init__( self.switch_on_total_max = 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( - name_prefix, self.effects_per_switch_on, 'per_switch_on' + def transform_data(self, name_prefix: str = '') -> None: + 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 = flow_system.fit_effects_to_model_coords( - name_prefix, self.effects_per_running_hour, 'per_running_hour' + self.effects_per_running_hour = self._fit_effect_coords( + prefix=name_prefix, + effect_values=self.effects_per_running_hour, + suffix='per_running_hour', ) - self.consecutive_on_hours_min = 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 = 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 = 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 = 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 = 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 = 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 = 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 9ddf46d31..21383913b 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 = {} @@ -95,6 +96,7 @@ def __init__(self, flow_system: FlowSystem, normalize_weights: bool): self.submodels: Submodels = Submodels({}) def do_modeling(self): + # 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) @@ -264,21 +266,96 @@ 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, 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. + + 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 + """ + self._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: NumericOrBool | None, dims: Collection[FlowSystemDimensions] | None = None + ) -> xr.DataArray | 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 | 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: + 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. @@ -856,6 +933,7 @@ 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 def _plausibility_checks(self) -> None: """This function is used to do some basic plausibility checks for each Element during initialization. @@ -1417,7 +1495,12 @@ def hours_per_step(self): return self._model.hours_per_step def _do_modeling(self): - """Called at the end of initialization. Override in subclasses to create variables and constraints.""" + """ + Override in subclasses to create variables, constraints, and submodels. + + This method is called during __init__. Create all nested submodels first + (so their variables exist), then create constraints that reference those variables. + """ pass