06-763 / L13

Lecture 13: Scientific machine learning: PINNs, neural ODEs and neural DAEs

Week 7, Machine learning and deep learning

Systems and Toolchains for AI Engineers

Systems and Toolchains for AI Engineers
06-763 / L13

What today is about

  1. Why physics: a bioreactor model that fits, then predicts the impossible
  2. Gen 0 to Gen 2: from surrogates to scientific machine learning
  3. PINNs: physics in the loss
  4. From RNNs to neural ODEs: learning the vector field
  5. Physics-enforced machine learning: hard constraints, and neural DAEs
  6. Sequential or simultaneous: two ways to train a dynamic model
  7. Training neural DAEs: SiNDAE, and monoclonal antibodies
  8. Constraints inside the network: projection layers
  9. Not only neural networks
Systems and Toolchains for AI Engineers
06-763 / L13

Why physics

Systems and Toolchains for AI Engineers
06-763 / L13

Why physics, a fed-batch bioreactor

  • Fed-batch: start with cells and medium, feed substrate during the run, harvest at the end
  • Nothing leaves before harvest, so the volume grows
  • The feed sets how much food the cells see, so it steers growth and product
  • Used for: monoclonal antibodies (mAbs) and other therapeutic proteins, recombinant proteins in E. coli, baker's yeast, penicillin (summary)
  • Why it matters: fed-batch is still "the most commonly utilized process type for biomanufacturing" (Bioprocess Biosyst. Eng., 2024)
  • CHO (Chinese hamster ovary) cells make about 80% of U.S.-licensed mAb products (mAbs, 2025)
Systems and Toolchains for AI Engineers
06-763 / L13

Why physics, the balances term by term

  • Biomass

  • Product

  • Substrate

  • Volume

: concentrations (g/L). : volume (L). : feed rate (L/h). : substrate in the feed (g/L). , : yields (g per g). : specific growth rate (1/h), the kinetics. Not known: we learn it.

Systems and Toolchains for AI Engineers
06-763 / L13

Why physics, why learn the kinetics?

  • Balances: conservation of mass, certain
  • Rates: uncertain. Even a simplified CHO network has 30 reaction rates, to , each with its own law and parameters
  • in our reactor is one of these rates: learn it from data

Simplified CHO metabolic network, Fig. 5 of Wang, Harcum and Xie (2025), CC BY 4.0.

Systems and Toolchains for AI Engineers
06-763 / L13

Why physics, inference on unseen data

  • Three models, all trained on the same three batches (0 to 40 h), predict a new batch for 60 h
  • Purely data-driven, no physics: substrate down to −0.43 g/L, biomass off by 0.43 g/L on average
  • Physics-informed, the balances with a learned : closer, but still −0.37 g/L
  • Physics-enforced, a neural DAE: holds, lowest value 0.05 g/L
Systems and Toolchains for AI Engineers
06-763 / L13

Machine learning in science and engineering, Gen 0 and Gen 1

Gen 0: surrogate models

  • "Often, the codes are computationally expensive to run, and a common objective of an experiment is to fit a cheaper predictor of the output to the data" (Sacks et al., 1989)
  • Goal: a cheap stand-in for an expensive model

Gen 1: machine learning with intention

  • The surrogate also says how uncertain it is, and that decides where to sample next (Jones, Schonlau and Welch, 1998)
  • This is the idea behind Bayesian optimization: a cheap surrogate with its uncertainty chooses the next expensive experiment or simulation

Neither generation knows any physics: everything comes from the data.

Systems and Toolchains for AI Engineers
06-763 / L13

Machine learning in science and engineering, Gen 2

First-principles
Data-driven (ML)
Gen 0, Gen 1
Minimal data
Extensive data
No ML
No physics
Challenges
unknown physics
time to impact
Challenges
no gained insight
cost to impact

Moderate data
Physics + ML
Hybrid
Gen 2

Gen 2: guaranteed constraints and interpretability. Missing: tools to train, optimize, quantify uncertainty.

Schematic: V. Alves, after Fig. 1 of Shah et al. (2025), Comput. Chem. Eng.

Systems and Toolchains for AI Engineers
06-763 / L13

Scientific machine learning

Scientific machine learning (SciML) combines mechanistic, first-principles models with machine learning: the network learns only what the physics does not already say.

Pendulum: unknown friction

Heat conduction: unknown source

CSTR: unknown kinetics

Tank: unknown valve law

Black: conservation laws, known. Red: a closure term, a job for scientific machine learning.

Systems and Toolchains for AI Engineers
06-763 / L13

Scientific machine learning, four ways to add the physics

PINN

Physics in the loss: a penalty, small, not zero.

Neural ODE

Network inside the equations. Bounds not enforced.

Neural DAE

One optimization problem: constraints enforced.

Projection layer

Constraints in the last layer, one output at a time.

Same ingredients each time: data, a network, and the physics in green. Where the green sits decides what it guarantees.

Systems and Toolchains for AI Engineers
06-763 / L13

PINNs

Systems and Toolchains for AI Engineers
06-763 / L13

Physics-informed neural networks (PINNs), introductory example

A spring-mass with a damper that can be dragged and released

PINN (physics-informed neural network): a network trained to fit the data and, at the same time, to make the residual of the governing equation small (Raissi, Perdikaris and Karniadakis, 2019).

  • Drag the green mass, let go
  • Our case: kg, N·s/m, N/m, from m
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, the equation of motion

: displacement from rest (m), up is positive. : velocity (m/s). : acceleration (m/s²). : mass (kg). : spring stiffness (N/m). : damping coefficient (N·s/m).

  • Newton's second law:

  • Spring (Hooke's law) pulls back toward rest:

  • Damper (viscous friction) opposes the motion:

  • Move both forces to the left:

  • Released from rest at 1 m: m,

Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, the task

  • We measure at 10 times, all in the first 0.36 s: about one oscillation
  • We want for the whole second, including to s, where there are no measurements
  • We know the equation of motion, and , ,
  • Question: can a network use that knowledge where it has no data?
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, a plain network first

  • A network : input the time, output the displacement

  • Trained on the data alone, by mean squared error over the points:

  • It fits the 10 points

  • After the last one it has nothing to go on

Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, physics in the loss

→

Physics + ML

→

data loss

physics loss

Physics is favored, not enforced: a penalty, made small but not zero

Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, collocation points

Collocation point: a time at which we ask the equation to hold. It needs no measurement, so we can place as many as we like, anywhere.

Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, the loss term by term

L(θ) = 1N Σ (xNN(ti) − xi)2 + λ 1M Σ (m x″NN + μ x′NN + k xNN)2
Fit the data: the 10 measurements, as before
Obey the equation: its residual at the M = 40 collocation points
λ = 0.0001 weighs physics against data: the residual is in newtons, hundreds of times larger than x
Derivatives x″ and x′ of the network with respect to its input t, by automatic differentiation (Lecture 11)
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs in JAX, the network

def init(key, sizes):
    params = []
    for n_in, n_out in zip(sizes[:-1], sizes[1:]):
        key, sub = jax.random.split(key)
        W = jax.random.normal(sub, (n_in, n_out)) * np.sqrt(1 / n_in)
        params.append((W, jnp.zeros(n_out)))
    return params

def net(params, t):
    h = t[:, None]
    for W, b in params[:-1]:
        h = jnp.tanh(h @ W + b)
    W, b = params[-1]
    return (h @ W + b)[:, 0]

def x_nn(params, t):
    return net(params, jnp.array([t]))[0]

params0 = init(jax.random.PRNGKey(0), [1, 32, 32, 32, 1])
Weights: one (W, b) per layer, random, scaled by 1/√(inputs) so tanh does not saturate (LeCun et al., 1998)
Forward pass: times in, three tanh layers, displacements out. tanh is smooth, so x″ exists
One time in, one x out: the form jax.grad needs
1 input, 3 × 32 hidden units, 1 output
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs in JAX, the physics

dx = jax.grad(x_nn, argnums=1)      # dx/dt
ddx = jax.grad(dx, argnums=1)       # d2x/dt2

def residual(params, ts):
    x = jax.vmap(lambda t: x_nn(params, t))(ts)
    v = jax.vmap(lambda t: dx(params, t))(ts)
    a = jax.vmap(lambda t: ddx(params, t))(ts)
    return m * a + mu * v + k * x

def data_loss(params):
    return jnp.mean((net(params, jnp.array(t_data)) - x_data) ** 2)

def physics_loss(params):
    return jnp.mean(residual(params, jnp.array(t_phys)) ** 2)

lam = 1e-4

def pinn_loss(params):
    return data_loss(params) + lam * physics_loss(params)
argnums=1: differentiate with respect to t, not the weights. Twice: x″
vmap: evaluate at all collocation points at once. Returns m x″ + μ x′ + k x
Mean squared error on the 10 measurements; mean squared residual at the 40 collocation points
The PINN loss: data plus λ times physics
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs in JAX, training

def train(loss, steps=30_000):
    opt = optax.adam(1e-3)
    params = params0
    state = opt.init(params)

    @jax.jit
    def step(params, state):
        grads = jax.grad(loss)(params)
        updates, state = opt.update(grads, state, params)
        return optax.apply_updates(params, updates), state

    for _ in range(steps):
        params, state = step(params, state)
    return params

p_nn = train(data_loss)
p_pinn = train(pinn_loss)
Adam, learning rate 0.001, from the same starting weights
One step: gradient of the loss with respect to the weights, then Adam's update. jit compiles it
30,000 steps, a few seconds for each network
Same network, same steps: only the loss differs
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, training

A plain network and a PINN trained on ten points of a damped spring-mass, with three hanging masses
  • Play: the curves train, then the three masses move with what each model predicts
  • First 1,000 steps: both fit the 10 points; after 0.36 s both are wrong (error over the second about 0.36 m)
  • Steps 1,300 to 8,500: the physics term pulls the PINN onto the oscillation, error 0.34 m → 0.01 m
  • Plain network: never improves after the data, about 0.45 m
  • Residual at the collocation points: 109 N → 0.7 N
Systems and Toolchains for AI Engineers
06-763 / L13

PINNs, result and limits

0.54 m
plain network: error after the data
0.001 m
PINN: error after the data
163 N
plain network: residual of the equation
5.3 N
PINN: residual, small but not zero

Error: RMSE against the exact solution, 0.36 to 1 s. Residual: RMS over 0 to 1 s. Same experiment as Ben Moseley's blog post.

  • A penalty: the equation holds approximately, and only where you evaluate it
  • λ trades physics against data. A different λ gives a different model
  • Training can fail on harder problems (Krishnapriyan et al., 2021); the two loss terms compete (Wang et al., 2021)
  • One trained PINN is one solution , for this release from 1 m with these , , . Release from 0.5 m and you must train again: the network learned a trajectory, not the law of motion
Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs

Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs, the residual update

  • Residual networks and recurrent networks build a long computation from many small updates of a hidden state :

  • An RNN over time: the same at every step, steps of a fixed size :

  • For the spring-mass, : the network learns how the state jumps from one sample to the next

  • "These iterative updates can be seen as an Euler discretization of a continuous transformation" (Chen et al., 2018)

Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs, two pictures

Recurrent network: a fixed step Δt, one jump at a time
h₀
→
h₀ + Δt f(h₀; θ)
→
h₁
→
h₁ + Δt f(h₁; θ)
→
h₂
→
h₂ + Δt f(h₂; θ)
→
h₃
The same network f, with the same weights θ, at every step. It learns the jump for this Δt.
Neural ODE: a vector field f, integrated by an ODE solver
h(t₀)
→
ODE solver: dh/dt = f(h, t; θ)
→
h(t), any t
The network learns the arrows f. The solver takes as many steps as it needs, and returns h at any time.
Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs, the limit

  • Read the hidden state as a function of time: , with

  • The recurrent update, divided by the step:

  • Let . The left side tends to the derivative of :

  • "In the limit, we parameterize the continuous dynamics of hidden units using an ordinary differential equation (ODE) specified by a neural network" (Chen et al., 2018)

Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs, shrinking the step

What it shows: the spring-mass, with the true . Orange: a recurrent network stepping at a time. Blue: the ODE. Press Play: as shrinks, the orange steps fall onto the blue curve. With a large step they even spiral out.

Euler steps of a spring-mass on its phase plane converging to the exact trajectory
Systems and Toolchains for AI Engineers
06-763 / L13

From RNNs to neural ODEs, what is learned

A neural ODE learns the right-hand side : the vector field of the system. A recurrent network learns one time-discrete realization of it, tied to its step .

Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs, training through the solver

Neural ODE: a differential equation whose right-hand side is a neural network, . An ODE solver computes the output.

network weights θ: the same f(h, t; θ) in every solver step
h(t₀)
solver step
solver step
solver step
h(tᵢ)
loss L: compare h(tᵢ) with the data
∂L/∂θ: back through every solver step, by reverse-mode automatic differentiation (Lecture 11)
a differentiable integrator: every step is made of JAX operations
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs, Patrick Kidger's work

  • Oxford thesis, On Neural Differential Equations (2021): networks and differential equations are "two sides of the same coin"
  • Neural CDEs: "the continuous time analogue of an RNN", for irregular time series (2020)
  • At Google X: Diffrax, differentiable ODE solvers in JAX, and Equinox, neural networks in JAX
  • Both are what the bioreactor below is built with
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs for the bioreactor, the network

mu_scale = 0.3                                 # 1/h: sets the size of mu
x_typical = jnp.array([5.0, 1.0, 10.0, 3.0])   # typical X, P, S, V

class GrowthRate(eqx.Module):
    mlp: eqx.nn.MLP

    def __call__(self, x):
        out = self.mlp(x / x_typical)[0]       # any real number
        return mu_scale * jax.nn.softplus(out) # mu >= 0, of order mu_scale
Two scales. Each state is divided by a typical value, so the network sees numbers near 1. mu_scale sets the size of μ
The unknown term: a small Equinox network (Equinox: neural networks in JAX)
softplus(out) is always positive and about 0.7 at the start, so μ starts near 0.2 1/h, a sensible growth rate. Nothing keeps S ≥ 0

Known balances with one learned term: also called a universal differential equation (UDE), "defined in full or part by a universal approximator" (Rackauckas et al., 2020). Neural ODE, UDE and hybrid model are used interchangeably.

Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs for the bioreactor, the balances

def balances(t, x, args):
    mu_net, Sf = args
    X, P, S, V = x
    mu = mu_net(x)                               # the learned term
    dX = mu * X - F / V * X
    dP = Ypx * mu * X - F / V * P
    dS = F / V * (Sf - S) - mu * X / Yxs
    dV = F
    return jnp.array([dX, dP, dS, dV])
The right-hand side f(t, x) in the form Diffrax wants. args carries the network and the feed concentration Sf
The network gives μ from the four states
The four balances: dX/dt = μX − (F/V)X, ..., dV/dt = F
The four derivatives, stacked
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs in Diffrax, the solve

def simulate(mu_net, x0, ts):
    sol = diffrax.diffeqsolve(
        diffrax.ODETerm(balances),
        diffrax.Tsit5(),
        t0=0.0,
        t1=ts[-1],
        dt0=0.1,
        y0=x0,
        args=(mu_net, x0[2]),
        saveat=diffrax.SaveAt(ts=ts),
        stepsize_controller=diffrax.PIDController(
            rtol=1e-6,
            atol=1e-8,
        ),
    )
    return sol.ys
The right-hand side: our balances
Tsitouras' fifth-order Runge-Kutta method
From 0 to the last sample, first step 0.1 h, from the initial state. args passes the network to balances
Return the states at the sample times; adapt the step to keep the error below the tolerances
  • Every operation is JAX, so jax.grad of anything computed from sol.ys goes back through the solver
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs for the bioreactor, training

def loss(mu_net):
    pred = jax.vmap(simulate, in_axes=(None, 0, None))(mu_net, x0s, ts)
    return jnp.mean(((pred - ys) / y_scale) ** 2)

mu_net = GrowthRate(eqx.nn.MLP(
    in_size=4,
    out_size=1,
    width_size=16,
    depth=2,
    activation=jnp.tanh,
    key=jax.random.PRNGKey(0),
))
opt = optax.adam(optax.cosine_decay_schedule(1e-2, steps, alpha=0.05))
state = opt.init(eqx.filter(mu_net, eqx.is_array))

@eqx.filter_jit
def step(mu_net, state):
    value, grads = eqx.filter_value_and_grad(loss)(mu_net)
    updates, state = opt.update(grads, state, mu_net)
    return eqx.apply_updates(mu_net, updates), state, value

for _ in range(10_000):
    mu_net, state, value = step(mu_net, state)
Simulate all three batches, compare with their 31 samples, each state divided by its spread
4 inputs, two hidden layers of 16 tanh units, 1 output
Adam, with a learning rate that decays from 0.01
Sequential: every step simulates, then updates θ. 10,000 steps, 36 s
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs for the bioreactor, the fit

  • Dots: three batches, 31 noisy samples of each state. Lines: the trained neural ODE, simulated from each batch's start
  • Error against the true model: 0.047 g/L in and 0.094 g/L in , within the measurement noise (0.05 and 0.5 g/L)
Systems and Toolchains for AI Engineers
06-763 / L13

Neural ODEs for the bioreactor, the learned growth rate

Takeaway: the network learned well where the data were. It was never told that cells cannot grow without substrate, so at it still gives 1/h, and the balances drive below zero.

Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning

Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning, Gen 2

Gen 2 promised guaranteed constraints. A PINN informs the training with the physics; we want the physics enforced.

Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning, from a penalty to a constraint

Informed (PINN)

small, traded against the data

→

Enforced

to solver tolerance, every bound held

Why enforce?

  • Conservation: mass and energy balances must close
  • Safety: a reactor temperature or a pressure must stay below its limit at every instant
  • Physical bounds: concentrations cannot be negative, the we just saw
  • Quality and regulation: a product specification must be met, not approached
  • No λ to tune, and the guarantee holds outside the training data too
Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning, path constraints

Path constraint: an inequality on the states that must hold at every instant, for all .

  • "Problems of dynamic optimisation with inequality path constraints are common in industrial plants" (Souza et al., 2006)
  • Batch reactor: limits on cooling duty and "an adiabatic temperature rise constraint that enforces safe operation even if cooling is lost" (Biegler, 2017)
  • Plate reactor start-up: "a temperature path constraint, which is straight forward to enforce using a simultaneous method" (Haugwitz et al., 2007)
  • Our bioreactor: for all
Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning, DAEs

DAE (differential-algebraic equation): differential equations plus algebraic equations that hold at every instant.

  • : differential states: levels, concentrations, temperatures
  • : algebraic variables: flows, pressures, rates, fixed at every instant by
  • Process models are full of them: flash equilibrium, flows that balance at a node, a valve equation

Tank manifold, redrawn from Lueg et al. (2026), Fig. 1 (CC BY 4.0). Levels x0 to x3, flows y0 to y4; the algebraic equation x0 = x1 fixes the flow y2.

Systems and Toolchains for AI Engineers
06-763 / L13

Physics-enforced machine learning, neural DAEs

→

  • Same data, same physics. The unknown part, here the damping, is the network output ; the equation is a constraint, not a term in the loss

Neural DAE: a DAE in which some unknown terms are neural networks, trained with the simultaneous approach (Lueg, Alves, Schicksnus, Kitchin, Laird and Biegler, 2026, open access).

Physics is now enforced, rather than informed: the equation is a constraint of the optimization problem

Systems and Toolchains for AI Engineers
06-763 / L13

Sequential or simultaneous

Systems and Toolchains for AI Engineers
06-763 / L13

Training a dynamic model, the problem

  • Fit a dynamic model to data: an optimization problem with the model as a constraint

  • : measurements. : parameters, or network weights

  • A dynamic optimization problem. Two classic ways to solve it (Biegler, 2007):

    • Sequential: simulate, then optimize
    • Simultaneous: discretize, then optimize
Systems and Toolchains for AI Engineers
06-763 / L13

Sequential approach, simulate then optimize

Sequential approach (single shooting): guess , simulate with an ODE solver, compute the loss and its gradient, update , repeat.

1. Guess

The loss over : one narrow valley

2. Simulate


3. Loss and gradient


and its gradient

4. Update

(BFGS: a quasi-Newton step)

The loss after each update

Repeat: every iterate is a full ODE simulation, costly!
Systems and Toolchains for AI Engineers
06-763 / L13

From sequential to simultaneous

One simulated curve for the current parameters Time cut into finite elements Every state at every point is an unknown, started on the data Residuals of the model equations at each point The solved NLP: points on a trajectory that obeys the model

Discretize first, then optimize: the states become variables, the ODE becomes algebraic equations.

Systems and Toolchains for AI Engineers
06-763 / L13

Simultaneous approach, discretize then optimize

Simultaneous approach (collocation): every state at every time point is an unknown, the ODE becomes algebraic equations, and one optimization problem solves for the states and together.

Lecture 9: a constrained optimization problem

Here is every state value and

Pyomo.DAEPyomo.DAE writes the model in continuous time and discretizes it
IPOPTIPOPT solves the NLP
Systems and Toolchains for AI Engineers
06-763 / L13

Simultaneous approach, polynomials on finite elements

  • In each element the state is a polynomial (here a cubic) through four points; their values are variables of the NLP
  • At each collocation point the polynomial's slope must equal the slope the model gives, : this is how the ODE becomes algebraic equations
  • Gray: the solution of the ODE. Element by element, the polynomials approximate it
Systems and Toolchains for AI Engineers
06-763 / L13

Sequential or simultaneous, iteration by iteration

Estimate the damping and stiffness from 51 noisy samples over 2 s, starting from , . True: ,

Iterates of single shooting and of the simultaneous approach on the same estimation problem
Systems and Toolchains for AI Engineers
06-763 / L13

Sequential or simultaneous, the result

Sequential: 37 BFGS iterations

, SSE 5.23: a local minimum

Simultaneous: 29 interior-point iterations

, SSE 0.036: the truth gives 0.038

A slightly wrong oscillation fits the data badly everywhere. The simultaneous states stay pinned to the data while the model equations are repaired.

Systems and Toolchains for AI Engineers
06-763 / L13

Sequential or simultaneous, side by side

Sequential
Simultaneous
at every iterate
Model equations hold
only at convergence
the parameters θ
Unknowns
θ and every state at every point
can stall in a local minimum
Oscillating or unstable dynamics
states pinned to the data
penalties, checked afterwards
Path constraints, S ≥ 0
rows of the NLP, enforced
rewritten as ODEs first
DAEs
written at the collocation points
adaptive, error-controlled
Time steps
a fixed grid, chosen beforehand
mini-batches on GPUs
Big data, big networks
the NLP grows with the data
Systems and Toolchains for AI Engineers
06-763 / L13

Training neural DAEs

Systems and Toolchains for AI Engineers
06-763 / L13

Neural DAEs, the training problem

  • Every trajectory + regularization
  • Differential equations
  • Algebraic equations
  • Inequalities, such as S ≥ 0
  • The network: the unknown terms
  • Initial conditions
  • Every line at every collocation point of every batch → the physics holds exactly at the solution
  • One for all batches → one network, reusable for new starting charges and feed profiles
Systems and Toolchains for AI Engineers
06-763 / L13

Neural DAEs, the fed-batch reactor again

Inference: use the trained model to predict a case it has never seen, with the network weights fixed.

  • Same three training batches for both models
  • New batch: 0.2 g/L of cells, 7 g/L of substrate, 0.9 L, for 60 h
  • Neural ODE: integrate the learned model (sequential)
  • Neural DAE: solve the balances as an NLP, with (simultaneous)
Systems and Toolchains for AI Engineers
06-763 / L13

Neural DAEs, inference on a new batch

Substrate in a new fed-batch run: the true model, a neural ODE and a neural DAE, with three vessels
Systems and Toolchains for AI Engineers
06-763 / L13

Neural DAEs, the result

0.07 g/L
lowest substrate, true mechanistic model
−0.37 g/L
sequential approach / neural ODE: below zero for 22 of 60 h
0.05 g/L
neural DAE: never below zero

Feasible is not the same as accurate: the substrate error (RMSE) is 0.19 g/L for the neural ODE and 0.30 g/L for the neural DAE.

Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE

SiNDAE: Simultaneous Neural Differential-Algebraic Systems of Equations. A Python package that learns unknown terms in dynamical systems from noisy time series while fully satisfying the constraints. pip install sindae, free and open source: code, documentation

Physics modeling
PyomoPyomo: write the model's equations and constraints
Pyomo.DAEPyomo.DAE: derivatives in time, and their discretization
Machine learning
JAXJAX: the network and its derivatives
ONNXONNX: a file format to export the trained network
OMLTOMLT: puts trained networks into Pyomo models
Constrained optimization
IPOPTIPOPT: interior-point solver for the NLP
POUNCEPOUNCE: interior-point solver, no license needed
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, the problem in code: the class

class FedBatchBioreactorProblem(ProblemDefinition):
    def __init__(
        self,
        params,
        ics,
        input_dim,
        z_dim,
        t_span,
        nfe,
        ncp,
        obs_times=None,
        obs_values=None,
        obs_dim=None,
    ):
        super().__init__(ics, input_dim, z_dim, t_span, nfe, ncp,
                         obs_times, obs_values, obs_dim)
        self.params = params
A problem is a subclass of SiNDAE's ProblemDefinition
Its arguments: the initial charge of each batch, 4 states and 1 learned term, the time span, the finite elements and collocation points, and the measurements
ProblemDefinition stores them
Plus our parameters: feed rate and yields
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, the problem in code: the variables

    def build_trajectory(self, block, traj_idx):
        p = self.params
        t0 = self.t_span[0]
        x0 = self.ics[traj_idx]
        Sf = float(x0[2])                # feed concentration = initial substrate
        Feed, Ypx, Yxs = p["Feed"], p["Ypx"], p["Yxs"]

        block.t = dae.ContinuousSet(bounds=self.t_span)
        block.x = pyo.Var(
            block.t,
            range(self.input_dim),
            domain=pyo.NonNegativeReals,
            initialize=1.0,
        )
        block.z = pyo.Var(block.t, range(self.z_dim), initialize=0.1)
        block.dxdt = dae.DerivativeVar(block.x, wrt=block.t)
Called once per batch: that batch's initial charge, feed and yields
Time, continuous: Pyomo.DAE discretizes it later
The states X, P, S, V at every time, non-negative: here is S ≥ 0
z = μ: a free variable the network will set. dx/dt: the derivatives of the states
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, the problem in code: the balances

        @block.Constraint(block.t, range(self.input_dim))
        def diffeq(b, t, s):
            mu = b.z[t, 0]               # growth rate, learned by the network
            X, P, S, V = b.x[t, 0], b.x[t, 1], b.x[t, 2], b.x[t, 3]
            if s == 0:
                return b.dxdt[t, 0] == mu * X - Feed * (X / V)
            elif s == 1:
                return b.dxdt[t, 1] == Ypx * mu * X - Feed * (P / V)
            elif s == 2:
                return b.dxdt[t, 2] == Feed * (Sf - S) / V - mu * (X / Yxs)
            else:
                return b.dxdt[t, 3] == Feed

        for j in range(self.input_dim):
            block.x[t0, j].fix(float(x0[j]))

    def get_input_vars(self, block, t):
        return [block.x[t, j] for j in range(self.input_dim)]

    def get_output_vars(self, block, t):
        return [block.z[t, 0]]
A constraint for every time t and every state s: these are rows of the NLP
The four balances, one per branch. No formula for μ
Fix the initial charge of the batch
The network: the states in, z = μ out
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, the network and the training stages

mlp = SimpleMLP(
    in_size=INPUT_DIM,
    out_size=Z_DIM,
    widths=[20, 20],
    activations=[jax.nn.softplus] * 2,
    key=jax.random.PRNGKey(SEED),
)
smoother_config = SmootherConfig(smooth_coef=10.0)
pretrain_config = PretrainConfig(
    epochs=200,
    batch_size=32,
    reg_coef=1e-3,
)
simul_config = SimultaneousConfig(
    use_gbm=True,
    reg_coef=1e-3,
)
solver_options = SolverConfig(
    tol=1e-6,
    max_iter=1000,
)
The network: 4 states in, μ out, two hidden layers of 20 softplus units, smooth as the solver needs
Stage 1, smoother: smooth trajectories through the data, μ free
Stage 2, pretraining: fit the network to the smoother's states and μ, with Adam
Stage 3, simultaneous: one NLP with the weights, every state and the balances. The network is evaluated in JAX, so the solver approximates second derivatives (L-BFGS)
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, training

problem = FedBatchBioreactorProblem(
    params=FB_PARAMS,
    ics=BATCH_ICS,
    input_dim=INPUT_DIM,
    z_dim=Z_DIM,
    t_span=T_SPAN,
    nfe=NFE_TRAIN,
    ncp=NCP_TRAIN,
    obs_dim=OBS_DIM,
    obs_times=obs_times,
    obs_values=obs_values,
)
The three batches over 40 h, cut into 40 finite elements of 3 collocation points
The measurements: 31 samples of the 4 states in each batch
model = HybridDAE(
    method="simultaneous",
    nlp_solver="pounce",
    net=mlp,
    smoother=smoother_config,
    pretrain=pretrain_config,
    train=simul_config,
    solver_options=solver_options,
    unfix_io=True,
)
model.fit(problem)
Simultaneous training, solved by POUNCE, with the network and the three stages
fit: smoother, pretraining, then the full NLP. About 30 s
Systems and Toolchains for AI Engineers
06-763 / L13

SiNDAE, prediction

new_problem = FedBatchBioreactorProblem(
    params=FB_PARAMS,
    ics=np.array([[0.20, 0.0, 7.0, 0.90]]),
    input_dim=INPUT_DIM,
    z_dim=Z_DIM,
    t_span=(0.0, 60.0),
    nfe=45,
    ncp=3,
    obs_dim=OBS_DIM,
)
prediction = model.predict(
    new_problem,
    slack_coef=1e-5,
)
The same class, a new batch: 0.2 g/L of cells, 7 g/L of substrate, 0.9 L
60 h, 20 h past the training batches, and no measurements
predict: the network fixed, the DAE solved with S ≥ 0 held. slack_coef: the price on moving μ off the network where a constraint would break
Systems and Toolchains for AI Engineers
06-763 / L13

Monoclonal antibodies, why glycosylation

  • Monoclonal antibodies (mAbs): medicines for cancers, autoimmune and infectious diseases, mostly made by CHO cells in fed-batch reactors
  • Glycosylation: sugar chains (glycans) attached to the antibody inside the cell, in the Golgi
  • The glycans are a critical quality attribute (CQA): "the terminal sugars of Fc glycans have been shown to be critical for safety or efficacy" (Reusch and Tejada, 2015)
  • They also change how long the drug stays in the body, and immune reactions to it
Systems and Toolchains for AI Engineers
06-763 / L13

Monoclonal antibodies, a hybrid glycosylation model

  • 9 states: viable cells, glucose, glutamine, lactate, ammonia, antibody, and three glycoforms G0F, G1F, G2F
  • Mechanistic: six culture balances, and three glycan equations with the rates above
  • Algebraic: nucleotide sugar donors (NSD) at quasi-steady state, set by glucose at every instant: a DAE
  • Learned: the growth rate, , 9 → 16 → 16 → 1, tanh
  • Why : it is never measured, and its usual law (Monod, inhibited by lactate and ammonia) is an assumption
Systems and Toolchains for AI Engineers
06-763 / L13

Monoclonal antibodies, the multiscale equations

Process scale
the bioreactor, 6 ODEs

Cell scale
sugar donors, algebraic

Golgi scale
glycans, 3 ODEs

, . Red: the learned term. The algebraic sugar donors make the whole model a DAE.

Systems and Toolchains for AI Engineers
06-763 / L13

Monoclonal antibodies, results

  • Data at two scales, one model: process measurements (blue box) and product-quality measurements, the glycan fractions (orange box)
  • Synthetic data: the mechanistic model plus noise, half of the samples kept, three runs of 10 days
  • Blue: true model. Green: neural DAE. Red dots: noisy data. Limits: 3 of more than 33 glycoforms; quasi-steady sugar donors; the learned oscillates
Systems and Toolchains for AI Engineers
06-763 / L13

Constraints inside the network

Systems and Toolchains for AI Engineers
06-763 / L13

Projection layers

Projection layer: a last layer that moves the network's raw output to the closest point that satisfies the constraints.

  • Input: the hot-water flow . Outputs: the two outlet temperatures
  • First law: heat given up = heat taken,
  • Linear in the outputs, . A plain network can still break it: energy from nothing
flow ṁh
network
weights θ
raw ỹ
projection, fixed
violation v = aᵀỹ − b
y = ỹ − a v / (aᵀa)
y
loss against the data
gradient ∂L/∂θ flows back through the projection, a linear map with no weights
Systems and Toolchains for AI Engineers
06-763 / L13

Projection layers, training

Training a network with a projection layer on a heat exchanger's energy balance
  • Noisy data: the measurements break the energy balance by 5.4 kW on average
  • Without the layer: the predictions create or destroy up to 10 kW
  • With the layer: kW, at every epoch
  • Error in the outlet temperatures: 0.49 → 0.33 °C: the layer removes the noise that breaks the balance
Systems and Toolchains for AI Engineers
06-763 / L13

Four approaches side by side

PINN

Physics in: the loss

Learns: one solution

Trained with: Adam

physics favored, not enforced

Neural ODE

Physics in: the right-hand side

Learns: the vector field

Trained with: unconstrained solver + Adam (sequential)

balances held; bounds not

Neural DAE

Physics in: the constraints of one NLP

Learns: the unknown terms

Trained with: nonlinear constrained optimization solver (simultaneous)

physics enforced, training and prediction

Projection layer

Physics in: the last layer

Learns: a static map

Trained with: Adam

physics enforced: linear balances; nonlinear ones in new work (KKT-Hardnet)

Systems and Toolchains for AI Engineers
06-763 / L13

Four approaches, as optimization problems

PINN: physics in the objective

Variables: the weights θ
Constraints: none
Solver: Adam, unconstrained
The physics is one more term to make small: the residual r is traded against the data through λ.
Neural ODE: physics inside a simulation

Variables: the weights θ
Constraints: none; the ODE is solved inside the objective
Solver: ODE solver + Adam (sequential)
The balances hold at every iterate, because every iterate is a simulation. Bounds such as S ≥ 0 are not in the problem.
Neural DAE: physics as constraints

Variables: θ, and every state at every collocation point
Constraints: the equations and bounds, at every point
Solver: nonlinear constrained optimization (simultaneous)
The physics holds at the solution, to solver tolerance, in training and in prediction.
Projection layer: physics in the last layer

Variables: the weights θ
Constraints: none outside; the layer solves a small problem in closed form
Solver: Adam, unconstrained
Every output satisfies Ay = b exactly. Nonlinear balances need the layer to solve a nonlinear problem (KKT-Hardnet).
Systems and Toolchains for AI Engineers
06-763 / L13

Not only neural networks

The learned term can be any model that fits data. Four families used with physics:

A Gaussian process: a mean curve with a shaded uncertainty band through six points

Gaussian processes

The unknown right-hand side of a fed-batch mAb culture model, with its uncertainty (Cruz-Bournazou et al., 2022)
Sparse regression: coefficients of eight candidate terms, two nonzero

Sparse regression (SINDy)

Learns the equation: keeps the few terms of a library that reproduce the data (Brunton, Proctor and Kutz, 2016)
An expression tree for mu max times S divided by K plus S

Symbolic regression

A search over formulas returns a closed form, such as μmax S / (K + S) (PySR, Cranmer, 2023)
A decision tree that splits on substrate and biomass to give a growth rate

Tree ensembles

A random forest sets the coefficients of a mAb process's mass balances (Nemoto et al., 2025); OMLT puts trees into optimization models too

The physics sits in the loss, the right-hand side or the constraints, whatever the learner.

Systems and Toolchains for AI Engineers
06-763 / L13

Limitations and trade-offs

  • PINNs: a penalty; can fail on stiff or multiscale problems; one solution per condition. Good with few data and a known equation
  • Neural ODEs (sequential): bounds not enforced; local minima on oscillating dynamics. Scale to big data and GPUs
  • Neural DAEs (simultaneous): the NLP grows with network and data; smooth activations; a fixed time grid. Use when constraints must hold exactly: safety, quality, balances, DAEs
  • Projection layers: one prediction at a time; inequalities and nonlinear constraints cost a solve per pass
  • All hybrids: the mechanistic part must be right; feasible is not the same as accurate
Systems and Toolchains for AI Engineers
06-763 / L13

Worked examples

Every code block on today's slides is a cell in these two notebooks, in the same order:

  • l13-pinn-jax.ipynb: the spring-mass PINN from scratch in JAX, against a plain network
  • l13-neural-dae.ipynb: the fed-batch bioreactor as a neural ODE in Diffrax and as a neural DAE in SiNDAE
Systems and Toolchains for AI Engineers
06-763 / L13

Before next time

  • The miniproject is due Friday 10-09
  • Practice module for this session: on the course site
  • Run both worked examples; the neural DAE one installs SiNDAE and Diffrax with pip install sindae diffrax
Systems and Toolchains for AI Engineers

Why physics 12, PINNs 16, RNN to neural ODE 16, physics enforced 7, seq vs sim 10, neural DAEs + SiNDAE + mAb 16, layers 8, beyond networks 2, close 5. Long deck: skip the code slides live if late, they are in the notebooks.

One bioreactor at both ends. One spring-mass in the middle.

The drawing is from my seminar. mAbs come back at the end with glycosylation.

Mass balances: accumulation = in - out + made - used. Every term is bookkeeping except mu.

Left: my multiscale picture, process to cell to Golgi. Right: what one cell looks like inside, simplified.

The gray band is time the training batches never reached. Today: how to get from the purple curve to the blue one.

Plots from my seminar. Gen 1: exploration, the green diamond is where the uncertainty is largest.

Click order: arrows, the two challenges, then the hybrid box, then what Gen 2 promises.

CSTR: continuous stirred-tank reactor (Bequette's exothermic example), the spike is ignition. Same pattern as mu in the bioreactor.

Let a student drag it. About 3 oscillations per second; use slow x4.

Gravity only shifts the rest position, so x is measured from rest and gravity drops out.

The gray band is all the data there is. Everything to its right is extrapolation.

Same layout comes back for neural DAEs. Raissi, Perdikaris, Karniadakis 2019.

The physics loss only sees the green points. Between them the equation is not checked at all.

Residual is zero for any function that obeys the equation.

From scratch on purpose, no Flax or Equinox here. Same code as l13-pinn-jax.ipynb, section 2.

Lecture 11's training loop, unchanged. The only new thing is the loss.

One Play: training, then the masses. Watch the orange mass stop oscillating after 0.36 s.

A PINN learns one solution. From here on we learn the right-hand side, which works for any starting point.

Euler: new state = old state + step size x slope at the old state.

Top: what an RNN computes. Bottom: what a neural ODE computes. Next slide: why the bottom is the limit of the top.

Same f on both sides. The recurrent network is the forward Euler discretization of this ODE.

Left: phase plane, position against velocity; the arrows are f. Right: x(t). Errors 0.85, 0.25, 0.106, 0.049, 0.024 m: halving the step halves the error.

Redrawn after Chen et al. Fig. 1. Same arrows work for any starting state and any sampling times.

Forward: simulate. Backward: the gradient goes through the solver. Over an RNN: data at any times, the solver picks the steps, the state can be physical.

Remember "simultaneous": it is the next section.

Algebraic equations are constraints that hold all the time: exactly what we want to enforce.

Compare with the PINN slide: the physics moved from the loss to the constraints.

Biegler 2017: sequential approaches are known to fail on unstable dynamic systems. Here: stuck at SSE 5.2.

Five clicks. Frame 4: red bars are the residuals of the trapezoid rule at the starting point. Frame 5: the POUNCE solution.

x_jk: the state in element j at point k. The polynomial's slope at a point is a weighted sum of its point values, so each line is algebraic in the unknowns. Details: Biegler 2007, and his open lecture slides on collocation. Next slide: what these polynomials look like.

Same data, start, model. Sequential: Runge-Kutta inside BFGS. Simultaneous: trapezoid on 200 steps, POUNCE, 404 unknowns. Iteration 0: simultaneous on the data, ODE broken.

From k = 250, 350 or 450 both methods find mu 0.98, k 400: the sequential failure is about the start.

Green: the better side for that row. Multiple shooting sits between the two.

The paper's Eq. 4. Smooth activations (tanh): the solver uses second derivatives of the network.

Play. The neural ODE vessel drains through the floor at 25 h. The neural DAE sits on the floor.

As in the SiNDAE fed-batch example, with F P / V in the product balance. This and the next two slides are one cell of l13-neural-dae.ipynb.

CQA: must stay within a range for the product to be released. Process scale is mechanistic; Golgi scale is poorly known.

Proof of concept from my group: model after Kotidis et al. 2019. Precursor -> G0F (fucosylation) -> G1F -> G2F (galactosylation).

Six culture ODEs: viable cells, glucose, glutamine, lactate, ammonia, antibody. Glycan rates from Kotidis et al. 2019.

Counterflow, water on both sides: hot in at 90 C, cold in at 20 C and 1 kg/s. a = (m_h, m_c) changes with the input; the projection stays linear in y~. Any linear balance Ay = b works the same way. It is part of the network in training and in prediction, not a fix-up afterwards.

Play: the pulse runs forward through the layers, the red one back through the projection.

Green: where the physics sits in each formulation.

Click through. Same data-fit term every time; what changes is where the physics sits and what the solver must handle.

GPs also come with physics built into the prior: Raissi and Karniadakis 2017 learn the parameters of linear ODEs this way. SINDy and symbolic regression return equations a person can read.