Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Piecewise inequalities — per-tuple sign

add_piecewise_formulation accepts an optional third tuple element, "<=" or ">=", that marks one expression as bounded by the piecewise curve instead of pinned to it:

m.add_piecewise_formulation(
    (fuel,  y_pts, "<="),   # bounded above by the curve
    (power, x_pts),         # pinned to the curve
)

This notebook walks through the geometry, the curvature × sign matching that lets method="auto" skip MIP machinery entirely, and the feasible regions produced by each method (LP, SOS2, incremental). For the formulation math see the reference page.

Key points

Tuple form	Behaviour
(expr, breaks)	Pinned: expr lies exactly on the curve.
(expr, breaks, "<=")	Bounded above: expr ≤ f(other tuples).
(expr, breaks, ">=")	Bounded below: expr ≥ f(other tuples).

Currently at most one tuple may carry a non-equality sign, and 3+ tuples must all be equality. Multi-bounded and N≥3 inequality cases aren’t supported yet — if you have a concrete use case, please open an issue at PyPSA/linopy#issues so we can scope it properly.

import matplotlib.pyplot as plt
import numpy as np

import linopy

Setup — a concave curve

We use a concave, monotonically increasing curve. With a tuple bounded <=, the LP method is applicable (concave + <= is a tight relaxation).

x_pts = np.array([0.0, 10.0, 20.0, 30.0])
y_pts = np.array([0.0, 20.0, 30.0, 35.0])  # slopes 2, 1, 0.5 (concave)

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(x_pts, y_pts, "o-", color="C0", lw=2)
ax.set(xlabel="power", ylabel="fuel", title="Concave reference curve f(x)")
ax.grid(alpha=0.3)
plt.tight_layout()

Three methods, identical feasible region

With one tuple bounded <= and our concave curve, the three methods give the same feasible region within [x_0, x_n]:

• ``method=”lp”`` — tangent lines + domain bounds. No auxiliary variables.

• ``method=”sos2”`` — lambdas + SOS2 + split link (pinned equality, bounded signed). Solver picks the piece.

• ``method=”incremental”`` — delta fractions + binaries + split link. Same mathematics, MIP encoding instead of SOS2.

method="auto" dispatches to "lp" whenever applicable — it’s always preferable because it’s pure LP.

Let’s verify they produce the same solution at power=15.

def solve(method, power_val):
    m = linopy.Model()
    power = m.add_variables(lower=0, upper=30, name="power")
    fuel = m.add_variables(lower=0, upper=40, name="fuel")
    m.add_piecewise_formulation(
        (fuel, y_pts, "<="),  # bounded
        (power, x_pts),  # pinned
        method=method,
    )
    m.add_constraints(power == power_val)
    m.add_objective(-fuel)  # maximise fuel to push against the bound
    m.solve(output_flag=False)
    return float(m.solution["fuel"]), list(m.variables), list(m.constraints)


for method in ["lp", "sos2", "incremental"]:
    fuel_val, vars_, cons_ = solve(method, 15)
    print(f"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}")

Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-ve6z9phj.lp
Reading time = 0.00 seconds
obj: 6 rows, 2 columns, 9 nonzeros

/tmp/ipykernel_2119/3538038423.py:5: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the per-tuple sign convention, active+sign semantics) may be refined in minor releases.  Please share your use cases or concerns at https://github.com/PyPSA/linopy/issues — your feedback shapes what stabilises.  This warning fires once per session; silence entirely with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
  m.add_piecewise_formulation(

lp          : fuel=25.00  vars=['power', 'fuel']  cons=['pwl0_chord', 'pwl0_domain_lo', 'pwl0_domain_hi', 'con0']
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-i_r1ijfh.lp
Reading time = 0.00 seconds
obj: 4 rows, 6 columns, 13 nonzeros

Dual values of MILP couldn't be parsed

sos2        : fuel=25.00  vars=['power', 'fuel', 'pwl0_lambda']  cons=['pwl0_convex', 'pwl0_link', 'pwl0_output_link', 'con0']
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-_7pdyslf.lp
Reading time = 0.00 seconds
obj: 10 rows, 8 columns, 23 nonzeros

Dual values of MILP couldn't be parsed

incremental : fuel=25.00  vars=['power', 'fuel', 'pwl0_delta', 'pwl0_order_binary']  cons=['pwl0_delta_bound', 'pwl0_fill_order', 'pwl0_binary_order', 'pwl0_link', 'pwl0_output_link', 'con0']

All three give fuel=25 at power=15 (which is f(15) exactly) — the math is equivalent. The LP method is strictly cheaper: no auxiliary variables, just three chord constraints and two domain bounds.

The SOS2 and incremental methods create lambdas (or deltas + binaries) and split the link into a pinned-equality constraint plus a signed bounded link — but the feasible region is the same.

Visualising the feasible region

The feasible region for (power, fuel) with fuel bounded <= is the hypograph of f restricted to the curve’s x-domain:

\[\{ (x, y) : x_0 \le x \le x_n,\ y \le f(x) \}\]

We colour green feasible points, red infeasible ones. Three test points:

• (15, 15) — inside the curve, 15 ≤ f(15)=25 ✓

• (15, 25) — on the curve ✓

• (15, 29) — above f(15), should be infeasible ✗

• (35, 20) — power beyond domain, infeasible ✗

def in_hypograph(px, py):
    if px < x_pts[0] or px > x_pts[-1]:
        return False
    return py <= np.interp(px, x_pts, y_pts)


xx, yy = np.meshgrid(np.linspace(-2, 38, 200), np.linspace(-5, 45, 200))
region = np.vectorize(in_hypograph)(xx, yy)

test_points = [(15, 15), (15, 25), (15, 29), (35, 20)]

fig, ax = plt.subplots(figsize=(6, 5))
ax.contourf(xx, yy, region, levels=[0.5, 1], colors=["lightsteelblue"], alpha=0.5)
ax.plot(x_pts, y_pts, "o-", color="C0", lw=2, label="f(x)")
for px, py in test_points:
    feas = in_hypograph(px, py)
    ax.scatter(
        [px], [py], color="green" if feas else "red", zorder=5, s=80, edgecolors="black"
    )
    ax.annotate(f"({px}, {py})", (px, py), textcoords="offset points", xytext=(8, 5))
ax.set(
    xlabel="power",
    ylabel="fuel",
    title="sign='<=' feasible region — hypograph of f(x) on [x_0, x_n]",
)
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()

When is LP the right choice?

tangent_lines imposes the intersection of chord inequalities. Whether that intersection matches the true hypograph/epigraph of f depends on the curvature × sign combination:

curvature	bounded <=	bounded >=
concave	hypograph (exact ✓)	wrong region — requires y ≥ max_k chord_k(x) > f(x)
convex	wrong region — requires y ≤ min_k chord_k(x) < f(x)	epigraph (exact ✓)
linear	exact	exact
mixed (non-convex)	convex hull of f (wrong for exact hypograph)	concave hull of f (wrong for exact epigraph)

In the ✗ cases, tangent lines do not give a loose relaxation — they give a strictly wrong feasible region that rejects points satisfying the true constraint. Example: for a concave f with y ≥ f(x), the chord of any piece extrapolated over another piece’s x-range lies above f, so y ≥ max_k chord_k(x) forbids y = f(x) itself.

method="auto" dispatches to LP only in the two exact cases (concave + <= or convex + >=). For the other combinations it falls back to SOS2 or incremental, which encode the hypograph/epigraph exactly via discrete piece selection.

method="lp" explicitly forces LP and raises on a mismatched curvature rather than silently producing a wrong feasible region.

For non-convex curves with either sign, the only exact option is a piecewise formulation. That’s what the bounded-tuple path does internally: it falls back to SOS2/incremental with the sign on the bounded link. No relaxation, no wrong bounds.

# 1. Non-convex curve: auto falls back (LP relaxation would be loose)
x_nc = [0, 10, 20, 30]
y_nc = [0, 20, 10, 30]  # slopes change sign → mixed convexity

m1 = linopy.Model()
x1 = m1.add_variables(lower=0, upper=30, name="x")
y1 = m1.add_variables(lower=0, upper=40, name="y")
f1 = m1.add_piecewise_formulation((y1, y_nc, "<="), (x1, x_nc))
print(f"non-convex + '<=' → {f1.method}")

# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP
x_cc = [0, 10, 20, 30]
y_cc = [0, 20, 30, 35]  # concave

m2 = linopy.Model()
x2 = m2.add_variables(lower=0, upper=30, name="x")
y2 = m2.add_variables(lower=0, upper=40, name="y")
f2 = m2.add_piecewise_formulation((y2, y_cc, ">="), (x2, x_cc))
print(f"concave + '>='    → {f2.method}")

# 3. Explicit method="lp" with mismatched curvature raises
try:
    m3 = linopy.Model()
    x3 = m3.add_variables(lower=0, upper=30, name="x")
    y3 = m3.add_variables(lower=0, upper=40, name="y")
    m3.add_piecewise_formulation((y3, y_cc, ">="), (x3, x_cc), method="lp")
except ValueError as e:
    print(f"lp(concave, '>=')  → raises: {e}")

non-convex + '<=' → sos2
concave + '>='    → incremental
lp(concave, '>=')  → raises: method='lp' is not applicable: sign='>=' needs convex/linear curvature, got 'concave'. Use method='auto'.

Summary

• Default is all-equality: every tuple lies on the curve.

• Append "<=" or ">=" as a third tuple element to mark one expression as bounded by the curve.

• method="auto" picks the most efficient formulation: LP for matching-curvature 2-variable inequalities, otherwise SOS2 or incremental.

• At most one tuple may be bounded; with 3+ tuples all must be equality. Multi-bounded and N≥3 inequality use cases — please open an issue at PyPSA/linopy#issues so we can scope them.