feat(piecewise): add sign parameter and LP method to add_piecewise_formulation

[FBumann]

Draft — targets feat/piecewise-api-refactor (PR #638). Merge into that branch, not master.

What this PR adds on top of #638

PR #638 establishes add_piecewise_formulation as the stateless construction entry point, covering equality formulations (SOS2 / incremental / disjunctive). Inequality bounds are left to the user via the bare tangent_lines helper, which returns chord expressions but adds no variables, no domain bounds, and no convexity check.

This PR closes those gaps by folding inequality support into add_piecewise_formulation itself:

API additions

• sign: Literal["==", "<=", ">="], default "=" — selects equality (current behaviour) or a one-sided bound.
• method: Literal["sos2", "incremental", "auto", "lp"] — new "lp" method uses tangent lines directly, no auxiliary variables.

Semantic additions

• First-tuple convention: when sign != "=", the first tuple is the output bounded by the sign; all other tuples are inputs forced to equality on the curve. Only one variable carries the sign — prevents ambiguous cases (dominated region, mixed per-variable signs) from slipping in through the main API.
• Automatic domain bounds (x_0 ≤ x ≤ x_n) added by the LP path. tangent_lines alone doesn't add these, and chord extrapolation past the breakpoints produces surprising feasible regions; the LP method in add_piecewise_formulation closes that footgun.
• Convexity check (_detect_convexity) — classifies the curve as convex / concave / linear / mixed and requires the sign to match curvature for LP to be a valid (tight) bound. Mismatched combinations raise rather than silently producing a wrong feasible region.
• Auto-dispatch: method="auto" picks "lp" for 2-variable inequalities when curvature matches sign (concave + <= or convex + >=). Otherwise falls back to SOS2 (non-monotonic / mixed) or incremental (monotonic), always with the sign applied to the output link.

Implementation details

• Split link constraint only when sign != "=". Equality keeps the single stacked link as in feat: Unify piecewise API behind add_piecewise_formulation (sign + LP dispatch) #638 (no regression). Inequality splits into one stacked equality link for inputs plus one single-row signed link for the output.
• New constraint suffixes: PWL_LP_SUFFIX, PWL_DOMAIN_LO_SUFFIX, PWL_DOMAIN_HI_SUFFIX, PWL_OUTPUT_LINK_SUFFIX.
• Uses existing SIGNS / EQUAL / LESS_EQUAL / GREATER_EQUAL from linopy.constants, plus sign_replace_dict to normalise "==" → "=".

API before vs after

# --- PR #638 state ---

# Equality
m.add_piecewise_formulation((x, xp), (y, yp))

# Inequality: user assembles manually; no domain bounds; no convexity check
t = linopy.tangent_lines(x, xp, yp)        # bare chord expressions
m.add_constraints(y <= t)                  # extrapolates past [x_0, x_n]
# (user is on their own if convexity is wrong — silently gets a loose/wrong region)


# --- This PR ---

# Equality (unchanged)
m.add_piecewise_formulation((x, xp), (y, yp))

# Inequality: single call, auto dispatch, domain bounds + convexity check
m.add_piecewise_formulation((y, yp), (x, xp), sign="<=")
# → concave curve → LP (no aux vars, domain bounds added)
# → non-concave  → SOS2/incremental, exact bound via output-link split

# tangent_lines remains public as the bare escape hatch for users who want
# chord expressions without the LP dispatch / checks / domain bounds.

Why the sign × curvature check raises instead of relaxing

Tangent lines impose the intersection of chord inequalities. Whether that equals the true hypograph/epigraph depends on curvature:

curvature	sign="<="	sign=">="
concave	hypograph (exact ✓)	y ≥ max_k chord_k(x) > f(x) — rejects y = f(x) ✗
convex	y ≤ min_k chord_k(x) < f(x) — rejects y = f(x) ✗	epigraph (exact ✓)
linear / mixed	tight / convex hull	tight / concave hull

In the ✗ cases, the feasible region is a strict subset of the true hypograph/epigraph — it silently rejects valid points. That's wrong, not merely loose, so we raise on explicit method="lp" and fall back to an exact MIP method for method="auto". Users who want the chord expressions anyway can still call linopy.tangent_lines() directly.

Notebook

examples/piecewise-inequality-bounds.ipynb:

• Mathematical formulations (equality, inequality, LP, incremental)
• 2-D feasibility plot
• 3-D feasibility ribbon for 3-variable inequality (three viewpoints)
• First-tuple convention + why we restrict to it
• Convexity × sign table with the wrong-region warning
• Worked fallback example for non-convex curves

Test plan

• 14 new tests: all sign values, all methods × curvature, LP vs SOS2 vs incremental solution consistency, error cases
• 135 piecewise tests pass
• Notebook runs end-to-end
• Review with @FabianHofmann @coroa

🤖 Generated with Claude Code

[FBumann]

🤖 Claude Code review

Overall this is a well-motivated feature with a clean public API and a genuinely good notebook. But there are two real correctness bugs that let the LP path silently produce wrong feasible regions, plus several lesser issues.

🔴 Critical — correctness bugs

1. _detect_convexity flips sign when x is strictly decreasing

linopy/piecewise.py:571-605 assumes strictly increasing x_points but only strict monotonicity is checked upstream. With decreasing x, dy/dx picks up an extra sign, so a concave graph is classified as convex (and vice-versa).

Repro:

m.add_piecewise_formulation(
    (y, [35, 30, 20, 0]),     # same concave graph as the PR example,
    (x, [30, 20, 10, 0]),     # just passed with x decreasing
    sign=">=",                  # docs call this case "wrong, not loose"
)
# auto picks method="lp"; solving gives y=30 at x=15
# correct SOS2 solution: y=25 = f(15)

Auto dispatches LP on a case the PR explicitly says must fall back, and the resulting constraints cut the curve. Explicit method=\"lp\" has the same hole — validation happily accepts it.

Fix options: (a) flip y by sign(dx) before classifying; (b) require dx > 0 strictly; (c) sort before classifying. Add a test with decreasing x.

2. _add_lp produces wrong constraints with NaN-padded per-entity breakpoints

_add_lp (piecewise.py:1278-1313) never consumes bp_mask. For multi-entity inputs where some entities have fewer breakpoints (NaN tail), NaN slopes/intercepts end up as 0 coefficients → constraint y ≤ 0 for the padded segment.

Repro:

bp_y = pd.DataFrame([[0,20,30,35],[0,10,20,nan]], index=['a','b'])
bp_x = pd.DataFrame([[0,10,20,30],[0,5,15,nan]],  index=['a','b'])
# LP:  y_b forced to 0 (f_b(10) = 15 expected)
# SOS2: y_b = 15 ✓

LP's siblings (_add_sos2, _add_incremental, _add_disjunctive) all respect bp_mask via lambda_mask / delta_mask. LP needs either (a) a constraint mask on the chord line, or (b) skip padded segments, or (c) reject per-entity masked input upfront.

🟠 Docstring / API contradictions

3. Wrong sign doc in the Parameters section

add_piecewise_formulation contradicts itself:

• Prose (piecewise.py:706-713): "The sign parameter follows the first-tuple convention… the first tuple's expression is bounded above…"
• Parameters (piecewise.py:735-737): "Constraint sign applied to the last tuple's link constraint."

Implementation matches the prose (first tuple = output). The Parameters block is stale.

4. Silently accepted sign aliases aren't in Literal[…]

sign = sign_replace_dict.get(sign, sign) (piecewise.py:757) accepts \"=\", \">\", \"<\", but the type hint is Literal[\"==\", \"<=\", \">=\"]. Either widen the Literal or drop the normalization.

🟠 Semantic gaps

5. active + sign != EQUAL is under-specified and untested

With active=0:

• sign=\"==\": output forced to 0 (desired "off" behavior). ✓
• sign=\"<=\": output is only forced to ≤ 0, not = 0. Likely surprising — a deactivated unit can still take any non-positive output.

No test combines active with a non-equality sign. Either document the asymmetry or add a both-sided constraint when combined. At minimum, a test.

6. Per-entity convexity classification is global

_detect_convexity aggregates across all entities. Safe for the concave+convex mix (\"mixed\" → fall back), but brittle if entities share a label while one has decreasing x (see bug 1). Either iterate per-entity or at least comment that convexity must be homogeneous.

🟡 Duplication / refactor

7. _add_lp duplicates tangent_lines

_add_lp (piecewise.py:1292-1304) and tangent_lines (piecewise.py:443-505) compute slopes / x_base / y_base / intercepts identically. Since tangent_lines is the documented escape hatch, _add_lp should call it:

def _add_lp(model, name, x_expr, y_expr, x_points, y_points, sign):
    tangents = tangent_lines(x_expr, x_points, y_points)
    _add_signed_link(model, y_expr, tangents, sign, f\"{name}{PWL_LP_SUFFIX}\")
    x_min = x_points.min(dim=BREAKPOINT_DIM)
    x_max = x_points.max(dim=BREAKPOINT_DIM)
    model.add_constraints(x_expr >= x_min, name=f\"{name}{PWL_DOMAIN_LO_SUFFIX}\")
    model.add_constraints(x_expr <= x_max, name=f\"{name}{PWL_DOMAIN_HI_SUFFIX}\")

Fixing mask handling once in tangent_lines benefits both call sites.

8. _add_sos2 / _add_incremental signatures ballooned

Both now take input_expr, output_expr, stacked_bp, input_stacked_bp, output_bp, bp_mask, link_dim, rhs, sign (+ active). Nine-ten params with Optional pairs that must stay consistent. Either package into a small _PwlLinks dataclass, or split so each function builds both links from bp_list internally.

9. _add_continuous has an EQUAL / non-EQUAL split that duplicates structure

Lines 993-1014 branch on sign == EQUAL to build either a single stacked link or a split output+input. Both paths then flow into the same _add_sos2 / _add_incremental which again branch internally on input_expr is None. The split-vs-single distinction crosses three function boundaries — unify or push up.

10. Auto-dispatch is a 5-deep nested if

piecewise.py:919-939. Flatten with early return and a helper _can_use_lp(lin_exprs, bp_list, sign, active) -> bool.

🟡 Housekeeping

11. Stale # type: ignore reintroduced

linopy/expressions.py:2359 re-adds a # type: ignore that the base branch removed in f0267f5 fix: remove unused type: ignore on merge() cls assignment. Will vanish on rebase.

12. Growing suffix-constant surface

_lp, _domain_lo, _domain_hi, _output_link — four new module-level constants for names not reused anywhere else. Consider inlining the f-strings or bundling.

13. PWL_LP_SUFFIX = \"_lp\" is generic

pwl0_lp doesn't identify which direction. Last chance to rename — _chord or _tangent (matches tangent_lines) would be clearer.

🟡 Test coverage gaps

• No test with decreasing x (would have caught bug 1).
• No test with per-entity masked breakpoints + LP (would have caught bug 2).
• No test solving a domain-bound violation (fix x beyond x_max, check infeasibility).
• No test with active + sign != EQUAL (issue 5).
• No test for LP on a linear curve (docs say it's accepted).
• No test for LP on N-D variables (extra dim beyond breakpoint dim).
• test_lp_consistency_with_sos2 checks only the max side — stronger if it also solved with a negative-f objective, testing both sides of the inequality.

✅ Things I liked

• _add_signed_link avoids if/elif ladders at every call site.
• "Raise, don't relax" stance on method=\"lp\" with mismatched curvature is correct and well-justified in both docstring and notebook.
• Reuse of SIGNS / sign_replace_dict from constants.py rather than re-parsing.
• First-tuple convention is the right call — per-tuple signs would create ambiguous cases.
• The notebook's treatment of why concave + >= is wrong (not loose) is the best explanation of chord-intersection semantics I've seen.

Suggested minimum before merging

1. Fix _detect_convexity order-dependence (bug 1) + test.
2. Fix LP NaN-mask handling (bug 2) — ideally by delegating _add_lp to tangent_lines (issue 7) + test.
3. Fix the sign Parameters docstring (issue 3).
4. Decide & document active + non-EQUAL sign semantics (issue 5).
5. Rebase to drop the stale expressions.py diff (issue 11).

Issues 8-10, 12-13 are refactor hygiene — nice to do now before the shape freezes, but not merge blockers.

---

🤖 Review generated with Claude Code (Opus 4.7)

[FBumann]

🤖 Claude Code review — round 2 (holistic, combined with #638)

Reviewing this PR as a unit with #638, against master. Stats: ~5.5k lines added / 2.3k removed across piecewise.py, tests, two notebooks, and the RST guide.

Big picture

The combined diff turns add_piecewise_formulation(...) into the single entry point for piecewise work — equality, inequality, N-variable linking, per-entity breakpoints, unit-commitment gating, and auto method dispatch all live behind one call. tangent_lines() survives as the bare "give me chord expressions" escape hatch. That is a coherent story, and the layered split between #638 (stateless construction + N-variable) and #663 (sign / LP) reads well as two landmarks.

Code-wise, after the last round of refactors the module is in good shape: _try_lp + _build_links + _resolve_sos2_vs_incremental + _PwlLinks cleanly separate concerns; _detect_convexity is now 8 lines and x-direction-invariant; _add_lp delegates to tangent_lines with a proper per-segment mask. The two correctness bugs from round 1 are genuinely gone.

🟠 The main gap: docs didn't keep up with the code

doc/piecewise-linear-constraints.rst is stale. It still reads like the #638-only state:

• method="auto" docstring (line 50) says "auto, sos2, or incremental" — "lp" is missing.
• The sign parameter is nowhere in the RST.
• The "When to Use What" table (lines 93-128) still positions add_piecewise_formulation as equality-only and tangent_lines as inequality-only. After this PR that's no longer true — add_piecewise_formulation(..., sign="<=") does inequality with domain bounds and curvature check baked in.
• Worst: line 127 says "For other combinations the bound is valid but loose (a relaxation)." — the opposite of the "wrong, not loose" conclusion the PR body and notebook argue at length. This one contradicts shipped code.
• "Generated variables and constraints" (lines 372-388) doesn't list the new suffixes (_chord, _domain_lo, _domain_hi, _output_link) and lists _x_link / _fill which aren't even correct for the current code (which uses _link / _fill_order). The latter is pre-existing staleness from feat: Unify piecewise API behind add_piecewise_formulation (sign + LP dispatch) #638, but it shows nobody has swept this page recently.
• No mention of PiecewiseFormulation.convexity on the result object.

This is the single biggest blocker for me. A user reading the docs after merge will be misled. Either update this RST page in this PR, or carve it out as an explicit follow-up ticket in the description.

Release notes — one bullet for all of #638 + all of #663 is hard to scan. The sign / LP feature deserves its own bullet ("Add one-sided bounds via sign parameter; pure-LP tangent dispatch for convex/concave cases with automatic domain bounds and convexity check"). The active parameter also vanished from the notes during the rewrite — it was in the pre-#638 entry, got folded into the add_piecewise_formulation mega-bullet, and a reader scanning for unit-commitment support won't find the keyword.

🟠 Notebook coordination

You now ship two piecewise notebooks:

• piecewise-linear-constraints.ipynb (expanded in feat: Unify piecewise API behind add_piecewise_formulation (sign + LP dispatch) #638) — 8 examples. Example 4 is "Tangent lines — Concave efficiency bound" and uses the bare tangent_lines() + add_constraints(fuel <= t).
• piecewise-inequality-bounds.ipynb (new in this PR) — covers exactly the same use case through the new sign="<=" path.

After this PR the bare-tangent_lines recipe is arguably the worse one (no domain bounds, no convexity check). Options:

1. Preferred: fold the new notebook into example 4 of the main notebook (or add an "Inequality bounds" section), delete the standalone. Users get one tutorial, not two.
2. Keep both but add a one-line "this shows the low-level building block; see the other notebook for the ergonomic API" in example 4, and a reciprocal link in the new notebook.

Right now a reader has no hint that two equivalent recipes exist.

🟡 netCDF round-trip drops convexity

before save: incremental  concave
after load:  incremental  None

linopy/io.py serialises method / variable_names / constraint_names but not convexity, so the attribute silently becomes None on reload. Either persist it or document that it's recomputed-only (it actually could just be recomputed on load from the saved constraints, though that's extra plumbing). Simpler: add it to the dict.

🟡 PR description nits

• "sign: Literal["==", "<=", ">="], default "="" — default should be "==" to match the Literal. Minor.
• Test plan line still says "135 piecewise tests", but you've added more since; actual count is higher.
• The description doesn't mention two shipped features: PiecewiseFormulation.convexity metadata + LP skip-reason logging. Both are user-visible — worth a line.
• Consider adding a short "tangent_lines vs sign="<=" — when to reach for which" note. Right now a first-time reader sees both and has no guidance.

🟡 Smaller code-level observations (all post-refactor)

1. _PwlLinks.eq_bp == stacked_bp when sign == EQUAL. Same object reference, just aliased. Fine, but either null out eq_bp in the EQUAL branch and have consumers fall back to stacked_bp, or add a one-line comment noting the intentional aliasing. Today it looks like two independent fields.

2. _add_lp's per-segment validity mask is computed inline (piecewise.py:1441-1449). SOS2 / incremental / disjunctive all get their masks via links.bp_mask from _build_links. A small helper _segment_validity_mask(bp) used by both _add_lp and _add_incremental's delta_mask construction would make the parallel obvious and keep the "mask everywhere" discipline centralised.

3. tangent_lines still says in its docstring "Breakpoint x-coordinates (must be strictly increasing)". After round 1 you made _detect_convexity direction-invariant, but tangent_lines itself doesn't enforce or rely on increasing order — yet the docstring says it's required. Either add validation or soften the docstring to "strictly monotonic".

4. _detect_convexity and tangent_lines both recompute slopes from scratch. Considering _add_lp now delegates to tangent_lines, and LP dispatch always calls _detect_convexity right before, the slope computation happens twice per call. Not a performance issue — but _detect_convexity(slopes) taking the pre-computed slopes DataArray (or exposing a _segment_slopes(x, y) shared helper) would be the cleanest shape. Low priority.

5. Is tangent_lines still the right public name? After this PR, it's the "advanced, no-checks, no-domain-bounds" escape hatch. The docstring still recommends it for the common case, which now has a better path. Consider a docstring note like "For most users, prefer add_piecewise_formulation(..., sign="<=") — this helper is the low-level building block."

✅ What I'd call shipped-quality

• The first-tuple convention — rejecting per-tuple signs prevents ambiguous cases and keeps the math clean. The notebook explains it well.
• _detect_convexity after the simplification + direction fix is a good example of "do one thing, 8 lines".
• _PwlLinks dataclass is the right call over the original 9-parameter signatures.
• LP skip-reason logging is a nice usability touch (user can tell why auto didn't pick LP).
• Error messages are now terse and actionable.

Minimum before merging

1. Update doc/piecewise-linear-constraints.rst to document sign, method="lp", the first-tuple convention, and fix the "loose" vs "wrong" statement. Also correct the stale constraint-suffix listings.
2. Split the release-notes bullet; make sure sign / LP / active / convexity are each findable.
3. Decide on the two-notebook story (merge, or cross-link).
4. Persist convexity in to_netcdf so round-trip is lossless.
5. Fix the default-"=" typo and refresh test count in the PR description.

Code changes I'd defer to a follow-up: the small _PwlLinks.eq_bp tidy, the _segment_validity_mask helper, the tangent_lines docstring tweak. None block merge.

---

🤖 Review generated with Claude Code (Opus 4.7)