Read this post on tubhyam.dev for interactive components, animations, and the best reading experience.

A chemical reactor is a dynamical system. Raw materials flow in, chemical reactions transform them, and products flow out. The state of the reactor — temperatures, concentrations, pressures — evolves continuously in time according to coupled differential equations. If you know the rate laws, the heat transfer coefficients, and the flow patterns, you can write down these equations and solve them numerically.

The problem is that you rarely know all the parameters. Rate constants depend on catalyst condition. Heat transfer coefficients change as fouling accumulates on surfaces. Flow patterns shift when internals degrade. A first-principles model calibrated on commissioning data drifts out of accuracy within months.

Neural ODEs offer a way to learn the dynamics directly from operating data — and to do it in a physically consistent way.

The Classical Approach

A continuous stirred-tank reactor (CSTR) with a single exothermic reaction obeys:

$$\frac{dC_A}{dt} = \frac{F}{V}(C_{A0} - C_A) - k_0 e^{-E_a/RT} C_A^n$$

$$\frac{dT}{dt} = \frac{F}{V}(T_0 - T) + \frac{(-\Delta H_r)}{\rho C_p} k_0 e^{-E_a/RT} C_A^n - \frac{UA}{\rho C_p V}(T - T_c)$$

These are elegant equations with clear physical meaning. Each term corresponds to a specific physical process: flow, reaction, heat generation, heat removal. The parameters — $k_0$, $E_a$, $n$, $U$, $\Delta H_r$ — have units and physical interpretations.

Let me unpack these equations term by term, because understanding the structure is essential for understanding where the neural network enters.

The mass balance (first equation) says: the rate of change of reactant concentration equals the rate of fresh feed entering ($F/V \cdot C_{A0}$) minus the rate of material leaving ($F/V \cdot C_A$) minus the rate consumed by reaction ($k_0 e^{-E_a/RT} C_A^n$). The Arrhenius term $k_0 e^{-E_a/RT}$ is the temperature-dependent rate constant — it encodes the exponential sensitivity of reaction rate to temperature. A 10 K increase in temperature can double the reaction rate.

The energy balance (second equation) has three terms: enthalpy of the feed stream ($F/V \cdot (T_0 - T)$), heat generated by the exothermic reaction ($(-\Delta H_r)/(\rho C_p) \cdot r$), and heat removed by the cooling jacket ($UA/(\rho C_p V) \cdot (T - T_c)$). The coupling between these equations creates the reactor's characteristic nonlinear dynamics — temperature affects reaction rate, which generates heat, which changes temperature.

The model works perfectly on day one and degrades steadily. Recalibration requires taking the reactor offline, running test batches, and re-fitting parameters — a process that costs tens of thousands of dollars in lost production. You need continuous adaptation.

Neural ODEs: The Idea

A neural ODE replaces (or augments) the right-hand side of the differential equation with a neural network:

$$\frac{d\mathbf{x}}{dt} = f_\theta(\mathbf{x}, t, \mathbf{u})$$

where $\mathbf{x}$ is the state vector, $\mathbf{u}$ is the input (feed conditions, coolant temperature), and $f_\theta$ is a neural network parameterized by $\theta$. The key insight from Chen et al. (2018): you can train this by backpropagating through the ODE solver using the adjoint method.

$$\frac{d\mathbf{a}}{dt} = -\mathbf{a}^T \frac{\partial f_\theta}{\partial \mathbf{x}} \quad \text{(adjoint equation)}$$

The adjoint equation runs backward in time, computing gradients without storing intermediate states. Memory cost is $O(1)$ in the number of solver steps — you can integrate over arbitrarily long time horizons without running out of GPU memory.

Why This Matters for Reactors

Chemical reactor dynamics unfold over hours to days. A batch reactor might run for 8 hours with measurements every 30 seconds — that's 960 time steps. Standard backpropagation through an RNN would require storing all 960 hidden states. The adjoint method makes this tractable: solve forward, solve the adjoint backward, compute gradients, update parameters.

Pure Neural vs. Physics-Informed

The naive approach is to replace the entire right-hand side with a neural network. This works for interpolation — predicting behavior within the training distribution — but fails catastrophically for extrapolation:

The physics-informed version is slightly worse at interpolation (the constraints limit expressiveness) but 17x better at extrapolation. More importantly, it never produces physically impossible states — temperatures don't diverge, concentrations don't go negative, energy is conserved.

The Hybrid Architecture

The architecture that works is not "neural network replaces physics" but neural network corrects physics. Keep the known structure — mass balances, energy balances, thermodynamic constraints — and use the neural network only for the uncertain terms:

The neural network $r_\theta$ learns the effective reaction rate — absorbing catalyst deactivation, mixing imperfections, and side reactions that the first-principles model doesn't capture. The network $q_\theta$ learns the effective heat transfer — absorbing fouling, flow maldistribution, and ambient losses.

Explore the CSTR dynamics yourself — adjust the coolant temperature, dilution rate, and activation energy to see how operating conditions affect the phase portrait. Toggle "neural ODE" to see how the learned correction shifts the trajectory:

Key Design Principle

The neural network should learn the residual between the first-principles model and reality — not the entire dynamics. This residual is typically much simpler than the full dynamics: it's smooth, low-dimensional, and bounded. A small network (2-3 layers, 64-128 hidden units) suffices, reducing overfitting risk and training time.

The Full Mathematical Specification

To make the hybrid approach concrete, here are the complete CSTR balance equations with the neural correction terms written out explicitly. A real reactor model involves more terms than the simplified version above.

Full mass balance for species $A$ in a CSTR with volume $V$, volumetric flow rate $F$, and $N_r$ reactions:

$$\frac{dC_A}{dt} = \frac{F}{V}(C_{A0} - C_A) - \sum_{j=1}^{N_r} \nu_{A,j} \cdot k_{0,j} \exp!\left(\frac{-E_{a,j}}{RT}\right) \prod_i C_i^{n_{i,j}} + \underbrace{\delta r_\theta(C_A, C_B, T, t)}_{\text{neural correction}}$$

The neural correction $\delta r_\theta$ accounts for everything the mechanistic model misses: catalyst deactivation (which reduces the effective $k_0$ over weeks), competitive adsorption on the catalyst surface, mass transfer limitations from poor mixing, and side reactions that produce byproducts not tracked in the simplified model.

Full energy balance including the cooling jacket dynamics:

$$\frac{dT}{dt} = \frac{F}{V}(T_0 - T) + \sum_{j=1}^{N_r} \frac{(-\Delta H_{r,j})}{\rho C_p} r_j - \frac{UA}{\rho C_p V}(T - T_c) + \underbrace{\delta q_\theta(T, T_c, C_A, F)}_{\text{neural correction}}$$

$$\frac{dT_c}{dt} = \frac{F_c}{V_c}(T_{c,\text{in}} - T_c) + \frac{UA}{\rho_c C_{p,c} V_c}(T - T_c)$$

The cooling jacket equation is typically left purely mechanistic — the heat transfer coefficient $U$ is the parameter that drifts most as fouling builds up, and the neural correction $\delta q_\theta$ absorbs this drift. The neural term learns to output an additive correction to the heat transfer rate, effectively replacing the uncertain $UA$ product with a data-driven estimate that adapts online.

The key constraint: $\delta r_\theta$ and $\delta q_\theta$ are initialized near zero (via small weight initialization or a final tanh with a small gain). At deployment time, the physics model alone provides a reasonable baseline; the neural corrections only activate when the data shows systematic deviations from the first-principles prediction.

Training Pipeline

Taking a hybrid neural ODE from research prototype to deployed digital twin requires a careful multi-stage pipeline. Each stage has a specific purpose — skipping stages leads to models that work in notebooks but fail in production.

Stage 1: Simulate. Generate synthetic training data from the first-principles CSTR model across a wide range of operating conditions — varying feed concentration (0.5-5.0 mol/L), feed temperature (280-340 K), coolant temperature (280-320 K), and flow rate (0.5-5.0 L/min). Include simulated disturbances: step changes, ramp changes, and sinusoidal perturbations. This gives the neural network a broad base of physically consistent trajectories to learn from.

Stage 2: Pretrain. Train the hybrid model on synthetic data with the neural correction terms active. The physics branch provides the correct structure; the neural branch learns to produce near-zero corrections (since the synthetic data came from the same physics model). This stage teaches the neural network the format of corrections — what scale they operate on, how they should interact with the physics terms — without yet learning the real-world residual.

Stage 3: Fine-tune on plant data. Transfer the pretrained model and fine-tune on real operating data (typically 2-4 weeks of continuous operation). The physics parameters can optionally be re-estimated during this stage. The neural corrections now learn the actual plant-model mismatch: catalyst aging, fouling, unmeasured disturbances. Learning rate is reduced 10x from pretraining.

Stage 4: Deploy. Export the model for edge inference. The ODE solver runs on the plant's control system (typically a PLC or DCS with a Python runtime). Inference must complete within the control cycle time — usually 1-5 seconds for a CSTR.

Stage 5: Online adaptation. Continuously update the neural branch as new plant data arrives. This is the stage that keeps the model accurate over months without manual recalibration.

Training: What the Loss Function Needs

The loss function for a reactor digital twin has three components:

$$\mathcal{L} = \underbrace{\mathcal{L}{\text{data}}}{\text{match measurements}} + \lambda_1 \underbrace{\mathcal{L}{\text{physics}}}{\text{conservation laws}} + \lambda_2 \underbrace{\mathcal{L}{\text{stability}}}{\text{bounded dynamics}}$$

Data loss is standard: MSE between predicted and measured state trajectories. But measurements are sparse — temperature every 30 seconds, composition every 15 minutes (requires lab analysis). The ODE solver fills in between measurements, but the gradient signal is weak when observations are far apart.

Physics loss enforces conservation laws that hold regardless of operating conditions:

$$\mathcal{L}{\text{physics}} = \left| \frac{d(\text{total mass})}{dt} - (\dot{m}{\text{in}} - \dot{m}_{\text{out}}) \right|^2$$

Stability loss prevents the neural network from learning dynamics that diverge:

$$\mathcal{L}{\text{stability}} = \max\left(0, \frac{\partial f\theta}{\partial \mathbf{x}} - \epsilon\right)$$

This is a soft constraint on the Jacobian eigenvalues — the learned dynamics must be locally dissipative. Without this, the neural ODE can find solutions that fit the training data but blow up under slightly different initial conditions.

Why Reactors Are Not Weather

Weather forecasting with neural ODEs (Pangu-Weather, GraphCast) has attracted enormous attention. Chemical reactor modeling seems similar — both are dynamical systems governed by PDEs. But the operating regime is fundamentally different:

Reactors have active control inputs that change the dynamics. The operator adjusts feed rates, setpoints, and coolant flow in response to the model's predictions. This creates a feedback loop: the model influences the data it will be trained on. Weather forecasting has no such loop — the atmosphere doesn't respond to forecasts.

Reactors also have hard safety constraints. A temperature prediction that's wrong by 20 K might trigger a thermal runaway. The model must be conservative — it's better to predict slightly worse performance than to miss a dangerous excursion. This asymmetric loss structure doesn't exist in weather.

Let me make these differences concrete with specific examples that distinguish reactor modeling from weather modeling:

Coolant control loop (manipulated variable). A CSTR has a cooling jacket with a controllable coolant flow rate $F_c$. The plant's PID controller adjusts $F_c$ based on the temperature setpoint and the model's predicted trajectory. If the model predicts a temperature rise, the controller increases $F_c$ preemptively. This means the model's output directly affects the system's future state — a closed-loop coupling that weather models never face. The atmosphere does not adjust wind patterns because a forecast predicted rain.

Safety constraint ($T < T_{\text{runaway}}$). Exothermic reactors have a critical temperature above which the reaction rate increases faster than heat removal — this is thermal runaway. For a typical organic oxidation, the runaway temperature might be $T_{\text{runaway}} = 420$ K. The model must satisfy $T(t) < T_{\text{runaway}}$ for all $t$ — not just on average, not just 95% of the time, but always. A weather model that predicts 30°C when the actual temperature is 32°C is slightly wrong. A reactor model that predicts 410 K when the actual temperature is 425 K has missed a potentially catastrophic excursion. This asymmetry demands conservative prediction — overestimating temperature is always safer than underestimating it.

Feedback from product quality. In a real plant, downstream quality measurements (e.g., product purity from a gas chromatograph, measured every 15-30 minutes) feed back to the reactor's operating conditions. If purity drops, the operator may increase residence time (decrease $F$), raise reactor temperature, or increase catalyst regeneration frequency. Each of these actions changes the reactor dynamics. The model must predict not just the current trajectory but how the trajectory will change in response to corrective actions that haven't happened yet. Weather prediction is purely open-loop — there is no atmospheric control room adjusting the jet stream.

Unmeasured disturbances. Reactor feed composition varies with upstream processes — the concentration of impurities might shift by 5-10% over a day. Ambient temperature changes affect heat loss through the reactor walls. Catalyst activity declines continuously as active sites are poisoned or sintered. None of these disturbances are directly measured; they manifest as slow drift in the model residuals. The neural correction terms $r_\theta$ and $q_\theta$ must learn to infer these hidden disturbances from their effect on the measured states.

Online Adaptation

The real value of a neural ODE reactor model is continuous online learning. As new operating data arrives, the neural network parameters update to track slow drifts in reactor behavior:

The neural ODE maintains <1.5 K accuracy over 6 months with no manual recalibration. The classical model degrades to 11.7 K — nearly 10x worse. The key is that only the neural network parameters ($r_\theta$, $q_\theta$) are updated online; the physics structure (mass balances, energy balances) is fixed. This prevents the model from learning spurious dynamics during transients.

The Adaptation Algorithm in Detail

The online update follows a specific protocol designed for safety and stability:

Gradient update rule. Every 100 new measurements (approximately every 50 minutes at 30-second sampling), the neural branch parameters are updated:

$$\theta_{t+1} = \theta_t - \eta \nabla_\theta \mathcal{L}_{\text{physics+data}}(\theta_t)$$

The loss $\mathcal{L}_{\text{physics+data}}$ is computed over a sliding window of the most recent 200 measurements (approximately 100 minutes of operation). This window is long enough to capture slow dynamics but short enough to respond to step changes in operating conditions.

Frozen vs. adaptive parameters. The physics branch parameters ($F$, $V$, $C_{A0}$, $T_0$) are treated as known constants from the plant's instrumentation — they do not update. The neural branch parameters ($r_\theta$, $q_\theta$ — approximately 8,000 weights for a 2-layer, 64-hidden-unit MLP per correction term) are the only parameters that adapt. This separation is critical: if the physics parameters were also adaptive, the model could "explain" a temperature drift by changing the reactor volume $V$ instead of learning the actual cause (fouling on the heat transfer surface). The physics provides structural constraints that prevent the neural network from learning physically nonsensical explanations.

Learning rate schedule. The adaptation learning rate follows cosine decay from an initial value of $\eta_0 = 10^{-4}$ to a minimum of $\eta_{\min} = 10^{-6}$ over a 30-day cycle, then resets:

$$\eta_t = \eta_{\min} + \frac{1}{2}(\eta_0 - \eta_{\min})\left(1 + \cos\left(\frac{\pi \cdot (t \bmod T)}{T}\right)\right)$$

where $T$ is the cycle length (30 days of measurements). The cosine schedule allows larger updates at the start of each cycle (when drift may have accumulated) and smaller updates at the end (when the model should be tracking smoothly). The 30-day cycle aligns with typical catalyst deactivation timescales.

Safety guard. Before applying any parameter update, the system runs a forward simulation over the next 30 minutes using the proposed new parameters. If the predicted trajectory violates the safety constraint ($T > T_{\text{runaway}} - 10$ K margin), the update is rejected and the previous parameters are retained. This rollback mechanism ensures that the online learning never produces a model that could miss a thermal runaway event.

Practical Impact

A reactor digital twin that stays accurate for 6+ months without manual recalibration eliminates the most expensive part of model maintenance in process industries. Recalibration campaigns typically cost $50-200K in lost production and engineering time. Continuous neural adaptation replaces this with a compute cost of ~$10/month on a single GPU.

What's Next

The ReactorTwin project implements this architecture with PyTorch's torchdiffeq solver. Current work focuses on multi-reactor networks — cascaded CSTRs and tubular reactors where the outlet of one unit feeds the inlet of the next. The challenge is that neural corrections in upstream reactors propagate through the entire network, requiring careful gradient management to avoid instability.

Originally published at tubhyam.dev/blog/neural-odes-for-reactor-modeling