Design: Model.restrict / advance — slice a model along a coordinate and condition on boundary values (rolling horizon) #768

Describe the feature you'd like to see

A principled API for slicing a model along a coordinate dimension and conditioning it on fixed boundary values — making rolling-horizon optimization a pure linopy pattern:

rm = m.restrict(snapshot=window)                    # pure slice; raises on dangling boundary terms
rm = m.restrict(snapshot=window, fix=prev_solution) # slice + fold out-of-window terms into rhs

rm.advance(snapshot=next_window, fix=prev_solution) # move the window in place (values only)

Motivation

Rolling-horizon workflows (e.g. PyPSA) repeatedly solve overlapping windows of one long-horizon problem. With the persistent solver interface (#718) and the in-place update API (#727), a prototype on a PyPSA-Eur 2040 network (426k vars / 1.28M cons per 168h window, stepping 24h) reaches a 1.76x end-to-end speedup on warm iterations with bitwise-identical objectives vs rebuilding the model each window — but currently requires hand-written, formulation-specific update code. restrict/advance would make this generic.

The mathematical operation is: restrict the model to the window's variables/constraints, and for constraint rows referencing variables outside the window (e.g. storage energy balance at the window's first snapshot referencing the previous snapshot's state of charge), fold those terms into the rhs: rhs -= coeff * fixed_value.

Semantics

Model.restrict(fix=None, **selection) — creates structure:

• slices every container along the selected dims; containers without the dim are kept whole
• detects dangling boundary terms (terms whose variable labels fall outside the selection — cheap on the flat/CSR representation)
• without fix: raises, listing the offending constraints (fail fast — pure selection is just the degenerate case where nothing dangles)
• with fix: folds dangling terms into the rhs using the actual matrix coefficient, then masks them
• returns a window-sized model (own labels, own solve)

restricted.advance(fix=..., **selection) — mutates values only:

• shifts the coordinate (cf. Model.assign_coords, Model.assign_coords: coordinate reassignment on an existing model #767), overwrites rhs/bounds/objective coefficients from the parent slice, re-folds the boundary rhs with the new fix values
• guarantees labels, active set, container order, and sparsity stay identical — so solver.update(..., ignore_dims={...}) stays on the in-place path (no rebuild)
• raises if the new window is not structurally identical (caller falls back to restrict + rebuild)

In short: restrict is the constructor (rebuild-scale, once), advance is the mutation (sub-second, every iteration).

Why this shape (alternatives ruled out)

The solver column/row space is defined by the dense positions of active labels; the persistent diff forces a full rebuild whenever the active set changes (STRUCTURAL_LABELS) and tolerates only coordinate shifts under ignore_dims plus value changes. This rules out:

• Full-size view, out-of-window variables fixed via bounds (the Pyomo/JuMP/gurobipy pattern): right update protocol, wrong size. Micro-experiment: HiGHS presolve does not collapse a mostly-fixed model — a full model with 90% of columns fixed solved ~13–16x slower than a genuinely window-sized LP.
• Materialized sub-model per window: relabeling triggers STRUCTURAL_LABELS every window — "rebuild, but smaller".
• Mask-shifting (set out-of-window labels to -1 per window): same rebuild trigger; would need new add/remove-row diff capability against a live solver (position renumbering on delete makes this a correctness minefield).

The only shape that gets both window-sized solve cost and in-place updates is a window-sized active set whose coordinate shifts under ignore_dims, with boundary terms folded into the rhs — exactly the prototype's pattern, promoted to a first-class API. Notably, no ecosystem peer auto-folds boundary terms; they all rely on fixed columns.

Design notes / open questions

1. Not Model.sel: .sel everywhere in xarray and linopy means pure selection — result data are unchanged subsets of the parent. Folding rewrites rhs and yields a different optimization problem, so it gets a different verb. Pure selection is the no-dangling-terms degenerate case of restrict; a separate strict Model.sel is unnecessary.
2. fix naming clash — resolved via Variable.fix: move from equality-constraint to bound-collapse semantics #769: Variable.fix() moves to bound-collapse semantics, so "fix" means hold a variable at a value everywhere — bound collapse in a live model, fold-out (elimination, the strongest form) in restrict/advance. Renaming the kwarg instead was ruled out: boundary= is wrong for the Benders use (see decomposition section below), and fixed_values= would concede the standard term to the non-standard equality-constraint semantics.
3. Folding coefficient semantics: generic folding uses the real matrix coefficient. PyPSA's own per-window builds fold the storage initial state with coefficient 1 instead of the standing-loss coefficient used on interior rows — so restrict reproduces a slightly different (arguably more correct) model than PyPSA's hand-rolled rolling horizon. Must be documented, not mimicked.
4. Fold order: fold before masking/sanitizing (the fold needs the original coefficient and the fixed value; sanitize_zeros must not null coefficients first). Folding must not introduce NaN rhs (auto_mask treats NaN as a mask trigger).
5. Solution write-back: the restricted model carries the current window's coords, so assign_result lands on the right index per window; mapping back to the parent horizon is the caller's concern (or a later convenience).

Relation to decomposition methods

restrict(fix=...) is not just a rolling-horizon helper — algebraically it is the subproblem-conditioning primitive of decomposition methods. Partition variables into in-window $x_W$ and out-of-window $x_B$; folding at the actual matrix coefficient forms

$$\min ; c_W^\top x_W \quad \text{s.t.} \quad A_{WW}, x_W \sim b_W - A_{WB}, \bar{x}_B$$

which is exactly the Benders subproblem at a fixed complicating point $\bar x_B$. And advance — structure frozen, only values change, persistent solver stays on the in-place path with dual-simplex warm starts — is exactly the per-iteration contract of every cutting-plane scheme. The methods differ only in what drives the fix values:

Rolling horizon is literally the cut-free forward pass of SDDP. The rhs-side methods fit restrict/advance directly; the objective-side methods are already served by the update API (#727) and need nothing from restrict.

Decomposition algorithms themselves (master problem, cut management, feasibility rays, convergence loops) stay out of linopy. But three design decisions are cheap now and breaking later, and make restrict the building block for user-side Benders/SDDP:

1. fix folds out any covered variable, not only slicing-induced dangling terms. The canonical capacity-expansion split fixes investment variables (no snapshot dim) — under the "containers without the dim are kept whole" rule they would stay in the subproblem as free columns. Specifying fix as "variables in fix are folded out at these values" (a superset of the dangling rule, same fold code path) makes m.restrict(snapshot=horizon, fix=capacities) the Benders operational subproblem, even if v1 only implements the dangling case.
2. Keep the fold ledger. The cut gradient w.r.t. the fixed values is $-A_{WB}^\top \lambda$ — boundary-row duals times the folded coefficients, i.e. exactly the (row, fixed_var, coeff) triples the fold computes (and which advance needs anyway for re-folding). Exposing them after solve (e.g. rm.fixed_sensitivity, a Series over fixed-variable labels) makes user-side Benders ~20 lines of linopy; without it, users must re-derive $A_{WB}$ from the parent model.
3. advance with selection optional: a pure re-fix at unchanged coords (rm.advance(fix=new_master_sol)) must be legal — that is the Benders iteration.

This framing sharpens two open questions above: the fix naming clash (point 2) is resolved on its own terms in #769 — "fixing complicating variables" is the standard decomposition term, and boundary= would be wrong for the Benders use; and "folding uses the real matrix coefficient" (point 3) is not a PyPSA-compatibility footnote but a correctness requirement — cut gradients computed from anything else would be wrong.

Relation to other issues

• Model.assign_coords: coordinate reassignment on an existing model #767 (Model.assign_coords) is the standalone primitive advance builds on.
• Variable.fix: move from equality-constraint to bound-collapse semantics #769 (Variable.fix → bound-collapse semantics) clears the fix= kwarg naming clash (open question 2).
• feat(persistent): in-place solver updates (Solver framework + HiGHS/Gurobi/Xpress/Mosek) #718 / Variable.update() / Constraint.update() as canonical mutation API #727 provide the persistent diff/update machinery this rides on — advance needs no new diff features.

Prototype and benchmark: dev-scripts/persistent-solver/ on the feat/solver-update branch (inplace_update.py, pypsa-rh-eur-inplace.py, verify-window-update.py).