Introduction

What the Tennessee Eastman Process is

The Tennessee Eastman Process is the chemical engineering profession's standard open benchmark for process control and for fault detection. Downs and Vogel published it in 1993 as a challenge problem: a model of an industrial chemical plant whose components are identified only by the letters A to H, distributed as Fortran so that competing control schemes and competing detectors could be compared on identical dynamics rather than on identical prose descriptions. Thirty years of published results rest on that code.

The plant makes two liquid products, G and H, from four gaseous feeds, A, C, D and E, with an inert B and a liquid byproduct F. Four irreversible, exothermic gas-phase reactions run in the reactor's vapour space (teprob.f:503-528):

1:  A + C + D -> G          2:  A + C + E -> H
3:  A + E     -> F          4:  3 D       -> 2 F

B takes part in none of them. It arrives with the mixed feed and leaves only through the purge, which is the entire reason the plant has a purge. Around those reactions sit five unit operations, described one page each later in this book: a feed mixing zone, the reactor with its cooling coil and agitator, a condenser and vapour-liquid separator, a compressor with a recycle loop, and a steam stripper. The whole plant is fifty integrated states (teprob.f:24-26), forty-one measurements, twelve manipulated variables, and twenty programmed disturbances (teprob.f:340, which loops DO 500 I=1,20).

Those five units, the four feeds and the thirteen streams between them are the whole plant, and it is worth seeing them laid out once before reading anything else. The diagram below was drawn from teprob.f rather than from the paper's figure, so every stream carries both of the numbers it goes by: the FTM index the source uses, and, where one exists, the stream number Downs and Vogel print. Those two disagree on three of the four feeds, which is the most common way to wire up a reimplementation incorrectly.

Two properties make it a good benchmark and a hard one. It is open-loop unstable, so a run with the valves held still does not merely drift, it trips: this port measures that trip at 3.060 simulated hours on reactor pressure (LOG entry B-0041). And its composition analysers report their previous sample rather than the current one (teprob.f:711-761), so the measurements a controller sees carry real dead time.

What this port is

tepsim is a pure-Rust reimplementation of that Fortran, ported from the original source rather than from any later reimplementation. The reference material is the Braatz group release, vendored unmodified under reference/: teprob.f, 1594 lines holding the plant model and eight utility routines, and temain_mod.f, 1413 lines holding the Euler driver and the decentralised PI control suite.

Two design commitments shape everything. The first is that the port is checked against the Fortran continuously rather than at the end: a development-only crate, tepsim-oracle, compiles the unmodified Fortran with gfortran and links it into the test binary, so a Rust test can force both implementations into the same state and compare them in the same process. The second is determinism. The core crate is no_std, forbids unsafe, uses no f32, no SIMD, no clock and no randomness outside the model's own generator, and answers exp, pow and ln from a vendored pure-Rust libm, so that the same scenario produces the same bits on x86-64, on aarch64 and in a browser.

Every place the port knowingly departs from the original is written down, with its class, its measured effect and the test that measured it, in the quirk and delta register.

The headline result

A complete 48-hour closed-loop run, 172,800 integrator steps at a one-second step, is bit-identical to the Fortran in all 41 measurements and all 12 manipulated variables, at every one of those steps, when both are given the same exp and pow.

That measurement is from the LOG.org entry for B-0041, "Phase 4 acceptance: the driver's full 48-hour run", with gfortran 15.2.0:

The two rows are the two halves of one claim. Under libm-system, where the transcendental functions are the ones gfortran itself calls, the port and the original agree exactly: not to a tolerance, to the bit. Under the vendored libm that actually ships, the only difference is transcendental rounding, and after two simulated days it has not reached a tenth of a billionth of any instrument's noise standard deviation.

The size of that transcendental difference is measured too. The first transcendental call in the model is the Antoine vapour pressure at teprob.f:485 and teprob.f:488; over the whole range of arguments this model reaches, the vendored libm and gfortran's disagree on 9.945% of them, by exactly one ULP each time (LOG entry B-0018). Everything downstream of a condensible partial pressure therefore carries about 1.1e-16 of relative difference that no care in the algebra removes. Rather than accept a tolerance and lose the sharper claim, every differential test runs twice, once against each libm, and the second run is held to zero ULP.

Since B-0052 that claim attaches to the public API and not only to the test harness: tepsim::Simulation reproduces the validated loop bit for bit over the nominal scenario and IDV(1), IDV(6), IDV(13) and IDV(20), on all 60 samples and all 53 channels, in both libm configurations.

Run it in the browser first

TEP Studio is this simulator compiled to WebAssembly, running entirely in your browser with no server and nothing to install. Start the plant, switch a disturbance on, and watch the controllers fight it. Every scenario is shareable as a URL, because the whole scenario is in the link.

It is the same code the rest of this book documents, not a demonstration model: the page prints a determinism digest that matches the one the native build computes, and apps/studio/node/deployed.test.mjs asserts that on the module that actually ships.

How to read this book

Getting started is the working Rust API and the command line tool, and is the fastest way to a CSV of plant data. The Python package is the same ground for the tepsim wheel, with installation, a quickstart and the API reference, and the four tutorials are the narrative around the executed notebooks in notebooks/. The plant and the right-hand side explain how the fifty states, the thirteen internal streams and the six hundred lines of TEFUNC are organised here. The five unit-operation pages carry the governing equations, the variable tables and the teprob.f line ranges they were derived from. Validation is the ten-tier ladder and the numbers actually measured at each rung. Status says what is finished and what is not, which for a port in progress is the page to read second.

Every claim about the model in this book cites the teprob.f or temain_mod.f line range it came from, so a sceptical reader can check it against the vendored source. Every validation number cites the LOG.org iteration that measured it.

Status

Everything planned is built. Every phase has landed, every validation tier has a harness, and the four decisions that were once open have been taken. The backlog is 85 items done, none open and none blocked.

That makes this page shorter than it used to be, and its job is now to say where the limits are rather than where the edge of construction is. Numbers here come from the current-state block of BACKLOG.org and the closing entries of LOG.org, and each names the iteration that measured it.

What exists

The plant model, the control layer, the public API, the command line, the Python wheel and the browser app are all complete and validated. tepsim-core ports the whole of TEFUNC (teprob.f:196-816); tepsim-control ports the twenty control loops and the driver's own scheduling; tepsim is the API most callers want; tep drives it from a terminal, including tep dataset for generating d00-d21 shaped files; the wheel carries no C dependency; and TEP Studio runs the whole simulator in a browser tab with no server.

The differential harness is tepsim-oracle, which compiles and links the unmodified Fortran, and tepsim-stats, which implements every statistic the ladder needs in Rust with known-answer tests, so no part of the validation depends on numpy or scipy.

All ten tiers have a harness and all ten have run. Tier 8, differential fuzzing with shrinking, and Tier 9, cross-platform determinism, were the last two and landed together. The validation chapters are generated from the suite's own output for the tiers that have a generator; the narrative in Validation is hand-written and transcribed from LOG.org, and each number there names the iteration that measured it so the transcription can be checked.

What is genuinely still missing

These are limits, not unfinished work, and most of them are limits of the machine this was built on rather than of the code.

Tier 9 has no x86-64 leg and no real-browser leg. Six committed digests are identical on aarch64 and on wasm32 under Node, across three build profiles, and the browser app's transport path reproduces them independently. Nobody has run them on x86-64, Windows or aarch64 Linux, and Node is not a browser. The table is committed constants rather than a value computed twice in one process, precisely so that running cargo xtask tier9 on another machine completes the claim with no code change.

Tier 8 has one open counterexample. In five million tuples, fuzz#863105 misses the 1e-12 gate at 4.607e-12 of the scale of its terms. It is recorded and attributed rather than fixed, because there is nothing to fix: it is an accepted one-ULP exp difference amplified by 973 simulated hours of IDV(13)'s kinetic drift, and it is bit-identical under libm-system.

Two IDV faults miss the Tier 5 margin under the vendored libm. IDV(14) and IDV(19), on exactly the valves those faults stick, and bit-identical under the platform libm. The cause is teprob.f:801, a discontinuous branch on a floating-point comparison. Those valves are judged on their distribution rather than on a mean, because a series of plateaux has no meaningful centre. No margin was widened.

Apache Arrow and Parquet sinks are not implemented. The recorder sinks that exist are Columnar, Csv, Ring, Decimating and Selecting, and Simulation::run_into streams into one rather than collecting. Arrow's dependency cost was judged not to be worth paying until someone has a dataset large enough that CSV is the bottleneck; the right home would be tepsim-cli, never tepsim, which is no_std and compiles to wasm32 under a size budget.

A historian export does not exist. Everything is wide tabular. A long-format tag,timestamp,value,quality record would need a caller-supplied epoch, since the core may not read a clock, a decision about what quality code a dead-time analyser and a frozen channel deserve, and units promoted from prose in measurements.rs to data on a channel.

The step size is not adapted. Three integrators exist as of B-0053: fixed step explicit Euler, classical RK4, and Dormand-Prince 5(4) with an embedded error estimate. Only Euler reproduces the original, and it is the default. A variable step changes when the discrete phases run, which is a decision about fidelity rather than about numerics, so it has not been taken.

That comparison produced a result about the original rather than about the port. RK4 and Dormand-Prince agree with each other to 1.5e-6 while both differ from Euler by about 1.1e-2, so the published Tennessee Eastman data carries roughly one percent of integration error against an accurate solution of the same equations. Reproducing that is the point, but it means "the TEP" names a particular discretisation and not only a set of differential equations.

The two Class C quirks, and what a default Scenario now does

Both were signed off on 2026-08-28, so a default Scenario fixes them and Scenario::faithful() reproduces them.

A trip ends the run (delta D-007, teprob.f:807-811). The original freezes the plant and keeps reporting, which is where four of the forty-four published files get their frozen tails: 1,832 rows in total, and 363 of d06.dat's 480. The fix is pure truncation, since every sample before the trip is bit-identical either way, and the argument for keeping the freeze was that it preserved an option it does not preserve, because the plant cannot be restarted in either case.

The driver does not force IDV(12) at hour eight (delta D-011, temain_mod.f:366-368). Tier 7 established that the published files were generated with that line replaced rather than kept: every dNN_te except d12_te sits at the nominal operating point straight across row 160.

Every comparison against the Fortran or against published data runs Scenario::faithful(), and tier5::run_port pins it so no differential can lose it. From the command line the flags are tep run --faithful, or --force-idv12 and --freeze-on-trip individually.

Getting started

This page is the Rust API and the command line tool. If you want to drive the simulator from Python, The Python package is the equivalent page for the tepsim wheel, and the tutorials and the notebooks they link to are all Python.

Building

The workspace pins its toolchain in rust-toolchain.toml, so a rustup install picks the right compiler on its own. Nothing in the shipping crates needs Fortran; gfortran is required only to build tepsim-oracle, the development-only differential harness, which sits behind the oracle feature and is never a dependency of tepsim.

$ cargo build --release
$ cargo xtask ci          # fmt, clippy, tests, docs, cargo-deny

The library

The whole public surface is three types. A Scenario describes a run, a Simulation performs one, and a Run holds the output as columns.

use tepsim::{Scenario, Simulation};

// Twenty-four hours with disturbance IDV(4), the reactor cooling water
// inlet temperature step.
let run = Simulation::new(Scenario::fault(4).with_hours(24.0)).run();

// One measurement across the run, one-based as XMEAS(n) is.
let reactor_pressure = run.measurement(7);

// One manipulated variable, one-based as XMV(n) is.
let coolant_valve = run.manipulated(10);

assert!(run.outcome.is_completed());

Scenario is a plain Copy struct, so a caller can build one, tweak it, and keep both. Its builders are with_seed, with_hours, sampling_every, with_fault and open_loop, and the two constructors are Scenario::baseline() and Scenario::fault(n).

Three of those defaults deserve a sentence. Open loop is not a useful operating mode, it is a diagnostic one: the plant trips on reactor pressure after about three hours, and the difference between the two settings is the clearest single statement of what the control layer does.

The other two are the Class C quirks, signed off on 2026-08-28 and fixed by default: a trip ends the run (delta D-007) and the driver does not force IDV(12) (delta D-011). Both are one call away. Scenario::faithful() gives the original's behaviour on both, and it is what every comparison against the Fortran or against published data runs. See the delta register for the evidence behind each.

A Run holds samples, each carrying step, hours, measurements (XMEAS(1..41)), manipulated (XMV(1..12)) and labels. Run::column, Run::columns, Run::measurement and Run::manipulated reshape it into the columns a statistic or a detector wants, and Sample::row gives all 53 channels in one array, measurements first, which is the layout every downstream consumer uses. channel_names() returns matching names in the same order.

Run::outcome is how a run ended: Completed, Tripped with the step, the hour and the first shutdown condition that fired, or SolveFailed with the step at which a temperature solve failed to converge. A trip ends the run by default, so a tripped run is shorter than its scenario planned. teprob.f:807-811 instead freezes the plant and keeps reporting, which is what Scenario::faithful() gives and what the constant tails in d06 and d18 are.

Labels is ground truth: which disturbances were active at that instant, and how long each had been. The original records nothing of the sort, so every detection-delay figure in the literature is computed against an onset the author assumed.

Stepping a run by hand

Simulation::run() is the whole scenario in one call. For an online detector, a live display or a reinforcement learning loop, Simulation::step() advances one integrator step and hands back a Sample on the steps where one is due.

use tepsim::{Scenario, Simulation};

let mut sim = Simulation::new(Scenario::baseline().with_hours(2.0));
while let Some(sample) = sim.step() {
    println!("{:.3} h  pressure {:.1}", sample.hours, sample.measurements[6]);
}

The order inside a step is temain_mod.f's: force IDV(12) if it is due, run the controllers on the previous step's measurements, integrate, then clamp. That one plant step of dead time in every loop is delta D-010, and getting it backwards leaves XMEAS(14) 23% out after four hours.

The command line

$ tep run --fault 4 --hours 24 --labels

tep has four subcommands, run, dataset, faults and help, and run takes these flags:

tep faults prints the twenty disturbances with the shape of each and the teprob.f line it acts on; the same table appears in The twenty disturbances.

Output is CSV on stdout, with progress and the outcome on stderr, so that redirecting stdout to a file gives a clean file. The header is step,hours followed by the 53 channel names, plus fault and hours_since_onset when --labels is given. Values are written with seventeen significant digits, which round-trips an f64 exactly, so a CSV written here is reproducible rather than approximately reproducible.

A trip is reported on stderr and exits successfully, because a trip is a result rather than a malfunction; the message says whether the run ended there or the plant froze and carried on. Only a failed temperature solve is an error exit.

Generating a dataset

$ tep dataset --out ./data
$ tep dataset --set all --faults 0,4,6 --format csv --out ./small
$ tep dataset --list

tep dataset writes files with the same geometry as the published d00 through d21: 960 rows of 52 columns for the _te half, 480 (500 for d00) for the training half, at the seeds teprob.f:1187-1256 records, in the same fixed sixteen-column layout the shipped files use. --set chooses te (the default), training or all; --faults takes all or a comma list; --format is dat or csv, the latter at full precision with a header.

It is emphatically not the published data, and it says so on stderr every run. Four things stand in the way and none is a matter of trying harder: the toolchain that made the shipped files is unrecorded and its exp was not this one's, the output was rounded to five significant figures before anyone saw it, IDV(21) is not in this revision of the model so d21 cannot be generated at all, and the protocol behind the twenty-two training files is written down nowhere. Tier 7 measures the gap rather than closing it.

Reproducibility

A run is a pure function of its Scenario. There is no clock, no thread-local state and no global, so the same scenario gives the same bits on x86-64, aarch64 and wasm32. That is what makes a recorded dataset reproducible from its description rather than from a file, and it is the property Tier 9 will assert across the CI matrix.

The Python package

The simulator ships as a Python package called tepsim. It is not a reimplementation and not a wrapper around a subprocess: it is the same Rust core the rest of this book documents, compiled into an extension module with PyO3. There is no C in the build, no Fortran, and no runtime dependency except NumPy.

A run started from Python is bit-identical to the same run started from Rust or from the browser, because all three call the same code with the same scenario. That is what makes the worked examples in the notebooks usable as evidence rather than as illustrations.

Installing it

It is not on PyPI. The version is still 0.0.0 and no release tag has been pushed, so nothing has ever been published. The publish job in .github/workflows/wheels.yml is wired up and gated on a v* tag, and the name is free, so pip install tepsim will be the line once there is a release. It is not the line today.

Install it from the repository instead:

$ pip install "git+https://github.com/jkitchin/tep-rust#subdirectory=crates/tepsim-py"

Two things about that command, both of which will bite otherwise.

It builds from source, so you need a Rust toolchain. There is no published wheel to fall back on, so pip clones the repository, invokes maturin, and compiles the whole core in release mode. That takes a couple of minutes on a warm machine and it fails partway through if cargo is not on your PATH. Install Rust from rustup.rs first. NumPy is pulled in automatically as a declared dependency, so nothing else is needed.

The quotes are not optional. Most shells treat # as the start of a comment, and without the quotes the #subdirectory= fragment is silently dropped, at which point pip tries to build the workspace root and fails with an error that does not mention the fragment at all.

The result is a cp39-abi3 wheel, so one build covers every GIL-enabled CPython from 3.9 upwards. Free-threaded interpreters have no stable ABI to target and need a version-specific build; 3.14t is the earliest that works, because PyO3 0.29 does not support the free-threaded build of anything below 3.14.

From a checkout

If you have the repository already, cargo xtask python is the supported path. It builds a release wheel, creates a throwaway virtualenv at .xtask-python/venv, installs the wheel and NumPy into it, proves that import tepsim resolves inside that virtualenv rather than somewhere else on the machine, and runs the binding's pytest suite against it.

$ cargo xtask python
$ .xtask-python/venv/bin/python -c "import tepsim; print(tepsim.__version__)"

That interpreter is the one the book's own tests use to check the quickstart transcript below. Note that cargo xtask python deletes and rebuilds the virtualenv every time it runs, so anything else installed into it, jupyter and matplotlib for the notebooks in particular, has to be reinstalled afterwards.

To build a wheel without installing it, maturin build --release -m crates/tepsim-py/Cargo.toml writes one and prints the path.

Quickstart

from concurrent.futures import ThreadPoolExecutor

import numpy as np

import tepsim as tep

print("tepsim %s: XMEAS(1..%d), XMV(1..%d), %d channels, IDV(1..%d)"
      % (tep.__version__, tep.MEASUREMENTS, tep.MANIPULATED, tep.CHANNELS,
         tep.DISTURBANCES))

# A Scenario says what to simulate, a Simulation does it, a Run holds what
# came out. That is the whole API.
run = tep.Simulation(tep.Scenario.baseline(seed=42, hours=48)).run()
print()
print("run:      %r" % run)
print("matrix:   %s %s" % (run.to_numpy().shape, run.to_numpy().dtype))
print("outcome:  %s" % run.outcome)
print("XMEAS(7): mean %.2f kPa over %.0f h"
      % (run.measurement(7).mean(), run.hours[-1]))

# The twenty disturbances say what they do, not only what the original header
# called them. Five of the published descriptions are the word "Unknown".
print()
print("the five the original leaves unexplained")
for fault in tep.faults():
    if fault.published == "Unknown":
        print("  IDV(%2d) %-9s %s" % (fault.index, fault.shape, fault.effect))

# Ground truth travels with the data, which the original records nowhere.
faulted = tep.Simulation(tep.Scenario.fault(1, hours=8)).run()
labels = faulted.labels()
print()
print("IDV(1) over 8 h: active at the last sample %s, %.2f h since onset"
      % (labels["active"][-1, 0], labels["since_onset"][-1, 0]))

# run() releases the GIL, so an ensemble is a thread pool and nothing has to be
# pickled.
sims = [tep.Simulation(tep.Scenario.fault(n, hours=8)) for n in range(1, 21)]
with ThreadPoolExecutor() as pool:
    runs = list(pool.map(tep.Simulation.run, sims))
print()
print("twenty 8-hour faulted runs: %d completed, %d tripped"
      % (sum(r.outcome == "completed" for r in runs),
         sum(r.outcome == "tripped" for r in runs)))

# A run is a pure function of its scenario, and a scenario is one line of text
# that parses back to an equal scenario.
scenario = tep.Scenario.fault(4, hours=8)
print()
print("digest:   %s" % scenario.digest)
print("text:     %s" % scenario.to_text())
print("parses back equal:      %s"
      % (tep.Scenario.from_text(scenario.to_text()) == scenario))
print("two runs bit-identical: %s"
      % np.array_equal(tep.Simulation(scenario).run().to_numpy(),
                       tep.Simulation(scenario).run().to_numpy()))

tepsim 0.0.0: XMEAS(1..41), XMV(1..12), 53 channels, IDV(1..20)

run:      <tepsim.Run 960 samples x 53 channels, 48.0 h, completed>
matrix:   (960, 53) float64
outcome:  completed
XMEAS(7): mean 2706.16 kPa over 48 h

the five the original leaves unexplained
  IDV(16) random    enables walk channel 9, the stripper steam valve capacity
  IDV(17) random    enables spike channel 10, the reactor coolant duty
  IDV(18) random    enables spike channel 11, the condenser coolant duty
  IDV(19) sticking  sticks valves 5, 7, 8 and 9; touches no equation in the model
  IDV(20) random    enables spike channel 12, the reactor outlet flow

IDV(1) over 8 h: active at the last sample True, 8.00 h since onset

twenty 8-hour faulted runs: 19 completed, 1 tripped

digest:   7f9accc04bb61e7f
text:     tepsim.scenario.v1;seed=4651207995;hours=8;step=2.777777777777778e-4;every=180;faults=4;controlled=1;idv12=0;trip=1;continuous=0;integrator=euler;events=
parses back equal:      True
two runs bit-identical: True

The API

Three classes carry the whole surface, and they divide the way the problem does. A Scenario is a description and holds no state. A Simulation is the machinery. A Run is the result.

Module level

channel_names() returns the 53 names in row order, so a CSV header and a matrix column cannot disagree about which is which. faults() returns the twenty Fault records.

Scenario

Immutable and cheap to copy, so a caller can build one, derive variants from it, and keep all of them.

There are four constructors: the bare Scenario(...), Scenario.baseline(...), Scenario.fault(n, ...) and Scenario.from_text(text). All but the last take the same keyword arguments, which are seed, hours, step_hours, sample_every, controlled, driver_forces_idv12 and trip_ends_the_run, plus faults on the bare one.

The two defaults worth knowing are that a trip ends the run and that the driver does not force IDV(12). Both are Class C quirks of the original that this port fixes by default and can reproduce on request; see the delta register.

Simulation

Simulation(scenario) holds a plant, a controller stack and an integrator state, all of it owned. The original keeps its whole working set in six Fortran COMMON blocks, which allows exactly one simulation per process and no reentrancy at all. This allows as many as there are threads.

run() runs the whole scenario and returns a Run. It runs a copy, so calling it twice gives two equal runs rather than one run and one empty one, and it never raises for a plant that misbehaves: a trip or a failed temperature solve is reported through Run.outcome, because a run that ended early is data.

run() releases the GIL for the entire integration, which is where all the time goes, so an ensemble is a ThreadPoolExecutor and not a ProcessPoolExecutor. Nothing has to be pickled, and the plant does no I/O and touches no Python object while the GIL is down.

Run

Every array is a read-only view over one buffer, which is filled once when the run finishes and moved into NumPy rather than copied. to_numpy() returns the same object every call, and measurement, manipulated, column and columns are strided views into it. Call .copy() for something writable, or numpy.ascontiguousarray if contiguity matters. columns() is a dict in channel order, which makes it a pandas.DataFrame constructor argument as it stands.

labels() is the part with no counterpart in the original. A published Tennessee Eastman dataset is a matrix and a filename, so every detection-delay figure in the literature is computed against whatever onset its author assumed. Here the onset is in the data.

Fault

faults() returns twenty of these. published is the description at teprob.f:172-191 verbatim, five of which say only "Unknown". effect is what the source actually does, which for those five is perfectly explicit: only the physical interpretation was withheld. shape is 'step', 'random' or 'sticking', line names the teprob.f line the fault acts on, and affects_the_plant is False for the three sticking faults, which touch no equation in the model at all. The same table is in The twenty disturbances.

What Scenario(...) cannot say

Three of a scenario's eleven fields are out of the constructors' reach: the event schedule, the continuous extension that allows a fault at a fraction of its full strength, and the choice of integrator. Scenario.from_text is the way to them from Python, and since the text is the scenario's real serialisation, the thing the digest is taken over and the browser app puts in a URL fragment, building through it is building in the form the run will be described by.

text = tep.Scenario.baseline(hours=18).to_text()
scheduled = tep.Scenario.from_text(text.replace("events=", "events=6:start:4,12:stop:4"))

repr() knows about this. A scenario the constructor can express is printed as a constructor call, because that is what is readable at a prompt, and anything else is printed as Scenario.from_text(...). Both round-trip, so eval(repr(s)) == s holds for either shape. The alternative, printing only the constructor's arguments, produced valid Python that evaluated to a different scenario with a different digest and said nothing about it.

notebooks/04-custom-scenarios.ipynb works through all three fields, with a helper that rebuilds the text line from a dictionary rather than patching it, so that a mistyped field name is an error rather than a silently different run.

What is not in the package

tepsim-stats, the crate the repository's own validation ladder runs its statistics on, is development-only. cargo xtask ci asserts that no shipped crate so much as names it, and it is therefore not reachable from Python. That is deliberate: the ladder's statistics have to be free of any dependency the port itself does not have, so they can never be the reason a validation number changed.

For monitoring work the notebooks use notebooks/pcamon.py instead, which is PCA, both monitoring statistics and their control limits in NumPy and the Python standard library alone, with no SciPy and no scikit-learn. Where the two could differ they agree: 02-fault-detection-pca.ipynb checks every printed digit of its detector against the Rust one.

Types

The package ships py.typed, and _tepsim.pyi beside __init__.py describes the compiled half, so mypy and an editor see the real signatures rather than Any.

Worked examples

The four notebooks in notebooks/ are where the package is actually used. They are executed and committed with their outputs and plots, and rendered copies are published beside this book:

• Getting started
• Fault detection with PCA
• The three hard faults
• Custom scenarios

Those four links resolve on the published site, where .github/workflows/pages.yml renders the notebooks beside the book with jupyter nbconvert. In a local mdbook build they point at files that are not there; the sources are notebooks/*.ipynb in the repository, and notebooks/README.md says how to run them.

The tutorials in this book are the narrative around them.

The twenty disturbances

IDV(1..20) is a bare integer array in the original, and the header at teprob.f:172-191 names each entry in prose. Nothing there connects a name to the line that implements it, and five of the twenty are called only "Unknown". This page makes the connection: every fault names the teprob.f line it acts on, where it is injected, and what shape it has.

There are twenty, not the twenty-one of the later literature. teprob.f:340 is DO 500 I=1,20. The extra one comes from later versions of the model, and this port follows the vendored source.

Three shapes, and the third is not a plant disturbance at all

Step faults change a feed condition the moment they are switched on and hold it. Seven of the twenty, at teprob.f:407-414, 567 and 568.

Random faults enable a walk channel through IDVWLK (teprob.f:347-358). Ten of the twenty, of which the last three drive spike trains rather than walks.

Sticking faults do not touch the model. They set IVST (teprob.f:793-798), which widens the dead band a valve command must cross before the valve follows it. Three of the twenty.

That third kind matters more than its size. A sticking fault is not a disturbance to the plant; it is a disturbance to the controller's authority over the plant. In an open-loop run, where the command never moves, it does nothing whatever, and a scenario engine that treated it as a plant fault would report an injected disturbance with no effect and look broken. Tier 4 confirmed this from a direction that could not have been arranged: over a four-hour run, IDV(14), IDV(15) and IDV(19) report worst errors identical to the nominal case to every digit, because with the command held still their trajectory is the nominal trajectory.

The table

The channel column is not decoration. A test asserts that this table agrees with the code that maps IDV flags to channel flags, that every one of the twelve channels is driven by exactly one fault, that the spiking flag is set exactly for channels 10 and above, and that exactly IDV(14), IDV(15) and IDV(19) fail to reach the plant. Two statements of one fact, which is the point.

The five "Unknown" faults are not unknown

The header calls IDV(16) through IDV(20) unknown, and every paper on TEP repeats it. The source is perfectly explicit about what they do; only their physical interpretation was withheld. They enter the model at these points:

So IDV(19) is a sticking fault and the other four are not, which the shared label hides. Three of the four are the spike channels, which is why they are reported in the literature as the hardest to detect: they are intermittent rather than sustained.

Nine walks and three spike trains

The twelve channels are not the same kind of object, which the shared array names hide (teprob.f:340-406).

Channels 1 to 9 are random walks. When one runs out, teprob.f:359-371 evaluates the old segment at its endpoint, takes the value and the slope there, and builds a segment that continues smoothly from them. TESUB5 (teprob.f:1506-1537) chooses the next knot value and slope from the uniform generator and fits a Hermite cubic; TESUB8 (teprob.f:1300-1359) evaluates the cubic at the current time.

Channels 10 to 12 are spike trains, and teprob.f:372-396 gives them their own rule. They alternate between two states. Dwelling: the channel sits at zero for a randomly drawn interval, its segment being a parabola rising from zero, and the dwell ends when the value reaches 0.1. Spiking: once the value exceeds 0.1, the channel is given a cubic through the current value and slope that lasts exactly 0.1 hours, with coefficients that drive it hard and then back down. So a spike channel is off, off, off, then briefly on.

The flag scales the dwell, not the schedule. CDIST(I) = IDVWLK(I) / h^2 at teprob.f:391 is the only place the flag enters a spike channel. With the disturbance off it is zero, the parabola is flat at zero, the channel never reaches 0.1, and it dwells forever, drawing a fresh interval each time. With it on, the parabola climbs and the channel eventually spikes. Either way it keeps drawing at the same rate.

That last property is load-bearing and is easy to optimise away. IDVWLK multiplies the two endpoint draws at teprob.f:1529-1530 but not the duration draw at teprob.f:1528, so an inactive walk channel still consumes all three draws. Skipping the two endpoint draws when the flag is zero produces identical segment values, because an inactive channel lands on SZERO with zero slope either way, and leaves the generator two steps behind. Every subsequent draw in the run then differs: the noise, the other channels, everything. That mutation was implemented deliberately to check the tests have teeth, and it was caught on the stream position rather than on any value. It is the entire argument for Tier 3 demonstrated at Tier 1.

Using them

$ tep faults                       # the table above, from the source
$ tep run --fault 4 --hours 24 --labels

or, from Rust, Scenario::fault(4) and Scenario::baseline().with_fault(4). Faults can be combined; the flags are independent.

One caveat, though it is no longer the default. The original driver switches IDV(12) on at eight hours whatever the scenario asked for (temain_mod.f:366-368). That is delta D-011, and it is off here: Tier 7 established that the published files were generated with that line replaced rather than kept, so a request for IDV(4) gets IDV(4) and nothing else. Scenario::faithful() and the driver_forces_idv12 field turn it back on, and the ground-truth labels record it either way, so a run carrying it is visibly fault-free for eight hours and then not.

A first run, and the 53 channels

The worked example is notebook 1, Getting started. It runs the plant, plots what normal operation looks like, walks the twenty faults, injects one, reads the ground truth back, checks reproducibility from a seed, and finishes with a trip. Source: notebooks/01-getting-started.ipynb.

The whole library is three objects. A Scenario says what to simulate, a Simulation does it, and a Run holds what came out. The notebook uses all three in its first cell; this page is the part of the story that is worth stating once in prose rather than re-deriving from a plot.

The examples are Python, and The Python package is how to install it. The simulator itself is Rust, and its Rust API is Getting started, but a run started from either binding is bit-identical, because both call the same core.

Why forty samples in two hours, and not 7200

The integrator takes one step per simulated second, so two hours is 7200 steps. A sample is written every 180 of them, which is the three-minute spacing temain_mod.f:401 writes at and the spacing of the published d00 through d21 files. Two hours is therefore forty rows, and a 48-hour run is 960.

The first sample of a run is at 0.0497 hours rather than at 0.05, and that is not an off-by-one. run.steps[0] on that row is 180, and run.hours[0] is the time at which step 180 began, which is 179 seconds. The simulated clock is advanced at the end of a step, after the row has been written, because that is the order temain_mod.f writes in, and a row that carried the post-step time would be labelled with a clock the plant had not reached when it was measured. The last row of a two-hour run is at 1.9997 hours for the same reason. If you want the initial condition itself, it is the model's nominal state and not a row of the run.

What the 53 channels are

A row is the plant as an operator sees it: 41 measurements and then the 12 valve positions the controllers are holding. run.to_numpy() returns exactly that as one (n_samples, 53) array, measurements first, and channel_names() returns names in the same order, so a CSV header and a matrix column can never disagree about which is which.

The 41 measurements split into two groups that behave quite differently.

XMEAS(1) through XMEAS(22) are continuous instruments: flows, pressures, levels, temperatures, the compressor work. They are read every step and carry Gaussian measurement noise whose standard deviation is the XNS table in teprob.f.

XMEAS(23) through XMEAS(41) are the three gas chromatographs, and they are not instruments in the same sense at all. They sample on a schedule, take time to run, and then hold the answer until the next result arrives. The notebook plots this, and it is the single most surprising thing about the data the first time you see it: a composition channel is a staircase where a flow meter is a curve. The reactor feed and purge analysers run every 0.1 hours and the product analyser every 0.25, so at the three-minute output cadence each answer appears twice or five times. A detector that treats those repeats as independent observations is counting the same measurement several times, and that is worth knowing before it produces a p-value.

Both 0.1 literals in teprob.f are single precision, so the gas interval is really 0.10000000149011612 and a step landing on exactly 0.1 does not sample. That is faithfully reproduced here; see the delta register.

XMV(1) through XMV(12) are the manipulated variables. Eleven of them move. XMV(12), the agitator speed, sits at 50 forever, because the published control scheme never touches it. A statistical model fitted to all 53 channels has to cope with that constant column, and the detector tutorial shows what a PCA model has to do about it.

Getting at the data

run.to_numpy() is the record itself, one read-only (n_samples, 53) array of float64. Three views on top of it cover most uses. run.measurement(n) and run.manipulated(n) take the one-based indices of the original, so run.measurement(7) is XMEAS(7) and needs no mental arithmetic. run.column(i) takes the zero-based row position, and run.columns() returns all 53 at once keyed by name, which is the shape a covariance matrix, a chart or a pandas.DataFrame wants. Alongside them, run.hours and run.steps are the simulated time and the integrator step each row was written at.

Every one of those is a view into the same buffer rather than a copy, so .copy() is needed before writing into one.

The digest

scenario.digest is a content hash over everything the run's output depends on: the seed, the duration, the step, the cadence, the disturbances, the control mode, the quirk flags, the schedule and the integrator. Two scenarios that describe the same experiment produce the same sixteen characters, and two that differ in any respect do not.

It is worth putting in a filename or a file header, because it turns "this is the fault 4 data, I think" into something checkable. scenario.to_text() writes the whole scenario out in one line that Scenario.from_text reads back, so a dataset can carry its own description rather than a memory of one. The notebook checks both properties rather than asserting them: the same scenario run twice is bit-identical, a different seed is not, and the text round-trips.

What a trip looks like

The plant has eight shutdown conditions, on reactor pressure, reactor level, reactor temperature, separator level and stripper level. Cross one and the simulation is over. run() never raises for a plant that misbehaves, because a run that ended early is data, and throwing it away would hide the difference between a port that trips where the original does and one that does not. The outcome, the step, the hour and the condition that fired are all on the Run.

The notebook's survey of all twenty faults at 48 hours is the useful version of this: at the default seed exactly one of them trips, IDV(6), the total loss of the A feed. That is a property of the seed as much as of the fault, which is the kind of thing worth measuring before designing an experiment around it.

The same run from the command line

$ tep run --hours 2

writes the same 40 rows as CSV on stdout, with seventeen significant digits, so the file round-trips an f64 exactly.

Injecting a fault, and finding it in the data

The worked example is the second half of notebook 1, Getting started, which switches IDV(4) on eight hours into a 24-hour run, plots the reactor temperature and the cooling water valve side by side, and reads the ground truth back. Source: notebooks/01-getting-started.ipynb.

A disturbance is one line: Scenario.fault(n) is the baseline with IDV(n) switched on for the whole run, one-based exactly as the Fortran's IDV(n) is. That is what the original does, and it is the right thing when what you want is a faulted record. It is the wrong thing when what you want is to watch the fault arrive, and arrival is what a detector is judged on, so the notebook schedules it instead: IDV(4) off for eight hours and on afterwards, which is the layout the published test sets use.

IDV(4) is a five degree step in the reactor cooling water inlet temperature. The interesting part is that it is nearly invisible where you would first look.

The fault is in the valve, not in the temperature

Over the fourteen hours after the fault settles, the mean reactor temperature is 120.399542 degrees on the fault-free run and 120.399514 degrees on the faulted one. The difference is 28 millionths of a degree, against a fault-free standard deviation of 0.018728 degrees: two thousandths of one standard deviation. To any instrument, and to any detector watching that channel, nothing happened.

The cooling water valve tells the other half of the story. It sits at 41.103% open without the fault and 44.868% with it, a shift of 3.766 percentage points of valve travel. That is the whole of IDV(4): the cooling water arrives hotter, the temperature controller notices immediately, and it opens the valve until the temperature comes back. The controller has converted a disturbance in a measured variable into a disturbance in a manipulated one, which is what a controller is for.

The lesson generalises well past this fault, and it is the single most important thing to understand before building a monitor for a plant under closed-loop control. The evidence of a disturbance often sits in the manipulated variables rather than in the measurements, because the controller has been busy cleaning up the measurements. Both halves are in the array run.to_numpy() returns, and the detector in the next tutorial uses all 53 channels for exactly this reason.

Ground truth

run.labels() records what was actually true at each instant. It returns two (n_samples, 20) arrays indexed by IDV(n) - 1: active says whether that disturbance was on at that sample, and since_onset says how many hours it had been on, nan where it was not on at all.

The original records nothing of the sort. A published Tennessee Eastman dataset is a matrix and a filename, so every detection delay in the literature is measured against an onset its author knew from the experimental protocol rather than from the data. That is fine until two papers disagree about where sample 160 falls. Here the onset is in the data, and finding it is int(np.argmax(labels["active"][:, 3])) rather than a constant somebody has to remember.

The precision is deliberate too. A fault live from the first step reports 0.04972222222222225 hours at the first sample, not "about 0.05": the label carries the same simulated clock the row does, so it is 179 seconds. That is what lets a delay of one sample be distinguished from a delay of zero.

The disturbance you did not ask for, and no longer get

Pass driver_forces_idv12=True and run for longer than eight hours, and a second fault appears at hour eight without being asked for. That is not a bug in the scheduler. It is temain_mod.f:366-368, which switches IDV(12) on at eight hours whatever the scenario said, and both IDV(4) and IDV(12) act on cooling water, so a run nominally labelled IDV(4) was really the two together.

It was the default here until 2026-08-28, on the belief that reproducing d01 through d21 required it. Tier 7 showed the opposite: every dNN_te file except d12_te sits at the nominal operating point straight across row 160, which is hour eight, so the published files were made with that line replaced. The default now follows the files.

It is delta D-011 in the register, and it used to be the single most common way a comparison against the published files went quietly wrong, because it was the default. It is not any more: a request for IDV(4) gets IDV(4), and driver_forces_idv12=True is how to ask for the driver's version. The Rust equivalent is Scenario::faithful(), which sets that quirk and the freeze-on-trip one together, and the command line spelling is tep run --force-idv12. Say which one you used when you report numbers. The labels make the difference visible either way, which is the point of recording ground truth rather than assuming it.

From the command line

$ tep run --fault 4 --hours 8 --labels

--labels adds the fault and hours_since_onset columns to the CSV, so the ground truth travels with the data.

Building a detector, and measuring it

The worked example is notebook 2, Fault detection with PCA, which fits the model, plots both statistics against their limits on fault-free and faulted records, scores eight disturbances, and then diagnoses the one limit that does not work. Its sequel, notebook 3, The three hard faults, measures the benchmark's detectability floor three separate ways. Sources: notebooks/02-fault-detection-pca.ipynb and notebooks/03-hard-faults.ipynb.

This is the canonical Tennessee Eastman monitoring experiment. Fit a principal component model to a record of the fault-free plant, watch two statistics on new data, and raise an alarm when either leaves its control limit. Hotelling's T-squared measures distance from the training mean inside the subspace the model retained. The squared prediction error, written SPE or Q, measures how much of an observation the model could not reconstruct at all.

The detector is notebooks/pcamon.py, about four hundred lines of NumPy and the Python standard library with no SciPy and no scikit-learn: the eigendecomposition is numpy.linalg.eigh, the normal quantile is statistics.NormalDist, and the F quantile the T-squared limit needs is a continued-fraction incomplete beta with a bisection on top.

That file, rather than a crate, is deliberate. The repository has a Rust implementation of exactly this detector in tepsim-stats, and this page used to show it. tepsim-stats is development-only, and cargo xtask ci asserts that no shipped crate so much as names it, so a reader who installed tepsim and followed the page could not build what the page was teaching. Everything the notebooks do runs against the package you can install.

The two implementations agree. Notebook 2 opens by reproducing the Rust transcript, and every printed digit matches: the same 33 components, the same limits to three decimals, and the same four detection rates and delays, from a hand-written cyclic Jacobi sweep on one side and LAPACK on the other. What the notebook measures after that is therefore a property of the plant and the method, not of the linear algebra.

The measurement is the point

Any detector can be made to look good by reporting only the faults it catches. Every detection rate in the notebook is reported next to the false alarm rate on held-out fault-free data the model never saw, and one of the two statistics comes out badly.

Three numbers do the work, and each has a trap in it.

The fault detection rate is the fraction of post-onset samples that raised an alarm. It is a rate over samples, not a per-run yes or no, so a detector that catches a fault and then loses it scores badly. That is the intent: a statistic that drops back inside its limit while the fault is still running is flickering, not detecting. Its complement, the missed detection rate, is what the tables in the literature report.

The false alarm rate is the same fraction over the pre-onset samples, and it is the number that makes a detection rate meaningful. A detector with a 20% false alarm rate that achieves a 20% detection rate has detected nothing.

The detection delay is the number of samples from the onset to the first run of three consecutive alarms. The persistence requirement is not decoration. With a run length of one the delay is just the first alarm after the onset, and on a detector with any false alarm rate at all that is mostly luck: at a 3% false alarm rate the first post-onset sample alarms by chance one time in thirty, and calling that a delay of zero flatters the detector. The literature uses three and six and does not agree, so the run length travels with the number. Delays are in samples, and one sample is three minutes.

What the model does with a constant column

pcamon.fit standardises each column to zero mean and unit sample variance before forming the correlation matrix. Standardisation rather than mere centring is not optional here: reactor pressure lives near 2705 kPa and a composition is a percentage, so a covariance model would be a model of the pressure and nothing else.

How many components to keep is passed as a named rule rather than a bare k, because two detectors that retain different numbers of components are different detectors and a result that does not say which rule produced it cannot be reproduced. At 90% of the variance the notebook keeps 33 of 52 components. That is two thirds of them, and it is worth pausing on: the largest eigenvalue is 5.910, only 11.4% of the total of 52, so this plant has no small handful of dominant directions and a monitoring scheme on it works in a fairly high-dimensional retained subspace.

XMV(12), the agitator speed, never moves, so its standard deviation is zero and it cannot be standardised. Rather than divide by zero or quietly drop the column, pcamon records its index, zeroes its row and column of the correlation matrix, and says so when asked. A silent drop here is how a 53-variable model becomes a 52-variable model that nobody can reproduce.

The number the detector would rather you did not see

Both limits are drawn at 99% confidence, so the nominal false alarm rate is 1%. On a fresh 48-hour fault-free run the model never saw, T-squared alarms on 32 of 960 samples, a rate of 0.0333, which is high but recognisable. SPE alarms on 174 of 960, a rate of 0.1812. That is a factor of eighteen.

It is not an arithmetic error. pcamon.spe_limit computes the Jackson-Mudholkar expression exactly as stated. It is the assumption underneath the expression failing: the limit is derived for residuals that are normal and independent, and the plant's are neither. Several feed conditions are driven by slow random walks that never stop, so a 48-hour record wanders somewhere a 25-hour training record did not go, and when it does, every sample in the excursion alarms together. The notebook's plot shows exactly that shape, with the SPE alarms arriving in long blocks rather than as isolated points.

The tempting response is to lower the confidence level until the number looks right. The notebook does the experiment instead, holding everything else fixed and lengthening the training record:

That settles it. With enough fault-free training data both statistics land on the nominal 1%, so Jackson-Mudholkar was never the problem. A 25-hour record simply does not contain the tails of a process whose feed conditions are driven by random walks that never stop, and the limit it produces is therefore too tight.

This has a direct consequence for the literature, and not a comfortable one. The published training file d00 holds 500 samples, which at a three-minute interval is exactly 25 hours: the first row of that table. Every SPE false alarm rate reported for static PCA on the published Tennessee Eastman data inherits this.

Estimating the limit empirically from the training residuals is the usual advice, and the notebook checks that too rather than assuming it. It does not help: the empirical SPE limit gives 0.1448 against the analytic 0.1812, and the empirical T-squared limit is markedly worse than the analytic one at 0.1042 against 0.0333. Both are drawn from the same 25 hours that were too short in the first place, and neither can know about the excursions it never saw.

The fix is more data, or a method that models the serial correlation rather than assuming it away. Dynamic PCA, which augments each observation with lagged copies of itself, and canonical variate analysis are the two the literature reaches for, and Russell, Chiang and Braatz's 2000 comparison of both against static PCA is on exactly these files.

The floor under the benchmark

Notebook 3 is about the best known empirical result on this problem: IDV(3), IDV(9) and IDV(15) are effectively undetectable by these methods. It measures that on the published d00 through d21 files, on simulated runs at the seeds those files record, and over a ten-seed ensemble, and the three measurements agree.

The plainest form of the result is distributional. The median post-onset T-squared for those three faults is within one percent of the fault-free median, and their median SPE within four percent, on quantities whose fault-free values span two orders of magnitude. No threshold separates them because there is nothing to separate.

That is worth stating carefully, because it reads like a criticism of PCA and is not one. The plant really is running normally under those three disturbances. A detector that alarmed on them would be reporting a fault with no consequence, and the correct behaviour for a monitoring scheme is what it does. What the result measures is that this benchmark has a detectability floor set by the process rather than by the method, which is part of what makes it a good benchmark: any paper claiming to detect all twenty is claiming something about its false alarm rate that it has probably not measured.

A fault that arrives, and clears

The worked example is notebook 4, Custom scenarios, which schedules IDV(4) between hours six and twelve and plots the valve moving and coming back, composes two faults, sweeps a fault's magnitude from zero to one, and compares the three integrators. Source: notebooks/04-custom-scenarios.ipynb.

Every published Tennessee Eastman dataset is the same shape of experiment: set some IDV flags before the run, then leave them. That is all the original admits, which is why the literature's fault onsets are all in the same place and why almost nobody studies a fault that arrives, persists, and then goes away.

A Scenario here carries a schedule of events, each at a time and each doing one thing, so a disturbance can arrive at a stated hour, several can arrive independently of each other, and any of them can clear. Two further things the original cannot express sit beside it: a fault applied at a fraction of its full strength, and a choice of integrator. All three live in the scenario's canonical text form, which is where Python reaches them, and none of it would be worth much if the description of a run could not travel with the run.

The schedule is text, and that is the point

The Python constructors cover eight of a scenario's eleven fields. The other three, events, continuous and integrator, are reached by rendering a scenario to text, editing a field, and parsing it back. That sounds like a workaround and is closer to the opposite. The text is the scenario's real serialisation: it is what the digest is taken over, what the browser app puts in a URL fragment, and what the Rust and wasm sides read and write, so building a schedule through it is building it in the form the run will be described by.

An event is time:verb followed by the verb's own fields, so an events field of 6:start:4,12:stop:4 switches IDV(4) on at hour six and off at hour twelve. The notebook's helper rebuilds the whole line from a dictionary rather than patching it with a string replacement, so a mistyped field name is a KeyError and not a silently different run.

An event is applied on the first step whose time is at or past its own. The window is half-open at the start and closed at the end, so an event is applied exactly once however the step size divides its time, and an event at time zero is applied on the first step rather than never. Two events at the same instant keep the order they were written, because stop then start on one fault is a different scenario from the reverse and the digest has to be able to tell them apart.

What the trace shows

The reactor cooling water valve sits at 41.063% open before the fault. Half an hour after it arrives the valve is at 45.492%, it holds near 44.5% for the six hours the fault is live, and half an hour after it clears it is back at 41.045%, within a fiftieth of a percentage point of where it started. The controller does not know a fault happened in either direction; it is rejecting a disturbance both times.

That return is worth noticing rather than assuming. It is this loop doing its job on a step disturbance that was removed, and it is not guaranteed in general: a fault that shifts an inventory leaves the plant somewhere else even after it clears, and the plant after a fault is not the plant before it, because the controllers have moved. The recovery half of the problem is expressible here and barely studied, which is most of the reason the schedule exists.

The labels follow the schedule exactly. Once the fault clears, since_onset goes back to nan, which means the ground truth describes the plant's current condition rather than its history. A study of recovery has to look at transitions in active, not at that column.

Composition

Several disturbances can be described independently, each with its own arrival time. The original allows more than one IDV flag to be set, but not when each arrives, and interactions are usually where the interesting behaviour is.

IDV(1) steps the A/C feed ratio and IDV(8) enables random walks on the A and B feed compositions. Individually the plant absorbs both: reactor pressure means of 2709.58 and 2702.19 kPa against a fault-free 2706.89. Together the mean is 2731.92 with a standard deviation of 113.07, where adding the two individual shifts would predict 2704.87. The effect of the pair is not the sum of the effects of the parts, which is the whole reason to be able to compose them.

Half a fault

teprob.f:341-346 opens TEFUNC by forcing every IDV to zero or one, so the original has exactly two states per disturbance. The disturbances are then used multiplicatively, XST(1,4) = TESUB8(1,TIME) - IDV(1)*0.03 at teprob.f:407 being the pattern, so a magnitude between zero and one scales a fault smoothly and needs no other change to the model.

That is the continuous extension, and it is off by default. With it on, a sweep of IDV(4) from magnitude 0 to 1 moves the mean cooling water valve from 41.023% to 44.779% in steps of about 0.9 percentage points, linear in the magnitude to within 0.042 percentage points of valve travel, while the reactor temperature it is defending stays at 120.400 degrees at every magnitude. A detection threshold study is then a sweep over one number rather than an argument about what "harder" means.

Two things are worth saying plainly. A run with the extension on is not comparable to any published dataset, and none of the validation ladder applies to it. And a fractional magnitude without the extension is refused rather than rounded, with an error naming the line of Fortran that makes it impossible. Silently turning a request for half a fault into a whole one would produce a run that does not match its own description, which is exactly what the content hash exists to prevent. A magnitude of exactly 1.0 is bit-identical to the faithful path, which is asserted, so the extension can be left on for a study that mixes full and partial faults.

The trap in refining the step

The notebook also compares Euler, RK4 and Dormand-Prince, and the result is about the original rather than about the port: RK4 and Dormand-Prince agree with each other to about 2e-6 relative while both differ from Euler by about 1.5e-2 on the worst channel. Two independent methods agreeing that closely and disagreeing with the third is what convergence looks like. Euler is not the accurate choice, it is the faithful one, and everything the validation ladder claims is a claim about Euler.

The obvious next move, keeping Euler and halving the step, does not do what you expect. Halving the step moves the answer by more than changing the method does, and it does not converge as the step shrinks. That is neither a bug nor stiffness: the disturbance walks and the measurement noise advance once per step, so a run at half the step draws twice as many random numbers and is a different realisation of the stochastic forcing. The step size is part of the disturbance model in this plant and not only a numerical parameter. Change the method to integrate the same realisation more accurately; change the step only when you mean to change the noise.

The description travels with the run

A scenario, schedule included, writes out as one line of text, that line parses back to an equal scenario, and the digest is unchanged across the round trip. A dataset generated from a scenario can carry the line in its header, and a reader can reconstruct the exact run rather than the description of a run. The same line is what TEP Studio puts in a URL fragment, so a scheduled experiment can be handed to somebody as a link.

The format is versioned and its parser is strict: a missing field, an unknown field, a value out of range or a version this build does not know is an error that says what was wrong, rather than a default quietly substituted.

repr() handles the two shapes separately, and knowing why is useful. A scenario the constructors can express prints as a constructor call, because that is what is readable at a prompt. Anything else prints as Scenario.from_text(...), which round-trips by construction. Both satisfy eval(repr(s)) == s. The obvious implementation, printing only the constructor's arguments, produced valid Python that evaluated to a different scenario with a different digest and said nothing about it.

A schedule holds up to 32 events, which is far more than any experiment in the literature uses. The published datasets have exactly one event each: the fault, switched on before the run.

The plant

This page is the vocabulary the rest of the process chapters use: the eight components, the fifty states, and the thirteen internal streams. All three are recovered from teprob.f rather than from the 1993 paper, because on two of the three the paper and the source disagree.

The flowsheet is the map for all of it. The mixing zone, the reactor, the separator and the stripper are the four vessels that hold state, and the state table below is their contents; the condenser and the compressor hold none. Every solid line is one of the thirteen FTM entries the stream table enumerates, and every tag is an XMEAS or XMV index whose teprob.f line is cited on the instrumentation page. Nothing on the drawing was taken from the paper's figure. Where a stream carries two numbers they are both shown, because that disagreement is the subject of the third section of this page.

Three things on it are worth reading twice, because each is a place the source and the received description of the plant part company. The mixed A and C feed, stream 4, does not reach the mixing zone at all: it enters the stripper base as the stripping gas (teprob.f:614-662). The pressure published as stripper pressure, XMEAS(16), is the mixing zone's PTV (teprob.f:694), because the model carries no separate stripper vapour space. And FTM(12), the liquid that fails to strip out and falls into the stripper sump, has no number in the paper at all.

Components

Eight, A through H. A, B and C are non-condensible and are treated as ideal gases throughout (teprob.f:478). D through H are condensible and get an Antoine vapour pressure (teprob.f:484). B is the inert: it appears in none of the four reactions, arrives with the mixed feed, and leaves only through the purge.

A, B and C have no real liquid density correlation. AD is 1.0 with BD and CD zero for all three (teprob.f:973, 983, 993), so they contribute a flat, temperature-independent term. That is a placeholder keeping the mixing rule finite rather than a fitted number, because the model never puts them in a liquid phase in quantity.

The fifty states

The original carries the state in a bare YY(50) and unpacks it by index arithmetic inside TEFUNC (teprob.f:417-440). Recovering that mapping is the first act of the port, and here it becomes typed structure, pinned by a test that reads the corresponding COMMON/TEPROC/ variables back out of the Fortran rather than trusting a comment.

The eight slots do not mean the same thing in every vessel, and this is the part that is easy to get wrong. For the reactor and the separator the array is split by phase: slots 1 to 3 are the vapour holdups of A, B and C, and slots 4 to 8 are the liquid holdups of D through H. UCLR(1..3) is set to zero at teprob.f:420-421 because the non-condensibles never form a liquid, and UCVR(4..8) does not come from the state at all: it is derived from the vapour-liquid equilibrium later in the same call (teprob.f:500-501). For the stripper all eight slots are liquid, and for the mixing zone all eight are vapour.

The four temperatures are state, not derived quantities

TESUB2 takes its temperature argument as both the initial guess and the result (teprob.f:1432, 1438), and the four call sites at teprob.f:460-465 pass TCR, TCS, TCC and TCV straight out of COMMON. Every evaluation therefore starts its Newton solves from the previous evaluation's answers, and since the iteration stops on a step below 1e-12 the converged value depends on where it started.

That is not a detail. B-0015 measured the cost of getting it wrong: seeding the solves from a different point on the nominal trajectory moves up to 21 of the 50 derivatives. A port that solved from a fixed guess would be tidier and would not be bit-exact. The warm-start temperatures are carried explicitly here, and B-0034 found the same thing again from the other end, where a trajectory started from the nominal literals instead of from the values TEINIT's own evaluation leaves behind is a different trajectory rather than a rounding of the same one.

Those four values are asserted against the oracle bit for bit; they are from the LOG.org entry for B-0052.

The thirteen internal streams

The Fortran's stream indices are not the stream numbers in the paper. FTM(1) is the D feed, which Downs and Vogel call stream 2. FTM(3) is the A feed, which they call stream 1. Nothing in the source says so, and every reimplementation of TEP has to rediscover it; getting it wrong produces a plant that runs, looks plausible, and is wired up incorrectly.

The mapping was established from the source, not from the paper: teprob.f:565 drives FTM(1) from valve 1 and XMV(1) is documented as "D Feed Flow (stream 2)"; teprob.f:567 gates FTM(3) on IDV(6), documented as "A Feed Loss (Stream 1)"; teprob.f:688 reports FTM(10) as XMEAS(10), "Purge Rate (stream 9)"; teprob.f:683 reports FTM(9) as XMEAS(5), "Recycle Flow (stream 8)"; and so on for the remaining six.

Streams 6 and 7 are the same fluid. teprob.f:656-661 copies flow, enthalpy, temperature, composition and component flows from 6 to 7 wholesale, with no mixing, no pressure drop and no heat loss. Stream 7 exists so that the reactor's balance at teprob.f:763-772 can name its own inlet.

Vessel volumes

Four vessels, four fixed total volumes, all four written in the original as single-precision literals (teprob.f:1118-1121):

The reactor and the separator hold two phases, so their vapour space is whatever the liquid does not occupy. The stripper is treated as liquid only and the mixing zone as vapour only.

Precision is a property of each literal

The line above is not pedantry. 182 of the assignments in TEINIT are single-precision literals, and a literal written without a D suffix is stored by gfortran as a single-precision value widened to double, which differs from the decimal number by up to about 6e-8 relative. Since the port must reproduce the original's arithmetic bit for bit, every constant has to be transcribed according to the suffix on its own line.

The original is not consistent, so the precision cannot be inferred from elsewhere in the file. teprob.f:1411 writes 273.15 and teprob.f:594 writes 273.15D0. The 1.8 at teprob.f:790 and 792 is single while every other occurrence in the file (1396, 1404, 1464, 1471) is 1.8D0. The gas constant at teprob.f:475 is RG=998.9, single, and it multiplies six of the eight partial pressures in every vessel.

The canary the constants table was built around is XMW(2), which must come out 25.399999618530273 and not 25.4.

The right-hand side

TEFUNC occupies teprob.f:196-816, six hundred lines that present themselves as a derivative evaluation. They are not one. Understanding why is the single most consequential structural decision in this port, so it comes before the unit operations rather than after them.

Three phases, because TEFUNC is not a pure function of (t, y)

Inside what looks like a right-hand side, the original also advances the disturbance random walks, draws measurement noise, ticks the three sampled analysers, and latches the valve commands. That is harmless for the fixed-step Euler integrator the original uses, which evaluates the right-hand side exactly once per step. It is wrong for anything else: an RK4 step would advance the walks four times and draw four sets of noise.

The impure work does not sit in one place, which is the part that is not obvious until you look. Some of it must happen before the derivative and some can only happen after it. The walks are read at teprob.f:407-416, so they must be advanced first. The measurement vector is assembled at teprob.f:679-701 out of flows the evaluation computes, so noise cannot be added until afterwards. One impure call cannot sit on both sides.

The port therefore splits the routine in three:

The pure phase hands back the signals it computed alongside the derivative, so the post-phase does not have to re-run the model to find out what to add noise to.

The valve latch is hoisted, and the hoist is proved mechanically

teprob.f:793-798 sets IVST from IDV and 799-804 latches VCV from XMV, at the very end of the routine. That block reads only XMV, VST, IVST, IDV and TIME, and nothing in 345-792 writes any of them, so moving it into the pre-phase changes no number.

That is a claim about four hundred and fifty lines of Fortran, which is too large to check by eye, so it is checked by machine instead: a dedicated oracle test drives both orderings and asserts they agree. The latch shares the DO 9020 loop with the valve derivative at teprob.f:805, so the port splits that loop, sending the latch to the pre-phase and leaving YP(I+38) in the pure one.

What each stage of the pure phase does

The pure phase runs in the original's order, and each block is a module in tepsim-core with its own page or section in this chapter.

The fifty balances

Everything above exists to feed teprob.f:762-811. Four vessels, each with eight component balances and one energy balance, plus two cooling-water wall temperatures and twelve valve lags.

\[ \frac{dn_i}{dt} = \sum_{\text{in}} \dot n_i - \sum_{\text{out}} \dot n_i + r_i \]

\[ \frac{dE}{dt} = \sum_{\text{in}} h F - \sum_{\text{out}} h F + Q \]

The reactor is the only vessel with a reaction term. The cooling water walls follow

\[ \frac{dT_w}{dt} = \frac{F_w \times 500.53 \times (T_{in} - T_w) - Q \times 10^6 / 1.8}{H_w} \]

where 500.53 converts a cooling water flow to a heat capacity rate and the factor \(10^6 / 1.8\) undoes the scaling the enthalpy correlations carry (teprob.f:789-792). The 1.8 on those two lines is single precision and is the only such occurrence in the file, which is why it has to be read off the line rather than inferred from teprob.f:1396, 1404, 1464 or 1471, where it is written 1.8D0.

Each valve is a first-order lag toward its latched command (teprob.f:805):

\[ \frac{dv_i}{dt} = \frac{c_i - v_i}{\tau_i} \]

A shutdown freezes the plant

teprob.f:807-811 zeroes all fifty derivatives whenever any shutdown condition holds. That does not stop the plant, it freezes it: the state stops moving, the clock keeps running, and nothing in the original says so.

PLAN.org classes this as a Class C quirk, "behaviour-defining and benchmark-relevant", so the fix needed a measured delta and a sign-off. Both arrived on 2026-08-28, and a default Scenario now ends the run at the trip. It is delta D-007. The port reports the trip and its cause either way, rather than leaving a caller to infer a freeze from a vector of zeros.

Scenario::faithful() reproduces the freeze, and Tier 2 needs it to: the adversarial sampling pool contains states that trip, and a port that did not freeze would disagree with the oracle on all fifty components for every one of them. So does any comparison against d06 or d18, whose published files are between 45% and 76% frozen tail.

Determinism, and the two libm builds

tepsim-core is no_std, forbids unsafe, and contains no f32, no SIMD, no rayon, no reordered reductions and no source of time or randomness outside the model's own generator. exp, pow and ln come from a vendored pure-Rust libm so that the answer does not depend on the host's C library.

That choice costs something measurable, and the project measures it rather than hoping. Over the range of Antoine arguments this model reaches, the vendored libm and gfortran's disagree on 9.945% of them, by exactly one ULP (LOG entry B-0018). So a differential against the default build can only assert a tolerance. Every Tier 2 comparison therefore runs a second time under a libm-system feature, where exp is the one gfortran calls, and is held to zero ULP there. Both runs are in the CI gate. Without the second, a reassociation worth one or two ULP would pass silently.

Integer powers are a related trap and are not routed through pow at all. gfortran expands X**4 into multiplications, and the shape of the expansion is load-bearing. Measured over 200,000 values with this project's pinned flags (LOG entry B-0023): (x*x)*(x*x) matches gfortran on 200,000 of 200,000, ((x*x)*x)*x on 132,040, and pow(x, 4.0) on 99,523. So it is binary exponentiation, squaring twice, and the two plausible alternatives are each wrong a third to half of the time.

The mixing zone

The feed mixing zone is a single vapour volume in which the three pure feeds, the compressor recycle and the stripper overhead combine before entering the reactor. It is the simplest of the four vessels: one phase, a fixed volume, no reaction, no heat duty.

Source: teprob.f:428 and 434 for its states, teprob.f:465-466 for its temperature, teprob.f:492 for its pressure, teprob.f:576-579 for its outlet flow, and teprob.f:783-788 for its energy balance.

Equations

All eight components are vapour and the volume is fixed at VTV, so the equilibrium block needs only the ideal gas law applied to the mixture (teprob.f:492):

\[ P_v = \frac{N_v R T_K}{V_{tv}} \]

The temperature comes from the specific internal energy, not from a state directly. The block totals the holdups, forms the mole fractions and the energy per mole, and then solves for whatever temperature makes the mixture's specific internal energy equal to it (teprob.f:465-466):

\[ N = \sum_i n_i, \qquad x_i = \frac{n_i}{N}, \qquad e = \frac{E}{N} \]

That solve is Newton's method in TESUB2 (teprob.f:1415-1442), warm-started from the previous evaluation's answer, which is why TCV is state rather than a derived quantity. See The plant.

Flow out of the mixing zone is not valve-driven. It is a square-root resistance across the pressure difference to the reactor, converted from mass to moles by the stream's mean molecular weight (teprob.f:576-579):

\[ F_6 = \frac{1937.6 \, \sqrt{\max(P_v - P_r,\, 0)}}{\overline{M}_6} \]

The clamp at zero is what stops a reversed pressure gradient from producing a NaN out of the square root, and it is reachable from an adversarial state though not from the nominal trajectory.

The component and energy balances have five inlets and one outlet (teprob.f:762-770 and 783-788):

\[ \frac{dn_i}{dt} = \dot n_{i,1} + \dot n_{i,2} + \dot n_{i,3}

• \dot n_{i,5} + \dot n_{i,9} - \dot n_{i,6} \]

\[ \frac{dE}{dt} = h_1 F_1 + h_2 F_2 + h_3 F_3 + h_5 F_5 + h_9 F_9 - h_6 F_6 \]

There is no Q term: the mixing zone is adiabatic.

Variables

Three things the source settles

The mixed A/C feed does not pass through here. Stream 4 goes directly to the stripper, and it appears in the stripper's energy balance at teprob.f:778-782 and in the stripper's feed at teprob.f:614-662, never in the mixing zone's. The three feeds that do enter are streams 1, 2 and 3, the D, E and A feeds.

Stream 7 is an alias of stream 6, made in the stripper block. teprob.f:656-661 copies flow, enthalpy, temperature, composition and component flows from 6 to 7 wholesale. There is no mixing, no pressure drop and no heat loss between them; stream 7 exists so that the reactor's balance at teprob.f:763-772 can name its own inlet.

The pressure published as stripper pressure is this vessel's. teprob.f:694 reports PTV as XMEAS(16), which Downs and Vogel's measurement table names stripper pressure. The model carries no separate stripper vapour space: the stripper's overhead discharges into the mixing zone as stream 5, and PTV is the pressure of that shared vapour node. The mixing zone's own outlet flow is reported as XMEAS(6), the reactor feed rate (teprob.f:684).

The reactor

The reactor is a two-phase vessel with an internal cooling coil and an agitator. Four exothermic gas-phase reactions run in its vapour space, its liquid holdup sets how much of the coil is wetted, and it is the only vessel in the plant with a reaction term in its balances.

Source: teprob.f:473-502 for the vapour-liquid equilibrium, teprob.f:503-528 for the kinetics, teprob.f:663-673 for the coil, and teprob.f:762-772 for the balances.

Vapour-liquid equilibrium

The vapour space is what the liquid does not occupy (teprob.f:473):

\[ V_{vr} = V_{tr} - V_{lr} \]

A, B and C are non-condensible and are treated as ideal gases, so their partial pressures come from the holdup directly (teprob.f:478-483):

\[ p_i = \frac{n_i R T_K}{V_{vr}}, \qquad i \in \{A, B, C\} \]

D through H are condensible, so their partial pressures come from Raoult's law with an Antoine vapour pressure in degrees Celsius (teprob.f:484-491):

\[ p_i = x_i \exp\!\left(A_i + \frac{B_i}{T_c + C_i}\right), \qquad i \in \{D \ldots H\} \]

The total is the sum of all eight, the vapour composition is \(y_i = p_i / P\) (teprob.f:493-496), and the vapour holdup follows from the ideal gas law applied to the mixture (teprob.f:497 and 500):

\[ N_v = \frac{P V_{vr}}{R T_K}, \qquad n_i = N_v y_i \]

UCVR is both an input and an output

teprob.f:418 fills UCVR(1..3) from the state and teprob.f:500 fills UCVR(4..8) from the equilibrium computed here. It is one Fortran array written by two different mechanisms, and the halves are not interchangeable: the non-condensibles are integrated, the condensibles are derived. Read in that order it also explains why the A, B and C partial pressures are computed first, at teprob.f:478-483: they need UCVR before it is overwritten.

This is where bit equality with gfortran ends

DEXP at teprob.f:485 and teprob.f:488 is the model's first transcendental call. The port answers it from the vendored pure-Rust libm rather than the platform's, for the determinism reason set out in The right-hand side. Measured over the whole Antoine range this model reaches, the two disagree on 9.945% of arguments, by exactly one ULP.

Everything downstream of a condensible partial pressure therefore carries about 1.1e-16 of relative difference from the Fortran that no care in the algebra removes. The bit-exactness claim does not disappear, it moves: under the libm-system feature the transcendental is the one gfortran calls, and the port is bit-identical again, so the algebra is still held to zero ULP rather than to a tolerance.

Kinetics

1:  A + C + D -> G          2:  A + C + E -> H
3:  A + E     -> F          4:  3 D       -> 2 F

Rates 1 and 2 are Arrhenius in reactor temperature with fractional pressure orders on A and C, multiplied by a disturbance drift factor (teprob.f:503-504, 508-511):

\[ r_1 = f_1 \, e^{\,a_1 - E_1/T_K} \; p_A^{1.1544} \, p_C^{0.3735} \, p_D \, V_{vr} \]

\[ r_2 = f_2 \, e^{\,a_2 - E_2/T_K} \; p_A^{1.1544} \, p_C^{0.3735} \, p_E \, V_{vr} \]

Both are guarded: teprob.f:507 requires \(p_A > 0\) and \(p_C > 0\), and sets \(r_1 = r_2 = 0\) otherwise. Rates 3 and 4 are first order in each reactant, and rate 4 shares rate 3's exponential rather than having one of its own (teprob.f:505-506, 516-517):

\[ r_3 = e^{\,a_3 - E_3/T_K} \, p_A \, p_E \, V_{vr}, \qquad r_4 = 0.767488334 \; e^{\,a_3 - E_3/T_K} \, p_A \, p_D \, V_{vr} \]

All four are multiplied by the vapour volume at teprob.f:518-520, so a rate is an extent in moles per hour rather than a volumetric rate. Net production per species follows the stoichiometry (teprob.f:521-527), and the heat release comes from reactions 1 and 2 only (teprob.f:528):

\[ Q_{rxn} = r_1 h_1 + r_2 h_2 \]

with \(h_1 = 0.06899381054\) and \(h_2 = 0.05\) (teprob.f:1122-1123). Reactions 3 and 4 contribute no heat: the original simply does not include them in RH, and HTR(3) is declared but never assigned or read.

R1F and R2F are two different quantities under one name

They arrive from TESUB8(7) and TESUB8(8) at teprob.f:415-416 as the IDV(13) kinetics-drift multipliers, are consumed at teprob.f:503-504, and are then reassigned in place at teprob.f:508-509 to hold the fractional pressure powers. The two meanings share nothing but the storage.

Reading teprob.f:510 as though R1F were still the drift factor gives a plausible and completely wrong rate law, so the port gives the two roles separate names. That is delta D-002: no numerical effect, and the only defence against a misreading no test would catch, because a wrong-but-consistent reading still reproduces itself.

CRXR(2) is never assigned

Seven of the eight slots are written at teprob.f:521-527. CRXR(2), the inert, is not, and it is read anyway at teprob.f:763. It works because COMMON is zero-initialised and nothing ever writes it, so B's net production is zero by static initialisation rather than by statement. That is delta D-003, class A: the value is right, the mechanism is an accident. Here the slot is explicitly zero and a test asserts the oracle agrees.

Precision hazards in this range

Every literal except 0.767488334D0 and 1.5D0 is single precision: the three pre-exponentials, the three activation energies, the gas constant 1.987, and both fractional exponents.

Worse, 40000.0/1.987 at teprob.f:503 is a quotient of two single-precision literals, so Fortran evaluates the division itself in single precision. Widening the operands first and dividing in double is wrong by 4e-9 relative, inside a DEXP argument.

The cooling coil

The coil's effective area ramps with liquid level, because the coil is only wetted over part of its height. VLR/7.8 is the level as a percentage, and the ramp is piecewise linear between 10% and 50% (teprob.f:663-669):

\[ \lambda = \begin{cases} 1 & \ell > 50 \\ 0 & \ell < 10 \\ 0.025\,\ell - 0.25 & \text{otherwise} \end{cases} \qquad \ell = V_{lr} / 7.8 \]

The overall coefficient is quadratic in agitator speed (teprob.f:670-671), and the duty is the coefficient times the driving temperature difference, scaled by a disturbance drift factor (teprob.f:672-673):

\[ U A_r = \lambda \left(-0.5\,\omega^2 + 2.75\,\omega - 2.5\right) \times 855490 \times 10^{-6} \]

\[ Q_r = U A_r \, (T_w - T_c) \, (1 - 0.35 \, d_{10}) \]

That parabola peaks at \(\omega = 2.75\), which is above the agitator's whole range: teprob.f:575 puts it between 1.5 and 2.5. So the coefficient rises monotonically with speed everywhere the plant can go, from 0.5 to 1.25 times the scale, and the falling half of the parabola is unreachable. Its roots are at 1.149 and 4.351, both outside that range too, so the coefficient never reaches zero from the agitator alone. The model is a fit, not a mechanism.

The ramp does not quite meet its flat sections

0.025 at teprob.f:668 is single precision, so it is stored as 0.02500000037252903 and the ramp misses both of its endpoints:

Both breakpoint comparisons are strict, so a level of exactly 10 or exactly 50 takes the ramp, and UARLEV is discontinuous by those amounts as the level crosses either one. This is faithful reproduction rather than a delta: it is what the original computes, the gaps are eight orders below the quantity itself, and the coefficient they scale is an empirical fit in the first place. It is written down because "the ramp meets the flat sections" is the obvious assumption, it is false, and a test asserting it would fail for a reason that looks like a porting error.

Balances

Inlet is stream 7, outlet is stream 8, and the reaction term is the only one of its kind in the plant (teprob.f:762-772):

\[ \frac{dn_i}{dt} = \dot n_{i,7} - \dot n_{i,8} + r_i, \qquad \frac{dE}{dt} = h_7 F_7 - h_8 F_8 + Q_{rxn} + Q_r \]

QUR is heat removed, so it enters positive here because \(U A_r (T_w - T_c)\) is already negative when the coil is cooling. The reactor's cooling water wall temperature is YP(37) (teprob.f:789-790).

Variables

Three of the eight shutdown conditions belong to this vessel: reactor pressure above 3000 kPa gauge (teprob.f:703), liquid volume outside 2 to 24 cubic metres (teprob.f:704-705), and temperature above 175 degrees Celsius (teprob.f:706). See Instrumentation.

The condenser and separator

The condenser cools the reactor effluent and the separator splits it into a vapour and a liquid. The vapour leaves through two paths, the compressor recycle and the purge; the liquid goes to the stripper. The compressor and its recycle valve belong here too, because their whole job is deciding how much of the separator's vapour goes back around the loop.

Source: teprob.f:473-502 for the equilibrium, teprob.f:585-601 for the flow network and the compressor, teprob.f:674-676 for the condenser duty, and teprob.f:773-777 for the balances.

Vapour-liquid equilibrium

The separator's equilibrium has exactly the shape of the reactor's and shares its code path. The vapour space is VTS less the liquid volume (teprob.f:474); A, B and C get ideal gas partial pressures (teprob.f:481); D through H get Raoult's law with an Antoine vapour pressure (teprob.f:488-490); the composition is \(y_i = p_i / P\) (teprob.f:495); and the vapour holdup comes back out of the ideal gas law at teprob.f:498 and 501.

The one number worth keeping in mind is that PTS floors around 811 mmHg over the whole sampled domain. That matters for the purge clamp below.

The condenser

A smooth saturating function of reactor outlet flow, approaching 0.404655 as the flow grows (teprob.f:674), with the duty taken against the stream 8 temperature rather than the separator's own, and scaled by a disturbance drift factor (teprob.f:675-676):

\[ U A_s = 0.404655 \left(1 - \frac{1}{1 + (F_8/3528.73)^4}\right) \]

\[ Q_s = U A_s \, (T_{ws} - T_{st,8}) \, (1 - 0.25 \, d_{11}) \]

**2 and **4 are integer powers and must not go through pow

gfortran expands an integer exponent into multiplications rather than calling libm, and the shape of that expansion is load-bearing. Measured over 200,000 values with this project's pinned flags:

So it is binary exponentiation, squaring twice, and the two plausible alternatives are each wrong about a third and a half of the time. X**2 is x*x on all 200,000, which is the only thing it could be.

Three ways out

The underflow to the stripper is valve-lagged, and is the simple case (teprob.f:570):

\[ F_{11} = \frac{v_7 R_7}{100} \]

The purge is pressure-driven to atmosphere through valve 6 (teprob.f:585-588):

\[ F_{10} = \frac{v_6 \times 0.151169 \times \sqrt{\max(P_s - 760,\, 0)}}{\overline{M}_{10}} \]

That clamp at zero is the one clamp in the whole flow network that cannot be reached. PTS is a sum of eight partial pressures in a vessel that always holds material, and it floors around 811 mmHg against a threshold of 760, so no trajectory state, no random perturbation and no adversarial boundary takes that branch. It is therefore covered by a unit test at a composition chosen to reach it, rather than by the differential, on the principle that a branch no test enters is indistinguishable from a branch that is wrong.

The recycle goes through the compressor. The operating point on the curve is the pressure ratio between the mixing zone and the separator, clamped at both ends (teprob.f:589-591), and the machine is fixed-speed with a cubic pressure-ratio curve (teprob.f:592-593). The recycle valve then bleeds flow back and the result has a floor (teprob.f:596-599):

\[ \dot m = \max\!\left( F_{\max}\left(1 + \frac{1 - r^3}{1.197}\right) - v_5 \times 53.349 \times \sqrt{\max(P_v - P_s,\, 0)}, \; 10^{-3} \right) \]

with \(r = \mathrm{clamp}(P_v / P_s,\, 1,\, 1.3)\), \(F_{\max} = 280275\) (teprob.f:1170) and the ratio ceiling CPPRMX at 1.3 (teprob.f:1171). The floor at \(10^{-3}\) exists so that the division at teprob.f:600-601 cannot blow up.

PR**3 at teprob.f:593 is an integer power, for the same reason **4 is above, and is written out as three multiplications rather than routed through pow.

The compressor work appears twice: as an enthalpy bump on the recycle stream, and as measurement 20 (teprob.f:594-595, 601, 699):

\[ W = \dot m \, (T_{cs} + 273.15) \times 1.8 \times 10^{-6} \times 1.9872 \times \frac{P_v - P_s}{\overline{M}_9 P_s} \]

Note that the 273.15 on that line is written 273.15D0, double precision, unlike the one in TESUB2 at teprob.f:1411. The original is not consistent about that constant and each occurrence has to be read off its own line.

HST(10) = HST(9) is a snapshot, not an alias

teprob.f:562 copies the separator vapour enthalpy into the purge. Both streams leave the separator vapour space, so at that moment they are the same fluid at the same temperature and the copy is exact.

Then teprob.f:601 adds the compressor work to HST(9). The recycle gains the work; the purge does not, because it was copied first. Reading line 562 as an alias rather than as a copy would give the purge a share of compressor work it never receives, and the two lines are seventy apart, so the ordering is easy to miss. Both energy balances that read stream 9 come after line 601 and therefore see the bumped value.

Balances

One inlet, three outlets, and the condenser duty (teprob.f:762-770 and 773-777):

\[ \frac{dn_i}{dt} = \dot n_{i,8} - \dot n_{i,9} - \dot n_{i,10} - \dot n_{i,11} \]

\[ \frac{dE}{dt} = h_8 F_8 - h_9 F_9 - h_{10} F_{10} - h_{11} F_{11} + Q_s \]

The condenser cooling water wall temperature is YP(38) (teprob.f:791-792).

Variables

Two of the eight shutdown conditions belong to this vessel: separator liquid volume above 12 or below 1 cubic metre (teprob.f:707-708).

The stripper

Steam strips light components out of the separator liquid. What leaves overhead rejoins the mixing zone as stream 5; what does not becomes the stripper's own liquid, stream 12, which leaves as the product, stream 13. The vessel is treated as liquid only: it has no vapour space of its own in the model.

Source: teprob.f:614-662 for the column, teprob.f:677-678 for the reboiler, and teprob.f:778-782 for the energy balance.

Equations

The feed is the mixed A/C feed plus the separator underflow (teprob.f:635-639). Note that stream 4 arrives here directly rather than through the mixing zone:

\[ f_i = \dot n_{i,4} + \dot n_{i,11} \]

A vapour-to-liquid ratio sets how hard the column strips, scaled by a temperature factor (teprob.f:622):

\[ \Lambda = \frac{F_4}{F_{11}} \, \tau(T_c) \]

Each condensible then strips according to a Langmuir-shaped saturating function of that ratio (teprob.f:623-627):

\[ s_i = \frac{k_i \Lambda}{1 + k_i \Lambda}, \qquad i \in \{D \ldots H\} \]

and the split is simply (teprob.f:643-644)

\[ \dot n_{i,5} = s_i f_i, \qquad \dot n_{i,12} = f_i - \dot n_{i,5} \]

Both product streams leave at the stripper's own temperature (teprob.f:652-653), and their enthalpies are taken on different bases: stream 5 with ITY = 1, the vapour basis, and stream 12 with ITY = 0, the liquid one (teprob.f:654-655).

The temperature factor has a pole at 177 C

\[ \tau(T) = \begin{cases} T - 120.262 & T > 170 \\ 0.1 & T < 5.292 \\ \dfrac{363.744}{177 - T} - 2.22579488 & \text{otherwise} \end{cases} \]

from teprob.f:615-621. The middle branch diverges at 177 C, which is inside the range the two outer branches leave for it only if TCC exceeds 170, and it does not: the T > 170 branch takes over first. So the pole is unreachable by seven degrees, and the two branches are continuous to within 0.1% at 170.

The adversarial state catalogue built for Tier 2 places a state at 176 C anyway, to sit near the pole and confirm that it stays on the linear branch. That state is coverage of the guard, not of the pole.

FTM(11) > 0.1 switches the whole block

Below that threshold the column is not really running, and teprob.f:629-633 substitutes five fixed stripping factors (0.9999, 0.999, 0.999, 0.99, 0.98) rather than evaluating the correlation. The reason is visible in the arithmetic: \(\Lambda = F_4 / F_{11}\) diverges as \(F_{11} \to 0\).

Both sides are covered by the adversarial catalogue, which places a state exactly on FTM(11) = 0.1. Since the test at teprob.f:614 is .GT., that state takes the fixed-factor branch.

SFR(1..3) are never recomputed

teprob.f:623-627 and 629-633 both write slots 4 through 8 only. Slots 1, 2 and 3 are set once in TEINIT (teprob.f:1126-1128) and are read at teprob.f:643 on every evaluation, so A, B and C strip at a fixed 99.5%, 99.1% and 99.0% no matter what the column is doing.

That is the intended physics rather than an oversight: the non-condensibles are gases, they leave overhead essentially completely, and no temperature or flow ratio in the plant's range would change that. It is worth stating because the loop at teprob.f:643 runs I=1,8 and looks as though all eight factors come from the branch above it.

The reboiler

Steam is at 100 C, so above that there is nothing to transfer and the original sets the duty to zero rather than letting it go negative (teprob.f:677-678):

\[ Q_c = \begin{cases} U A_c \, (100 - T_c) & T_c < 100 \\ 0 & \text{otherwise} \end{cases} \]

The coefficient UAC is not computed here. It is a valve-lagged capacity with a disturbance drift factor, set at teprob.f:572:

\[ U A_c = \frac{v_9 R_9 (1 + d_9)}{100} \]

which is the point at which IDV(16), published as "Unknown", enters the model.

The nominal trajectory sits near 65 C, so the cutoff is the branch at risk of never being exercised. B-0021 measured 300 of 300 nominal states below 100.

The reactor inlet is an alias, and this is where it is made

teprob.f:656-661 copies flow, enthalpy, temperature, composition and component flows from stream 6 to stream 7 wholesale. There is no mixing, no pressure drop and no heat loss between them: stream 7 exists so that the reactor's balance at teprob.f:763-772 can name its own inlet. It lives in the stripper block for no reason other than that is where the original put it.

Balances, and an asymmetry between them

The component balances are a straight pass-through of the two streams the column produced (teprob.f:762-770), while the energy balance names the column's inputs instead (teprob.f:778-782):

\[ \frac{dn_i}{dt} = \dot n_{i,12} - \dot n_{i,13} \]

\[ \frac{dE}{dt} = h_4 F_4 + h_{11} F_{11} - h_5 F_5 - h_{13} F_{13} + Q_c \]

That is what the source does, and the two forms describe the same vessel: the component split at teprob.f:643-644 conserves moles by construction, so streams 4 and 11 in, less stream 5 out, is stream 12.

Variables

Two of the eight shutdown conditions belong to this vessel: stripper liquid volume above 8 or below 1 cubic metre (teprob.f:709-710). Its level measurement is also the odd one out among the three, because its span is the vessel volume VTC itself rather than a separately hard-coded range (teprob.f:693).

Instrumentation

Forty-one measurements come out of the plant, and they are not all the same kind of thing. XMEAS(1..22) are continuous instruments, read every step, computed at teprob.f:679-701 and given additive noise at teprob.f:711-735. XMEAS(23..41) are composition analysers, read on a schedule and reporting the composition from their previous sample (teprob.f:736-761).

Nothing in the continuous block is physics. Every one of the twenty-two lines takes a quantity the model has already computed and converts it into the unit an operator's instrument would read. Getting a conversion wrong changes no state and no derivative; it changes only what the controller sees, which is worse, because the plant then runs correctly and is controlled wrongly.

The conversion factors

So FTM * 0.359 / 35.3145 is lbmol/h to standard cubic metres per hour, and (P - 760)/760 * 101.325 is mmHg absolute to kPa gauge. Levels are reported as a percentage of a span, and the three vessels do not do it the same way: the reactor and the separator carry hard-coded ranges (teprob.f:686, 690) while the stripper's span is the vessel volume VTC itself (teprob.f:693).

XMEAS(20) is assigned twice

      XMEAS(20)=CPDH*0.0003927D6
      XMEAS(20)=CPDH*0.29307D3

at teprob.f:698-699. The first is dead, and the two factors are not equal: 392.7 against 293.07, a third apart. So this is not a harmless duplicate but a superseded conversion, and a port that took the first line would report compressor work 34% high. That is delta D-006.

The shutdown detector

Eight limits, checked at teprob.f:702-710:

The original records only that something tripped, in a single integer ISD. This port reports which, because "the plant tripped" without a reason is nearly useless to a caller and the information is free.

All eight comparisons are strict, so a state exactly on a limit does not trip. That matters for how the adversarial sampling pool is built: states placed on the limits exercise the not-tripped side, and the tripping side needs states past them. Both were built.

Two of the eight are phrased in terms of XMEAS rather than the underlying quantity. teprob.f:703 tests the converted reactor pressure against 3000 kPa gauge, and teprob.f:706 tests XMEAS(9), which is TCR unconverted. Testing PTR against an equivalent mmHg threshold instead would be arithmetically different in the last bits.

What a trip does is described in The right-hand side: it freezes all fifty derivatives (teprob.f:807-811), which is delta D-007.

Noise and dead time

Noise is drawn by TESUB6 (teprob.f:1538-1546), twelve uniform draws summed and scaled, and is skipped entirely at TIME = 0 and on a tripped plant (teprob.f:711). Only the continuous noise is skipped, though. The analyser blocks at teprob.f:744-761 have no such guard, so a tripped plant still draws: 258 draws in a tripped evaluation against 522 in a healthy one, measured in B-0027. A port that silenced everything on a trip would leave the generator 264 steps behind and desynchronise every later draw.

The dead time is a latch, and the order of two lines makes it:

      XMEAS(I)=XDEL(I)
      CALL TESUB6(XNS(I),XMNS)
      XMEAS(I)=XMEAS(I)+XMNS
      XDEL(I)=XCMP(I)

The reported value is taken from the store before the store is updated. Swapping those two lines gives an analyser with no dead time at all, which produces entirely plausible numbers and a plant that is much easier to control than the real one.

Three further details are about when rather than about what. The schedules advance from their own previous value rather than from the current time (TGAS = TGAS + 0.1 at teprob.f:751), so a step arriving late does not shift the schedule. Both 0.1 literals, at teprob.f:741 and 751, are single precision, so the gas interval is 0.10000000149011612 and a step landing on exactly 0.1 does not sample; 0.25 is exactly representable, so the product analyser is unaffected, which is precisely the kind of inconsistency that has to be read off the line rather than inferred from its neighbour. And at TIME = 0 the analysers are primed rather than sampled (teprob.f:736-743): the store and the reported value are both set to the current composition, with no noise and no draw.

The 53 channels

Measurements and manipulated variables together are the 53 columns every downstream consumer sees, in this order. The measurement names follow Downs and Vogel's Table 4 and the manipulated ones their Table 3.

XMEAS(23..28) and XMEAS(29..36) are the two gas analysers, on a 0.1 hour schedule; XMEAS(37..41) is the product analyser, on 0.25 hours.

XMV(12), the agitator, is never written by any controller the driver calls, so in a closed-loop run it has zero variance. That is a fact about the control scheme rather than about the plant, and it is the reason the Tier 5 harness had to learn to report a degenerate ensemble as NaN rather than as a number.

Validation

This is the part of the project that determines whether anyone trusts the result, so it gets the most care. The strategy is a ladder: prove the pieces exactly, prove the derivative to near machine precision, prove the stochastic call order exactly, then prove the long-run behaviour statistically and, finally, prove it on the downstream task people actually care about.

Every number on this page was measured, and each one names the LOG.org iteration that measured it. None of them is a target, a plan, or a rounded recollection. The project's operating rule is to record numbers rather than verdicts, precisely so that a degradation inside tolerance is still visible: if a maximum relative error moves from 3e-14 to 8e-13, both pass a 1e-12 gate and something has broken, and the only place that is visible is the logged history.

The figures are held to the same rule. Each one is drawn by cargo xtask validate from the run's own output, committed alongside the generated chapters, and captioned with the condition that would make it false. Nothing in one is placed by hand: a marker is orange because its value is on the wrong side of a gate, never because a test's name was recognised, so a genuine regression and a deliberate positive control are drawn identically and the caption is what tells them apart.

All numbers below were produced with gfortran 15.2.0 and the pinned rustc 1.97.1. The oracle's compiler flags are fixed in build.rs and asserted by a test; changing them invalidates every Tier 1 and Tier 2 number here, so it is a logged re-baseline rather than a casual edit.

The oracle harness

tepsim-oracle is a development-only crate whose build.rs compiles the unmodified teprob.f and temain_mod.f with gfortran and links them through a small C shim exposing TEINIT, TEFUNC, each TESUBn, and read and write access to every COMMON block. A Rust test can therefore set the Fortran into an arbitrary state, call it, and compare against the Rust implementation in the same process.

That crate is never a dependency of tepsim, tepsim-py or tepsim-wasm. It runs on Linux and macOS runners where gfortran is available, and building it first is what converts the whole port from "read carefully and hope" into a differential testing exercise.

Alongside it sits a cheaper check that needs no Fortran at all. cargo xtask fidelity runs the port forward 100 steps from the nominal state and diffs states, derivatives, measurements and the generator word against a golden trace committed to the repository. It takes about a second and runs at the top of every session. Since B-0026 it has reported 100 of 100 steps diffed, worst 3.521e-14 at YP(12) step 77 against a 1e-12 gate, and that same number has now been recorded unchanged for eleven consecutive iterations, most recently in the entry for B-0052.

Tier 1: the utility routines, exactly

TESUB1 (enthalpy), TESUB2 (temperature from enthalpy by Newton), TESUB3 (heat capacity) and TESUB4 (liquid density) are swept over a simplex grid, ten million random Dirichlet samples, and a boundary pool, at every temperature in the physical range, for each of the three ITY modes. The gate PLAN.org sets is a maximum relative error below 1e-13 with a ULP histogram reported rather than a pass or fail.

The measured result is not "inside 1e-13". It is zero.

That is 59,924,940 evaluations across TESUB1 and TESUB3 with zero differing bits, p50 = p90 = p99 = p100 = 0 in all six sweeps. A separate test asserts bit equality directly, because a 1e-13 threshold would let a drift to 1e-15 pass unnoticed forever.

The 1e-13 gate is therefore not what is holding the line, and the entry for B-0009 says why that is the right expectation: both routines are straight-line arithmetic over constants already proved bit-identical, so once the association and the literal precisions are right there is nothing left to differ by. Any future Tier 1 routine that lands at 1e-15 rather than at zero has something wrong with it that a tolerance would hide.

The sabotage check. Substituting a double-precision 273.15 for the widened f32 at teprob.f:1411, and changing nothing else, gives a maximum relative error of 1.597e-5 at grid#704 T=21.875, 103,098,852,352 ULP (B-0009). That is the size of the error a single mis-transcribed literal produces, and it is why constants in this project are transcribed and asserted rather than retyped.

TESUB2 is the Newton solve, and it is measured from two starting strategies (B-0011). A warm start, which is what every call site actually does, converges in one step and recovers the answer exactly: round-trip error 0e0. A cold start from the far end of the range is what exercises the iteration, and its worst round trip is 5.12e-13 C, at dirichlet#720561 T=171.86, face#1309309 T=146.56 and face#96 T=144.14. Over 59,924,940 solves there were zero differing bits and zero abandoned iterations, which is the measurement that makes delta D-001's effect zero on the physical domain.

TESUB7, the random generator, must be exact and is (B-0005): 10^7 draws from the compiled-in seed with the XOR fold and the final state both matching; 10^6 draws each from five dataset seeds; draw-by-draw comparison over 200,000 draws for four seeds with every draw and every intermediate state bit-identical; and interleaved output modes over 50,000 draws on an irregular pattern. Its exactness rests on reproducing the rounding rather than removing it: the product exceeds 2^53 on 0.7716 of draws, against 0.7728 predicted from the arithmetic, so a "fixed" integer recurrence diverges at draw 0.

TESUB5 and TESUB6 are exact in both libm configurations, since neither touches a transcendental, and their draw counts were recovered independently by stepping a port-side generator from the word before each call to the word after: 3 draws for TESUB5, 12 for TESUB6, on every case and for both flag values (B-0028).

Tier 2: single-step derivative equivalence

Both implementations are forced into an identical state, all fifty states, the twelve manipulated variables, the twenty IDV flags, the full walk state, the generator word and the nineteen held analyser readings, evaluated once, and compared on all fifty components. Sampling is from three pools: states along the nominal closed-loop trajectory, random perturbations of those states scaled across several orders of magnitude, and adversarial states placed deliberately at every discontinuity and clamp in the model.

The tolerance is relative to the scale of the terms, not to the result, and that is a decision with a measurement behind it rather than a relaxation. A balance is inflow minus outflow, and near steady state those nearly agree. YP(2), the inert's reactor balance, is a difference of two flows around 660 whose result is a few parts in ten thousand of either. One ULP of difference in each term, which is all the vendored libm costs, is 1e-16 of the terms and 1e-4 of the result. Measured against the result, 28 of the 50 components exceed 1e-12 while the whole right-hand side is bit-identical to gfortran under libm-system. Each balance therefore reports the magnitude of its largest term alongside its value, and the gate is the error over that.

Acceptance, from B-0026: all fifty components, 2,412 running states, all three pools, worst 6.093e-14 at YP(7), perturbed#300, against a 1e-12 gate.

The whole tier fits in one picture. Every comparison it makes, run twice, once against each libm, plotted at its own maximum relative error:

28 of 50 components cancel by more than 100x, and the 22 that do not are exactly 0.000e0: bit-identical rather than merely inside the gate. The acceptance test asserts that contrast and fails if the count of heavily cancelling components drops below twenty, so the reasoning behind the decision expires loudly if the model ever stops cancelling.

Tier 3: RNG call-order equivalence

Both sides are instrumented to emit every draw, and the traces are diffed. This test exists specifically because it is the one the existing Python port would fail: its documented divergence is attributed to "the exact sequence of calls differs due to implementation details", which is exactly the defect a trace diff catches on the first run instead of after a 48-hour statistical comparison.

From B-0029, draw counts at four points in a run, with exact agreement at every one:

Scalings in one real evaluation: 30 signed and 432 unit. The 30 is 9 x 3 + 3, the walk advance; the 432 is 36 compositions at twelve draws each. Worst evaluation is 462 draws against a trace buffer capacity of 4096. The instrument is checked against the generator word rather than only against the values, which covers completeness, ordering and fidelity at once, and it holds over 113,088 draws.

Two independent methods agree on the counts. B-0027 measured them with no instrumentation at all, by stepping a port-side generator from the word before a call to the word after, and the trace lengths match that census exactly. Either method alone would be a claim; the two agreeing is evidence.

Tier 4: trajectory equivalence, diagnostic

Tier 4 is a diagnostic, not a gate. Long-horizon divergence between two programs that use different exp implementations is expected. PLAN.org asks for two things: that the error stay below the corresponding measurement noise standard deviation XNS(i) for at least the first several hours, and that the onset of divergence be explained by showing that switching libm moves it.

Open loop, nominal, 8 simulated hours (28,800 steps), from B-0034:

All 21 scenarios at 4 hours each stay within XNS for the whole run in both configurations, and every one is exactly 0.00e0 on the platform libm. Worst error by scenario on the vendored libm at 4 hours: nominal 2.45e-5, IDV(10) 3.66e-8, IDV(8) 1.42e-8, IDV(13) 1.05e-8, and the rest below 1e-8.

The explanation is better than the one asked for. With identical transcendentals the two trajectories are not merely close, they are bit-identical for 28,800 steps on all 41 measurements. So the divergence is transcendental rounding and nothing else, and the residual under the vendored libm is twenty thousand times below what the instruments could resolve after eight hours.

Closed loop, 48 hours, 172,800 steps, nominal scenario with the driver's forced IDV(12) from hour eight, from B-0041:

Bit-identical for all 172,800 closed-loop steps under the platform libm. Not one measurement, not one valve, not one step. Under the vendored libm the error never reaches a tenth of a billionth of any instrument's noise over two simulated days.

Compare the open-loop figure above: nominal at 8 hours ends at 5.149e-5 of XNS, six orders of magnitude worse. Closing the loop suppresses the amplification, which is what a controller pulling toward a setpoint should do and is worth having measured rather than assumed.

The same run measures what the control layer is actually buying. Open loop from the same start, with the valves held, trips at step 11,017 (3.060 h) on reactor pressure high. So "48 hours without tripping" is a statement about the control layer, not about a placid plant. B-0052 reproduced that trip from the public API and the CLI, at exactly the same step, without either being told about the other.

Tier 5: statistical equivalence, the real gate

Tier 5 tests equivalence rather than difference, because a failure to reject the null of no difference is not evidence of equivalence and two one-sided tests are. Every statistic ships in Rust, in tepsim-stats, with known-answer tests: Welch's t and the TOST wrapper, the two-sample Kolmogorov-Smirnov statistic, energy distance, autocorrelation, Welch power spectra and the Pearson correlation matrix. Nothing calls out to numpy or scipy.

The harness

From B-0047a: a 48-hour run produces 960 samples of 53 variables, a port run takes 766 ms in release, and a full battery of 2100 runs takes about 27 minutes per source. Over three scenarios the two sources agree bit-identically (0) under the platform libm and to 2.056e-13 under the vendored one, while five nominal seeds spread by 1.366e1, which is the scale that makes those first two numbers meaningful.

The same entry checked that the scenarios are not vacuous. All twenty disturbances move the plant within two hours; the smallest departures are IDV(3) at 4.076e-3, IDV(9) at 3.531e-3 and IDV(15) at 9.839e-3, and the largest are IDV(6) at 3.467e0 and IDV(1) at 2.456e0.

The battery

Four of the six statistics have margins that are measured rather than chosen (B-0047b). The reference's seeds are split in half, the statistic is computed against itself over twenty deterministic splits, and the cross-source value is treated as one more draw from that null, giving p = (1 + #{within >= cross}) / (K + 1) gated at 0.05. With twenty splits the cross-source value fails exactly when it is the strict maximum of the twenty-one.

The smoke battery in the CI gate is 3 scenarios, 4 seeds, 2 hours, 12 runs per source, 14 seconds. The full battery is partial: 3 of 21 scenarios at 100 seeds by 48 hours, stopped by direction after about 25 minutes.

Every calibrated statistic passed at p = 1.0000 on all three: the cross-source value was never the maximum of its null. The correlation matrix's Frobenius distance is three to five orders of magnitude inside the within-source spread, which matters because that matrix is exactly what a PCA-based detector consumes.

That relationship, rather than any single number, is what Tier 5 claims, so it is worth drawing:

That plot answers "are the two sources closer to each other than the reference is to itself", which is the calibrated question. It does not answer "how much room is left", and neither does a TOST verdict: a p-value reports that a test did not reject, not by what distance. The margin is a stated quantity, a tenth of the reference's standard deviation for each variable, so the distance can simply be drawn.

TOST power against battery size, measured in the same entry:

A power above 1 means the battery is underpowered for the margin, not that the sources differ. PLAN.org sets the mean margin at a tenth of the pooled standard deviation, and a disturbed plant has a much larger pooled spread than a quiet one, so the same absolute run-to-run wander sits comfortably inside the margin under IDV(2) and outside it at the nominal operating point. At the nominal plant the margin needs about 255 seeds, not 100. The variables that run out of power are the manipulated ones: an integrating controller's output has a random-walk component, so its mean over 48 hours varies between seeds by much more than a tenth of its own within-run spread. That is a fact about the plant, not about the port, and the battery reports it as undecided rather than as a difference while still asserting on the measured gap.

Physics invariants

These are the only tests in the whole ladder that can catch an error the port faithfully inherited, because every other tier compares against the Fortran and an error in the Fortran is invisible to all of them. From B-0046:

All four reactions balance exactly with the published molecular weights, `2 + 28

• 32 = 62, 2 + 28 + 46 = 76, 2 + 46 = 48and3 x 32 = 2 x 48`, which is why I-1 is an equality rather than a tolerance.

The invariants' teeth were checked by mutation, and one result is worth repeating because it is not obvious. Deleting the reaction term from the reactor's component balance leaves the total mass balance passing at 2e-16. That is not a weak test, it is a true fact about the invariant: I-2 is the molecular-weight-weighted sum of I-4, and I-1 says the reaction is mass-neutral, so the reaction term cancels out of I-2 exactly. Total mass conservation is blind to stoichiometry by construction. The unweighted per-component balance, I-4, sees it immediately. An invariant that is a sum of other invariants is weaker than the set, and how much weaker is not obvious from reading it.

Tiers 6 through 10

These have not run. They are listed here so that the shape of the remaining claim is visible rather than implied.

Tier 6, downstream-task equivalence, is the operational definition of "practically equivalent" for this project's research audience: train a detector suite on Fortran d00 data, evaluate on Fortran and on Rust test data, then reverse it, and compare detection rate, false alarm rate and detection delay. The claim to be able to make is that cross-source performance matches within-source performance within its own run-to-run variability. Backlog item B-0050.

Tier 7, published dataset reproduction, attempts direct reproduction of the bundled d00 through d21 files under the documented generation protocol, reporting per-file agreement with the Tier 5 machinery. The groundwork exists: teprob.f:1187-1256 carries fifty-four generator words in comments, one per published dataset, and B-0047a transcribed and asserted them. Three facts about that table are worth knowing in advance. Thirty-four of the fifty-four exceed 2^32, which is not a transcription error because TESUB7 reduces any seed on the first draw. Twenty-seven are even, and a multiplicative generator modulo a power of two keeps the factors of two its seed has, so those runs have a shorter period and low bits that never move. Reproducing the published files means doing the same rather than fixing it. Backlog item B-0051.

Tier 8 is differential fuzzing, Tier 9 is cross-platform determinism by golden BLAKE3 digest including wasm in a real browser, and Tier 10 requires that every quirk fix ship with a measured delta from the full Tier 5 battery with the fix on and off. None has started.

Two bugs that only long runs found

Worth recording because they are the argument for running the battery long rather than often (B-0047b).

TRCN, the Tier 3 trace counter, is a Fortran INTEGER that nothing clears between evaluations. At about 264 draws per step over 172,800 steps it passes 2^31 after roughly fifty runs, goes negative, passes the capacity guard, and writes outside the array. The symptom was a SIGSEGV deep inside the Fortran, tens of millions of steps into the battery, with nothing nearby to suggest why: the same seed ran perfectly in isolation, and starting at seed 40 moved the crash to seed 90, which is what identified it as cumulative rather than data-dependent. Nothing in Tiers 1 to 4 was ever affected, because no shorter run approached the overflow.

welch_t computed (n - 1) on a summary with no observations. In debug that panics; in release it wraps and returns a degrees-of-freedom figure that looks like a number. The case is XMV(12), the agitator, which no controller ever writes, so every run has zero variance and the log-variance ensemble is empty.

Both are the same lesson. A validation harness is code, it has bugs, and the bugs it has are the ones only its longest runs reach.

Validation, measured

This page is generated. cargo xtask validate wrote it from commit 03ec3d6-dirty. Every number on it was captured from that run's own output. To change what it says, change what the suite measures and run the command again.

The narrative version of this material, with the reasoning behind each tier and the history of what it caught, is in Validation. This section is the other half: the numbers, written by the command that ran the suite, from the suite's own output. Nothing here was transcribed.

A tier with no chapter has not been generated. That is stated rather than left to be inferred from an absence, because a missing page and a page nobody updated look the same from the table of contents.

Toolchain on the machine that ran this: rustc 1.97.1 (8bab26f4f 2026-07-14), gfortran 15.2.0.

Fidelity preflight

cargo xtask fidelity runs the port forward from the nominal state and diffs states, derivatives, measurements and the generator word against a golden oracle trace committed to the repository. It needs no Fortran toolchain, so it runs everywhere in about a second.

What the long tiers look like

Tiers 4 and 5 run but write no chapter yet, for the reason GENERATED_TIERS gives in xtask. Their figures are here, each drawn by the run named under it.

The ladder

This run selected tier(s) [5] and wrote chapter(s) for []. Tiers 4 to 7, 9 and 10 run but do not write a chapter yet; tier 8 has no harness. The delta index is generated separately, by cargo xtask deltas.

Tier 1: the utility routines

This page is generated. cargo xtask validate --tiers 1,2,3 --smoke wrote it from commit b6ef4ee-dirty. Every number on it was captured from that run's own output. To change what it says, change what the suite measures and run the command again.

TESUB1 (enthalpy), TESUB2 (temperature from enthalpy by Newton), TESUB3 (heat capacity) and TESUB4 (liquid density) are swept against the Fortran over a simplex grid, a Dirichlet sample and a boundary pool, at every temperature in the physical range, for each of the three ITY modes. PLAN.org sets the gate at a maximum relative error below 1e-13, with a ULP histogram reported rather than a verdict.

Reduced volume. This run passed --smoke, so the sweeps are the short ones the CI gate uses rather than the full ones PLAN.org specifies. The case counts in the tables below are what actually ran. Drop --smoke for the gate volume.

Produced with rustc 1.97.1 (8bab26f4f 2026-07-14), gfortran 15.2.0. The oracle's compiler flags are fixed in crates/tepsim-oracle/build.rs and asserted by a test; changing them invalidates every number on this page, which is why it is a logged re-baseline and not an edit.

What ran

2 test binaries: 7 test(s) passed, 0 failed, 0 ignored.

Figures

Measurements

15 block(s), lifted from the transcripts below. The columns are whatever fields the run printed, so a new field in the reporter becomes a new column here rather than data this page drops.

The test column matters as much as the numbers, because not every row is the port being measured. Some tests deliberately mis-type a constant, or solve from the wrong guess, to show what that would cost; a row from one of those is supposed to be enormous, and its test name says so. The what column carries a tier1 prefix in every tier, because it is the shared comparison reporter's own label rather than a claim about which tier printed it.

Transcripts

Each block below is a test binary's own output, verbatim, with the command that produced it. The summary above is derived from these; they are not derived from it.

tier1_enthalpy, vendored libm

cargo test -p tepsim-oracle --features oracle --release --test tier1_enthalpy -- --nocapture --test-threads 1

running 5 tests
test both_routines_are_bit_identical_to_the_fortran ... sweep: SMOKE, 7950 cases. The 9987490-case gate is `cargo xtask validate --tiers 1`.
ok
test reading_the_offset_as_double_precision_would_fail_the_gate_by_orders ... tier1 TESUB1 ity=2, offset as f64
  cases          : 7950
  max rel err    : 1.597e-5 at grid#704 T=21.875
  max ulp        : 103098852352 at grid#704 T=21.875
  ulp percentiles: p50>=16 p90>=16 p99>=16 p100>=16
  ulp histogram  : >=16:7950
  non-finite     : 0 seen, 0 mismatched
ok
test tesub1_matches_the_fortran_over_the_full_sweep ... sweep: SMOKE, 7950 cases. The 9987490-case gate is `cargo xtask validate --tiers 1`.
tier1 TESUB1 ity=0
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
tier1 TESUB1 ity=1
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
tier1 TESUB1 ity=2
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
ok
test tesub3_matches_the_fortran_over_the_full_sweep ... sweep: SMOKE, 7950 cases. The 9987490-case gate is `cargo xtask validate --tiers 1`.
tier1 TESUB3 ity=0
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
tier1 TESUB3 ity=1
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
tier1 TESUB3 ity=2
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
ok
test tesub4_matches_the_fortran_over_the_full_sweep ... sweep: SMOKE, 7950 cases. The 9987490-case gate is `cargo xtask validate --tiers 1`.
tier1 TESUB4
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched  [0.0 s]
ok

test result: ok. 5 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.01s

tier1_temperature, vendored libm

cargo test -p tepsim-oracle --features oracle --release --test tier1_temperature -- --nocapture --test-threads 1

running 2 tests
test tesub2_matches_the_fortran_and_the_silent_failure_never_fires ... sweep: SMOKE, 7950 cases. The 9987490-case gate is `cargo xtask validate --tiers 1`.
tier1 TESUB2 ity=0 start=warm
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 0e0 C
  elapsed        : 0.0 s
tier1 TESUB2 ity=0 start=cold
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 5.684341886080802e-14 C at grid#3637 T=131.25
  elapsed        : 0.0 s
tier1 TESUB2 ity=1 start=warm
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 0e0 C
  elapsed        : 0.0 s
tier1 TESUB2 ity=1 start=cold
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 1.9895196601282805e-13 C at grid#4617 T=170
  elapsed        : 0.0 s
tier1 TESUB2 ity=2 start=warm
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 0e0 C
  elapsed        : 0.0 s
tier1 TESUB2 ity=2 start=cold
  cases          : 7950
  max rel err    : 0.000e0 at grid#0 T=0
  max ulp        : 0 at grid#0 T=0
  ulp percentiles: p50=0 p90=0 p99=0 p100=0
  ulp histogram  : 0:7950
  non-finite     : 0 seen, 0 mismatched
  abandoned      : 0 (delta D-001)
  round trip     : max |solved - true| = 5.115907697472721e-13 C at face#96 T=144.14305056234923
  elapsed        : 0.0 s
ok
test the_fortran_silently_returns_the_guess_where_the_port_reports_failure ... D-001 demonstrated:
  Fortran: returned 120.4 silently
  port:    Newton did not converge in 100 iterations from a guess of 120.4 C: reached -48598544002334720 C with a final step of 24299272000832196, which is not below 1e-12
ok

test result: ok. 2 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.03s

Tier 2: the plant model

This page is generated. cargo xtask validate --tiers 1,2,3 --smoke wrote it from commit b6ef4ee-dirty. Every number on it was captured from that run's own output. To change what it says, change what the suite measures and run the command again.

Both implementations are forced into an identical state, evaluated once, and compared on all fifty derivative components. Sampling is from three pools: states along the nominal closed-loop trajectory, random perturbations of those states, and adversarial states placed at every discontinuity and clamp in the model.

Every comparison runs twice. The vendored libm disagrees with gfortran's by an ULP on about a tenth of exp and pow calls, so the default build can only be held to 1e-12; the libm-system build removes the transcendental from the comparison and is held to bit equality. Both runs are below, and the libm column says which is which.

The tolerance is relative to the scale of the terms rather than to the result. A balance is inflow minus outflow, and near steady state those nearly cancel, so an error that is 1e-16 of either term can be 1e-4 of their difference.

Reduced volume. This run passed --smoke, so the sweeps are the short ones the CI gate uses rather than the full ones PLAN.org specifies. The case counts in the tables below are what actually ran. Drop --smoke for the gate volume.

Produced with rustc 1.97.1 (8bab26f4f 2026-07-14), gfortran 15.2.0. The oracle's compiler flags are fixed in crates/tepsim-oracle/build.rs and asserted by a test; changing them invalidates every number on this page, which is why it is a logged re-baseline and not an edit.

What ran

20 test binaries: 94 test(s) passed, 0 failed, 0 ignored.

Figures

Measurements

143 block(s), lifted from the transcripts below. The columns are whatever fields the run printed, so a new field in the reporter becomes a new column here rather than data this page drops.

The test column matters as much as the numbers, because not every row is the port being measured. Some tests deliberately mis-type a constant, or solve from the wrong guess, to show what that would cost; a row from one of those is supposed to be enormous, and its test name says so. The what column carries a tier1 prefix in every tier, because it is the shared comparison reporter's own label rather than a claim about which tier printed it.