Design: Aspen PIMS, and delta-base planning beyond refining

Design: Aspen PIMS, and delta-base planning beyond refining#

Status: design proposal. Nothing described here is implemented. No module difflow.planning.pims exists. This document records the shape of the feature, the decisions already taken, and the one thing that has to happen before any code is written.

It covers two plays that share a mechanism and almost nothing else. The first — connecting to PIMS, which is the bulk of this document — is an integration into a market with a thirty-year incumbent, and it deliberately keeps difflow subordinate: we supply vectors, the customer’s PIMS keeps the model. The second is the observation that the exclusions which force that subordination are refinery-specific and stop binding the moment the domain changes; it starts at Beyond refining. Between them sits the digital twin, which both plays want and neither has.

Table of contents#

Summary
Why there is something to connect
What this is not
Architecture: three tiers
The hard part: basis and naming
Provenance: what difflow can ship that PIMS cannot
API sketch
Testing without a licence
Phase 0: the gate
Decisions taken, and alternatives rejected
Beyond refining: where the exclusions stop binding
The digital twin, and where the loop breaks
Structural gaps both of those depend on
Risks
Where it lands

Summary#

This section and the nine that follow are the PIMS bridge. The domain-generality and digital-twin arguments begin at Beyond refining.

difflow supplies unit submodel delta vectors to a PIMS model the customer keeps and continues to run. It does not replace the planning LP, the recursion, the assay library or anything else PIMS owns.

The deliverable in one line: a PIMS submodel table, generated by AD from a difflow flowsheet, carrying the trust radius it was validated over and its finite-difference cross-check — and, read in the other direction, a report saying which entries of the incumbent PIMS vectors have gone stale and what the plan does when they are corrected.

Why there is something to connect#

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 \(J\) is built 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 — and why a planning model is usually running on submodels that describe the plant as it was, not as it is.

A difflow 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, and it costs a small constant multiple of one model evaluation regardless of \(n\) (difflow.planning.benchmark measures it; the ratio is locked in as a regression test).

That asymmetry — \(O(1)\) against \(O(n)\) — is the whole basis of the feature. It is not “difflow is a better planner.” It is “the expensive step in your existing planner is free here.”

What this is not#

difflow.planning already declines two categories of work, for two entirely different reasons, and the PIMS integration inherits both boundaries.

Bilinearity is a mathematical exclusion. A stream quality multiplied by a stream flow stays nonconvex however good the unit linearisation is (Haverly, ACM SIGMAP Bulletin 25, 1978). Better delta vectors do not touch it: the nonconvexity lives in the blending structure, not in the unit submodel. Every guarantee in difflow.planning chains off the subproblem being an LP — each trust-region subproblem solved to global optimality, the Eason–Biegler filter convergence theory, duals reading as prices — and a bilinear term voids all three silently while the code still returns a number. Pools and property recursion stay on the PIMS side.

Assays, blending correlations and scheduling are a scope exclusion. A crude assay library is curated proprietary data; octane, RVP and cloud-point rules are empirical fits. None of it is derivable from a flowsheet, which is difflow’s only claim to novelty here. Reimplementing it would mean competing with thirty years of accumulated domain content on its own terms, and losing.

Also out: replacing the PIMS LP, and reverse-engineering the PIMS model database. Integration stays on the documented spreadsheet interchange path.

Architecture: three tiers#

Ordered by how much works without an Aspen licence.

Tier 1 — Export#

This is a new writer on existing machinery, not a new exporter. difflow.planning.export already defines DeltaVector / DeltaVectorSet and writes them as write_json, write_csv, write_lp and write_mps. PIMS becomes a fifth writer:

write_pims(dvs: DeltaVectorSet, mapping: PIMSMapping, path) -> None

DeltaVector already carries nearly everything the table needs — u0, y0, J, lb/ub, radius, tr_lo/tr_hi, mode, phase, and per-lever u_units/y_units. It also carries scaled_J, the dimensionless Jacobian, whose own docstring calls it “the form a shift-vector table is normally read in”. That is the PIMS shift vector in all but name and file format.

So Tier 1 adds exactly two things that do not already exist: the basis conversion (below), and the tag mapping. Everything else is plumbing an existing record into a spreadsheet layout.

Block.from_flowsheet(fs, u=[...], y=[...]) is the front of the same pipeline — the bridge from a difflow flowsheet to a planning block, with lever keys in _apply_params notation and feed levers under a feed: prefix. A PIMS export is then: flowsheet → Block.from_flowsheet → linearize_block → DeltaVectorSet.from_block → write_pims.

Output is a spreadsheet, one sheet per table. Buildable and testable today with no Aspen software present.

Tier 2 — Import, and the staleness report#

read_pims_model(path) -> list[Block] parses submodel tables back into affine Blocks. Note what this changes: difflow.planning.export is one-way today — there is no importer at all, for any format. A PIMS reader would be the first, so it is a change to the module’s contract and not merely another file parser. It earns that only because of what it unlocks:

a staleness report. Incumbent PIMS vector against fresh AD vector, entry by entry, ranked by economic impact rather than by relative error — “six of these forty-two entries move the plan; here they are, and here is what the plan does when they are corrected.” A large relative error on an entry the LP never touches is not news; a small one on an active lever is.

Secondary, nearly free once the importer exists: run difflow’s guarded trust-region loop and plan_sensitivity over the PIMS structure. PIMS cannot offer that itself, because its submodels are not differentiable.

Tier 3 — Live coupling#

Writing into a running model, or driving whatever automation surface the installed PIMS version exposes. Deferred: it needs a licence and an install to scope responsibly. Follow the existing precedent for third-party bridges — difflow.dwsim_import is pythonnet-gated (checked against DWSIM 9.0.5), difflow.pyglenn_import and difflow.cantera_import are optional imports — and keep it behind a lazy import so the core path stays dependency-free.

The hard part: basis and naming#

Writing tables is a week of work. This section is the feature.

Basis conversion. PIMS submodel yields are per unit of feed on a weight or volume basis and are expected to close. difflow blocks emit absolute molar flows. So J must be divided through by the feed rate and pushed onto the target basis via molecular weights and densities — and the exporter must assert closure, not hope for it. A yield set that silently fails to sum to one produces a planning model that is wrong in a way no reviewer will catch by reading it.

Property rows are not free outputs. Rows in a PIMS submodel that carry recursed structural qualities are not the same kind of object as a difflow block output. Only the subset of y_names that maps onto a genuine yield or property row can be exported. Anything else — a spec variable such as T_colfeed, which limits but is not priced — has no PIMS home and must raise, never be quietly dropped.

Naming. PIMS row and column tags are short and structured. The mapping is explicit, declared per block, and validated on export. No automatic mangling of difflow names into tags.

Shift ranges. PIMS shift levers have modeller-set ranges and are linear over them by construction; difflow’s trust radius is adaptive per cycle. Export commits to a range, and records which one.

Provenance: what difflow can ship that PIMS cannot#

Every exported vector carries:

the trust radius it was validated over;
the check_delta_vectors residual against central differences at that radius — the same \(2n\)-evaluation method the vector replaces;
the linearisation point u0, and the phase regime recorded at it.

Three of those four are already fields on DeltaVector (radius, tr_lo/ tr_hi, u0, phase); only the residual would be added. check_delta_vectors already exists, so this is cheap. It matters because a delta vector’s accuracy is local and the receiving LP has no way to know where it stops. A vector that arrives stating its own domain of validity is a different object from one that arrives as a block of numbers.

A block that declares a phase_fn raises PhaseBoundaryWarning when the linearisation straddles a phase appearance or disappearance. Such a vector is meaningless rather than merely inaccurate — the underlying function is not differentiable there — and must never be exported silently.

API sketch#

Illustrative, not settled; the table layout in Phase 0 will move it.

from difflow.planning import Block, linearize_block
from difflow.planning.export import DeltaVectorSet
from difflow.planning.pims import PIMSMapping, write_pims, read_pims_model

blk = Block.from_flowsheet(fs, u=["deethanizer.reflux_ratio", "feed:feed.T"],
                           y=["ngl.F_C2", "residue.F_total"])
dvs = DeltaVectorSet.from_block(blk, linearize_block(blk))

mapping = PIMSMapping(
    submodel="SDEETH",              # the PIMS submodel this block feeds
    feed="feed_F",                  # which lever is the feed yields are per
    basis="weight",                 # "weight" | "volume"
    tags={"ngl.F_C2": "...", "residue.F_total": "..."},   # difflow -> PIMS tag
    mw={...}, density={...},        # what the conversion needs
)

write_pims(dvs, mapping, "deethanizer.xlsx")   # asserts closure, or raises

incumbent = read_pims_model("plant_model.xlsx")          # -> list[Block]
report = compare_vectors(incumbent["SDEETH"], dvs, prices=prices)
print(report.summary())        # entries ranked by effect on the plan

Testing without a licence#

No PIMS in CI, so the test strategy carries the weight:

Round trip against a synthetic fixture workbook in the PIMS table layout, committed to tests/fixtures/.
Closure assertions on every basis conversion, including deliberately non-closing inputs that must raise.
The golden test: the exported base-plus-shift model reproduces the block’s nonlinear response inside the stated radius, to the stated tolerance. That is the property being sold, so it is the property under test.
Refusal tests: an unmapped output raises; a phase-straddling linearisation raises.
Any live-PIMS test marked optional and licence-gated.

Phase 0: the gate#

Obtain one real PIMS submodel table, even redacted, before writing code.

Tag conventions, the layout of base and shift columns, and which rows are yields versus recursed properties are facts that cannot be derived from this side. Guessing produces an exporter that looks right and is wrong — the failure mode this whole design is organised against.

Pick the pilot unit at the same time: one whose difflow block maps naturally onto a PIMS submodel (a fractionator, a simple conversion unit), not one whose outputs are mostly recursed qualities.

Decisions taken, and alternatives rejected#

Integrate at the submodel, not the LP. difflow supplies vectors; PIMS keeps the model. Rejected: emitting a whole planning model, which would require the pooling and assay machinery that is out of scope by design.

Spreadsheet interchange, not the model database. Rejected: parsing the proprietary model store — fragile, and legally murkier.

A module, not a plugin. src/difflow/planning/pims.py, alongside lp.py. difflow.plugins is for unit operations; difflow.planning is deliberately not one, and the PIMS bridge sits at the same level.

The LP solver question is orthogonal, and was investigated separately. The planning LP stays on HiGHS via scipy.optimize.linprog / milp. What difflow exports is J from linearize.py; what validates it is the nonlinear block. Neither touches the LP solver, so no backend decision can block this work.

Two incidental findings from that investigation are worth recording here, because they are not obvious from reading the code and they constrain any future solver swap:

LPSolution.duals and LPSolution.active_bounds have no consumers anywhere in the repository. plan_sensitivity derives its multipliers by AD on the caller’s nonlinear blocks (sensitivity.py:304), not from LP duals — deliberately, and consistent with the invariant that violations are scored from the nonlinear model and never from LP slacks. The LP is a throwaway subproblem; its duals are shadow prices of a Taylor model, not of the plant.
If an interior-point LP backend is ever adopted, crossover must be enabled. The load-bearing consumer of vertex structure is not any bound-reporting helper but _hit_trust_region (planner.py:709), which gates radius expansion at planner.py:656. An IPM returning the analytic centre of a tied optimal face never reports at_boundary, so the radius can only shrink and the loop crawls to radius_min reporting convergence — a silent convergence-rate regression. A snap-to-bound tolerance does not address this: on a tie the analytic centre is \(O(1)\) from the bound, not \(10^{-9}\).

Beyond refining: where the exclusions stop binding#

What this is not declines two categories of work for two different reasons, and the difference between those reasons is the whole of this section. Bilinearity is a mathematical exclusion; assays, blending correlations and scheduling are a scope exclusion. Only the first is a statement about delta-base planning. The second is a statement about refining — and it is the one that makes difflow subordinate here.

Nothing in the mechanism is refinery-specific. A trust-region SLP over AD delta vectors is generic; Baker and Lasdon happened to write it up at Exxon. What is refinery-specific is who holds the scarce asset. In a refinery it is curated crude assays and empirical property correlations, accumulated over decades, and difflow does not have them and should not try to. Change domain and the scarce asset changes with it: where a planning model’s submodels are physics, the thing nobody has is a rigorous model you can differentiate, and that is the only asset difflow holds.

difflow already carries plugins for five domains of that kind — difflow_gas, difflow_power, difflow_cc, difflow_bio, difflow_ree — which is not the same as five planning problems. Candidates, and the reason each is or is not one:

Gas transmission — the strongest candidate. Multi-period nomination planning over difflow_gas. The submodels are Weymouth pipes and compressor stations, network_residuals is a single traceable definition of the equation set, and the tear solve is differentiated implicitly, so the reduced Jacobian is already there. Line-pack is genuine inventory: a state that must carry from one period to the next, which is exactly the structure a planning model has and the reference chain does not.
REE separation — the clean scope-exclusion inversion. Campaign planning across ore lots whose feed composition changes lot to lot. There is no assay library to lose to, because the assay library would be difflow_ree’s database, and n_stages is already a continuous traceable decision.
Carbon capture fleet retrofit. Planning against a carbon price, where the question being asked is a price sensitivity — which is what plan_sensitivity and price_switch_point return and no commercial planning system can, because its submodels are not differentiable.
Power dispatch — a non-candidate, named so the mistake is not made. difflow_power solves AC-OPF directly through its own interior-point NLP. A delta-base layer earns nothing where the full nonlinear program is already tractable; linearising it would be strictly worse than the solver that is there. Delta-base planning earns its keep when the rigorous model is too expensive or too structured to put inside the optimiser whole, not when it fits.

Two cautions, both load-bearing.

Bilinearity does not disappear by leaving refining. Any domain where a quality rides on a flow into a pool carries the Haverly nonconvexity. Gas heating value and Wobbe index blending have it outright; a raffinate blend has it. Every guarantee in difflow.planning chains off the subproblem being an LP, and a bilinear term voids all of them while the code still returns a number. The mathematical exclusion survives the domain change intact — so checking for pooling structure is part of choosing a domain, not something to discover later.

No incumbent cuts both ways. Outside refining there is no PIMS to integrate with, so no Phase 0 gate blocks the work — and equally no customer with a planning model already running, a modelling group that understands shift vectors, and a budget line for keeping them fresh. This is a harder sell and an easier build. The PIMS bridge is the reverse.

The digital twin, and where the loop breaks#

The ambition is that a flowsheet is not a one-off source of linearisations but a model kept current against operating data, so that a plan is built on the plant as it is rather than as it was commissioned. That is the same complaint the opening argument makes about annually refreshed delta vectors, pushed one step further: refresh the model, and the vectors follow.

Every piece of it exists.

Step	Where
Is today’s data consistent with the model?	`difflow.reconciliation.monitor`
Bad sensor, or model drift?	`MonitorResult.diagnose` — blame concentration, `reconciliation/monitoring.py`
Re-estimate the drifted parameter over a window	`reconcile_multi`, holding a parameter common across data sets
Track it online through the dynamics	`difflow.mhe`; `MHEResult.parameters` already returns the shape `Block.theta` takes
Refresh the delta vectors	nothing — `jax.jacobian` follows the parameter change for free
Size the constraint margin for what the estimate does not know	`constraint_backoff`, from the estimate’s Σ_θ

That fifth row is the point. Refreshing the vectors after a parameter update is the step that costs O(n) simulator perturbations in a commercial system, and here it is not a step at all.

Where it breaks. run_modifier_adaptation(planner, plant_fns) (planning/modifiers.py:178) requires plant callables, and update_modifiers obtains the gradient correction λ by running linearize_block on the plant function (modifiers.py:122). A historian supplies neither: it supplies noisy records at whichever operating points the plant happened to visit, with no gradient anywhere. The docstring already declares the seam — “in a real deployment they would come from a gradient estimator, and this function accepts whatever plant_fn provides” — but no gradient estimator exists, and the driver loop hard-requires callables. As shipped, modifier adaptation is an RTO study against a simulated plant. It is not yet a data-driven twin.

The way through is to split the mismatch, because the two halves want entirely different machinery and only one of them is hard.

Parametric drift needs no modifiers at all. Re-estimate θ and let AD carry the change onto J. This moves the model rather than papering over it, it is the case difflow is already fully equipped for, and — by diagnose()’s own account — it is the case that can actually be detected from routine monitoring. Shipping this loop is assembly, not research.

Structural mismatch is the one that genuinely needs λ from data, and a historian is close to the worst possible source for it: steady-state records cluster at the operating points the plant actually holds, which is precisely where the gradient is least identifiable. Recovering λ needs either deliberate excitation — a dither or a designed step, and difflow.estimation.design already does experiment design — or a gradient estimated over a window of recent operating points, carrying its own covariance. The second is a research question, not plumbing, and it should be gated accordingly.

Two things to record before either is built:

Whatever supplies λ must supply its uncertainty. update_modifiers filters at a fixed gain of 0.5 precisely because an unfiltered gradient correction from noisy data oscillates. A constant is a placeholder; from a real estimator the gain should follow the estimate’s covariance, the same way constraint_backoff sizes κσ from Σ_θ.
A twin whose parameters move invalidates every vector already exported. That is an argument for Tier 2: the staleness report described under Architecture is the only artefact here that says which entries actually moved the plan, and a live twin is what makes that question recur rather than arise once.

Structural gaps both of those depend on#

Neither gap below touches the PIMS bridge — Tier 1 exports one block’s vectors and PIMS owns the commercial structure. They bite only on the play where difflow is itself the planner.

There are no commercial variables. build_lp creates columns for block inputs, block outputs, elastic spec slacks and piecewise SOS2 weights (planning/assemble.py:120) and nothing else. A planning model needs purchases at a tiered price, sales against a contract, transport arcs, inventory. Some of that can be faked with an identity-function Block, but such a block carries a Jacobian, a trust region and a phase regime it has no use for, and its lb/ub become the only place a contract limit can live. That is a design decision worth taking deliberately rather than discovering.

Periods are replicated, not linked. two_plant_chain(horizon=4) names blocks ngl@t0 … ngl@t3 and couples them only through a shared CO₂ cap (planning/chain.py:341). That is a horizon built to make the AD scaling argument measurable, which is what it was for — not a multi-period planning model, which couples periods through inventory.

But a Link is output-to-input, and Network._topological_order (planning/network.py:118) rejects only cycles — a self-recycle, or a loop among blocks. A forward link from tank@t0.level to tank@t1.level_in is an ordinary DAG edge and is legal today. Multi-period inventory may therefore need no new machinery at all, only a demonstration — and finding that out is the cheapest item in this document. It should be settled before anything else here is designed, because the answer changes the size of both plays.

Risks#

Risk	Severity	Mitigation
Real submodels don’t map cleanly onto difflow blocks (feed basis, recursed property rows)	Could kill the feature	Phase 0; pick the pilot unit for mapping cleanliness
Silent basis-conversion error	High — wrong plan, no visible symptom	Closure assertions; golden test against the nonlinear model
Table layout guessed wrong	Medium — rework	Phase 0 gate
Vector exported outside its validity domain	Medium	Radius and FD residual travel with every vector; phase warnings raise
Customer expects a full PIMS replacement	Medium — expectation, not code	Scope section, stated up front

Risks carried by the two other plays rather than by the bridge:

Risk	Severity	Mitigation
Bilinearity assumed to be refinery-specific, and is not	High — voids every LP guarantee while still returning a number	Pooling structure is out of scope in every domain; check for quality-times-flow when choosing one
Outside refining there is no incumbent to integrate with, and no customer already running a planning model	Medium — a harder sell, not a harder build	Choose a domain where difflow already owns the physics plugin and a real planning question is being asked
Structural mismatch needs plant gradients a historian cannot supply	Medium — confines the twin to parametric drift	Ship the parametric loop, which is assembly; gate the structural one on a gradient estimator that reports its own covariance
A delta vector exported before a twin updates θ is silently stale	Medium	Tier 2 staleness report; a live twin makes this recur rather than arise once

Where it lands#

Path	Contents
`src/difflow/planning/pims.py`	`PIMSMapping`, `write_pims`, `read_pims_model`, `compare_vectors` — a writer on `export.DeltaVectorSet`, kept out of `export.py` because the importer and the basis layer have no business there
`tests/test_planning_pims.py`	Round trip, closure, golden, refusal tests
`tests/fixtures/`	Synthetic PIMS-layout workbook
`docs/planning.md`	Cross-reference once implemented
`pyproject.toml`	New `pims` extra (`openpyxl`) — there is no Excel dependency in the project today

That table is the PIMS bridge only. The other two plays land in difflow.planning proper — new LP columns in assemble.py, a period-linking demonstration on the existing Link machinery, a data-driven plant_fn seam in modifiers.py — with docs/planning.md and a worked example notebook rather than this document. None of it is scoped until one of them is chosen, and the multi-period question should be answered first because it changes the size of the rest.

Related reading: docs/planning.md, difflow.planning.export (the writers this extends), Block.from_flowsheet, examples/30_delta_base_planning.ipynb, difflow.planning.chain.two_plant_chain, docs/data-reconciliation.md and docs/moving-horizon-estimation.md (the twin’s existing halves), examples/29_model_updating.ipynb.