feat: add sign parameter to add_piecewise_formulation

[FBumann]

Adds "==" (default), "<=", and ">=" as a uniform modifier on the link constraint. Works for N-variable formulations — sign applies symmetrically to all expressions, bounding the operating point above/below the curve.

For 2-var convex/concave cases, users can still use tangent_lines() for pure LP. This adds the inequality option for non-convex cases and for N-variable feasible-region modelling.

Closes # (if applicable).

Changes proposed in this Pull Request

Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• Unit tests for new features were added (if applicable).
• A note for the release notes doc/release_notes.rst of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.

[FBumann]

🤖 Snippet drafted by Claude Code

Clarification on tangent_lines domain behaviour (surfaced during notebook discussion):

Looking at the current tangent_lines implementation:

def tangent_lines(x, x_points, y_points):
    # ... computes slopes, intercepts
    return slopes * _to_linexpr(x) + intercepts

It returns a LinearExpression (one chord per segment) — no domain constraints are added, and the function has no way to add them because it doesn't take a model. When the user writes m.add_constraints(fuel <= t), they get one linear inequality per segment, e.g. fuel <= 0.5·power + 20, extending indefinitely in both directions.

So in the visualisation, the feasible region continues past x_0 and x_n along the outer chord extrapolations. The curve's "endpoints" aren't enforced unless the user explicitly adds variable bounds.

In the pre-refactor code, _add_pwl_lp used to include explicit PWL_LP_DOMAIN_SUFFIX constraints (x >= x_min, x <= x_max). That safety net is gone in the new tangent_lines — it's intentionally a "bare math" utility.

Why this matters for the sign discussion

None of the add_piecewise_formulation(sign=...) options recover this behaviour. SOS2/incremental lambdas are in [0, 1] with sum(lambda) = 1, so the interpolated point is always a convex combination of breakpoints — it can never leave [x_0, x_n] × [y_min, y_max]. So:

At x > x_n	tangent_lines	(==, <=)	(>=, <=)	(<=, <=)
Behaviour	fuel ≤ extrapolated chord	infeasible	fuel ≤ y_n	infeasible

tangent_lines is fundamentally about trends extending past the breakpoints; the piecewise formulations are about hard operating regions bounded by breakpoints.

This nudges toward folding tangent_lines into add_piecewise_formulation

Using tangent_lines cleanly actually requires the user to:

1. Call linopy.tangent_lines(x, x_points, y_points) to get the expression
2. Add m.add_constraints(y <= t) — one constraint per segment
3. Add m.add_constraints(x >= x_min) and m.add_constraints(x <= x_max) to prevent chord extrapolation

That's three add_* calls and multiple constraints. In other words, tangent_lines + downstream constraints is itself a formulation — just a pure-LP one. The name "tangent_lines" (returning an expression) understates what the user ends up adding to the model, and leaves domain handling to the user.

If we fold this into add_piecewise_formulation, the user gets:

• A single entry point with consistent semantics across all methods
• Automatic domain bounds (closing the footgun)
• A PiecewiseFormulation return object for inspection, tracking, and unified naming
• Internal dispatch: if the function is convex/concave and the user asks for an inequality, use the pure-LP tangent-lines approach (no aux variables). Otherwise fall back to SOS2/incremental with <= on the link.

Options

1. Leave as-is, document the footgun — current state
2. Fold tangent_lines into add_piecewise_formulation as an internal LP method. Inequality signs (sign="<="/">=") dispatch to the pure-LP formulation when the function is convex/concave, else to SOS2/incremental. tangent_lines either becomes internal or is kept public as the "bare math" escape hatch.

Option (2) feels like the most coherent API — the user's modelling intent ("y bounded by f(x)") always goes through one function, and the implementation picks the most efficient formulation. Happy to prototype if that's the direction.

[FBumann]

A exploratory artifact, superseeded by #663