refac: Refactor piecewise API to a stateless construction layer + tangent_lines utility

[FBumann]

Summary

Follows up on the discussion in #602. Replaces the descriptor pattern (PiecewiseExpression, PiecewiseConstraintDescriptor, piecewise()) with a stateless construction layer. Establishes a clean separation between piecewise equality constraints and mathematically very different tangent-line inequality bounds.

Final API

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

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

# With breakpoints() factory
m.add_piecewise_constraints(
    (power, linopy.breakpoints([0, 30, 60, 100])),
    (fuel,  linopy.breakpoints([0, 36, 84, 170])),
)

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

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

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

# Per-entity breakpoints (different curves per generator)
m.add_piecewise_constraints(
    (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")),
)

# Unit commitment with active parameter
m.add_piecewise_constraints(
    (power, x_pts),
    (fuel,  y_pts),
    active=commit,  # binary variable gates the PWL
)

# Inequality bounds (pure LP, no aux variables) — separate utility
t = linopy.tangent_lines(power, x_pts, y_pts)
m.add_constraints(fuel <= t)   # upper bound (concave f)
m.add_constraints(fuel >= t)   # lower bound (convex f)

Key changes

• Remove piecewise() function and descriptor classes — no intermediate state objects
• Tuple-based add_piecewise_constraints — each (expression, breakpoints) pair is a positional arg. 2-var and N-var are the same pattern.
• tangent_lines() utility — computes tangent-line LinearExpression per segment. Pure LP, no auxiliary variables. Separate from piecewise equality.
• sign parameter removed — piecewise only does equality linking. Inequality via tangent_lines + regular add_constraints.
• N-variable support restored — link 3+ expressions through shared lambda weights (CHP use case from feat: add piecewise linear constraint API (SOS2, incremental, disjunctive) #576)
• Variable names as link coords — the _pwl_var dimension shows power, fuel, heat instead of 0, 1, 2
• Clean operator overloads — Variable.__le__, __ge__, __eq__ simplified (no more piecewise dispatch)

Design principles

• API surface stays minimal: Model, Variables, Constraints, Expression — no new user-facing types
• Piecewise is a construction layer that produces regular linopy objects
• Equality linking (core formulation) is cleanly separated from inequality bounds (tangent lines)
• tangent_lines creates zero variables — it returns a LinearExpression with one tangent line per segment

Notebook examples

#	Example	Feature demonstrated
1	Gas turbine	SOS2 formulation
2	Coal plant	Incremental formulation
3	Diesel generator	Disjunctive (disconnected segments)
4	Concave bound	tangent_lines for inequality
5	Slopes mode	breakpoints(slopes=...)
6	Unit commitment	active parameter
7	CHP plant	N-variable linking (3 vars)
8	Generator fleet	Per-entity breakpoints

Other API designs considered

A) Descriptor pattern (PR #602, current 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_constraints(
    exprs={"power": power, "fuel": fuel, "heat": heat},
    breakpoints=bp,
)

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

C) Keyword-only 2-variable + dict N-variable (two separate forms)

m.add_piecewise_constraints(x=power, y=fuel, x_points=xp, y_points=yp)
m.add_piecewise_constraints(exprs={...}, breakpoints=bp)

Considered: Clear but maintains two separate code paths. The 2-var form doesn't generalize.

D) List of expressions + breakpoints array (name-matching)

m.add_piecewise_constraints(expressions=[power, fuel], breakpoints=bp)

Rejected: Magic name matching — breaks if variable names don't match breakpoint coordinates.

E) Chosen: Tuple-based

m.add_piecewise_constraints(
    (power, [0, 30, 60, 100]),
    (fuel,  [0, 40, 85, 160]),
)

Selected: Each pair collocates expression + breakpoints. No string keys, no separate dict, no name matching. 2-var and N-var are the same pattern. Minimal, explicit, composable.

Test plan

• All 107 piecewise tests pass
• Full test suite passes
• Notebook executes cleanly
• Ruff lint + format clean
• Review with @FabianHofmann and @coroa

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