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

[FabianHofmann]

Note

AI assisted writing

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:

Method	Subproblem per iteration	What changes
Rolling horizon	restrict(snapshot=W, fix=prev)	rhs (boundary state) — myopic, no feedback
Benders / L-shaped	restrict(..., fix=master_sol)	rhs — master feeds back via cuts
SDDP	rolling forward pass + Benders backward pass	rhs both ways
Lagrangian / DW pricing	block subproblem	objective coeffs (multipliers/prices)

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).

[FBumann]

Some thoughts of me through discussion with AI. Maybe I got something wrong, but I think naming and api are important here.

Note

AI generated

On naming and API shape

fold is more accurate than restrict — it names the load-bearing operation (absorbing boundary terms into the RHS) rather than the slicing, which users already understand from .sel. refold is the natural companion: same operation, new values, cheaper execution.

The two-method split reflects an implementation boundary, not a conceptual one. From the user’s perspective both calls express identical intent — “give me this window with these boundary values folded in.” The difference is purely whether the solver can warm-start.

That suggests a single fold that detects structural identity internally and takes the appropriate path:

rm = m.fold(snapshot=window, fix=prev_solution) # first call: builds structure
rm = m.fold(snapshot=next_window, fix=new_solution) # subsequent: mutates in place

The performance contract — first call expensive, subsequent calls cheap — is surfaced not through a separate method name but through a property:

rm.folded # True if the last call took the in-place path

This keeps the API to one concept and one method, makes the common rolling-horizon loop uniform across iterations, and leaves the solver optimization as an implementation detail where it belongs.

[FabianHofmann]

I like the name fold but I would not apply it to both kinds of backends but rather keep the distinction. One is cold-starting the other is guaranteed to warm-start. Both stay lean without auto-detection on identity. Users decide where to place which.

We could rename to

1. fold == always cold-start and (slow, creates new model instance)
2. refold == always warm-start (fast, repurpose initialized model)

Also, the two functions will have different return types. fold returns a new model instance; refold modifies inplace and returns None.

[FBumann]

True. Then we definitely should keep 2 methods.

[FBumann]

I think this is an amazing addition in general, really leveraging the dimension aware foundation of linopy to enable things much more complex in other tools

[odow]

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

Do you have the instance? The HiGHS folks would be interested

[coroa]

Something that does not show up here yet, but is important for later use:

in decomposition problems as you noted above you will have multiple restricted/folded models at the same time that you want to keep persisted and update. So you want the folding process to be able to split one linopy model into multiple ones living in separate threads (or at least the solver processes living in separate threads).

Maybe as a simple example think capacity variables in one master and two dispatch blocks for 2-weeks in the sub problem in benders. you would want to be able to advance/refold the subproblems separately

[coroa]

The OptiGraph/Plasmo.jl framework allows effectively labelling parts of the optimization problem into different nodes and then use graph partitioning algorithms to split it up. In the linopy space this would be something like a map on the variables to optigraph node names, but this might be a diversion from what is already proposed here. https://arxiv.org/abs/2501.02098 and https://arxiv.org/abs/2511.14966 for how it connects to multi machine stuff.

[FBumann]

My first intuition would be to not pull such things into linopy, but rather provide the plumbing to enable this.

Like with a sync solving, we only provide the primitives to start multiple solves in different threads, but opening those threads is the users responsibility, as there is not one right way to do it.

[FabianHofmann]

Ideally we generalize the fold/refold enough to allow fancy stuff like graph partitioning. If we need another data handling feature in linopy that has a clear scope and is not too opinionated, I would be down for it. Probably good to align with the jump approach and move diverging things into standalone packages.