fix: compatibility with linopy v1 arithmetic semantics

docs/user-guide/linopy-v1-compatibility.md

@@ -0,0 +1,186 @@
+# linopy v1 Semantics Compatibility
+
+!!! info "Who this is for"
+    This is a **contributor / maintainer** guide. flixopt end users never
+    interact with linopy directly, so nothing in your models or scripts changes.
+    This page documents how flixopt's constraint-construction code stays correct
+    under linopy's upcoming *v1 arithmetic convention* and how to keep new code
+    compatible.
+
+---
+
+## Background
+
+[linopy PR #717](https://github.com/PyPSA/linopy/pull/717) introduces a global
+switch:
+
+```python
+import linopy
+linopy.options["semantics"] = "legacy"  # current default
+linopy.options["semantics"] = "v1"      # opt-in; becomes default before linopy 1.0
+```
+
+Legacy linopy silently *guessed* in three situations where a guess can quietly
+change a model. **v1** replaces each silent guess with an explicit rule and
+**raises** where legacy would have guessed, so a wrong model surfaces as a build
+error instead of a wrong number at solve time:
+
+| Legacy behaviour | Under v1 |
+|---|---|
+| Coordinate mismatch on a shared dim → aligned by position / inner join | **Raises** — align by label explicitly |
+| A conflicting auxiliary (non-dimension) coordinate → silently dropped | **Raises** — drop or relabel explicitly |
+| `NaN` in a constant / a shifted-in variable → filled with `0` (or `1` for `/`) | **Raises** / absence propagates |
+
+flixopt builds and solves **identically** under both settings. The full test
+suite is verified under both `"legacy"` and `"v1"`.
+
+---
+
+## The three patterns flixopt uses
+
+### 1. Adjacent-step constraints → `shift(-1)` + a legacy-only slice
+
+Constraints that relate `x[t+1]` to `x[t]` (storage balance, on/off duration
+tracking, state transitions, level tracking, linked-period sizing) used to slice
+the same variable two ways:
+
+```python
+x.isel(time=slice(1, None))    # x[t+1] — labels t[1:]
+x.isel(time=slice(None, -1))   # x[t]   — labels t[:-1]
+```
+
+The two slices carry **different labels** on `time`, so combining them is a
+coordinate mismatch. Legacy aligned them by position; v1 raises.
+
+**Fix:** express "next step" with `.shift({dim: -1})` — the v1-native form — and
+pair it with a legacy-only boundary slice:
+
+```python
+# Before (relies on legacy positional alignment):
+x.isel(time=slice(1, None)) - x.isel(time=slice(None, -1))
+
+# After (identical under both semantics):
+x.shift(time=-1).isel(time=slice(None, -1)) - x.isel(time=slice(None, -1))
+#                ^^^^^^^^^^^^^^^^^^^^^^^^^^^ legacy-only (see below)
+```
+
+`x.shift({dim: -1})` places `x[t+1]` at position `t`. Under **v1** the
+out-of-range boundary slot is marked *absent*, so that constraint row drops out
+on its own — the shift alone is correct. Under **legacy** `.shift` instead
+*fills* the boundary slot with `0`, leaving a spurious row, so the trailing
+`.isel({dim: slice(None, -1)})` drops it. With the slice, both conventions
+produce **identical** constraints.
+
+!!! note "The `.isel(...)` is a legacy-only shim"
+    Pure `x - x.shift(...)` on its own is **not** identical across conventions —
+    that is exactly why the slice is there. But it is needed *only* for legacy:
+    once legacy support is dropped, every `.isel({dim: slice(None, -1)})` on a
+    shifted operand can be removed, leaving bare `.shift(...)` constraints (v1
+    masks the boundary natively). flixopt marks each such slice as legacy-only
+    (see the module note in `modeling.py`) so the future cleanup is mechanical.
+
+### 2. Scalar-endpoint constraints → `drop=True`
+
+A scalar integer index keeps the indexed coordinate as a **leftover scalar
+aux-coord**:
+
+```python
+x.isel(time=0)    # drops the `time` *dimension*, keeps `time` = first timestamp
+x.isel(time=-1)   # keeps `time` = last timestamp
+```
+
+Cyclic / initial / final constraints compare two such endpoints:
+
+```python
+x.isel(time=0) == x.isel(time=-1)   # `time` = first vs last → conflict
+```
+
+Legacy silently drops the conflicting coord; v1 raises. **Fix:** pass
+`drop=True` so the meaningless scalar coord never appears:
+
+```python
+x.isel(time=0, drop=True) == x.isel(time=-1, drop=True)
+```
+
+This applies to any scalar `.isel(...)` whose result is combined with another
+operand — `time`, `scenario`, `cluster_boundary`, and to a value *derived* from
+such a select (e.g. `previous_status.isel(time=-1, drop=True)`).
+
+### 3. Constraint mutation → build once, or `_set_constraint_lhs`
+
+The `Constraint.lhs` / `.rhs` **setters** are deprecated in linopy 0.8 (they were
+the only mutation API in 0.7). Two cases:
+
+**Prefer building the constraint in one pass** rather than creating it and then
+mutating it. Where a term is only conditionally present, assemble the full LHS
+first (this is what the transmission-loss and bus-balance constraints do now):
+
+```python
+# Before:
+con = model.add_constraints(a == b)
+if extra:
+    con.lhs += term          # deprecated setter
+
+# After:
+lhs = a
+if extra:
+    lhs = lhs + term
+con = model.add_constraints(lhs == b)
+```
+
+**For a genuine incremental accumulator** (e.g. the effect-share totals, which
+are grown across the whole model build) mutation is unavoidable. Use the
+`_set_constraint_lhs` shim, which works on **both** linopy versions — 0.8+ has
+`Constraint.update(lhs=...)` (and deprecates the setter), while 0.7 has only the
+setter (not deprecated there):
+
+```python
+from flixopt.modeling import _set_constraint_lhs
+
+_set_constraint_lhs(con, con.lhs - share)   # update() on >=0.8, setter on <0.8
+```
+
+---
+
+## Writing v1-safe constraints (checklist)
+
+- Relating `x[t+1]` to `x[t]`? → write the next step as `x.shift({dim: -1})` and
+  add the legacy-only `.isel({dim: slice(None, -1)})`; never combine raw
+  `isel(slice(1, None))` with `isel(slice(None, -1))`.
+- Scalar `.isel(dim=k)` whose result feeds arithmetic or a comparison? → add
+  `drop=True`.
+- Mutating an existing constraint? → `con.update(lhs=...)`, not `con.lhs = ...`.
+- Passing a plain numpy array / list into linopy arithmetic? → wrap it in a
+  named `DataArray` so it aligns by label, not by size.
+- Never introduce a user-supplied `NaN` as a stand-in for "absent" — express
+  absence with `mask=` / `.where(...)` on the variable.
+
+---
+
+## Testing under v1
+
+Run the suite with the option flipped (e.g. via a `conftest`/plugin that sets
+`linopy.options["semantics"] = "v1"` at configure time), and turn the legacy
+warning into an error to surface every remaining site under legacy:
+
+```python
+import warnings
+from linopy import LinopySemanticsWarning
+warnings.filterwarnings("error", category=LinopySemanticsWarning)
+```
+
+A green run under **both** `"legacy"` and `"v1"` is the compatibility contract.
+
+---
+
+## Known non-issue: cyclic-cluster SOC degeneracy
+
+A handful of clustered-storage tests assert an exact absolute charge-state
+trajectory. In `cluster_mode='cyclic'` the initial/final charge-state
+constraints are intentionally skipped, so the **absolute** SOC level of each
+representative cluster is unpinned — only the per-cluster *delta* is fixed. The
+problem is therefore degenerate in the SOC level: multiple equally-optimal
+solutions exist, and different solver/linopy builds may return a different (but
+equally valid) vertex. Objective and all flow rates are unaffected. This is a
+test-robustness matter, **independent of the semantics setting** (it reproduces
+identically under `"legacy"` and `"v1"`).

flixopt/components.py

@@ -843,15 +843,14 @@ def create_transmission_equation(self, name: str, in_flow: Flow, out_flow: Flow)
         """Creates an Equation for the Transmission efficiency and adds it to the model"""
         # eq: out(t) + on(t)*loss_abs(t) = in(t)*(1 - loss_rel(t))
         rel_losses = 0 if self.element.relative_losses is None else self.element.relative_losses
-        con_transmission = self.add_constraints(
-            out_flow.submodel.flow_rate == in_flow.submodel.flow_rate * (1 - rel_losses),
-            short_name=name,
-        )
-
+        lhs = out_flow.submodel.flow_rate
         if (self.element.absolute_losses is not None) and np.any(self.element.absolute_losses != 0):
-            con_transmission.lhs += in_flow.submodel.status.status * self.element.absolute_losses
+            lhs = lhs + in_flow.submodel.status.status * self.element.absolute_losses
 
-        return con_transmission
+        return self.add_constraints(
+            lhs == in_flow.submodel.flow_rate * (1 - rel_losses),
+            short_name=name,
+        )
 
 
 class LinearConverterModel(ComponentModel):
@@ -974,7 +973,7 @@ def _add_cluster_cyclic_constraint(self):
         """For 'cyclic' cluster mode: each cluster's start equals its end."""
         if self._model.flow_system.clusters is not None and self.element.cluster_mode == 'cyclic':
             self.add_constraints(
-                self.charge_state.isel(time=0) == self.charge_state.isel(time=-2),
+                self.charge_state.isel(time=0, drop=True) == self.charge_state.isel(time=-2, drop=True),
                 short_name='cluster_cyclic',
             )
 
@@ -1018,7 +1017,7 @@ def _add_initial_final_constraints(self):
         if self.element.initial_charge_state is not None:
             if isinstance(self.element.initial_charge_state, str):
                 self.add_constraints(
-                    self.charge_state.isel(time=0) == self.charge_state.isel(time=-1),
+                    self.charge_state.isel(time=0, drop=True) == self.charge_state.isel(time=-1, drop=True),
                     short_name='initial_charge_state',
                 )
             else:
@@ -1072,8 +1071,10 @@ def _build_energy_balance_lhs(self):
         eff_charge = self.element.eta_charge
         eff_discharge = self.element.eta_discharge
 
+        # charge_state[t+1] via shift(-1); the `.isel(time=slice(None, -1))` on the
+        # shifted term is a legacy-only boundary drop (v1 masks it natively).
         return (
-            charge_state.isel(time=slice(1, None))
+            charge_state.shift(time=-1).isel(time=slice(None, -1))
             - charge_state.isel(time=slice(None, -1)) * ((1 - rel_loss) ** timestep_duration)
             - charge_rate * eff_charge * timestep_duration
             + discharge_rate * timestep_duration / eff_discharge
@@ -1415,7 +1416,8 @@ def _add_intercluster_linking(self) -> None:
         # 6. Add cyclic or initial constraint
         if self.element.cluster_mode == 'intercluster_cyclic':
             self.add_constraints(
-                soc_boundary.isel(cluster_boundary=0) == soc_boundary.isel(cluster_boundary=n_original_clusters),
+                soc_boundary.isel(cluster_boundary=0, drop=True)
+                == soc_boundary.isel(cluster_boundary=n_original_clusters, drop=True),
                 short_name='cyclic',
             )
         else:
@@ -1425,13 +1427,13 @@ def _add_intercluster_linking(self) -> None:
                 if isinstance(initial, str):
                     # 'equals_final' means cyclic
                     self.add_constraints(
-                        soc_boundary.isel(cluster_boundary=0)
-                        == soc_boundary.isel(cluster_boundary=n_original_clusters),
+                        soc_boundary.isel(cluster_boundary=0, drop=True)
+                        == soc_boundary.isel(cluster_boundary=n_original_clusters, drop=True),
                         short_name='initial_SOC_boundary',
                     )
                 else:
                     self.add_constraints(
-                        soc_boundary.isel(cluster_boundary=0) == initial,
+                        soc_boundary.isel(cluster_boundary=0, drop=True) == initial,
                         short_name='initial_SOC_boundary',
                     )
 
@@ -1482,7 +1484,7 @@ def _compute_delta_soc(self, n_clusters: int, timesteps_per_cluster: int) -> xr.
             DataArray with 'cluster' dimension containing delta_SOC for each cluster.
         """
         # With 2D structure: result already has cluster dimension
-        return self.charge_state.isel(time=-1) - self.charge_state.isel(time=0)
+        return self.charge_state.isel(time=-1, drop=True) - self.charge_state.isel(time=0, drop=True)
 
     def _add_linking_constraints(
         self,
@@ -1595,7 +1597,7 @@ def _add_combined_bound_constraints(
 
         for sample_name, offset in zip(['start', 'mid', 'end'], sample_offsets, strict=False):
             # With 2D structure: select time offset, then reorder by cluster_assignments
-            cs_at_offset = charge_state.isel(time=offset)  # Shape: (cluster, ...)
+            cs_at_offset = charge_state.isel(time=offset, drop=True)  # Shape: (cluster, ...)
             # Reorder to original_cluster order using cluster_assignments indexer
             cs_t = cs_at_offset.isel(cluster=cluster_assignments)
             # Suppress xarray warning about index loss - we immediately assign new coords anyway
@@ -1607,7 +1609,7 @@ def _add_combined_bound_constraints(
             # Apply decay factor (1-loss)^hours to SOC_boundary per Eq. 9
             # Convert timestep offset to hours using cumulative duration for non-uniform timesteps
             if cumulative_hours is not None:
-                hours_offset = cumulative_hours.isel(time=offset)
+                hours_offset = cumulative_hours.isel(time=offset, drop=True)
             else:
                 hours_offset = offset * mean_timestep_duration
             decay_t = (1 - rel_loss) ** hours_offset

flixopt/elements.py

@@ -1012,9 +1012,9 @@ def _do_modeling(self):
             self.register_variable(flow.submodel.flow_rate, flow.label_full)
         inputs = sum([flow.submodel.flow_rate for flow in self.element.inputs.values()])
         outputs = sum([flow.submodel.flow_rate for flow in self.element.outputs.values()])
-        eq_bus_balance = self.add_constraints(inputs == outputs, short_name='balance')
+        balance = inputs - outputs
 
-        # Add virtual supply/demand to balance and penalty if needed
+        # Add virtual supply/demand to the balance and penalty if needed
         if self.element.allows_imbalance:
             imbalance_penalty = self.element.imbalance_penalty_per_flow_hour * self._model.timestep_duration
 
@@ -1033,7 +1033,7 @@ def _do_modeling(self):
             )
 
             # Σ(inflows) + virtual_supply = Σ(outflows) + virtual_demand
-            eq_bus_balance.lhs += self.virtual_supply - self.virtual_demand
+            balance = balance + self.virtual_supply - self.virtual_demand
 
             # Add penalty shares as temporal effects (time-dependent)
             from .effects import PENALTY_EFFECT_LABEL
@@ -1045,6 +1045,8 @@ def _do_modeling(self):
                 target='temporal',
             )
 
+        self.add_constraints(balance == 0, short_name='balance')
+
     def results_structure(self):
         inputs = [flow.submodel.flow_rate.name for flow in self.element.inputs.values()]
         outputs = [flow.submodel.flow_rate.name for flow in self.element.outputs.values()]

flixopt/features.py

@@ -10,7 +10,7 @@
 import linopy
 import numpy as np
 
-from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilities
+from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilities, _set_constraint_lhs
 from .structure import FlowSystemModel, Submodel, VariableCategory
 
 if TYPE_CHECKING:
@@ -91,8 +91,10 @@ def _create_variables_and_constraints(self):
 
         if self.parameters.linked_periods is not None:
             masked_size = self.size.where(self.parameters.linked_periods, drop=True)
+            # size[p] == size[p+1]; size[p+1] via shift(-1). The `.isel(period=slice(None, -1))`
+            # on the shifted term is a legacy-only boundary drop (v1 masks it natively).
             self.add_constraints(
-                masked_size.isel(period=slice(None, -1)) == masked_size.isel(period=slice(1, None)),
+                masked_size.isel(period=slice(None, -1)) == masked_size.shift(period=-1).isel(period=slice(None, -1)),
                 short_name='linked_periods',
             )
 
@@ -236,7 +238,9 @@ def _do_modeling(self):
             )
 
             # Determine previous_state: None means relaxed (no constraint at t=0)
-            previous_state = self._previous_status.isel(time=-1) if self._previous_status is not None else None
+            previous_state = (
+                self._previous_status.isel(time=-1, drop=True) if self._previous_status is not None else None
+            )
 
             BoundingPatterns.state_transition_bounds(
                 self,
@@ -295,7 +299,7 @@ def _add_cluster_cyclic_constraint(self):
         """For 'cyclic' cluster mode: each cluster's start status equals its end status."""
         if self._model.flow_system.clusters is not None and self.parameters.cluster_mode == 'cyclic':
             self.add_constraints(
-                self.status.isel(time=0) == self.status.isel(time=-1),
+                self.status.isel(time=0, drop=True) == self.status.isel(time=-1, drop=True),
                 short_name='cluster_cyclic',
             )
 
@@ -664,7 +668,9 @@ def _do_modeling(self):
             # Add it to the total (cluster_weight handles cluster representation, defaults to 1.0)
             # Sum over all temporal dimensions (time, and cluster if present)
             weighted_per_timestep = self.total_per_timestep * self._model.weights.get('cluster', 1.0)
-            self._eq_total.lhs -= weighted_per_timestep.sum(dim=self._model.temporal_dims)
+            _set_constraint_lhs(
+                self._eq_total, self._eq_total.lhs - weighted_per_timestep.sum(dim=self._model.temporal_dims)
+            )
 
     def add_share(
         self,
@@ -694,7 +700,7 @@ def add_share(
                 raise ValueError('Cannot add share with scenario-dim to a model without scenario-dim')
 
         if name in self.shares:
-            self.share_constraints[name].lhs -= expression
+            _set_constraint_lhs(self.share_constraints[name], self.share_constraints[name].lhs - expression)
         else:
             # Temporal shares (with 'time' dim) are segment totals that need division
             category = VariableCategory.SHARE if 'time' in dims else None
@@ -710,6 +716,6 @@ def add_share(
             )
 
             if 'time' not in dims:
-                self._eq_total.lhs -= self.shares[name]
+                _set_constraint_lhs(self._eq_total, self._eq_total.lhs - self.shares[name])
             else:
-                self._eq_total_per_timestep.lhs -= self.shares[name]
+                _set_constraint_lhs(self._eq_total_per_timestep, self._eq_total_per_timestep.lhs - self.shares[name])