Tutorial: a reactor–mixer plant with hidden conditions¶

In this tutorial you build a small plant from the process library, run it through a four-day production plan, and find out why the product misses its plan even though every controller holds its setpoint. The files are in examples/reactor_mixer/ in the repository.

You will learn how to:

• connect units through their ports;
• make one flow follow another with a ratio station;
• use a reactor whose product depends on states nobody measures;
• describe a production plan and the disturbances nobody plans;
• compare what the historian records with the true values.

The plant¶

  A ── FIC-101 ──────────────┐
         │ FT-101 × ratio    ├──> TIC-201 reactor ── product (q) ──┐
  B ── FIC-102 (CAS) ────────┘    (lab AT-201)                     ├──> MX-301 mixer ──> AT-301
  C ── FIC-103 ────────────────────────────── material C ──────────┘

Materials A and B react in a stirred reactor, TIC-201. A mixer, MX-301, blends the reactor product with a third material, C. The operators set four values: the flow of A, the ratio of B to A, the reactor temperature and the flow of C. Everything else follows from the plant.

The whole plant is examples/reactor_mixer/plant.yaml. The steps below go through it one section at a time.

Step 1: the feed flow loops¶

seed: 1
dt: 1s
library: process@1

units:
  - {id: FIC-101, template: flow_loop, sp: 8, K: 0.2, Kp: 2, span: [0, 20], tuning: simc}              # A
  - {id: FIC-102, template: flow_loop, mode: CAS, sp: 2, K: 0.05, Kp: 0.5, span: [0, 5], tuning: simc} # B
  - {id: FIC-103, template: flow_loop, sp: 5, K: 0.1, Kp: 1, span: [0, 10], tuning: simc}              # C

Each flow_loop is a flow controller, a valve with a first-order-plus-dead-time flow response, and a flowmeter. The unit id gives the ISA tags: FIC-101 measures its flow as FT-101 and moves the valve FV-101. K is the flow per % of valve opening, chosen so that each valve works near mid-travel, and tuning: simc computes the controller gain and integral time from that valve model. Kp is how much an upstream pressure change moves the flow; step 5 gives each loop a slowly varying pressure.

The base step dt is 1 s because flow loops are fast, with a 5 s time constant. With dt: 2s the flow control gets visibly noisier, and at 5 s it oscillates.

Step 2: ratio control for B¶

The B flow must stay in proportion to A. A ratio station multiplies the measured A flow by a ratio setpoint, and the B controller follows the result as a remote setpoint:

custom:
  - {id: FFC-102.ratio, op: schedule, value: 0.25, unit: "1"}
  - {id: FFC-102.out, op: product, inputs: {in0: FT-101, in1: FFC-102.ratio}, unit: m3/h}

exogenous:
  - {target: FIC-102.rsp, from: FFC-102.out}

custom adds core operators directly: a schedule holds a value that plans and events can change, and a product multiplies its inputs. Every loop template has an optional rsp port, its remote setpoint. Binding it with from, and starting the loop in mode: CAS as in step 1, makes FIC-102 follow the ratio station. When an operator raises the A flow, the B flow rises with it.

Step 3: the reactor¶

units:
  - {id: TIC-201, template: cstr, sp: 80}

exogenous:
  - {target: TIC-201.flow_a, from: FIC-101.cv}
  - {target: TIC-201.flow_b, from: FIC-102.cv}

cstr is a well-mixed reactor with a temperature loop. Its feeds are the true flows of the two loops (<loop>.cv), not their meter readings. Inside, A and B react on a catalyst:

• The fraction of B left in the reactor follows a mass balance: the feed brings B in, and the outflow and the reaction take it out. It depends on the feed ratio, the residence time V/F, the temperature and the catalyst.
• The reaction rate rises with the temperature (an Arrhenius law) and falls when an impurity poisons the catalyst or when the catalyst ages.
• The reaction heat warms the reactor, and a larger feed cools it. The temperature controller compensates both, so the heating valve TV-201 carries a trace of the reaction.
• The product property q depends on the fraction of B left, the temperature, the impurity and the catalyst's activity. The product leaves after mixing in the reactor and a short outlet line.
• The lab, AT-201, takes a sample every 4 h and reports it 1 h later.

The equations are in the docstring of _cstr_plant in src/homeostat/domains/process/templates.py, and the process library reference lists the template's parameters. None of these states is measured: they appear in run.truth and never in run.observed.

You can look at them right away. The plant alone runs the grade A recipe (8 m³/h of A, B/A = 0.25, 80 °C):

import yaml
import homeostat

plant = yaml.safe_load(open("examples/reactor_mixer/plant.yaml"))
run = homeostat.simulate({**plant, "duration": "10min"})
run.truth[["TIC-201.x_B", "TIC-201.q_product", "TV-201"]].iloc[0]

This takes about 45 s, almost all of it spent settling the plant into its steady state before t = 0. With fresh catalyst and 1 ppm of impurity in A, the fraction of B left is 0.0777, the product property is 1.102, and the heating valve is 65 % open.

Step 4: the mixer and material C¶

units:
  - {id: MX-301, template: mixer, volume: 1, span: [0, 4]}

exogenous:
  - {target: MX-301.flow_a, from: TIC-201.flow}       # the reactor outflow
  - {target: MX-301.prop_a, from: TIC-201.q_product}
  - {target: MX-301.flow_b, from: FIC-103.cv}
  - target: MX-301.prop_b                              # property of material C
    source: {kind: ou, mean: 3.0, std: 0.05, tau: 12h}

The mixer averages the two properties by flow, then the mixture travels through a 1 m³ outlet pipe in plug flow, about 4 minutes at 15 m³/h. An online analyzer, AT-301, takes a sample every 10 minutes and reports it 5 minutes later. Instrument noise is a percentage of the span, so span: [0, 4] keeps the analyzer's noise in proportion to a property near 1.7.

Step 5: disturbances and hidden conditions¶

custom:
  - {id: A.lot_impurity, op: schedule, value: 1.0, unit: ppm}
  - {id: A.impurity, op: sum, weights: [1, 1], inputs: {in0: A.lot_impurity, in1: A.impurity_drift}, unit: ppm}

exogenous:
  - {target: TIC-201.impurity, from: A.impurity}
  - target: A.impurity_drift
    unit: ppm
    source: {kind: ou, mean: 0, std: 0.08, tau: 8h}
  - target: TIC-201.feed_temp
    source: {kind: ou, mean: 25, std: 2, tau: 6h}
  - target: FIC-101.pressure
    source: {kind: ou, mean: 0, std: 0.1, tau: 20min}
  # and the same for FIC-102 and FIC-103

degradation:
  - {target: TIC-201, kind: catalyst_decay, rate: 0.02}

The impurity in feed A is the level of the current lot, which changes when a new lot arrives, plus a slow drift. The catalyst loses about 2 % of its activity per day. The feed temperature and the upstream pressures are slow random disturbances: Ornstein–Uhlenbeck sources, given by their mean, standard deviation and correlation time. Nobody controls or records any of these.

Step 6: the production plan¶

examples/reactor_mixer/plan.yaml moves the plant through three grades in four days and brings in two new lots of A:

duration: 4d

regimes:
  grade_A: {FIC-101.sp: 8, FFC-102.ratio: 0.25, TIC-201.sp: 80, FIC-103.sp: 5}
  grade_B: {FIC-101.sp: 8, FFC-102.ratio: 0.375, TIC-201.sp: 84, FIC-103.sp: 4}
  grade_C: {FIC-101.sp: 10, FFC-102.ratio: 0.15, TIC-201.sp: 78, FIC-103.sp: 6}

plan:
  start: grade_A
  production:
    - {at: 18h, to: grade_B, over: 2h}
    - {at: 42h, to: grade_C, over: 3h}
    - {at: 66h, to: grade_A, over: 2h}

interventions:
  - {at: 30h, target: A.lot_impurity, value: 1.4, id: lot-2}
  - {at: 78h, target: A.lot_impurity, value: 0.7, id: lot-3}

output:
  every: 1min

A regime names the values of one grade, and a transition ramps every value to the next grade over over. The lot changes are ordinary events on the A.lot_impurity schedule; they are recorded in the ground truth but not in the historian. With fresh catalyst and 1 ppm of impurity, the three recipes settle at:

Grade C needs a much lower ratio than a proportional guess suggests. It runs cooler and faster, and both leave more B unreacted: at B/A = 0.2 it would make almost the grade A product (1.056).

Step 7: run it¶

uv run python examples/reactor_mixer/run.py

examples/reactor_mixer/run.py merges the two files, simulates 345,600 steps, writes CSV files to examples/reactor_mixer/output/, and prints the last 12 h of each grade. It takes about five minutes; the first 45 s settle the plant into the grade A steady state.

Simulated 4 days in 284 s.
Wrote 18 historian tags and 86 true signals to examples/reactor_mixer/output.

grade hours  q_true  q_lab  impurity  activity  B_left  heating  mixer_true  mixer_AT301
    A  0-18   1.096  1.128      0.82     0.990  0.0784     63.0       1.699        1.692
    B 20-42   1.491  1.432      1.05     0.969  0.1026     69.3       1.881        1.875
    C 45-66   0.903  0.902      1.25     0.949  0.0653     67.5       1.666        1.673
    A 68-96   0.974  1.006      0.65     0.924  0.0793     65.3       1.674        1.675

q_true is the product property and q_lab the lab's reading of it. impurity is the impurity in the reactor, lower than in feed A because feed B dilutes it. B_left is the fraction of B left unreacted, and heating the heating valve's opening in %. The last two columns are the mixer outlet's true property and its analyzer reading.

Step 8: read the results¶

In Python, the same run gives:

run.observed            # the historian: 18 tags every minute
run.truth               # all 86 true signals, including the hidden ones
run.meta["events"]      # every change, planned or not, with the value before and after
run.meta["timeline"]    # when each grade and transition started and ended

The controllers do their job: every flow and the reactor temperature stay on their setpoints. The product does not stay on plan:

• Grade A, the first time: on plan.
• Grade B: the lot that arrives at 30 h carries 1.4 ppm. The impurity poisons the catalyst, less B reacts, and q ends the grade 6 % above plan without any change of setpoint.
• Grade C: the same lot is still in use, and q stays 5 % above plan.
• Grade A, the second time: a cleaner lot (0.7 ppm) arrives at 78 h, and the catalyst has lost 7.6 % of its activity. The recipe that made 1.096 on the first day now makes 0.974.

The lab follows these shifts hours late and with noise. An operator sees only the lab, the flows, the temperature and the heating valve. The heating valve holds a mixed clue: when less B reacts, less reaction heat is released and the controller opens the valve further to hold the temperature, but the feed temperature moves the valve just as much. Separating the two is exactly the kind of problem a soft sensor has to solve.

What to try next¶

• Build a soft sensor. Predict TIC-201.q_product from run.observed. The true property is in run.truth every minute, while the lab gives about two dozen values.
• Change the hidden conditions. Move the lot changes in plan.yaml, change their impurity, or change the catalyst's decay rate, and see how far the product moves from plan.
• Change the chemistry. The reactor's kinetic parameters (k_ref, E_k, poison, b, gamma, c) are template parameters: set them on the TIC-201 unit in plant.yaml.

Limitations¶

• Speed. The flow loops need a 1 s step, and the engine takes about 0.7 ms per step for this plant, so four days take about five minutes. Settling into the steady state alone takes about 45 s. A compiled step is planned.
• Noise in the B flow. The ratio station passes the A flowmeter's noise into the B setpoint, and the fast B loop amplifies it: the B valve moves by about ±10 %. Real ratio stations usually filter their input.
• The ratio setpoint is not in the historian. FFC-102.ratio is a custom schedule, not a loop setpoint, so the historian does not record it. Its changes are in run.meta["events"].