feat: Unify piecewise API behind add_piecewise_formulation (sign + LP dispatch)

[FBumann]

Summary

Follows up on the discussion in #602. Replaces the descriptor pattern (PiecewiseExpression, PiecewiseConstraintDescriptor, piecewise()) with a single, stateless construction layer: add_piecewise_formulation.

add_piecewise_formulation becomes the one entry point for piecewise work — equality, one-sided inequality, N-variable linking, per-entity breakpoints, unit-commitment gating, and automatic method dispatch all live behind one call. tangent_lines() survives as a low-level helper for users who want to build chord expressions manually.

The API

Each (expression, breakpoints) tuple links a variable to its breakpoints. All tuples share interpolation weights, coupling them on the same curve segment.

Equality (default)

# 2-variable: fuel = f(power)
m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel,  [0, 36, 84, 170]),
)

# N-variable (CHP plant — same pattern, more tuples)
m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel,  [0, 40, 85, 160]),
    (heat,  [0, 25, 55, 95]),
)

# Disjunctive (disconnected operating regions)
m.add_piecewise_formulation(
    (power, linopy.segments([(0, 0), (50, 80)])),
    (cost,  linopy.segments([(0, 0), (125, 200)])),
)

# Per-entity breakpoints (different curves per generator)
m.add_piecewise_formulation(
    (power, linopy.breakpoints({"gas": [0, 30, 60, 100], "coal": [0, 50, 100, 150]}, dim="gen")),
    (fuel,  linopy.breakpoints({"gas": [0, 40, 90, 180], "coal": [0, 55, 130, 225]}, dim="gen")),
)

# Breakpoints from slopes
m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel,  linopy.breakpoints(slopes=[1.2, 1.4, 1.7], x_points=[0, 30, 60, 100], y0=0)),
)

Inequality bounds — sign="<=" / ">="

The sign parameter applies the bound to the first tuple (first-tuple convention); remaining tuples are inputs pinned to the curve.

# fuel ≤ f(power).  "auto" picks the cheapest correct formulation — a pure
# LP chord formulation when the curvature matches the sign (concave + "<=",
# or convex + ">="), SOS2/incremental otherwise.
m.add_piecewise_formulation(
    (fuel,  [0, 20, 30, 35]),   # bounded output listed FIRST
    (power, [0, 10, 20, 30]),   # input, always equality
    sign="<=",
)

On a matching-curvature curve this dispatches to method="lp": one chord inequality per segment plus a domain bound x_min ≤ x ≤ x_max, and no auxiliary variables. Mismatched curvature (convex + "<=", or concave + ">=") describes the wrong region (not just a looser one) — method="auto" detects this and falls back to SOS2/incremental with an explanatory info log; method="lp" raises with a clear error.

Inequalities with 3+ variables cant go lp. They are always sos/...

Unit-commitment gating — active

# active=1: power in [30, 100], fuel = f(power)
# active=0: power = 0, fuel = 0
m.add_piecewise_formulation(
    (power, [30, 60, 100]),
    (fuel,  [40, 90, 170]),
    active=commit,  # binary variable
)

Works with SOS2, incremental, and disjunctive. With sign != "==", deactivation only pushes the signed bound to 0 — set lower=0 on naturally non-negative outputs (fuel, cost, heat) and the combination forces y = 0 automatically.

Low-level: tangent_lines

# Returns a LinearExpression with one chord per segment — no variables.
# Most users should prefer add_piecewise_formulation(..., sign="<="),
# which builds on this helper and adds domain bounds + curvature checks.
t = linopy.tangent_lines(power, x_pts, y_pts)
m.add_constraints(fuel <= t)

Key changes

New public API

• add_piecewise_formulation — one entry point for all piecewise constructions. Returns a PiecewiseFormulation object carrying .method (resolved formulation), .convexity (when well-defined), .variables, .constraints.
• sign parameter — "==" (default) / "<=" / ">=", with the first-tuple convention.
• method parameter — "auto" (default), "sos2", "incremental", "lp".
• active parameter — binary gating for unit commitment.
• linopy.breakpoints() / linopy.segments() factories.
• linopy.slopes_to_points() utility.
• linopy.tangent_lines() low-level helper.

Auto-dispatch (method="auto")

1. 2-variable inequality on a matching-curvature curve → lp (chord + domain bounds, no aux vars)
2. Strictly monotonic breakpoints → incremental (deltas + binaries)
3. Otherwise → sos2 (lambdas + SOS2 constraint)
4. Disjunctive (segments) → sos2 with binary segment selection

The resolved method is logged at INFO; when LP is skipped, the log explains why (curvature, NaN layout, tuple count, or active).

Correctness properties

• _detect_convexity is direction-invariant — ascending and descending x give the same label for the same graph.
• NaN-padded per-entity breakpoints are masked everywhere, including the LP chord constraints (padded segments don't create spurious y ≤ 0 constraints).
• method="lp" validates curvature + sign + tuple count + active absence upfront; mismatch raises with a clear fix pointer to method="auto".
• PiecewiseFormulation.method and .convexity persist across netCDF round-trip.

Architecture (internal)

• _PwlLinks dataclass carries the equality-side and signed-side link expressions; _build_links() is the single place where the sign split happens.
• _try_lp() / _resolve_sos2_vs_incremental() flatten the dispatch — _add_continuous is ~10 lines of real logic.
• _add_sos2 / _add_incremental / _add_disjunctive take 4–5 arguments (down from 9–11).
• _add_lp delegates chord construction to the public tangent_lines helper.

Removed / deprecated

• PiecewiseExpression, PiecewiseConstraintDescriptor, linopy.piecewise() — descriptor pattern gone.
• Old add_piecewise_constraints shape that dispatched on Variable.__le__/__ge__/__eq__ return types.
• Variable operator overloads (__le__, __ge__, __eq__) simplified — no more piecewise dispatch.

Design principles

• API surface stays minimal: Model, Variables, Constraints, Expression, plus the single PiecewiseFormulation return type. No user-facing intermediate state.
• Piecewise is a construction layer that produces regular linopy objects.
• One entry point for equality and inequality — sign selects the semantics, method selects the cost.
• Variable names as link coords: the internal _pwl_var dimension shows power, fuel, heat instead of 0, 1, 2.

Alternative designs considered

A) Descriptor pattern (PR #602, previous master)

pw = linopy.piecewise(x, x_pts, y_pts)
m.add_piecewise_constraints(pw == y)

Rejected: Introduces PiecewiseExpression and PiecewiseConstraintDescriptor as user-facing state. Breaks the "only Variables, Constraints, Expressions" principle. Structurally limited to 2-variable x → y.

B) Dict of expressions + shared breakpoints DataArray

m.add_piecewise_formulation(exprs={"power": power, "fuel": fuel}, breakpoints=bp)

Considered: Supports N-variable. But requires string keys to match breakpoint coordinates — an indirection layer.

C) Per-tuple signs

m.add_piecewise_formulation((power, xp, "<="), (fuel, yp, "<="))

Rejected: Creates ambiguous cases (what if tuples disagree on sign? what if only some have it?). Global sign= + the first-tuple convention keeps the math clean.

D) Last-tuple convention (output listed last)

m.add_piecewise_formulation((power, xp), (fuel, yp), sign="<=")  # fuel bounded

Rejected: The y = f(x) convention reads left-to-right; listing the bounded output first ("first expression ≤ f(the rest)") is more natural. First-tuple convention chosen after discussion.

E) Chosen: tuple-based, global sign, first-tuple convention

m.add_piecewise_formulation(
    (fuel,  y_pts),    # bounded output, listed first
    (power, x_pts),    # input
    sign="<=",
)

Notebook examples

examples/piecewise-linear-constraints.ipynb (main tutorial, 7 sections):

#	Section	Feature
1	Getting started	Basic 2-variable equality + method inspection via pwf
2	Picking a method	auto vs sos2 vs incremental give the same optimum
3	Disjunctive segments	segments() for gaps in operation
4	Inequality bounds	sign="<=" (auto-dispatches LP)
5	Unit commitment	active parameter
6	N-variable linking	CHP plant (3 variables)
7	Per-entity breakpoints	Fleet with different curves

examples/piecewise-inequality-bounds.ipynb (dedicated deep-dive on sign / LP):

• Mathematical formulation (equality, inequality, LP, incremental)
• 2D hypograph feasibility plot
• 3D ribbon for 3-variable case (multiple viewpoints)
• Convexity × sign rule table, "wrong region" analysis
• Auto-dispatch fallback demonstration

Test plan

• All 148 piecewise tests pass (incl. 11 direct unit tests on _detect_convexity)
• Full test suite passes
• Both notebooks execute cleanly
• Ruff lint + format clean, mypy clean on touched files
• Round-trip tested: PiecewiseFormulation.method and .convexity persist across to_netcdf / read_netcdf
• Review

🤖 Generated with Claude Code

[FBumann]

@FabianHofmann This is pretty urgent to me, as the current piecewise API is already on master and should NOT be included in the next release in my opinion.
But i understand that the CSR PR is more important atm.

The current code can be reorganized a bit, but id like to hear your thoughts about the API first.

[FBumann]

edit: did some refactoring anyway, grouped in #641

[FBumann]

N-variable disjunctive + follow-up notes

Done: Refactored _add_disjunctive to use the same stacked N-variable pattern as _add_continuous. The 2-variable restriction is removed — disjunctive now works with any number of (expression, breakpoints) pairs, using a single _link constraint with a _pwl_var dimension (same as continuous SOS2/incremental).

[FBumann]

Additional follow-ups

Single-variable support: All formulation methods (SOS2, incremental, disjunctive) should support a single (expression, breakpoints) pair. Use case: constraining a variable to lie on a set of discrete segments or breakpoints without linking to another variable. The stacking logic should handle len(pairs) == 1 naturally — just relax the len(pairs) < 2 check.

Segments with 2+ breakpoints: Currently segments() only supports pairs of points (start, end) per segment. It should support segments with 2 or more breakpoints, enabling piecewise-within-segment curves (e.g. a nonlinear segment approximated by multiple linear pieces within each disconnected region. Currently, a user could express this with "touching" segments, but its mathematically less efficient).

[FBumann]

Comparison with JuMP's piecewise linear formulations

JuMP's approach (via PiecewiseLinearOpt.jl)

JuMP's core piecewiselinear() returns a single variable representing the PWL output. All auxiliary variables/constraints are added internally. The package offers 8 formulation methods that all share the same convex-combination (λ) core — they differ only in how SOS2 adjacency is enforced:

Method	Binary vars	Notes	linopy
:SOS2	0	Delegates to solver's native SOS2 branching	✅ _add_sos2
:CC (Convex Combination)	O(n)	Same as SOS2 but adjacency enforced with explicit binaries instead of SOS2 set	✅ same formulation as _add_sos2but with our internal sos_reformulation
:MC (Multiple Choice)	O(n)	Disaggregated segment copies with binary selection	✅ essentially _add_disjunctive — and we extend it with gap support via segments()
:Incremental	O(n)	Delta chaining	✅ _add_incremental
:Logarithmic (default)	O(log n)	Gray code encoding	❌ not yet
:DisaggLogarithmic	O(log n)	Disaggregated + Gray codes	❌ not yet
:ZigZag	O(log n)	Zig-zag codes	❌ not yet
:ZigZagInteger	1 general int	Most compact	❌ not yet

The default is :Logarithmic — ceil(log2(n-1)) binary variables, significantly more efficient for branch-and-bound with many breakpoints.

JuMP also supports bivariate PWL via triangulation (UnionJack, K1, Stencil, BestFit patterns) for z = f(x, y) with independent inputs.

Reference: Huchette & Vielma, "Nonconvex piecewise linear functions: Advanced formulations and simple modeling tools" (2017).

Two separate concerns: piecewise structure vs. SOS2 enforcement

JuMP mixes these into a single method parameter. We can do better by separating them:

Piecewise formulation structure (what add_piecewise_formulation handles):

Method	What it does	linopy
SOS2	λ weights + sum=1 + linking constraints. Delegates adjacency enforcement to add_sos_constraints.	✅ _add_sos2
Incremental	Completely different structure — δ deltas + binary chaining. No λ, no SOS2 at all.	✅ _add_incremental
Disjunctive	Per-segment λ + binary segment selection. Uses SOS2 within each segment, but segment decomposition is piecewise logic.	✅ _add_disjunctive

SOS2 adjacency enforcement (what add_sos_constraints / sos_reformulation.py handles — orthogonal to piecewise):

Method	Binary vars	Notes	linopy
Native	0	Solver handles branching	✅ add_sos_constraints
CC / Big-M	O(n)	x[i] <= M[i] * (z[i-1] + z[i])	✅ sos_reformulation.py (auto via reformulate_sos=True)
Logarithmic (Gray code)	O(log n)	Best trade-off for many breakpoints	❌ not yet, high value follow op ⭐
DisaggLogarithmic	O(log n)	Disaggregated + Gray codes	❌ not yet
ZigZag	O(log n)	Different code pattern, similar performance	❌ not yet
ZigZagInteger	1 general int	Most compact	❌ not yet

This separation means:

• Piecewise stays simple: 3 structural methods, already complete
• SOS2 reformulations benefit all SOS constraints in the model, not just piecewise
• The two choices compose: e.g. disjunctive piecewise + logarithmic SOS2 enforcement

JuMP bundles all 8 as piecewise methods. We separate structure from enforcement — cleaner architecture, broader applicability.

What linopy has on top of JuMP

Feature	Description
Disconnected segments	segments() factory + _add_disjunctive for forbidden operating zones with gaps (e.g. off or 50–80 MW). JuMP's MC handles contiguous segments only.
N-variable linking	Link 3+ expressions through shared λ weights (CHP: power/fuel/heat). JuMP's piecewiselinear() is 1→1 only.
active parameter (unit commitment)	Binary variable gates the entire PWL — when active=0, all aux variables are zero. In JuMP this requires manual formulation.
tangent_lines() utility	Pure LP inequality bounds (no auxiliary variables). Returns a LinearExpression per segment for use with regular add_constraints.
Per-entity breakpoints	Different curves per generator via breakpoints({"gas": [...], "coal": [...]}, dim="gen") with automatic xarray broadcasting.
Auto method selection	method="auto" detects monotonicity and chooses incremental (faster) vs SOS2 (flexible). JuMP requires explicit method choice.

Why our architecture doesn't block any of these

The design in this PR is method-agnostic by construction:

1. The dispatch is a simple if/else on a string in _add_continuous. The 3 structural methods are already complete.

2. All formulation methods reduce to the same primitives: add_variables(continuous/binary/integer) + add_constraints + optionally add_sos_constraints. Every JuMP method uses exactly these building blocks.

3. Result tracking is formulation-blind: PiecewiseFormulation snapshots which variable/constraint names were added via before/after diffing. It doesn't care how they were structured — any new method's artifacts get captured automatically.

4. SOS2 reformulations are additive: Adding logarithmic/zigzag to sos_reformulation.py requires zero changes to the piecewise layer — it just gets faster SOS2 enforcement for free.

5. No formulation-specific assumptions leak into the public API. The user writes method="sos2" and gets back the same PiecewiseFormulation with .variables and .constraints, regardless of how SOS2 is enforced.

Our N-variable linking vs. JuMP's triangulation

These solve different problems (complementary, not competing):

	Our N-variable (1D curve)	Triangulation (2D+ mesh)
Inputs	1 hidden parameter (e.g. load)	2+ independent variables (x, y)
Breakpoints	1D sequence	2D grid, triangulated into simplices
λ weights	2 adjacent points (SOS2 along 1 axis)	3 vertices of one triangle
Degrees of freedom	1	2+
Use case	Coupled outputs of one parameter (CHP, efficiency curves)	z = f(x, y) with independent inputs

Our 3-variable (power, fuel, heat) links three outputs coupled by a single operating parameter — you can't set power=60 and heat=25 independently; the curve determines heat given power. Triangulation handles z = f(x, y) where x and y move independently (e.g., cost as a function of temperature AND load).

Follow up: Triangulation

Triangulation is mathematically quite a bit more complex and needs different data. It takes a mesh instead of a 1D breakpoint sequence. (mesh). Its clearly a follow up, probably as a new Model.add_*()* method that can reuse the PiecewiseFormulation return type, name tracking, and model primitives.

[FBumann]

@coroa Can we discuss this together with @FabianHofmann ?

[coroa]

@coroa Can we discuss this together with @FabianHofmann ?

@FabianHofmann is on vacation this week. Can we set up something next week? Or is pressing more urgent than that?

[FBumann]

@coroa Can we discuss this together with @FabianHofmann ?

@FabianHofmann is on vacation this week. Can we set up something next week? Or is pressing more urgent than that?

Ah ok i missed that. No, next week is fine :) Could you setup a meeting?

[FBumann]

From call:

• Add a sign parameter to add_piecewise_constraints()
• add_piecewise_constraints() should forward to tangent lines if possible
• Add new Slopes object, drop slopes support from breakpoints
• Maybe add a reference variable in add_piecewise_constraints() so other variables/slopes can refer to it (especially for tangent lines)

[FBumann]

@coroa I thought a lot about the inequality in non lp piecewise relationships.
I added a notebook and plot about it in #662 (5992788). Essentially, the sign needs to be "translated" to behave like tangent lines. And we need to make a decision about the convention, because givin all options to the user makes a pretty bad UX whie "allowing" cases which probably noone wants.
Here is a visualization for feasible regions for a 2-variable case with SOS2/lambda/delta formulations based on which sign is used for which variable:

Compared to tangent line one:

We can see that we need to fix one of the variables to the slope! TO be precise, all but one variable needs to be fixed!

[FBumann]

@coroa What do you think about the new api?
I think its quite good and clear now. The testing is also already quite good imo.

One thing we should discuss is the notebooks and docs. If you have the time to read it, feedback and potentially changes would be great!
The one who implemented it usually has a different angle on the docs than the reviewer/user.
https://linopy--638.org.readthedocs.build/en/638/piecewise-linear-constraints.html

[FBumann]

@coroa One API detail im not sure about yet, is where to but the sign.
It could be per tuple as a third optional parameter:

m.add_piecewise_formulation(
   (power, xp, "<="),
   (fuel, yp)
)

with defaulting to '==' where not specified.

I decided against this api for simplicity, as using more than i single equality produces mathematically weird bounds as far as i understand it (only relevant for the 3 var case). ALso, for the 2 variable LP case, it doestn make any sense to have more than a single sign. So there is a lot of API surface that would lead to errors we need to catch.

We could however also just disallow the 3 variable case with inequalities for now, as i don`t understand the use case that well. This would allow the arguably more obvious notation of having the sing right next to the expression...

With the EvolvingAPIWarning however, im fine with keeping it and adjusting it if needed.

Looking forward, this API seems better suited for triangular bounding, where we could actually allow multiple inequalities. But the implementation of this is not part of this PR and quite complex i think, so we can just change it as soon as this potential feature lands

[FBumann]

One thing I didn't cover yet, is a dedicated slopes object. Should we add it? Or can we do it in a follow up?