Delta-Base Planning

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#

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$:

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:

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#

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:

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):

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)

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:

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:

Linearise every block at the current point — the delta vectors.
Solve the LP inside a trust region to get a proposal.
Evaluate the caller’s own nonlinear blocks at that proposal.
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:

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. 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:

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 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:

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:

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:

Block(..., phase_fn=lambda u, th: jnp.atleast_1d(info["V_frac"]),
      phase_names=("V_frac",), phase_bounds=(0.0, 1.0))

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:

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.

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:

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.

[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:

>>> 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.

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{split}\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}\end{split}\]

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):

\[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:

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:

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$:

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.

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.

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:

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:

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:

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:

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:

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.

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. The details, including why yields are the levers and cut points the outputs, are in 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.

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

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:

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

And from the shell, via the difflow entry point:

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
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
Haverly, C. A. Studies of the behavior of recursion for the pooling problem. ACM SIGMAP Bulletin 25, 19–28, 1978.