Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Piecewise Linear Constraints Tutorial#
add_piecewise_formulation links variables through shared breakpoint weights. Every section below stacks one feature on top of a small shared dispatch pattern — if you want the math, see the reference page. For inequality bounds and the LP chord formulation in depth, see the inequality bounds notebook.
The baseline we extend:
m.add_piecewise_formulation(
(power, [0, 30, 60, 100]),
(fuel, [0, 36, 84, 170]),
)
[1]:
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import linopy
time = pd.Index([1, 2, 3], name="time")
def plot_curve(
bp_x, bp_y, operating_x, operating_y, *, ax=None, xlabel="power", ylabel="fuel"
):
"""PWL curve with solver's operating points overlaid."""
ax = ax if ax is not None else plt.subplots(figsize=(4.5, 3.5))[1]
ax.plot(bp_x, bp_y, "o-", color="C0", label="breakpoints")
ax.plot(operating_x, operating_y, "D", color="C1", ms=10, label="solved", alpha=0.8)
ax.set(xlabel=xlabel, ylabel=ylabel)
ax.legend()
return ax
1. Getting started#
A gas turbine with a convex heat rate. Each (variable, 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 curve.
[2]:
demand = xr.DataArray([50, 80, 30], coords=[time])
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=100, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
x_pts = [0, 30, 60, 100]
y_pts = [0, 36, 84, 170]
pwf = m.add_piecewise_formulation((power, x_pts), (fuel, y_pts))
m.add_constraints(power == demand, name="demand")
m.add_objective(fuel.sum())
m.solve(reformulate_sos="auto", output_flag=False)
print(pwf) # inspect the auto-resolved method
m.solution[["power", "fuel"]].to_pandas()
/tmp/ipykernel_2227/4097225550.py:9: 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)`.
pwf = m.add_piecewise_formulation((power, x_pts), (fuel, y_pts))
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-0vct90zz.lp
Reading time = 0.00 seconds
obj: 30 rows, 24 columns, 69 nonzeros
Dual values of MILP couldn't be parsed
PiecewiseFormulation `pwl0` [time: 3] — incremental, concave
Variables:
* pwl0_delta (time, _breakpoint_piece)
* pwl0_order_binary (time, _breakpoint_piece)
Constraints:
* pwl0_delta_bound (time, _breakpoint_piece)
* pwl0_fill_order (time, _breakpoint_piece)
* pwl0_binary_order (time, _breakpoint_piece)
* pwl0_link (_breakpoint, _pwl_var, time)
[2]:
| power | fuel | |
|---|---|---|
| time | ||
| 1 | 50.0 | 68.0 |
| 2 | 80.0 | 127.0 |
| 3 | 30.0 | 36.0 |
[3]:
plot_curve(x_pts, y_pts, m.solution["power"].values, m.solution["fuel"].values);
2. Picking a method#
method="auto" (default) picks the cheapest correct formulation based on sign, curvature and breakpoint layout. The explicit options — "sos2", "incremental", "lp" — give the same optimum on equality cases where they all apply, so the choice is about cost (auxiliary variables, solver capability), not correctness.
method |
needs |
creates |
|---|---|---|
|
SOS2-capable solver |
lambdas (continuous) |
|
MIP solver, strictly monotonic breakpoints |
deltas (continuous) + binaries |
|
any LP solver |
no variables — requires |
Below: all applicable methods yield the same fuel dispatch on this convex curve.
[4]:
def solve_method(method):
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=100, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
m.add_piecewise_formulation((power, x_pts), (fuel, y_pts), method=method)
m.add_constraints(power == demand, name="demand")
m.add_objective(fuel.sum())
m.solve(reformulate_sos="auto", output_flag=False)
return m.solution["fuel"].to_pandas()
pd.DataFrame({m: solve_method(m) for m in ["auto", "sos2", "incremental"]})
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-qg4vv4g_.lp
Reading time = 0.00 seconds
obj: 30 rows, 24 columns, 69 nonzeros
Dual values of MILP couldn't be parsed
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-nwuo5w7o.lp
Reading time = 0.00 seconds
obj: 12 rows, 18 columns, 39 nonzeros
Dual values of MILP couldn't be parsed
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-o4p6bhpx.lp
Reading time = 0.00 seconds
obj: 30 rows, 24 columns, 69 nonzeros
Dual values of MILP couldn't be parsed
[4]:
| auto | sos2 | incremental | |
|---|---|---|---|
| time | |||
| 1 | 68.0 | 68.0 | 68.0 |
| 2 | 127.0 | 127.0 | 127.0 |
| 3 | 36.0 | 36.0 | 36.0 |
3. Disjunctive segments — gaps in the operating range#
When operating regions are disconnected (a diesel generator that is either off or running in [50, 80] MW, never in between), use segments() instead of breakpoints(). A binary picks which segment is active; inside it SOS2 interpolates as usual.
[5]:
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=80, coords=[time])
cost = m.add_variables(name="cost", lower=0, coords=[time])
backup = m.add_variables(name="backup", lower=0, coords=[time])
m.add_piecewise_formulation(
(power, linopy.segments([(0, 0), (50, 80)])), # two disjoint segments
(cost, linopy.segments([(0, 0), (125, 200)])),
)
m.add_constraints(power + backup == xr.DataArray([15, 60, 75], coords=[time]))
m.add_objective(cost.sum() + 10 * backup.sum())
m.solve(reformulate_sos="auto", output_flag=False)
m.solution[["power", "cost", "backup"]].to_pandas()
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-wtnsh1di.lp
Reading time = 0.00 seconds
obj: 18 rows, 27 columns, 48 nonzeros
Dual values of MILP couldn't be parsed
[5]:
| power | cost | backup | |
|---|---|---|---|
| time | |||
| 1 | 0.0 | 0.0 | 15.0 |
| 2 | 60.0 | 150.0 | 0.0 |
| 3 | 75.0 | 187.5 | 0.0 |
At t=1 the 15 MW demand falls in the forbidden zone; the unit sits at 0 and backup fills the gap.
4. Inequality bounds — per-tuple sign#
Append a third tuple element ("<=" or ">=") to mark a single expression as bounded by the piecewise curve instead of pinned to it. The other tuples stay on the curve. The 2-variable hypograph (y ≤ f(x)) and epigraph (y ≥ f(x)) are the canonical cases.
On a concave curve with <= (or convex with >=), method="auto" dispatches to a pure LP chord formulation — no binaries, no SOS2.
At most one tuple may carry a non-equality sign. With 3 or more tuples, all signs must be "=="; the multi-input bounded case is reserved for a future bivariate / triangulated piecewise API.
See the inequality bounds notebook for mismatched curvature, auto-dispatch fallbacks, and more geometry.
[6]:
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=120, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
# concave curve: diminishing marginal fuel per MW
x_pts = [0, 50, 90, 120]
y_pts = [0, 40, 80, 120]
pwf = m.add_piecewise_formulation(
(fuel, x_pts, "<="), # bounded above by the curve
(power, y_pts), # pinned to the curve
)
m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))
m.add_objective(-fuel.sum()) # push fuel against the bound
m.solve(reformulate_sos="auto", output_flag=False)
print(f"resolved method={pwf.method}, curvature={pwf.convexity}")
m.solution[["power", "fuel"]].to_pandas()
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-uqqr9tzd.lp
Reading time = 0.00 seconds
obj: 18 rows, 6 columns, 27 nonzeros
resolved method=lp, curvature=concave
[6]:
| power | fuel | |
|---|---|---|
| time | ||
| 1 | 30.0 | 37.5 |
| 2 | 80.0 | 90.0 |
| 3 | 100.0 | 105.0 |
[7]:
# x_pts are fuel breakpoints, y_pts are power breakpoints — swap so power is on the x-axis
plot_curve(y_pts, x_pts, m.solution["power"].values, m.solution["fuel"].values);
5. Unit commitment — active#
A binary variable gates the whole formulation. active=0 forces the PWL variables (and thus all linked outputs) to zero. Combined with the natural lower=0 on cost/fuel/heat, this gives a clean on/off coupling:
active=1: the unit operates in its full range, outputs tied to the curve.active=0:power = 0,fuel = 0.
[8]:
m = linopy.Model()
p_min, p_max = 30, 100
power = m.add_variables(name="power", lower=0, upper=p_max, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
backup = m.add_variables(name="backup", lower=0, coords=[time])
commit = m.add_variables(name="commit", binary=True, coords=[time])
x_pts = [p_min, 60, p_max]
y_pts = [40, 90, 170]
m.add_piecewise_formulation(
(power, x_pts),
(fuel, y_pts),
active=commit,
)
# demand below p_min at t=1 — commit must be 0 and backup covers it
m.add_constraints(power + backup == xr.DataArray([15, 80, 40], coords=[time]))
m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())
m.solve(reformulate_sos="auto", output_flag=False)
m.solution[["commit", "power", "fuel", "backup"]].to_pandas()
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-nmqh3fil.lp
Reading time = 0.00 seconds
obj: 27 rows, 24 columns, 66 nonzeros
Dual values of MILP couldn't be parsed
[8]:
| commit | power | fuel | backup | |
|---|---|---|---|---|
| time | ||||
| 1 | 0.0 | 0.0 | 0.000000 | 15.0 |
| 2 | 1.0 | 80.0 | 130.000000 | 0.0 |
| 3 | 1.0 | 40.0 | 56.666667 | 0.0 |
[9]:
plot_curve(x_pts, y_pts, m.solution["power"].values, m.solution["fuel"].values);
6. N-variable linking — CHP plant#
More than two variables can share the same interpolation — useful for combined heat-and-power plants where power, fuel and heat are all functions of a single operating point.
[10]:
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=100, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
heat = m.add_variables(name="heat", lower=0, coords=[time])
x_pts = [0, 30, 60, 100]
y_pts = [0, 40, 85, 160]
z_pts = [0, 25, 55, 95]
m.add_piecewise_formulation(
(power, x_pts),
(fuel, y_pts),
(heat, z_pts),
)
m.add_constraints(fuel == xr.DataArray([20, 100, 160], coords=[time]))
m.add_objective(power.sum())
m.solve(reformulate_sos="auto", output_flag=False)
m.solution[["power", "fuel", "heat"]].to_pandas().round(2)
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-en9sh8dm.lp
Reading time = 0.00 seconds
obj: 33 rows, 27 columns, 81 nonzeros
Dual values of MILP couldn't be parsed
[10]:
| power | fuel | heat | |
|---|---|---|---|
| time | |||
| 1 | 15.0 | 20.0 | 12.5 |
| 2 | 68.0 | 100.0 | 63.0 |
| 3 | 100.0 | 160.0 | 95.0 |
[11]:
fig, axes = plt.subplots(1, 2, figsize=(8, 3))
plot_curve(
x_pts, y_pts, m.solution["power"].values, m.solution["fuel"].values, ax=axes[0]
)
plot_curve(
x_pts,
z_pts,
m.solution["power"].values,
m.solution["heat"].values,
ylabel="heat",
ax=axes[1],
);
7. Per-entity breakpoints — a fleet of generators#
Pass a dict to breakpoints() with entity names as keys for different curves per entity. Ragged lengths are NaN-padded automatically, and breakpoints broadcast over any remaining dimensions (here, time).
[12]:
gens = pd.Index(["gas", "coal"], name="gen")
x_gen = linopy.breakpoints(
{"gas": [0, 30, 60, 100], "coal": [0, 50, 100, 150]}, dim="gen"
)
y_gen = linopy.breakpoints(
{"gas": [0, 40, 90, 180], "coal": [0, 55, 130, 225]}, dim="gen"
)
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=150, coords=[gens, time])
fuel = m.add_variables(name="fuel", lower=0, coords=[gens, time])
m.add_piecewise_formulation((power, x_gen), (fuel, y_gen))
m.add_constraints(power.sum("gen") == xr.DataArray([80, 120, 50], coords=[time]))
m.add_objective(fuel.sum())
m.solve(reformulate_sos="auto", output_flag=False)
m.solution[["power", "fuel"]].to_dataframe()
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-4_fk7uqe.lp
Reading time = 0.00 seconds
obj: 57 rows, 48 columns, 138 nonzeros
Dual values of MILP couldn't be parsed
[12]:
| power | fuel | ||
|---|---|---|---|
| gen | time | ||
| gas | 1 | 30.0 | 40.0 |
| 2 | 30.0 | 40.0 | |
| 3 | 0.0 | 0.0 | |
| coal | 1 | 50.0 | 55.0 |
| 2 | 90.0 | 115.0 | |
| 3 | 50.0 | 55.0 |
8. Specifying with slopes — Slopes#
When marginal costs (slopes) are more natural than absolute y-values, wrap them in linopy.Slopes. The x grid is borrowed from the sibling tuple — no need to repeat it. Same curve as section 1:
[13]:
m = linopy.Model()
power = m.add_variables(name="power", lower=0, upper=100, coords=[time])
fuel = m.add_variables(name="fuel", lower=0, coords=[time])
m.add_piecewise_formulation(
(power, [0, 30, 60, 100]),
(fuel, linopy.Slopes([1.2, 1.6, 2.15], y0=0)),
)
m.add_constraints(power == demand, name="demand")
m.add_objective(fuel.sum())
m.solve(reformulate_sos="auto", output_flag=False)
m.solution[["power", "fuel"]].to_pandas()
/tmp/ipykernel_2227/3771982636.py:7: EvolvingAPIWarning: piecewise: Slopes is a new API; the constructor signature and the dispatch rules for inheriting an x grid from sibling tuples may be refined in minor releases.
(fuel, linopy.Slopes([1.2, 1.6, 2.15], y0=0)),
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-nbofw_sr.lp
Reading time = 0.00 seconds
obj: 30 rows, 24 columns, 69 nonzeros
Dual values of MILP couldn't be parsed
[13]:
| power | fuel | |
|---|---|---|
| time | ||
| 1 | 50.0 | 68.0 |
| 2 | 80.0 | 127.0 |
| 3 | 30.0 | 36.0 |