From cdc2df6369ba8f102c25c429eff5d05d08870540 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:17:05 +0200
Subject: [PATCH 01/29] docs(piecewise): tighten language and fix small
 inaccuracies

- "breakpoint" is now framed as one knot defined by the i-th entry of
  each tuple, so the N-tuple linking case isn't (x, y)-specific.
- "pinned" is reframed as a sign role, with a brief breakdown of what
  it actually constrains: joint-on-curve only with 2+ pinned tuples;
  with one bounded + one pinned, the pinned axis collapses to a
  [x_min, x_max] domain box (LP enforces it directly, SOS2/incremental
  via the weight link).
- Disjunctive auxiliary-variables cell corrected to
  "Continuous + binary + SOS2".
- SOS2/disjunctive solver-requirement cells now mention the Big-M
  reformulation path and link to :ref:`sos-reformulation`.
- LP domain bound math uses x_min/x_max (descending grids accepted).
- Per-tuple-sign formulation math switched to a method-agnostic
  W_j(weights, B) so the section covers both SOS2 (lambda) and
  incremental (delta) accurately.
- "Two factories" lead-in widened to three building blocks so Slopes
  isn't hidden in the code block.
- Quick Start inequality bounds heat (a curtailable output) instead of
  fuel, matching the doc's own "choice of bounded tuple" guidance.
- Tutorials: add output_flag=False to m.solve(...) calls so HiGHS's
  banner/progress doesn't clutter the nbsphinx-rendered output;
  linopy's INFO logs (including auto-dispatch resolution) are kept.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        | 126 +++++++++++++-------
 examples/piecewise-inequality-bounds.ipynb  |   2 +-
 examples/piecewise-linear-constraints.ipynb |  46 +++----
 3 files changed, 105 insertions(+), 69 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index e364988c..c96ad54a 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -9,7 +9,11 @@ production functions within a linear programming framework.
 
 **Terminology used in this page:**
 
-- **breakpoint** — an :math:`(x, y)` knot where the slope can change.
+- **breakpoint** — a knot on the piecewise curve where the slope can
+  change.  Each tuple supplies a 1-D breakpoint array; the ``i``-th
+  entries across all tuples, taken together, define one knot.  For two
+  tuples this is the usual :math:`(x, y)` pair; for three or more, an
+  ``N``-dim knot.
 - **piece** — a linear part between two adjacent breakpoints on a single
   connected curve.  ``n`` breakpoints define ``n − 1`` pieces.
 - **segment** — a *disjoint* operating region in the disjunctive
@@ -45,19 +49,21 @@ Quick Start
 
 .. code-block:: python
 
-    # 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).
+    # heat <= f(power).  "auto" picks the cheapest correct formulation:
+    # pure LP (chord constraints) when curvature matches the sign,
+    # SOS2/incremental otherwise.  Bound a curtailable output so
+    # undershooting the curve is physically realisable — see *Choice of
+    # bounded tuple* below.
     m.add_piecewise_formulation(
-        (fuel, [0, 20, 30, 35], "<="),  # bounded by the curve
-        (power, [0, 10, 20, 30]),  # pinned to the curve
+        (heat, [0, 20, 30, 35], "<="),
+        (power, [0, 10, 20, 30]),
     )
 
-Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with its
-breakpoint values, and optionally marks it as bounded by the curve (``"<="``
-or ``">="``) instead of pinned to it.  All tuples share interpolation weights,
-so at any feasible point every variable corresponds to the *same* point on
-the piecewise curve.
+Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with
+its breakpoint values.  The optional sign (default ``"=="``) is ``"<="``
+or ``">="`` to mark that expression as bounded by the curve.  With every
+sign ``"=="``, all tuples land on the same point of the piecewise curve
+— see *Per-tuple sign* below for the geometry of the inequality cases.
 
 
 API
@@ -69,7 +75,7 @@ API
 .. code-block:: python
 
     m.add_piecewise_formulation(
-        (expr1, breakpoints1),  # pinned (sign defaults to "==")
+        (expr1, breakpoints1),  # sign defaults to "==" (pinned role)
         (expr2, breakpoints2, "<="),  # or with an explicit sign
         ...,
         method="auto",  # "auto", "sos2", "incremental", or "lp"
@@ -77,53 +83,72 @@ API
         name=None,  # base name for generated variables/constraints
     )
 
-Creates auxiliary variables and constraints that enforce either a joint
-equality (all tuples on the curve, the default) or a one-sided bound
-(at most one tuple bounded by the curve, the rest pinned).
+Adds constraints — and, depending on the resolved method, auxiliary
+variables — for either an all-equality joint (every tuple at the same
+point on the curve, the default) or a one-sided bound on a single
+tuple.  The pure-LP path adds chord and domain constraints only; SOS2,
+incremental, and disjunctive also add interpolation weights and/or
+binaries (see *Formulation Methods* below).
 
-``breakpoints`` and ``segments``
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Breakpoint inputs
+~~~~~~~~~~~~~~~~~
 
-Two factories with distinct geometric meaning:
+Three building blocks with distinct geometric meaning — two factories
+and one class:
 
 - ``breakpoints()`` — values along a single **connected** curve.  Linear
   pieces between adjacent breakpoints are interpolated continuously.
 - ``segments()`` — **disjoint** operating regions with gaps between them
   (e.g. forbidden zones).  Builds a 2-D array consumed by the
   *disjunctive* formulation, where exactly one region is active at a time.
+- :class:`~linopy.Slopes` — per-piece slopes plus an initial ``y0``,
+  deferred until an x grid is supplied.  Inside
+  ``add_piecewise_formulation`` the x grid is borrowed from a sibling
+  tuple; standalone, call :meth:`~linopy.Slopes.to_breakpoints`.
 
 .. code-block:: python
 
     linopy.breakpoints([0, 50, 100])  # connected
     linopy.breakpoints({"gen1": [0, 50], "gen2": [0, 80]}, dim="gen")  # per-entity
-    linopy.Slopes(
-        [1.2, 1.4], y0=0
-    )  # from slopes (deferred — pairs with a sibling tuple)
     linopy.segments([(0, 10), (50, 100)])  # two disjoint regions
     linopy.segments({"gen1": [(0, 10)], "gen2": [(0, 80)]}, dim="gen")
+    linopy.Slopes([1.2, 1.4], y0=0)  # deferred — pairs with a sibling tuple
 
 
 Per-tuple sign — equality vs inequality
 ----------------------------------------
 
-By default each tuple's expression is **pinned** to the piecewise curve.
-Pass a third tuple element (``"<="`` or ``">="``) to mark a single
-expression as **bounded** by the curve — it can undershoot (``"<="``) or
-overshoot (``">="``) the interpolated value, while every other tuple
-stays pinned.
+Each tuple's optional third element is a sign:
+
+- ``"=="`` (default) — **pinned**: the tuple enters as an equality.
+- ``"<="`` / ``">="`` — **bounded**: the expression undershoots /
+  overshoots the curve.
+
+These describe the tuple's *role*, not a geometric property of the
+variable.  What "pinned" actually constrains depends on the count:
+
+- **All tuples pinned (default).**  Shared interpolation weights put
+  the joint :math:`(e_1, \ldots, e_N)` exactly on the curve.
+- **One bounded + one pinned (2 tuples).**  The joint :math:`(x, y)`
+  lies in the hypograph / epigraph (see *Geometry* below).  A single
+  coordinate can't locate a curve point, so the pinned axis's marginal
+  feasible set is just :math:`[x_{\min}, x_{\max}]` — same in LP
+  (enforced directly) and SOS2/incremental (enforced via the weight
+  link).
+- **One bounded + 3+ tuples.**  Not supported (see restrictions below).
 
 .. code-block:: python
 
-    # Joint equality (default): both expressions on the curve.
+    # All-equality: joint (x, y) on the curve.
     m.add_piecewise_formulation((y, y_pts), (x, x_pts))
 
-    # Bounded above: y <= f(x), x pinned.
+    # Bounded: joint (x, y) in the hypograph — y ≤ f(x), x ∈ [x_min, x_max].
     m.add_piecewise_formulation((y, y_pts, "<="), (x, x_pts))
 
-    # Bounded below: y >= f(x), x pinned.
+    # Bounded: joint (x, y) in the epigraph — y ≥ f(x), x ∈ [x_min, x_max].
     m.add_piecewise_formulation((y, y_pts, ">="), (x, x_pts))
 
-    # 3-variable equality (CHP heat/power/fuel): all three on one curve.
+    # 3-variable all-equality (CHP): joint (power, fuel, heat) on the curve.
     m.add_piecewise_formulation((power, p_pts), (fuel, f_pts), (heat, h_pts))
 
 **Restrictions (current):**
@@ -137,16 +162,20 @@ https://github.com/PyPSA/linopy/issues so we can scope it properly.
 
 **Formulation.**  For methods that introduce shared interpolation
 weights (SOS2 and incremental — see below), only the link constraint
-between the weights and the bounded expression changes.  Pinned tuples
-:math:`j` keep the equality, and the bounded tuple :math:`b` flips to
-the requested sign:
+between the weights and the bounded expression changes.  Write the
+method-specific weighted sum of breakpoints for tuple :math:`j` as
+:math:`W_j(\text{weights}, B)` — the explicit form is :math:`\sum_i
+\lambda_i B_{j,i}` for SOS2 and :math:`B_{j,0} + \sum_i \delta_i (B_{j,i}
+- B_{j,i-1})` for incremental (see the method sections below).  Pinned
+tuples :math:`j` keep the equality, and the bounded tuple :math:`b` flips
+to the requested sign:
 
 .. math::
 
-   &e_j = \sum_{i=0}^{n} \lambda_i \, B_{j,i}
+   &e_j = W_j(\text{weights}, B)
    \quad \text{(pinned, } j \ne b \text{)}
 
-   &e_b \ \text{sign}\ \sum_{i=0}^{n} \lambda_i \, B_{b,i}
+   &e_b \ \text{sign}\ W_b(\text{weights}, B)
    \quad \text{(bounded)}
 
 Internally this shows up as a stacked ``*_link`` equality covering the
@@ -330,11 +359,15 @@ formulation based on ``sign``, curvature and breakpoint layout:
   no auxiliary variables)
 - **All breakpoints monotonic** → ``incremental``
 - **Otherwise** → ``sos2``
-- **Disjunctive (segments)** → always ``sos2`` with binary segment selection
+- **Disjunctive (segments)** → SOS2 applied per segment with binary
+  segment selection (the disjunctive formulation in the table below).
 
 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.
+``.method`` (and ``.convexity`` when well-defined).  Disjunctive
+formulations report ``method="sos2"`` even though their structure is the
+per-segment variant — the table below treats it as a separate column for
+clarity.  An ``INFO``-level log line explains the resolution whenever
+``method="auto"`` is in play.
 
 At-a-glance comparison:
 
@@ -376,7 +409,7 @@ At-a-glance comparison:
      - **None**
      - Continuous + binary
      - Continuous + SOS2
-     - Binary + SOS2
+     - Continuous + binary + SOS2
    * - ``active=`` supported
      - No
      - Yes
@@ -385,8 +418,8 @@ At-a-glance comparison:
    * - Solver requirement
      - **Any LP solver**
      - MIP-capable
-     - SOS2-capable
-     - SOS2 + MIP
+     - SOS2-capable (or MIP via :ref:`Big-M reformulation <sos-reformulation>`)
+     - SOS2 + MIP (or MIP via :ref:`Big-M reformulation <sos-reformulation>`)
 
 LP (chord-line) Formulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -400,12 +433,15 @@ no MIP relaxation:
    &y \ \text{sign}\ m_k \cdot x + c_k
    \quad \forall\ \text{pieces } k
 
-   &x_0 \le x \le x_n
+   &x_{\min} \le x \le x_{\max}
 
 where :math:`m_k = (y_{k+1} - y_k)/(x_{k+1} - x_k)` and
-:math:`c_k = y_k - m_k\, x_k`.  For concave :math:`f` with ``sign="<="``,
-the intersection of all chord inequalities equals the hypograph of
-:math:`f` on its domain.
+:math:`c_k = y_k - m_k\, x_k`.  The domain bound uses
+:math:`x_{\min}` and :math:`x_{\max}` rather than the first/last
+breakpoint so that descending x grids work too — strictly-monotonic
+breakpoints are accepted in either order.  For concave :math:`f` with
+``sign="<="``, the intersection of all chord inequalities equals the
+hypograph of :math:`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
diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index a79c612d..bb597e5b 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -98,7 +98,7 @@
     }
    },
    "outputs": [],
-   "source": "def solve(method, power_val):\n    m = linopy.Model()\n    power = m.add_variables(lower=0, upper=30, name=\"power\")\n    fuel = m.add_variables(lower=0, upper=40, name=\"fuel\")\n    m.add_piecewise_formulation(\n        (fuel, y_pts, \"<=\"),  # bounded\n        (power, x_pts),  # pinned\n        method=method,\n    )\n    m.add_constraints(power == power_val)\n    m.add_objective(-fuel)  # maximise fuel to push against the bound\n    m.solve()\n    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n\n\nfor method in [\"lp\", \"sos2\", \"incremental\"]:\n    fuel_val, vars_, cons_ = solve(method, 15)\n    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
+   "source": "def solve(method, power_val):\n    m = linopy.Model()\n    power = m.add_variables(lower=0, upper=30, name=\"power\")\n    fuel = m.add_variables(lower=0, upper=40, name=\"fuel\")\n    m.add_piecewise_formulation(\n        (fuel, y_pts, \"<=\"),  # bounded\n        (power, x_pts),  # pinned\n        method=method,\n    )\n    m.add_constraints(power == power_val)\n    m.add_objective(-fuel)  # maximise fuel to push against the bound\n    m.solve(output_flag=False)\n    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n\n\nfor method in [\"lp\", \"sos2\", \"incremental\"]:\n    fuel_val, vars_, cons_ = solve(method, 15)\n    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
   },
   {
    "cell_type": "markdown",
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 392ca8f1..978caa49 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "# Piecewise Linear Constraints Tutorial\n",
     "\n",
-    "`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](piecewise-linear-constraints). For inequality bounds and the LP chord formulation in depth, see the [inequality bounds notebook](piecewise-inequality-bounds-tutorial).\n",
+    "`add_piecewise_formulation` links variables through shared breakpoint weights. Every section below stacks one feature on top of a small shared dispatch pattern \u2014 if you want the math, see the [reference page](piecewise-linear-constraints). For inequality bounds and the LP chord formulation in depth, see the [inequality bounds notebook](piecewise-inequality-bounds-tutorial).\n",
     "\n",
     "The baseline we extend:\n",
     "\n",
@@ -81,7 +81,7 @@
     "pwf = m.add_piecewise_formulation((power, x_pts), (fuel, y_pts))\n",
     "m.add_constraints(power == demand, name=\"demand\")\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
@@ -107,13 +107,13 @@
    "source": [
     "## 2. Picking a method\n",
     "\n",
-    "`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.\n",
+    "`method=\"auto\"` (default) picks the cheapest correct formulation based on sign, curvature and breakpoint layout. The explicit options \u2014 `\"sos2\"`, `\"incremental\"`, `\"lp\"` \u2014 give the same optimum on equality cases where they all apply, so the choice is about **cost** (auxiliary variables, solver capability), not correctness.\n",
     "\n",
     "| method | needs | creates |\n",
     "|---|---|---|\n",
     "| `sos2` | SOS2-capable solver | lambdas (continuous) |\n",
     "| `incremental` | MIP solver, strictly monotonic breakpoints | deltas (continuous) + binaries |\n",
-    "| `lp` | any LP solver | no variables — requires `sign != \"==\"`, 2 tuples, matching curvature |\n",
+    "| `lp` | any LP solver | no variables \u2014 requires `sign != \"==\"`, 2 tuples, matching curvature |\n",
     "\n",
     "Below: all applicable methods yield the same fuel dispatch on this convex curve."
    ]
@@ -136,7 +136,7 @@
     "    m.add_piecewise_formulation((power, x_pts), (fuel, y_pts), method=method)\n",
     "    m.add_constraints(power == demand, name=\"demand\")\n",
     "    m.add_objective(fuel.sum())\n",
-    "    m.solve(reformulate_sos=\"auto\")\n",
+    "    m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "    return m.solution[\"fuel\"].to_pandas()\n",
     "\n",
     "\n",
@@ -147,7 +147,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 3. Disjunctive segments — gaps in the operating range\n",
+    "## 3. Disjunctive segments \u2014 gaps in the operating range\n",
     "\n",
     "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."
    ]
@@ -175,7 +175,7 @@
     ")\n",
     "m.add_constraints(power + backup == xr.DataArray([15, 60, 75], coords=[time]))\n",
     "m.add_objective(cost.sum() + 10 * backup.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"cost\", \"backup\"]].to_pandas()"
    ]
   },
@@ -190,11 +190,11 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 4. Inequality bounds — per-tuple sign\n",
+    "## 4. Inequality bounds \u2014 per-tuple sign\n",
     "\n",
-    "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.\n",
+    "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 \u2264 f(x)`) and epigraph (`y \u2265 f(x)`) are the canonical cases.\n",
     "\n",
-    "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
+    "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation \u2014 no binaries, no SOS2.\n",
     "\n",
     "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.\n",
     "\n",
@@ -225,7 +225,7 @@
     ")\n",
     "m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))\n",
     "m.add_objective(-fuel.sum())  # push fuel against the bound\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
@@ -237,7 +237,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# x_pts are fuel breakpoints, y_pts are power breakpoints — swap so power is on the x-axis\n",
+    "# x_pts are fuel breakpoints, y_pts are power breakpoints \u2014 swap so power is on the x-axis\n",
     "plot_curve(y_pts, x_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
    ]
   },
@@ -245,7 +245,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 5. Unit commitment — `active`\n",
+    "## 5. Unit commitment \u2014 `active`\n",
     "\n",
     "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:\n",
     "\n",
@@ -279,10 +279,10 @@
     "    (fuel, y_pts),\n",
     "    active=commit,\n",
     ")\n",
-    "# demand below p_min at t=1 — commit must be 0 and backup covers it\n",
+    "# demand below p_min at t=1 \u2014 commit must be 0 and backup covers it\n",
     "m.add_constraints(power + backup == xr.DataArray([15, 80, 40], coords=[time]))\n",
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"commit\", \"power\", \"fuel\", \"backup\"]].to_pandas()"
    ]
   },
@@ -299,9 +299,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 6. N-variable linking — CHP plant\n",
+    "## 6. N-variable linking \u2014 CHP plant\n",
     "\n",
-    "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."
+    "More than two variables can share the same interpolation \u2014 useful for combined heat-and-power plants where power, fuel and heat are all functions of a single operating point."
    ]
   },
   {
@@ -330,7 +330,7 @@
     ")\n",
     "m.add_constraints(fuel == xr.DataArray([20, 100, 160], coords=[time]))\n",
     "m.add_objective(power.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\", \"heat\"]].to_pandas().round(2)"
    ]
   },
@@ -358,7 +358,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 7. Per-entity breakpoints — a fleet of generators\n",
+    "## 7. Per-entity breakpoints \u2014 a fleet of generators\n",
     "\n",
     "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`)."
    ]
@@ -388,7 +388,7 @@
     "m.add_piecewise_formulation((power, x_gen), (fuel, y_gen))\n",
     "m.add_constraints(power.sum(\"gen\") == xr.DataArray([80, 120, 50], coords=[time]))\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\"]].to_dataframe()"
    ]
   },
@@ -396,9 +396,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 8. Specifying with slopes — `Slopes`\n",
+    "## 8. Specifying with slopes \u2014 `Slopes`\n",
     "\n",
-    "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:"
+    "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 \u2014 no need to repeat it. Same curve as section 1:"
    ]
   },
   {
@@ -417,7 +417,7 @@
     ")\n",
     "m.add_constraints(power == demand, name=\"demand\")\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\")\n",
+    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
    ]

From e7cd756e0d58e19ea133f7f011ac5fb8820a37b7 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:22:39 +0200
Subject: [PATCH 02/29] docs(piecewise): rename tutorial heading to "Creating
 Piecewise Linear Constraints"

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 978caa49..af79fee1 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Piecewise Linear Constraints Tutorial\n",
+    "# Creating Piecewise Linear Constraints\n",
     "\n",
     "`add_piecewise_formulation` links variables through shared breakpoint weights. Every section below stacks one feature on top of a small shared dispatch pattern \u2014 if you want the math, see the [reference page](piecewise-linear-constraints). For inequality bounds and the LP chord formulation in depth, see the [inequality bounds notebook](piecewise-inequality-bounds-tutorial).\n",
     "\n",

From c5104be9911104589bfd4e47b03385960cdf0552 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:23:01 +0200
Subject: [PATCH 03/29] docs(piecewise): rename inequality-bounds tutorial
 heading to match

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index bb597e5b..76766687 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Piecewise inequalities \u2014 per-tuple sign\n",
+    "# Creating Piecewise Inequality Bounds\n",
     "\n",
     "`add_piecewise_formulation` accepts an optional third tuple element, `\"<=\"` or `\">=\"`, that marks one expression as **bounded** by the piecewise curve instead of pinned to it:\n",
     "\n",

From 716eb714539b11b44089eff8a1f87e3b51b6f769 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:33:01 +0200
Subject: [PATCH 04/29] docs(piecewise): polish both tutorials for readability
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Tut 1 intro now opens with a one-line description plus an 8-section
  roadmap, replacing the jargon-y "stacks one feature on top of a
  small shared dispatch pattern" line.
- Tut 1 §3 explains the degenerate (0, 0) "off" segment instead of
  leaving it for the reader to puzzle out.
- Tut 1 §4 inequality example now bounds `heat` (a curtailable output)
  on a concave curve, matching the rst doc's "choice of bounded tuple"
  guidance.  Variables use intuitive `power_pts` / `heat_pts` names so
  the plot cell no longer needs a "swap to put power on the x-axis"
  comment.
- Tut 1 ends with a "When to use what" table cross-linking to the rst
  reference page and to the inequality-bounds tutorial.
- Tut 2 intro leads with the one-sided-bound motivation and the
  pure-LP pay-off before the API snippet, and the "Tuple roles" table
  matches the rst's precision: with one bounded + one equality tuple,
  the equality tuple's marginal feasible set is just its breakpoint
  domain (not "lies exactly on the curve").
- Tut 2's summary cross-links back to the linear-constraints tutorial.
- Both notebooks: removed remaining "pinned" wording where it implied
  geometric on-curve placement instead of a sign role.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb  | 60 ++++++++++---
 examples/piecewise-linear-constraints.ipynb | 97 ++++++++++++++-------
 2 files changed, 113 insertions(+), 44 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 76766687..682a3656 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -6,24 +6,26 @@
    "source": [
     "# Creating Piecewise Inequality Bounds\n",
     "\n",
-    "`add_piecewise_formulation` accepts an optional third tuple element, `\"<=\"` or `\">=\"`, that marks one expression as **bounded** by the piecewise curve instead of pinned to it:\n",
+    "When you only need a one-sided bound from a piecewise curve \u2014 `y \u2264 f(x)` for a concave upper envelope, `y \u2265 f(x)` for a convex lower envelope \u2014 `add_piecewise_formulation` accepts an optional sign as the third tuple element:\n",
     "\n",
     "```python\n",
     "m.add_piecewise_formulation(\n",
     "    (fuel,  y_pts, \"<=\"),   # bounded above by the curve\n",
-    "    (power, x_pts),         # pinned to the curve\n",
+    "    (power, x_pts),         # equality role\n",
     ")\n",
     "```\n",
     "\n",
-    "This notebook walks through the geometry, the curvature \u00d7 sign matching that lets `method=\"auto\"` skip MIP machinery entirely, and the feasible regions produced by each method (LP, SOS2, incremental). For the formulation math see the [reference page](piecewise-linear-constraints).\n",
+    "The pay-off is a pure-LP encoding when the curve's curvature matches the sign \u2014 no SOS2, no binaries. This notebook covers the geometry of the feasible region, the curvature \u00d7 sign combinations that unlock the LP path, and what happens when they don't match.\n",
     "\n",
-    "## Key points\n",
+    "For the formulation math see the [reference page](piecewise-linear-constraints); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial).\n",
     "\n",
-    "| Tuple form | Behaviour |\n",
-    "|---|---|\n",
-    "| `(expr, breaks)` | Pinned: `expr` lies exactly on the curve. |\n",
-    "| `(expr, breaks, \"<=\")` | Bounded above: `expr \u2264 f(other tuples)`. |\n",
-    "| `(expr, breaks, \">=\")` | Bounded below: `expr \u2265 f(other tuples)`. |\n",
+    "## Tuple roles\n",
+    "\n",
+    "| Tuple form | Role | What it constrains |\n",
+    "|---|---|---|\n",
+    "| `(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]` \u2014 one coordinate alone can't locate a curve point. |\n",
+    "| `(expr, breaks, \"<=\")` | bounded above | `expr \u2264 f(other tuples)`. |\n",
+    "| `(expr, breaks, \">=\")` | bounded below | `expr \u2265 f(other tuples)`. |\n",
     "\n",
     "Currently at most one tuple may carry a non-equality sign, and 3+ tuples must all be equality. Multi-bounded and N\u22653 inequality cases aren't supported yet \u2014 if you have a concrete use case, please open an issue at https://github.com/PyPSA/linopy/issues so we can scope it properly."
    ]
@@ -80,7 +82,7 @@
     "With one tuple bounded `<=` and our concave curve, the three methods give the **same** feasible region within `[x_0, x_n]`:\n",
     "\n",
     "- **`method=\"lp\"`** \u2014 tangent lines + domain bounds. No auxiliary variables.\n",
-    "- **`method=\"sos2\"`** \u2014 lambdas + SOS2 + split link (pinned equality, bounded signed). Solver picks the piece.\n",
+    "- **`method=\"sos2\"`** \u2014 lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
     "- **`method=\"incremental\"`** \u2014 delta fractions + binaries + split link. Same mathematics, MIP encoding instead of SOS2.\n",
     "\n",
     "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable \u2014 it's always preferable because it's pure LP.\n",
@@ -98,12 +100,35 @@
     }
    },
    "outputs": [],
-   "source": "def solve(method, power_val):\n    m = linopy.Model()\n    power = m.add_variables(lower=0, upper=30, name=\"power\")\n    fuel = m.add_variables(lower=0, upper=40, name=\"fuel\")\n    m.add_piecewise_formulation(\n        (fuel, y_pts, \"<=\"),  # bounded\n        (power, x_pts),  # pinned\n        method=method,\n    )\n    m.add_constraints(power == power_val)\n    m.add_objective(-fuel)  # maximise fuel to push against the bound\n    m.solve(output_flag=False)\n    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n\n\nfor method in [\"lp\", \"sos2\", \"incremental\"]:\n    fuel_val, vars_, cons_ = solve(method, 15)\n    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
+   "source": [
+    "def solve(method, power_val):\n",
+    "    m = linopy.Model()\n",
+    "    power = m.add_variables(lower=0, upper=30, name=\"power\")\n",
+    "    fuel = m.add_variables(lower=0, upper=40, name=\"fuel\")\n",
+    "    m.add_piecewise_formulation(\n",
+    "        (fuel, y_pts, \"<=\"),  # bounded above by the curve\n",
+    "        (power, x_pts),  # equality role (domain-bounded to [0, 30])\n",
+    "        method=method,\n",
+    "    )\n",
+    "    m.add_constraints(power == power_val)\n",
+    "    m.add_objective(-fuel)  # maximise fuel to push against the bound\n",
+    "    m.solve(output_flag=False)\n",
+    "    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n",
+    "\n",
+    "\n",
+    "for method in [\"lp\", \"sos2\", \"incremental\"]:\n",
+    "    fuel_val, vars_, cons_ = solve(method, 15)\n",
+    "    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
+   ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "All three give `fuel=25` at `power=15` (which is `f(15)` exactly) \u2014 the math is equivalent. The LP method is strictly cheaper: no auxiliary variables, just three chord constraints and two domain bounds.\n\nThe SOS2 and incremental methods create lambdas (or deltas + binaries) and split the link into a pinned-equality constraint plus a signed bounded link \u2014 but the feasible region is the same."
+   "source": [
+    "All three give `fuel=25` at `power=15` (which is `f(15)` exactly) \u2014 the math is equivalent. The LP method is strictly cheaper: no auxiliary variables, just three chord constraints and two domain bounds.\n",
+    "\n",
+    "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 \u2014 but the feasible region is the same."
+   ]
   },
   {
    "cell_type": "markdown",
@@ -190,7 +215,16 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Summary\n\n- Default is all-equality: every tuple lies on the curve.\n- Append `\"<=\"` or `\">=\"` as a third tuple element to mark one expression as bounded by the curve.\n- `method=\"auto\"` picks the most efficient formulation: LP for matching-curvature 2-variable inequalities, otherwise SOS2 or incremental.\n- At most one tuple may be bounded; with 3+ tuples all must be equality. Multi-bounded and N\u22653 inequality use cases \u2014 please open an issue at https://github.com/PyPSA/linopy/issues so we can scope them."
+   "source": [
+    "## Summary\n",
+    "\n",
+    "- **One bounded tuple + a 2-variable formulation** gives a hypograph (`<=`) or epigraph (`>=`) feasible region.\n",
+    "- **Curvature \u00d7 sign matching** \u2014 concave + `<=` or convex + `>=` \u2014 lets `method=\"auto\"` skip MIP entirely. Mismatched combinations fall back to SOS2/incremental with a signed link.\n",
+    "- **`method=\"lp\"` is strict** \u2014 it raises on a mismatched curvature rather than silently encoding the wrong region.\n",
+    "- At most one tuple may carry a non-`==` sign, and 3+ tuples must all be `==`. Multi-bounded / N\u22653 inequalities \u2014 open an issue at https://github.com/PyPSA/linopy/issues.\n",
+    "\n",
+    "**See also**: [reference page](piecewise-linear-constraints) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial) for all-equality, unit commitment, CHP, fleets, slopes."
+   ]
   }
  ],
  "metadata": {
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index af79fee1..3a962f44 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -6,16 +6,27 @@
    "source": [
     "# Creating Piecewise Linear Constraints\n",
     "\n",
-    "`add_piecewise_formulation` links variables through shared breakpoint weights. Every section below stacks one feature on top of a small shared dispatch pattern \u2014 if you want the math, see the [reference page](piecewise-linear-constraints). For inequality bounds and the LP chord formulation in depth, see the [inequality bounds notebook](piecewise-inequality-bounds-tutorial).\n",
-    "\n",
-    "The baseline we extend:\n",
+    "`add_piecewise_formulation` links variables through a shared piecewise-linear curve. Pair each variable with its breakpoint values; the solver puts every variable on the *same* point of the curve at every feasible solution.\n",
     "\n",
     "```python\n",
     "m.add_piecewise_formulation(\n",
     "    (power, [0, 30, 60, 100]),\n",
     "    (fuel,  [0, 36, 84, 170]),\n",
     ")\n",
-    "```"
+    "```\n",
+    "\n",
+    "This tutorial walks through the main features of `add_piecewise_formulation`. For the formulation math see the [reference page](piecewise-linear-constraints); for the inequality variant in depth see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial).\n",
+    "\n",
+    "**Roadmap**\n",
+    "\n",
+    "1. Getting started — the basic 2-variable equality.\n",
+    "2. Picking a method — `\"auto\"`, `\"sos2\"`, `\"incremental\"`, `\"lp\"`.\n",
+    "3. Disjunctive segments — disconnected operating regions with `segments()`.\n",
+    "4. Inequality bounds — `<=` / `>=` per-tuple sign.\n",
+    "5. Unit commitment — gating with `active=...`.\n",
+    "6. N-variable linking — CHP plants and beyond.\n",
+    "7. Per-entity breakpoints — fleets with different curves.\n",
+    "8. Specifying with slopes — `linopy.Slopes`."
    ]
   },
   {
@@ -107,13 +118,13 @@
    "source": [
     "## 2. Picking a method\n",
     "\n",
-    "`method=\"auto\"` (default) picks the cheapest correct formulation based on sign, curvature and breakpoint layout. The explicit options \u2014 `\"sos2\"`, `\"incremental\"`, `\"lp\"` \u2014 give the same optimum on equality cases where they all apply, so the choice is about **cost** (auxiliary variables, solver capability), not correctness.\n",
+    "`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.\n",
     "\n",
     "| method | needs | creates |\n",
     "|---|---|---|\n",
     "| `sos2` | SOS2-capable solver | lambdas (continuous) |\n",
     "| `incremental` | MIP solver, strictly monotonic breakpoints | deltas (continuous) + binaries |\n",
-    "| `lp` | any LP solver | no variables \u2014 requires `sign != \"==\"`, 2 tuples, matching curvature |\n",
+    "| `lp` | any LP solver | no variables — requires `sign != \"==\"`, 2 tuples, matching curvature |\n",
     "\n",
     "Below: all applicable methods yield the same fuel dispatch on this convex curve."
    ]
@@ -147,9 +158,11 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 3. Disjunctive segments \u2014 gaps in the operating range\n",
+    "## 3. Disjunctive segments — gaps in the operating range\n",
+    "\n",
+    "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.\n",
     "\n",
-    "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."
+    "Below the first segment is `(0, 0)` — a degenerate \"off\" state where both endpoints sit at the origin, so the unit produces no power and incurs no cost. The second segment is the active range, with cost rising from 125 to 200 over 50–80 MW."
    ]
   },
   {
@@ -190,15 +203,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 4. Inequality bounds \u2014 per-tuple sign\n",
+    "## 4. Inequality bounds — per-tuple sign\n",
     "\n",
-    "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 \u2264 f(x)`) and epigraph (`y \u2265 f(x)`) are the canonical cases.\n",
+    "Append a third tuple element (`\"<=\"` or `\">=\"`) to mark a single expression as **bounded** by the curve instead of entering as an equality. The 2-variable hypograph (`y ≤ f(x)`) and epigraph (`y ≥ f(x)`) are the canonical cases.\n",
     "\n",
-    "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation \u2014 no binaries, no SOS2.\n",
+    "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "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.\n",
-    "\n",
-    "See the [inequality bounds notebook](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
+    "For physical realism, bound a *curtailable* output (heat rejection, emissions after post-treatment, curtailable power) rather than a consumption-side variable like fuel intake — otherwise the formulation admits operating points the plant can't realise. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
    ]
   },
   {
@@ -214,21 +225,21 @@
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=120, coords=[time])\n",
-    "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
+    "heat = m.add_variables(name=\"heat\", lower=0, coords=[time])\n",
     "\n",
-    "# concave curve: diminishing marginal fuel per MW\n",
-    "x_pts = [0, 50, 90, 120]\n",
-    "y_pts = [0, 40, 80, 120]\n",
+    "# Concave heat-extraction curve: diminishing marginal heat per MW of power\n",
+    "power_pts = [0, 50, 90, 120]\n",
+    "heat_pts = [0, 40, 80, 120]\n",
     "pwf = m.add_piecewise_formulation(\n",
-    "    (fuel, x_pts, \"<=\"),  # bounded above by the curve\n",
-    "    (power, y_pts),  # pinned to the curve\n",
+    "    (heat, heat_pts, \"<=\"),  # heat ≤ f(power) — curtailable output\n",
+    "    (power, power_pts),  # power ∈ [0, 120]\n",
     ")\n",
     "m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))\n",
-    "m.add_objective(-fuel.sum())  # push fuel against the bound\n",
+    "m.add_objective(-heat.sum())  # push heat against the bound\n",
     "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
-    "m.solution[[\"power\", \"fuel\"]].to_pandas()"
+    "m.solution[[\"power\", \"heat\"]].to_pandas()"
    ]
   },
   {
@@ -237,15 +248,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# x_pts are fuel breakpoints, y_pts are power breakpoints \u2014 swap so power is on the x-axis\n",
-    "plot_curve(y_pts, x_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
+    "plot_curve(\n",
+    "    power_pts,\n",
+    "    heat_pts,\n",
+    "    m.solution[\"power\"].values,\n",
+    "    m.solution[\"heat\"].values,\n",
+    "    ylabel=\"heat\",\n",
+    ");"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 5. Unit commitment \u2014 `active`\n",
+    "## 5. Unit commitment — `active`\n",
     "\n",
     "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:\n",
     "\n",
@@ -279,7 +295,7 @@
     "    (fuel, y_pts),\n",
     "    active=commit,\n",
     ")\n",
-    "# demand below p_min at t=1 \u2014 commit must be 0 and backup covers it\n",
+    "# demand below p_min at t=1 — commit must be 0 and backup covers it\n",
     "m.add_constraints(power + backup == xr.DataArray([15, 80, 40], coords=[time]))\n",
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
     "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
@@ -299,9 +315,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 6. N-variable linking \u2014 CHP plant\n",
+    "## 6. N-variable linking — CHP plant\n",
     "\n",
-    "More than two variables can share the same interpolation \u2014 useful for combined heat-and-power plants where power, fuel and heat are all functions of a single operating point."
+    "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."
    ]
   },
   {
@@ -358,7 +374,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 7. Per-entity breakpoints \u2014 a fleet of generators\n",
+    "## 7. Per-entity breakpoints — a fleet of generators\n",
     "\n",
     "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`)."
    ]
@@ -396,9 +412,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 8. Specifying with slopes \u2014 `Slopes`\n",
+    "## 8. Specifying with slopes — `Slopes`\n",
     "\n",
-    "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 \u2014 no need to repeat it. Same curve as section 1:"
+    "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:"
    ]
   },
   {
@@ -421,6 +437,25 @@
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## When to use what\n",
+    "\n",
+    "| Pattern | API |\n",
+    "|---|---|\n",
+    "| `y` is a function of `x` | `(x, x_pts), (y, y_pts)` — all-equality |\n",
+    "| `y` bounded by `f(x)` on a convex/concave curve | `(y, y_pts, \"<=\"` or `\">=\"), (x, x_pts)` — LP if curvature matches |\n",
+    "| Disconnected operating regions | `linopy.segments(...)` per tuple |\n",
+    "| Unit on/off coupling | `active=binary_var` |\n",
+    "| Multiple synchronized outputs (e.g. CHP) | 3+ tuples, all `\"==\"` |\n",
+    "| Different curves per entity | `linopy.breakpoints({...}, dim=...)` |\n",
+    "| Slopes more natural than absolute y-values | `linopy.Slopes(...)` |\n",
+    "\n",
+    "For the formulation math, see the [reference page](piecewise-linear-constraints). For inequality bounds in depth, see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial)."
+   ]
   }
  ],
  "metadata": {

From c5055abe31a75a45f3bf78aa1c23ab71c008a084 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:34:08 +0200
Subject: [PATCH 05/29] docs(piecewise): make solver_name='highs' explicit in
 tutorials

Other linopy tutorials (create-a-model.ipynb, manipulating-models.ipynb,
etc.) explicitly pass solver_name='highs' to m.solve(...). The piecewise
tutorials relied on the default-available-solver fallback, which made
the HiGHS-specific output_flag kwarg look mysterious.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb  | 48 ++++++++++-----------
 examples/piecewise-linear-constraints.ipynb | 16 +++----
 2 files changed, 32 insertions(+), 32 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 682a3656..9ae689be 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "# Creating Piecewise Inequality Bounds\n",
     "\n",
-    "When you only need a one-sided bound from a piecewise curve \u2014 `y \u2264 f(x)` for a concave upper envelope, `y \u2265 f(x)` for a convex lower envelope \u2014 `add_piecewise_formulation` accepts an optional sign as the third tuple element:\n",
+    "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:\n",
     "\n",
     "```python\n",
     "m.add_piecewise_formulation(\n",
@@ -15,7 +15,7 @@
     ")\n",
     "```\n",
     "\n",
-    "The pay-off is a pure-LP encoding when the curve's curvature matches the sign \u2014 no SOS2, no binaries. This notebook covers the geometry of the feasible region, the curvature \u00d7 sign combinations that unlock the LP path, and what happens when they don't match.\n",
+    "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.\n",
     "\n",
     "For the formulation math see the [reference page](piecewise-linear-constraints); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial).\n",
     "\n",
@@ -23,11 +23,11 @@
     "\n",
     "| Tuple form | Role | What it constrains |\n",
     "|---|---|---|\n",
-    "| `(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]` \u2014 one coordinate alone can't locate a curve point. |\n",
-    "| `(expr, breaks, \"<=\")` | bounded above | `expr \u2264 f(other tuples)`. |\n",
-    "| `(expr, breaks, \">=\")` | bounded below | `expr \u2265 f(other tuples)`. |\n",
+    "| `(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. |\n",
+    "| `(expr, breaks, \"<=\")` | bounded above | `expr ≤ f(other tuples)`. |\n",
+    "| `(expr, breaks, \">=\")` | bounded below | `expr ≥ f(other tuples)`. |\n",
     "\n",
-    "Currently at most one tuple may carry a non-equality sign, and 3+ tuples must all be equality. Multi-bounded and N\u22653 inequality cases aren't supported yet \u2014 if you have a concrete use case, please open an issue at https://github.com/PyPSA/linopy/issues so we can scope it properly."
+    "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 https://github.com/PyPSA/linopy/issues so we can scope it properly."
    ]
   },
   {
@@ -50,7 +50,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Setup \u2014 a concave curve\n\nWe use a concave, monotonically increasing curve. With a tuple bounded `<=`, the LP method is applicable (concave + `<=` is a tight relaxation)."
+   "source": "## Setup — a concave curve\n\nWe use a concave, monotonically increasing curve. With a tuple bounded `<=`, the LP method is applicable (concave + `<=` is a tight relaxation)."
   },
   {
    "cell_type": "code",
@@ -81,11 +81,11 @@
     "\n",
     "With one tuple bounded `<=` and our concave curve, the three methods give the **same** feasible region within `[x_0, x_n]`:\n",
     "\n",
-    "- **`method=\"lp\"`** \u2014 tangent lines + domain bounds. No auxiliary variables.\n",
-    "- **`method=\"sos2\"`** \u2014 lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
-    "- **`method=\"incremental\"`** \u2014 delta fractions + binaries + split link. Same mathematics, MIP encoding instead of SOS2.\n",
+    "- **`method=\"lp\"`** — tangent lines + domain bounds. No auxiliary variables.\n",
+    "- **`method=\"sos2\"`** — lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
+    "- **`method=\"incremental\"`** — delta fractions + binaries + split link. Same mathematics, MIP encoding instead of SOS2.\n",
     "\n",
-    "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable \u2014 it's always preferable because it's pure LP.\n",
+    "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable — it's always preferable because it's pure LP.\n",
     "\n",
     "Let's verify they produce the same solution at `power=15`."
    ]
@@ -112,7 +112,7 @@
     "    )\n",
     "    m.add_constraints(power == power_val)\n",
     "    m.add_objective(-fuel)  # maximise fuel to push against the bound\n",
-    "    m.solve(output_flag=False)\n",
+    "    m.solve(solver_name=\"highs\", output_flag=False)\n",
     "    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n",
     "\n",
     "\n",
@@ -125,15 +125,15 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "All three give `fuel=25` at `power=15` (which is `f(15)` exactly) \u2014 the math is equivalent. The LP method is strictly cheaper: no auxiliary variables, just three chord constraints and two domain bounds.\n",
+    "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.\n",
     "\n",
-    "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 \u2014 but the feasible region is the same."
+    "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."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Visualising the feasible region\n\nThe feasible region for `(power, fuel)` with `fuel` bounded `<=` is the **hypograph** of `f` restricted to the curve's x-domain:\n\n$$\\{ (x, y) : x_0 \\le x \\le x_n,\\ y \\le f(x) \\}$$\n\nWe colour green feasible points, red infeasible ones. Three test points:\n\n- `(15, 15)` \u2014 inside the curve, `15 \u2264 f(15)=25` \u2713\n- `(15, 25)` \u2014 on the curve \u2713\n- `(15, 29)` \u2014 above `f(15)`, should be infeasible \u2717\n- `(35, 20)` \u2014 power beyond domain, infeasible \u2717"
+   "source": "## Visualising the feasible region\n\nThe feasible region for `(power, fuel)` with `fuel` bounded `<=` is the **hypograph** of `f` restricted to the curve's x-domain:\n\n$$\\{ (x, y) : x_0 \\le x \\le x_n,\\ y \\le f(x) \\}$$\n\nWe colour green feasible points, red infeasible ones. Three test points:\n\n- `(15, 15)` — inside the curve, `15 ≤ f(15)=25` ✓\n- `(15, 25)` — on the curve ✓\n- `(15, 29)` — above `f(15)`, should be infeasible ✗\n- `(35, 20)` — power beyond domain, infeasible ✗"
   },
   {
    "cell_type": "code",
@@ -169,7 +169,7 @@
     "ax.set(\n",
     "    xlabel=\"power\",\n",
     "    ylabel=\"fuel\",\n",
-    "    title=\"sign='<=' feasible region \u2014 hypograph of f(x) on [x_0, x_n]\",\n",
+    "    title=\"sign='<=' feasible region — hypograph of f(x) on [x_0, x_n]\",\n",
     ")\n",
     "ax.grid(alpha=0.3)\n",
     "ax.legend()\n",
@@ -182,16 +182,16 @@
    "source": [
     "## When is LP the right choice?\n",
     "\n",
-    "`tangent_lines` imposes the **intersection** of chord inequalities. Whether that intersection matches the true hypograph/epigraph of `f` depends on the curvature \u00d7 sign combination:\n",
+    "`tangent_lines` imposes the **intersection** of chord inequalities. Whether that intersection matches the true hypograph/epigraph of `f` depends on the curvature × sign combination:\n",
     "\n",
     "| curvature | bounded `<=` | bounded `>=` |\n",
     "|-----------|--------------|--------------|\n",
-    "| **concave** | **hypograph (exact \u2713)** | **wrong region** \u2014 requires `y \u2265 max_k chord_k(x) > f(x)` |\n",
-    "| **convex**  | **wrong region** \u2014 requires `y \u2264 min_k chord_k(x) < f(x)` | **epigraph (exact \u2713)** |\n",
+    "| **concave** | **hypograph (exact ✓)** | **wrong region** — requires `y ≥ max_k chord_k(x) > f(x)` |\n",
+    "| **convex**  | **wrong region** — requires `y ≤ min_k chord_k(x) < f(x)` | **epigraph (exact ✓)** |\n",
     "| linear      | exact | exact |\n",
     "| mixed (non-convex) | convex hull of `f` (wrong for exact hypograph) | concave hull of `f` (wrong for exact epigraph) |\n",
     "\n",
-    "In the \u2717 cases, tangent lines do **not** give a loose relaxation \u2014 they give a **strictly wrong feasible region** that rejects points satisfying the true constraint. Example: for a concave `f` with `y \u2265 f(x)`, the chord of any piece extrapolated over another piece's x-range lies *above* `f`, so `y \u2265 max_k chord_k(x)` forbids `y = f(x)` itself.\n",
+    "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.\n",
     "\n",
     "`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.\n",
     "\n",
@@ -210,7 +210,7 @@
     }
    },
    "outputs": [],
-   "source": "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\nx_nc = [0, 10, 20, 30]\ny_nc = [0, 20, 10, 30]  # slopes change sign \u2192 mixed convexity\n\nm1 = linopy.Model()\nx1 = m1.add_variables(lower=0, upper=30, name=\"x\")\ny1 = m1.add_variables(lower=0, upper=40, name=\"y\")\nf1 = m1.add_piecewise_formulation((y1, y_nc, \"<=\"), (x1, x_nc))\nprint(f\"non-convex + '<=' \u2192 {f1.method}\")\n\n# 2. Concave curve + sign='>=': LP would be loose \u2192 auto falls back to MIP\nx_cc = [0, 10, 20, 30]\ny_cc = [0, 20, 30, 35]  # concave\n\nm2 = linopy.Model()\nx2 = m2.add_variables(lower=0, upper=30, name=\"x\")\ny2 = m2.add_variables(lower=0, upper=40, name=\"y\")\nf2 = m2.add_piecewise_formulation((y2, y_cc, \">=\"), (x2, x_cc))\nprint(f\"concave + '>='    \u2192 {f2.method}\")\n\n# 3. Explicit method=\"lp\" with mismatched curvature raises\ntry:\n    m3 = linopy.Model()\n    x3 = m3.add_variables(lower=0, upper=30, name=\"x\")\n    y3 = m3.add_variables(lower=0, upper=40, name=\"y\")\n    m3.add_piecewise_formulation((y3, y_cc, \">=\"), (x3, x_cc), method=\"lp\")\nexcept ValueError as e:\n    print(f\"lp(concave, '>=')  \u2192 raises: {e}\")"
+   "source": "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\nx_nc = [0, 10, 20, 30]\ny_nc = [0, 20, 10, 30]  # slopes change sign → mixed convexity\n\nm1 = linopy.Model()\nx1 = m1.add_variables(lower=0, upper=30, name=\"x\")\ny1 = m1.add_variables(lower=0, upper=40, name=\"y\")\nf1 = m1.add_piecewise_formulation((y1, y_nc, \"<=\"), (x1, x_nc))\nprint(f\"non-convex + '<=' → {f1.method}\")\n\n# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP\nx_cc = [0, 10, 20, 30]\ny_cc = [0, 20, 30, 35]  # concave\n\nm2 = linopy.Model()\nx2 = m2.add_variables(lower=0, upper=30, name=\"x\")\ny2 = m2.add_variables(lower=0, upper=40, name=\"y\")\nf2 = m2.add_piecewise_formulation((y2, y_cc, \">=\"), (x2, x_cc))\nprint(f\"concave + '>='    → {f2.method}\")\n\n# 3. Explicit method=\"lp\" with mismatched curvature raises\ntry:\n    m3 = linopy.Model()\n    x3 = m3.add_variables(lower=0, upper=30, name=\"x\")\n    y3 = m3.add_variables(lower=0, upper=40, name=\"y\")\n    m3.add_piecewise_formulation((y3, y_cc, \">=\"), (x3, x_cc), method=\"lp\")\nexcept ValueError as e:\n    print(f\"lp(concave, '>=')  → raises: {e}\")"
   },
   {
    "cell_type": "markdown",
@@ -219,9 +219,9 @@
     "## Summary\n",
     "\n",
     "- **One bounded tuple + a 2-variable formulation** gives a hypograph (`<=`) or epigraph (`>=`) feasible region.\n",
-    "- **Curvature \u00d7 sign matching** \u2014 concave + `<=` or convex + `>=` \u2014 lets `method=\"auto\"` skip MIP entirely. Mismatched combinations fall back to SOS2/incremental with a signed link.\n",
-    "- **`method=\"lp\"` is strict** \u2014 it raises on a mismatched curvature rather than silently encoding the wrong region.\n",
-    "- At most one tuple may carry a non-`==` sign, and 3+ tuples must all be `==`. Multi-bounded / N\u22653 inequalities \u2014 open an issue at https://github.com/PyPSA/linopy/issues.\n",
+    "- **Curvature × sign matching** — concave + `<=` or convex + `>=` — lets `method=\"auto\"` skip MIP entirely. Mismatched combinations fall back to SOS2/incremental with a signed link.\n",
+    "- **`method=\"lp\"` is strict** — it raises on a mismatched curvature rather than silently encoding the wrong region.\n",
+    "- At most one tuple may carry a non-`==` sign, and 3+ tuples must all be `==`. Multi-bounded / N≥3 inequalities — open an issue at https://github.com/PyPSA/linopy/issues.\n",
     "\n",
     "**See also**: [reference page](piecewise-linear-constraints) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial) for all-equality, unit commitment, CHP, fleets, slopes."
    ]
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 3a962f44..8ef6ce9f 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -92,7 +92,7 @@
     "pwf = m.add_piecewise_formulation((power, x_pts), (fuel, y_pts))\n",
     "m.add_constraints(power == demand, name=\"demand\")\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
@@ -147,7 +147,7 @@
     "    m.add_piecewise_formulation((power, x_pts), (fuel, y_pts), method=method)\n",
     "    m.add_constraints(power == demand, name=\"demand\")\n",
     "    m.add_objective(fuel.sum())\n",
-    "    m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "    m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "    return m.solution[\"fuel\"].to_pandas()\n",
     "\n",
     "\n",
@@ -188,7 +188,7 @@
     ")\n",
     "m.add_constraints(power + backup == xr.DataArray([15, 60, 75], coords=[time]))\n",
     "m.add_objective(cost.sum() + 10 * backup.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"cost\", \"backup\"]].to_pandas()"
    ]
   },
@@ -236,7 +236,7 @@
     ")\n",
     "m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))\n",
     "m.add_objective(-heat.sum())  # push heat against the bound\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"heat\"]].to_pandas()"
@@ -298,7 +298,7 @@
     "# demand below p_min at t=1 — commit must be 0 and backup covers it\n",
     "m.add_constraints(power + backup == xr.DataArray([15, 80, 40], coords=[time]))\n",
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"commit\", \"power\", \"fuel\", \"backup\"]].to_pandas()"
    ]
   },
@@ -346,7 +346,7 @@
     ")\n",
     "m.add_constraints(fuel == xr.DataArray([20, 100, 160], coords=[time]))\n",
     "m.add_objective(power.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\", \"heat\"]].to_pandas().round(2)"
    ]
   },
@@ -404,7 +404,7 @@
     "m.add_piecewise_formulation((power, x_gen), (fuel, y_gen))\n",
     "m.add_constraints(power.sum(\"gen\") == xr.DataArray([80, 120, 50], coords=[time]))\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\"]].to_dataframe()"
    ]
   },
@@ -433,7 +433,7 @@
     ")\n",
     "m.add_constraints(power == demand, name=\"demand\")\n",
     "m.add_objective(fuel.sum())\n",
-    "m.solve(reformulate_sos=\"auto\", output_flag=False)\n",
+    "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
    ]

From 29cc65b71f60d3951488b65611c3824925fcc0d2 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:37:56 +0200
Subject: [PATCH 06/29] =?UTF-8?q?docs(piecewise):=20reframe=20inequality-b?=
 =?UTF-8?q?ounds=20tutorial=20as=20power=20=E2=89=A4=20f(fuel)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The previous framing — fuel <= f(power) with objective -fuel — is
logically off: bounding fuel from above lets the solver choose any
non-negative fuel, including zero (since fuel is consumption-side).
The artificial maximise-fuel objective hid this but the underlying
relation didn't make physical sense.

Power ≤ f(fuel) is the classic production-function bound: f maps fuel
input to the maximum power the unit can deliver, and the unit can
always run below that (output is curtailable).  Maximising power
against the bound is now a sensible LP objective rather than a
contrivance.

- Reframed cells 0, 2, 3, 4, 5, 6, 7, 8, 10 around fuel_pts /
  power_pts; axis labels and the hypograph plot now read (fuel, power)
  with f(fuel) as the production curve.
- Cell 5's objective is now `-power` (maximise power up to the curve
  bound), matching real LPs.
- Test points and curvature × sign analysis unchanged — only the
  physical interpretation flipped.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb | 114 +++++++++++++++------
 1 file changed, 81 insertions(+), 33 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 9ae689be..1e4d04a3 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -10,8 +10,8 @@
     "\n",
     "```python\n",
     "m.add_piecewise_formulation(\n",
-    "    (fuel,  y_pts, \"<=\"),   # bounded above by the curve\n",
-    "    (power, x_pts),         # equality role\n",
+    "    (power, power_pts, \"<=\"),  # power ≤ f(fuel) — output curtailable\n",
+    "    (fuel,  fuel_pts),         # equality role\n",
     ")\n",
     "```\n",
     "\n",
@@ -50,7 +50,11 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Setup — a concave curve\n\nWe use a concave, monotonically increasing curve. With a tuple bounded `<=`, the LP method is applicable (concave + `<=` is a tight relaxation)."
+   "source": [
+    "## Setup — a concave production curve\n",
+    "\n",
+    "A concave, monotonically increasing curve maps fuel input to the maximum power the plant can deliver.  Bounding `power` by this curve with `<=` says the unit may run *at or below* the maximum — power output is curtailable for a given fuel draw.  Concave + `<=` is exactly the combination that lets the LP method apply."
+   ]
   },
   {
    "cell_type": "code",
@@ -63,12 +67,12 @@
    },
    "outputs": [],
    "source": [
-    "x_pts = np.array([0.0, 10.0, 20.0, 30.0])\n",
-    "y_pts = np.array([0.0, 20.0, 30.0, 35.0])  # slopes 2, 1, 0.5 (concave)\n",
+    "fuel_pts = np.array([0.0, 10.0, 20.0, 30.0])\n",
+    "power_pts = np.array([0.0, 20.0, 30.0, 35.0])  # slopes 2, 1, 0.5 (concave)\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(5, 4))\n",
-    "ax.plot(x_pts, y_pts, \"o-\", color=\"C0\", lw=2)\n",
-    "ax.set(xlabel=\"power\", ylabel=\"fuel\", title=\"Concave reference curve f(x)\")\n",
+    "ax.plot(fuel_pts, power_pts, \"o-\", color=\"C0\", lw=2)\n",
+    "ax.set(xlabel=\"fuel\", ylabel=\"power\", title=\"Concave production curve f(fuel)\")\n",
     "ax.grid(alpha=0.3)\n",
     "plt.tight_layout()"
    ]
@@ -79,7 +83,7 @@
    "source": [
     "## Three methods, identical feasible region\n",
     "\n",
-    "With one tuple bounded `<=` and our concave curve, the three methods give the **same** feasible region within `[x_0, x_n]`:\n",
+    "With `power` bounded `<=` and our concave curve, the three methods give the **same** feasible region for `fuel ∈ [0, 30]`:\n",
     "\n",
     "- **`method=\"lp\"`** — tangent lines + domain bounds. No auxiliary variables.\n",
     "- **`method=\"sos2\"`** — lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
@@ -87,7 +91,7 @@
     "\n",
     "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable — it's always preferable because it's pure LP.\n",
     "\n",
-    "Let's verify they produce the same solution at `power=15`."
+    "Let's verify they produce the same solution at `fuel=15`."
    ]
   },
   {
@@ -101,31 +105,31 @@
    },
    "outputs": [],
    "source": [
-    "def solve(method, power_val):\n",
+    "def solve(method, fuel_val):\n",
     "    m = linopy.Model()\n",
-    "    power = m.add_variables(lower=0, upper=30, name=\"power\")\n",
-    "    fuel = m.add_variables(lower=0, upper=40, name=\"fuel\")\n",
+    "    fuel = m.add_variables(lower=0, upper=30, name=\"fuel\")\n",
+    "    power = m.add_variables(lower=0, upper=40, name=\"power\")\n",
     "    m.add_piecewise_formulation(\n",
-    "        (fuel, y_pts, \"<=\"),  # bounded above by the curve\n",
-    "        (power, x_pts),  # equality role (domain-bounded to [0, 30])\n",
+    "        (power, power_pts, \"<=\"),  # power ≤ f(fuel) — curtailable output\n",
+    "        (fuel, fuel_pts),  # equality role (domain-bounded to [0, 30])\n",
     "        method=method,\n",
     "    )\n",
-    "    m.add_constraints(power == power_val)\n",
-    "    m.add_objective(-fuel)  # maximise fuel to push against the bound\n",
+    "    m.add_constraints(fuel == fuel_val)\n",
+    "    m.add_objective(-power)  # maximise power output against the upper bound\n",
     "    m.solve(solver_name=\"highs\", output_flag=False)\n",
-    "    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n",
+    "    return float(m.solution[\"power\"]), list(m.variables), list(m.constraints)\n",
     "\n",
     "\n",
     "for method in [\"lp\", \"sos2\", \"incremental\"]:\n",
-    "    fuel_val, vars_, cons_ = solve(method, 15)\n",
-    "    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
+    "    power_val, vars_, cons_ = solve(method, 15)\n",
+    "    print(f\"{method:12}: power={power_val:.2f}  vars={vars_}  cons={cons_}\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "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.\n",
+    "All three give `power=25` at `fuel=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.\n",
     "\n",
     "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."
    ]
@@ -133,7 +137,20 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": "## Visualising the feasible region\n\nThe feasible region for `(power, fuel)` with `fuel` bounded `<=` is the **hypograph** of `f` restricted to the curve's x-domain:\n\n$$\\{ (x, y) : x_0 \\le x \\le x_n,\\ y \\le f(x) \\}$$\n\nWe colour green feasible points, red infeasible ones. Three test points:\n\n- `(15, 15)` — inside the curve, `15 ≤ f(15)=25` ✓\n- `(15, 25)` — on the curve ✓\n- `(15, 29)` — above `f(15)`, should be infeasible ✗\n- `(35, 20)` — power beyond domain, infeasible ✗"
+   "source": [
+    "## Visualising the feasible region\n",
+    "\n",
+    "The feasible region for `(fuel, power)` with `power` bounded `<=` is the **hypograph** of `f` restricted to the fuel domain:\n",
+    "\n",
+    "$$\\{ (\\mathrm{fuel}, \\mathrm{power}) : 0 \\le \\mathrm{fuel} \\le 30,\\ \\mathrm{power} \\le f(\\mathrm{fuel}) \\}$$\n",
+    "\n",
+    "Below we colour feasible points green, infeasible ones red:\n",
+    "\n",
+    "- `(15, 15)` — under the curve, `15 ≤ f(15)=25` ✓\n",
+    "- `(15, 25)` — on the curve ✓\n",
+    "- `(15, 29)` — above `f(15)`, infeasible ✗\n",
+    "- `(35, 20)` — fuel beyond domain, infeasible ✗"
+   ]
   },
   {
    "cell_type": "code",
@@ -146,10 +163,10 @@
    },
    "outputs": [],
    "source": [
-    "def in_hypograph(px, py):\n",
-    "    if px < x_pts[0] or px > x_pts[-1]:\n",
+    "def in_hypograph(fx, py):\n",
+    "    if fx < fuel_pts[0] or fx > fuel_pts[-1]:\n",
     "        return False\n",
-    "    return py <= np.interp(px, x_pts, y_pts)\n",
+    "    return py <= np.interp(fx, fuel_pts, power_pts)\n",
     "\n",
     "\n",
     "xx, yy = np.meshgrid(np.linspace(-2, 38, 200), np.linspace(-5, 45, 200))\n",
@@ -159,17 +176,17 @@
     "\n",
     "fig, ax = plt.subplots(figsize=(6, 5))\n",
     "ax.contourf(xx, yy, region, levels=[0.5, 1], colors=[\"lightsteelblue\"], alpha=0.5)\n",
-    "ax.plot(x_pts, y_pts, \"o-\", color=\"C0\", lw=2, label=\"f(x)\")\n",
-    "for px, py in test_points:\n",
-    "    feas = in_hypograph(px, py)\n",
+    "ax.plot(fuel_pts, power_pts, \"o-\", color=\"C0\", lw=2, label=\"f(fuel)\")\n",
+    "for fx, py in test_points:\n",
+    "    feas = in_hypograph(fx, py)\n",
     "    ax.scatter(\n",
-    "        [px], [py], color=\"green\" if feas else \"red\", zorder=5, s=80, edgecolors=\"black\"\n",
+    "        [fx], [py], color=\"green\" if feas else \"red\", zorder=5, s=80, edgecolors=\"black\"\n",
     "    )\n",
-    "    ax.annotate(f\"({px}, {py})\", (px, py), textcoords=\"offset points\", xytext=(8, 5))\n",
+    "    ax.annotate(f\"({fx}, {py})\", (fx, py), textcoords=\"offset points\", xytext=(8, 5))\n",
     "ax.set(\n",
-    "    xlabel=\"power\",\n",
-    "    ylabel=\"fuel\",\n",
-    "    title=\"sign='<=' feasible region — hypograph of f(x) on [x_0, x_n]\",\n",
+    "    xlabel=\"fuel\",\n",
+    "    ylabel=\"power\",\n",
+    "    title=\"sign='<=' feasible region — hypograph of f(fuel) on [0, 30]\",\n",
     ")\n",
     "ax.grid(alpha=0.3)\n",
     "ax.legend()\n",
@@ -210,7 +227,38 @@
     }
    },
    "outputs": [],
-   "source": "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\nx_nc = [0, 10, 20, 30]\ny_nc = [0, 20, 10, 30]  # slopes change sign → mixed convexity\n\nm1 = linopy.Model()\nx1 = m1.add_variables(lower=0, upper=30, name=\"x\")\ny1 = m1.add_variables(lower=0, upper=40, name=\"y\")\nf1 = m1.add_piecewise_formulation((y1, y_nc, \"<=\"), (x1, x_nc))\nprint(f\"non-convex + '<=' → {f1.method}\")\n\n# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP\nx_cc = [0, 10, 20, 30]\ny_cc = [0, 20, 30, 35]  # concave\n\nm2 = linopy.Model()\nx2 = m2.add_variables(lower=0, upper=30, name=\"x\")\ny2 = m2.add_variables(lower=0, upper=40, name=\"y\")\nf2 = m2.add_piecewise_formulation((y2, y_cc, \">=\"), (x2, x_cc))\nprint(f\"concave + '>='    → {f2.method}\")\n\n# 3. Explicit method=\"lp\" with mismatched curvature raises\ntry:\n    m3 = linopy.Model()\n    x3 = m3.add_variables(lower=0, upper=30, name=\"x\")\n    y3 = m3.add_variables(lower=0, upper=40, name=\"y\")\n    m3.add_piecewise_formulation((y3, y_cc, \">=\"), (x3, x_cc), method=\"lp\")\nexcept ValueError as e:\n    print(f\"lp(concave, '>=')  → raises: {e}\")"
+   "source": [
+    "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\n",
+    "fuel_nc = [0, 10, 20, 30]\n",
+    "power_nc = [0, 20, 10, 30]  # slopes change sign → mixed convexity\n",
+    "\n",
+    "m1 = linopy.Model()\n",
+    "fuel1 = m1.add_variables(lower=0, upper=30, name=\"fuel\")\n",
+    "power1 = m1.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "f1 = m1.add_piecewise_formulation((power1, power_nc, \"<=\"), (fuel1, fuel_nc))\n",
+    "print(f\"non-convex + '<=' → {f1.method}\")\n",
+    "\n",
+    "# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP\n",
+    "fuel_cc = [0, 10, 20, 30]\n",
+    "power_cc = [0, 20, 30, 35]  # concave\n",
+    "\n",
+    "m2 = linopy.Model()\n",
+    "fuel2 = m2.add_variables(lower=0, upper=30, name=\"fuel\")\n",
+    "power2 = m2.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "f2 = m2.add_piecewise_formulation((power2, power_cc, \">=\"), (fuel2, fuel_cc))\n",
+    "print(f\"concave + '>='    → {f2.method}\")\n",
+    "\n",
+    "# 3. Explicit method=\"lp\" with mismatched curvature raises\n",
+    "try:\n",
+    "    m3 = linopy.Model()\n",
+    "    fuel3 = m3.add_variables(lower=0, upper=30, name=\"fuel\")\n",
+    "    power3 = m3.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "    m3.add_piecewise_formulation(\n",
+    "        (power3, power_cc, \">=\"), (fuel3, fuel_cc), method=\"lp\"\n",
+    "    )\n",
+    "except ValueError as e:\n",
+    "    print(f\"lp(concave, '>=')  → raises: {e}\")"
+   ]
   },
   {
    "cell_type": "markdown",

From 9829e8e89f1155218b2042adab1fb3b046b3b06c Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:40:26 +0200
Subject: [PATCH 07/29] =?UTF-8?q?docs(piecewise):=20reframe=20tutorial=201?=
 =?UTF-8?q?=20=C2=A74=20as=20fuel=20=E2=89=A5=20f(power)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Previous version bounded `heat` by a curve labelled "concave" — but the
slopes (0.8, 1.0, 1.33) are increasing, so the curve is actually convex,
and `heat <= f(power)` on a convex curve doesn't dispatch to LP.

Switched to a convex heat-rate curve with `fuel >= f(power)`: over-
fuelling is physically admissible (waste heat) but wasteful, so
minimising fuel pulls the operating point onto the curve.  This is the
epigraph half of the LP-applicable region, complementing the inequality-
bounds tutorial's hypograph example (concave + `<=`).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 28 ++++++++-------------
 1 file changed, 11 insertions(+), 17 deletions(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 8ef6ce9f..f8c31f65 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -209,7 +209,7 @@
     "\n",
     "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "For physical realism, bound a *curtailable* output (heat rejection, emissions after post-treatment, curtailable power) rather than a consumption-side variable like fuel intake — otherwise the formulation admits operating points the plant can't realise. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
+    "Below: a convex heat-rate curve, `fuel ≥ f(power)`. Over-fuelling is admissible but wasteful, so minimising fuel pulls the operating point down onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for the dual case (concave + `<=`), mismatched curvature, and more geometry."
    ]
   },
   {
@@ -224,22 +224,22 @@
    "outputs": [],
    "source": [
     "m = linopy.Model()\n",
-    "power = m.add_variables(name=\"power\", lower=0, upper=120, coords=[time])\n",
-    "heat = m.add_variables(name=\"heat\", lower=0, coords=[time])\n",
+    "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
+    "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
     "\n",
-    "# Concave heat-extraction curve: diminishing marginal heat per MW of power\n",
-    "power_pts = [0, 50, 90, 120]\n",
-    "heat_pts = [0, 40, 80, 120]\n",
+    "# Convex heat-rate curve: fuel rises faster than linearly with power\n",
+    "power_pts = [0, 30, 60, 100]\n",
+    "fuel_pts = [0, 36, 84, 170]  # slopes 1.2, 1.6, 2.15 (convex)\n",
     "pwf = m.add_piecewise_formulation(\n",
-    "    (heat, heat_pts, \"<=\"),  # heat ≤ f(power) — curtailable output\n",
-    "    (power, power_pts),  # power ∈ [0, 120]\n",
+    "    (fuel, fuel_pts, \">=\"),  # fuel ≥ f(power) — over-fuelling allowed\n",
+    "    (power, power_pts),  # equality role\n",
     ")\n",
     "m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))\n",
-    "m.add_objective(-heat.sum())  # push heat against the bound\n",
+    "m.add_objective(fuel.sum())  # minimise fuel against the lower bound\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
-    "m.solution[[\"power\", \"heat\"]].to_pandas()"
+    "m.solution[[\"power\", \"fuel\"]].to_pandas()"
    ]
   },
   {
@@ -248,13 +248,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_curve(\n",
-    "    power_pts,\n",
-    "    heat_pts,\n",
-    "    m.solution[\"power\"].values,\n",
-    "    m.solution[\"heat\"].values,\n",
-    "    ylabel=\"heat\",\n",
-    ");"
+    "plot_curve(power_pts, fuel_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
    ]
   },
   {

From d7592f3a77b70caa2d8f7622cdfc20f72c80e075 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:42:34 +0200
Subject: [PATCH 08/29] =?UTF-8?q?docs(piecewise):=20use=20fuel=20=E2=89=A5?=
 =?UTF-8?q?=20f(power)=20in=20the=20rst=20Quick=20Start?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Matches the tutorial 1 §4 reframe: convex heat-rate curve with
``fuel >= f(power)``, where over-fuelling is admissible but wasteful so
minimisation pulls the operating point onto the curve.  Replaces the
heat <= f(power) example (which used a curve that was actually convex,
not concave, so it wouldn't have dispatched to LP).

Generalised the "Choice of bounded tuple" guidance to cover both signs:

- ``"<="`` for a controllable dissipation path (curtailment, post-
  treatment).
- ``">="`` for an input whose over-supply is admissible but wasteful
  (fuel, raw materials).

The wrong-direction warning now spells out both anti-patterns: ``"<="``
on fuel, ``">="`` on a non-curtailable output.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 40 ++++++++++++++++------------
 1 file changed, 23 insertions(+), 17 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index c96ad54a..09afcff1 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -49,14 +49,14 @@ Quick Start
 
 .. code-block:: python
 
-    # heat <= f(power).  "auto" picks the cheapest correct formulation:
-    # pure LP (chord constraints) when curvature matches the sign,
-    # SOS2/incremental otherwise.  Bound a curtailable output so
-    # undershooting the curve is physically realisable — see *Choice of
-    # bounded tuple* below.
+    # fuel >= f(power) on a convex heat-rate curve.  Over-fuelling is
+    # admissible but wasteful, so minimising fuel pulls the operating
+    # point onto the curve.  "auto" picks the cheapest correct
+    # formulation: pure LP (chord constraints) when curvature matches
+    # the sign — here convex + ">="  — and SOS2/incremental otherwise.
     m.add_piecewise_formulation(
-        (heat, [0, 20, 30, 35], "<="),
-        (power, [0, 10, 20, 30]),
+        (fuel, [0, 36, 84, 170], ">="),
+        (power, [0, 30, 60, 100]),
     )
 
 Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with
@@ -196,16 +196,22 @@ sign + curvature (convex + ``"<="``, or concave + ``">="``) describes a
 *non-convex* region — ``method="auto"`` falls back to SOS2/incremental
 and ``method="lp"`` raises.
 
-**Choice of bounded tuple.**  The bounded 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.  Marking a consumption-side variable such as fuel
-intake as bounded yields a valid but **loose** formulation: the
-characteristic curve fixes fuel draw at a given load, so ``"<="`` 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.
+**Choice of bounded tuple and sign.**  Pick the sign matching the
+physically admissible direction for that expression:
+
+- ``"<="`` for a quantity with a controllable *dissipation path* — heat
+  rejection via cooling tower (*thermal curtailment*), electrical
+  curtailment, emissions after post-treatment — so undershooting the
+  curve is realisable.
+- ``">="`` for an *input* whose over-supply is admissible but wasteful —
+  fuel, raw materials — so overshooting the curve is realisable
+  (objective pressure then pulls the operating point onto the curve).
+
+The wrong direction (``"<="`` on fuel, ``">="`` on a non-curtailable
+output) yields a valid but **loose** formulation that admits operating
+points the plant cannot physically realise; an objective rewarding the
+wrong direction may then find a non-physical optimum — safe only when
+no such objective pressure exists.
 
 **When is a one-sided bound wanted?**
 

From 31f31f00651c1e125c586d87dfa3f1dfa1918631 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:44:44 +0200
Subject: [PATCH 09/29] =?UTF-8?q?docs(piecewise):=20switch=20Quick=20Start?=
 =?UTF-8?q?=20+=20tutorial=20=C2=A74=20to=20power=20=E2=89=A4=20f(fuel)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Reuses §1's exact breakpoints — only the sign on one tuple changes —
so the reader sees "the same curve, now with an inequality" rather than
a new curve.  Reads as a production function: power output is bounded
by what the fuel input can support, and the unit may run below the
maximum (curtailable output).

Same physical meaning as the previous fuel ≥ f(power) framing (the two
forms describe the same feasible region under inversion), but more
intuitive and continuous with the rest of the tutorial.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        | 15 ++++++++-------
 examples/piecewise-linear-constraints.ipynb | 21 ++++++++++++---------
 2 files changed, 20 insertions(+), 16 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 09afcff1..98187ee0 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -49,14 +49,15 @@ Quick Start
 
 .. code-block:: python
 
-    # fuel >= f(power) on a convex heat-rate curve.  Over-fuelling is
-    # admissible but wasteful, so minimising fuel pulls the operating
-    # point onto the curve.  "auto" picks the cheapest correct
-    # formulation: pure LP (chord constraints) when curvature matches
-    # the sign — here convex + ">="  — and SOS2/incremental otherwise.
+    # power <= f(fuel) on the same curve as above, but read as a
+    # production function: power output is bounded by what the fuel
+    # input can support, and the unit may run below the max
+    # (curtailable output).  "auto" picks the cheapest correct
+    # formulation: pure LP (chord constraints) here, since concave +
+    # "<=" is LP-applicable; SOS2/incremental otherwise.
     m.add_piecewise_formulation(
-        (fuel, [0, 36, 84, 170], ">="),
-        (power, [0, 30, 60, 100]),
+        (power, [0, 30, 60, 100], "<="),
+        (fuel, [0, 36, 84, 170]),
     )
 
 Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index f8c31f65..d54bd2c1 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -209,7 +209,7 @@
     "\n",
     "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "Below: a convex heat-rate curve, `fuel ≥ f(power)`. Over-fuelling is admissible but wasteful, so minimising fuel pulls the operating point down onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for the dual case (concave + `<=`), mismatched curvature, and more geometry."
+    "Below: the same curve as §1, but read as a production function `power ≤ f(fuel)` — the unit *can* run below the maximum power for a given fuel draw (output curtailable). Minimising fuel against a power demand pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
    ]
   },
   {
@@ -227,15 +227,13 @@
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
     "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
     "\n",
-    "# Convex heat-rate curve: fuel rises faster than linearly with power\n",
-    "power_pts = [0, 30, 60, 100]\n",
-    "fuel_pts = [0, 36, 84, 170]  # slopes 1.2, 1.6, 2.15 (convex)\n",
+    "# Same curve as §1, now read as power = g(fuel) (concave production function)\n",
     "pwf = m.add_piecewise_formulation(\n",
-    "    (fuel, fuel_pts, \">=\"),  # fuel ≥ f(power) — over-fuelling allowed\n",
-    "    (power, power_pts),  # equality role\n",
+    "    (power, [0, 30, 60, 100], \"<=\"),  # power ≤ f(fuel) — curtailable output\n",
+    "    (fuel, [0, 36, 84, 170]),  # equality role\n",
     ")\n",
-    "m.add_constraints(power == xr.DataArray([30, 80, 100], coords=[time]))\n",
-    "m.add_objective(fuel.sum())  # minimise fuel against the lower bound\n",
+    "m.add_constraints(power == demand)\n",
+    "m.add_objective(fuel.sum())  # minimise fuel against the bound\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
@@ -248,7 +246,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_curve(power_pts, fuel_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
+    "plot_curve(\n",
+    "    [0, 30, 60, 100],\n",
+    "    [0, 36, 84, 170],\n",
+    "    m.solution[\"power\"].values,\n",
+    "    m.solution[\"fuel\"].values,\n",
+    ");"
    ]
   },
   {

From f516fc9e802aa63c80aa57a7655aeb28cd70b76a Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:46:26 +0200
Subject: [PATCH 10/29] docs(piecewise): align tutorial 2 breakpoints with the
 rest
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Inequality-bounds tutorial now uses the same curve as §1/§4 of the
linear-constraints tutorial and the rst Quick Start:

  fuel_pts  = [0, 36, 84, 170]
  power_pts = [0, 30, 60, 100]

This concave production curve has f(60) = 45, so the verification cell
checks power = 45 at fuel = 60.  The hypograph-visualisation test
points are scaled to the new axes: (60, 30) under, (60, 45) on,
(60, 55) above, (180, 50) beyond domain.  The non-convex fallback
example also uses the [0, 170] fuel domain.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb | 59 +++++++++++-----------
 1 file changed, 29 insertions(+), 30 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 1e4d04a3..db767677 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -67,8 +67,8 @@
    },
    "outputs": [],
    "source": [
-    "fuel_pts = np.array([0.0, 10.0, 20.0, 30.0])\n",
-    "power_pts = np.array([0.0, 20.0, 30.0, 35.0])  # slopes 2, 1, 0.5 (concave)\n",
+    "fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])\n",
+    "power_pts = np.array([0.0, 30.0, 60.0, 100.0])  # slopes 0.83, 0.62, 0.46 (concave)\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(5, 4))\n",
     "ax.plot(fuel_pts, power_pts, \"o-\", color=\"C0\", lw=2)\n",
@@ -83,7 +83,7 @@
    "source": [
     "## Three methods, identical feasible region\n",
     "\n",
-    "With `power` bounded `<=` and our concave curve, the three methods give the **same** feasible region for `fuel ∈ [0, 30]`:\n",
+    "With `power` bounded `<=` and our concave curve, the three methods give the **same** feasible region for `fuel ∈ [0, 170]`:\n",
     "\n",
     "- **`method=\"lp\"`** — tangent lines + domain bounds. No auxiliary variables.\n",
     "- **`method=\"sos2\"`** — lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
@@ -91,7 +91,7 @@
     "\n",
     "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable — it's always preferable because it's pure LP.\n",
     "\n",
-    "Let's verify they produce the same solution at `fuel=15`."
+    "Let's verify they produce the same solution at `fuel=60`, where `f(60)=45`."
    ]
   },
   {
@@ -107,11 +107,11 @@
    "source": [
     "def solve(method, fuel_val):\n",
     "    m = linopy.Model()\n",
-    "    fuel = m.add_variables(lower=0, upper=30, name=\"fuel\")\n",
-    "    power = m.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "    fuel = m.add_variables(lower=0, upper=170, name=\"fuel\")\n",
+    "    power = m.add_variables(lower=0, upper=100, name=\"power\")\n",
     "    m.add_piecewise_formulation(\n",
     "        (power, power_pts, \"<=\"),  # power ≤ f(fuel) — curtailable output\n",
-    "        (fuel, fuel_pts),  # equality role (domain-bounded to [0, 30])\n",
+    "        (fuel, fuel_pts),  # equality role (domain-bounded to [0, 170])\n",
     "        method=method,\n",
     "    )\n",
     "    m.add_constraints(fuel == fuel_val)\n",
@@ -121,7 +121,7 @@
     "\n",
     "\n",
     "for method in [\"lp\", \"sos2\", \"incremental\"]:\n",
-    "    power_val, vars_, cons_ = solve(method, 15)\n",
+    "    power_val, vars_, cons_ = solve(method, 60)\n",
     "    print(f\"{method:12}: power={power_val:.2f}  vars={vars_}  cons={cons_}\")"
    ]
   },
@@ -129,7 +129,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "All three give `power=25` at `fuel=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.\n",
+    "All three give `power=45` at `fuel=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.\n",
     "\n",
     "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."
    ]
@@ -142,14 +142,14 @@
     "\n",
     "The feasible region for `(fuel, power)` with `power` bounded `<=` is the **hypograph** of `f` restricted to the fuel domain:\n",
     "\n",
-    "$$\\{ (\\mathrm{fuel}, \\mathrm{power}) : 0 \\le \\mathrm{fuel} \\le 30,\\ \\mathrm{power} \\le f(\\mathrm{fuel}) \\}$$\n",
+    "$$\\{ (\\mathrm{fuel}, \\mathrm{power}) : 0 \\le \\mathrm{fuel} \\le 170,\\ \\mathrm{power} \\le f(\\mathrm{fuel}) \\}$$\n",
     "\n",
     "Below we colour feasible points green, infeasible ones red:\n",
     "\n",
-    "- `(15, 15)` — under the curve, `15 ≤ f(15)=25` ✓\n",
-    "- `(15, 25)` — on the curve ✓\n",
-    "- `(15, 29)` — above `f(15)`, infeasible ✗\n",
-    "- `(35, 20)` — fuel beyond domain, infeasible ✗"
+    "- `(60, 30)` — under the curve, `30 ≤ f(60)=45` ✓\n",
+    "- `(60, 45)` — on the curve ✓\n",
+    "- `(60, 55)` — above `f(60)`, infeasible ✗\n",
+    "- `(180, 50)` — fuel beyond domain, infeasible ✗"
    ]
   },
   {
@@ -169,10 +169,10 @@
     "    return py <= np.interp(fx, fuel_pts, power_pts)\n",
     "\n",
     "\n",
-    "xx, yy = np.meshgrid(np.linspace(-2, 38, 200), np.linspace(-5, 45, 200))\n",
+    "xx, yy = np.meshgrid(np.linspace(-10, 200, 200), np.linspace(-10, 120, 200))\n",
     "region = np.vectorize(in_hypograph)(xx, yy)\n",
     "\n",
-    "test_points = [(15, 15), (15, 25), (15, 29), (35, 20)]\n",
+    "test_points = [(60, 30), (60, 45), (60, 55), (180, 50)]\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(6, 5))\n",
     "ax.contourf(xx, yy, region, levels=[0.5, 1], colors=[\"lightsteelblue\"], alpha=0.5)\n",
@@ -186,7 +186,7 @@
     "ax.set(\n",
     "    xlabel=\"fuel\",\n",
     "    ylabel=\"power\",\n",
-    "    title=\"sign='<=' feasible region — hypograph of f(fuel) on [0, 30]\",\n",
+    "    title=\"sign='<=' feasible region — hypograph of f(fuel) on [0, 170]\",\n",
     ")\n",
     "ax.grid(alpha=0.3)\n",
     "ax.legend()\n",
@@ -229,32 +229,31 @@
    "outputs": [],
    "source": [
     "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\n",
-    "fuel_nc = [0, 10, 20, 30]\n",
-    "power_nc = [0, 20, 10, 30]  # slopes change sign → mixed convexity\n",
+    "fuel_nc = [0, 50, 100, 170]\n",
+    "power_nc = [0, 60, 30, 100]  # slopes change sign → mixed convexity\n",
     "\n",
     "m1 = linopy.Model()\n",
-    "fuel1 = m1.add_variables(lower=0, upper=30, name=\"fuel\")\n",
-    "power1 = m1.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "fuel1 = m1.add_variables(lower=0, upper=170, name=\"fuel\")\n",
+    "power1 = m1.add_variables(lower=0, upper=100, name=\"power\")\n",
     "f1 = m1.add_piecewise_formulation((power1, power_nc, \"<=\"), (fuel1, fuel_nc))\n",
     "print(f\"non-convex + '<=' → {f1.method}\")\n",
     "\n",
     "# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP\n",
-    "fuel_cc = [0, 10, 20, 30]\n",
-    "power_cc = [0, 20, 30, 35]  # concave\n",
-    "\n",
     "m2 = linopy.Model()\n",
-    "fuel2 = m2.add_variables(lower=0, upper=30, name=\"fuel\")\n",
-    "power2 = m2.add_variables(lower=0, upper=40, name=\"power\")\n",
-    "f2 = m2.add_piecewise_formulation((power2, power_cc, \">=\"), (fuel2, fuel_cc))\n",
+    "fuel2 = m2.add_variables(lower=0, upper=170, name=\"fuel\")\n",
+    "power2 = m2.add_variables(lower=0, upper=100, name=\"power\")\n",
+    "f2 = m2.add_piecewise_formulation(\n",
+    "    (power2, list(power_pts), \">=\"), (fuel2, list(fuel_pts))\n",
+    ")\n",
     "print(f\"concave + '>='    → {f2.method}\")\n",
     "\n",
     "# 3. Explicit method=\"lp\" with mismatched curvature raises\n",
     "try:\n",
     "    m3 = linopy.Model()\n",
-    "    fuel3 = m3.add_variables(lower=0, upper=30, name=\"fuel\")\n",
-    "    power3 = m3.add_variables(lower=0, upper=40, name=\"power\")\n",
+    "    fuel3 = m3.add_variables(lower=0, upper=170, name=\"fuel\")\n",
+    "    power3 = m3.add_variables(lower=0, upper=100, name=\"power\")\n",
     "    m3.add_piecewise_formulation(\n",
-    "        (power3, power_cc, \">=\"), (fuel3, fuel_cc), method=\"lp\"\n",
+    "        (power3, list(power_pts), \">=\"), (fuel3, list(fuel_pts)), method=\"lp\"\n",
     "    )\n",
     "except ValueError as e:\n",
     "    print(f\"lp(concave, '>=')  → raises: {e}\")"

From dbcea0bc0a066065864ec29f5ca03a6fe55aa1fd Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:50:27 +0200
Subject: [PATCH 11/29] docs(piecewise): silence EvolvingAPIWarning in tutorial
 setup cells
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The warning fires once per session on first use of
add_piecewise_formulation / Slopes / tangent_lines.  Useful in user
code, but in a tutorial it's a distracting wall of yellow above the
first solve.  Silenced in the imports cell of each notebook via
warnings.filterwarnings — the canonical recipe from the warning
message itself.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb  |   7 +-
 examples/piecewise-linear-constraints.ipynb | 151 ++++++++++++--------
 2 files changed, 94 insertions(+), 64 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index db767677..79c497be 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -41,10 +41,15 @@
    },
    "outputs": [],
    "source": [
+    "import warnings\n",
+    "\n",
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
-    "import linopy"
+    "import linopy\n",
+    "\n",
+    "# Silence the evolving-API warning for cleaner tutorial output.\n",
+    "warnings.filterwarnings(\"ignore\", category=linopy.EvolvingAPIWarning)"
    ]
   },
   {
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index d54bd2c1..0b3d737e 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -31,21 +31,24 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:58.302751Z",
-     "start_time": "2026-04-22T23:31:58.299283Z"
+     "end_time": "2026-05-11T12:48:45.162889Z",
+     "start_time": "2026-05-11T12:48:40.871619Z"
     }
    },
-   "outputs": [],
    "source": [
+    "import warnings\n",
+    "\n",
     "import matplotlib.pyplot as plt\n",
     "import pandas as pd\n",
     "import xarray as xr\n",
     "\n",
     "import linopy\n",
     "\n",
+    "# Silence the evolving-API warning for cleaner tutorial output.\n",
+    "warnings.filterwarnings(\"ignore\", category=linopy.EvolvingAPIWarning)\n",
+    "\n",
     "time = pd.Index([1, 2, 3], name=\"time\")\n",
     "\n",
     "\n",
@@ -59,7 +62,9 @@
     "    ax.set(xlabel=xlabel, ylabel=ylabel)\n",
     "    ax.legend()\n",
     "    return ax"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -72,14 +77,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:58.464773Z",
-     "start_time": "2026-04-22T23:31:58.310016Z"
+     "end_time": "2026-05-11T12:48:45.817102Z",
+     "start_time": "2026-05-11T12:48:45.170166Z"
     }
    },
-   "outputs": [],
    "source": [
     "demand = xr.DataArray([50, 80, 30], coords=[time])\n",
     "\n",
@@ -96,21 +99,23 @@
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:58.532078Z",
-     "start_time": "2026-04-22T23:31:58.473509Z"
+     "end_time": "2026-05-11T12:48:45.925332Z",
+     "start_time": "2026-05-11T12:48:45.825967Z"
     }
    },
-   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -131,14 +136,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:58.952185Z",
-     "start_time": "2026-04-22T23:31:58.537015Z"
+     "end_time": "2026-05-11T12:48:46.288470Z",
+     "start_time": "2026-05-11T12:48:45.983599Z"
     }
    },
-   "outputs": [],
    "source": [
     "def solve_method(method):\n",
     "    m = linopy.Model()\n",
@@ -152,7 +155,9 @@
     "\n",
     "\n",
     "pd.DataFrame({m: solve_method(m) for m in [\"auto\", \"sos2\", \"incremental\"]})"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -167,15 +172,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:59.092539Z",
-     "start_time": "2026-04-22T23:31:58.956054Z"
-    },
-    "scrolled": true
+     "end_time": "2026-05-11T12:48:46.422121Z",
+     "start_time": "2026-05-11T12:48:46.294571Z"
+    }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=80, coords=[time])\n",
@@ -190,7 +193,9 @@
     "m.add_objective(cost.sum() + 10 * backup.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"cost\", \"backup\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -214,14 +219,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:59.210868Z",
-     "start_time": "2026-04-22T23:31:59.098774Z"
+     "end_time": "2026-05-11T12:48:46.488654Z",
+     "start_time": "2026-05-11T12:48:46.427536Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -238,13 +241,18 @@
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-05-11T12:48:46.540310Z",
+     "start_time": "2026-05-11T12:48:46.501095Z"
+    }
+   },
    "source": [
     "plot_curve(\n",
     "    [0, 30, 60, 100],\n",
@@ -252,7 +260,9 @@
     "    m.solution[\"power\"].values,\n",
     "    m.solution[\"fuel\"].values,\n",
     ");"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -268,14 +278,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:59.422636Z",
-     "start_time": "2026-04-22T23:31:59.232150Z"
+     "end_time": "2026-05-11T12:48:46.678189Z",
+     "start_time": "2026-05-11T12:48:46.544709Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "p_min, p_max = 30, 100\n",
@@ -297,16 +305,23 @@
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"commit\", \"power\", \"fuel\", \"backup\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-05-11T12:48:46.722833Z",
+     "start_time": "2026-05-11T12:48:46.683526Z"
+    }
+   },
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -319,14 +334,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:59.598540Z",
-     "start_time": "2026-04-22T23:31:59.433551Z"
+     "end_time": "2026-05-11T12:48:46.835780Z",
+     "start_time": "2026-05-11T12:48:46.732959Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -345,13 +358,18 @@
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\", \"heat\"]].to_pandas().round(2)"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-05-11T12:48:46.911022Z",
+     "start_time": "2026-05-11T12:48:46.849158Z"
+    }
+   },
    "source": [
     "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
     "plot_curve(\n",
@@ -365,7 +383,9 @@
     "    ylabel=\"heat\",\n",
     "    ax=axes[1],\n",
     ");"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -378,14 +398,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-04-22T23:31:59.801734Z",
-     "start_time": "2026-04-22T23:31:59.606692Z"
+     "end_time": "2026-05-11T12:48:47.052842Z",
+     "start_time": "2026-05-11T12:48:46.923409Z"
     }
    },
-   "outputs": [],
    "source": [
     "gens = pd.Index([\"gas\", \"coal\"], name=\"gen\")\n",
     "x_gen = linopy.breakpoints(\n",
@@ -403,7 +421,9 @@
     "m.add_objective(fuel.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\"]].to_dataframe()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -416,9 +436,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-05-11T12:48:47.150365Z",
+     "start_time": "2026-05-11T12:48:47.058209Z"
+    }
+   },
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -433,7 +456,9 @@
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",

From 3d84a908b7afe442e2b6e61d4ecc5175c401cccf Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 14:59:51 +0200
Subject: [PATCH 12/29] =?UTF-8?q?docs(piecewise):=20switch=20unified=20ine?=
 =?UTF-8?q?quality=20example=20to=20fuel=20=E2=89=A5=20f(power)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Heat-rate framing — `fuel = f(power)` — is native to power-systems
modelling, so the inequality variant reads naturally as `fuel ≥
f(power)`: the curve is the design minimum, over-fuelling is admissible
but wasteful, and minimising fuel pulls the operating point onto the
curve.  Same feasible region as the previous `power ≤ f(fuel)` form
under inversion; the LP path uses the convex + `>=` (epigraph) half
instead of concave + `<=` (hypograph).

- rst Quick Start: bound fuel ≥ f(power), same breakpoints as §1.
- Tutorial 1 §4: same.
- Tutorial 2: setup curve labelled "convex heat-rate curve f(power)";
  solve checks fuel=84 at power=60; visualisation is now the epigraph
  with test points (60, 100), (60, 84), (60, 70), (120, 100); fallback
  cell exercises the dual mismatched-curvature cases (convex + `<=`).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        |  15 +-
 examples/piecewise-inequality-bounds.ipynb  | 104 +++++++-------
 examples/piecewise-linear-constraints.ipynb | 144 ++++++++++----------
 3 files changed, 131 insertions(+), 132 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 98187ee0..ece5c823 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -49,15 +49,14 @@ Quick Start
 
 .. code-block:: python
 
-    # power <= f(fuel) on the same curve as above, but read as a
-    # production function: power output is bounded by what the fuel
-    # input can support, and the unit may run below the max
-    # (curtailable output).  "auto" picks the cheapest correct
-    # formulation: pure LP (chord constraints) here, since concave +
-    # "<=" is LP-applicable; SOS2/incremental otherwise.
+    # fuel >= f(power) on the same heat-rate curve as above.  Over-
+    # fuelling is physically admissible but wasteful, so minimising
+    # fuel pulls the operating point onto the curve.  "auto" picks the
+    # cheapest correct formulation: pure LP (chord constraints) here,
+    # since convex + ">=" is LP-applicable; SOS2/incremental otherwise.
     m.add_piecewise_formulation(
-        (power, [0, 30, 60, 100], "<="),
-        (fuel, [0, 36, 84, 170]),
+        (fuel, [0, 36, 84, 170], ">="),
+        (power, [0, 30, 60, 100]),
     )
 
 Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with
diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 79c497be..728f1ee8 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -10,8 +10,8 @@
     "\n",
     "```python\n",
     "m.add_piecewise_formulation(\n",
-    "    (power, power_pts, \"<=\"),  # power ≤ f(fuel) — output curtailable\n",
-    "    (fuel,  fuel_pts),         # equality role\n",
+    "    (fuel,  fuel_pts, \">=\"),   # fuel ≥ f(power) — over-fuelling admissible\n",
+    "    (power, power_pts),        # equality role\n",
     ")\n",
     "```\n",
     "\n",
@@ -56,9 +56,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Setup — a concave production curve\n",
+    "## Setup — a convex heat-rate curve\n",
     "\n",
-    "A concave, monotonically increasing curve maps fuel input to the maximum power the plant can deliver.  Bounding `power` by this curve with `<=` says the unit may run *at or below* the maximum — power output is curtailable for a given fuel draw.  Concave + `<=` is exactly the combination that lets the LP method apply."
+    "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."
    ]
   },
   {
@@ -72,12 +72,12 @@
    },
    "outputs": [],
    "source": [
-    "fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])\n",
-    "power_pts = np.array([0.0, 30.0, 60.0, 100.0])  # slopes 0.83, 0.62, 0.46 (concave)\n",
+    "power_pts = np.array([0.0, 30.0, 60.0, 100.0])\n",
+    "fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])  # slopes 1.2, 1.6, 2.15 (convex)\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(5, 4))\n",
-    "ax.plot(fuel_pts, power_pts, \"o-\", color=\"C0\", lw=2)\n",
-    "ax.set(xlabel=\"fuel\", ylabel=\"power\", title=\"Concave production curve f(fuel)\")\n",
+    "ax.plot(power_pts, fuel_pts, \"o-\", color=\"C0\", lw=2)\n",
+    "ax.set(xlabel=\"power\", ylabel=\"fuel\", title=\"Convex heat-rate curve f(power)\")\n",
     "ax.grid(alpha=0.3)\n",
     "plt.tight_layout()"
    ]
@@ -88,7 +88,7 @@
    "source": [
     "## Three methods, identical feasible region\n",
     "\n",
-    "With `power` bounded `<=` and our concave curve, the three methods give the **same** feasible region for `fuel ∈ [0, 170]`:\n",
+    "With `fuel` bounded `>=` and our convex curve, the three methods give the **same** feasible region for `power ∈ [0, 100]`:\n",
     "\n",
     "- **`method=\"lp\"`** — tangent lines + domain bounds. No auxiliary variables.\n",
     "- **`method=\"sos2\"`** — lambdas + SOS2 + a split link: equality for the equality-signed tuple, signed for the bounded one. Solver picks the piece.\n",
@@ -96,7 +96,7 @@
     "\n",
     "`method=\"auto\"` dispatches to `\"lp\"` whenever applicable — it's always preferable because it's pure LP.\n",
     "\n",
-    "Let's verify they produce the same solution at `fuel=60`, where `f(60)=45`."
+    "Let's verify they produce the same solution at `power=60`, where `f(60)=84`."
    ]
   },
   {
@@ -110,31 +110,31 @@
    },
    "outputs": [],
    "source": [
-    "def solve(method, fuel_val):\n",
+    "def solve(method, power_val):\n",
     "    m = linopy.Model()\n",
-    "    fuel = m.add_variables(lower=0, upper=170, name=\"fuel\")\n",
     "    power = m.add_variables(lower=0, upper=100, name=\"power\")\n",
+    "    fuel = m.add_variables(lower=0, upper=200, name=\"fuel\")\n",
     "    m.add_piecewise_formulation(\n",
-    "        (power, power_pts, \"<=\"),  # power ≤ f(fuel) — curtailable output\n",
-    "        (fuel, fuel_pts),  # equality role (domain-bounded to [0, 170])\n",
+    "        (fuel, fuel_pts, \">=\"),  # fuel ≥ f(power) — over-fuelling admissible\n",
+    "        (power, power_pts),  # equality role (domain-bounded to [0, 100])\n",
     "        method=method,\n",
     "    )\n",
-    "    m.add_constraints(fuel == fuel_val)\n",
-    "    m.add_objective(-power)  # maximise power output against the upper bound\n",
+    "    m.add_constraints(power == power_val)\n",
+    "    m.add_objective(fuel)  # minimise fuel against the lower bound\n",
     "    m.solve(solver_name=\"highs\", output_flag=False)\n",
-    "    return float(m.solution[\"power\"]), list(m.variables), list(m.constraints)\n",
+    "    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n",
     "\n",
     "\n",
     "for method in [\"lp\", \"sos2\", \"incremental\"]:\n",
-    "    power_val, vars_, cons_ = solve(method, 60)\n",
-    "    print(f\"{method:12}: power={power_val:.2f}  vars={vars_}  cons={cons_}\")"
+    "    fuel_val, vars_, cons_ = solve(method, 60)\n",
+    "    print(f\"{method:12}: fuel={fuel_val:.2f}  vars={vars_}  cons={cons_}\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "All three give `power=45` at `fuel=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.\n",
+    "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.\n",
     "\n",
     "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."
    ]
@@ -145,16 +145,16 @@
    "source": [
     "## Visualising the feasible region\n",
     "\n",
-    "The feasible region for `(fuel, power)` with `power` bounded `<=` is the **hypograph** of `f` restricted to the fuel domain:\n",
+    "The feasible region for `(power, fuel)` with `fuel` bounded `>=` is the **epigraph** of `f` restricted to the power domain:\n",
     "\n",
-    "$$\\{ (\\mathrm{fuel}, \\mathrm{power}) : 0 \\le \\mathrm{fuel} \\le 170,\\ \\mathrm{power} \\le f(\\mathrm{fuel}) \\}$$\n",
+    "$$\\{ (\\mathrm{power}, \\mathrm{fuel}) : 0 \\le \\mathrm{power} \\le 100,\\ \\mathrm{fuel} \\ge f(\\mathrm{power}) \\}$$\n",
     "\n",
     "Below we colour feasible points green, infeasible ones red:\n",
     "\n",
-    "- `(60, 30)` — under the curve, `30 ≤ f(60)=45` ✓\n",
-    "- `(60, 45)` — on the curve ✓\n",
-    "- `(60, 55)` — above `f(60)`, infeasible ✗\n",
-    "- `(180, 50)` — fuel beyond domain, infeasible ✗"
+    "- `(60, 100)` — above the curve, `100 ≥ f(60)=84` ✓\n",
+    "- `(60, 84)` — on the curve ✓\n",
+    "- `(60, 70)` — below `f(60)`, infeasible ✗\n",
+    "- `(120, 100)` — power beyond domain, infeasible ✗"
    ]
   },
   {
@@ -168,30 +168,30 @@
    },
    "outputs": [],
    "source": [
-    "def in_hypograph(fx, py):\n",
-    "    if fx < fuel_pts[0] or fx > fuel_pts[-1]:\n",
+    "def in_epigraph(px, fy):\n",
+    "    if px < power_pts[0] or px > power_pts[-1]:\n",
     "        return False\n",
-    "    return py <= np.interp(fx, fuel_pts, power_pts)\n",
+    "    return fy >= np.interp(px, power_pts, fuel_pts)\n",
     "\n",
     "\n",
-    "xx, yy = np.meshgrid(np.linspace(-10, 200, 200), np.linspace(-10, 120, 200))\n",
-    "region = np.vectorize(in_hypograph)(xx, yy)\n",
+    "xx, yy = np.meshgrid(np.linspace(-10, 130, 200), np.linspace(-10, 200, 200))\n",
+    "region = np.vectorize(in_epigraph)(xx, yy)\n",
     "\n",
-    "test_points = [(60, 30), (60, 45), (60, 55), (180, 50)]\n",
+    "test_points = [(60, 100), (60, 84), (60, 70), (120, 100)]\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(6, 5))\n",
     "ax.contourf(xx, yy, region, levels=[0.5, 1], colors=[\"lightsteelblue\"], alpha=0.5)\n",
-    "ax.plot(fuel_pts, power_pts, \"o-\", color=\"C0\", lw=2, label=\"f(fuel)\")\n",
-    "for fx, py in test_points:\n",
-    "    feas = in_hypograph(fx, py)\n",
+    "ax.plot(power_pts, fuel_pts, \"o-\", color=\"C0\", lw=2, label=\"f(power)\")\n",
+    "for px, fy in test_points:\n",
+    "    feas = in_epigraph(px, fy)\n",
     "    ax.scatter(\n",
-    "        [fx], [py], color=\"green\" if feas else \"red\", zorder=5, s=80, edgecolors=\"black\"\n",
+    "        [px], [fy], color=\"green\" if feas else \"red\", zorder=5, s=80, edgecolors=\"black\"\n",
     "    )\n",
-    "    ax.annotate(f\"({fx}, {py})\", (fx, py), textcoords=\"offset points\", xytext=(8, 5))\n",
+    "    ax.annotate(f\"({px}, {fy})\", (px, fy), textcoords=\"offset points\", xytext=(8, 5))\n",
     "ax.set(\n",
-    "    xlabel=\"fuel\",\n",
-    "    ylabel=\"power\",\n",
-    "    title=\"sign='<=' feasible region — hypograph of f(fuel) on [0, 170]\",\n",
+    "    xlabel=\"power\",\n",
+    "    ylabel=\"fuel\",\n",
+    "    title=\"sign='>=' feasible region — epigraph of f(power) on [0, 100]\",\n",
     ")\n",
     "ax.grid(alpha=0.3)\n",
     "ax.legend()\n",
@@ -234,34 +234,34 @@
    "outputs": [],
    "source": [
     "# 1. Non-convex curve: auto falls back (LP relaxation would be loose)\n",
-    "fuel_nc = [0, 50, 100, 170]\n",
-    "power_nc = [0, 60, 30, 100]  # slopes change sign → mixed convexity\n",
+    "power_nc = [0, 30, 60, 100]\n",
+    "fuel_nc = [0, 50, 30, 80]  # slopes change sign → mixed convexity\n",
     "\n",
     "m1 = linopy.Model()\n",
-    "fuel1 = m1.add_variables(lower=0, upper=170, name=\"fuel\")\n",
     "power1 = m1.add_variables(lower=0, upper=100, name=\"power\")\n",
-    "f1 = m1.add_piecewise_formulation((power1, power_nc, \"<=\"), (fuel1, fuel_nc))\n",
-    "print(f\"non-convex + '<=' → {f1.method}\")\n",
+    "fuel1 = m1.add_variables(lower=0, upper=200, name=\"fuel\")\n",
+    "f1 = m1.add_piecewise_formulation((fuel1, fuel_nc, \">=\"), (power1, power_nc))\n",
+    "print(f\"non-convex + '>=' → {f1.method}\")\n",
     "\n",
-    "# 2. Concave curve + sign='>=': LP would be loose → auto falls back to MIP\n",
+    "# 2. Convex curve + sign='<=': LP would be loose → auto falls back to MIP\n",
     "m2 = linopy.Model()\n",
-    "fuel2 = m2.add_variables(lower=0, upper=170, name=\"fuel\")\n",
     "power2 = m2.add_variables(lower=0, upper=100, name=\"power\")\n",
+    "fuel2 = m2.add_variables(lower=0, upper=200, name=\"fuel\")\n",
     "f2 = m2.add_piecewise_formulation(\n",
-    "    (power2, list(power_pts), \">=\"), (fuel2, list(fuel_pts))\n",
+    "    (fuel2, list(fuel_pts), \"<=\"), (power2, list(power_pts))\n",
     ")\n",
-    "print(f\"concave + '>='    → {f2.method}\")\n",
+    "print(f\"convex + '<='     → {f2.method}\")\n",
     "\n",
     "# 3. Explicit method=\"lp\" with mismatched curvature raises\n",
     "try:\n",
     "    m3 = linopy.Model()\n",
-    "    fuel3 = m3.add_variables(lower=0, upper=170, name=\"fuel\")\n",
     "    power3 = m3.add_variables(lower=0, upper=100, name=\"power\")\n",
+    "    fuel3 = m3.add_variables(lower=0, upper=200, name=\"fuel\")\n",
     "    m3.add_piecewise_formulation(\n",
-    "        (power3, list(power_pts), \">=\"), (fuel3, list(fuel_pts)), method=\"lp\"\n",
+    "        (fuel3, list(fuel_pts), \"<=\"), (power3, list(power_pts)), method=\"lp\"\n",
     "    )\n",
     "except ValueError as e:\n",
-    "    print(f\"lp(concave, '>=')  → raises: {e}\")"
+    "    print(f\"lp(convex, '<=')  → raises: {e}\")"
    ]
   },
   {
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 0b3d737e..3648411e 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -31,12 +31,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:45.162889Z",
-     "start_time": "2026-05-11T12:48:40.871619Z"
+     "end_time": "2026-05-11T12:50:48.865977Z",
+     "start_time": "2026-05-11T12:50:48.861703Z"
     }
    },
+   "outputs": [],
    "source": [
     "import warnings\n",
     "\n",
@@ -62,9 +64,7 @@
     "    ax.set(xlabel=xlabel, ylabel=ylabel)\n",
     "    ax.legend()\n",
     "    return ax"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -77,12 +77,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:45.817102Z",
-     "start_time": "2026-05-11T12:48:45.170166Z"
+     "end_time": "2026-05-11T12:51:46.043507Z",
+     "start_time": "2026-05-11T12:51:45.930380Z"
     }
    },
+   "outputs": [],
    "source": [
     "demand = xr.DataArray([50, 80, 30], coords=[time])\n",
     "\n",
@@ -99,23 +101,21 @@
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:45.925332Z",
-     "start_time": "2026-05-11T12:48:45.825967Z"
+     "end_time": "2026-05-11T12:50:49.027645Z",
+     "start_time": "2026-05-11T12:50:48.982161Z"
     }
    },
+   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -136,12 +136,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.288470Z",
-     "start_time": "2026-05-11T12:48:45.983599Z"
+     "end_time": "2026-05-11T12:50:49.330961Z",
+     "start_time": "2026-05-11T12:50:49.039795Z"
     }
    },
+   "outputs": [],
    "source": [
     "def solve_method(method):\n",
     "    m = linopy.Model()\n",
@@ -155,9 +157,7 @@
     "\n",
     "\n",
     "pd.DataFrame({m: solve_method(m) for m in [\"auto\", \"sos2\", \"incremental\"]})"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -172,13 +172,15 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
-    "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.422121Z",
-     "start_time": "2026-05-11T12:48:46.294571Z"
-    }
+     "end_time": "2026-05-11T12:50:49.487115Z",
+     "start_time": "2026-05-11T12:50:49.347614Z"
+    },
+    "scrolled": true
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=80, coords=[time])\n",
@@ -193,9 +195,7 @@
     "m.add_objective(cost.sum() + 10 * backup.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"cost\", \"backup\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -214,45 +214,47 @@
     "\n",
     "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "Below: the same curve as §1, but read as a production function `power ≤ f(fuel)` — the unit *can* run below the maximum power for a given fuel draw (output curtailable). Minimising fuel against a power demand pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
+    "Below: the same heat-rate curve as §1, now read as `fuel ≥ f(power)` — over-fuelling is admissible but wasteful (the curve is the design minimum), so an objective that minimises fuel pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
    ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.488654Z",
-     "start_time": "2026-05-11T12:48:46.427536Z"
+     "end_time": "2026-05-11T12:50:49.556562Z",
+     "start_time": "2026-05-11T12:50:49.492505Z"
     }
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
     "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
     "\n",
-    "# Same curve as §1, now read as power = g(fuel) (concave production function)\n",
+    "# Same convex heat-rate curve as §1, now bounded with \">=\"\n",
     "pwf = m.add_piecewise_formulation(\n",
-    "    (power, [0, 30, 60, 100], \"<=\"),  # power ≤ f(fuel) — curtailable output\n",
-    "    (fuel, [0, 36, 84, 170]),  # equality role\n",
+    "    (fuel, [0, 36, 84, 170], \">=\"),  # fuel ≥ f(power) — over-fuelling allowed\n",
+    "    (power, [0, 30, 60, 100]),  # equality role\n",
     ")\n",
     "m.add_constraints(power == demand)\n",
-    "m.add_objective(fuel.sum())  # minimise fuel against the bound\n",
+    "m.add_objective(fuel.sum())  # minimise fuel against the lower bound\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.540310Z",
-     "start_time": "2026-05-11T12:48:46.501095Z"
+     "end_time": "2026-05-11T12:50:49.603928Z",
+     "start_time": "2026-05-11T12:50:49.566231Z"
     }
    },
+   "outputs": [],
    "source": [
     "plot_curve(\n",
     "    [0, 30, 60, 100],\n",
@@ -260,9 +262,7 @@
     "    m.solution[\"power\"].values,\n",
     "    m.solution[\"fuel\"].values,\n",
     ");"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -278,12 +278,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.678189Z",
-     "start_time": "2026-05-11T12:48:46.544709Z"
+     "end_time": "2026-05-11T12:50:49.746059Z",
+     "start_time": "2026-05-11T12:50:49.613733Z"
     }
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "p_min, p_max = 30, 100\n",
@@ -305,23 +307,21 @@
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"commit\", \"power\", \"fuel\", \"backup\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.722833Z",
-     "start_time": "2026-05-11T12:48:46.683526Z"
+     "end_time": "2026-05-11T12:50:49.792782Z",
+     "start_time": "2026-05-11T12:50:49.751516Z"
     }
    },
+   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -334,12 +334,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.835780Z",
-     "start_time": "2026-05-11T12:48:46.732959Z"
+     "end_time": "2026-05-11T12:50:49.914868Z",
+     "start_time": "2026-05-11T12:50:49.802514Z"
     }
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -358,18 +360,18 @@
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\", \"heat\"]].to_pandas().round(2)"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:46.911022Z",
-     "start_time": "2026-05-11T12:48:46.849158Z"
+     "end_time": "2026-05-11T12:50:49.993753Z",
+     "start_time": "2026-05-11T12:50:49.922140Z"
     }
    },
+   "outputs": [],
    "source": [
     "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
     "plot_curve(\n",
@@ -383,9 +385,7 @@
     "    ylabel=\"heat\",\n",
     "    ax=axes[1],\n",
     ");"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -398,12 +398,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:47.052842Z",
-     "start_time": "2026-05-11T12:48:46.923409Z"
+     "end_time": "2026-05-11T12:50:50.141147Z",
+     "start_time": "2026-05-11T12:50:50.003565Z"
     }
    },
+   "outputs": [],
    "source": [
     "gens = pd.Index([\"gas\", \"coal\"], name=\"gen\")\n",
     "x_gen = linopy.breakpoints(\n",
@@ -421,9 +423,7 @@
     "m.add_objective(fuel.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\"]].to_dataframe()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -436,12 +436,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:48:47.150365Z",
-     "start_time": "2026-05-11T12:48:47.058209Z"
+     "end_time": "2026-05-11T12:50:50.244428Z",
+     "start_time": "2026-05-11T12:50:50.151699Z"
     }
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -456,9 +458,7 @@
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",

From 825af889145dc27a1cd474b5324cd658cd36f7bf Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 15:04:43 +0200
Subject: [PATCH 13/29] docs(piecewise): add reformulate_sos="auto" to tutorial
 2's solve

HiGHS doesn't natively support SOS constraints, so the sos2 and
incremental cases of the three-methods comparison need
reformulate_sos="auto" (matching tutorial 1's pattern).  Without it,
m.solve(solver_name="highs", ...) raises ValueError on the SOS2 path.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 728f1ee8..0c47ffae 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -121,7 +121,7 @@
     "    )\n",
     "    m.add_constraints(power == power_val)\n",
     "    m.add_objective(fuel)  # minimise fuel against the lower bound\n",
-    "    m.solve(solver_name=\"highs\", output_flag=False)\n",
+    "    m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "    return float(m.solution[\"fuel\"]), list(m.variables), list(m.constraints)\n",
     "\n",
     "\n",

From 1e691ea8c2893ef965dc4ae6c59a9727e1655a7f Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 15:49:23 +0200
Subject: [PATCH 14/29] docs(piecewise): restructure rst + tutorial flow tweaks
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

RST restructure
- Moved "Breakpoint Construction" ahead of "Per-tuple sign".  Readers
  looking up "how do I supply breakpoints" no longer have to scroll
  through the inequality section first.
- Split "Per-tuple sign" — previously a chapter masquerading as a
  section, covering six concepts under one heading — into five proper
  subsections: Roles and restrictions, Geometry, Choice of bounded
  tuple and sign, When is a one-sided bound wanted?, Formulation math.
- Folded the duplicate "Breakpoint inputs" overview (under API) into
  the lead-in of "Breakpoint Construction"; the three building blocks
  (breakpoints, segments, Slopes) are now introduced once.
- Promoted tangent_lines to its own "Chord expressions as a building
  block" subsection under LP — previously a half-hidden paragraph.
- Moved the disjunctive caveat ("method returns 'sos2' but the table
  treats it separately") into a `.. note::` directly under the
  comparison table where the apparent contradiction is.
- Shortened the terminology block at the top to one paragraph (it
  previously defined "segment" hundreds of lines before the disjunctive
  section).
- Centralised N≥3 sign restriction — N-variable linking section now
  cross-references Per-tuple sign instead of redefining it.
- Removed the duplicated active+sign warning under "Advanced Features"
  — single source of truth now in "Per-tuple sign / Formulation math".

Tutorial flow
- Tut 1 §2: trimmed the upfront method-comparison table to a forward
  pointer — the table previously claimed "lp requires sign != ==" and
  "matching curvature" before §3 and §4 had tutorialised those terms.
- Tut 2: dropped the standalone setup plot that showed only the curve.
  The hypograph/epigraph visualisation later in the notebook already
  shows the curve, with the operating points overlaid — saves a
  figure and tightens the narrative.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        | 438 ++++++++++----------
 examples/piecewise-inequality-bounds.ipynb  |  12 +-
 examples/piecewise-linear-constraints.ipynb |  10 +-
 3 files changed, 223 insertions(+), 237 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index ece5c823..6e558772 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -3,23 +3,15 @@
 Piecewise Linear Constraints
 ============================
 
-Piecewise linear (PWL) constraints approximate nonlinear functions as connected
-linear pieces, allowing you to model cost curves, efficiency curves, or
-production functions within a linear programming framework.
-
-**Terminology used in this page:**
-
-- **breakpoint** — a knot on the piecewise curve where the slope can
-  change.  Each tuple supplies a 1-D breakpoint array; the ``i``-th
-  entries across all tuples, taken together, define one knot.  For two
-  tuples this is the usual :math:`(x, y)` pair; for three or more, an
-  ``N``-dim knot.
-- **piece** — a linear part between two adjacent breakpoints on a single
-  connected curve.  ``n`` breakpoints define ``n − 1`` pieces.
-- **segment** — a *disjoint* operating region in the disjunctive
-  formulation, built via the :func:`~linopy.segments` factory.  Within
-  one segment the curve is itself piecewise-linear (made of pieces);
-  between segments there are gaps.
+Piecewise linear (PWL) constraints approximate nonlinear functions as
+connected linear pieces, allowing you to model cost curves, efficiency
+curves, or production functions within a linear programming framework.
+
+Throughout this page: a **breakpoint** is a knot where the slope can
+change; a **piece** is the linear part between adjacent breakpoints; a
+**segment** is a disjoint operating region in the disjunctive
+formulation.  Per-tuple breakpoint arrays are paired by index — the
+``i``-th entries across all tuples together define one knot.
 
 .. contents::
    :local:
@@ -75,7 +67,7 @@ API
 .. code-block:: python
 
     m.add_piecewise_formulation(
-        (expr1, breakpoints1),  # sign defaults to "==" (pinned role)
+        (expr1, breakpoints1),  # sign defaults to "==" (equality role)
         (expr2, breakpoints2, "<="),  # or with an explicit sign
         ...,
         method="auto",  # "auto", "sos2", "incremental", or "lp"
@@ -90,166 +82,20 @@ tuple.  The pure-LP path adds chord and domain constraints only; SOS2,
 incremental, and disjunctive also add interpolation weights and/or
 binaries (see *Formulation Methods* below).
 
-Breakpoint inputs
-~~~~~~~~~~~~~~~~~
-
-Three building blocks with distinct geometric meaning — two factories
-and one class:
-
-- ``breakpoints()`` — values along a single **connected** curve.  Linear
-  pieces between adjacent breakpoints are interpolated continuously.
-- ``segments()`` — **disjoint** operating regions with gaps between them
-  (e.g. forbidden zones).  Builds a 2-D array consumed by the
-  *disjunctive* formulation, where exactly one region is active at a time.
-- :class:`~linopy.Slopes` — per-piece slopes plus an initial ``y0``,
-  deferred until an x grid is supplied.  Inside
-  ``add_piecewise_formulation`` the x grid is borrowed from a sibling
-  tuple; standalone, call :meth:`~linopy.Slopes.to_breakpoints`.
-
-.. code-block:: python
-
-    linopy.breakpoints([0, 50, 100])  # connected
-    linopy.breakpoints({"gen1": [0, 50], "gen2": [0, 80]}, dim="gen")  # per-entity
-    linopy.segments([(0, 10), (50, 100)])  # two disjoint regions
-    linopy.segments({"gen1": [(0, 10)], "gen2": [(0, 80)]}, dim="gen")
-    linopy.Slopes([1.2, 1.4], y0=0)  # deferred — pairs with a sibling tuple
-
-
-Per-tuple sign — equality vs inequality
-----------------------------------------
-
-Each tuple's optional third element is a sign:
-
-- ``"=="`` (default) — **pinned**: the tuple enters as an equality.
-- ``"<="`` / ``">="`` — **bounded**: the expression undershoots /
-  overshoots the curve.
-
-These describe the tuple's *role*, not a geometric property of the
-variable.  What "pinned" actually constrains depends on the count:
-
-- **All tuples pinned (default).**  Shared interpolation weights put
-  the joint :math:`(e_1, \ldots, e_N)` exactly on the curve.
-- **One bounded + one pinned (2 tuples).**  The joint :math:`(x, y)`
-  lies in the hypograph / epigraph (see *Geometry* below).  A single
-  coordinate can't locate a curve point, so the pinned axis's marginal
-  feasible set is just :math:`[x_{\min}, x_{\max}]` — same in LP
-  (enforced directly) and SOS2/incremental (enforced via the weight
-  link).
-- **One bounded + 3+ tuples.**  Not supported (see restrictions below).
-
-.. code-block:: python
-
-    # All-equality: joint (x, y) on the curve.
-    m.add_piecewise_formulation((y, y_pts), (x, x_pts))
-
-    # Bounded: joint (x, y) in the hypograph — y ≤ f(x), x ∈ [x_min, x_max].
-    m.add_piecewise_formulation((y, y_pts, "<="), (x, x_pts))
-
-    # Bounded: joint (x, y) in the epigraph — y ≥ f(x), x ∈ [x_min, x_max].
-    m.add_piecewise_formulation((y, y_pts, ">="), (x, x_pts))
-
-    # 3-variable all-equality (CHP): joint (power, fuel, heat) on the curve.
-    m.add_piecewise_formulation((power, p_pts), (fuel, f_pts), (heat, h_pts))
-
-**Restrictions (current):**
-
-- At most one tuple may carry a non-equality sign — a single bounded side.
-- With **3 or more** tuples, all signs must be ``"=="``.
-
-Multi-bounded and N≥3-inequality use cases aren't supported yet.  If you
-have a concrete use case, please open an issue at
-https://github.com/PyPSA/linopy/issues so we can scope it properly.
-
-**Formulation.**  For methods that introduce shared interpolation
-weights (SOS2 and incremental — see below), only the link constraint
-between the weights and the bounded expression changes.  Write the
-method-specific weighted sum of breakpoints for tuple :math:`j` as
-:math:`W_j(\text{weights}, B)` — the explicit form is :math:`\sum_i
-\lambda_i B_{j,i}` for SOS2 and :math:`B_{j,0} + \sum_i \delta_i (B_{j,i}
-- B_{j,i-1})` for incremental (see the method sections below).  Pinned
-tuples :math:`j` keep the equality, and the bounded tuple :math:`b` flips
-to the requested sign:
-
-.. math::
-
-   &e_j = W_j(\text{weights}, B)
-   \quad \text{(pinned, } j \ne b \text{)}
-
-   &e_b \ \text{sign}\ W_b(\text{weights}, B)
-   \quad \text{(bounded)}
-
-Internally this shows up as a stacked ``*_link`` equality covering the
-pinned tuples plus a separate signed ``*_output_link`` for the bounded
-tuple.  The ``method="lp"`` path encodes the same one-sided semantics
-without weights — see the LP section below.
-
-**Geometry.**  For 2 variables with ``sign="<="`` on a concave curve
-:math:`f`, the feasible region is the **hypograph** of :math:`f` on its
-domain:
-
-.. math::
-
-   \{ (x, y) \ :\ x_0 \le x \le x_n,\ y \le f(x) \}.
-
-For convex :math:`f` with ``sign=">="`` it is the **epigraph**.  Mismatched
-sign + curvature (convex + ``"<="``, or concave + ``">="``) describes a
-*non-convex* region — ``method="auto"`` falls back to SOS2/incremental
-and ``method="lp"`` raises.
-
-**Choice of bounded tuple and sign.**  Pick the sign matching the
-physically admissible direction for that expression:
-
-- ``"<="`` for a quantity with a controllable *dissipation path* — heat
-  rejection via cooling tower (*thermal curtailment*), electrical
-  curtailment, emissions after post-treatment — so undershooting the
-  curve is realisable.
-- ``">="`` for an *input* whose over-supply is admissible but wasteful —
-  fuel, raw materials — so overshooting the curve is realisable
-  (objective pressure then pulls the operating point onto the curve).
-
-The wrong direction (``"<="`` on fuel, ``">="`` on a non-curtailable
-output) yields a valid but **loose** formulation that admits operating
-points the plant cannot physically realise; an objective rewarding the
-wrong direction may then find a non-physical optimum — safe only when
-no such objective pressure exists.
-
-**When is a one-sided bound wanted?**
-
-For *continuous* curves, the main reason to reach for ``"<="`` / ``">="``
-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.  Inequality formulations also tighten the LP relaxation
-of SOS2/incremental, which can reduce branch-and-bound work even when
-LP itself is not applicable.
-
-For *disjunctive* curves (``segments(...)``), the per-tuple 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).
-
-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 link, which gives a valid but much
-more expensive model.  In that case prefer ``"=="`` unless you genuinely
-need the one-sided semantics.  See the
-:doc:`piecewise-inequality-bounds-tutorial` notebook for a full
-walkthrough.
-
-.. warning::
-
-   With a bounded tuple and ``active=0``, the output is only forced to
-   ``0`` on the signed side — the complementary bound still comes from
-   the output variable's own lower/upper 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
 -----------------------
 
+Three building blocks provide breakpoint data:
+
+- :func:`~linopy.breakpoints` — values along a single **connected**
+  curve.
+- :func:`~linopy.segments` — **disjoint** operating regions with gaps
+  between them (e.g. forbidden zones); selects the disjunctive
+  formulation automatically.
+- :class:`~linopy.Slopes` — per-piece slopes plus an initial ``y0``,
+  deferred until an x grid is supplied by a sibling tuple.
+
 From lists
 ~~~~~~~~~~
 
@@ -349,10 +195,151 @@ Link any number of variables through shared breakpoints (joint equality):
     )
 
 All variables are symmetric here; every feasible point is the same
-``λ``-weighted combination of breakpoints across all three.  With 3 or
-more tuples, only ``"=="`` signs are accepted — bounding one expression
-by a multi-input curve isn't supported yet; see the per-tuple sign
-section above for the issue link.
+``λ``-weighted combination of breakpoints across all three.  Sign
+restrictions apply (see *Per-tuple sign* below) — for ``N ≥ 3`` tuples
+all signs must be ``"=="``.
+
+
+Per-tuple sign — equality vs inequality
+----------------------------------------
+
+Roles and restrictions
+~~~~~~~~~~~~~~~~~~~~~~
+
+Each tuple's optional third element is a sign:
+
+- ``"=="`` (default) — **equality role**: the tuple enters as an
+  equality.
+- ``"<="`` / ``">="`` — **bounded**: the expression undershoots /
+  overshoots the curve.
+
+**Restrictions (current):**
+
+- At most one tuple may carry a non-equality sign — a single bounded side.
+- With **3 or more** tuples, all signs must be ``"=="``.
+
+Multi-bounded and N≥3-inequality use cases aren't supported yet.  If
+you have a concrete use case, please open an issue at
+https://github.com/PyPSA/linopy/issues so we can scope it properly.
+
+Geometry
+~~~~~~~~
+
+What the formulation actually constrains depends on the tuple count and
+signs:
+
+- **All-equality (default).**  Shared interpolation weights put the
+  joint :math:`(e_1, \ldots, e_N)` exactly on the curve.
+- **One bounded + one equality (2 tuples).**  The joint :math:`(x, y)`
+  lies in the **hypograph** (``"<="`` on a concave :math:`f`) or
+  **epigraph** (``">="`` on a convex :math:`f`):
+
+  .. math::
+
+     \{ (x, y) \ :\ x_{\min} \le x \le x_{\max},\ y \le f(x) \}
+     \qquad \text{(hypograph)}
+
+  The equality axis is just confined to its breakpoint domain
+  :math:`[x_{\min}, x_{\max}]` — a single coordinate can't locate a
+  curve point.  Same projection in LP (enforced directly) and
+  SOS2/incremental (enforced via the weight link).
+- **Mismatched sign + curvature** (convex + ``"<="``, or concave +
+  ``">="``) describes a *non-convex* region — ``method="auto"`` falls
+  back to SOS2/incremental and ``method="lp"`` raises.
+
+.. code-block:: python
+
+    # All-equality: joint (x, y) on the curve.
+    m.add_piecewise_formulation((y, y_pts), (x, x_pts))
+
+    # Bounded: joint (x, y) in the hypograph — y ≤ f(x), x ∈ [x_min, x_max].
+    m.add_piecewise_formulation((y, y_pts, "<="), (x, x_pts))
+
+    # Bounded: joint (x, y) in the epigraph — y ≥ f(x), x ∈ [x_min, x_max].
+    m.add_piecewise_formulation((y, y_pts, ">="), (x, x_pts))
+
+    # 3-variable all-equality (CHP): joint (power, fuel, heat) on the curve.
+    m.add_piecewise_formulation((power, p_pts), (fuel, f_pts), (heat, h_pts))
+
+Choice of bounded tuple and sign
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Pick the sign matching the physically admissible direction for that
+expression:
+
+- ``"<="`` for a quantity with a controllable *dissipation path* — heat
+  rejection via cooling tower (*thermal curtailment*), electrical
+  curtailment, emissions after post-treatment — so undershooting the
+  curve is realisable.
+- ``">="`` for an *input* whose over-supply is admissible but wasteful —
+  fuel, raw materials — so overshooting the curve is realisable
+  (objective pressure then pulls the operating point onto the curve).
+
+The wrong direction (``"<="`` on fuel, ``">="`` on a non-curtailable
+output) yields a valid but **loose** formulation that admits operating
+points the plant cannot physically realise; an objective rewarding the
+wrong direction may then find a non-physical optimum — safe only when
+no such objective pressure exists.
+
+When is a one-sided bound wanted?
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For *continuous* curves, the main reason to reach for ``"<="`` /
+``">="`` 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.  Inequality formulations also tighten
+the LP relaxation of SOS2/incremental, which can reduce branch-and-
+bound work even when LP itself is not applicable.
+
+For *disjunctive* curves (``segments(...)``), the per-tuple 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).
+
+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 link, which gives a valid but much
+more expensive model.  In that case prefer ``"=="`` unless you
+genuinely need the one-sided semantics.  See the
+:doc:`piecewise-inequality-bounds-tutorial` notebook for a full
+walkthrough.
+
+Formulation math
+~~~~~~~~~~~~~~~~
+
+For methods that introduce shared interpolation weights (SOS2 and
+incremental — see below), only the link constraint between the weights
+and the bounded expression changes.  Write the method-specific weighted
+sum of breakpoints for tuple :math:`j` as
+:math:`W_j(\text{weights}, B)` — the explicit form is
+:math:`\sum_i \lambda_i B_{j,i}` for SOS2 and
+:math:`B_{j,0} + \sum_i \delta_i (B_{j,i} - B_{j,i-1})` for incremental
+(see the method sections below).  Equality-signed tuples :math:`j` keep
+the equality, and the bounded tuple :math:`b` flips to the requested
+sign:
+
+.. math::
+
+   &e_j = W_j(\text{weights}, B)
+   \quad \text{(equality, } j \ne b \text{)}
+
+   &e_b \ \text{sign}\ W_b(\text{weights}, B)
+   \quad \text{(bounded)}
+
+Internally this shows up as a stacked ``*_link`` equality covering the
+equality-signed tuples plus a separate signed ``*_output_link`` for
+the bounded tuple.  The ``method="lp"`` path encodes the same one-sided
+semantics without weights — see the LP section below.
+
+.. warning::
+
+   With a bounded tuple and ``active=0``, the output is only forced to
+   ``0`` on the signed side — the complementary bound still comes from
+   the output variable's own lower/upper 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.
 
 
 Formulation Methods
@@ -368,11 +355,9 @@ formulation based on ``sign``, curvature and breakpoint layout:
 - **Disjunctive (segments)** → SOS2 applied per segment with binary
   segment selection (the disjunctive formulation in the table below).
 
-The resolved choice is exposed on the returned ``PiecewiseFormulation`` via
-``.method`` (and ``.convexity`` when well-defined).  Disjunctive
-formulations report ``method="sos2"`` even though their structure is the
-per-segment variant — the table below treats it as a separate column for
-clarity.  An ``INFO``-level log line explains the resolution whenever
+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.
 
 At-a-glance comparison:
@@ -427,12 +412,19 @@ At-a-glance comparison:
      - SOS2-capable (or MIP via :ref:`Big-M reformulation <sos-reformulation>`)
      - SOS2 + MIP (or MIP via :ref:`Big-M reformulation <sos-reformulation>`)
 
-LP (chord-line) Formulation
+.. note::
+
+   Disjunctive formulations report ``method="sos2"`` (the underlying
+   per-segment encoding uses SOS2), but the table treats them as a
+   separate column because the per-segment binaries change the
+   auxiliary-variable structure and solver requirements.
+
+LP (chord-line) formulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-For **2-variable inequality** on a **convex** or **concave** curve.  Adds one
-chord inequality per piece plus a domain bound — no auxiliary variables and
-no MIP relaxation:
+For **2-variable inequality** on a **convex** or **concave** curve.
+Adds one chord inequality per piece plus a domain bound — no auxiliary
+variables and no MIP relaxation:
 
 .. math::
 
@@ -442,18 +434,18 @@ no MIP relaxation:
    &x_{\min} \le x \le x_{\max}
 
 where :math:`m_k = (y_{k+1} - y_k)/(x_{k+1} - x_k)` and
-:math:`c_k = y_k - m_k\, x_k`.  The domain bound uses
-:math:`x_{\min}` and :math:`x_{\max}` rather than the first/last
-breakpoint so that descending x grids work too — strictly-monotonic
-breakpoints are accepted in either order.  For concave :math:`f` with
-``sign="<="``, the intersection of all chord inequalities equals the
-hypograph of :math:`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
-:doc:`piecewise-inequality-bounds-tutorial`).  ``method="auto"`` detects
-this and falls back; ``method="lp"`` raises.
+:math:`c_k = y_k - m_k\, x_k`.  The domain bound uses :math:`x_{\min}`
+and :math:`x_{\max}` rather than the first/last breakpoint so that
+descending x grids work too — strictly-monotonic breakpoints are
+accepted in either order.  For concave :math:`f` with ``sign="<="``,
+the intersection of all chord inequalities equals the hypograph of
+:math:`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 :doc:`piecewise-inequality-bounds-tutorial`).
+``method="auto"`` detects this and falls back; ``method="lp"`` raises.
 
 .. code-block:: python
 
@@ -463,18 +455,26 @@ this and falls back; ``method="lp"`` raises.
     # Or explicitly:
     m.add_piecewise_formulation((y, yp, "<="), (x, xp), method="lp")
 
-**Not supported with** ``method="lp"``: all-equality, more than 2 tuples,
-and ``active``.  ``method="auto"`` falls back to SOS2/incremental in all
-three cases.
+**Not supported with** ``method="lp"``: all-equality, 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-piece
-chord as a :class:`~linopy.expressions.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.
+Chord expressions as a building block
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The underlying chord expressions are also exposed as a standalone
+helper, :func:`~linopy.tangent_lines`, which returns the per-piece
+chord as a :class:`~linopy.expressions.LinearExpression` with no
+variables created.  Use it directly when you want to compose the chord
+bound with other constraints by hand, without the domain bound that
+``method="lp"`` adds automatically:
+
+.. code-block:: python
 
-Incremental (Delta) Formulation
+    chord = linopy.tangent_lines(x, x_pts, y_pts)
+    m.add_constraints(y <= chord + slack)
+
+Incremental (Delta) formulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The default MIP encoding when ``method="auto"`` is in play and breakpoints
@@ -491,8 +491,8 @@ binary indicators :math:`z_i`:
    &e_j = B_{j,0} + \sum_{i=1}^{n} \delta_i \, (B_{j,i} - B_{j,i-1})
 
 With a bounded tuple, the link to that tuple's expression flips to the
-requested sign while the pinned tuples keep the equality above (see
-the *Per-tuple sign* section's *Formulation* block).
+requested sign while the equality-signed tuples keep the equality above
+(see the *Formulation math* block in *Per-tuple sign*).
 
 .. code-block:: python
 
@@ -500,7 +500,7 @@ the *Per-tuple sign* section's *Formulation* block).
 
 **Limitation:** breakpoint sequences must be strictly monotonic.
 
-SOS2 (Convex Combination)
+SOS2 (Convex combination)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Fallback when breakpoints aren't strictly monotonic (the only case
@@ -536,7 +536,7 @@ above.
 
     m.add_piecewise_formulation((power, xp), (fuel, yp), method="sos2")
 
-Disjunctive (Disaggregated Convex Combination)
+Disjunctive (Disaggregated convex combination)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 For **disconnected segments** (gaps between operating regions).  Binary
@@ -583,15 +583,9 @@ are forced to zero:
 - ``commit=1``: power operates in [30, 100], fuel = f(power)
 - ``commit=0``: power = 0, fuel = 0
 
-Not supported with ``method="lp"``.
-
-.. note::
-
-   With a bounded tuple, 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 per-tuple sign section above for details.
+Not supported with ``method="lp"``.  For bounded-tuple semantics under
+``active=0`` see the warning in *Per-tuple sign — Formulation math*
+above.
 
 Auto-broadcasting
 ~~~~~~~~~~~~~~~~~
diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 0c47ffae..260ff56e 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -58,7 +58,9 @@
    "source": [
     "## Setup — a convex heat-rate curve\n",
     "\n",
-    "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."
+    "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.\n",
+    "\n",
+    "The breakpoint arrays:"
    ]
   },
   {
@@ -73,13 +75,7 @@
    "outputs": [],
    "source": [
     "power_pts = np.array([0.0, 30.0, 60.0, 100.0])\n",
-    "fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])  # slopes 1.2, 1.6, 2.15 (convex)\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(5, 4))\n",
-    "ax.plot(power_pts, fuel_pts, \"o-\", color=\"C0\", lw=2)\n",
-    "ax.set(xlabel=\"power\", ylabel=\"fuel\", title=\"Convex heat-rate curve f(power)\")\n",
-    "ax.grid(alpha=0.3)\n",
-    "plt.tight_layout()"
+    "fuel_pts = np.array([0.0, 36.0, 84.0, 170.0])  # slopes 1.2, 1.6, 2.15 (convex)"
    ]
   },
   {
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 3648411e..d93c54e3 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -123,15 +123,11 @@
    "source": [
     "## 2. Picking a method\n",
     "\n",
-    "`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.\n",
+    "`method=\"auto\"` (default) picks the cheapest correct formulation based on sign, curvature and breakpoint layout. The explicit options are `\"sos2\"`, `\"incremental\"`, and `\"lp\"`; the choice is about **cost** (auxiliary variables, solver capability), not correctness — on cases where they all apply they give the same optimum.\n",
     "\n",
-    "| method | needs | creates |\n",
-    "|---|---|---|\n",
-    "| `sos2` | SOS2-capable solver | lambdas (continuous) |\n",
-    "| `incremental` | MIP solver, strictly monotonic breakpoints | deltas (continuous) + binaries |\n",
-    "| `lp` | any LP solver | no variables — requires `sign != \"==\"`, 2 tuples, matching curvature |\n",
+    "For now: a quick sanity check that all applicable methods yield the same fuel dispatch on the convex curve from §1.\n",
     "\n",
-    "Below: all applicable methods yield the same fuel dispatch on this convex curve."
+    "A full comparison — when each method dispatches, what sign/curvature/breakpoint patterns each requires — lives in §3 (disjunctive), §4 (inequalities), and the [reference page's \"Formulation Methods\" section](piecewise-linear-constraints)."
    ]
   },
   {

From e706734972dfbb0f6694c4112e1ca8ab29e01099 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 16:11:26 +0200
Subject: [PATCH 15/29] docs(piecewise): drop "Formulation math" subsection,
 fold warning into Active
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The W_j(weights, B) abstraction it introduced was immediately re-
derived in each method section below (SOS2 with λ, incremental with
δ), so it added a layer for the reader to climb through without
delivering content the method sections didn't already.  The
"*_link" / "*_output_link" naming detail contradicted our own guidance
elsewhere ("exact name suffixes are an implementation detail and may
evolve").  The "equality keeps equality, bounded flips sign" idea is
already conveyed by Roles & restrictions and Geometry above.

The one load-bearing nugget was the active=0 warning.  Moved it into
the Active parameter (unit commitment) subsection under Advanced
Features, where readers setting `active=...` actually land — that's
where the gotcha needs to be visible.

Also dropped the now-stale cross-reference from the Incremental method
section.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 53 ++++++----------------------
 1 file changed, 11 insertions(+), 42 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 6e558772..d4163b8b 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -305,43 +305,6 @@ genuinely need the one-sided semantics.  See the
 :doc:`piecewise-inequality-bounds-tutorial` notebook for a full
 walkthrough.
 
-Formulation math
-~~~~~~~~~~~~~~~~
-
-For methods that introduce shared interpolation weights (SOS2 and
-incremental — see below), only the link constraint between the weights
-and the bounded expression changes.  Write the method-specific weighted
-sum of breakpoints for tuple :math:`j` as
-:math:`W_j(\text{weights}, B)` — the explicit form is
-:math:`\sum_i \lambda_i B_{j,i}` for SOS2 and
-:math:`B_{j,0} + \sum_i \delta_i (B_{j,i} - B_{j,i-1})` for incremental
-(see the method sections below).  Equality-signed tuples :math:`j` keep
-the equality, and the bounded tuple :math:`b` flips to the requested
-sign:
-
-.. math::
-
-   &e_j = W_j(\text{weights}, B)
-   \quad \text{(equality, } j \ne b \text{)}
-
-   &e_b \ \text{sign}\ W_b(\text{weights}, B)
-   \quad \text{(bounded)}
-
-Internally this shows up as a stacked ``*_link`` equality covering the
-equality-signed tuples plus a separate signed ``*_output_link`` for
-the bounded tuple.  The ``method="lp"`` path encodes the same one-sided
-semantics without weights — see the LP section below.
-
-.. warning::
-
-   With a bounded tuple and ``active=0``, the output is only forced to
-   ``0`` on the signed side — the complementary bound still comes from
-   the output variable's own lower/upper 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.
-
-
 Formulation Methods
 -------------------
 
@@ -491,8 +454,7 @@ binary indicators :math:`z_i`:
    &e_j = B_{j,0} + \sum_{i=1}^{n} \delta_i \, (B_{j,i} - B_{j,i-1})
 
 With a bounded tuple, the link to that tuple's expression flips to the
-requested sign while the equality-signed tuples keep the equality above
-(see the *Formulation math* block in *Per-tuple sign*).
+requested sign while the equality-signed tuples keep the equality above.
 
 .. code-block:: python
 
@@ -583,9 +545,16 @@ are forced to zero:
 - ``commit=1``: power operates in [30, 100], fuel = f(power)
 - ``commit=0``: power = 0, fuel = 0
 
-Not supported with ``method="lp"``.  For bounded-tuple semantics under
-``active=0`` see the warning in *Per-tuple sign — Formulation math*
-above.
+Not supported with ``method="lp"``.
+
+.. warning::
+
+   With a bounded tuple and ``active=0``, the output is only forced to
+   ``0`` on the signed side — the complementary bound still comes from
+   the output variable's own lower/upper 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.
 
 Auto-broadcasting
 ~~~~~~~~~~~~~~~~~

From c2fe7f6ce8d4b4becec2d65b14149c225a18c1d5 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 16:15:33 +0200
Subject: [PATCH 16/29] docs(piecewise): explain why active is not supported
 with method="lp"

"Not supported with method='lp'" reads as a missing feature; one-line
addition explains it's structural (gating needs a binary) and points
to the auto-dispatch and to tangent_lines for manual gating.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index d4163b8b..8bfd669d 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -545,7 +545,9 @@ are forced to zero:
 - ``commit=1``: power operates in [30, 100], fuel = f(power)
 - ``commit=0``: power = 0, fuel = 0
 
-Not supported with ``method="lp"``.
+Not supported with ``method="lp"`` (gating needs a binary).  Use
+``method="auto"``, or *Chord expressions as a building block* for
+manual gating.
 
 .. warning::
 

From 04fccec40f159478baf7213b14d58ad0b722c73a Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 16:21:06 +0200
Subject: [PATCH 17/29] ci: trigger docs build


From 644bb50442207561c2aaeb2231e110b56aa6c59d Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:37:35 +0200
Subject: [PATCH 18/29] =?UTF-8?q?docs(piecewise):=20replace=20=C2=A73=20(0?=
 =?UTF-8?q?,0)=20"off"=20segment=20with=20discrete-sizing=20example?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The previous example used `segments([(0, 0), (50, 80)])` — the (0,0)
"off" segment was a hack to encode on/off behaviour, which is exactly
what `active=...` (§5) is for.  Using disjunctive for on/off
duplicates §5 and teaches a bad pattern.

New example shows what disjunctive is actually for: equipment without
a continuous design space — gas turbines in three commercial classes
(small / medium / large), each with its own non-overlapping power band
and 2-piece heat-rate curve.  The formulation picks the appropriate
class per timestep based on demand.

Added a forward pointer to §5 for the on/off case.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 30 ++++++++++++---------
 1 file changed, 18 insertions(+), 12 deletions(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index d93c54e3..d999ad42 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -159,11 +159,15 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 3. Disjunctive segments — gaps in the operating range\n",
+    "## 3. Disjunctive segments — choose one of several operating modes\n",
     "\n",
-    "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.\n",
+    "Some equipment doesn't have a continuous design space.  Gas turbines come in commercial classes (small / medium / large), pumps in discrete VSD steps, batteries in catalogued module sizes — each option has its own operating range and performance curve, and the model has to pick one.\n",
     "\n",
-    "Below the first segment is `(0, 0)` — a degenerate \"off\" state where both endpoints sit at the origin, so the unit produces no power and incurs no cost. The second segment is the active range, with cost rising from 125 to 200 over 50–80 MW."
+    "`segments()` models this directly: one segment per option.  A binary picks exactly one; SOS2 interpolates within it as usual.\n",
+    "\n",
+    "Below: a gas turbine offered in three classes, with non-overlapping power ranges and class-specific heat-rate curves.  The formulation picks the appropriate class per timestep based on demand.\n",
+    "\n",
+    "(For a single on/off gate on one curve — \"the unit is either off or running in its operating range\" — use `active=...` instead; see §5.)"
    ]
   },
   {
@@ -179,25 +183,27 @@
    "outputs": [],
    "source": [
     "m = linopy.Model()\n",
-    "power = m.add_variables(name=\"power\", lower=0, upper=80, coords=[time])\n",
-    "cost = m.add_variables(name=\"cost\", lower=0, coords=[time])\n",
-    "backup = m.add_variables(name=\"backup\", lower=0, coords=[time])\n",
+    "power = m.add_variables(name=\"power\", lower=0, upper=200, coords=[time])\n",
+    "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
     "\n",
+    "# Three commercially-available turbine classes, each with a 2-piece\n",
+    "# heat-rate curve over its operating band.\n",
     "m.add_piecewise_formulation(\n",
-    "    (power, linopy.segments([(0, 0), (50, 80)])),  # two disjoint segments\n",
-    "    (cost, linopy.segments([(0, 0), (125, 200)])),\n",
+    "    (power, linopy.segments([(10, 20, 30), (40, 60, 80), (100, 150, 200)])),\n",
+    "    (fuel, linopy.segments([(12, 22, 35), (40, 56, 80), (105, 150, 215)])),\n",
     ")\n",
-    "m.add_constraints(power + backup == xr.DataArray([15, 60, 75], coords=[time]))\n",
-    "m.add_objective(cost.sum() + 10 * backup.sum())\n",
+    "demand = xr.DataArray([25, 60, 150], coords=[time])\n",
+    "m.add_constraints(power == demand)\n",
+    "m.add_objective(fuel.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
-    "m.solution[[\"power\", \"cost\", \"backup\"]].to_pandas()"
+    "m.solution[[\"power\", \"fuel\"]].to_pandas()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "At *t=1* the 15 MW demand falls in the forbidden zone; the unit sits at 0 and backup fills the gap."
+    "At each timestep the formulation picks the unique class whose band contains the demand: small (10–30 MW) at *t=1*, medium (40–80 MW) at *t=2*, large (100–200 MW) at *t=3*.  Fuel is interpolated on the chosen class's heat-rate curve."
    ]
   },
   {

From 7935b4c5dcfdc34f76d41f140a77a72721ad9181 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:38:14 +0200
Subject: [PATCH 19/29] docs(piecewise): replace (0,0) segment in rst example
 with discrete-sizing
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Matches the tutorial §3 rewrite — uses three turbine classes with
non-overlapping operating bands as the canonical disjunctive example,
and adds a one-line steer away from the (0,0) anti-pattern toward
`active=...` for on/off gating.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 8bfd669d..31d2eace 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -169,17 +169,25 @@ over remaining dimensions (e.g. ``time``).
 Disjunctive segments
 ~~~~~~~~~~~~~~~~~~~~
 
-For disconnected operating regions (e.g. forbidden zones), use ``segments()``:
+For equipment with disconnected operating regions — discrete commercial
+sizes, multi-speed pumps, forbidden zones — use ``segments()``.  Each
+segment is one option's (range, curve); a binary picks exactly one.
 
 .. code-block:: python
 
+    # Three turbine classes with non-overlapping operating bands.
     m.add_piecewise_formulation(
-        (power, linopy.segments([(0, 0), (50, 80)])),
-        (cost, linopy.segments([(0, 0), (125, 200)])),
+        (power, linopy.segments([(10, 20, 30), (40, 60, 80), (100, 150, 200)])),
+        (fuel, linopy.segments([(12, 22, 35), (40, 56, 80), (105, 150, 215)])),
     )
 
-The disjunctive formulation is selected automatically when breakpoints have a
-segment dimension.  A bounded tuple (``"<="`` / ``">="``) also works here.
+The disjunctive formulation is selected automatically when breakpoints
+have a segment dimension.  A bounded tuple (``"<="`` / ``">="``) also
+works here.
+
+For a single on/off gate on one continuous curve, prefer ``active=...``
+(see *Advanced Features* below) — using a degenerate ``(0, 0)`` segment
+to encode "off" mixes the disjunctive concept with on/off logic.
 
 N-variable linking
 ~~~~~~~~~~~~~~~~~~

From bf32c46825fe6f5b78f42bafce92dbfb657b63bc Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:46:31 +0200
Subject: [PATCH 20/29] docs(piecewise): use pump VSD as the disjunctive
 example
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The turbine-classes framing was claiming "investment" but the
formulation puts the binary on the time dimension by default — so
the model was actually re-deciding the class each timestep, not
making a one-time sizing call.  Either we'd need a contrived
constraint to make the binary time-invariant, or we'd be using an
"investment" narrative on per-period selection.

Pump VSD is honest per-period selection: a pump with three stepped
speed settings, each covering a different flow band, switches speed
between dispatch periods as a matter of normal operation.  No
constraints to bolt on, no narrative-vs-math gap.

Also added "switchable combustion cycles" and "allowed bands around
forbidden vibration zones" as alternative real use cases in the
intro, so readers see disjunctive is for genuine multi-mode
equipment.

Updated both the rst Disjunctive-segments code block and the
tutorial §3 (markdown intro, code, explanatory cell) to use the
pump example.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        | 13 ++++-----
 examples/piecewise-linear-constraints.ipynb | 30 ++++++++++-----------
 2 files changed, 22 insertions(+), 21 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 31d2eace..6ad98798 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -169,16 +169,17 @@ over remaining dimensions (e.g. ``time``).
 Disjunctive segments
 ~~~~~~~~~~~~~~~~~~~~
 
-For equipment with disconnected operating regions — discrete commercial
-sizes, multi-speed pumps, forbidden zones — use ``segments()``.  Each
-segment is one option's (range, curve); a binary picks exactly one.
+For equipment with disconnected operating modes — stepped pump speeds,
+switchable combustion cycles, allowed bands around forbidden vibration
+zones — use ``segments()``.  Each segment is one mode's (range, curve);
+a binary picks exactly one per operating point.
 
 .. code-block:: python
 
-    # Three turbine classes with non-overlapping operating bands.
+    # Pump VSD with three stepped speed settings.
     m.add_piecewise_formulation(
-        (power, linopy.segments([(10, 20, 30), (40, 60, 80), (100, 150, 200)])),
-        (fuel, linopy.segments([(12, 22, 35), (40, 56, 80), (105, 150, 215)])),
+        (flow, linopy.segments([(5, 15, 25), (30, 50, 70), (80, 115, 150)])),
+        (power, linopy.segments([(1, 3.5, 7), (10, 19, 32), (40, 70, 115)])),
     )
 
 The disjunctive formulation is selected automatically when breakpoints
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index d999ad42..e872f93d 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -161,13 +161,13 @@
    "source": [
     "## 3. Disjunctive segments — choose one of several operating modes\n",
     "\n",
-    "Some equipment doesn't have a continuous design space.  Gas turbines come in commercial classes (small / medium / large), pumps in discrete VSD steps, batteries in catalogued module sizes — each option has its own operating range and performance curve, and the model has to pick one.\n",
+    "Some equipment has **discrete operating modes** rather than a continuous design space.  A pump with a stepped variable-speed drive (VSD) has three speed settings, each covering a different flow band with its own flow-vs-power curve; a gas turbine might switch between open-cycle and combined-cycle; a synchronous machine has allowed operating bands around forbidden vibration zones.  At each operating point the model picks one mode.\n",
     "\n",
-    "`segments()` models this directly: one segment per option.  A binary picks exactly one; SOS2 interpolates within it as usual.\n",
+    "`segments()` models this directly: one segment per mode, a binary picks exactly one, and SOS2 interpolates within it.\n",
     "\n",
-    "Below: a gas turbine offered in three classes, with non-overlapping power ranges and class-specific heat-rate curves.  The formulation picks the appropriate class per timestep based on demand.\n",
+    "Below: a pump with three stepped speed settings, each with its own flow band and a flow-vs-power curve (power rises faster than linearly with flow).  Per timestep the formulation picks the speed that covers the required flow.\n",
     "\n",
-    "(For a single on/off gate on one curve — \"the unit is either off or running in its operating range\" — use `active=...` instead; see §5.)"
+    "(For an on/off gate on a single continuous curve — \"the unit is either off or running\" — use `active=...` instead; see §5.)"
    ]
   },
   {
@@ -183,27 +183,27 @@
    "outputs": [],
    "source": [
     "m = linopy.Model()\n",
-    "power = m.add_variables(name=\"power\", lower=0, upper=200, coords=[time])\n",
-    "fuel = m.add_variables(name=\"fuel\", lower=0, coords=[time])\n",
+    "flow = m.add_variables(name=\"flow\", lower=0, upper=150, coords=[time])\n",
+    "power = m.add_variables(name=\"power\", lower=0, coords=[time])\n",
     "\n",
-    "# Three commercially-available turbine classes, each with a 2-piece\n",
-    "# heat-rate curve over its operating band.\n",
+    "# Pump VSD with three stepped speed settings.  Each step covers a\n",
+    "# different flow band and has a mildly convex flow-vs-power curve.\n",
     "m.add_piecewise_formulation(\n",
-    "    (power, linopy.segments([(10, 20, 30), (40, 60, 80), (100, 150, 200)])),\n",
-    "    (fuel, linopy.segments([(12, 22, 35), (40, 56, 80), (105, 150, 215)])),\n",
+    "    (flow, linopy.segments([(5, 15, 25), (30, 50, 70), (80, 115, 150)])),\n",
+    "    (power, linopy.segments([(1, 3.5, 7), (10, 19, 32), (40, 70, 115)])),\n",
     ")\n",
-    "demand = xr.DataArray([25, 60, 150], coords=[time])\n",
-    "m.add_constraints(power == demand)\n",
-    "m.add_objective(fuel.sum())\n",
+    "demand = xr.DataArray([15, 50, 115], coords=[time])\n",
+    "m.add_constraints(flow == demand)\n",
+    "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
-    "m.solution[[\"power\", \"fuel\"]].to_pandas()"
+    "m.solution[[\"flow\", \"power\"]].to_pandas()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "At each timestep the formulation picks the unique class whose band contains the demand: small (10–30 MW) at *t=1*, medium (40–80 MW) at *t=2*, large (100–200 MW) at *t=3*.  Fuel is interpolated on the chosen class's heat-rate curve."
+    "At each timestep the formulation picks the speed setting whose flow band covers demand: low (5–25 m³/h) at *t=1*, medium (30–70) at *t=2*, high (80–150) at *t=3*.  Power is interpolated on the chosen setting's curve."
    ]
   },
   {

From 6438db037901966feeef5e2d0effa58b70794142 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:51:39 +0200
Subject: [PATCH 21/29] docs(piecewise): two pumps with two bands each, demand
 catches the gap
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Two pumps in parallel, each with two operating bands (low 5–25 m³/h,
high 40–100 m³/h) and a forbidden zone in between.  Demand profile
[30, 75, 150] makes every timestep single-pump-infeasible:

- t=1, demand=30: lands in the single-pump gap (25, 40); both pumps
  run in low band, splitting the load.
- t=2, demand=75: too much for low+low (max 50), too little for
  high+high (min 80); the low pump tops out at 25, the high pump
  covers the remaining 50.
- t=3, demand=150: exceeds a single pump's maximum (100); both pumps
  run in high band.

Each segment is now a single piece (2 breakpoints), so the example is
clean and the load-split arithmetic is explicit.  Rst snippet updated
to match.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst        |  11 +-
 examples/piecewise-linear-constraints.ipynb | 111 ++++++++++----------
 2 files changed, 63 insertions(+), 59 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 6ad98798..4a537634 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -169,17 +169,18 @@ over remaining dimensions (e.g. ``time``).
 Disjunctive segments
 ~~~~~~~~~~~~~~~~~~~~
 
-For equipment with disconnected operating modes — stepped pump speeds,
+For equipment with disconnected operating bands — stepped pump speeds,
 switchable combustion cycles, allowed bands around forbidden vibration
-zones — use ``segments()``.  Each segment is one mode's (range, curve);
+zones — use ``segments()``.  Each segment is one band's (range, curve);
 a binary picks exactly one per operating point.
 
 .. code-block:: python
 
-    # Pump VSD with three stepped speed settings.
+    # Stepped pump with two speed bands; the gap between them is a
+    # forbidden zone.
     m.add_piecewise_formulation(
-        (flow, linopy.segments([(5, 15, 25), (30, 50, 70), (80, 115, 150)])),
-        (power, linopy.segments([(1, 3.5, 7), (10, 19, 32), (40, 70, 115)])),
+        (flow, linopy.segments([(5, 25), (40, 100)])),
+        (power, linopy.segments([(1, 7), (15, 50)])),
     )
 
 The disjunctive formulation is selected automatically when breakpoints
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index e872f93d..66c69755 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -31,14 +31,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:48.865977Z",
-     "start_time": "2026-05-11T12:50:48.861703Z"
+     "end_time": "2026-05-11T17:47:02.640180Z",
+     "start_time": "2026-05-11T17:47:02.630870Z"
     }
    },
-   "outputs": [],
    "source": [
     "import warnings\n",
     "\n",
@@ -64,7 +62,9 @@
     "    ax.set(xlabel=xlabel, ylabel=ylabel)\n",
     "    ax.legend()\n",
     "    return ax"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -77,14 +77,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:51:46.043507Z",
-     "start_time": "2026-05-11T12:51:45.930380Z"
+     "end_time": "2026-05-11T17:47:02.757350Z",
+     "start_time": "2026-05-11T17:47:02.645672Z"
     }
    },
-   "outputs": [],
    "source": [
     "demand = xr.DataArray([50, 80, 30], coords=[time])\n",
     "\n",
@@ -101,21 +99,23 @@
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.027645Z",
-     "start_time": "2026-05-11T12:50:48.982161Z"
+     "end_time": "2026-05-11T17:47:02.817859Z",
+     "start_time": "2026-05-11T17:47:02.763695Z"
     }
    },
-   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -132,14 +132,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.330961Z",
-     "start_time": "2026-05-11T12:50:49.039795Z"
+     "end_time": "2026-05-11T17:47:03.115144Z",
+     "start_time": "2026-05-11T17:47:02.832347Z"
     }
    },
-   "outputs": [],
    "source": [
     "def solve_method(method):\n",
     "    m = linopy.Model()\n",
@@ -153,57 +151,62 @@
     "\n",
     "\n",
     "pd.DataFrame({m: solve_method(m) for m in [\"auto\", \"sos2\", \"incremental\"]})"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 3. Disjunctive segments — choose one of several operating modes\n",
+    "## 3. Disjunctive segments — discrete operating bands\n",
     "\n",
-    "Some equipment has **discrete operating modes** rather than a continuous design space.  A pump with a stepped variable-speed drive (VSD) has three speed settings, each covering a different flow band with its own flow-vs-power curve; a gas turbine might switch between open-cycle and combined-cycle; a synchronous machine has allowed operating bands around forbidden vibration zones.  At each operating point the model picks one mode.\n",
+    "Some equipment has **disjoint operating ranges** rather than a continuous one.  A stepped pump has two speed bands with a forbidden zone between them — the pump physically can't operate in that gap.  `segments()` models this directly: one segment per band, a binary picks exactly one per operating point.\n",
     "\n",
-    "`segments()` models this directly: one segment per mode, a binary picks exactly one, and SOS2 interpolates within it.\n",
+    "Below: two pumps in parallel, each with a low band (5–25 m³/h) and a high band (40–100 m³/h).  Demands that land in the single-pump gap or above its maximum force the optimiser to combine bands across the two pumps.\n",
     "\n",
-    "Below: a pump with three stepped speed settings, each with its own flow band and a flow-vs-power curve (power rises faster than linearly with flow).  Per timestep the formulation picks the speed that covers the required flow.\n",
-    "\n",
-    "(For an on/off gate on a single continuous curve — \"the unit is either off or running\" — use `active=...` instead; see §5.)"
+    "(For an on/off gate on a single continuous curve, use `active=...` instead; see §5.)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.487115Z",
-     "start_time": "2026-05-11T12:50:49.347614Z"
-    },
-    "scrolled": true
+     "end_time": "2026-05-11T17:47:03.268700Z",
+     "start_time": "2026-05-11T17:47:03.129554Z"
+    }
    },
-   "outputs": [],
    "source": [
+    "pumps = pd.Index([\"p1\", \"p2\"], name=\"pump\")\n",
+    "\n",
     "m = linopy.Model()\n",
-    "flow = m.add_variables(name=\"flow\", lower=0, upper=150, coords=[time])\n",
-    "power = m.add_variables(name=\"power\", lower=0, coords=[time])\n",
+    "flow = m.add_variables(name=\"flow\", lower=0, upper=100, coords=[pumps, time])\n",
+    "power = m.add_variables(name=\"power\", lower=0, coords=[pumps, time])\n",
     "\n",
-    "# Pump VSD with three stepped speed settings.  Each step covers a\n",
-    "# different flow band and has a mildly convex flow-vs-power curve.\n",
+    "# Each pump has two operating bands; the gap between them is a forbidden zone.\n",
     "m.add_piecewise_formulation(\n",
-    "    (flow, linopy.segments([(5, 15, 25), (30, 50, 70), (80, 115, 150)])),\n",
-    "    (power, linopy.segments([(1, 3.5, 7), (10, 19, 32), (40, 70, 115)])),\n",
+    "    (flow, linopy.segments([(5, 25), (40, 100)])),\n",
+    "    (power, linopy.segments([(1, 7), (15, 50)])),\n",
     ")\n",
-    "demand = xr.DataArray([15, 50, 115], coords=[time])\n",
-    "m.add_constraints(flow == demand)\n",
+    "demand = xr.DataArray([30, 75, 150], coords=[time])\n",
+    "m.add_constraints(flow.sum(\"pump\") == demand)\n",
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"flow\", \"power\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "At each timestep the formulation picks the speed setting whose flow band covers demand: low (5–25 m³/h) at *t=1*, medium (30–70) at *t=2*, high (80–150) at *t=3*.  Power is interpolated on the chosen setting's curve."
+    "Every timestep is single-pump-infeasible:\n",
+    "\n",
+    "- *t=1*, demand=30: in the single-pump gap (25, 40).  Both pumps run in low band, splitting the load.\n",
+    "- *t=2*, demand=75: too much for low+low (max 50), too little for high+high (min 80).  The low pump tops out at 25 m³/h; the high pump covers the remaining 50.\n",
+    "- *t=3*, demand=150: above a single pump's maximum (100).  Both pumps run in high band."
    ]
   },
   {
@@ -221,14 +224,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.556562Z",
-     "start_time": "2026-05-11T12:50:49.492505Z"
+     "end_time": "2026-05-11T17:47:03.568401Z",
+     "start_time": "2026-05-11T17:47:03.284318Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -245,14 +246,16 @@
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.603928Z",
+     "end_time": "2026-05-11T17:47:03.583650Z",
      "start_time": "2026-05-11T12:50:49.566231Z"
     }
    },
@@ -283,7 +286,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.746059Z",
+     "end_time": "2026-05-11T17:47:03.584848Z",
      "start_time": "2026-05-11T12:50:49.613733Z"
     }
    },
@@ -316,7 +319,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.792782Z",
+     "end_time": "2026-05-11T17:47:03.585384Z",
      "start_time": "2026-05-11T12:50:49.751516Z"
     }
    },
@@ -339,7 +342,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.914868Z",
+     "end_time": "2026-05-11T17:47:03.585496Z",
      "start_time": "2026-05-11T12:50:49.802514Z"
     }
    },
@@ -369,7 +372,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:49.993753Z",
+     "end_time": "2026-05-11T17:47:03.590541Z",
      "start_time": "2026-05-11T12:50:49.922140Z"
     }
    },
@@ -403,7 +406,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:50.141147Z",
+     "end_time": "2026-05-11T17:47:03.590784Z",
      "start_time": "2026-05-11T12:50:50.003565Z"
     }
    },
@@ -441,7 +444,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T12:50:50.244428Z",
+     "end_time": "2026-05-11T17:47:03.590905Z",
      "start_time": "2026-05-11T12:50:50.151699Z"
     }
    },

From c0c4cb4fa572787d51368c0f27f6a8baed965654 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:53:57 +0200
Subject: [PATCH 22/29] =?UTF-8?q?docs(piecewise):=20flatten=20=C2=A73=20ou?=
 =?UTF-8?q?tput=20to=20a=203-column=20DataFrame?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

power_p1 / power_p2 / flow (total) per timestep — the asymmetric split
at t=2 is visible directly in the powers (one pump at 7 kW, the other
at ~21 kW) instead of buried in a multi-index DataFrame.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 85 +++++++++++----------
 1 file changed, 46 insertions(+), 39 deletions(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 66c69755..98f7596a 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -31,12 +31,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:02.640180Z",
-     "start_time": "2026-05-11T17:47:02.630870Z"
+     "end_time": "2026-05-11T17:51:41.226105Z",
+     "start_time": "2026-05-11T17:51:41.221083Z"
     }
    },
+   "outputs": [],
    "source": [
     "import warnings\n",
     "\n",
@@ -62,9 +64,7 @@
     "    ax.set(xlabel=xlabel, ylabel=ylabel)\n",
     "    ax.legend()\n",
     "    return ax"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -77,12 +77,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:02.757350Z",
-     "start_time": "2026-05-11T17:47:02.645672Z"
+     "end_time": "2026-05-11T17:51:41.334560Z",
+     "start_time": "2026-05-11T17:51:41.229084Z"
     }
    },
+   "outputs": [],
    "source": [
     "demand = xr.DataArray([50, 80, 30], coords=[time])\n",
     "\n",
@@ -99,23 +101,21 @@
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:02.817859Z",
-     "start_time": "2026-05-11T17:47:02.763695Z"
+     "end_time": "2026-05-11T17:51:41.388313Z",
+     "start_time": "2026-05-11T17:51:41.338111Z"
     }
    },
+   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -132,12 +132,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.115144Z",
-     "start_time": "2026-05-11T17:47:02.832347Z"
+     "end_time": "2026-05-11T17:51:41.686866Z",
+     "start_time": "2026-05-11T17:51:41.401345Z"
     }
    },
+   "outputs": [],
    "source": [
     "def solve_method(method):\n",
     "    m = linopy.Model()\n",
@@ -151,9 +153,7 @@
     "\n",
     "\n",
     "pd.DataFrame({m: solve_method(m) for m in [\"auto\", \"sos2\", \"incremental\"]})"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
@@ -170,13 +170,15 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
-    "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.268700Z",
-     "start_time": "2026-05-11T17:47:03.129554Z"
-    }
+     "end_time": "2026-05-11T17:52:01.049181Z",
+     "start_time": "2026-05-11T17:52:00.881744Z"
+    },
+    "scrolled": true
    },
+   "outputs": [],
    "source": [
     "pumps = pd.Index([\"p1\", \"p2\"], name=\"pump\")\n",
     "\n",
@@ -193,10 +195,15 @@
     "m.add_constraints(flow.sum(\"pump\") == demand)\n",
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
-    "m.solution[[\"flow\", \"power\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+    "\n",
+    "pd.DataFrame(\n",
+    "    {\n",
+    "        \"power_p1\": m.solution[\"power\"].sel(pump=\"p1\").to_pandas(),\n",
+    "        \"power_p2\": m.solution[\"power\"].sel(pump=\"p2\").to_pandas(),\n",
+    "        \"flow\": m.solution[\"flow\"].sum(\"pump\").to_pandas(),\n",
+    "    }\n",
+    ")"
+   ]
   },
   {
    "cell_type": "markdown",
@@ -224,12 +231,14 @@
   },
   {
    "cell_type": "code",
+   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.568401Z",
+     "end_time": "2026-05-11T17:51:41.867742Z",
      "start_time": "2026-05-11T17:47:03.284318Z"
     }
    },
+   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -246,16 +255,14 @@
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ],
-   "outputs": [],
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.583650Z",
+     "end_time": "2026-05-11T17:51:41.868298Z",
      "start_time": "2026-05-11T12:50:49.566231Z"
     }
    },
@@ -286,7 +293,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.584848Z",
+     "end_time": "2026-05-11T17:51:41.874712Z",
      "start_time": "2026-05-11T12:50:49.613733Z"
     }
    },
@@ -319,7 +326,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.585384Z",
+     "end_time": "2026-05-11T17:51:41.875307Z",
      "start_time": "2026-05-11T12:50:49.751516Z"
     }
    },
@@ -342,7 +349,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.585496Z",
+     "end_time": "2026-05-11T17:51:41.875528Z",
      "start_time": "2026-05-11T12:50:49.802514Z"
     }
    },
@@ -372,7 +379,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.590541Z",
+     "end_time": "2026-05-11T17:51:41.877675Z",
      "start_time": "2026-05-11T12:50:49.922140Z"
     }
    },
@@ -406,7 +413,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.590784Z",
+     "end_time": "2026-05-11T17:51:41.881102Z",
      "start_time": "2026-05-11T12:50:50.003565Z"
     }
    },
@@ -444,7 +451,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:47:03.590905Z",
+     "end_time": "2026-05-11T17:51:41.882210Z",
      "start_time": "2026-05-11T12:50:50.151699Z"
     }
    },

From 262895241bc0b6431202b5d374415acb6c8635b6 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 19:58:01 +0200
Subject: [PATCH 23/29] =?UTF-8?q?docs(piecewise):=20show=20per-pump=20flow?=
 =?UTF-8?q?=20and=20power=20in=20=C2=A73=20output?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

flat 4-column DataFrame (flow_p1, flow_p2, power_p1, power_p2) per
timestep, via unstack on the pump dimension.  The asymmetric splits
are visible directly: at t=2 one pump runs at flow=25 (low band), the
other at flow=50 (high band).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 97 ++++++++++-----------
 1 file changed, 47 insertions(+), 50 deletions(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 98f7596a..d374d540 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -31,14 +31,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.226105Z",
-     "start_time": "2026-05-11T17:51:41.221083Z"
+     "end_time": "2026-05-11T17:55:29.339781Z",
+     "start_time": "2026-05-11T17:55:29.335906Z"
     }
    },
-   "outputs": [],
    "source": [
     "import warnings\n",
     "\n",
@@ -64,7 +62,9 @@
     "    ax.set(xlabel=xlabel, ylabel=ylabel)\n",
     "    ax.legend()\n",
     "    return ax"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -77,14 +77,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.334560Z",
-     "start_time": "2026-05-11T17:51:41.229084Z"
+     "end_time": "2026-05-11T17:55:29.448389Z",
+     "start_time": "2026-05-11T17:55:29.351171Z"
     }
    },
-   "outputs": [],
    "source": [
     "demand = xr.DataArray([50, 80, 30], coords=[time])\n",
     "\n",
@@ -101,21 +99,23 @@
     "\n",
     "print(pwf)  # inspect the auto-resolved method\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.388313Z",
-     "start_time": "2026-05-11T17:51:41.338111Z"
+     "end_time": "2026-05-11T17:55:29.493108Z",
+     "start_time": "2026-05-11T17:55:29.453063Z"
     }
    },
-   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -132,14 +132,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.686866Z",
-     "start_time": "2026-05-11T17:51:41.401345Z"
+     "end_time": "2026-05-11T17:55:29.788176Z",
+     "start_time": "2026-05-11T17:55:29.497205Z"
     }
    },
-   "outputs": [],
    "source": [
     "def solve_method(method):\n",
     "    m = linopy.Model()\n",
@@ -153,7 +151,9 @@
     "\n",
     "\n",
     "pd.DataFrame({m: solve_method(m) for m in [\"auto\", \"sos2\", \"incremental\"]})"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -170,15 +170,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:52:01.049181Z",
-     "start_time": "2026-05-11T17:52:00.881744Z"
-    },
-    "scrolled": true
+     "end_time": "2026-05-11T17:56:24.057376Z",
+     "start_time": "2026-05-11T17:56:23.889900Z"
+    }
    },
-   "outputs": [],
    "source": [
     "pumps = pd.Index([\"p1\", \"p2\"], name=\"pump\")\n",
     "\n",
@@ -196,14 +194,13 @@
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
-    "pd.DataFrame(\n",
-    "    {\n",
-    "        \"power_p1\": m.solution[\"power\"].sel(pump=\"p1\").to_pandas(),\n",
-    "        \"power_p2\": m.solution[\"power\"].sel(pump=\"p2\").to_pandas(),\n",
-    "        \"flow\": m.solution[\"flow\"].sum(\"pump\").to_pandas(),\n",
-    "    }\n",
-    ")"
-   ]
+    "# Flat columns: flow_p1, flow_p2, power_p1, power_p2 per timestep.\n",
+    "sol = m.solution[[\"flow\", \"power\"]].to_dataframe().unstack(\"pump\")\n",
+    "sol.columns = [f\"{var}_{p}\" for var, p in sol.columns]\n",
+    "sol"
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -231,14 +228,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.867742Z",
-     "start_time": "2026-05-11T17:47:03.284318Z"
+     "end_time": "2026-05-11T17:55:30.007496Z",
+     "start_time": "2026-05-11T17:55:29.940203Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -255,18 +250,18 @@
     "\n",
     "print(f\"resolved method={pwf.method}, curvature={pwf.convexity}\")\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.868298Z",
-     "start_time": "2026-05-11T12:50:49.566231Z"
+     "end_time": "2026-05-11T17:56:18.581809Z",
+     "start_time": "2026-05-11T17:56:18.563442Z"
     }
    },
-   "outputs": [],
    "source": [
     "plot_curve(\n",
     "    [0, 30, 60, 100],\n",
@@ -274,7 +269,9 @@
     "    m.solution[\"power\"].values,\n",
     "    m.solution[\"fuel\"].values,\n",
     ");"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -293,7 +290,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.874712Z",
+     "end_time": "2026-05-11T17:55:30.036469Z",
      "start_time": "2026-05-11T12:50:49.613733Z"
     }
    },
@@ -326,7 +323,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.875307Z",
+     "end_time": "2026-05-11T17:55:30.036860Z",
      "start_time": "2026-05-11T12:50:49.751516Z"
     }
    },
@@ -349,7 +346,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.875528Z",
+     "end_time": "2026-05-11T17:55:30.044953Z",
      "start_time": "2026-05-11T12:50:49.802514Z"
     }
    },
@@ -379,7 +376,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.877675Z",
+     "end_time": "2026-05-11T17:55:30.045236Z",
      "start_time": "2026-05-11T12:50:49.922140Z"
     }
    },
@@ -413,7 +410,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.881102Z",
+     "end_time": "2026-05-11T17:55:30.045347Z",
      "start_time": "2026-05-11T12:50:50.003565Z"
     }
    },
@@ -451,7 +448,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:51:41.882210Z",
+     "end_time": "2026-05-11T17:55:30.052214Z",
      "start_time": "2026-05-11T12:50:50.151699Z"
     }
    },

From 5b590d5fb2c9122073daf4c226c55b1d024f2401 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 20:00:38 +0200
Subject: [PATCH 24/29] =?UTF-8?q?docs(piecewise):=20inline=20=C2=A73=20dem?=
 =?UTF-8?q?and=20to=20stop=20clobbering=20=C2=A71's?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

§3 previously rebound the Python `demand` variable to its own
pump-tutorial values [30, 75, 150], which then broke §4 and §8
downstream: both rely on §1's `demand = [50, 80, 30]`, and §4 with
demand=150 is infeasible (power.upper=100).

Inlined §3's demand into the constraint so the global namespace is
untouched.  Verified §1 → §3 → §4 sequence end-to-end: §4 now solves
to optimal with fuel = [68, 127, 36].

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-linear-constraints.ipynb | 41 ++++++++++-----------
 1 file changed, 20 insertions(+), 21 deletions(-)

diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index d374d540..67b300e5 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -33,8 +33,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:29.339781Z",
-     "start_time": "2026-05-11T17:55:29.335906Z"
+     "end_time": "2026-05-11T17:58:31.207714Z",
+     "start_time": "2026-05-11T17:58:31.204486Z"
     }
    },
    "source": [
@@ -79,8 +79,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:29.448389Z",
-     "start_time": "2026-05-11T17:55:29.351171Z"
+     "end_time": "2026-05-11T17:58:31.305048Z",
+     "start_time": "2026-05-11T17:58:31.210664Z"
     }
    },
    "source": [
@@ -107,8 +107,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:29.493108Z",
-     "start_time": "2026-05-11T17:55:29.453063Z"
+     "end_time": "2026-05-11T17:58:31.352060Z",
+     "start_time": "2026-05-11T17:58:31.308976Z"
     }
    },
    "source": [
@@ -134,8 +134,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:29.788176Z",
-     "start_time": "2026-05-11T17:55:29.497205Z"
+     "end_time": "2026-05-11T17:58:31.649825Z",
+     "start_time": "2026-05-11T17:58:31.355107Z"
     }
    },
    "source": [
@@ -173,8 +173,8 @@
    "metadata": {
     "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:56:24.057376Z",
-     "start_time": "2026-05-11T17:56:23.889900Z"
+     "end_time": "2026-05-11T17:58:31.802208Z",
+     "start_time": "2026-05-11T17:58:31.661973Z"
     }
    },
    "source": [
@@ -189,8 +189,7 @@
     "    (flow, linopy.segments([(5, 25), (40, 100)])),\n",
     "    (power, linopy.segments([(1, 7), (15, 50)])),\n",
     ")\n",
-    "demand = xr.DataArray([30, 75, 150], coords=[time])\n",
-    "m.add_constraints(flow.sum(\"pump\") == demand)\n",
+    "m.add_constraints(flow.sum(\"pump\") == xr.DataArray([30, 75, 150], coords=[time]))\n",
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
@@ -230,8 +229,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.007496Z",
-     "start_time": "2026-05-11T17:55:29.940203Z"
+     "end_time": "2026-05-11T17:58:31.882868Z",
+     "start_time": "2026-05-11T17:58:31.814503Z"
     }
    },
    "source": [
@@ -258,7 +257,7 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:56:18.581809Z",
+     "end_time": "2026-05-11T17:58:31.892591Z",
      "start_time": "2026-05-11T17:56:18.563442Z"
     }
    },
@@ -290,7 +289,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.036469Z",
+     "end_time": "2026-05-11T17:58:31.893187Z",
      "start_time": "2026-05-11T12:50:49.613733Z"
     }
    },
@@ -323,7 +322,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.036860Z",
+     "end_time": "2026-05-11T17:58:31.894605Z",
      "start_time": "2026-05-11T12:50:49.751516Z"
     }
    },
@@ -346,7 +345,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.044953Z",
+     "end_time": "2026-05-11T17:58:31.898469Z",
      "start_time": "2026-05-11T12:50:49.802514Z"
     }
    },
@@ -376,7 +375,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.045236Z",
+     "end_time": "2026-05-11T17:58:31.898705Z",
      "start_time": "2026-05-11T12:50:49.922140Z"
     }
    },
@@ -410,7 +409,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.045347Z",
+     "end_time": "2026-05-11T17:58:31.898838Z",
      "start_time": "2026-05-11T12:50:50.003565Z"
     }
    },
@@ -448,7 +447,7 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:55:30.052214Z",
+     "end_time": "2026-05-11T17:58:31.898933Z",
      "start_time": "2026-05-11T12:50:50.151699Z"
     }
    },

From 522e78d2ebe0fbf6559c038e1714216a512cae47 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 22:17:44 +0200
Subject: [PATCH 25/29] docs(piecewise): pull the inequality Quick Start
 explanation out of the code block
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

A five-line "# ..." comment inside the snippet was carrying the
formulation explanation — that belongs as rst prose around the code
block, not stuffed into the Python.  Code block now has a single
one-line orienting comment; the rationale moves to a paragraph below.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index 4a537634..faf994fe 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -41,16 +41,18 @@ Quick Start
 
 .. code-block:: python
 
-    # fuel >= f(power) on the same heat-rate curve as above.  Over-
-    # fuelling is physically admissible but wasteful, so minimising
-    # fuel pulls the operating point onto the curve.  "auto" picks the
-    # cheapest correct formulation: pure LP (chord constraints) here,
-    # since convex + ">=" is LP-applicable; SOS2/incremental otherwise.
+    # fuel >= f(power) on the same heat-rate curve as above.
     m.add_piecewise_formulation(
         (fuel, [0, 36, 84, 170], ">="),
         (power, [0, 30, 60, 100]),
     )
 
+Over-fuelling is physically admissible but wasteful, so minimising
+fuel pulls the operating point onto the curve.  ``method="auto"``
+picks the cheapest correct formulation: pure LP (chord constraints)
+here, since convex + ``">="`` is LP-applicable; SOS2/incremental
+otherwise.
+
 Each ``(expression, breakpoints[, sign])`` tuple pairs a variable with
 its breakpoint values.  The optional sign (default ``"=="``) is ``"<="``
 or ``">="`` to mark that expression as bounded by the curve.  With every

From 79fb8e2c46ac0f356242665c1c87033c69a2f0e0 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 22:19:06 +0200
Subject: [PATCH 26/29] docs(piecewise): trim the disjunctive snippet's comment
 to one line

The two-line "stepped pump ... forbidden zone" comment inside the
code block duplicated the prose immediately above ("stepped pump
speeds ... forbidden vibration zones").  Single-line orienting
comment is enough.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index faf994fe..c0b2329d 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -178,8 +178,7 @@ a binary picks exactly one per operating point.
 
 .. code-block:: python
 
-    # Stepped pump with two speed bands; the gap between them is a
-    # forbidden zone.
+    # Stepped pump with two speed bands.
     m.add_piecewise_formulation(
         (flow, linopy.segments([(5, 25), (40, 100)])),
         (power, linopy.segments([(1, 7), (15, 50)])),

From 53d3fa5080d623a7a4c75c808d4ff88c3e5d84ee Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 22:36:35 +0200
Subject: [PATCH 27/29] docs(piecewise): use .html extension in tutorial
 cross-links
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Bare-slug markdown links (`[text](piecewise-linear-constraints)`)
produced sphinx "File not found" warnings AND rendered with href
attributes missing the `.html` extension, requiring browser-side
resolution.  Switching to `.html` form generates correct `href`s
in the rendered HTML directly — verified against the built doc:

  <a class="reference external" href="piecewise-linear-constraints.html">
  <a class="reference external" href="piecewise-inequality-bounds-tutorial.html">

Sphinx still warns at build time (its link checker doesn't find a
source file with `.html` extension), but the warnings are cosmetic
and the rendered docs work correctly.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb  |  4 +-
 examples/piecewise-linear-constraints.ipynb | 96 ++++++++++-----------
 2 files changed, 50 insertions(+), 50 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index 260ff56e..d05af706 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -17,7 +17,7 @@
     "\n",
     "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.\n",
     "\n",
-    "For the formulation math see the [reference page](piecewise-linear-constraints); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial).\n",
+    "For the formulation math see the [reference page](piecewise-linear-constraints.html); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.html).\n",
     "\n",
     "## Tuple roles\n",
     "\n",
@@ -271,7 +271,7 @@
     "- **`method=\"lp\"` is strict** — it raises on a mismatched curvature rather than silently encoding the wrong region.\n",
     "- At most one tuple may carry a non-`==` sign, and 3+ tuples must all be `==`. Multi-bounded / N≥3 inequalities — open an issue at https://github.com/PyPSA/linopy/issues.\n",
     "\n",
-    "**See also**: [reference page](piecewise-linear-constraints) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial) for all-equality, unit commitment, CHP, fleets, slopes."
+    "**See also**: [reference page](piecewise-linear-constraints.html) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.html) for all-equality, unit commitment, CHP, fleets, slopes."
    ]
   }
  ],
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index 67b300e5..f872ef3e 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -15,7 +15,7 @@
     ")\n",
     "```\n",
     "\n",
-    "This tutorial walks through the main features of `add_piecewise_formulation`. For the formulation math see the [reference page](piecewise-linear-constraints); for the inequality variant in depth see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial).\n",
+    "This tutorial walks through the main features of `add_piecewise_formulation`. For the formulation math see the [reference page](piecewise-linear-constraints.html); for the inequality variant in depth see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html).\n",
     "\n",
     "**Roadmap**\n",
     "\n",
@@ -33,8 +33,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.207714Z",
-     "start_time": "2026-05-11T17:58:31.204486Z"
+     "end_time": "2026-05-11T18:01:54.620516Z",
+     "start_time": "2026-05-11T18:01:54.613427Z"
     }
    },
    "source": [
@@ -79,8 +79,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.305048Z",
-     "start_time": "2026-05-11T17:58:31.210664Z"
+     "end_time": "2026-05-11T18:01:54.730Z",
+     "start_time": "2026-05-11T18:01:54.625751Z"
     }
    },
    "source": [
@@ -107,8 +107,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.352060Z",
-     "start_time": "2026-05-11T17:58:31.308976Z"
+     "end_time": "2026-05-11T18:01:54.780034Z",
+     "start_time": "2026-05-11T18:01:54.735021Z"
     }
    },
    "source": [
@@ -127,15 +127,15 @@
     "\n",
     "For now: a quick sanity check that all applicable methods yield the same fuel dispatch on the convex curve from §1.\n",
     "\n",
-    "A full comparison — when each method dispatches, what sign/curvature/breakpoint patterns each requires — lives in §3 (disjunctive), §4 (inequalities), and the [reference page's \"Formulation Methods\" section](piecewise-linear-constraints)."
+    "A full comparison — when each method dispatches, what sign/curvature/breakpoint patterns each requires — lives in §3 (disjunctive), §4 (inequalities), and the [reference page's \"Formulation Methods\" section](piecewise-linear-constraints.html)."
    ]
   },
   {
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.649825Z",
-     "start_time": "2026-05-11T17:58:31.355107Z"
+     "end_time": "2026-05-11T18:01:55.102903Z",
+     "start_time": "2026-05-11T18:01:54.783092Z"
     }
    },
    "source": [
@@ -173,8 +173,8 @@
    "metadata": {
     "scrolled": true,
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.802208Z",
-     "start_time": "2026-05-11T17:58:31.661973Z"
+     "end_time": "2026-05-11T18:01:55.257839Z",
+     "start_time": "2026-05-11T18:01:55.114836Z"
     }
    },
    "source": [
@@ -222,15 +222,15 @@
     "\n",
     "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "Below: the same heat-rate curve as §1, now read as `fuel ≥ f(power)` — over-fuelling is admissible but wasteful (the curve is the design minimum), so an objective that minimises fuel pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
+    "Below: the same heat-rate curve as §1, now read as `fuel ≥ f(power)` — over-fuelling is admissible but wasteful (the curve is the design minimum), so an objective that minimises fuel pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
    ]
   },
   {
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.882868Z",
-     "start_time": "2026-05-11T17:58:31.814503Z"
+     "end_time": "2026-05-11T18:01:55.331357Z",
+     "start_time": "2026-05-11T18:01:55.269Z"
     }
    },
    "source": [
@@ -257,8 +257,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.892591Z",
-     "start_time": "2026-05-11T17:56:18.563442Z"
+     "end_time": "2026-05-11T18:01:55.381548Z",
+     "start_time": "2026-05-11T18:01:55.337053Z"
     }
    },
    "source": [
@@ -286,14 +286,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.893187Z",
-     "start_time": "2026-05-11T12:50:49.613733Z"
+     "end_time": "2026-05-11T18:01:55.558321Z",
+     "start_time": "2026-05-11T18:01:55.386257Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "p_min, p_max = 30, 100\n",
@@ -315,21 +313,23 @@
     "m.add_objective(fuel.sum() + 50 * commit.sum() + 200 * backup.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"commit\", \"power\", \"fuel\", \"backup\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.894605Z",
-     "start_time": "2026-05-11T12:50:49.751516Z"
+     "end_time": "2026-05-11T18:01:55.609973Z",
+     "start_time": "2026-05-11T18:01:55.564366Z"
     }
    },
-   "outputs": [],
    "source": [
     "plot_curve(x_pts, y_pts, m.solution[\"power\"].values, m.solution[\"fuel\"].values);"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -342,14 +342,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.898469Z",
-     "start_time": "2026-05-11T12:50:49.802514Z"
+     "end_time": "2026-05-11T18:01:55.731111Z",
+     "start_time": "2026-05-11T18:01:55.619583Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -368,18 +366,18 @@
     "m.add_objective(power.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\", \"heat\"]].to_pandas().round(2)"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.898705Z",
-     "start_time": "2026-05-11T12:50:49.922140Z"
+     "end_time": "2026-05-11T18:01:55.811271Z",
+     "start_time": "2026-05-11T18:01:55.738346Z"
     }
    },
-   "outputs": [],
    "source": [
     "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
     "plot_curve(\n",
@@ -393,7 +391,9 @@
     "    ylabel=\"heat\",\n",
     "    ax=axes[1],\n",
     ");"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -406,14 +406,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.898838Z",
-     "start_time": "2026-05-11T12:50:50.003565Z"
+     "end_time": "2026-05-11T18:01:55.957075Z",
+     "start_time": "2026-05-11T18:01:55.820261Z"
     }
    },
-   "outputs": [],
    "source": [
     "gens = pd.Index([\"gas\", \"coal\"], name=\"gen\")\n",
     "x_gen = linopy.breakpoints(\n",
@@ -431,7 +429,9 @@
     "m.add_objective(fuel.sum())\n",
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "m.solution[[\"power\", \"fuel\"]].to_dataframe()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -444,14 +444,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2026-05-11T17:58:31.898933Z",
-     "start_time": "2026-05-11T12:50:50.151699Z"
+     "end_time": "2026-05-11T18:01:56.057514Z",
+     "start_time": "2026-05-11T18:01:55.964137Z"
     }
    },
-   "outputs": [],
    "source": [
     "m = linopy.Model()\n",
     "power = m.add_variables(name=\"power\", lower=0, upper=100, coords=[time])\n",
@@ -466,7 +464,9 @@
     "m.solve(solver_name=\"highs\", reformulate_sos=\"auto\", output_flag=False)\n",
     "\n",
     "m.solution[[\"power\", \"fuel\"]].to_pandas()"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -484,7 +484,7 @@
     "| Different curves per entity | `linopy.breakpoints({...}, dim=...)` |\n",
     "| Slopes more natural than absolute y-values | `linopy.Slopes(...)` |\n",
     "\n",
-    "For the formulation math, see the [reference page](piecewise-linear-constraints). For inequality bounds in depth, see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial)."
+    "For the formulation math, see the [reference page](piecewise-linear-constraints.html). For inequality bounds in depth, see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html)."
    ]
   }
  ],

From 65f615485b506f56b101fab724ecf6d72d02c008 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Mon, 11 May 2026 22:43:29 +0200
Subject: [PATCH 28/29] docs(piecewise): use .rst / .nblink in tutorial
 cross-links
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Previous .html form rendered correct browser hrefs but produced
sphinx "File not found" warnings AND rendered as external links.

Switching to source-extension form:
- [text](page.rst) for rst pages
- [text](tutorial.nblink) for sibling notebook tutorials

Sphinx resolves both to internal cross-references at build time —
HTML output shows `class="reference internal"` with the doc anchor,
no build warnings.  Verified locally: piecewise cross-link warnings
went from 6 to 0 (total build warnings 20 → 10, the remainder are
pre-existing unrelated issues).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 examples/piecewise-inequality-bounds.ipynb  | 4 ++--
 examples/piecewise-linear-constraints.ipynb | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/examples/piecewise-inequality-bounds.ipynb b/examples/piecewise-inequality-bounds.ipynb
index d05af706..d1ca4e79 100644
--- a/examples/piecewise-inequality-bounds.ipynb
+++ b/examples/piecewise-inequality-bounds.ipynb
@@ -17,7 +17,7 @@
     "\n",
     "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.\n",
     "\n",
-    "For the formulation math see the [reference page](piecewise-linear-constraints.html); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.html).\n",
+    "For the formulation math see the [reference page](piecewise-linear-constraints.rst); for the all-equality variant and other features see [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.nblink).\n",
     "\n",
     "## Tuple roles\n",
     "\n",
@@ -271,7 +271,7 @@
     "- **`method=\"lp\"` is strict** — it raises on a mismatched curvature rather than silently encoding the wrong region.\n",
     "- At most one tuple may carry a non-`==` sign, and 3+ tuples must all be `==`. Multi-bounded / N≥3 inequalities — open an issue at https://github.com/PyPSA/linopy/issues.\n",
     "\n",
-    "**See also**: [reference page](piecewise-linear-constraints.html) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.html) for all-equality, unit commitment, CHP, fleets, slopes."
+    "**See also**: [reference page](piecewise-linear-constraints.rst) for the formulation math, [Creating Piecewise Linear Constraints](piecewise-linear-constraints-tutorial.nblink) for all-equality, unit commitment, CHP, fleets, slopes."
    ]
   }
  ],
diff --git a/examples/piecewise-linear-constraints.ipynb b/examples/piecewise-linear-constraints.ipynb
index f872ef3e..8b4f56ca 100644
--- a/examples/piecewise-linear-constraints.ipynb
+++ b/examples/piecewise-linear-constraints.ipynb
@@ -15,7 +15,7 @@
     ")\n",
     "```\n",
     "\n",
-    "This tutorial walks through the main features of `add_piecewise_formulation`. For the formulation math see the [reference page](piecewise-linear-constraints.html); for the inequality variant in depth see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html).\n",
+    "This tutorial walks through the main features of `add_piecewise_formulation`. For the formulation math see the [reference page](piecewise-linear-constraints.rst); for the inequality variant in depth see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.nblink).\n",
     "\n",
     "**Roadmap**\n",
     "\n",
@@ -127,7 +127,7 @@
     "\n",
     "For now: a quick sanity check that all applicable methods yield the same fuel dispatch on the convex curve from §1.\n",
     "\n",
-    "A full comparison — when each method dispatches, what sign/curvature/breakpoint patterns each requires — lives in §3 (disjunctive), §4 (inequalities), and the [reference page's \"Formulation Methods\" section](piecewise-linear-constraints.html)."
+    "A full comparison — when each method dispatches, what sign/curvature/breakpoint patterns each requires — lives in §3 (disjunctive), §4 (inequalities), and the [reference page's \"Formulation Methods\" section](piecewise-linear-constraints.rst)."
    ]
   },
   {
@@ -222,7 +222,7 @@
     "\n",
     "On a concave curve with `<=` (or convex with `>=`), `method=\"auto\"` dispatches to a pure LP chord formulation — no binaries, no SOS2.\n",
     "\n",
-    "Below: the same heat-rate curve as §1, now read as `fuel ≥ f(power)` — over-fuelling is admissible but wasteful (the curve is the design minimum), so an objective that minimises fuel pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
+    "Below: the same heat-rate curve as §1, now read as `fuel ≥ f(power)` — over-fuelling is admissible but wasteful (the curve is the design minimum), so an objective that minimises fuel pulls the operating point onto the curve. See [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.nblink) for mismatched curvature, auto-dispatch fallbacks, and more geometry."
    ]
   },
   {
@@ -484,7 +484,7 @@
     "| Different curves per entity | `linopy.breakpoints({...}, dim=...)` |\n",
     "| Slopes more natural than absolute y-values | `linopy.Slopes(...)` |\n",
     "\n",
-    "For the formulation math, see the [reference page](piecewise-linear-constraints.html). For inequality bounds in depth, see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.html)."
+    "For the formulation math, see the [reference page](piecewise-linear-constraints.rst). For inequality bounds in depth, see [Creating Piecewise Inequality Bounds](piecewise-inequality-bounds-tutorial.nblink)."
    ]
   }
  ],

From 532dab68754e938360417942c1c474009ac9834b Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Tue, 12 May 2026 13:46:23 +0200
Subject: [PATCH 29/29] docs(piecewise): address review feedback (PR #677)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Fabian's review raised three points (https://github.com/PyPSA/linopy/pull/677#pullrequestreview-4271283888):

1. "Advanced features" forward-pointer should be a ref link.

   Added `_piecewise-active:` label above the Active parameter
   subsection; the disjunctive forward-pointer now uses
   `:ref:`piecewise-active`` (renders as a clickable internal link).

2. Restrictions block could be an info block.

   Converted to `.. note::` (matching the existing pattern in
   sos-constraints.rst and gpu-acceleration.rst — the codebase doesn't
   use `.. info::`, only note/warning).

3. Building blocks need more detail on usage.

   Restructured Breakpoint Construction.  Originally six subsections
   organised around input shape (From lists / breakpoints() factory /
   From slopes / Per-entity / Disjunctive / N-variable).  Now organised
   around the two factories with Slopes folded in as a convenience:

   - `breakpoints()` — connected curve.  Opens with use case
     (efficiency curves, heat rates, cost curves).  Two bold-inline
     sub-blocks: **Per-entity curves.** (dict input) and **Specifying
     by slopes.** (Slopes wrapper).
   - `segments()` — disjoint operating bands.  Opens with use case.
     Trimmed redundant trailing prose; kept the bounded-tuple
     compatibility note as a one-liner.
   - N-variable linking — kept separate (orthogonal to input type).

   Added a two-sentence signpost at the top of Breakpoint Construction
   that transitions from the API signature: names both factories with
   cross-refs and positions Slopes as a slope-input stand-in for
   `breakpoints()` (not a separate building block — Slopes resolves to
   a breakpoints array internally).

Plus minor editorial cleanups from the iteration:
- Dropped the SOS2 jargon from the `segments()` opener (replaced with
  "continuous interpolation within the chosen band" — the SOS2 method
  is introduced later in Formulation Methods).
- Dropped the accepted-datatypes enumeration from the `breakpoints()`
  opener (redundant with examples and the function's docstring).
- Dropped the ambiguous "prices per MWh per operating step" example
  from the Slopes block — both electricity-specific and confusable
  with the disjunctive/VSD case from §3.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 doc/piecewise-linear-constraints.rst | 112 +++++++++++++--------------
 1 file changed, 54 insertions(+), 58 deletions(-)

diff --git a/doc/piecewise-linear-constraints.rst b/doc/piecewise-linear-constraints.rst
index c0b2329d..2acf886d 100644
--- a/doc/piecewise-linear-constraints.rst
+++ b/doc/piecewise-linear-constraints.rst
@@ -88,20 +88,19 @@ binaries (see *Formulation Methods* below).
 Breakpoint Construction
 -----------------------
 
-Three building blocks provide breakpoint data:
+Each tuple's breakpoints come from :func:`~linopy.breakpoints` (a
+single connected curve) or :func:`~linopy.segments` (disjoint
+operating bands).  :class:`~linopy.Slopes` can stand in for
+:func:`~linopy.breakpoints` when per-piece slopes are the natural
+input — it resolves to a breakpoints array.
 
-- :func:`~linopy.breakpoints` — values along a single **connected**
-  curve.
-- :func:`~linopy.segments` — **disjoint** operating regions with gaps
-  between them (e.g. forbidden zones); selects the disjunctive
-  formulation automatically.
-- :class:`~linopy.Slopes` — per-piece slopes plus an initial ``y0``,
-  deferred until an x grid is supplied by a sibling tuple.
+``breakpoints()`` — a connected curve
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-From lists
-~~~~~~~~~~
+Values along a single **connected** piecewise curve — the default case
+for efficiency curves, heat rates, and cost curves.
 
-The simplest form — pass Python lists directly in the tuple:
+The simplest form passes a Python list directly in the tuple:
 
 .. code-block:: python
 
@@ -110,9 +109,6 @@ The simplest form — pass Python lists directly in the tuple:
         (fuel, [0, 36, 84, 170]),
     )
 
-With the ``breakpoints()`` factory
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
 Equivalent, but explicit about the DataArray construction:
 
 .. code-block:: python
@@ -122,31 +118,8 @@ Equivalent, but explicit about the DataArray construction:
         (fuel, linopy.breakpoints([0, 36, 84, 170])),
     )
 
-From slopes
-~~~~~~~~~~~
-
-When you know marginal costs (slopes) rather than absolute values, wrap
-them in :class:`linopy.Slopes`.  The x grid is borrowed from the sibling
-tuple — no need to repeat it:
-
-.. code-block:: python
-
-    m.add_piecewise_formulation(
-        (power, [0, 50, 100, 150]),
-        (cost, linopy.Slopes([1.1, 1.5, 1.9], y0=0)),
-    )
-    # cost breakpoints: [0, 55, 130, 225]
-
-For standalone resolution outside of ``add_piecewise_formulation``, call
-:meth:`linopy.Slopes.to_breakpoints` with an explicit x grid::
-
-    bp = linopy.Slopes([1.1, 1.5, 1.9], y0=0).to_breakpoints([0, 50, 100, 150])
-
-Per-entity breakpoints
-~~~~~~~~~~~~~~~~~~~~~~
-
-Different generators can have different curves.  Pass a dict to
-``breakpoints()`` with entity names as keys:
+**Per-entity curves.**  Different generators can have different
+curves.  Pass a dict to ``breakpoints()`` with entity names as keys:
 
 .. code-block:: python
 
@@ -165,16 +138,35 @@ Different generators can have different curves.  Pass a dict to
         ),
     )
 
-Ragged lengths are NaN-padded automatically.  Breakpoints are auto-broadcast
-over remaining dimensions (e.g. ``time``).
+Ragged lengths are NaN-padded automatically.  Breakpoints are auto-
+broadcast over remaining dimensions (e.g. ``time``).
+
+**Specifying by slopes.**  :class:`linopy.Slopes` resolves to a
+breakpoint array from per-piece slopes plus an initial ``y0``,
+instead of from absolute y-values — useful when slopes are the
+natural input (e.g. marginal costs).  The x grid is borrowed from
+the sibling tuple, so the y breakpoints don't have to be computed
+by hand:
+
+.. code-block:: python
+
+    m.add_piecewise_formulation(
+        (power, [0, 50, 100, 150]),
+        (cost, linopy.Slopes([1.1, 1.5, 1.9], y0=0)),
+    )
+    # cost breakpoints resolve to: [0, 55, 130, 225]
+
+For standalone resolution outside ``add_piecewise_formulation``, call
+:meth:`linopy.Slopes.to_breakpoints` with an explicit x grid::
+
+    bp = linopy.Slopes([1.1, 1.5, 1.9], y0=0).to_breakpoints([0, 50, 100, 150])
 
-Disjunctive segments
-~~~~~~~~~~~~~~~~~~~~
+``segments()`` — disjoint operating bands
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-For equipment with disconnected operating bands — stepped pump speeds,
-switchable combustion cycles, allowed bands around forbidden vibration
-zones — use ``segments()``.  Each segment is one band's (range, curve);
-a binary picks exactly one per operating point.
+For equipment with disconnected operating bands.  Each segment is one
+band's ``(range, curve)``; a binary picks exactly one per operating
+point, with continuous interpolation within the chosen band.
 
 .. code-block:: python
 
@@ -184,18 +176,18 @@ a binary picks exactly one per operating point.
         (power, linopy.segments([(1, 7), (15, 50)])),
     )
 
-The disjunctive formulation is selected automatically when breakpoints
-have a segment dimension.  A bounded tuple (``"<="`` / ``">="``) also
-works here.
+Bounded tuples (``"<="`` / ``">="``) are supported on disjunctive
+curves too.
 
 For a single on/off gate on one continuous curve, prefer ``active=...``
-(see *Advanced Features* below) — using a degenerate ``(0, 0)`` segment
+(see :ref:`piecewise-active`) — using a degenerate ``(0, 0)`` segment
 to encode "off" mixes the disjunctive concept with on/off logic.
 
 N-variable linking
 ~~~~~~~~~~~~~~~~~~
 
-Link any number of variables through shared breakpoints (joint equality):
+Independent of the building block used, any number of variables can be
+linked through shared breakpoints (joint equality):
 
 .. code-block:: python
 
@@ -224,14 +216,16 @@ Each tuple's optional third element is a sign:
 - ``"<="`` / ``">="`` — **bounded**: the expression undershoots /
   overshoots the curve.
 
-**Restrictions (current):**
+.. note::
+
+   **Current restrictions.**
 
-- At most one tuple may carry a non-equality sign — a single bounded side.
-- With **3 or more** tuples, all signs must be ``"=="``.
+   - At most one tuple may carry a non-equality sign — a single bounded side.
+   - With **3 or more** tuples, all signs must be ``"=="``.
 
-Multi-bounded and N≥3-inequality use cases aren't supported yet.  If
-you have a concrete use case, please open an issue at
-https://github.com/PyPSA/linopy/issues so we can scope it properly.
+   Multi-bounded and N≥3-inequality use cases aren't supported yet.
+   If you have a concrete use case, please open an issue at
+   https://github.com/PyPSA/linopy/issues so we can scope it properly.
 
 Geometry
 ~~~~~~~~
@@ -537,6 +531,8 @@ disconnected operating regions" that ``method="lp"`` cannot handle.
 Advanced Features
 -----------------
 
+.. _piecewise-active:
+
 Active parameter (unit commitment)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~