flixOpt
diff --git a/‎.github/workflows/docs.yaml‎
Lines changed: 9 additions & 0 deletions b/‎.github/workflows/docs.yaml‎
Lines changed: 9 additions & 0 deletions
diff --git a/‎CHANGELOG.md‎
Lines changed: 126 additions & 6 deletions b/‎CHANGELOG.md‎
Lines changed: 126 additions & 6 deletions
diff --git a/‎CITATION.cff‎
Lines changed: 2 additions & 2 deletions b/‎CITATION.cff‎
Lines changed: 2 additions & 2 deletions
diff --git a/‎README.md‎
Lines changed: 1 addition & 1 deletion b/‎README.md‎
Lines changed: 1 addition & 1 deletion
diff --git a/‎docs/notebooks/08c-clustering.ipynb‎
Lines changed: 1 addition & 1 deletion b/‎docs/notebooks/08c-clustering.ipynb‎
Lines changed: 1 addition & 1 deletion
diff --git a/‎docs/notebooks/08d-clustering-multiperiod.ipynb‎
Lines changed: 1 addition & 0 deletions b/‎docs/notebooks/08d-clustering-multiperiod.ipynb‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎flixopt/comparison.py‎
Lines changed: 77 additions & 14 deletions b/‎flixopt/comparison.py‎
Lines changed: 77 additions & 14 deletions
diff --git a/‎flixopt/components.py‎
Lines changed: 14 additions & 4 deletions b/‎flixopt/components.py‎
Lines changed: 14 additions & 4 deletions
@@ -12,6 +12,15 @@ on:
       - 'docs/**'
       - 'mkdocs.yml'
   workflow_dispatch:
+    inputs:
+      deploy:
+        description: 'Deploy docs to GitHub Pages'
+        type: boolean
+        default: false
+      version:
+        description: 'Version to deploy (e.g., v6.0.0)'
+        type: string
+        required: false
   workflow_call:
     inputs:
       deploy:
 
@@ -52,7 +52,70 @@ If upgrading from v2.x, see the [v3.0.0 release notes](https://github.com/flixOp
 
 Until here -->
 
-## [6.0.0] - Upcoming
+## [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.
+
+### 🐛 Fixed
+
+- **Comparison Coordinates**: Fixed `component` coordinate becoming `(case, contributor)` shaped after concatenation in `Comparison` class. Non-index coordinates are now properly merged before concat in `solution`, `inputs`, and all statistics properties. Added warning when coordinate mappings conflict (#599)
+
+### 📝 Docs
+
+- **Docs Workflow**: Added `workflow_dispatch` inputs for manual docs deployment with version selection (#599)
+
+### 👷 Development
+
+- Updated dev dependencies to newer versions
+---
+
+## [6.0.1] - 2026-02-04
+
+**Summary**: Bugfix release addressing clustering issues with multi-period systems and ExtremeConfig.
+
+### 🐛 Fixed
+
+- **Multi-period clustering with ExtremeConfig** - Fixed `ValueError: cannot reshape array` when clustering multi-period or multi-scenario systems with `ExtremeConfig`. The fix uses pandas `.unstack()` instead of manual reshape for robustness.
+- **Consistent cluster count validation** - Added validation to detect inconsistent cluster counts across periods/scenarios, providing clear error messages.
+
+### 💥 Breaking Changes
+
+- **ExtremeConfig method restriction for multi-period systems** - When using `ExtremeConfig` with multi-period or multi-scenario systems, only `method='replace'` is now allowed. Using `method='new_cluster'` or `method='append'` will raise a `ValueError`. This works around a tsam bug where these methods can produce inconsistent cluster counts across slices.
+
+### 📦 Dependencies
+
+- Excluded tsam 3.1.0 from compatible versions due to clustering bug.
+
+---
+
+## [6.0.0] - 2026-02-03
 
 **Summary**: Major release featuring tsam v3 migration, complete rewrite of the clustering/aggregation system, 2-3x faster I/O for large systems, new `plotly` plotting accessor, FlowSystem comparison tools, and removal of deprecated v5.0 classes.
 
@@ -226,12 +289,12 @@ comp = fx.Comparison([fs_base, fs_modified])
 comp = fx.Comparison([fs1, fs2, fs3], names=['baseline', 'low_cost', 'high_eff'])
 
 # Side-by-side plots (auto-facets by 'case' dimension)
-comp.statistics.plot.balance('Heat')
-comp.statistics.flow_rates.plotly.line()
+comp.stats.plot.balance('Heat')
+comp.stats.flow_rates.plotly.line()
 
 # Access combined data with 'case' dimension
 comp.solution  # xr.Dataset
-comp.statistics.flow_rates  # xr.Dataset
+comp.stats.flow_rates  # xr.Dataset
 
 # Compute differences relative to a reference case
 comp.diff()  # vs first case
@@ -262,6 +325,58 @@ flow_system.topology.set_component_colors('turbo', overwrite=False)  # Only unse
 
 `Component.inputs`, `Component.outputs`, and `Component.flows` now use `FlowContainer` (dict-like) with dual access by index or label: `inputs[0]` or `inputs['Q_th']`.
 
+#### `before_solve` Callback
+
+New callback parameter for `optimize()` and `rolling_horizon()` allows adding custom constraints before solving:
+
+```python
+def add_constraints(fs):
+    model = fs.model
+    boiler = model.variables['Boiler(Q_th)|flow_rate']
+    model.add_constraints(boiler >= 10, name='min_boiler')
+
+flow_system.optimize(solver, before_solve=add_constraints)
+
+# Works with rolling_horizon too
+flow_system.optimize.rolling_horizon(
+    solver,
+    horizon=168,
+    before_solve=add_constraints
+)
+```
+
+#### `cluster_mode` for StatusParameters
+
+New parameter to control status behavior at cluster boundaries:
+
+```python
+fx.StatusParameters(
+    ...,
+    cluster_mode='relaxed',  # Default: no constraint at boundaries, prevents phantom startups
+    # cluster_mode='cyclic',  # Each cluster's final status equals its initial status
+)
+```
+
+#### Comparison Class Enhancements
+
+- **`Comparison.inputs`**: Compare inputs across FlowSystems for easy side-by-side input parameter comparison
+- **`data_only` parameter**: Get data without generating plots in Comparison methods
+- **`threshold` parameter**: Filter small values when comparing
+
+#### Plotting Enhancements
+
+- **`threshold` parameter**: Added to all plotting methods to filter values below a threshold (default: `1e-5`)
+- **`round_decimals` parameter**: Control decimal precision in `balance()`, `carrier_balance()`, and `storage()` plots
+- **`flow_colors` property**: Map flows to their component's colors for consistent visualization
+
+#### `FlowSystem.from_old_dataset()`
+
+New method for loading datasets saved with older flixopt versions:
+
+```python
+fs = fx.FlowSystem.from_old_dataset(old_dataset)
+```
+
 ### 💥 Breaking Changes
 
 #### tsam v3 Migration
@@ -306,17 +421,22 @@ fs.transform.cluster(
 
 - `FlowSystem.weights` returns `dict[str, xr.DataArray]` (unit weights instead of `1.0` float fallback)
 - `FlowSystemDimensions` type now includes `'cluster'`
-- `statistics.plot.balance()`, `carrier_balance()`, and `storage()` now use `xarray_plotly.fast_bar()` internally (styled stacked areas for better performance)
+- `stats.plot.balance()`, `carrier_balance()`, and `storage()` now use `xarray_plotly.fast_bar()` internally (styled stacked areas for better performance)
+- `stats.plot.carrier_balance()` now combines inputs and outputs to show net flow per component, and aggregates per component by default
 
 ### 🗑️ Deprecated
 
 The following items are deprecated and will be removed in **v7.0.0**:
 
+**Accessor renamed:**
+
+- `flow_system.statistics` → Use `flow_system.stats` (shorter, more convenient)
+
 **Classes** (use FlowSystem methods instead):
 
 - `Optimization` class → Use `flow_system.optimize(solver)`
 - `SegmentedOptimization` class → Use `flow_system.optimize.rolling_horizon()`
-- `Results` class → Use `flow_system.solution` and `flow_system.statistics`
+- `Results` class → Use `flow_system.solution` and `flow_system.stats`
 - `SegmentedResults` class → Use segment FlowSystems directly
 
 **FlowSystem methods** (use `transform` or `topology` accessor instead):
 
@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: "If you use this software, please cite it as below and consider citing the related publication."
 type: software
 title: "flixopt"
-version: 6.0.0rc17
-date-released: 2026-02-02
+version: 6.0.2
+date-released: 2026-02-05
 url: "https://github.com/flixOpt/flixopt"
 repository-code: "https://github.com/flixOpt/flixopt"
 license: MIT
 
@@ -54,7 +54,7 @@ flow_system.optimize(fx.solvers.HighsSolver())
 
 # 3. Analyze results
 flow_system.solution        # Raw xarray Dataset
-flow_system.statistics      # Convenient analysis accessor
+flow_system.stats           # Convenient analysis accessor
 ```
 
 **Get started with real examples:**
 
@@ -585,7 +585,7 @@
    "id": "37",
    "metadata": {},
    "source": [
-    "## API Reference\n",
+    "<cell_type>markdown</cell_type>## API Reference\n",
     "\n",
     "### `transform.cluster()` Parameters\n",
     "\n",
 
@@ -556,6 +556,7 @@
     "fs = fs.transform.isel(time=slice(0, 168))  # First 168 timesteps\n",
     "\n",
     "# Cluster (applies per period/scenario)\n",
+    "# Note: For multi-period systems, only method='replace' is supported\n",
     "fs_clustered = fs.transform.cluster(\n",
     "    n_clusters=10,\n",
     "    cluster_duration='1D',\n",
 
@@ -29,6 +29,69 @@
 _CASE_SLOTS = frozenset(slot for slots in SLOT_ORDERS.values() for slot in slots)
 
 
+def _extract_nonindex_coords(datasets: list[xr.Dataset]) -> tuple[list[xr.Dataset], dict[str, tuple[str, dict]]]:
+    """Extract and merge non-index coords, returning cleaned datasets and merged mappings.
+
+    Non-index coords (like `component` on `contributor` dim) cause concat conflicts.
+    This extracts them, merges the mappings, and returns datasets without them.
+    """
+    if not datasets:
+        return datasets, {}
+
+    # Find non-index coords and collect mappings
+    merged: dict[str, tuple[str, dict]] = {}
+    coords_to_drop: set[str] = set()
+
+    for ds in datasets:
+        for name, coord in ds.coords.items():
+            if len(coord.dims) != 1:
+                continue
+            dim = coord.dims[0]
+            if dim == name or dim not in ds.coords:
+                continue
+
+            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,
+                )
+                continue
+
+            for dv, cv in zip(ds.coords[dim].values, coord.values, strict=False):
+                if dv not in merged[name][1]:
+                    merged[name][1][dv] = cv
+                elif merged[name][1][dv] != cv:
+                    warnings.warn(
+                        f"Coordinate '{name}' has conflicting values for dim value '{dv}': "
+                        f"'{merged[name][1][dv]}' vs '{cv}'. Keeping first value.",
+                        stacklevel=4,
+                    )
+
+    # Drop these coords from datasets
+    if coords_to_drop:
+        datasets = [ds.drop_vars(coords_to_drop, errors='ignore') for ds in datasets]
+
+    return datasets, merged
+
+
+def _apply_merged_coords(ds: xr.Dataset, merged: dict[str, tuple[str, dict]]) -> xr.Dataset:
+    """Apply merged coord mappings to concatenated dataset."""
+    if not merged:
+        return ds
+
+    new_coords = {}
+    for name, (dim, mapping) in merged.items():
+        if dim not in ds.dims:
+            continue
+        new_coords[name] = (dim, [mapping.get(dv, dv) for dv in ds.coords[dim].values])
+
+    return ds.assign_coords(new_coords)
+
+
 def _apply_slot_defaults(plotly_kwargs: dict, defaults: dict[str, str | None]) -> None:
     """Apply default slot assignments to plotly kwargs.
 
@@ -256,12 +319,10 @@ def solution(self) -> xr.Dataset:
             self._require_solutions()
             datasets = [fs.solution for fs in self._systems]
             self._warn_mismatched_dimensions(datasets)
-            self._solution = xr.concat(
-                [ds.expand_dims(case=[name]) for ds, name in zip(datasets, self._names, strict=True)],
-                dim='case',
-                join='outer',
-                fill_value=float('nan'),
-            )
+            expanded = [ds.expand_dims(case=[name]) for ds, name in zip(datasets, self._names, strict=True)]
+            expanded, merged_coords = _extract_nonindex_coords(expanded)
+            result = xr.concat(expanded, dim='case', join='outer', coords='minimal', fill_value=float('nan'))
+            self._solution = _apply_merged_coords(result, merged_coords)
         return self._solution
 
     @property
@@ -324,12 +385,10 @@ def inputs(self) -> xr.Dataset:
         if self._inputs is None:
             datasets = [fs.to_dataset(include_solution=False) for fs in self._systems]
             self._warn_mismatched_dimensions(datasets)
-            self._inputs = xr.concat(
-                [ds.expand_dims(case=[name]) for ds, name in zip(datasets, self._names, strict=True)],
-                dim='case',
-                join='outer',
-                fill_value=float('nan'),
-            )
+            expanded = [ds.expand_dims(case=[name]) for ds, name in zip(datasets, self._names, strict=True)]
+            expanded, merged_coords = _extract_nonindex_coords(expanded)
+            result = xr.concat(expanded, dim='case', join='outer', coords='minimal', fill_value=float('nan'))
+            self._inputs = _apply_merged_coords(result, merged_coords)
         return self._inputs
 
 
@@ -374,7 +433,9 @@ def _concat_property(self, prop_name: str) -> xr.Dataset:
                 continue
         if not datasets:
             return xr.Dataset()
-        return xr.concat(datasets, dim='case', join='outer', fill_value=float('nan'))
+        datasets, merged_coords = _extract_nonindex_coords(datasets)
+        result = xr.concat(datasets, dim='case', join='outer', coords='minimal', fill_value=float('nan'))
+        return _apply_merged_coords(result, merged_coords)
 
     def _merge_dict_property(self, prop_name: str) -> dict[str, str]:
         """Merge a dict property from all cases (later cases override)."""
@@ -528,7 +589,9 @@ def _combine_data(self, method_name: str, *args, **kwargs) -> tuple[xr.Dataset,
         if not datasets:
             return xr.Dataset(), ''
 
-        return xr.concat(datasets, dim='case', join='outer', fill_value=float('nan')), title
+        datasets, merged_coords = _extract_nonindex_coords(datasets)
+        combined = xr.concat(datasets, dim='case', join='outer', coords='minimal', fill_value=float('nan'))
+        return _apply_merged_coords(combined, merged_coords), title
 
     def _finalize(self, ds: xr.Dataset, fig, show: bool | None) -> PlotResult:
         """Handle show and return PlotResult."""
 
@@ -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)