PyPSA
diff --git a/‎dev-scripts/piecewise-feasibility-tests-walkthrough.ipynb‎
Lines changed: 348 additions & 0 deletions b/‎dev-scripts/piecewise-feasibility-tests-walkthrough.ipynb‎
Lines changed: 348 additions & 0 deletions
diff --git a/‎doc/api.rst‎
Lines changed: 4 additions & 2 deletions b/‎doc/api.rst‎
Lines changed: 4 additions & 2 deletions
diff --git a/‎doc/index.rst‎
Lines changed: 0 additions & 1 deletion b/‎doc/index.rst‎
Lines changed: 0 additions & 1 deletion
diff --git a/‎doc/piecewise-inequality-bounds-tutorial.nblink‎
Lines changed: 3 additions & 0 deletions b/‎doc/piecewise-inequality-bounds-tutorial.nblink‎
Lines changed: 3 additions & 0 deletions
@@ -0,0 +1,348 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# `test_piecewise_feasibility.py` — visual walkthrough\n",
+    "\n",
+    "**Purpose:** document what each test class in `test/test_piecewise_feasibility.py` actually probes, with pictures. Intended as review aid for the PR — **not** merged into master.\n",
+    "\n",
+    "The test file stress-tests the claim that `add_piecewise_formulation(sign=\"<=\"/\">=\")` yields the **same feasible region** for `(x, y)` regardless of which method (`lp` / `sos2` / `incremental`) dispatches the formulation, on curves where all three are applicable.\n",
+    "\n",
+    "Four test classes:\n",
+    "\n",
+    "| class | what it probes | scope |\n",
+    "|---|---|---|\n",
+    "| `TestRotatedObjective` | support-function equivalence — 16 rotation directions | the strong test |\n",
+    "| `TestDomainBoundary` | `x` outside `[x_min, x_max]` is infeasible | LP explicit vs SOS2 implicit |\n",
+    "| `TestPointwiseInfeasibility` | `y` just past `f(x)` is infeasible | targeted sanity check |\n",
+    "| `TestNVariableInequality` | 3-variable: first tuple bounded, rest equality | SOS2 vs incremental only |\n",
+    "\n",
+    "Below: one visualization per class.\n",
+    "\n",
+    "*Run this notebook from the repository root so that `from test.test_piecewise_feasibility import ...` resolves.*"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:47.376525Z",
+     "start_time": "2026-04-23T08:00:46.142492Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from test.test_piecewise_feasibility import (\n",
+    "    CURVES,\n",
+    "    Y_HI,\n",
+    "    Y_LO,\n",
+    "    Curve,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Shared primitive: draw the curve and its feasible region"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:47.384959Z",
+     "start_time": "2026-04-23T08:00:47.381361Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def draw_curve_and_region(ax, curve: Curve, *, shade: bool = True) -> None:\n",
+    "    \"\"\"Plot breakpoints + shade the feasible region (hypograph or epigraph).\"\"\"\n",
+    "    xs = np.array(curve.x_pts)\n",
+    "    ys = np.array(curve.y_pts)\n",
+    "    ax.plot(xs, ys, \"o-\", color=\"C0\", lw=2, label=\"breakpoints\")\n",
+    "\n",
+    "    if shade:\n",
+    "        if curve.sign == \"<=\":\n",
+    "            ax.fill_between(\n",
+    "                xs,\n",
+    "                np.full_like(ys, Y_LO),\n",
+    "                ys,\n",
+    "                alpha=0.15,\n",
+    "                color=\"C0\",\n",
+    "                label=f\"feasible: y {curve.sign} f(x)\",\n",
+    "            )\n",
+    "        else:\n",
+    "            ax.fill_between(\n",
+    "                xs,\n",
+    "                ys,\n",
+    "                np.full_like(ys, Y_HI),\n",
+    "                alpha=0.15,\n",
+    "                color=\"C0\",\n",
+    "                label=f\"feasible: y {curve.sign} f(x)\",\n",
+    "            )\n",
+    "\n",
+    "    pad_x = 0.15 * (xs.max() - xs.min())\n",
+    "    pad_y = 0.15 * (ys.max() - ys.min()) + 1\n",
+    "    ax.set_xlim(xs.min() - pad_x, xs.max() + pad_x)\n",
+    "    ax.set_ylim(ys.min() - pad_y, ys.max() + pad_y)\n",
+    "    ax.set_xlabel(\"x\")\n",
+    "    ax.set_ylabel(\"y\")\n",
+    "    ax.grid(alpha=0.3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## `TestRotatedObjective` — the strong test\n",
+    "\n",
+    "For every direction `(α, β)` on the unit circle, minimize `α·x + β·y` under the PWL. The answer is the **support function** of the feasible region in direction `(α, β)` — and for a convex region, the support function uniquely determines the region. If LP and SOS2/incremental give the same support-function value for 16 directions, their feasible regions are identical.\n",
+    "\n",
+    "Each red dot below is the extreme point the solver lands at for one direction. The arrows show the objective-push direction. A failure would manifest as one method's dot landing at a different vertex than the oracle's."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:47.829483Z",
+     "start_time": "2026-04-23T08:00:47.388542Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def panel_rotated_objective(ax, curve: Curve, n_dirs: int = 16) -> None:\n",
+    "    draw_curve_and_region(ax, curve)\n",
+    "    xs, ys = np.array(curve.x_pts), np.array(curve.y_pts)\n",
+    "    cx = 0.5 * (xs.min() + xs.max())\n",
+    "    cy = 0.5 * (ys.min() + ys.max())\n",
+    "    arrow_len = 0.25 * min(xs.max() - xs.min(), (ys.max() - ys.min()) + 5)\n",
+    "\n",
+    "    for i in range(n_dirs):\n",
+    "        theta = 2 * np.pi * i / n_dirs\n",
+    "        alpha, beta = np.cos(theta), np.sin(theta)\n",
+    "        ax.annotate(\n",
+    "            \"\",\n",
+    "            xytext=(cx, cy),\n",
+    "            xy=(cx + arrow_len * alpha, cy + arrow_len * beta),\n",
+    "            arrowprops=dict(arrowstyle=\"->\", color=\"C3\", alpha=0.4, lw=1),\n",
+    "        )\n",
+    "        # Oracle extreme point in this direction\n",
+    "        verts = curve.vertices()\n",
+    "        extreme = min(verts, key=lambda v: alpha * v[0] + beta * v[1])\n",
+    "        ax.plot(*extreme, \"o\", color=\"C3\", ms=4, alpha=0.7)\n",
+    "\n",
+    "    ax.plot([], [], \"o\", color=\"C3\", alpha=0.7, label=f\"{n_dirs} extreme points\")\n",
+    "    ax.legend(loc=\"upper left\", fontsize=8)\n",
+    "    ax.set_title(f\"{curve.name}  (sign={curve.sign})\")\n",
+    "\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
+    "panel_rotated_objective(axes[0], CURVES[0])  # concave-smooth\n",
+    "panel_rotated_objective(axes[1], CURVES[2])  # convex-steep\n",
+    "panel_rotated_objective(axes[2], CURVES[5])  # two-segment\n",
+    "fig.suptitle(\n",
+    "    \"TestRotatedObjective — support function sampled at 16 directions\", fontsize=12\n",
+    ")\n",
+    "plt.tight_layout();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice the **dots cluster at the curve breakpoints** (top edges) and at the **bottom corners** `(x_min, Y_LO)`, `(x_max, Y_LO)`. That's because the feasible region is a polygon: linear objectives always attain their optimum at a vertex.\n",
+    "\n",
+    "The 288 pytest items (6 curves × 3 methods × 16 directions) check that all three methods land at the same extreme point for every direction."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## `TestDomainBoundary` — enforce `x ∈ [x_min, x_max]`\n",
+    "\n",
+    "LP enforces this with an explicit constraint; SOS2/incremental enforce it implicitly via `sum(λ) = 1`. Two different implementations of the same bound — worth a direct probe."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:48.103641Z",
+     "start_time": "2026-04-23T08:00:47.835275Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def panel_domain_boundary(ax, curve: Curve) -> None:\n",
+    "    draw_curve_and_region(ax, curve)\n",
+    "    xs = np.array(curve.x_pts)\n",
+    "    y_span = ax.get_ylim()\n",
+    "    ax.axvline(xs[0], color=\"C2\", lw=1.5, label=f\"x_min={xs[0]}\")\n",
+    "    ax.axvline(xs[-1], color=\"C2\", lw=1.5, label=f\"x_max={xs[-1]}\")\n",
+    "    ax.axvline(xs[0] - 1, color=\"C3\", lw=1.5, ls=\"--\")\n",
+    "    ax.axvline(xs[-1] + 1, color=\"C3\", lw=1.5, ls=\"--\")\n",
+    "    yy = y_span[1] - 0.12 * (y_span[1] - y_span[0])\n",
+    "    ax.text(\n",
+    "        xs[0] - 1, yy, \"INFEASIBLE\\n(x < x_min)\", ha=\"center\", fontsize=8, color=\"C3\"\n",
+    "    )\n",
+    "    ax.text(\n",
+    "        xs[-1] + 1, yy, \"INFEASIBLE\\n(x > x_max)\", ha=\"center\", fontsize=8, color=\"C3\"\n",
+    "    )\n",
+    "    ax.legend(loc=\"lower center\", fontsize=7)\n",
+    "    ax.set_title(f\"{curve.name} — domain probe\")\n",
+    "\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
+    "panel_domain_boundary(axes[0], CURVES[0])  # concave-smooth\n",
+    "panel_domain_boundary(axes[1], CURVES[1])  # concave-shifted (negative domain)\n",
+    "panel_domain_boundary(axes[2], CURVES[5])  # two-segment\n",
+    "fig.suptitle(\n",
+    "    \"TestDomainBoundary — x outside the breakpoint range is infeasible\", fontsize=12\n",
+    ")\n",
+    "plt.tight_layout();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## `TestPointwiseInfeasibility` — y just past the curve\n",
+    "\n",
+    "Rotated objectives probe *extremes*; this test specifically nudges `y` past `f(x)` by a small margin (`0.01`) and asserts infeasibility. Catches NaN-mask or off-by-one-segment bugs that might accidentally allow slack."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:48.366674Z",
+     "start_time": "2026-04-23T08:00:48.112127Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def panel_pointwise(ax, curve: Curve) -> None:\n",
+    "    draw_curve_and_region(ax, curve)\n",
+    "    xs = np.array(curve.x_pts)\n",
+    "    x_mid = 0.5 * (xs[0] + xs[-1])\n",
+    "    fx = curve.f(x_mid)\n",
+    "    y_bad = fx + 0.01 if curve.sign == \"<=\" else fx - 0.01\n",
+    "    ax.plot(x_mid, fx, \"o\", color=\"C2\", ms=9, label=f\"on curve: f({x_mid:g})={fx:g}\")\n",
+    "    ax.plot(\n",
+    "        x_mid, y_bad, \"x\", color=\"C3\", ms=14, mew=3, label=f\"infeasible: y={y_bad:g}\"\n",
+    "    )\n",
+    "    ax.legend(loc=\"lower right\", fontsize=7)\n",
+    "    ax.set_title(f\"{curve.name} — nudge past f(x)\")\n",
+    "\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
+    "panel_pointwise(axes[0], CURVES[0])  # concave-smooth, sign=\"<=\"\n",
+    "panel_pointwise(axes[1], CURVES[2])  # convex-steep, sign=\">=\"\n",
+    "panel_pointwise(axes[2], CURVES[4])  # linear-gte\n",
+    "fig.suptitle(\n",
+    "    \"TestPointwiseInfeasibility — y past the curve by 0.01 in the sign direction\",\n",
+    "    fontsize=12,\n",
+    ")\n",
+    "plt.tight_layout();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## `TestNVariableInequality` — 3-variable sign split\n",
+    "\n",
+    "With three tuples `(fuel, power, heat)` and `sign=\"<=\"`:\n",
+    "- `fuel` (the **first** tuple) is **bounded above** by its interpolated value,\n",
+    "- `power` and `heat` (remaining tuples) are **forced to equality** — pinned on the curve.\n",
+    "\n",
+    "LP doesn't support N > 2 tuples, so this class compares SOS2 vs incremental only. The 3D plot shows the CHP curve and the 7 test points (one per `power_fix`) that both methods must agree on."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2026-04-23T08:00:48.489668Z",
+     "start_time": "2026-04-23T08:00:48.371526Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "bp = {\n",
+    "    \"power\": np.array([0, 30, 60, 100]),\n",
+    "    \"fuel\": np.array([0, 40, 85, 160]),\n",
+    "    \"heat\": np.array([0, 25, 55, 95]),\n",
+    "}\n",
+    "\n",
+    "fig = plt.figure(figsize=(9, 6.5))\n",
+    "ax = fig.add_subplot(projection=\"3d\")\n",
+    "ax.plot(\n",
+    "    bp[\"power\"], bp[\"fuel\"], bp[\"heat\"], \"o-\", color=\"C0\", lw=2, label=\"CHP breakpoints\"\n",
+    ")\n",
+    "\n",
+    "for p in [0, 15, 30, 45, 60, 80, 100]:\n",
+    "    f = np.interp(p, bp[\"power\"], bp[\"fuel\"])\n",
+    "    h = np.interp(p, bp[\"power\"], bp[\"heat\"])\n",
+    "    ax.plot([p], [f], [h], \"o\", color=\"C3\", ms=7)\n",
+    "    # drop to base plane\n",
+    "    ax.plot([p, p], [f, 0], [h, h], color=\"C3\", alpha=0.3, lw=0.8)\n",
+    "\n",
+    "ax.set_xlabel(\"power\")\n",
+    "ax.set_ylabel(\"fuel\")\n",
+    "ax.set_zlabel(\"heat\")\n",
+    "ax.plot(\n",
+    "    [],\n",
+    "    [],\n",
+    "    \"o\",\n",
+    "    color=\"C3\",\n",
+    "    label=\"7 test points — power pinned,\\nfuel at upper bound, heat on curve\",\n",
+    ")\n",
+    "ax.set_title('TestNVariableInequality — CHP curve (sign=\"<=\")')\n",
+    "ax.legend(loc=\"upper left\", fontsize=8);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## What a failing test would tell you\n",
+    "\n",
+    "- **Rotated objective fails**: the methods disagree on the feasible region in some direction. The failure message includes the attained `(x, y)` point — you'd see which extreme point one method landed at that the others didn't.\n",
+    "- **Domain boundary fails**: one method lets `x` escape `[x_min, x_max]`. LP path most likely: the domain-bound constraint was dropped. SOS2 path: the `sum(λ) = 1` constraint was weakened.\n",
+    "- **Pointwise infeasibility fails**: one method accepts a point past the curve. Most often a NaN-mask bug in per-entity formulations, or a wrong segment getting picked.\n",
+    "- **N-variable fails**: the sign split went wrong — either an input leaked into the signed link or the first-tuple convention is misrouting.\n",
+    "\n",
+    "All 356 pytest items are currently green at `TOL = 1e-5`."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.13.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
@@ -18,10 +18,12 @@ Creating a model
     model.Model.add_variables
     model.Model.add_constraints
     model.Model.add_objective
-    model.Model.add_piecewise_constraints
-    piecewise.piecewise
+    model.Model.add_piecewise_formulation
+    piecewise.PiecewiseFormulation
     piecewise.breakpoints
     piecewise.segments
+    piecewise.slopes_to_points
+    piecewise.tangent_lines
     model.Model.linexpr
     model.Model.remove_constraints
     model.Model.copy
 
@@ -114,7 +114,6 @@ This package is published under MIT license.
    coordinate-alignment
    sos-constraints
    piecewise-linear-constraints
-   piecewise-linear-constraints-tutorial
    manipulating-models
    testing-framework
    transport-tutorial
 
@@ -0,0 +1,3 @@
+{
+    "path": "../examples/piecewise-inequality-bounds.ipynb"
+}
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,3 @@`
	`1`	`+{`
	`2`	`+ "path": "../examples/piecewise-inequality-bounds.ipynb"`
	`3`	`+}`