Note

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

Creating Piecewise Inequality Bounds

When you only need a one-sided bound from a piecewise curve — y ≤ f(x) for a concave upper envelope, y ≥ f(x) for a convex lower envelope — add_piecewise_formulation accepts an optional sign as the third tuple element:

m.add_piecewise_formulation(
    (fuel,  fuel_pts, ">="),   # fuel ≥ f(power) — over-fuelling admissible
    (power, power_pts),        # equality role
)

The pay-off is a pure-LP encoding when the curve’s curvature matches the sign — no SOS2, no binaries. This notebook covers the geometry of the feasible region, the curvature × sign combinations that unlock the LP path, and what happens when they don’t match.

For the formulation math see the reference page; for the all-equality variant and other features see Creating Piecewise Linear Constraints.

Tuple roles

Tuple form	Role	What it constrains
(expr, breaks)	== (equality)	With 2+ equality tuples sharing weights, the joint point lies on the curve. With 1 equality + 1 bounded, the equality tuple’s marginal feasible set is just its breakpoint domain [x_min, x_max] — one coordinate alone can’t locate a curve point.
(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 warnings

import matplotlib.pyplot as plt
import numpy as np

import linopy

# Silence the evolving-API warning for cleaner tutorial output.
warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)

Setup — a convex heat-rate curve

A convex, monotonically increasing curve maps power output to the fuel required (the classic heat-rate curve). Bounding fuel by this curve with >= says the unit must consume at least the design fuel for a given power output — over-fuelling is physically admissible but wasteful, so an objective that minimises fuel pulls the operating point onto the curve. Convex + >= is exactly the combination that lets the LP method apply.

The breakpoint arrays:

power_pts = np.array([0.0, 30.0, 60.0, 100.0])
fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])  # slopes 1.2, 1.6, 2.15 (convex)

Three methods, identical feasible region

With fuel bounded >= and our convex curve, the three methods give the same feasible region for power ∈ [0, 100]:

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

• ``method=”sos2”`` — lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. 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=60, where f(60)=84.

def solve(method, power_val):
    m = linopy.Model()
    power = m.add_variables(lower=0, upper=100, name="power")
    fuel = m.add_variables(lower=0, upper=200, name="fuel")
    m.add_piecewise_formulation(
        (fuel, fuel_pts, ">="),  # fuel ≥ f(power) — over-fuelling admissible
        (power, power_pts),  # equality role (domain-bounded to [0, 100])
        method=method,
    )
    m.add_constraints(power == power_val)
    m.add_objective(fuel)  # minimise fuel against the lower bound
    m.solve(solver_name="highs", reformulate_sos="auto", 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, 60)
    print(f"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}")

lp          : fuel=84.00  vars=['power', 'fuel']  cons=['pwl0_chord', 'pwl0_domain_lo', 'pwl0_domain_hi', 'con0']
sos2        : fuel=84.00  vars=['power', 'fuel', 'pwl0_lambda']  cons=['pwl0_convex', 'pwl0_link', 'pwl0_output_link', 'con0']
incremental : fuel=84.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=84 at power=60 (which is f(60) 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 an equality constraint for the equality-signed tuple plus a signed link for the bounded tuple — but the feasible region is the same.

Visualising the feasible region

The feasible region for (power, fuel) with fuel bounded >= is the epigraph of f restricted to the power domain:

\[\{ (\mathrm{power}, \mathrm{fuel}) : 0 \le \mathrm{power} \le 100,\ \mathrm{fuel} \ge f(\mathrm{power}) \}\]

Below we colour feasible points green, infeasible ones red:

• (60, 100) — above the curve, 100 ≥ f(60)=84 ✓

• (60, 84) — on the curve ✓

• (60, 70) — below f(60), infeasible ✗

• (120, 100) — power beyond domain, infeasible ✗

def in_epigraph(px, fy):
    if px < power_pts[0] or px > power_pts[-1]:
        return False
    return fy >= np.interp(px, power_pts, fuel_pts)


xx, yy = np.meshgrid(np.linspace(-10, 130, 200), np.linspace(-10, 200, 200))
region = np.vectorize(in_epigraph)(xx, yy)

test_points = [(60, 100), (60, 84), (60, 70), (120, 100)]

fig, ax = plt.subplots(figsize=(6, 5))
ax.contourf(xx, yy, region, levels=[0.5, 1], colors=["lightsteelblue"], alpha=0.5)
ax.plot(power_pts, fuel_pts, "o-", color="C0", lw=2, label="f(power)")
for px, fy in test_points:
    feas = in_epigraph(px, fy)
    ax.scatter(
        [px], [fy], color="green" if feas else "red", zorder=5, s=80, edgecolors="black"
    )
    ax.annotate(f"({px}, {fy})", (px, fy), textcoords="offset points", xytext=(8, 5))
ax.set(
    xlabel="power",
    ylabel="fuel",
    title="sign='>=' feasible region — epigraph of f(power) on [0, 100]",
)
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)
power_nc = [0, 30, 60, 100]
fuel_nc = [0, 50, 30, 80]  # slopes change sign → mixed convexity

m1 = linopy.Model()
power1 = m1.add_variables(lower=0, upper=100, name="power")
fuel1 = m1.add_variables(lower=0, upper=200, name="fuel")
f1 = m1.add_piecewise_formulation((fuel1, fuel_nc, ">="), (power1, power_nc))
print(f"non-convex + '>=' → {f1.method}")

# 2. Convex curve + sign='<=': LP would be loose → auto falls back to MIP
m2 = linopy.Model()
power2 = m2.add_variables(lower=0, upper=100, name="power")
fuel2 = m2.add_variables(lower=0, upper=200, name="fuel")
f2 = m2.add_piecewise_formulation(
    (fuel2, list(fuel_pts), "<="), (power2, list(power_pts))
)
print(f"convex + '<='     → {f2.method}")

# 3. Explicit method="lp" with mismatched curvature raises
try:
    m3 = linopy.Model()
    power3 = m3.add_variables(lower=0, upper=100, name="power")
    fuel3 = m3.add_variables(lower=0, upper=200, name="fuel")
    m3.add_piecewise_formulation(
        (fuel3, list(fuel_pts), "<="), (power3, list(power_pts)), method="lp"
    )
except ValueError as e:
    print(f"lp(convex, '<=')  → raises: {e}")

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

Summary

• One bounded tuple + a 2-variable formulation gives a hypograph (<=) or epigraph (>=) feasible region.

• Curvature × sign matching — concave + <= or convex + >= — lets method="auto" skip MIP entirely. Mismatched combinations fall back to SOS2/incremental with a signed link.

• ``method=”lp”`` is strict — it raises on a mismatched curvature rather than silently encoding the wrong region.

• At most one tuple may carry a non-== sign, and 3+ tuples must all be ==. Multi-bounded / N≥3 inequalities — open an issue at PyPSA/linopy#issues.

See also: reference page for the formulation math, Creating Piecewise Linear Constraints for all-equality, unit commitment, CHP, fleets, slopes.