# Delta-Base Planning

The `difflow.planning` module turns one or more differentiable flowsheets into a
linear (or mixed-integer linear) planning model, keeps it honest with a trust
region, and returns the sensitivity of the *plan* as well as the plan.

It is a module alongside `difflow.eo_solver`, `difflow.estimation` and
`difflow.uncertainty` — **not** a `difflow.plugins` entry point. That registry
is for unit operations.

## Table of Contents

1. [Why delta vectors, and why AD](#why-delta-vectors-and-why-ad)
2. [Quick start](#quick-start)
3. [Blocks, networks and links](#blocks-networks-and-links)
4. [Stating the problem, and the LP that gets solved](#stating-the-problem-and-the-lp-that-gets-solved)
5. [Drawing the flowsheet, the model and the region](#drawing-the-flowsheet-the-model-and-the-region)
6. [The trust-region loop](#the-trust-region-loop)
7. [Scoring: realised violation, not predicted](#scoring-realised-violation-not-predicted)
8. [Bang-bang levers and vertex seeding](#bang-bang-levers-and-vertex-seeding)
9. [Phase boundaries](#phase-boundaries)
10. [Large models: what degrades and what does not](#large-models-what-degrades-and-what-does-not)
11. [Sensitivity of the plan](#sensitivity-of-the-plan)
12. [Modifier adaptation](#modifier-adaptation)
13. [Which delta vectors are wrong: attribution from plant data](#which-delta-vectors-are-wrong-attribution-from-plant-data)
14. [Coefficient covariance and back-off](#coefficient-covariance-and-back-off)
15. [Piecewise-linear blocks and MILP](#piecewise-linear-blocks-and-milp)
16. [Second-order models: is a delta vector enough?](#second-order-models-is-a-delta-vector-enough)
17. [Multi-period planning and inventory](#multi-period-planning-and-inventory)
18. [Solving a quadratic subproblem](#solving-a-quadratic-subproblem)
19. [Feasibility restoration](#feasibility-restoration)
20. [Emitting Pyomo](#emitting-pyomo)
21. [From a flowsheet to a block](#from-a-flowsheet-to-a-block)
22. [A crude unit as a block](#a-crude-unit-as-a-block)
23. [Exporting delta vectors](#exporting-delta-vectors)
24. [What this module is not](#what-this-module-is-not)
25. [API summary](#api-summary)

---

## Why delta vectors, and why AD

Refinery and value-chain planning has been done for decades with linear programs
whose unit submodels are *base plus delta vectors* — a first-order Taylor
expansion

$$y \approx y_0 + J\,(u - u_0).$$

Aspen PIMS, Haverly GRTMPS, Honeywell RPMS and AVEVA Spiral Plan all work this
way; the refereed anchor is Baker and Lasdon, *Successive Linear Programming at
Exxon*, Management Science **31**(3), 1985.

In every one of those systems the delta vectors are generated by **perturbing a
rigorous simulator one variable at a time**. That is $O(n)$ in the number of
decisions, which is why the vectors are refreshed on the order of annually.

A flowsheet is a pure function with its flash, recycle and unit solves
embedded, so `jax.jacobian` returns the *reduced* input-output sensitivity
directly — already implicitly differentiated through those inner solves. That
reduced Jacobian **is** the delta vector. AD builds it in $\min(n_u, n_y)$
passes (forward mode per input, reverse mode per output) against the $2n_u$
evaluations of central differences, and exactly rather than to truncation
error; the gradient of a scalar — the planner's objective — is one reverse
pass, $O(1)$ in the number of decisions.

Measured on the two-plant chain in `difflow.planning.chain` (NGL recovery plus
gas-turbine power, coupled through the residue stream), with a scalar objective
and an $H$-period horizon so that $n = 5H$:

```python
from difflow.planning import scaling_study, format_scaling_table, planner_objective
from difflow.planning.chain import two_plant_chain

def make(n):
    problem = two_plant_chain(horizon=n // 5)
    return planner_objective(problem.planner()), problem.network.decision_start()

print(format_scaling_table(scaling_study(make, [5, 10, 20, 40, 80])))
```

The gradient costs a small constant multiple of one model evaluation — measured
between 0.7x and 2x across $n = 5 \ldots 80$, with no trend in $n$ — while
central differences cost exactly $2n$, tracking theory. (The ratio can dip below
one on a compiled model: XLA prunes parts of the forward program the gradient
does not need, so the two are not strictly nested.) The ratio between them therefore grows linearly in $n$,
and planning (horizon x units x decisions) is precisely where $n$ gets large.
`tests/test_planning.py::TestScaling::test_gradient_cost_ratio_scaling` locks
this in as a regression test.

Reverse mode is what delivers that scaling, and the mode is chosen from the
block's shape rather than hard-coded:

```python
from difflow.planning import choose_ad_mode
choose_ad_mode(n_u=80, n_y=1)    # 'rev'  — many decisions, few outputs
choose_ad_mode(n_u=1,  n_y=40)   # 'fwd'  — one lever, many reported outputs
```

## Quick start

```python
import jax.numpy as jnp
from difflow.planning import Block, Network, DeltaBasePlanner

def ngl_outputs(u):
    # any pure JAX callable: a flowsheet, a unit, an analytic model
    recovery, T_cold, split = u
    ngl_c2 = 12.0 * recovery * jnp.exp(-(T_cold - 218.0) / 60.0)
    residue_F = 90.0 - 20.0 * recovery + 0.4 * (T_cold - 218.0)
    T_colfeed = T_cold + 14.0 * (1.0 - split) + 6.0 * (recovery - 0.5)
    return jnp.array([ngl_c2, residue_F, T_colfeed])

def power_outputs(u):
    fuel_F, alloc = u
    burned = alloc * fuel_F
    power = 0.4 * burned * (1.0 - 0.3 * jnp.exp(-burned / 40.0))
    return jnp.array([power, 0.048 * burned])

ngl = Block(name="ngl", fn=ngl_outputs,
            u_names=["ethane_recovery", "T_coldbox", "split"],
            y_names=["NGL_C2", "residue_F", "T_colfeed"],
            lb=[0.30, 218.0, 0.0], ub=[0.98, 244.0, 1.0])

pwr = Block(name="power", fn=power_outputs,
            u_names=["fuel_F", "alloc"], y_names=["Power", "CO2"],
            lb=[0.0, 0.0], ub=[200.0, 1.0])

net = Network([ngl, pwr], links=[("ngl.residue_F", "power.fuel_F")])

prices = {"ngl.NGL_C2": 9.0, "power.Power": 55.0}      # linear in y: an LP
specs = [("ngl.T_colfeed", "<=", 236.0)]

planner = DeltaBasePlanner(net, prices=prices, specs=specs,
                           radius=0.3)                 # fraction of range

res = planner.solve()
res.plan                              # optimal decisions
res.delta_vectors                     # the J blocks actually used
# res.pyomo_model                     # the emitted Pyomo model (needs pyomo)
res.plan_sensitivity(wrt="prices")    # d(plan)/d(price)
print(res.summary())
```

Every variable is addressed by a qualified name, `"<block>.<variable>"`.

## Blocks, networks and links

A `Block` wraps a callable `u -> y`. It may also take parameters, in which case
its signature is `fn(u, theta)` and the parameters become available to the
sensitivity machinery. Set `jit=True` whenever the block is a real flowsheet:
the planner calls it once per trust-region cycle and once per AD pass, so
compiling it pays for itself immediately (roughly 18x on the reference chain).

A `Network` is a **DAG**. A link makes a downstream input equal to an upstream
output, so it stops being a free decision:

```python
net.decision_names      # ['ngl.ethane_recovery', ..., 'power.alloc']
net.is_linked("power.fuel_F")   # True
```

Blocks are linearised individually and the links become equality rows in the LP.
That is the structure a delta-base model actually has, and at first order it is
equivalent to differentiating the composed chain — LP elimination of the link
rows reproduces the chain rule. Keeping the blocks separate preserves the delta
vectors as an inspectable artefact, which is the thing planners audit.

Recycles **inside** a block are expected, and are exactly where difflow earns
its keep: the tear solve is differentiated implicitly, so the reduced Jacobian
comes back for free. Recycles **between** blocks are rejected with an error that
tells you to merge the loop into one block — that is, into one flowsheet.

## Stating the problem, and the LP that gets solved

Two reports write the problem out from the model that is actually solved, so a
statement cannot drift away from the code.

`DeltaBasePlanner.describe()` answers the three questions that come before any
result: what is being planned (the priced objective), what may be changed to get
it (the free decisions, their bounds, and how far one cycle may move them), and
what may not be violated (the links, the specs, and the acceptance test):

```python
print(planner.describe())
```

`LPModel.as_text()` writes the subproblem out algebraically, row by row. For
blocks `b`, links `s -> t` and specs `k`, with one slack `s_k >= 0` per elastic
spec, each cycle solves

```
max_x    sum_v c_v x_v  -  sum_k pi_k s_k              priced outputs, less slack
s.t.     y_b - J_b u_b  =  y0_b - J_b u0_b             model rows: the delta vectors
         u_t - y_s      =  0                           link rows: the network
         a_k' x - s_k  <=  r_k                         spec rows, elastic
         max(l_i, u0_i - D_i) <= u_i <= min(h_i, u0_i + D_i)   bounds ∩ trust region
         s_k >= 0,   D_i = radius * (h_i - l_i)
```

```python
state = res.state                    # the nonlinear network state at the plan
lp = planner.build_lp(planner.linearize(state), state, radius=0.25)
print(lp.as_text())          # or, after a solve, res.lp_model.as_text()
```

The model rows are the only place the flowsheet enters; eliminating the link
rows reproduces the chain rule. Everything nonlinear — the flash, the recycle,
the efficiency curve — lives outside the LP, in the acceptance test.

## Drawing the flowsheet, the model and the region

A planning model is worth looking at before it is solved. Five drawings, all in
`difflow.planning.diagram` and all exported from `difflow.planning`:

```python
from difflow.planning import (
    draw_chain, draw_planning_network, draw_delta_vectors, draw_taylor_model,
    draw_trust_region,
)

from difflow.planning import linearize_block
from difflow.planning.chain import two_plant_chain, ngl_block

problem = two_plant_chain()              # the reference chain behind draw_chain
chain_state = problem.network.evaluate(problem.network.decision_start())
chain_ngl = ngl_block()

draw_chain(chain_state.as_dict(), prices=problem.prices, specs=problem.specs)
draw_planning_network(problem.network, prices=problem.prices, specs=problem.specs)
draw_delta_vectors(linearize_block(chain_ngl), block=chain_ngl)
draw_taylor_model(chain_ngl, "T_coldbox", "residue_F", radius=0.25)
draw_trust_region(res, decisions=("ngl.ethane_recovery", "ngl.T_coldbox"))
```

| Drawing | Shows |
|---|---|
| `draw_chain` | the reference chain as a process flow diagram: units, streams, the decisions as levers, the priced streams, the spec, and which planning block each unit belongs to |
| `draw_planning_network` | any `Network` as the LP holds it: blocks, free decisions, links, priced outputs, specs |
| `draw_delta_vectors` | one block's `J`, shaded within each row because the rows carry different units |
| `draw_taylor_model` | the delta-vector prediction against the block along one decision, with the trust region marked — the picture of why the region exists |
| `draw_trust_region` | the accept/reject/shrink cycles over the *nonlinear* merit surface |

matplotlib is imported inside these functions, so a headless run pays nothing
for importing the module. `draw_trust_region` evaluates the merit on a grid, so
keep `grid` modest on an expensive flowsheet.

## The trust-region loop

Each cycle is:

1. Linearise every block at the current point — the delta vectors.
2. Solve the LP inside a trust region to get a proposal.
3. **Evaluate the caller's own nonlinear blocks at that proposal.**
4. Accept only if the realised merit improved as much as the LP promised;
   otherwise shrink the region and retry.

Step 3 is not optional. Without it — the "recursion" heuristic that commercial
planning systems use — repeated re-linearisation walks the iterate far outside
the region where any Taylor model is valid, and the LP confidently returns a
plan the real model does not support:

```python
guarded   = DeltaBasePlanner(net, prices=prices, specs=specs, radius=0.3).solve()
unguarded = DeltaBasePlanner(net, prices=prices, specs=specs, radius=0.3,
                             accept_test=False).solve()
guarded.converged      # True  — reaches the optimum
unguarded.converged    # False — oscillates, rho < 0 on alternate steps
```

`accept_test=False` exists so the difference can be measured, not so it can be
used.

The convergence theory is Eason and Biegler, *A trust region filter method for
glass box/black box optimization*, AIChE J **62**(9), 2016,
[doi:10.1002/aic.15325](https://doi.org/10.1002/aic.15325). Their first-order
consistency requirement — that the surrogate match the true model's value *and
gradient* at the trust-region centre — is satisfied **exactly** by an AD Taylor
model. A delta vector obtained by one-at-a-time finite differencing satisfies it
only to truncation error, which is why "recursion" has no comparable guarantee.

The subproblem is an LP, so its optimal value scales with the radius; the
planner uses predicted gain per unit radius as its first-order criticality
measure, confirmed on a probe radius that reuses the current linearisation and
therefore costs no model evaluations. Convergence near a smooth interior optimum
is linear, which is the known behaviour of successive linear programming; the
`radius_min` setting in `TrustRegionOptions` sets the final accuracy.

## Scoring: realised violation, not predicted

Specs are elastic by default: the LP gets a slack variable so it always has a
feasible point to report. But **realised economics are scored against the
nonlinear model**, never against the LP's own prediction:

```python
scored = planner.score(res.decisions)
scored["objective"]         # priced objective from the real blocks
scored["violations"]        # per-spec violation from the real blocks
scored["merit"]             # objective less the violation charge
```

This matters more than it looks. If a plan is scored on its own LP's slacks,
then a planner whose coefficients have gone stale scores *well* precisely by
running off-spec — it predicts feasibility it does not achieve, and is never
charged for the difference. In the prototype this let a static planner beat a
provably optimal oracle.

Back-off (below) is treated as a margin, not a promise: `Spec.violation`
measures against the stated right-hand side, so eating into the margin is not
scored as a violation.

### A block that cannot be evaluated

A block with an inner solve has operating points where the solve has no answer.
A distillation column can dry out, and a flash can lose a phase. The convention
is that such a block returns **NaN**. The planner treats any non-finite value
in the network state (`state_is_finite(state)` is false) as a model that cannot
be evaluated:

- Such a point scores merit `-inf` and violation `+inf` (`planner.score(u)["evaluable"]`
  is `False`).
- A proposal there is **rejected** and the radius shrinks. This holds even with
  `accept_test=False`, because a point that cannot be evaluated cannot be linearised.
- Restoration never takes such a point as its least-violating one.
- A start point that cannot be evaluated is an attempt with reason
  `"start_not_evaluable"`, and any evaluable seed beats it. If no start can be
  evaluated, `solve()` raises.

This is written down rather than left to IEEE arithmetic. A NaN *merit* already
fails `rho >= eta_accept`, but the other two cases do not fail on their own:

- A NaN *spec output* scores **zero** violation, because Python's `max(0.0, nan)`
  is `0.0`. Restoration would then take a failed solve for a feasible point.
- A NaN output that is neither priced nor in a spec leaves the merit finite. The
  main loop would accept the point and then try to linearise at it.

The tests are in `tests/test_planning_nonfinite.py`. The crude unit is the worked
case: [Planning with the crude unit](unit-operations-refinery.md#planning-with-the-crude-unit)
drives a column into the region where it does not converge and shows the plan
backing off.

## Bang-bang levers and vertex seeding

Levers like ethane recovery versus rejection sit at a bound and switch
discretely with prices. A single interior start converges to whichever corner it
happens to face, so the planner also seeds from bound vertices of the most
price-sensitive decisions. Ranking them costs one reverse-mode AD pass whatever
$n$ is:

```python
planner.vertex_seeds()        # the extra starting points
res.n_starts                  # how many were run
res.attempts                  # one PlanResult per start
```

Because the derivative of a bang-bang plan is zero — the corner does not move
under an infinitesimal price change, it *switches* at a finite one — the useful
question is where it switches:

```python
from difflow.planning import price_switch_point
price_switch_point(planner, "power.Power", 5.0, 60.0, tol=1.0)["price"]
```

## Phase boundaries

A delta vector computed across phase appearance or disappearance is
**meaningless**, not merely inaccurate: the function is not differentiable
there, and the Taylor model extrapolates a branch that has ceased to exist.
difflow's flash already guards its single-phase branch with safe inputs, so the
numbers keep coming — which is exactly why the planner has to say something.

Give a block a `phase_fn` and the planner warns when a proposal crosses a
regime:

<!-- doc-test: skip: template; `info` is the flash result computed inside the user's own block -->
```python
Block(..., phase_fn=lambda u, th: jnp.atleast_1d(info["V_frac"]),
      phase_names=("V_frac",), phase_bounds=(0.0, 1.0))
```

```text
PhaseBoundaryWarning: block 'ngl': indicator 'V_frac' crossed a phase boundary
between the linearisation point and the proposal (0 -> 0.4043, regime 0 -> 1).
The delta vectors either side describe different functions, so this
linearisation is not valid at the proposal — shrink the trust region or
re-centre inside one regime.
```

The messages are also collected on `res.phase_warnings`.

## Large models: what degrades and what does not

> **Does a very large flowsheet suffer "gradient collapse" — the vanishing-gradient
> failure familiar from RNNs and LSTMs?**
> No, and the reason is worth being precise about, because three *other* things
> do degrade with size and they have different remedies.

### The non-problem: small absolute sensitivities

Compose enough units and the sensitivity of a downstream product to an upstream
lever becomes numerically tiny. That is almost always **physics, not numerical
loss**. Run a chain of first-order stages deep enough to drive the outlet to
`1e-61` and the *relative* sensitivity is still exact:

| depth | `y0` | `dy/dk_1` | `dln y/dln k_1` | exact |
|---|---|---|---|---|
| 10  | 9.766e-04 | -4.883e-04 | -0.500000 | -0.5 |
| 100 | 7.889e-31 | -3.944e-31 | -0.500000 | -0.5 |
| 200 | 6.223e-61 | -3.112e-61 | -0.500000 | -0.5 |

Nothing was lost: `y` itself is `1e-61`. A lever twenty units upstream of a
product genuinely has little absolute leverage on it, and an LP that gives it
little weight is right. The LSTM analogy breaks on three counts — `difflow`
enables `float64` at import (~300 orders of headroom, not `float32`'s 38), a
flowsheet is tens of units deep rather than thousands of timesteps, and the
trust-region loop re-linearises every cycle instead of accumulating a product
over many gradient steps.

`composed_sensitivity` measures this directly, by **one AD pass over the
composed network** rather than by multiplying per-block Jacobians together, so
it does not itself suffer the rounding it is measuring:

```python
from difflow.planning import composed_sensitivity

S, outputs, decisions = composed_sensitivity(net)
# S[i, j] = fractional change in output i per full-bound-range move of decision j
```

An *exactly* zero entry is the meaningful one — it says no path from that
decision to that output survived the linearisation.

### Why the block decomposition helps

Linearising the whole plant as one `Block` does form the deep chain-rule
product, and its entries do collapse. Keeping it as linked blocks does not.
Measured on a chain of identical stages in mixed engineering units:

| depth | `cond(A_eq)` | worst block `cond(J)` | monolithic `min \|J\|` |
|---|---|---|---|
| 4   | 4.58e+02 | 319.6 | 1.43e+01 |
| 16  | 8.48e+02 | 319.6 | 4.61e-02 |
| 32  | 9.04e+02 | 319.6 | 2.20e-05 |
| 64  | 9.19e+02 | 319.6 | 5.00e-12 |
| 128 | 9.23e+02 | 319.6 | 2.58e-25 |

The LP's condition number **saturates** while the monolithic Jacobian collapses
to `1e-25`. Each block is linearised where the *network* puts it, so its
entries stay `O(1)` in its own units no matter how deep it sits; composition
lives in the link equality rows, where the solver eliminates with pivoting
rather than you forming a 128-term floating-point product. This is the
practical payoff of the rule in `Network` that inter-block recycles are
rejected and loops belong *inside* a block — but note the tension: a plant-wide
recycle forces you toward a bigger block, which is the column of that table
where sensitivities do collapse.

### The three things that actually degrade

`check_delta_health` reports all three. Nothing here raises during a solve; the
point is to make the failure visible before it is silently planned around.

```python
from difflow.planning import check_delta_health

report = check_delta_health(net)          # or a single Block
print(report.summary())
report.raise_on_error()                   # non-finite deltas only
report.warn()                             # as DeltaHealthWarning
```

or, including the assembled program's coefficient scaling:

```python
planner.check_health().summary()
```

**1. Dead levers — the one genuine analogue of lost information.** Every
`jnp.clip`, `jnp.minimum` and `jnp.where` sitting on an active spec contributes
an *exactly* zero column to `J`. The LP then correctly concludes the lever does
nothing and never moves it — not because it does not matter, but because the
linearisation cannot see that it does. This scales with the number of active
quality specs, which is to say it scales with model size.

```text
[warning] dead_lever: blend.butane: delta column is structurally zero (max
scaled sensitivity 0.000e+00); the LP will never move this lever. The lever is
interior, so the zero comes from a saturated expression downstream of it — a
clip, minimum or where on an active spec — rather than from the bound.
Re-centre inside the smooth region, or model the saturation explicitly with a
piecewise block rather than letting the clip hide it.
```

Companion findings: `dead_output` (an output that responds to nothing, so a
price or spec on it is being applied to a constant) and `no_influence` (a
decision that reaches no output anywhere in the network).

**2. Amplification, not attenuation.** A tear solve differentiated implicitly
returns `(I - A)^-1`, so a recycle of loop gain `g` multiplies sensitivities by
`1/(1 - g)`:

| loop gain | `x*` | `dx/dfeed` |
|---|---|---|
| 0.900 | 9.1743  | 9.1743  |
| 0.990 | 50.2513 | 50.2513 |
| 0.999 | 90.9918 | 90.9918 |

Recycle-to-extinction loops push `g` toward one, the delta vectors blow up, and
the trust region has to shrink to stay honest. **Large flowsheets fail by
exploding sensitivities far more often than by vanishing ones.** The
`amplifying` finding fires when a full trust-region step is predicted to change
an output by more than `AMPLIFY_TOL` times its own value, and distinguishes a
near-unity loop gain from the other common cause — bounds set far wider than
the range the lever is actually planned over. Outputs whose value is *zero* at
the linearisation point are excluded, because "a fraction of its own value" is
undefined there and a bang-bang lever at a corner routinely produces one.

**3. Scale spread.** Mixed engineering units — ppm against kbbl/d against
$/bbl — put entries spanning many orders of magnitude in one constraint matrix.
This is the mundane failure that actually stops a large planning model, and it
is a *units* problem, not a gradient problem. `check_lp_scaling` flags both the
whole matrix and individual rows whose small coefficients sit below the
solver's effective precision relative to their large ones. Nondimensionalise
the levers before blaming the gradients.

### What it finds on the reference chain

`two_plant_chain()` is small, and its one finding is the units problem:

```text
>>> problem.planner(radius=0.25).check_health().summary()
delta-vector health: 1 findings (0 error, 1 warning)
  [warning] scale_spread: network: constraint entries span 2.812e-08 to
  8.044e+01 (ratio 2.860e+09); the solver's pivot tolerances are being asked to
  separate signal from unit conversion. Re-scale the offending variables to
  comparable magnitudes.
```

`ngl.P_expander` is in pascals and operates at 2.5e6, so `dE_refrig/dP` is
`3.6e-08` per Pa, while `power.alloc` is dimensionless and gives
`d gas_sold/d alloc = 85.8`. Nine orders of magnitude in one matrix, entirely
from the choice of unit. Expressing the expander pressure in MPa removes it.

Nothing is dead and nothing amplifies: the composed sensitivities span
`0` to `2.2`, the deepest lever (`ngl.split`) still reaches an output at
`5.8e-02`, and both blocks are well conditioned. That is what a healthy model
looks like — the diagnostics are quiet until size makes them speak.

Every diagnostic works on the **scaled** Jacobian,
`Js[i, j] = J[i, j] * u_scale[j] / y_scale[i]` — the fractional change in output
`i` per full-bound-range move of input `j`. Comparing raw `J` entries across a
model in mixed units compares unit conversions. Linked inputs are scaled by
their value at the operating point instead, because the trust region never
steps them: the LP's link row ties them to the upstream output.

## Sensitivity of the plan

Because the blocks are differentiable, the planner returns the sensitivity of
the **plan itself**, not just the plan. This is what k_aug and sIPOPT provide
for a single NLP, and it is what makes capital planning tractable without
scenario enumeration. A commercial planning system structurally cannot provide
it, because its submodels are not differentiable.

```python
s = res.plan_sensitivity(wrt="prices")   # or wrt="theta"
s.d_plan                                 # d(u*)/d(theta)
s.d_objective                            # d(objective*)/d(theta)
s.multipliers                            # per active constraint
print(s.summary())
```

The machinery is the implicit function theorem on the KKT conditions at the
converged plan. With free decisions $u_F$ (those not at a bound), active
constraints $h_A$ and multipliers $\nu$:

$$\begin{bmatrix} \nabla^2_{uu} L & A^{\mathsf T} \\ A & 0 \end{bmatrix}
\begin{bmatrix} \mathrm{d}u_F/\mathrm{d}\theta \\ \mathrm{d}\nu/\mathrm{d}\theta \end{bmatrix}
= \begin{bmatrix} -\partial^2 L/\partial u\,\partial\theta \\ -\partial h_A/\partial\theta \end{bmatrix}$$

and the envelope theorem gives the objective sensitivity directly,

$$\frac{\mathrm{d}\phi^*}{\mathrm{d}\theta}
= \frac{\partial \phi}{\partial \theta} + \nu^{\mathsf T}\frac{\partial h_A}{\partial \theta},$$

which is exact even at a vertex where $\mathrm{d}u/\mathrm{d}\theta$ is zero.
Every derivative in that system comes from AD on your own blocks.

Rows for decisions pinned at a bound are exactly zero, and the result says so
in `s.note` rather than letting you read a vertex as if it were an interior
solution. A rank-deficient KKT matrix is reported through `s.degenerate`, and a
plan that is not actually stationary is called out in `s.note` too.

## Modifier adaptation

Under **structural** plant-model mismatch — the model has the wrong form, not
merely the wrong parameters — correcting only predicted *values* converges to
the model's optimum, not the plant's. Matching the plant's *gradients* fixes it
(Marchetti, Chachuat and Bonvin, I&ECR **48**(13), 2009,
[doi:10.1021/ie801352x](https://doi.org/10.1021/ie801352x)):

$$y_{\text{mod}}(u) = y_{\text{model}}(u) + \varepsilon + \lambda\,(u - u_{\text{ad}}).$$

In a delta-base LP the delta vectors **are** those gradients, so $\lambda$ is
added straight onto $J$ and the LP needs no change:

```python
from difflow.planning import run_modifier_adaptation

# The "plant": same inputs and outputs as the model block, but a different form
def real_ngl(u):
    return ngl_outputs(u) * jnp.array([0.85, 1.0, 1.0]) + jnp.array([0.0, 2.0, 0.0])

plant = {"ngl": real_ngl}
ma_planner = DeltaBasePlanner(net, prices=prices, specs=specs, radius=0.3)
ma = run_modifier_adaptation(ma_planner, plant, max_iter=6)
ma.history[-1]["plant_objective"]
print(ma.plan)

# The comparison the method exists for:
ma0_planner = DeltaBasePlanner(net, prices=prices, specs=specs, radius=0.3)
run_modifier_adaptation(ma0_planner, plant, use_gradients=False, max_iter=6)   # model's optimum
```

Both modifiers are first-order filtered, because an unfiltered gradient
correction from noisy plant data is a fast route to oscillation. When modifiers
are in force the *corrected* model is the planner's own model, so that — and not
the uncorrected blocks, and never the plant — is what the acceptance test is
judged against.

## Which delta vectors are wrong: attribution from plant data

`update_modifiers` takes the plant gradient from a callable. A running plant
is not a callable; what exists is a history of each block's inputs and some
of its measured outputs. `attribute_deltas` estimates the modifiers from that
history and, as importantly, says which of them the history can support:

```python
from difflow.planning import attribute_deltas

import numpy as np
from difflow.planning import attribute_deltas

# Routine plant history for the "ngl" block: inputs wander a little, and the
# plant sits 0.5 mol/s above what the model predicts for NGL_C2.
rng = np.random.default_rng(0)
n = 60
times = np.arange(n, dtype=float)
U = np.asarray(ngl.u0) + rng.normal(size=(n, 3)) * [0.03, 1.0, 0.05]
Y = np.array([ngl_outputs(jnp.asarray(u)) for u in U])
c2_meas = Y[:, 0] + 0.5 + rng.normal(scale=0.05, size=n)
resid_meas = Y[:, 1] * (1.0 + rng.normal(scale=0.01, size=n))

att = attribute_deltas(ngl, U, {"NGL_C2": c2_meas, "residue_F": resid_meas},
                       sigma_y={"NGL_C2": 0.05, "residue_F": 0.01},
                       t=times,
                       move={"ethane_recovery": 0.1, "T_coldbox": 2.0, "split": 0.1},
                       sigma_u={"T_coldbox": 0.1},        # errors in variables
                       log_outputs={"residue_F": 1e-3})   # relative errors
print(att.table())
ma_planner.modifiers[ngl.name] = att.to_modifiers()       # flagged terms only
att.exposure(res)           # level error x shadow price of its model row
```

For each measured output the residual $r = g(y_{\text{meas}}) -
g(y_{\text{model}}(u))$ is fitted by weighted least squares on a level, a
trend, and one slope per input scaled by that input's characteristic move.
Routine plant data are not a designed experiment, and the fit is built around
that:

- **Estimability is decided from the design, not the answer.** The scaled
  slope columns, with level and trend projected out, go through a
  column-pivoted QR; a slope is estimated only while
  $|R_{kk}|\cdot\text{materiality} \ge 1.9$. An input the operators held
  still is reported `not estimable` and left out. Expect most slopes to land
  there; the level is the reliable part.
- **Aliases are reported.** Inputs that moved together are only estimable as
  a combination; a significant estimate that absorbs a held-out input with
  $|A| > 0.3$ is reported as `combination`, not as a finding about one input.
- **Standard errors are inflated** by $\sqrt{\phi(1+\rho)/(1-\rho)}$, with
  $\phi = \max(1, \chi^2/\text{dof})$ and $\rho$ the lag-1 autocorrelation of
  the residual. Without it, slow drift produces a stream of false flags.
- **A Picard check separates a wrong delta from a wrong form.** Residual
  weight along a direction the design barely resolves would need an absurd
  slope to explain; the output is marked `structural`, and no affine modifier
  is the fix.

A term is flagged when it is estimable, $|z| > 3$ after inflation, and larger
than `materiality` (default: one `sigma_y`). The level refers to the latest
time and to `u_ref`, by default the mean input, where it is not aliased with
any held-out slope. Outputs that were not measured are listed in
`res.unobserved`: the data say nothing about them.

## Coefficient covariance and back-off

Delta vectors are functions of the model parameters. When those parameters come
from a reconciliation (`difflow.reconciliation`) or an estimation
(`difflow.estimation`) they carry a covariance, and that covariance propagates
through the block Jacobian onto the LP's coefficients — and hence onto the
constraint values the plan is built to respect.

Stale or uncertain parameters usually cost **feasibility**, not optimality: a
plan that predicts it sits exactly on a spec will spend roughly half its periods
on the wrong side of it. The remedy is a back-off sized by the propagated
uncertainty, $\kappa\sigma$ with $\sigma^2 = g^{\mathsf T}\Sigma_\theta g$:

```python
import jax.numpy as jnp
from difflow.planning import constraint_backoff, apply_backoff

# The blocks carry parameters (theta); here the reference chain's NGL block.
chain_planner = problem.planner(radius=0.25)
theta_covariance = jnp.diag(jnp.array([0.05, 0.03]) ** 2)   # from a fit

found = constraint_backoff(chain_planner, problem.network.decision_start(),
                           theta_covariance, ["ngl.feed_scale", "ngl.colfeed_rise"],
                           kappa=2.0)
print(found.summary())
chain_planner.specs = apply_backoff(chain_planner.specs, found)
```

This is a thin layer over `difflow.uncertainty.propagate_covariance`. No
commercial planning system does it.

## Piecewise-linear blocks and MILP

When a block's response to one lever is strongly curved over the whole operating
range, a piecewise model captures the range at once and turns the plan into a
MILP instead of a sequence of trust-region LPs. Building one is a batching
problem: `vmap` the block over every breakpoint in a single call, and `vmap` its
Jacobian too.

```python
from difflow.planning import PiecewiseSpec

pw_planner = DeltaBasePlanner(net, prices=prices, specs=specs, radius=0.1,
                              piecewise=[PiecewiseSpec("ngl", "T_coldbox",
                                                       n_points=9)])
pw_res = pw_planner.solve()
pw_res.lp_model.integer_cols     # the SOS2 interval binaries
pw_res.lp_model.sos2_sets        # emitted natively when you export to Pyomo
```

The response to the distinguished variable is piecewise-linear and *exact* at
the breakpoints; the response to the other inputs uses a single Jacobian taken
at the centre. That is exact when the block is separable and a first-order
approximation of the cross terms otherwise. The alternative would multiply the
SOS2 weights by the other inputs, which is bilinear — and bilinear is the line
this module does not cross.

## Second-order models: is a delta vector enough?

A delta vector is a first-order model. AD supplies the second order for about
what the first order costs: a Hessian-vector product is one `jax.jvp` through
`jax.grad`, so the exact Hessian of one scalar output costs $O(n_u)$ HVPs —
roughly what a *central-difference Jacobian alone* costs the systems that build
delta vectors by perturbation. The second-order model is available at the price
the incumbent already pays for its first-order one.

Where it helps, it helps enormously. Stepping a reduced AC power-flow model the
whole way from an incumbent generator schedule to the optimal one
(`tests/power/test_planning_opf.py`):

| step | true cost | linear error | quadratic error |
|---|---|---|---|
| 25% | 5379.35 | −7.86 | 0.00 |
| 50% | 5336.09 | −31.44 | +0.02 |
| 100% | 5296.69 | −125.73 | **+0.13** |

Three orders of magnitude, over a step no trust region would allow in one
cycle.

Whether to *use* it is a different question, and `difflow.planning.curvature`
answers it rather than assuming. Every guarantee in this module chains off the
subproblem being an LP: solved to global optimality, duals readable as prices,
Eason–Biegler filter convergence. A quadratic objective preserves those only
while its Hessian is definite in the direction of optimisation. An indefinite
one makes the subproblem a nonconvex QP and voids all three **while still
returning a number**.

```python
from difflow.planning import check_model_order

block = ngl      # any Block; here the quick-start NGL block
rep = check_model_order(block, "residue_F", radius=0.2, sense="min")
print(rep.summary())
rep.recommended        # 'linear' or 'quadratic'
rep.improvement        # error-reduction factor, inf when exact
rep.curvature          # Hessian, eigenvalues, definiteness verdict
rep.caveat             # why a good fit was still refused, or None
```

`check_model_order` recommends `"quadratic"` only when the model is *both*
materially more accurate over the step and convex in the direction of
optimisation. A perfect fit with an indefinite Hessian is refused, and
`rep.caveat` says what to do about it — Gauss-Newton, modified Cholesky, or a
damped BFGS update, which needs no second derivatives at all.

**Definiteness is a property of where you are, not of the model.** On the same
nine-bus network the reduced cost Hessian is positive definite at the incumbent
operating point and strongly indefinite at heavily loaded ones, and the
Hessians of the voltage and thermal limits are indefinite at nearly every point
sampled. So the check belongs at the linearisation point each cycle, not once
at commissioning. `test_definiteness_is_a_property_of_the_point` is the
regression.

Pass a mapping to take the curvature of the *priced objective* rather than of
one output:

```python
from difflow.planning import block_curvature
curv = block_curvature(block, {"NGL_C2": 12.0, "residue_F": -1.0})
curv.convex_for("max")
```

Nothing here changes the planner. These are diagnostics on a block, in the
spirit of `check_delta_vectors` and `check_delta_health`: they say what a
second-order subproblem would buy and what it would cost in guarantees, so the
decision to build one is taken on evidence.

## Multi-period planning and inventory

Periods are *replicated* by `two_plant_chain(horizon=n)`, which names blocks
`ngl@t0 … ngl@t3` and couples them only through a shared cap. That is a horizon
built to make the AD scaling argument measurable. A planning model couples
periods through **inventory**: what is not sold this period is still there next
period.

That needs no new machinery. A `Link` is output-to-input and the network rejects
only *cycles*, so a forward link from one period's tank level to the next is an
ordinary DAG edge:

```python
links = []
t = 1
links.append((f"tank@t{t-1}.level_out", f"tank@t{t}.level_in"))
```

`tests/test_planning_multiperiod.py` builds a four-period storage-arbitrage
model this way — make cheaply, hold, sell into a price spike — and it plans
correctly, holding inventory back and drawing the tank down into the spike.
The opening level is not constrained to zero; the first period simply uses a
block with no `level_in` input, so the model *starts* feasible.

Two traps, both found by building the model rather than by reading the code,
and both pinned as regressions.

**`Spec` is elastic by default, and that is wrong for a mass balance.** Elastic
slack is right for a commercial specification — you can ship off-spec at a cost
— and fiction for a physical one. With elastic inventory constraints the
planner reports a *higher* objective than the feasible plan earns, by running
the tank negative and selling from an empty vessel, and it converges and
reports that number without complaint. Physical balances must be
`elastic=False`.

**There is no feasibility restoration.** If an inelastic spec is violated at the
starting point, the LP is infeasible from the first cycle, and the planner's
response — shrink the radius — can only tighten it. The run ends at
`reason="lp_infeasible"` with every decision still on its start value. The
failure is reported rather than hidden, but the remedy is to start feasible;
a phase-1 restoration step is what would fix it properly.

## Solving a quadratic subproblem

Measuring the curvature is one thing; using it is another. `model_order`
switches the subproblem from an LP to a QP:

```python
qp_planner = DeltaBasePlanner(net, prices=prices, specs=specs, sense="max",
                              model_order="quadratic")
```

| | | |
|---|---|---|
| `"linear"` | delta vectors alone | the default, unchanged |
| `"quadratic"` | curvature everywhere, convexified where it points the wrong way | fewest iterations |
| `"auto"` | curvature only where it is already definite | the exact second-order model, or nothing |

The model stays a **QP, not a QCQP**. The objective is linear in the block
outputs, so the priced combination of one block's outputs has a single Hessian
and the whole correction collapses to one term per block:

$$p^{\mathsf T} y \approx p^{\mathsf T} y_0 + p^{\mathsf T} J\,\delta u
  + \tfrac12 \delta u^{\mathsf T} H_p\, \delta u, \qquad
  H_p = \sum_i p_i H_i.$$

Making the *rows* quadratic would turn each subproblem into a nonconvex QCQP,
and the subproblem solving to global optimality is what every guarantee here
rests on. Constraint rows stay first order.

### What it buys: termination

On the AC-OPF comparison (`tests/power/test_planning_opf.py`), from the
incumbent generator schedule on case9:

| `model_order` | AC cost | iterations | ended on |
|---|---|---|---|
| `"linear"` | 5296.69 | 40 | the iteration cap |
| `"auto"` | 5296.69 | 19 | its own radius test |
| `"quadratic"` | 5296.69 | 12 | its own radius test |

Same optimum, same AC feasibility. What changes is that the loop *stops*. A
first-order model of a curved objective promises a gain the blocks do not
deliver, so the radius ratchets down and the run exhausts its budget sitting
on the right answer without being able to say so. Iterations are the metric
that matters, not wall time: each one is an evaluation of the caller's blocks,
which on a real flowsheet is the entire cost.

### Convexification, and why it is reported

`convexify` clips the eigenvalues that point the wrong way for the sense being
solved and keeps the rest, so the directions the model curves in are
preserved and only the amount changes. A convexified model is **not** the true
second-order model. What keeps it honest is what keeps the first-order model
honest — the trust region, and an acceptance test against the caller's own
nonlinear blocks — and `qp.convexification` records every block it touched:

```python
lins = qp_planner.linearize(state)
lp, qp = qp_planner.build_subproblem(lins, state, radius=0.25)
qp.convexified                       # did anything have to be clipped?
for name, report in qp.convexification.items():
    print(name, report.summary())    # e.g. '3 of 5 eigenvalues clipped (worst -1.28e+03)'
```

Use `"auto"` to refuse the compromise instead: it takes curvature only where
the Hessian is already definite and falls back to the LP elsewhere, so the
model solved is always the real one.

Two limits worth stating. The QP is warm-started from the LP and falls back to
it whenever it fails to improve, so a quadratic subproblem can only match or
beat a linear one — that is what makes it safe to switch on, and it also means
the LP is still built and still solved every cycle. And a block with integer
columns (a piecewise SOS2 model) would make it a MIQP, which is out of scope:
those networks fall back to `"linear"` automatically.

## Feasibility restoration

An elastic spec cannot make the subproblem infeasible — that is what the slack
is for. An **inelastic** one can, and inelastic is the right setting for
anything physical.

Before restoration the loop's only response to an infeasible LP was to shrink
the radius, and shrinking a box that already excludes the feasible set
excludes it harder. The run ground down to `radius_min` and reported
`lp_infeasible` with every decision still on its starting value.

Restoration solves a phase-one problem instead — stop optimising, minimise
infeasibility:

$$\min\; \textstyle\sum a \quad\text{s.t.}\quad
  A_{eq} x = b_{eq},\;\; A_{ub} x - a \le b_{ub},\;\;
  l \le x \le u,\;\; a \ge 0.$$

Three choices in that statement are deliberate.

**The equality rows are not relaxed.** They are the model rows
($y = y_0 + J(u-u_0)$) and the link rows, and both are *definitional*: given
$u$ inside its bounds there is always a $y$ that satisfies them. An
equality-infeasible subproblem means the model itself is broken — contradictory
links, a degenerate block — and relaxing it would bury that under an artificial
variable. Only the spec rows, which are the caller's requirements rather than
the model's structure, get artificials.

**Restoration has its own trust region and its own acceptance test.** A
phase-one LP handed a big enough region will happily propose a point it
*predicts* feasible and the blocks are not: on the storage model this was built
against, doubling the radius took predicted violation to zero while the true
violation rose from 2.5 to 2.4. So a restoration step is kept only when the
caller's own blocks report less violation, and the region shrinks on one that
is not — the same rule the main acceptance test follows. The feasible point is
usually outside the initial box, and it is reached the way a trust-region
method reaches anything distant: as a sequence of accepted steps, re-centring
each time.

```
restoration 1: violation 2.500e+00 -> 1.580e+00 | accepted
restoration 2: violation 1.580e+00 -> 9.600e-02 | accepted
restoration 3: violation 9.600e-02 -> 1.393e+00 | rejected, shrink
restoration 4: violation 9.600e-02 -> 0.000e+00 | accepted
```

The optimisation radius is not the restoration radius. Searching for a
feasible point can leave the working region very small, and that smallness is
a fact about the search rather than about where the objective model is
trustworthy, so restoration hands back the radius it was called with.

**Anything that is not one of those two row kinds is an error.** `A_ub` holds
the caller's spec rows and the SOS2 adjacency rows of a piecewise block, and
restoration picks the first by a *positive* match on the row label — so
silence is the default for anything new. A third row kind that was the
caller's to relax would quietly be left hard, phase one would be unable to buy
down the violation on it, and the planner would report `restoration_failed` on
a recoverable problem with nothing to say why. So the taxonomy is asserted
total: an unrecognised row label raises, naming the row. Adding a row kind to
`assemble.py` means declaring which side it is on, in
`RELAXABLE_PREFIXES` or `STRUCTURAL_PREFIXES`.

Set `TrustRegionOptions(max_restoration=0)` to get the old behaviour; the
history records every restoration cycle with `Iteration.restoration` set, so
they are visible in the audit trail rather than hidden inside a solve.

## Emitting Pyomo

difflow is not short of solvers — `difflow.eo_solver` solves a flowsheet's
equations simultaneously with Newton, and the LP here is solved with HiGHS
through `scipy.optimize.linprog`. What `difflow.planning` does not own is a
*mathematical programming* stack: no branch-and-bound, no interior point, no
modelling language. `LPModel.to_pyomo()` emits a `ConcreteModel` so the plan
composes with the existing Pyomo/IDAES ecosystem instead of competing with it:

<!-- doc-test: skip: needs the optional pyomo package and a cbc solver binary -->
```python
model = res.pyomo_model          # or res.lp_model.to_pyomo()
import pyomo.environ as pyo
pyo.SolverFactory("cbc").solve(model)
```

Pyomo is an optional dependency (`pip install "difflow[planning]"`). Without it
the built-in HiGHS solve via SciPy is used for everything; only `to_pyomo()`
requires it.

The relation to the equation-oriented solver is worth stating, since the two are
cousins. Both assemble every unit's equations into one system rather than
marching unit by unit: `eo_solver` solves the *nonlinear* residuals `F(x) = 0`,
while the planning LP is the linearised, priced, bounded version of the same
assembly — the delta vectors are its model rows and the links its connectivity
rows, with an objective and specs a simulation does not have. An EO solve asks
"what does the plant do at these inputs?"; the LP asks "which inputs pay best?";
the trust-region loop alternates between the two, which is why its acceptance
test is a nonlinear evaluation.

## From a flowsheet to a block

Every example so far hand-wrote its `Block.fn`. For a real flowsheet you do not
have to: `Block.from_flowsheet` builds the callable, and the whole path stays
pure JAX, so `jax.jacobian` of it is the reduced input-output sensitivity of the
plant — implicitly differentiated through the recycle tear solve and every inner
unit solve.

```python
import jax.numpy as jnp
from difflow import (CSTR, CSTRParams, Flowsheet, IdealThermo, Mixer,
                     SpeciesData, Splitter, Unit, make_stream)
from difflow.planning import Block, check_delta_vectors

# A small flowsheet with a recycle: feed + recycle -> mixer -> CSTR -> splitter
thermo = IdealThermo({
    s: SpeciesData(s, MW=100.0, Cp_coeffs=(75.0, 0.0, 0.0, 0.0),
                   Hvap_coeffs=(35000.0, 0.38, 500.0),
                   antoine_coeffs=(10.0, 3000.0, -50.0))
    for s in ("A", "B")})

def rate_fn(C, T, params):
    return jnp.array([params["A"] * jnp.exp(-params["Ea"] / (8.314 * T)) * C["A"]])

cstr = CSTR(CSTRParams(V=jnp.array(1.5), rate_fn=rate_fn,
                       stoich=jnp.array([[-1.0], [+1.0]]),
                       rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
                       species_order=["A", "B"]),
            thermo=thermo, mode="isothermal")

fs = Flowsheet(species_order=["A", "B"], default_flow=1.0)
fs.add_feed("feed", make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0))
fs.add_unit(Unit("mix", Mixer(["A", "B"]), ["feed", "recycle"], ["mixed"]))
fs.add_unit(Unit("reactor", cstr, ["mixed"], ["product"], params={"T_spec": 350.0}))
fs.add_unit(Unit("split", Splitter(["A", "B"]), ["product"], ["purge", "recycled"],
                 params={"split_frac": 0.5}))
fs.add_recycle("recycled", "recycle")

blk = Block.from_flowsheet(
    fs,
    u=["reactor.V", "feed:feed.total_flow"],   # levers
    y=["purge.F_B", "purge.total_flow"],       # outputs
    name="plant", lb=[0.5, 5.0], ub=[5.0, 20.0],
    solve_kwargs={"tol": 1e-10, "max_iter": 200})

check_delta_vectors(blk)["passed"]             # AD vs central differences
```

Lever keys are `Flowsheet._apply_params` notation: `"<unit>.<param>"` for a unit
parameter, and `"feed:<stream>.<field>"` for a feed, where the field is `T`,
`P`, `total_flow`, `F_<species>` or `x_<species>`. Feed rate and composition are
the most common planning levers there are, so they are first-class: `total_flow`
scales the stream at constant composition, `x_<species>` moves a mole fraction
at constant total, and every `x_` in one call is applied together.

Output keys are `"<stream>.<quantity>"` with the same quantity vocabulary.
`u0` defaults to the flowsheet's *current* values, which is almost always the
base case you want, and the block's `metadata` records the original keys and the
units of every variable so the export below is self-describing.

Two things worth knowing. A flowsheet with a recycle is solved by fixed-point
iteration, and the Anderson/Wegstein accelerators are Python loops that cannot
be traced — so under `jax.jacobian` the solve routes automatically to the
optimistix fixed-point path, which carries an implicit-differentiation rule.
And `check_delta_vectors` is worth running before anything leaves the building:
it is `2 n_u` extra model evaluations against a Jacobian that costs `min(n_u, n_y)` AD passes, and
it is the cheapest way to find out that a lever does nothing.

## A crude unit as a block

`difflow_refinery.planning.cdu_block(unit, levers, outputs, rate=, T=, P=)` is a
ready-made block for a rigorous atmospheric crude column. Its levers and outputs
are named by meaning, for example `naphtha.yield`, `pa1.duty`, `crude.rate`,
`kero.tbp95`, `gap.kero_diesel` and `furnace.fired`. They are carried in planner
units: bbl/d, MW, °C for temperatures, K for temperature differences, and kg/h
for steam. The units are recorded in `metadata["u_units"]` and `["y_units"]`.

`product_value_block` and `link_cdu` give the downstream half of a
CDU → product value network. On the 30-stage test column:

- `check_delta_vectors` passes with an error of 5.6e-10 relative to the largest entry.
- `check_delta_health` is clean on the default outputs.
- A four-lever plan (crude rate, two yields, overflash) under a kero end-point
  spec terminates `stationary` in 6 iterations.
- A fresh column solve at that plan reproduces the planner's state to 4e-15.

The column has operating points where it does not converge, and the block
returns NaN there; see [A block that cannot be evaluated](#a-block-that-cannot-be-evaluated).
The details, including why yields are the levers and cut points the outputs, are in
[Planning with the crude unit](unit-operations-refinery.md#planning-with-the-crude-unit).
Example: `examples/35_refinery_cdu_planning.ipynb`.

## Exporting delta vectors

The person who owns the planning model is usually not the person who owns the
simulation. They want a table of coefficients and a base case, not a JAX
install. `difflow.planning.export` is that hand-off: one IR, `DeltaVectorSet`,
and several renderers over it — the same "structured IR plus renderers" shape as
`difflow.report`.

```python
from difflow.planning import (DeltaVectorSet, write_json, write_csv,
                              write_lp, write_mps, write_iterations_csv)

res = DeltaBasePlanner(net, prices=prices, specs=specs).solve()
dvs = DeltaVectorSet.from_result(res)

write_json(dvs, "plan.json")            # the lossless manifest
write_csv(dvs, "tables/")               # one shift-vector table per block
write_iterations_csv(res, "iters.csv")  # the trust-region audit trail
```

<!-- doc-test: skip: write_lp / write_mps go through Pyomo, an optional dependency -->
```python
write_mps(res.lp_model, "plan.mps")     # the assembled LP itself
```

The IR exists because the numbers alone are not enough. A Jacobian entry means
nothing without the base case it was taken at, the names and units of its rows
and columns, and the radius over which the first-order model was ever meant to
hold. `Linearization` carries the numbers, `Block` carries the names, bounds and
units, `PlanResult.radius` carries the validity — `from_result` joins the three,
and adds the links, prices, specs, LP duals, health findings and provenance.

| Renderer | Output |
|---|---|
| `write_json` | One self-describing manifest. Every field of the IR. |
| `write_csv` | `<block>_jacobian.csv` (outputs down the rows, levers across, base case in the margins) plus `base`, `bounds`, `links`, `specs`, `prices`, `duals`, `health` and `names` tables |
| `write_lp` / `write_mps` | The assembled LP in CPLEX LP or free MPS form, via Pyomo |
| `write_iterations_csv` | Radius, predicted vs realised merit, `rho`, accepted |

Some details that matter in practice:

- **Shadow prices.** `res.duals` and `res.solution` solve the final LP once and
  return the marginals by row and bound name. They ship inside the JSON and in
  `duals.csv`, because they are the first thing a planning engineer asks for.
- **Scaled coefficients.** `DeltaVector.scaled_J` is `scaled_jacobian` — the
  fractional change in an output per full-bound-range move of a lever. That is
  the dimensionless shift-vector form a planning table is normally read in, and
  `write_csv(..., scaled=True)` writes it instead of the raw matrix. Comparing
  raw `J` entries across a model in mixed units compares unit conversions.
- **Names.** LP and MPS files cannot hold `ngl.residue_F`, so every symbol is
  sanitised to `ngl_residue_F`. The mapping is recorded in `names.csv` and in
  the manifest's `meta["lp_symbols"]`, so the round trip is recoverable.
- **Sense.** Pyomo writes a minimisation; a maximising plan is emitted with
  negated costs and the constant term dropped. `meta["lp_sense"]` and
  `meta["lp_objective_offset"]` are what recover difflow's own objective.
- **Health travels with the coefficients.** `check_delta_health` findings are in
  the export, so a dead lever or a recycle loop gain near one is visible
  downstream and not just locally.
- **The export is one-way.** There is no importer for a foreign LP.

A single linearisation can be exported with no planner run at all, which is the
common case for a simulation engineer handing a table to a planning group:

```python
dvs = DeltaVectorSet.from_block(blk, radius=0.2)
write_csv(dvs, "tables/")
```

And from the shell, via the `difflow` entry point:

```bash
difflow plan-export plan.py --format csv -o tables/
difflow plan-export model.py -u reactor.V -u 'feed:F1.total_flow' \
    -y product.F_B --lb 0.5,5 --ub 5,20 --radius 0.2 -o plan.json
```

The source is a `.py` script — run, then searched for the first `PlanResult`,
`DeltaBasePlanner`, `Block` or `Flowsheet` — or a serialized flowsheet `.json`.
`--format lp|mps|iterations` needs a source that actually runs the planner,
since there is no assembled LP without one.

## What this module is not

**Not a commercial planning system.** Explicitly out of scope: crude assay
libraries, unit submodel libraries, empirical blending correlations (octane,
RVP, cloud point), scheduling, multi-site transport, reporting and audit trails.
Those are the real moat of a product like PIMS; they are not derivable from a
flowsheet, and pretending otherwise would make this module worse.

**Not a global optimiser for pooling and blending.** Bilinear terms — a stream
quality multiplied by a stream flow — stay nonconvex no matter how good the unit
linearisation is (Haverly, *ACM SIGMAP Bulletin* **25**, 1978). That needs
global optimisation, not better delta vectors. Pricing here is linear in the
block outputs, which is what keeps the problem an LP; put a bilinear pooling
term in and none of the guarantees above apply.

**Not a recycle solver between blocks.** Inter-block recycles are rejected;
merge them into one flowsheet, where difflow differentiates the tear solve for
you.

## API summary

| Object | Purpose |
|---|---|
| `Block` | A linearisable submodel: any pure JAX callable `u -> y` |
| `Network` | A DAG of blocks joined by output-to-input links |
| `DeltaBasePlanner` | The trust-region loop over the delta-base LP |
| `PlanResult` | Plan, delta vectors, LP/Pyomo model, history, sensitivity |
| `TrustRegionOptions` | Radius schedule, acceptance thresholds, iteration caps |
| `Spec` | A linear constraint, elastic by default, with optional back-off |
| `LPModel` / `LPSolution` | The assembled program and its solution |
| `Linearization` | One block's delta vectors, base point and phase regime |
| `Block.from_flowsheet` | The bridge from a difflow flowsheet to a planning block |
| `DeltaVector`, `DeltaVectorSet` | The export IR: coefficients plus everything needed to use them elsewhere |
| `write_json`, `write_csv` | The neutral interchange pair |
| `write_lp`, `write_mps` | The assembled LP in CPLEX LP / free MPS form |
| `write_iterations_csv` | The trust-region audit trail |
| `PlanResult.duals` | Shadow prices at the plan, by row and bound name |
| `difflow plan-export` | The same four formats from the shell |
| `DeltaBasePlanner.describe` | The problem statement: objective, decisions, bounds, links, specs |
| `LPModel.as_text` | The assembled program written out row by row |
| `draw_chain`, `draw_planning_network` | The flowsheet, and the network as the LP holds it |
| `draw_delta_vectors`, `draw_taylor_model`, `draw_trust_region` | The model, its locality, and the loop |
| `check_delta_vectors` | Verify AD deltas against central differences |
| `check_model_order` | Linear against quadratic, measured on the block itself |
| `model_order=` | `"linear"`, `"quadratic"` or `"auto"` subproblems |
| `QPModel`, `build_qp` | The quadratic subproblem, warm-started from the LP |
| `convexify`, `ConvexificationReport` | Spectral modification, and what it changed |
| `restoration_model` | The phase-one program for an infeasible subproblem |
| `restoration_violation` | Predicted infeasibility at a phase-one solution |
| `block_curvature`, `Curvature` | Exact Hessian by HVP, with its definiteness verdict |
| `hvp`, `hessian_of`, `block_hvp` | The second-order primitives |
| `ModelOrderReport` | Which model earns its cost here, and why not the other |
| `choose_ad_mode` | `jacrev` vs `jacfwd`, chosen by shape |
| `PhaseBoundaryWarning` | Raised when a proposal crosses a phase boundary |
| `check_delta_health` | Dead levers, recycle amplification, ill-conditioning |
| `check_lp_scaling` | Constraint-matrix scale spread (the units problem) |
| `composed_sensitivity` | End-to-end relative sensitivity, by one AD pass |
| `scaled_jacobian` | Delta vectors made dimensionless for comparison |
| `HealthReport`, `Finding` | Diagnostic results; `.summary()`, `.warn()` |
| `plan_sensitivity` | `d(plan)/d(price)`, `d(plan)/d(parameter)` |
| `price_switch_point` | The finite price at which a bang-bang lever flips |
| `Modifiers`, `run_modifier_adaptation` | Zeroth- and first-order plant corrections |
| `attribute_deltas`, `AttributionResult` | Those corrections estimated from plant history, with estimability, aliases and a structural check |
| `constraint_backoff`, `apply_backoff` | Coefficient covariance to spec margin |
| `PiecewiseSpec`, `sample_piecewise` | Batched SOS2 piecewise-linear blocks |
| `gradient_cost_ratio`, `scaling_study` | The AD-versus-perturbation measurement |
| `difflow.planning.chain` | The two-plant reference chain used above |

## References

- Baker, T. E. and Lasdon, L. S. *Successive Linear Programming at Exxon*.
  Management Science **31**(3), 264–274, 1985.
- Eason, J. P. and Biegler, L. T. *A trust region filter method for glass
  box/black box optimization*. AIChE Journal **62**(9), 3124–3136, 2016.
  [doi:10.1002/aic.15325](https://doi.org/10.1002/aic.15325)
- Marchetti, A., Chachuat, B. and Bonvin, D. *Modifier-Adaptation Methodology
  for Real-Time Optimization*. Ind. Eng. Chem. Res. **48**(13), 6022–6033, 2009.
  [doi:10.1021/ie801352x](https://doi.org/10.1021/ie801352x)
- Haverly, C. A. *Studies of the behavior of recursion for the pooling problem*.
  ACM SIGMAP Bulletin **25**, 19–28, 1978.