Piecewise Linear Constraints

Piecewise Linear Constraints#

Piecewise linear (PWL) constraints approximate nonlinear functions as connected linear segments, allowing you to model cost curves, efficiency curves, or production functions within a linear programming framework.

Quick Start#

Equality — lock variables onto the piecewise curve:

import linopy

m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=100)
fuel = m.add_variables(name="fuel")

# fuel = f(power) on the piecewise curve defined by these breakpoints
m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel, [0, 36, 84, 170]),
)

Inequality — bound one expression by the curve:

# fuel <= f(power).  "auto" picks the cheapest correct formulation
# (pure LP with chord constraints when the curve's curvature matches
# the requested sign; SOS2/incremental otherwise).
m.add_piecewise_formulation(
    (fuel, [0, 20, 30, 35]),  # bounded output listed FIRST
    (power, [0, 10, 20, 30]),  # input always on the curve
    sign="<=",
)

Each (expression, breakpoints) tuple pairs a variable with its breakpoint values. All tuples share interpolation weights, so at any feasible point every variable corresponds to the same point on the piecewise curve.

API#

add_piecewise_formulation#

m.add_piecewise_formulation(
    (expr1, breakpoints1),
    (expr2, breakpoints2),
    ...,
    sign="==",  # "==", "<=", or ">="
    method="auto",  # "auto", "sos2", "incremental", or "lp"
    active=None,  # binary variable to gate the constraint
    name=None,  # base name for generated variables/constraints
)

Creates auxiliary variables and constraints that enforce either an equality (sign="==", default) or a one-sided bound (sign="<=" / ">=") of the first expression by the piecewise function of the rest.

breakpoints and segments#

Factory functions that create DataArrays with the correct dimension names:

linopy.breakpoints([0, 50, 100])  # list
linopy.breakpoints({"gen1": [0, 50], "gen2": [0, 80]}, dim="gen")  # per-entity
linopy.breakpoints(slopes=[1.2, 1.4], x_points=[0, 30, 60], y0=0)  # from slopes
linopy.segments([(0, 10), (50, 100)])  # disjunctive
linopy.segments({"gen1": [(0, 10)], "gen2": [(0, 80)]}, dim="gen")

The sign parameter — equality vs inequality#

The sign argument of add_piecewise_formulation chooses whether all expressions are locked onto the curve or whether the first one is bounded:

  • sign="==" (default): every expression lies exactly on the piecewise curve — joint equality. All tuples are symmetric. The feasible region is a 1-D curve in N-space.

  • sign="<=": N−1 tuples are pinned to the curve (moving together along it), and the first tuple’s expression is bounded above by its interpolated value at that shared curve position — it can undershoot.

  • sign=">=": same split, but the first is bounded below (overshoots admitted).

Inequality relaxes one tuple’s curve-equality into a one-sided bound. The others keep moving together along the curve in lockstep — this is the first-tuple convention.

What this means geometrically.

For 2 variables ((y, yp), (x, xp), sign="<="), “N−1 pinned, 1 bounded” reduces to the familiar hypograph: x moves along the breakpoint axis, y ranges from its lower bound up to f(x).

For 3+ variables, the N−1 pinned tuples are jointly constrained to a single segment position on the piecewise curve. In a CHP characteristic (power, fuel, heat) with sign="<=", fuel and heat trace the curve simultaneously: specifying fuel = 85 determines the segment position, which in turn fixes heat to its curve value at that position. Assigning heat to any value inconsistent with the same segment renders the system infeasible.

The feasible region in N-space is a 2-dimensional manifold: the 1-D parametric curve at its upper boundary (for "<="), extended along the first tuple’s axis down to that variable’s lower bound.

Choice of bounded tuple. The first tuple should correspond to a quantity with a mechanism for below-curve operation — typically a controllable dissipation path: heat rejection via cooling tower (also called thermal curtailment), electrical curtailment, or emissions after post-treatment. Placing a consumption-side variable such as fuel intake in the bounded position yields a valid but loose formulation: the characteristic curve fixes fuel draw at a given load, so sign="<=" on fuel admits operating points the plant cannot physically realise. An objective that rewards lower fuel may then find a non-physical optimum — safe only when no such objective pressure exists.

Relatedly, inequality formulations can also be faster to solve: with 2 variables and matching curvature, method="auto" dispatches to the pure-LP chord formulation (no SOS2, no binaries). For N≥3 the solver still reaches for SOS2/incremental, but the relaxed feasible region often tightens the LP relaxation and reduces branch-and-bound work. Choose sign="==" when you want strict curve adherence (the tightest feasible region) and sign="<=" / ">=" when either the physics admits dissipation or the speedup is worth the relaxation.

When is a one-sided bound wanted?

For continuous curves, the main reason to reach for sign="<=" / ">=" is to unlock the LP chord formulation — no SOS2, no binaries, just pure LP. On a convex/concave curve with a matching sign, the chord inequalities are as tight as SOS2, so you get the same optimum with a cheaper model.

For disjunctive curves (segments(...)), sign is a first-class tool in its own right: disconnected operating regions with a bounded output, always exact regardless of segment curvature (see the disjunctive section below).

Beyond that: fuel-on-efficiency-envelope modelling (extra burn above the curve is admissible, cost is still bounded), emissions caps where the curve is itself a convex overestimator, or any situation where the curve bounds a variable that need not sit on it.

If the curvature doesn’t match the sign (convex + "<=", or concave + ">="), LP is not applicable — method="auto" falls back to SOS2/incremental with the signed output link, which gives a valid but much more expensive model. In that case prefer sign="==" unless you genuinely need the one-sided semantics; the equality formulation is typically simpler to reason about and no more expensive than the SOS2 inequality variant.

Math (2-variable ``sign=”<=”``, concave :math:`f`). The feasible region is the hypograph of \(f\) restricted to the breakpoint range:

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

For convex \(f\) with sign=">=", the feasible region is the epigraph. Mismatched sign+curvature (convex + <=, or concave + >=) describes a non-convex region — method="auto" will fall back to SOS2/incremental and method="lp" will raise. See the Piecewise inequalities — the sign parameter notebook for a full walkthrough.

Warning

With sign="<=" and active=0, the output is only bounded above by 0 — the lower side still comes from the output variable’s own lower bound. In the common case of non-negative outputs (fuel, cost, heat), set lower=0 on that variable: combined with the y ≤ 0 constraint from deactivation, this forces y = 0 automatically. See the docstring for the full recipe.

Breakpoint Construction#

From lists#

The simplest form — pass Python lists directly in the tuple:

m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel, [0, 36, 84, 170]),
)

With the breakpoints() factory#

Equivalent, but explicit about the DataArray construction:

m.add_piecewise_formulation(
    (power, linopy.breakpoints([0, 30, 60, 100])),
    (fuel, linopy.breakpoints([0, 36, 84, 170])),
)

From slopes#

When you know marginal costs (slopes) rather than absolute values:

m.add_piecewise_formulation(
    (power, [0, 50, 100, 150]),
    (
        cost,
        linopy.breakpoints(
            slopes=[1.1, 1.5, 1.9], x_points=[0, 50, 100, 150], y0=0
        ),
    ),
)
# cost breakpoints: [0, 55, 130, 225]

Per-entity breakpoints#

Different generators can have different curves. Pass a dict to breakpoints() with entity names as keys:

m.add_piecewise_formulation(
    (
        power,
        linopy.breakpoints(
            {"gas": [0, 30, 60, 100], "coal": [0, 50, 100, 150]}, dim="gen"
        ),
    ),
    (
        fuel,
        linopy.breakpoints(
            {"gas": [0, 40, 90, 180], "coal": [0, 55, 130, 225]}, dim="gen"
        ),
    ),
)

Ragged lengths are NaN-padded automatically. Breakpoints are auto-broadcast over remaining dimensions (e.g. time).

Disjunctive segments#

For disconnected operating regions (e.g. forbidden zones), use segments():

m.add_piecewise_formulation(
    (power, linopy.segments([(0, 0), (50, 80)])),
    (cost, linopy.segments([(0, 0), (125, 200)])),
)

The disjunctive formulation is selected automatically when breakpoints have a segment dimension. sign="<=" / ">=" also works here; the signed link is applied to the first tuple as usual.

N-variable linking#

Link any number of variables through shared breakpoints:

m.add_piecewise_formulation(
    (power, [0, 30, 60, 100]),
    (fuel, [0, 40, 85, 160]),
    (heat, [0, 25, 55, 95]),
)

With sign="==" (default) all variables are symmetric. With a non-equality sign the first tuple is the bounded output and the rest are forced to equality.

Formulation Methods#

Pass method="auto" (the default) and linopy picks the cheapest correct formulation based on sign, curvature and breakpoint layout:

  • 2-variable inequality on a convex/concave curve → lp (chord lines, no auxiliary variables)

  • All breakpoints monotonic → incremental

  • Otherwise → sos2

  • Disjunctive (segments) → always sos2 with binary segment selection

The resolved choice is exposed on the returned PiecewiseFormulation via .method (and .convexity when well-defined). An INFO-level log line explains the resolution whenever method="auto" is in play.

SOS2 (Convex Combination)#

Works for any breakpoint ordering. Introduces interpolation weights \(\lambda_i\) with an SOS2 adjacency constraint:

\[ \begin{align}\begin{aligned}&\sum_{i=0}^{n} \lambda_i = 1, \qquad \text{SOS2}(\lambda_0, \ldots, \lambda_n)\\&e_j = \sum_{i=0}^{n} \lambda_i \, B_{j,i} \quad \text{for each expression } j\end{aligned}\end{align} \]

The SOS2 constraint ensures at most two adjacent \(\lambda_i\) are non-zero, so every expression is interpolated within the same segment.

With sign != "==" the input tuples still use the equality above; the first tuple’s link is replaced by a one-sided e_1 \ \text{sign}\ \sum_i \lambda_i B_{1,i} constraint.

Note

SOS2 is handled via branch-and-bound, similar to integer variables. Prefer method="incremental" when breakpoints are monotonic.

m.add_piecewise_formulation((power, xp), (fuel, yp), method="sos2")

Incremental (Delta) Formulation#

For strictly monotonic breakpoints. Uses fill-fraction variables \(\delta_i\) with binary indicators \(z_i\):

\[ \begin{align}\begin{aligned}&\delta_i \in [0, 1], \quad z_i \in \{0, 1\}\\&\delta_{i+1} \le \delta_i, \quad z_{i+1} \le \delta_i, \quad \delta_i \le z_i\\&e_j = B_{j,0} + \sum_{i=1}^{n} \delta_i \, (B_{j,i} - B_{j,i-1})\end{aligned}\end{align} \]

With sign != "==" the same sign split as SOS2 applies: inputs use the equality above; the first tuple’s link uses the requested sign.

m.add_piecewise_formulation((power, xp), (fuel, yp), method="incremental")

Limitation: breakpoint sequences must be strictly monotonic.

LP (chord-line) Formulation#

For 2-variable inequality on a convex or concave curve. Adds one chord inequality per segment plus a domain bound — no auxiliary variables and no MIP relaxation:

\[ \begin{align}\begin{aligned}&y \ \text{sign}\ m_k \cdot x + c_k \quad \forall\ \text{segments } k\\&x_0 \le x \le x_n\end{aligned}\end{align} \]

where \(m_k = (y_{k+1} - y_k)/(x_{k+1} - x_k)\) and \(c_k = y_k - m_k\, x_k\). For concave \(f\) with sign="<=", the intersection of all chord inequalities equals the hypograph of \(f\) on its domain.

The LP dispatch requires curvature and sign to match: sign="<=" needs concave (or linear); sign=">=" needs convex (or linear). A mismatch is not just a loose bound — it describes the wrong region (see the Piecewise inequalities — the sign parameter). method="auto" detects this and falls back; method="lp" raises.

# y <= f(x) on a concave f — auto picks LP
m.add_piecewise_formulation((y, yp), (x, xp), sign="<=")

# Or explicitly:
m.add_piecewise_formulation((y, yp), (x, xp), sign="<=", method="lp")

Not supported with method="lp": sign="==", more than 2 tuples, and active. method="auto" falls back to SOS2/incremental in all three cases.

The underlying chord expressions are also exposed as a standalone helper, linopy.tangent_lines(x, x_pts, y_pts), which returns the per-segment chord as a LinearExpression with no variables created. Use it directly if you want to compose the chord bound with other constraints by hand, without the domain bound that method="lp" adds automatically.

Disjunctive (Disaggregated Convex Combination)#

For disconnected segments (gaps between operating regions). Binary indicators \(z_k\) select exactly one segment; SOS2 applies within it:

\[ \begin{align}\begin{aligned}&z_k \in \{0, 1\}, \quad \sum_{k} z_k = 1\\&\sum_{i} \lambda_{k,i} = z_k, \qquad e_j = \sum_{k} \sum_{i} \lambda_{k,i} \, B_{j,k,i}\end{aligned}\end{align} \]

No big-M constants are needed, giving a tight LP relaxation.

Disjunctive + ``sign``. sign="<=" / ">=" works here too, applied to the first tuple exactly as for the continuous methods. Because the disjunctive machinery already carries a per-segment binary, there is no curvature requirement on the segments — inequality is always exact on the hypograph (or epigraph) of the active segment, whatever its slope pattern. This makes disjunctive + sign a first-class tool for “bounded output on disconnected operating regions” that method="lp" cannot handle.

Advanced Features#

Active parameter (unit commitment)#

The active parameter gates the piecewise function with a binary variable. When active=0, all auxiliary variables (and thus the linked expressions) are forced to zero:

commit = m.add_variables(name="commit", binary=True, coords=[time])
m.add_piecewise_formulation(
    (power, [30, 60, 100]),
    (fuel, [40, 90, 170]),
    active=commit,
)

  • commit=1: power operates in [30, 100], fuel = f(power)

  • commit=0: power = 0, fuel = 0

Not supported with method="lp".

Note

With a non-equality sign, deactivation only pushes the signed bound to 0 — the complementary side comes from the output variable’s own lower/upper bound. Set lower=0 on naturally non-negative outputs (fuel, cost, heat) to pin the output to zero on deactivation. See the sign section above for details.

Auto-broadcasting#

Breakpoints are automatically broadcast to match expression dimensions — you don’t need expand_dims:

time = pd.Index([1, 2, 3], name="time")
x = m.add_variables(name="x", lower=0, upper=100, coords=[time])
y = m.add_variables(name="y", coords=[time])

# 1D breakpoints auto-expand to match x's time dimension
m.add_piecewise_formulation((x, [0, 50, 100]), (y, [0, 70, 150]))

NaN masking#

Trailing NaN values in breakpoints mask the corresponding lambda / delta variables (and, for LP, the corresponding chord constraints). This is useful for per-entity breakpoints with ragged lengths:

# gen1 has 3 breakpoints, gen2 has 2 (NaN-padded)
bp = linopy.breakpoints({"gen1": [0, 50, 100], "gen2": [0, 80]}, dim="gen")

Interior NaN values (gaps in the middle) are not supported and raise an error.

Generated variables and constraints#

Given a base name N (either user-supplied or auto-assigned like pwl0), each formulation creates a predictable set of names:

SOS2 (method="sos2"):

  • {N}_lambda — variable, interpolation weights

  • {N}_convex — constraint, sum(lambda) == 1 (or == active)

  • {N}_link — constraint, equality link (stacked inputs when sign != "==", all tuples when sign="==")

  • {N}_output_link — constraint, signed link on the first tuple (only when sign != "==" )

Incremental (method="incremental"):

  • {N}_delta — variable, fill fractions \(\delta_i\)

  • {N}_order_binary — variable, per-segment binaries \(z_i\)

  • {N}_delta_bound — constraint, \(\delta_i \le z_i\)

  • {N}_fill_order — constraint, \(\delta_{i+1} \le \delta_i\)

  • {N}_binary_order — constraint, \(z_{i+1} \le \delta_i\)

  • {N}_active_bound — constraint, \(\delta_i \le active\) (only when active is given)

  • {N}_link / {N}_output_link — same split as SOS2

LP (method="lp"):

  • {N}_chord — constraint, per-segment chord inequality

  • {N}_domain_lo, {N}_domain_hi — constraints, \(x_0 \le x \le x_n\)

  • no auxiliary variables

Disjunctive (segments(...) input):

  • {N}_segment_binary — variable, per-segment selectors \(z_k\)

  • {N}_select — constraint, sum(z_k) == 1 (or == active)

  • {N}_lambda — variable, within-segment weights

  • {N}_convex — constraint, per-segment \(\sum_i \lambda_{k,i} = z_k\)

  • {N}_link / {N}_output_link — same split as SOS2

See Also#