Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Piecewise inequalities — the sign parameter#
add_piecewise_formulation accepts a sign parameter to express one-sided bounds of the form y ≤ f(x) or y ≥ f(x):
m.add_piecewise_formulation(
(fuel, y_pts), # output — gets the sign
(power, x_pts), # input — always equality
sign="<=",
)
This notebook walks through the math, the first-tuple convention, and the feasible regions produced by each method (LP, SOS2, incremental).
Key points#
Behaviour |
|
|---|---|
|
All expressions lie exactly on the curve. |
|
First expression is bounded above by |
|
First expression is bounded below by |
First-tuple convention: only the first tuple’s variable gets the sign. All remaining tuples are equality (inputs on the curve). This restriction keeps the semantics unambiguous — it’s always “output sign function(inputs)”.
[1]:
import matplotlib.pyplot as plt
import numpy as np
import linopy
Mathematical formulation#
Equality (sign="==")#
For \(N\) expressions \(e_1, \dots, e_N\) with breakpoints \(B_{j,0}, \dots, B_{j,n}\) per expression \(j\), the SOS2 formulation introduces interpolation weights \(\lambda_i \in [0,1]\):
Every expression is tied to the same \(\lambda\) — they share a single point on the curve.
Inequality (sign="<=" or ">=", first-tuple convention)#
The first expression \(e_1\) is the output; the rest are inputs forced to equality:
Inputs \(e_2, \dots, e_N\) are pinned to the curve at a shared \(\lambda\); the output \(e_1\) is then bounded (above or below) by the interpolated value. The internal split is visible in the generated constraints: a single stacked *_link for inputs and a separate *_output_link carrying the sign.
LP method (2-variable inequality, convex/concave curve)#
For \(y \le f(x)\) on a concave \(f\) (or \(y \ge f(x)\) on convex), we add one tangent (chord) per segment \(k\):
where \(m_k = (y_{k+1}-y_k)/(x_{k+1}-x_k)\) and \(c_k = y_k - m_k x_k\). The intersection of all chord inequalities equals the hypograph within the x-domain. No auxiliary variables are created.
Incremental (delta) formulation#
An MIP alternative to SOS2 for strictly monotonic breakpoints, using fill fractions \(\delta_i \in [0,1]\) and binaries \(z_i\) per segment:
Same sign split as SOS2: inputs use equality, output uses the requested sign.
Setup — a concave curve#
We use a concave, monotonically increasing curve. With sign="<=", the LP method is applicable (concave + <= is a tight relaxation).
[2]:
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 sign="<=" 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 (input equality, output signed). Solver picks the segment.
``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.
[3]:
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), # output, signed
(power, x_pts), # input, ==
sign="<=",
method=method,
)
m.add_constraints(power == power_val)
m.add_objective(-fuel) # maximise fuel to push against the bound
m.solve()
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-ozec44xe.lp
Reading time = 0.00 seconds
obj: 6 rows, 2 columns, 9 nonzeros
Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 24.04 LTS")
CPU model: AMD EPYC 9R14, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 2 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 6 rows, 2 columns and 9 nonzeros (Min)
Model fingerprint: 0x5f732d5f
Model has 1 linear objective coefficients
Coefficient statistics:
Matrix range [5e-01, 2e+00]
Objective range [1e+00, 1e+00]
Bounds range [3e+01, 4e+01]
RHS range [1e+01, 3e+01]
Presolve removed 6 rows and 2 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration Objective Primal Inf. Dual Inf. Time
0 -2.5000000e+01 0.000000e+00 0.000000e+00 0s
Solved in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective -2.500000000e+01
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-srt_1grd.lp
Reading time = 0.00 seconds
obj: 4 rows, 6 columns, 13 nonzeros
Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 24.04 LTS")
CPU model: AMD EPYC 9R14, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 2 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 4 rows, 6 columns and 13 nonzeros (Min)
Model fingerprint: 0x9f7a8a10
Model has 1 linear objective coefficients
Model has 1 SOS constraint
Variable types: 6 continuous, 0 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 4e+01]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 4e+01]
RHS range [1e+00, 2e+01]
Presolve removed 4 rows and 6 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 2 available processors)
Solution count 1: -25
No other solutions better than -25
Optimal solution found (tolerance 1.00e-04)
Best objective -2.500000476837e+01, best bound -2.500000476837e+01, gap 0.0000%
/tmp/ipykernel_2080/3765156795.py:5: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
m.add_piecewise_formulation(
/tmp/ipykernel_2080/3765156795.py:5: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
m.add_piecewise_formulation(
Dual values of MILP couldn't be parsed
/tmp/ipykernel_2080/3765156795.py:5: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
m.add_piecewise_formulation(
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-w8qf3_wm.lp
Reading time = 0.00 seconds
obj: 10 rows, 8 columns, 23 nonzeros
Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 24.04 LTS")
CPU model: AMD EPYC 9R14, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 2 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 10 rows, 8 columns and 23 nonzeros (Min)
Model fingerprint: 0xa97916e1
Model has 1 linear objective coefficients
Variable types: 5 continuous, 3 integer (3 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 4e+01]
RHS range [2e+01, 2e+01]
Presolve removed 10 rows and 8 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 2 available processors)
Solution count 1: -25
No other solutions better than -25
Optimal solution found (tolerance 1.00e-04)
Best objective -2.500000000000e+01, best bound -2.500000000000e+01, gap 0.0000%
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 an input-equality constraint plus a signed output link — but the feasible region is the same.
Visualising the feasible region#
The feasible region for (power, fuel) with sign="<=" is the hypograph of f restricted to the curve’s x-domain:
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)— abovef(15), should be infeasible ✗(35, 20)— power beyond domain, infeasible ✗
[4]:
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()
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
x_pts_3d = np.array([0.0, 30.0, 60.0, 100.0]) # power z_pts_3d = np.array([0.0, 25.0, 55.0, 95.0]) # heat y_pts_3d = np.array([0.0, 40.0, 85.0, 160.0]) # fuel (output)
Dense parameterisation of the 1-D curve#
t_grid = np.linspace(0, len(x_pts_3d) - 1, 80) power_c = np.interp(t_grid, np.arange(len(x_pts_3d)), x_pts_3d) heat_c = np.interp(t_grid, np.arange(len(z_pts_3d)), z_pts_3d) fuel_c = np.interp(t_grid, np.arange(len(y_pts_3d)), y_pts_3d)
fig = plt.figure(figsize=(7, 6)) ax = fig.add_subplot(111, projection=”3d”)
Shaded ribbon: (power(t), heat(t), fuel) for fuel from 0 to f(t)#
polys = [] for i in range(len(t_grid) - 1): quad = [ (power_c[i], heat_c[i], 0.0), (power_c[i + 1], heat_c[i + 1], 0.0), (power_c[i + 1], heat_c[i + 1], fuel_c[i + 1]), (power_c[i], heat_c[i], fuel_c[i]), ] polys.append(quad) coll = Poly3DCollection(polys, facecolor=”lightsteelblue”, edgecolor=”none”, alpha=0.35) ax.add_collection3d(coll)
Upper boundary: the curve itself#
ax.plot(power_c, heat_c, fuel_c, color=”C0”, lw=2.5, label=”curve f(t)”) ax.scatter(x_pts_3d, z_pts_3d, y_pts_3d, color=”C0”, s=50)
Lower boundary at fuel = 0#
ax.plot(power_c, heat_c, 0 * fuel_c, color=”gray”, lw=1, linestyle=”:”, alpha=0.7, label=”projection in (power, heat)”)
ax.set(xlabel=”power”, ylabel=”heat”, zlabel=”fuel”, title=”sign=’<=’ feasible region for 3 variables”) ax.view_init(elev=20, azim=-70) ax.legend() plt.tight_layout()
[5]:
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
x_pts_3d = np.array([0.0, 30.0, 60.0, 100.0]) # power
z_pts_3d = np.array([0.0, 25.0, 55.0, 95.0]) # heat
y_pts_3d = np.array([0.0, 40.0, 85.0, 160.0]) # fuel (output)
# Dense parameterisation of the 1-D curve
t_grid = np.linspace(0, len(x_pts_3d) - 1, 80)
power_c = np.interp(t_grid, np.arange(len(x_pts_3d)), x_pts_3d)
heat_c = np.interp(t_grid, np.arange(len(z_pts_3d)), z_pts_3d)
fuel_c = np.interp(t_grid, np.arange(len(y_pts_3d)), y_pts_3d)
def draw_ribbon(ax, elev, azim, title):
# Shaded ribbon: (power(t), heat(t), fuel) for fuel from 0 to f(t)
polys = []
for i in range(len(t_grid) - 1):
quad = [
(power_c[i], heat_c[i], 0.0),
(power_c[i + 1], heat_c[i + 1], 0.0),
(power_c[i + 1], heat_c[i + 1], fuel_c[i + 1]),
(power_c[i], heat_c[i], fuel_c[i]),
]
polys.append(quad)
coll = Poly3DCollection(
polys, facecolor="lightsteelblue", edgecolor="none", alpha=0.35
)
ax.add_collection3d(coll)
# Upper boundary: the curve itself
ax.plot(power_c, heat_c, fuel_c, color="C0", lw=2.5)
ax.scatter(x_pts_3d, z_pts_3d, y_pts_3d, color="C0", s=40)
# (power, heat) projection at fuel=0
ax.plot(power_c, heat_c, 0 * fuel_c, color="gray", lw=1, linestyle=":", alpha=0.7)
ax.set(xlabel="power", ylabel="heat", zlabel="fuel", title=title)
ax.view_init(elev=elev, azim=azim)
fig = plt.figure(figsize=(13, 5))
ax1 = fig.add_subplot(131, projection="3d")
draw_ribbon(ax1, elev=20, azim=-70, title="perspective")
ax2 = fig.add_subplot(132, projection="3d")
draw_ribbon(ax2, elev=0, azim=-90, title="side (power vs fuel)")
ax3 = fig.add_subplot(133, projection="3d")
draw_ribbon(ax3, elev=90, azim=-90, title="top (power vs heat)")
plt.tight_layout()
The ribbon has the shape (input curve) × [0, output value]. Points above the curve in fuel are infeasible; points below are feasible, provided (power, heat) lies on the curve’s projection. If the user tried to set power=50, heat=20 (off-curve), the formulation would be infeasible — the inputs must be consistent with a shared \(\lambda\).
The first-tuple convention#
Why does only one variable get the sign? Because the math of “bound one output by a function of the others” has a single inequality direction. For 3+ variables:
m.add_piecewise_formulation(
(fuel, y_pts), # output — sign
(power, x_pts), # input — ==
(heat, z_pts), # input — ==
sign="<=",
)
reads as fuel ≤ g(power, heat) on the joint curve. All inputs must lie on the curve (equality); only the output is bounded.
Allowing arbitrary per-variable signs would open up cases like “fuel ≤ f(power) AND heat ≤ f(power)” which is a dominated region, not a hypograph — mathematically valid but rarely what users want. Restricting to one output keeps the semantics unambiguous.
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 |
|
|
|---|---|---|
concave |
hypograph (exact ✓) |
wrong region — requires |
convex |
wrong region — requires |
epigraph (exact ✓) |
linear |
exact |
exact |
mixed (non-convex) |
convex hull of |
concave hull of |
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 segment extrapolated over another segment’s x-range lies above f, so the constraint 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 segment 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 sign="<=" does internally: it falls back to SOS2/incremental with the sign on the output link. No relaxation, no wrong bounds.
For 3+ variables with inequality, LP is never applicable: tangent lines express one input → one output. With multiple inputs on a 1-D curve in N-D space, identifying which segment we’re on requires SOS2/binary. Auto dispatches to SOS2 here.
[6]:
# 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), sign="<=")
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), sign=">=")
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), sign=">=", method="lp")
except ValueError as e:
print(f"lp(concave, '>=') → raises: {e}")
non-convex + '<=' → sos2
concave + '>=' → incremental
lp(concave, '>=') → raises: method='lp' with sign='>=' requires convex or linear curvature; got 'concave'. Use method='auto'.
/tmp/ipykernel_2080/1973565292.py:8: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
f1 = m1.add_piecewise_formulation((y1, y_nc), (x1, x_nc), sign="<=")
/tmp/ipykernel_2080/1973565292.py:18: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
f2 = m2.add_piecewise_formulation((y2, y_cc), (x2, x_cc), sign=">=")
/tmp/ipykernel_2080/1973565292.py:26: EvolvingAPIWarning: piecewise: add_piecewise_formulation is a new API; some details (e.g. the sign/first-tuple 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. Silence with `warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.
m3.add_piecewise_formulation((y3, y_cc), (x3, x_cc), sign=">=", method="lp")
Summary#
Use
sign="="(default) for exact equality on the curve.Use
sign="<="/sign=">="for one-sided bounds on the first tuple’s expression.method="auto"picks the most efficient formulation: LP for convex/concave 2-variable inequalities, otherwise SOS2 or incremental.Only the first tuple gets the sign — all other tuples are always equality. This restriction keeps semantics unambiguous.