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

sign="==" (default)

All expressions lie exactly on the curve.

sign="<="

First expression is bounded above by f(rest).

sign=">="

First expression is bounded below by f(rest).

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]\):

\[\sum_{i=0}^{n} \lambda_i = 1, \qquad \text{SOS2}(\lambda_0, \dots, \lambda_n), \qquad e_j = \sum_{i=0}^{n} \lambda_i \, B_{j,i} \ \ \forall j.\]

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:

\[\sum_{i=0}^{n} \lambda_i = 1, \qquad \text{SOS2}(\lambda), \qquad e_j = \sum_{i=0}^{n} \lambda_i \, B_{j,i}\ \ \forall j \ge 2, \qquad e_1 \ \text{sign}\ \sum_{i=0}^{n} \lambda_i \, B_{1,i}.\]

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\):

\[y \le m_k \cdot x + c_k \ \ \forall k, \qquad x_0 \le x \le x_n,\]

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:

\[\delta_{i+1} \le \delta_i, \quad z_{i+1} \le \delta_i, \quad \delta_i \le z_i, \qquad e_j = B_{j,0} + \sum_i \delta_i\,(B_{j,i+1}-B_{j,i}).\]

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()
_images/piecewise-inequality-bounds-tutorial_4_0.svg

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:

\[\{ (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 ✗

[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()
_images/piecewise-inequality-bounds-tutorial_9_0.svg

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()
_images/piecewise-inequality-bounds-tutorial_11_0.svg

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

sign="<="

sign=">="

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 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.