Add mathematical correctness test suite

CHANGELOG.md

@@ -52,6 +52,33 @@ If upgrading from v2.x, see the [v3.0.0 release notes](https://github.com/flixOp
 
 Until here -->
 
+## [6.0.3] - Upcoming
+
+**Summary**: Bugfix release fixing `cluster_weight` loss during NetCDF roundtrip for manually constructed clustered FlowSystems.
+
+### 🐛 Fixed
+
+- **Clustering IO**: `cluster_weight` is now preserved during NetCDF roundtrip for manually constructed clustered FlowSystems (i.e. `FlowSystem(..., clusters=..., cluster_weight=...)`). Previously, `cluster_weight` was silently dropped to `None` during `save->reload->solve`, causing incorrect objective values. Systems created via `.transform.cluster()` were not affected.
+
+### 👷 Development
+
+- **New `test_math/` test suite**: Comprehensive mathematical correctness tests with exact, hand-calculated assertions. Each test runs in 3 IO modes (solve, save→reload→solve, solve→save→reload) via the `optimize` fixture:
+    - `test_flow.py` — flow bounds, merit order, relative min/max, on/off hours
+    - `test_flow_invest.py` — investment sizing, fixed-size, optional invest, piecewise invest
+    - `test_flow_status.py` — startup costs, switch-on/off constraints, status penalties
+    - `test_bus.py` — bus balance, excess/shortage penalties
+    - `test_effects.py` — effect aggregation, periodic/temporal effects, multi-effect objectives
+    - `test_components.py` — SourceAndSink, converters, links, combined heat-and-power
+    - `test_conversion.py` — linear converter balance, multi-input/output, efficiency
+    - `test_piecewise.py` — piecewise-linear efficiency, segment selection
+    - `test_storage.py` — charge/discharge, SOC tracking, final charge state, losses
+    - `test_multi_period.py` — period weights, invest across periods
+    - `test_scenarios.py` — scenario weights, scenario-independent flows
+    - `test_clustering.py` — exact per-timestep flow_rates, effects, and charge_state in clustered systems (incl. non-equal cluster weights to cover IO roundtrip)
+    - `test_validation.py` — plausibility checks and error messages
+
+---
+
 ## [6.0.2] - 2026-02-05
 
 **Summary**: Patch release which improves `Comparison` coordinate handling.

docs/notebooks/08d-clustering-multiperiod.ipynb

@@ -525,7 +525,7 @@
    "id": "29",
    "metadata": {},
    "source": [
-    "<cell_type>markdown</cell_type>## Summary\n",
+    "## Summary\n",
     "\n",
     "You learned how to:\n",
     "\n",

docs/notebooks/08f-clustering-segmentation.ipynb

@@ -527,7 +527,7 @@
    "id": "33",
    "metadata": {},
    "source": [
-    "<cell_type>markdown</cell_type>## API Reference\n",
+    "## API Reference\n",
     "\n",
     "### SegmentConfig Parameters\n",
     "\n",

flixopt/comparison.py

@@ -53,6 +53,14 @@ def _extract_nonindex_coords(datasets: list[xr.Dataset]) -> tuple[list[xr.Datase
             coords_to_drop.add(name)
             if name not in merged:
                 merged[name] = (dim, {})
+            elif merged[name][0] != dim:
+                warnings.warn(
+                    f"Coordinate '{name}' appears on different dims: "
+                    f"'{merged[name][0]}' vs '{dim}'. Dropping this coordinate.",
+                    stacklevel=4,
+                )
+                del merged[name]
+                continue
 
             for dv, cv in zip(ds.coords[dim].values, coord.values, strict=False):
                 if dv not in merged[name][1]:

flixopt/components.py

@@ -1144,8 +1144,13 @@ def _relative_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
             min_final_da = min_final_da.assign_coords(time=[timesteps_extra[-1]])
             min_bounds = xr.concat([rel_min, min_final_da], dim='time')
         else:
-            # Original is scalar - broadcast to full time range (constant value)
-            min_bounds = rel_min.expand_dims(time=timesteps_extra)
+            # Original is scalar - expand to regular timesteps, then concat with final value
+            regular_min = rel_min.expand_dims(time=timesteps_extra[:-1])
+            min_final_da = (
+                min_final_value.expand_dims('time') if 'time' not in min_final_value.dims else min_final_value
+            )
+            min_final_da = min_final_da.assign_coords(time=[timesteps_extra[-1]])
+            min_bounds = xr.concat([regular_min, min_final_da], dim='time')
 
         if 'time' in rel_max.dims:
             # Original has time dim - concat with final value
@@ -1155,8 +1160,13 @@ def _relative_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
             max_final_da = max_final_da.assign_coords(time=[timesteps_extra[-1]])
             max_bounds = xr.concat([rel_max, max_final_da], dim='time')
         else:
-            # Original is scalar - broadcast to full time range (constant value)
-            max_bounds = rel_max.expand_dims(time=timesteps_extra)
+            # Original is scalar - expand to regular timesteps, then concat with final value
+            regular_max = rel_max.expand_dims(time=timesteps_extra[:-1])
+            max_final_da = (
+                max_final_value.expand_dims('time') if 'time' not in max_final_value.dims else max_final_value
+            )
+            max_final_da = max_final_da.assign_coords(time=[timesteps_extra[-1]])
+            max_bounds = xr.concat([regular_max, max_final_da], dim='time')
 
         # Ensure both bounds have matching dimensions (broadcast once here,
         # so downstream code doesn't need to handle dimension mismatches)

flixopt/io.py

@@ -1631,15 +1631,12 @@ def _create_flow_system(
         # Extract cluster index if present (clustered FlowSystem)
         clusters = ds.indexes.get('cluster')
 
-        # For clustered datasets, cluster_weight is (cluster,) shaped - set separately
-        if clusters is not None:
-            cluster_weight_for_constructor = None
-        else:
-            cluster_weight_for_constructor = (
-                cls._resolve_dataarray_reference(reference_structure['cluster_weight'], arrays_dict)
-                if 'cluster_weight' in reference_structure
-                else None
-            )
+        # Resolve cluster_weight if present in reference structure
+        cluster_weight_for_constructor = (
+            cls._resolve_dataarray_reference(reference_structure['cluster_weight'], arrays_dict)
+            if 'cluster_weight' in reference_structure
+            else None
+        )
 
         # Resolve scenario_weights only if scenario dimension exists
         scenario_weights = None

pyproject.toml

@@ -176,7 +176,7 @@ extend-fixable = ["B"]   # Enable fix for flake8-bugbear (`B`), on top of any ru
 # Apply rule exceptions to specific files or directories
 [tool.ruff.lint.per-file-ignores]
 "tests/*.py" = ["S101"]  # Ignore assertions in test files
-"tests/test_integration.py" = ["N806"]  # Ignore NOT lowercase names in test files
+"tests/superseded/test_integration.py" = ["N806"]  # Ignore NOT lowercase names in test files
 "flixopt/linear_converters.py" = ["N803"]  # Parameters with NOT lowercase names
 
 [tool.ruff.format]
@@ -193,7 +193,7 @@ markers = [
     "examples: marks example tests (run only on releases)",
     "deprecated_api: marks tests using deprecated Optimization/Results API (remove in v6.0.0)",
 ]
-addopts = '-m "not examples"'  # Skip examples by default
+addopts = '-m "not examples" --ignore=tests/superseded'  # Skip examples and superseded tests by default
 
 # Warning filter configuration for pytest
 # Filters are processed in order; first match wins

tests/conftest.py

@@ -400,11 +400,8 @@ def gas_with_costs():
 # ============================================================================
 
 
-@pytest.fixture
-def simple_flow_system() -> fx.FlowSystem:
-    """
-    Create a simple energy system for testing
-    """
+def build_simple_flow_system() -> fx.FlowSystem:
+    """Create a simple energy system for testing (factory function)."""
     base_timesteps = pd.date_range('2020-01-01', periods=9, freq='h', name='time')
     timesteps_length = len(base_timesteps)
     base_thermal_load = LoadProfiles.thermal_simple(timesteps_length)
@@ -431,6 +428,12 @@ def simple_flow_system() -> fx.FlowSystem:
     return flow_system
 
 
+@pytest.fixture
+def simple_flow_system() -> fx.FlowSystem:
+    """Create a simple energy system for testing."""
+    return build_simple_flow_system()
+
+
 @pytest.fixture
 def simple_flow_system_scenarios() -> fx.FlowSystem:
     """

tests/test_flow_system_locking.py → tests/flow_system/test_flow_system_locking.py

@@ -12,7 +12,7 @@
 
 import flixopt as fx
 
-# Note: We use simple_flow_system fixture from conftest.py
+from ..conftest import build_simple_flow_system
 
 
 class TestIsLocked:
@@ -179,6 +179,14 @@ def test_reset_allows_reoptimization(self, simple_flow_system, highs_solver):
 class TestCopy:
     """Test the copy method."""
 
+    @pytest.fixture(scope='class')
+    def optimized_flow_system(self):
+        """Pre-optimized flow system shared across TestCopy (tests only work with copies)."""
+        fs = build_simple_flow_system()
+        solver = fx.solvers.HighsSolver(mip_gap=0, time_limit_seconds=300)
+        fs.optimize(solver)
+        return fs
+
     def test_copy_creates_new_instance(self, simple_flow_system):
         """Copy should create a new FlowSystem instance."""
         copy_fs = simple_flow_system.copy()
@@ -191,67 +199,58 @@ def test_copy_preserves_elements(self, simple_flow_system):
         assert set(copy_fs.components.keys()) == set(simple_flow_system.components.keys())
         assert set(copy_fs.buses.keys()) == set(simple_flow_system.buses.keys())
 
-    def test_copy_does_not_copy_solution(self, simple_flow_system, highs_solver):
+    def test_copy_does_not_copy_solution(self, optimized_flow_system):
         """Copy should not include the solution."""
-        simple_flow_system.optimize(highs_solver)
-        assert simple_flow_system.solution is not None
+        assert optimized_flow_system.solution is not None
 
-        copy_fs = simple_flow_system.copy()
+        copy_fs = optimized_flow_system.copy()
         assert copy_fs.solution is None
 
-    def test_copy_does_not_copy_model(self, simple_flow_system, highs_solver):
+    def test_copy_does_not_copy_model(self, optimized_flow_system):
         """Copy should not include the model."""
-        simple_flow_system.optimize(highs_solver)
-        assert simple_flow_system.model is not None
+        assert optimized_flow_system.model is not None
 
-        copy_fs = simple_flow_system.copy()
+        copy_fs = optimized_flow_system.copy()
         assert copy_fs.model is None
 
-    def test_copy_is_not_locked(self, simple_flow_system, highs_solver):
+    def test_copy_is_not_locked(self, optimized_flow_system):
         """Copy should not be locked even if original is."""
-        simple_flow_system.optimize(highs_solver)
-        assert simple_flow_system.is_locked is True
+        assert optimized_flow_system.is_locked is True
 
-        copy_fs = simple_flow_system.copy()
+        copy_fs = optimized_flow_system.copy()
         assert copy_fs.is_locked is False
 
-    def test_copy_can_be_modified(self, simple_flow_system, highs_solver):
+    def test_copy_can_be_modified(self, optimized_flow_system):
         """Copy should be modifiable even if original is locked."""
-        simple_flow_system.optimize(highs_solver)
-
-        copy_fs = simple_flow_system.copy()
+        copy_fs = optimized_flow_system.copy()
         new_bus = fx.Bus('NewBus')
         copy_fs.add_elements(new_bus)  # Should not raise
         assert 'NewBus' in copy_fs.buses
 
-    def test_copy_can_be_optimized_independently(self, simple_flow_system, highs_solver):
+    def test_copy_can_be_optimized_independently(self, optimized_flow_system):
         """Copy can be optimized independently of original."""
-        simple_flow_system.optimize(highs_solver)
-        original_cost = simple_flow_system.solution['costs'].item()
+        original_cost = optimized_flow_system.solution['costs'].item()
 
-        copy_fs = simple_flow_system.copy()
-        copy_fs.optimize(highs_solver)
+        copy_fs = optimized_flow_system.copy()
+        solver = fx.solvers.HighsSolver(mip_gap=0, time_limit_seconds=300)
+        copy_fs.optimize(solver)
 
         # Both should have solutions
-        assert simple_flow_system.solution is not None
+        assert optimized_flow_system.solution is not None
         assert copy_fs.solution is not None
 
         # Costs should be equal (same system)
         assert copy_fs.solution['costs'].item() == pytest.approx(original_cost)
 
-    def test_python_copy_uses_copy_method(self, simple_flow_system, highs_solver):
+    def test_python_copy_uses_copy_method(self, optimized_flow_system):
         """copy.copy() should use the custom copy method."""
-        simple_flow_system.optimize(highs_solver)
-
-        copy_fs = copy.copy(simple_flow_system)
+        copy_fs = copy.copy(optimized_flow_system)
         assert copy_fs.solution is None
         assert copy_fs.is_locked is False
 
-    def test_python_deepcopy_uses_copy_method(self, simple_flow_system, highs_solver):
+    def test_python_deepcopy_uses_copy_method(self, optimized_flow_system):
         """copy.deepcopy() should use the custom copy method."""
-        simple_flow_system.optimize(highs_solver)
-
-        copy_fs = copy.deepcopy(simple_flow_system)
+        copy_fs = copy.deepcopy(optimized_flow_system)
         assert copy_fs.solution is None
         assert copy_fs.is_locked is False
 

tests/test_io.py → tests/io/test_io.py

@@ -8,7 +8,7 @@
 
 import flixopt as fx
 
-from .conftest import (
+from ..conftest import (
     flow_system_base,
     flow_system_long,
     flow_system_segments_of_flows_2,

tests/test_io_conversion.py → tests/io/test_io_conversion.py

@@ -636,7 +636,7 @@ def test_load_old_results_from_resources(self):
 
         import flixopt as fx
 
-        resources_path = pathlib.Path(__file__).parent / 'ressources'
+        resources_path = pathlib.Path(__file__).parent.parent / 'ressources'
 
         # Load old results using new method
         fs = fx.FlowSystem.from_old_results(resources_path, 'Sim1')
@@ -655,7 +655,7 @@ def test_old_results_can_be_saved_new_format(self, tmp_path):
 
         import flixopt as fx
 
-        resources_path = pathlib.Path(__file__).parent / 'ressources'
+        resources_path = pathlib.Path(__file__).parent.parent / 'ressources'
 
         # Load old results
         fs = fx.FlowSystem.from_old_results(resources_path, 'Sim1')
@@ -674,7 +674,7 @@ def test_old_results_can_be_saved_new_format(self, tmp_path):
 class TestV4APIConversion:
     """Tests for converting v4 API result files to the new format."""
 
-    V4_API_PATH = pathlib.Path(__file__).parent / 'ressources' / 'v4-api'
+    V4_API_PATH = pathlib.Path(__file__).parent.parent / 'ressources' / 'v4-api'
 
     # All result names in the v4-api folder
     V4_RESULT_NAMES = [

tests/test_network_app.py → tests/plotting/test_network_app.py

@@ -2,7 +2,7 @@
 
 import flixopt as fx
 
-from .conftest import (
+from ..conftest import (
     flow_system_long,
     flow_system_segments_of_flows_2,
     simple_flow_system,

tests/superseded/__init__.py

@@ -0,0 +1,8 @@
+"""Superseded tests — replaced by tests/test_math/.
+
+These tests have been replaced by more thorough, analytically verified tests
+in tests/test_math/. They are kept temporarily for reference and will be
+deleted once confidence in the new test suite is established.
+
+All tests in this folder are skipped via pytestmark.
+"""

tests/superseded/math/__init__.py

@@ -0,0 +1,6 @@
+"""Model-building tests superseded by tests/test_math/.
+
+These tests verified linopy model structure (variables, constraints, bounds).
+They are implicitly covered by test_math: if solutions are mathematically correct,
+the model building must be correct.
+"""

tests/test_bus.py → tests/superseded/math/test_bus.py

@@ -1,6 +1,10 @@
+import pytest
+
 import flixopt as fx
 
-from .conftest import assert_conequal, assert_var_equal, create_linopy_model
+from ...conftest import assert_conequal, assert_var_equal, create_linopy_model
+
+pytestmark = pytest.mark.skip(reason='Superseded: model-building tests implicitly covered by tests/test_math/')
 
 
 class TestBusModel:

tests/test_component.py → tests/superseded/math/test_component.py

@@ -4,7 +4,7 @@
 import flixopt as fx
 import flixopt.elements
 
-from .conftest import (
+from ...conftest import (
     assert_almost_equal_numeric,
     assert_conequal,
     assert_dims_compatible,
@@ -13,6 +13,8 @@
     create_linopy_model,
 )
 
+pytestmark = pytest.mark.skip(reason='Superseded: model-building tests implicitly covered by tests/test_math/')
+
 
 class TestComponentModel:
     def test_flow_label_check(self):

tests/test_effect.py → tests/superseded/math/test_effect.py

@@ -4,13 +4,15 @@
 
 import flixopt as fx
 
-from .conftest import (
+from ...conftest import (
     assert_conequal,
     assert_sets_equal,
     assert_var_equal,
     create_linopy_model,
 )
 
+pytestmark = pytest.mark.skip(reason='Superseded: model-building tests implicitly covered by tests/test_math/')
+
 
 class TestEffectModel:
     """Test the FlowModel class."""