Physics-Informed Neural Networks (PINN) for ADM1¶
This page explains the PINN estimator from the ground up. It starts with the
idea, then builds up the exact mathematics, explains the ADM1-specific tricks
that make it work at all, and finally how it is implemented in this repository
(pinn.py, pinn_smoother.py).
In one sentence
A PINN is a neural network that reconstructs the digester's hidden state trajectory \(\hat{x}(t)\) by simultaneously fitting the few available sensor measurements and obeying the ADM1 differential equations.
Sections 1–3 are the intuition (no maths), section 4 is the formal loss function, section 5 explains why the stiff ADM1 needs special care, and section 6 maps every idea onto the actual code.
1. The problem this solves¶
A biogas digester is described by 41 internal ADM1 states: dissolved sugars and acids, particulate fractions, microbial biomass, ion/charge species and gas-phase pressures (see the ADM1da model). Online we can only measure a handful of them indirectly: biogas flow \(Q_\text{gas}\), methane \(Q_\text{ch4}\), pH, and a few more. State estimation is the task of reconstructing all 41 states over time from those few noisy signals.
There are two classic routes, each with a weakness:
- Pure mechanistic model (ADM1). Trustworthy structure, but it needs the
exact initial state and exact kinetic parameters. Small errors accumulate and the simulation drifts away from the real plant. - Pure machine learning. Flexible, but it needs a large labelled dataset of
true states to learn from, which on a real plant we never have, because the true states are exactly what we cannot measure.
A PINN combines the two: a neural network flexible enough to fit the sparse data, constrained by the ADM1 physics so that it needs far less data and cannot wander into physically impossible trajectories.
2. A neural network¶
- A neural network is a flexible mathematical function \(y = \text{NN}_\theta(x)\)
with many tunable numbers \(\theta\) (its weights). - Training means adjusting \(\theta\) so that a loss is minimised by gradient
descent. - The input is time \(t\) and the output is the state vector \(\hat{x}(t)\).
So the network is the estimated trajectory itself: ask it at any time \(t\) and it returns the estimated plant state at that instant.
Concretely it is a multilayer perceptron (MLP): a few Linear layers with a
smooth tanh nonlinearity in between. tanh is preferred in PINNs because its
derivatives are smooth.
3. What "physics-informed" means¶
A network trained on data alone would fit the measurement points but be free to do anything between and after them. A PINN adds a second requirement: the trajectory must satisfy the known governing equation
where \(f\) is the ADM1 right-hand side and \(u\) is the (known) substrate feed.
The trick that makes this checkable is automatic differentiation: because the network output \(\hat{x}(t)\) is an explicit function of the input \(t\), we can ask the deep-learning framework for its exact time derivative \(\tfrac{d\hat{x}}{dt}\). ADM1 gives us \(f(\hat{x}, u)\). Their mismatch,
is the physics residual → zero when the trajectory obeys the model.
Training then minimises a sum of three demands:
Gradient descent on \(\theta\) minimises \(L\) and feeds the updated weights back into the network.
The payoff of the physics term: it interpolates through gaps in the data, regularises against sensor noise, and — because the ODE is defined at every instant — lets the model forecast past the last measurement.
4. The mathematics¶
Let the network be \(\hat{x}_\theta(t)\). The estimator minimises
The weights \(\lambda_\text{phys}, \lambda_\text{prior}\) balance the three demands and are the main tuning knobs.
Data loss: fit the measurements. With the differentiable measurement map \(h\) (state \(\to\) sensor channels) and per-channel noise std \(\sigma\), over the observed cells \(\mathcal{O}\) (missing/gated readings are masked out):
Physics loss: obey the ODE. It measures how strongly the predicted trajectory violates the ODE. This is checked at \(N_c\) collocation points \(\{\tau_j\}\), self-chosen times placed across the window. ADM1 mixes very fast processes (acid–base reactions) with very slow ones (biomass). Without a fix the fast ones would dominate the loss, so the residual is divided per state by a typical scale \(s\) so that all count equally (details in section 5):
Prior loss: anchor the start. A weak boundary condition tying the trajectory to a physical prior state \(x_\text{prior}\) at \(t_0\):
Why collocation points matter
Because \(L_\text{data}\) lives only at the measurement times but \(L_\text{phys}\) can be evaluated at any collocation point, the two time ranges are decoupled. Put collocation points past the last measurement and the ODE alone carries the state forward: that is how the same fit forecasts.
5. Why ADM1 is hard¶
A textbook PINN fails on ADM1: the system is stiff (rates spanning many
orders of magnitude) and its states span from \(\sim 10^{-7}\) to tens, all of
which must stay non-negative. Six engineering choices in PinnSmoother make it
converge (it settles stably on a good solution instead of derailing).
These are not optional polish
Measured on the shipped benchmark: without the charge-balance solve of section 5.5 the fit never improves on its own starting point, at any learning rate, loss weighting or network size. The best data loss it ever reached equalled its value at epoch 0 — training only ever made the fit worse. See section 7 for the numbers.
1. Positive, well-scaled outputs¶
So far we said the network outputs \(\hat{x}(t)\). More precisely, its last layer first produces a raw, unbounded vector \(\text{raw}_\theta(t)\). The physical state is then formed from it as a log-deviation from a prior:
The final layer is zero-initialised, so at the start of training
\(\text{raw}_\theta \equiv 0\) and \(\hat{x}(t) = x_\text{prior}\) exactly. The
exponential guarantees \(\hat{x}(t) > 0\) for all time (raw is clamped to
\([-10, 10]\) for numerical safety).
2. Relative physics residual¶
Acid–base reactions in ADM1 run at rates \(\sim 10^{8}\), while biomass changes over days. Normalising the residual by state magnitude would let the fast terms swamp the loss. Instead the scale \(s\) from section 4 is not fixed: each equation is divided by its own current rate \(\lvert f_i \rvert\) (floored so near-equilibrium states with \(f\approx 0\) don't blow up). The result is a relative residual, so no single stiff equation dominates by sheer magnitude.
3. Quasi-steady gas¶
The four gas-phase pressures (especially the total pressure) are numerically
delicate: tiny changes can tip the solution over. Rather than let the network
predict them, it predicts only the 37 liquid states. The gas pressures are
computed from gas–liquid equilibrium at every evaluation (gas_equilibrium_torch).
The four gas ODEs are then satisfied automatically, so the physics loss only has
to check the 37 liquid states.
4. Prior anchoring and per-state weighting¶
The prior term from section 4 pulls all states toward their physical prior at \(t_0\). The data loss pulls the well-observed ones back off it, while the weakly-observed ones stay stuck there. Which states this affects therefore follows automatically from observability. To make the optimiser move such sluggish states instead of parking them at the prior, they can be up-weighted individually in \(L_\text{phys}\): a hand-set weight vector, not a learned parameter. Which states to up-weight is your choice (domain knowledge / diagnosis), not something the training decides.
The prior itself must be a physically reachable state. A componentwise
statistic (say the per-state median over a training set) generally is not: taking
each component's median independently breaks the identities that tie them
together. On ADM1 that identity is the charge balance, and the resulting "average"
state implies a pH far outside the physical range. Use a medoid — an actual
observed state — instead; PinnData does this for you (section 6).
5. Quasi-steady charge balance (solve_cation=True)¶
ADM1 gets pH from electroneutrality:
At a normal operating point fixed is a difference of terms of order
\(0.1\!-\!0.2\) kmol m⁻³ that cancels down to \(\sim\!10^{-6}\) — 5.3 decades of
cancellation. Since \(S_\mathrm{H} \approx K_w/\text{fixed}\), we effectively have
\(\text{pH} \approx -\log_{10} K_w + \log_{10}(\text{fixed})\): perturb any one ion
state by 1 % and fixed moves by orders of magnitude.
That is the same knife-edge as the biogas map in section 5.3, and it has the same consequence — pH responds with \(\sim\!5\cdot10^{4}\,\sigma\) to a 1 % state change, so its gradient dictates every optimiser step and the fit diverges.
The fix mirrors the gas solve: do not let the network guess the quantity the cancellation hangs on. The network predicts pH directly (a well-scaled \(\approx 7.5\) output whose sensitivity is 1:1) and \(S_\text{cat}\) is solved from the charge balance in closed form:
What this trades away
\(S_\text{cat}\) stops being a differential state and becomes an algebraic one
— the model turns from an ODE into a DAE, its trivial dilution equation
replaced by the electroneutrality constraint. PinnSmoother therefore masks
it out of both the physics residual and the prior anchor. This is deliberate:
\(S_\text{cat}\) is a bookkeeping charge rather than a measured species, and
enforcing its dilution equation is exactly what creates the ill-conditioning.
For the same reason the solved value may come out slightly negative — that is
a net anion excess, well within how ADM1 uses these two slots.
6. Bounded residuals (res_clip)¶
A squared residual has a gradient proportional to the residual, so one badly
conditioned channel dictates the whole step. res_clip switches the residual to
Huber: quadratic below the threshold \(\delta\), linear above it, so the
gradient is capped at \(2\delta\) while staying non-zero.
This is deliberately softer than a hard clamp (which zeroes the gradient beyond the threshold and leaves a channel that has already blown past it with no force pulling it back). It bounds the damage from an ill-conditioned channel but does not substitute for section 5.5 — it caps how far a bad step can go, not the conditioning that causes it.
6. How it is implemented here¶
The implementation has two layers: a textbook template and the production estimator built on top of it.
| Layer | File | Role |
|---|---|---|
ADM1PINN + PINNLoss |
deep_learning/pinn.py |
The building blocks: the MLP \(t \to x\) (tanh, configurable hidden layers, optional dropout) and a generic data + physics loss template. Readable reference implementation. |
PinnSmoother |
deep_learning/pinn_smoother.py |
The production estimator. Builds on ADM1PINN and adds the log-transform, the three-term loss, rate scaling, quasi-steady gas, forecasting, rolling updates and MC-Dropout uncertainty. |
The differentiable physics pieces it wires together all come from the base
pyadm1 package:
| Symbol | Provided by |
|---|---|
| \(f\) — ADM1 right-hand side | pyadm1.core.adm1_torch.adm1da_rhs_torch |
| gas–liquid equilibrium | pyadm1.core.adm1_torch.gas_equilibrium_torch |
| \(h\) — measurement map | deep_learning.observation_torch.TorchObservationModel |
| charge-balance inversion | deep_learning.charge_balance.solve_cation_for_ph |
| feed / parameters | Adm1TorchParams — the network input stays time only; a feed that changes inside the window enters through params_at |
Feeding it: the dataset adapter¶
PinnData (deep_learning/data_adapter.py) turns any dataset registered in
filter_tuning.datasets into what both PINN variants need. Three guarantees
matter for correctness:
- The same split as the filters. Both go through
EstimatorDataset.split_indices— stratified by the series label, grouped by series. Same(val_frac, seed)⇒ same series, which is what makes a network and a filter comparable. Freeze it withsave_split()/split_file=. - Train-only statistics. Feature normalisation,
x_ref,x_prior(as a
medoid) and the per-state scale are estimated on the training series alone. - Feed-matched physics.
smoother_inputs()attaches aparams_atclosure
over each series' own feed, so the ODE residual is evaluated against the plant that produced the data. This matters more than it sounds: an unfed parameter snapshot (q_ad = 0) models a closed batch reactor, and the influent term \(D_\text{in}\,s_\text{in}\) dominates the ADM1 right-hand side.
from pyadm1ode_estimation.estimation.deep_learning import PinnData
data = PinnData.build("benchmark", val_frac=0.2, seed=0)
inputs = data.smoother_inputs("val", days=5.0) # one payload per series
Batch, not recursive¶
Unlike the UKF (a recursive StateEstimator doing predict/update each step), the
PINN is a BatchEstimator: fit the whole window once, then query it.
from pyadm1ode_estimation.estimation.deep_learning import PinnSmoother
it = inputs[0]
smoother = PinnSmoother(
data.physics_params(), data.obs_model(), it.x_prior, it.x_scale,
quasi_steady_gas=True, # section 5.3 — slave the gas pressures
solve_cation=True, # section 5.5 — slave S_cation to a predicted pH
res_clip=3.0, # section 5.6 — Huber threshold in sigma
params_at=it.params_at, # the series' own, time-varying feed
)
smoother.fit(**it.fit_kwargs(), epochs=2000, lr=1e-3)
traj = smoother.estimate(it.obs_times) # (T, 41) states + std
estimate returns a TrajectoryEstimate (time, x_hat, std), the shared
output currency with the UKF, so the twin-experiment harness scores both
estimator families the same way.
restore_best=True (the default) returns the weights of the best-scoring epoch
rather than the last ones. The collocation fit is not monotone — it routinely
reaches its best trajectory and then walks away from it, so returning the final
weights throws away the answer it already found. It also makes the fit
monotone-safe: it can never return something worse than the prior.
Forecasting¶
Because the collocation window [t0, t1] is independent of the measurement
times, setting t1 past the last measurement makes the ADM1 ODE carry the
state into the data-free tail — a physics-driven forecast from the same fit.
Online operation¶
update(...) warm-starts from the current weights and optimiser state for
cheap incremental re-fits as new samples arrive:
- growing window: keep the whole history, anchor stays at the original
t0; or - sliding window: a fixed horizon
t0 = t1 − window, self-anchored at the
network's own current estimate at the window start (bounds compute on long runs).
Uncertainty¶
Optional MC-Dropout: set dropout > 0 and call estimate(..., mc_samples=k)
to get the mean and standard deviation over k stochastic forward passes.
7. Measured status¶
Numbers below are on the shipped benchmark's validation split (4 operating
modes), 5-day windows, reproducible via experiments/pinn_gate/.
The optimisation is solved. Data loss = mean squared standardised residual on the measured channels; the noise floor (what the true state scores) is 0.57.
| Configuration | Median best data loss |
|---|---|
| baseline | 3520 |
+ res_clip |
11.7 |
+ solve_cation (section 5.5) |
0.562 |
The measurement fit now reaches its theoretical optimum. Before the charge-balance solve, training never improved the fit at all — the best value it ever reached was its value at epoch 0.
The state accuracy is not. Median-over-41-states NRMSE against the trivial baseline of holding the prior constant:
| median NRMSE | |
|---|---|
hold x_prior (do nothing) |
25.1 % |
PINN A (best rate_floor) |
22.0 % — but only 2 of 4 modes |
| hold the true \(x(t_0)\) (reference) | 1.2 – 3.7 % |
A rate_floor sweep shows no robust optimum. Five sensors do not pin down 41
states, and the ADM1 term cannot close the gap because a per-window fit has no
way to learn the true kinetics (the benchmark perturbs them lognormally per
series, \(\sigma = 0.25\)). The remaining gap is observability, not architecture
— which is why network width, depth and activation have not been tuned: they are
demonstrably not the binding constraint.
Limits. MC-Dropout uncertainty is not calibrated. A from-scratch fit costs ~2.5 min per 5-day window, which caps how much hyperparameter search is affordable.
Where it sits. The UKF (see UKF in practice and SR-UKF performance) is recursive and cheap per step. Note that on this metric the UKF is itself worse than doing nothing in 3 of 4 modes (50.5 % vs. 34.2 % overall on the test split) — so report both baselines, or a result looks better than it is. The pre-trained observer clears the same bar in all four modes, because it learns across series what a single-window fit cannot. A covariance-intersection hybrid can fuse estimators.
Source files¶
pyadm1ode_estimation/estimation/deep_learning/pinn.py—ADM1PINN,PINNLosspyadm1ode_estimation/estimation/deep_learning/pinn_smoother.py—PinnSmootherpyadm1ode_estimation/estimation/deep_learning/charge_balance.py—solve_cation_for_ph,apply_phpyadm1ode_estimation/estimation/deep_learning/data_adapter.py—PinnData,FeatureSpecpyadm1ode_estimation/estimation/deep_learning/observation_torch.py—TorchObservationModelexperiments/pinn_gate/— the gate measurement and therate_floorsweep
References¶
- Raissi, M., Perdikaris, P. & Karniadakis, G. E. (2019). Physics-informed
neural networks. Journal of Computational Physics 378:686–707. - ADM1da model and state indices: ADM1da model.
API reference¶
pyadm1ode_estimation.estimation.deep_learning.pinn_smoother.PinnSmoother(params, obs, x_prior, x_scale=None, hidden_layers=(64, 64, 64), dropout=0.0, lambda_phys=1.0, lambda_prior=1.0, physics_scaling='rate', rate_floor=1.0, physics_weight=None, res_clip=None, restore_best=True, quasi_steady_gas=False, solve_cation=False, gas_n_iter=20, params_at=None, device='cpu', seed=0, dtype=torch.float64)
¶
Offline PINN state-trajectory smoother; see module docstring.
Methods:¶
estimate(times, mc_samples=1)
¶
Query the fitted trajectory (+ MC-Dropout uncertainty).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
times
|
Sequence[float]
|
Query times, shape |
required |
mc_samples
|
int
|
Number of stochastic forward passes for the MC-Dropout
uncertainty band. |
1
|
fit(obs_times, obs_values, t0, t1, n_collocation=200, epochs=2000, lr=0.001, grad_clip=1.0, verbose=False)
¶
Train the smoother over [t0, t1].
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
obs_times
|
Sequence[float]
|
Measurement times, shape |
required |
obs_values
|
ndarray | dict[str, Sequence[float]]
|
Either an array |
required |
t0, t1
|
Collocation window bounds [days] — the ODE residual is
enforced across all of |
required | |
n_collocation
|
int
|
Number of uniform physics-residual points in the window. |
200
|
epochs, lr
|
Adam training schedule. |
required |
Forecasting (physics-driven, still a per-window fit — no pre-training):
the data loss lives only at obs_times, so if t1 runs past
the last measurement the tail (max(obs_times), t1] has no data
and the ADM1 ODE alone (dx/dt = f(x)) carries the state forward
from the data-anchored end — an initial-value continuation. Pass the
historical measurements as obs_times and set t1 to the forecast
horizon; estimate(t) for t beyond the data is then a forecast
(put enough n_collocation points to cover the tail densely).
Returns:
| Type | Description |
|---|---|
dict[str, list[float]]
|
Loss history dict with keys |
update(obs_times, obs_values, t1, *, window=None, n_collocation=200, epochs=100, lr=0.0003, grad_clip=1.0, verbose=False)
¶
Warm-started incremental re-fit for rolling operation.
Reuses the current network weights and optimiser state instead of
retraining from scratch, so each new step (e.g. a new hourly sample)
converges in far fewer epochs. Requires a prior :meth:fit.
Window mode is chosen by window:
window is None— growing window:t0stays at the original :meth:fitt0(all measurements from the start are used and the original prior anchorsx̂(t0)). Use when you keep the whole history.windowa float — sliding window of that many days:t0 = t1 − window, only measurements in[t0, t1]are used, and the window start is anchored to the network's own current estimate att0(a self-consistent hand-off), not the stale initial prior. Use to cap compute / memory on a long run.
As in :meth:fit, t1 may run past the last measurement to forecast.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
obs_times, obs_values
|
measurements (the full history for a growing window; anything is fine, the window is applied here). |
required | |
t1
|
float
|
current window end [days] (also the forecast horizon). |
required |
window
|
float | None
|
sliding-window length [days], or |
None
|
epochs, lr
|
warm-start schedule (few epochs, smaller lr than fit). |
required |
Returns:
| Type | Description |
|---|---|
dict[str, list[float]]
|
Loss history dict ( |
pyadm1ode_estimation.estimation.deep_learning.pinn.ADM1PINN(input_dim=2, output_dim=41, hidden_layers=(64, 64, 64), dropout=0.0)
¶
Bases: Module
Physics-Informed Neural Network (PINN) for the ADM1 (AP 4.3).
Estimates the state vector x(t) from time t and the time-dependent input u(t) (substrate feed).
The output covers the full extended ADM1da state (41 slots: 12 soluble, 10 particulate PS/PF fractions, 7 biomass, 8 charge/ion states, 4 gas-phase slots). The 4 gas-phase states (p_gas_, p_total) are not optional: they couple back into diff_S_h2/ch4/co2 via the Rho_T_ transfer rates, so the ADM_ODE physics residual cannot be closed without them.
Initialize the PINN.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
input_dim
|
int
|
Input dimension (time t + u_dim), where u_dim is the number of substrate feed channels, i.e. input_dim = 1 + u_dim. |
2
|
output_dim
|
int
|
Dimension of the extended ADM1da state vector (41). |
41
|
hidden_layers
|
Sequence[int]
|
Number of neurons per hidden layer. |
(64, 64, 64)
|
dropout
|
float
|
Dropout probability after each hidden layer; |
0.0
|