diff --git a/CHANGELOG.md b/CHANGELOG.md
index 97486fb78..40e12ec42 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -51,6 +51,8 @@ This release introduces new model dimensions (periods and scenarios) for multi-p
 
 - **Period dimension**: Enables multi-period investment modeling with distinct decisions in each period for transformation pathway optimization
 - **Scenario dimension**: Supports stochastic modeling with weighted scenarios for robust decision-making under uncertainty (demand, prices, weather)
+    - Control variable independence across scenarios via `scenario_independent_sizes` and `scenario_independent_flow_rates` parameters
+    - By default, investment sizes are shared across scenarios while flow rates vary per scenario
 
 **Redesigned effect sharing system:**
 
diff --git a/docs/user-guide/mathematical-notation/dimensions.md b/docs/user-guide/mathematical-notation/dimensions.md
index e7104fa6a..d1bc99c8e 100644
--- a/docs/user-guide/mathematical-notation/dimensions.md
+++ b/docs/user-guide/mathematical-notation/dimensions.md
@@ -132,9 +132,9 @@ $$
 
 ### Shared Periodic Decisions: The Exception
 
-**Within a period, periodic (investment) decisions are shared across all scenarios:**
+**Investment decisions (sizes) can be shared across all scenarios:**
 
-If a period has multiple scenarios, periodic variables (e.g., component size) are **scenario-independent** but **temporal variables are scenario-specific**.
+By default, sizes (e.g., Storage capacity, Thermal power, ...) are **scenario-independent** but **flow_rates are scenario-specific**.
 
 **Example - Flow with investment:**
 
@@ -153,6 +153,28 @@ $$
 
 This reflects real-world investment under uncertainty: you build capacity once (periodic/investment decision), but it operates under different conditions (temporal/operational decisions per scenario).
 
+**Mathematical Flexibility:**
+
+Variables can be either scenario-independent or scenario-specific:
+
+| Variable Type | Scenario-Independent | Scenario-Specific |
+|---------------|---------------------|-------------------|
+| **Sizes** (e.g., $\text{P}$) | $\text{P}(y)$ - Single value per period | $\text{P}(y, s)$ - Different per scenario |
+| **Flow rates** (e.g., $p(\text{t}_i)$) | $p(\text{t}_i, y)$ - Same across scenarios | $p(\text{t}_i, y, s)$ - Different per scenario |
+
+**Use Cases:**
+
+*Investment problems (with InvestParameters):*
+- **Sizes shared** (default): Investment under uncertainty - build capacity that performs well across all scenarios
+- **Sizes vary**: Scenario-specific capacity planning where different investments can be made for each future
+- **Selected sizes shared**: Mix of shared critical infrastructure and scenario-specific optional/flexible capacity
+
+*Dispatch problems (fixed sizes, no investments):*
+- **Flow rates shared**: Robust dispatch - find a single operational strategy that works across all forecast scenarios (e.g., day-ahead unit commitment under demand/weather uncertainty)
+- **Flow rates vary** (default): Scenario-adaptive dispatch - optimize operations for each scenario's specific conditions (demand, weather, prices)
+
+For implementation details on controlling scenario independence, see the [`FlowSystem`][flixopt.flow_system.FlowSystem] API reference.
+
 ---
 
 ## Dimensional Impact on Objective Function
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 769391fb7..ad43c183b 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -60,6 +60,14 @@ class FlowSystem(Interface):
             If you use an array, take care that its long enough to cover all previous values!
         weights: The weights of each period and scenario. If None, all scenarios have the same weight (normalized to 1).
             Its recommended to normalize the weights to sum up to 1.
+        scenario_independent_sizes: Controls whether investment sizes are equalized across scenarios.
+            - True: All sizes are shared/equalized across scenarios
+            - False: All sizes are optimized separately per scenario
+            - list[str]: Only specified components (by label_full) are equalized across scenarios
+        scenario_independent_flow_rates: Controls whether flow rates are equalized across scenarios.
+            - True: All flow rates are shared/equalized across scenarios
+            - False: All flow rates are optimized separately per scenario
+            - list[str]: Only specified flows (by label_full) are equalized across scenarios
 
     Notes:
         - Creates an empty registry for components and buses, an empty EffectCollection, and a placeholder for a SystemModel.
@@ -75,6 +83,8 @@ def __init__(
         hours_of_last_timestep: float | None = None,
         hours_of_previous_timesteps: int | float | np.ndarray | None = None,
         weights: PeriodicDataUser | None = None,
+        scenario_independent_sizes: bool | list[str] = True,
+        scenario_independent_flow_rates: bool | list[str] = False,
     ):
         self.timesteps = self._validate_timesteps(timesteps)
         self.timesteps_extra = self._create_timesteps_with_extra(self.timesteps, hours_of_last_timestep)
@@ -104,6 +114,10 @@ def __init__(
 
         self._network_app = None
 
+        # Use properties to validate and store scenario dimension settings
+        self.scenario_independent_sizes = scenario_independent_sizes
+        self.scenario_independent_flow_rates = scenario_independent_flow_rates
+
     @staticmethod
     def _validate_timesteps(timesteps: pd.DatetimeIndex) -> pd.DatetimeIndex:
         """Validate timesteps format and rename if needed."""
@@ -268,6 +282,8 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
             else None,
             hours_of_last_timestep=reference_structure.get('hours_of_last_timestep'),
             hours_of_previous_timesteps=reference_structure.get('hours_of_previous_timesteps'),
+            scenario_independent_sizes=reference_structure.get('scenario_independent_sizes', True),
+            scenario_independent_flow_rates=reference_structure.get('scenario_independent_flow_rates', False),
         )
 
         # Restore components
@@ -762,6 +778,77 @@ def coords(self) -> dict[FlowSystemDimensions, pd.Index]:
     def used_in_calculation(self) -> bool:
         return self._used_in_calculation
 
+    def _validate_scenario_parameter(self, value: bool | list[str], param_name: str, element_type: str) -> None:
+        """
+        Validate scenario parameter value.
+
+        Args:
+            value: The value to validate
+            param_name: Name of the parameter (for error messages)
+            element_type: Type of elements expected in list (e.g., 'component label_full', 'flow label_full')
+
+        Raises:
+            TypeError: If value is not bool or list[str]
+            ValueError: If list contains non-string elements
+        """
+        if isinstance(value, bool):
+            return  # Valid
+        elif isinstance(value, list):
+            if not all(isinstance(item, str) for item in value):
+                raise ValueError(f'{param_name} list must contain only strings ({element_type} values)')
+        else:
+            raise TypeError(f'{param_name} must be bool or list[str], got {type(value).__name__}')
+
+    @property
+    def scenario_independent_sizes(self) -> bool | list[str]:
+        """
+        Controls whether investment sizes are equalized across scenarios.
+
+        Returns:
+            bool or list[str]: Configuration for scenario-independent sizing
+        """
+        return self._scenario_independent_sizes
+
+    @scenario_independent_sizes.setter
+    def scenario_independent_sizes(self, value: bool | list[str]) -> None:
+        """
+        Set whether investment sizes should be equalized across scenarios.
+
+        Args:
+            value: True (all equalized), False (all vary), or list of component label_full strings to equalize
+
+        Raises:
+            TypeError: If value is not bool or list[str]
+            ValueError: If list contains non-string elements
+        """
+        self._validate_scenario_parameter(value, 'scenario_independent_sizes', 'Element.label_full')
+        self._scenario_independent_sizes = value
+
+    @property
+    def scenario_independent_flow_rates(self) -> bool | list[str]:
+        """
+        Controls whether flow rates are equalized across scenarios.
+
+        Returns:
+            bool or list[str]: Configuration for scenario-independent flow rates
+        """
+        return self._scenario_independent_flow_rates
+
+    @scenario_independent_flow_rates.setter
+    def scenario_independent_flow_rates(self, value: bool | list[str]) -> None:
+        """
+        Set whether flow rates should be equalized across scenarios.
+
+        Args:
+            value: True (all equalized), False (all vary), or list of flow label_full strings to equalize
+
+        Raises:
+            TypeError: If value is not bool or list[str]
+            ValueError: If list contains non-string elements
+        """
+        self._validate_scenario_parameter(value, 'scenario_independent_flow_rates', 'Flow.label_full')
+        self._scenario_independent_flow_rates = value
+
     def sel(
         self,
         time: str | slice | list[str] | pd.Timestamp | pd.DatetimeIndex | None = None,
diff --git a/flixopt/structure.py b/flixopt/structure.py
index defccd4b8..72efc3df2 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -103,6 +103,55 @@ def do_modeling(self):
         for bus in self.flow_system.buses.values():
             bus.create_model(self)
 
+        # Add scenario equality constraints after all elements are modeled
+        self._add_scenario_equality_constraints()
+
+    def _add_scenario_equality_for_parameter_type(
+        self,
+        parameter_type: Literal['flow_rate', 'size'],
+        config: bool | list[str],
+    ):
+        """Add scenario equality constraints for a specific parameter type.
+
+        Args:
+            parameter_type: The type of parameter ('flow_rate' or 'size')
+            config: Configuration value (True = equalize all, False = equalize none, list = equalize these)
+        """
+        if config is False:
+            return  # All vary per scenario, no constraints needed
+
+        suffix = f'|{parameter_type}'
+        if config is True:
+            # All should be scenario-independent
+            vars_to_constrain = [var for var in self.variables if var.endswith(suffix)]
+        else:
+            # Only those in the list should be scenario-independent
+            all_vars = [var for var in self.variables if var.endswith(suffix)]
+            to_equalize = {f'{element}{suffix}' for element in config}
+            vars_to_constrain = [var for var in all_vars if var in to_equalize]
+
+        # Validate that all specified variables exist
+        missing_vars = [v for v in vars_to_constrain if v not in self.variables]
+        if missing_vars:
+            param_name = 'scenario_independent_sizes' if parameter_type == 'size' else 'scenario_independent_flow_rates'
+            raise ValueError(f'{param_name} contains invalid labels: {missing_vars}')
+
+        logger.debug(f'Adding scenario equality constraints for {len(vars_to_constrain)} {parameter_type} variables')
+        for var in vars_to_constrain:
+            self.add_constraints(
+                self.variables[var].isel(scenario=0) == self.variables[var].isel(scenario=slice(1, None)),
+                name=f'{var}|scenario_independent',
+            )
+
+    def _add_scenario_equality_constraints(self):
+        """Add equality constraints to equalize variables across scenarios based on FlowSystem configuration."""
+        # Only proceed if we have scenarios
+        if self.flow_system.scenarios is None or len(self.flow_system.scenarios) <= 1:
+            return
+
+        self._add_scenario_equality_for_parameter_type('flow_rate', self.flow_system.scenario_independent_flow_rates)
+        self._add_scenario_equality_for_parameter_type('size', self.flow_system.scenario_independent_sizes)
+
     @property
     def solution(self):
         solution = super().solution
diff --git a/tests/test_scenarios.py b/tests/test_scenarios.py
index 1ff9e9cea..3f0637c91 100644
--- a/tests/test_scenarios.py
+++ b/tests/test_scenarios.py
@@ -337,3 +337,356 @@ def test_scenarios_selection(flow_system_piecewise_conversion_scenarios):
     )  ## Account for rounding errors
 
     assert calc.results.solution.indexes['scenario'].equals(flow_system_full.scenarios[0:2])
+
+
+def test_sizes_per_scenario_default():
+    """Test that scenario_independent_sizes defaults to True (sizes equalized) and flow_rates to False (vary)."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios)
+
+    assert fs.scenario_independent_sizes is True
+    assert fs.scenario_independent_flow_rates is False
+
+
+def test_sizes_per_scenario_bool():
+    """Test scenario_independent_sizes with boolean values."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    # Test False (vary per scenario)
+    fs1 = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios, scenario_independent_sizes=False)
+    assert fs1.scenario_independent_sizes is False
+
+    # Test True (equalized across scenarios)
+    fs2 = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios, scenario_independent_sizes=True)
+    assert fs2.scenario_independent_sizes is True
+
+
+def test_sizes_per_scenario_list():
+    """Test scenario_independent_sizes with list of element labels."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=['solar->grid', 'battery->grid'],
+    )
+
+    assert fs.scenario_independent_sizes == ['solar->grid', 'battery->grid']
+
+
+def test_flow_rates_per_scenario_default():
+    """Test that scenario_independent_flow_rates defaults to False (flow rates vary by scenario)."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios)
+
+    assert fs.scenario_independent_flow_rates is False
+
+
+def test_flow_rates_per_scenario_bool():
+    """Test scenario_independent_flow_rates with boolean values."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    # Test False (vary per scenario)
+    fs1 = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios, scenario_independent_flow_rates=False)
+    assert fs1.scenario_independent_flow_rates is False
+
+    # Test True (equalized across scenarios)
+    fs2 = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios, scenario_independent_flow_rates=True)
+    assert fs2.scenario_independent_flow_rates is True
+
+
+def test_scenario_parameters_property_setters():
+    """Test that scenario parameters can be changed via property setters."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios)
+
+    # Change scenario_independent_sizes
+    fs.scenario_independent_sizes = True
+    assert fs.scenario_independent_sizes is True
+
+    fs.scenario_independent_sizes = ['component1', 'component2']
+    assert fs.scenario_independent_sizes == ['component1', 'component2']
+
+    # Change scenario_independent_flow_rates
+    fs.scenario_independent_flow_rates = True
+    assert fs.scenario_independent_flow_rates is True
+
+    fs.scenario_independent_flow_rates = ['flow1', 'flow2']
+    assert fs.scenario_independent_flow_rates == ['flow1', 'flow2']
+
+
+def test_scenario_parameters_validation():
+    """Test that scenario parameters are validated correctly."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(timesteps=timesteps, scenarios=scenarios)
+
+    # Test invalid type
+    with pytest.raises(TypeError, match='must be bool or list'):
+        fs.scenario_independent_sizes = 'invalid'
+
+    # Test invalid list content
+    with pytest.raises(ValueError, match='must contain only strings'):
+        fs.scenario_independent_sizes = [1, 2, 3]
+
+
+def test_size_equality_constraints():
+    """Test that size equality constraints are created when scenario_independent_sizes=True."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=True,  # Sizes should be equalized
+        scenario_independent_flow_rates=False,  # Flow rates can vary
+    )
+
+    bus = fx.Bus('grid')
+    source = fx.Source(
+        label='solar',
+        outputs=[
+            fx.Flow(
+                label='out',
+                bus='grid',
+                size=fx.InvestParameters(
+                    minimum_size=10,
+                    maximum_size=100,
+                    effects_of_investment_per_size={'cost': 100},
+                ),
+            )
+        ],
+    )
+
+    fs.add_elements(bus, source, fx.Effect('cost', 'Total cost', '€', is_objective=True))
+
+    calc = fx.FullCalculation('test', fs)
+    calc.do_modeling()
+
+    # Check that size equality constraint exists
+    constraint_names = [str(c) for c in calc.model.constraints]
+    size_constraints = [c for c in constraint_names if 'scenario_independent' in c and 'size' in c]
+
+    assert len(size_constraints) > 0, 'Size equality constraint should exist'
+
+
+def test_flow_rate_equality_constraints():
+    """Test that flow_rate equality constraints are created when scenario_independent_flow_rates=True."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=False,  # Sizes can vary
+        scenario_independent_flow_rates=True,  # Flow rates should be equalized
+    )
+
+    bus = fx.Bus('grid')
+    source = fx.Source(
+        label='solar',
+        outputs=[
+            fx.Flow(
+                label='out',
+                bus='grid',
+                size=fx.InvestParameters(
+                    minimum_size=10,
+                    maximum_size=100,
+                    effects_of_investment_per_size={'cost': 100},
+                ),
+            )
+        ],
+    )
+
+    fs.add_elements(bus, source, fx.Effect('cost', 'Total cost', '€', is_objective=True))
+
+    calc = fx.FullCalculation('test', fs)
+    calc.do_modeling()
+
+    # Check that flow_rate equality constraint exists
+    constraint_names = [str(c) for c in calc.model.constraints]
+    flow_rate_constraints = [c for c in constraint_names if 'scenario_independent' in c and 'flow_rate' in c]
+
+    assert len(flow_rate_constraints) > 0, 'Flow rate equality constraint should exist'
+
+
+def test_selective_scenario_independence():
+    """Test selective scenario independence with specific element lists."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=['solar(out)'],  # Only solar size is equalized
+        scenario_independent_flow_rates=['demand(in)'],  # Only demand flow_rate is equalized
+    )
+
+    bus = fx.Bus('grid')
+    source = fx.Source(
+        label='solar',
+        outputs=[
+            fx.Flow(
+                label='out',
+                bus='grid',
+                size=fx.InvestParameters(
+                    minimum_size=10, maximum_size=100, effects_of_investment_per_size={'cost': 100}
+                ),
+            )
+        ],
+    )
+    sink = fx.Sink(
+        label='demand',
+        inputs=[fx.Flow(label='in', bus='grid', size=50)],
+    )
+
+    fs.add_elements(bus, source, sink, fx.Effect('cost', 'Total cost', '€', is_objective=True))
+
+    calc = fx.FullCalculation('test', fs)
+    calc.do_modeling()
+
+    constraint_names = [str(c) for c in calc.model.constraints]
+
+    # Solar SHOULD have size constraints (it's in the list, so equalized)
+    solar_size_constraints = [c for c in constraint_names if 'solar(out)|size' in c and 'scenario_independent' in c]
+    assert len(solar_size_constraints) > 0
+
+    # Solar should NOT have flow_rate constraints (not in the list, so varies per scenario)
+    solar_flow_constraints = [
+        c for c in constraint_names if 'solar(out)|flow_rate' in c and 'scenario_independent' in c
+    ]
+    assert len(solar_flow_constraints) == 0
+
+    # Demand should NOT have size constraints (no InvestParameters, size is fixed)
+    demand_size_constraints = [c for c in constraint_names if 'demand(in)|size' in c and 'scenario_independent' in c]
+    assert len(demand_size_constraints) == 0
+
+    # Demand SHOULD have flow_rate constraints (it's in the list, so equalized)
+    demand_flow_constraints = [
+        c for c in constraint_names if 'demand(in)|flow_rate' in c and 'scenario_independent' in c
+    ]
+    assert len(demand_flow_constraints) > 0
+
+
+def test_scenario_parameters_io_persistence():
+    """Test that scenario_independent_sizes and scenario_independent_flow_rates persist through IO operations."""
+    import shutil
+    import tempfile
+
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    # Create FlowSystem with custom scenario parameters
+    fs_original = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=['solar(out)'],
+        scenario_independent_flow_rates=True,
+    )
+
+    bus = fx.Bus('grid')
+    source = fx.Source(
+        label='solar',
+        outputs=[
+            fx.Flow(
+                label='out',
+                bus='grid',
+                size=fx.InvestParameters(
+                    minimum_size=10, maximum_size=100, effects_of_investment_per_size={'cost': 100}
+                ),
+            )
+        ],
+    )
+
+    fs_original.add_elements(bus, source, fx.Effect('cost', 'Total cost', '€', is_objective=True))
+
+    # Save to dataset
+    fs_original.connect_and_transform()
+    ds = fs_original.to_dataset()
+
+    # Load from dataset
+    fs_loaded = fx.FlowSystem.from_dataset(ds)
+
+    # Verify parameters persisted
+    assert fs_loaded.scenario_independent_sizes == fs_original.scenario_independent_sizes
+    assert fs_loaded.scenario_independent_flow_rates == fs_original.scenario_independent_flow_rates
+
+
+def test_scenario_parameters_io_with_calculation():
+    """Test that scenario parameters persist through full calculation IO."""
+    import shutil
+    import tempfile
+
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['base', 'high'], name='scenario')
+
+    fs = fx.FlowSystem(
+        timesteps=timesteps,
+        scenarios=scenarios,
+        scenario_independent_sizes=True,
+        scenario_independent_flow_rates=['demand(in)'],
+    )
+
+    bus = fx.Bus('grid')
+    source = fx.Source(
+        label='solar',
+        outputs=[
+            fx.Flow(
+                label='out',
+                bus='grid',
+                size=fx.InvestParameters(
+                    minimum_size=10, maximum_size=100, effects_of_investment_per_size={'cost': 100}
+                ),
+            )
+        ],
+    )
+    sink = fx.Sink(
+        label='demand',
+        inputs=[fx.Flow(label='in', bus='grid', size=50)],
+    )
+
+    fs.add_elements(bus, source, sink, fx.Effect('cost', 'Total cost', '€', is_objective=True))
+
+    # Create temp directory for results
+    temp_dir = tempfile.mkdtemp()
+
+    try:
+        # Solve and save
+        calc = fx.FullCalculation('test_io', fs, folder=temp_dir)
+        calc.do_modeling()
+        calc.solve(fx.solvers.HighsSolver(mip_gap=0.01, time_limit_seconds=60))
+        calc.results.to_file()
+
+        # Load results
+        results = fx.results.CalculationResults.from_file(temp_dir, 'test_io')
+        fs_loaded = fx.FlowSystem.from_dataset(results.flow_system_data)
+
+        # Verify parameters persisted
+        assert fs_loaded.scenario_independent_sizes == fs.scenario_independent_sizes
+        assert fs_loaded.scenario_independent_flow_rates == fs.scenario_independent_flow_rates
+
+        # Verify constraints are recreated correctly
+        calc2 = fx.FullCalculation('test_io_2', fs_loaded, folder=temp_dir)
+        calc2.do_modeling()
+
+        constraint_names1 = [str(c) for c in calc.model.constraints]
+        constraint_names2 = [str(c) for c in calc2.model.constraints]
+
+        size_constraints1 = [c for c in constraint_names1 if 'scenario_independent' in c and 'size' in c]
+        size_constraints2 = [c for c in constraint_names2 if 'scenario_independent' in c and 'size' in c]
+
+        assert len(size_constraints1) == len(size_constraints2)
+
+    finally:
+        # Clean up
+        shutil.rmtree(temp_dir)