diff --git a/docs/user-guide/linopy-v1-compatibility.md b/docs/user-guide/linopy-v1-compatibility.md
new file mode 100644
index 000000000..44723ba9f
--- /dev/null
+++ b/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"`).
diff --git a/flixopt/components.py b/flixopt/components.py
index 06313d7f6..f219c9d7f 100644
--- a/flixopt/components.py
+++ b/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
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 446ef4bd7..d5478a76f 100644
--- a/flixopt/elements.py
+++ b/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()]
diff --git a/flixopt/features.py b/flixopt/features.py
index e85636435..8b9bc9c59 100644
--- a/flixopt/features.py
+++ b/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])
diff --git a/flixopt/modeling.py b/flixopt/modeling.py
index ff84c808f..676a78993 100644
--- a/flixopt/modeling.py
+++ b/flixopt/modeling.py
@@ -10,6 +10,16 @@
 
 logger = logging.getLogger('flixopt')
 
+# Adjacent-step constraints in this module relate x[t+1] to x[t]. The "next step"
+# is written `x.shift({dim: -1})` (x[t+1] placed at position t) — the v1-native
+# form: under the linopy v1 convention the out-of-range boundary slot is marked
+# *absent*, so that row drops from the constraint on its own. The trailing
+# `.isel({dim: slice(None, -1)})` paired with each such shift is a LEGACY-ONLY
+# shim: under the legacy convention `.shift` *fills* the boundary slot instead of
+# masking it, leaving a spurious row that the slice removes. Once legacy support
+# is dropped, every `.isel({dim: slice(None, -1)})` on a shifted operand here can
+# go and the constraints reduce to their bare `.shift(...)` form.
+
 
 def _scalar_safe_isel(data: xr.DataArray | Any, indexers: dict) -> xr.DataArray | Any:
     """Apply isel if data has the required dimensions, otherwise return data as-is.
@@ -76,6 +86,22 @@ def _scalar_safe_reduce(data: xr.DataArray | Any, dim: str, method: str = 'mean'
     return data
 
 
+def _set_constraint_lhs(constraint: linopy.Constraint, lhs) -> None:
+    """Replace a constraint's LHS, compatibly across linopy versions.
+
+    Incremental accumulators (e.g. share totals) must mutate a constraint after
+    it is created. linopy >= 0.8 exposes ``Constraint.update(lhs=...)`` and
+    deprecates the ``.lhs`` setter (flixopt escalates that deprecation to an
+    error); linopy < 0.8 has only the setter, which is not deprecated there.
+    Dispatch on whichever API the installed linopy provides so the same call
+    works — and warns on neither — under both.
+    """
+    if hasattr(constraint, 'update'):
+        constraint.update(lhs=lhs)
+    else:
+        constraint.lhs = lhs
+
+
 def _xr_allclose(a: xr.DataArray, b: xr.DataArray, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
     """Check if two DataArrays are element-wise equal within tolerance.
 
@@ -405,17 +431,17 @@ def consecutive_duration_tracking(
 
         # Forward constraint: duration[t+1] ≤ duration[t] + duration_per_step[t]
         constraints['forward'] = model.add_constraints(
-            duration.isel({duration_dim: slice(1, None)})
+            duration.shift({duration_dim: -1}).isel({duration_dim: slice(None, -1)})
             <= duration.isel({duration_dim: slice(None, -1)}) + duration_per_step.isel({duration_dim: slice(None, -1)}),
             name=f'{duration.name}|forward',
         )
 
         # Backward constraint: duration[t+1] ≥ duration[t] + duration_per_step[t] + (state[t+1] - 1) * M
         constraints['backward'] = model.add_constraints(
-            duration.isel({duration_dim: slice(1, None)})
+            duration.shift({duration_dim: -1}).isel({duration_dim: slice(None, -1)})
             >= duration.isel({duration_dim: slice(None, -1)})
             + duration_per_step.isel({duration_dim: slice(None, -1)})
-            + (state.isel({duration_dim: slice(1, None)}) - 1) * mega,
+            + (state.shift({duration_dim: -1}).isel({duration_dim: slice(None, -1)}) - 1) * mega,
             name=f'{duration.name}|backward',
         )
 
@@ -431,8 +457,11 @@ def consecutive_duration_tracking(
         # Minimum duration constraint if provided
         if minimum_duration is not None:
             constraints['lb'] = model.add_constraints(
-                duration
-                >= (state.isel({duration_dim: slice(None, -1)}) - state.isel({duration_dim: slice(1, None)}))
+                duration.isel({duration_dim: slice(None, -1)})
+                >= (
+                    state.isel({duration_dim: slice(None, -1)})
+                    - state.shift({duration_dim: -1}).isel({duration_dim: slice(None, -1)})
+                )
                 * _scalar_safe_isel(minimum_duration, {duration_dim: slice(None, -1)}),
                 name=f'{duration.name}|lb',
             )
@@ -714,8 +743,9 @@ def state_transition_bounds(
 
         # State transition constraints for t > 0
         transition = model.add_constraints(
-            activate.isel({coord: slice(1, None)}) - deactivate.isel({coord: slice(1, None)})
-            == state.isel({coord: slice(1, None)}) - state.isel({coord: slice(None, -1)}),
+            activate.shift({coord: -1}).isel({coord: slice(None, -1)})
+            - deactivate.shift({coord: -1}).isel({coord: slice(None, -1)})
+            == state.shift({coord: -1}).isel({coord: slice(None, -1)}) - state.isel({coord: slice(None, -1)}),
             name=f'{name}|transition',
         )
 
@@ -776,14 +806,26 @@ def continuous_transition_bounds(
 
         # Transition constraints for t > 0: continuous variable can only change when transitions occur
         transition_upper = model.add_constraints(
-            continuous_variable.isel({coord: slice(1, None)}) - continuous_variable.isel({coord: slice(None, -1)})
-            <= max_change * (activate.isel({coord: slice(1, None)}) + deactivate.isel({coord: slice(1, None)})),
+            continuous_variable.shift({coord: -1}).isel({coord: slice(None, -1)})
+            - continuous_variable.isel({coord: slice(None, -1)})
+            <= max_change
+            * (
+                activate.shift({coord: -1}).isel({coord: slice(None, -1)})
+                + deactivate.shift({coord: -1}).isel({coord: slice(None, -1)})
+            ),
             name=f'{name}|transition_ub',
         )
 
         transition_lower = model.add_constraints(
-            -(continuous_variable.isel({coord: slice(1, None)}) - continuous_variable.isel({coord: slice(None, -1)}))
-            <= max_change * (activate.isel({coord: slice(1, None)}) + deactivate.isel({coord: slice(1, None)})),
+            -(
+                continuous_variable.shift({coord: -1}).isel({coord: slice(None, -1)})
+                - continuous_variable.isel({coord: slice(None, -1)})
+            )
+            <= max_change
+            * (
+                activate.shift({coord: -1}).isel({coord: slice(None, -1)})
+                + deactivate.shift({coord: -1}).isel({coord: slice(None, -1)})
+            ),
             name=f'{name}|transition_lb',
         )
 
@@ -852,10 +894,10 @@ def link_changes_to_level_with_binaries(
 
         # 2. Transition periods: level[t] = level[t-1] + increase[t] - decrease[t] for t > 0
         transition_constraints = model.add_constraints(
-            level_variable.isel({coord: slice(1, None)})
+            level_variable.shift({coord: -1}).isel({coord: slice(None, -1)})
             == level_variable.isel({coord: slice(None, -1)})
-            + increase_variable.isel({coord: slice(1, None)})
-            - decrease_variable.isel({coord: slice(1, None)}),
+            + increase_variable.shift({coord: -1}).isel({coord: slice(None, -1)})
+            - decrease_variable.shift({coord: -1}).isel({coord: slice(None, -1)}),
             name=f'{name}|transitions',
         )
 
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 983681951..ab1cfd97c 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -273,7 +273,7 @@ def _add_scenario_equality_for_parameter_type(
         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)),
+                self.variables[var].isel(scenario=0, drop=True) == self.variables[var].isel(scenario=slice(1, None)),
                 name=f'{var}|scenario_independent',
             )
 
diff --git a/mkdocs.yml b/mkdocs.yml
index 4274cf97b..daf4cf033 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -51,6 +51,7 @@ nav:
       - Troubleshooting: user-guide/troubleshooting.md
       - Community: user-guide/support.md
     - Migration & Updates:
+      - linopy v1 Compatibility: user-guide/linopy-v1-compatibility.md
       - Migration Guide v7: user-guide/migration-guide-v7.md
       - Migration Guide v6: user-guide/migration-guide-v6.md
       - Migration Guide v5: user-guide/migration-guide-v5.md
diff --git a/tests/deprecated/test_flow.py b/tests/deprecated/test_flow.py
index 69922482a..b6c466482 100644
--- a/tests/deprecated/test_flow.py
+++ b/tests/deprecated/test_flow.py
@@ -703,7 +703,7 @@ def test_consecutive_on_hours(self, basic_flow_system_linopy_coords, coords_conf
 
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|forward'],
-            model.variables['Sink(Wärme)|uptime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|uptime'].shift(time=-1).isel(time=slice(None, -1))
             <= model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1)),
         )
@@ -711,10 +711,10 @@ def test_consecutive_on_hours(self, basic_flow_system_linopy_coords, coords_conf
         # eq: duration(t) >= duration(t - 1) + dt(t) + (On(t) - 1) * BIG
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|backward'],
-            model.variables['Sink(Wärme)|uptime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|uptime'].shift(time=-1).isel(time=slice(None, -1))
             >= model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1))
-            + (model.variables['Sink(Wärme)|status'].isel(time=slice(1, None)) - 1) * mega,
+            + (model.variables['Sink(Wärme)|status'].shift(time=-1).isel(time=slice(None, -1)) - 1) * mega,
         )
 
         assert_conequal(
@@ -725,10 +725,10 @@ def test_consecutive_on_hours(self, basic_flow_system_linopy_coords, coords_conf
 
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|lb'],
-            model.variables['Sink(Wärme)|uptime']
+            model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             >= (
                 model.variables['Sink(Wärme)|status'].isel(time=slice(None, -1))
-                - model.variables['Sink(Wärme)|status'].isel(time=slice(1, None))
+                - model.variables['Sink(Wärme)|status'].shift(time=-1).isel(time=slice(None, -1))
             )
             * 2,
         )
@@ -784,7 +784,7 @@ def test_consecutive_on_hours_previous(self, basic_flow_system_linopy_coords, co
 
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|forward'],
-            model.variables['Sink(Wärme)|uptime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|uptime'].shift(time=-1).isel(time=slice(None, -1))
             <= model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1)),
         )
@@ -792,10 +792,10 @@ def test_consecutive_on_hours_previous(self, basic_flow_system_linopy_coords, co
         # eq: duration(t) >= duration(t - 1) + dt(t) + (On(t) - 1) * BIG
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|backward'],
-            model.variables['Sink(Wärme)|uptime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|uptime'].shift(time=-1).isel(time=slice(None, -1))
             >= model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1))
-            + (model.variables['Sink(Wärme)|status'].isel(time=slice(1, None)) - 1) * mega,
+            + (model.variables['Sink(Wärme)|status'].shift(time=-1).isel(time=slice(None, -1)) - 1) * mega,
         )
 
         assert_conequal(
@@ -806,10 +806,10 @@ def test_consecutive_on_hours_previous(self, basic_flow_system_linopy_coords, co
 
         assert_conequal(
             model.constraints['Sink(Wärme)|uptime|lb'],
-            model.variables['Sink(Wärme)|uptime']
+            model.variables['Sink(Wärme)|uptime'].isel(time=slice(None, -1))
             >= (
                 model.variables['Sink(Wärme)|status'].isel(time=slice(None, -1))
-                - model.variables['Sink(Wärme)|status'].isel(time=slice(1, None))
+                - model.variables['Sink(Wärme)|status'].shift(time=-1).isel(time=slice(None, -1))
             )
             * 2,
         )
@@ -869,7 +869,7 @@ def test_consecutive_off_hours(self, basic_flow_system_linopy_coords, coords_con
 
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|forward'],
-            model.variables['Sink(Wärme)|downtime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|downtime'].shift(time=-1).isel(time=slice(None, -1))
             <= model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1)),
         )
@@ -877,10 +877,10 @@ def test_consecutive_off_hours(self, basic_flow_system_linopy_coords, coords_con
         # eq: duration(t) >= duration(t - 1) + dt(t) + (On(t) - 1) * BIG
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|backward'],
-            model.variables['Sink(Wärme)|downtime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|downtime'].shift(time=-1).isel(time=slice(None, -1))
             >= model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1))
-            + (model.variables['Sink(Wärme)|inactive'].isel(time=slice(1, None)) - 1) * mega,
+            + (model.variables['Sink(Wärme)|inactive'].shift(time=-1).isel(time=slice(None, -1)) - 1) * mega,
         )
 
         assert_conequal(
@@ -891,10 +891,10 @@ def test_consecutive_off_hours(self, basic_flow_system_linopy_coords, coords_con
 
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|lb'],
-            model.variables['Sink(Wärme)|downtime']
+            model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             >= (
                 model.variables['Sink(Wärme)|inactive'].isel(time=slice(None, -1))
-                - model.variables['Sink(Wärme)|inactive'].isel(time=slice(1, None))
+                - model.variables['Sink(Wärme)|inactive'].shift(time=-1).isel(time=slice(None, -1))
             )
             * 4,
         )
@@ -952,7 +952,7 @@ def test_consecutive_off_hours_previous(self, basic_flow_system_linopy_coords, c
 
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|forward'],
-            model.variables['Sink(Wärme)|downtime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|downtime'].shift(time=-1).isel(time=slice(None, -1))
             <= model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1)),
         )
@@ -960,10 +960,10 @@ def test_consecutive_off_hours_previous(self, basic_flow_system_linopy_coords, c
         # eq: duration(t) >= duration(t - 1) + dt(t) + (On(t) - 1) * BIG
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|backward'],
-            model.variables['Sink(Wärme)|downtime'].isel(time=slice(1, None))
+            model.variables['Sink(Wärme)|downtime'].shift(time=-1).isel(time=slice(None, -1))
             >= model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             + model.timestep_duration.isel(time=slice(None, -1))
-            + (model.variables['Sink(Wärme)|inactive'].isel(time=slice(1, None)) - 1) * mega,
+            + (model.variables['Sink(Wärme)|inactive'].shift(time=-1).isel(time=slice(None, -1)) - 1) * mega,
         )
 
         assert_conequal(
@@ -974,10 +974,10 @@ def test_consecutive_off_hours_previous(self, basic_flow_system_linopy_coords, c
 
         assert_conequal(
             model.constraints['Sink(Wärme)|downtime|lb'],
-            model.variables['Sink(Wärme)|downtime']
+            model.variables['Sink(Wärme)|downtime'].isel(time=slice(None, -1))
             >= (
                 model.variables['Sink(Wärme)|inactive'].isel(time=slice(None, -1))
-                - model.variables['Sink(Wärme)|inactive'].isel(time=slice(1, None))
+                - model.variables['Sink(Wärme)|inactive'].shift(time=-1).isel(time=slice(None, -1))
             )
             * 4,
         )
diff --git a/tests/deprecated/test_storage.py b/tests/deprecated/test_storage.py
index 3fd47fbf8..4bba38e83 100644
--- a/tests/deprecated/test_storage.py
+++ b/tests/deprecated/test_storage.py
@@ -71,7 +71,7 @@ def test_basic_storage(self, basic_flow_system_linopy_coords, coords_config):
         charge_state = model.variables['TestStorage|charge_state']
         assert_conequal(
             model.constraints['TestStorage|charge_state'],
-            charge_state.isel(time=slice(1, None))
+            charge_state.shift(time=-1).isel(time=slice(None, -1))
             == charge_state.isel(time=slice(None, -1))
             + model.variables['TestStorage(Q_th_in)|flow_rate'] * model.timestep_duration
             - model.variables['TestStorage(Q_th_out)|flow_rate'] * model.timestep_duration,
@@ -154,7 +154,7 @@ def test_lossy_storage(self, basic_flow_system_linopy_coords, coords_config):
 
         assert_conequal(
             model.constraints['TestStorage|charge_state'],
-            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 / eff_discharge * timestep_duration,
@@ -240,7 +240,7 @@ def test_charge_state_bounds(self, basic_flow_system_linopy_coords, coords_confi
         charge_state = model.variables['TestStorage|charge_state']
         assert_conequal(
             model.constraints['TestStorage|charge_state'],
-            charge_state.isel(time=slice(1, None))
+            charge_state.shift(time=-1).isel(time=slice(None, -1))
             == charge_state.isel(time=slice(None, -1))
             + model.variables['TestStorage(Q_th_in)|flow_rate'] * model.timestep_duration
             - model.variables['TestStorage(Q_th_out)|flow_rate'] * model.timestep_duration,
@@ -377,8 +377,8 @@ def test_storage_cyclic_initialization(self, basic_flow_system_linopy_coords, co
         # Check cyclic constraint formulation
         assert_conequal(
             model.constraints['CyclicStorage|initial_charge_state'],
-            model.variables['CyclicStorage|charge_state'].isel(time=0)
-            == model.variables['CyclicStorage|charge_state'].isel(time=-1),
+            model.variables['CyclicStorage|charge_state'].isel(time=0, drop=True)
+            == model.variables['CyclicStorage|charge_state'].isel(time=-1, drop=True),
         )
 
     @pytest.mark.parametrize(
diff --git a/tests/test_clustering/test_expansion_regression.py b/tests/test_clustering/test_expansion_regression.py
index 065c33ddd..a465e6850 100644
--- a/tests/test_clustering/test_expansion_regression.py
+++ b/tests/test_clustering/test_expansion_regression.py
@@ -93,8 +93,13 @@ def test_expanded_storage(self, system_with_storage, solver_fixture):
         fs_e = fs_c.transform.expand()
 
         sol = fs_e.solution
-        # Storage dispatch varies by solver — check charge_state is non-trivial
-        assert float(np.nansum(sol['S|charge_state'].values)) > 0
+        # The inter-cluster SOC level is unpinned (no explicit capacity bound), so its
+        # absolute value is degenerate and solver-dependent — an all-zero SOC is a valid
+        # optimum. Assert the determinate invariants instead: SOC stays finite and within
+        # its non-negative bound, and net discharge balances to ~0. See issue #733.
+        cs = sol['S|charge_state'].values
+        assert np.all(np.isfinite(cs))
+        assert float(np.nanmin(cs)) >= -1e-5
         # Net discharge should be ~0 (balanced storage)
         assert float(np.nansum(sol['S|netto_discharge'].values)) == pytest.approx(0, abs=1e-4)
 
diff --git a/tests/test_math/test_clustering.py b/tests/test_math/test_clustering.py
index aaa37923c..f65df33c9 100644
--- a/tests/test_math/test_clustering.py
+++ b/tests/test_math/test_clustering.py
@@ -340,7 +340,11 @@ def test_storage_cyclic_charge_discharge_pattern(self, optimize):
         assert charge_state.dims == ('cluster', 'time')
         cs_c0 = charge_state.isel(cluster=0).values[:5]
         cs_c1 = charge_state.isel(cluster=1).values[:5]
-        assert_allclose(cs_c0, [50, 50, 0, 50, 0], atol=1e-5)
-        assert_allclose(cs_c1, [100, 100, 50, 100, 50], atol=1e-5)
+        # In 'cyclic' cluster mode the absolute SOC level of each cluster is unpinned
+        # (only the per-step change is fixed by the charge/discharge flows), so the
+        # solution is degenerate in the level and solver-dependent — assert the
+        # determinate charge/discharge dynamics instead. See issue #733.
+        assert_allclose(np.diff(cs_c0), [0, -50, 50, -50], atol=1e-5)
+        assert_allclose(np.diff(cs_c1), [0, -50, 50, -50], atol=1e-5)
 
         assert_allclose(fs.solution['objective'].item(), 100.0, rtol=1e-5)