Making unknown physical constants trainable variables, recovering diffusivity, reaction rates, and material properties from 40 noisy sensor readings using physics-informed neural networks.
Most introductions to PDEs present what is called the forward problem: given a governing equation with fully known coefficients and initial/boundary conditions, compute the solution field. For the 1D heat equation, this reads: given thermal diffusivity $\alpha$, solve
for the temperature field $u(x,t)$. This is well-posed, and numerical methods, finite differences, finite elements, spectral methods, have been refined over decades to solve it efficiently and to guaranteed accuracy.
The inverse problem runs in the opposite direction. You have sensor measurements $\{u(x_i, t_i)\}_{i=1}^{N}$, possibly corrupted by noise, and the parameter $\alpha$ is unknown. You want to infer $\alpha$ from the data. This reversal is not merely a notational switch, it changes the problem's mathematical character entirely. Inverse problems are generically ill-posed in the sense of Hadamard: existence, uniqueness, and continuous dependence on data may all fail without additional regularization.
The inverse problem is often the real engineering question. Consider a few motivating examples:
In each case, the physics is not just background knowledge, it is the mechanism that makes parameter identification possible from data that would otherwise be severely underdetermined.
The key conceptual move in PINN-based inverse problems is simple but powerful: treat the unknown physical parameter as an additional trainable variable in PyTorch, optimized jointly with the network weights. The parameter is no longer a hyperparameter or an input, it is a learnable scalar that participates in the automatic differentiation graph alongside every weight and bias.
For the 1D heat equation with unknown diffusivity $\alpha^\star = 0.4$, initialized at $\alpha_{init} = 1.5$ (far from the truth), the network class becomes:
class InversePINN(nn.Module): def __init__(self, alpha_init=1.5): super().__init__() # Standard network layers for u(x, t) self.net = nn.Sequential( nn.Linear(2, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1) ) # KEY: alpha as a trainable parameter, stored in log-space # This guarantees alpha > 0 and improves the optimization landscape self.log_alpha = nn.Parameter( torch.tensor(np.log(alpha_init), dtype=torch.float32) ) def forward(self, x, t): inp = torch.cat([x, t], dim=1) return self.net(inp) @property def alpha(self): # Recover alpha from log-space during evaluation return torch.exp(self.log_alpha)
Storing log(α) rather than α directly serves two purposes. First, it enforces the constraint $\alpha > 0$ unconditionally, exponentiation always returns a positive value, eliminating the need for clipping or projection steps that would break gradient flow. Second, it improves the optimization landscape. Physical parameters often span several orders of magnitude (think permeabilities from $10^{-3}$ to $10^{3}$ m²/s), and working in log-space makes the gradient-based search scale-invariant. A step of $0.1$ in log-space corresponds to a multiplicative factor of $e^{0.1} \approx 1.105$ regardless of the current magnitude of $\alpha$.
More formally, if the loss $\mathcal{L}$ is a function of $\alpha$, the gradient with respect to $\log\alpha$ is
This natural scaling by $\alpha$ is the continuous analog of the relative update in Adam: the effective learning rate is proportional to the current parameter value, which prevents the optimizer from taking catastrophically large absolute steps when $\alpha$ is large and tiny steps when it is small.
The training data consists of $N = 40$ sensor observations drawn from the true solution $u^\star(x,t)$ with additive Gaussian noise at 2% relative amplitude:
where $\sigma_{u^\star}$ is the standard deviation of the true field. This noise level is representative of high-quality inline spectroscopic measurements in a pharmaceutical manufacturing context (typically $\pm 1$–$3\%$ relative error for NIR-based concentration estimates after multivariate calibration).
The combined loss has two terms that work in concert:
where $\theta$ denotes all network weights and biases, $\hat{u}$ is the network output, $\hat{u}_t$ and $\hat{u}_{xx}$ are computed via automatic differentiation, and $N_f$ is the number of collocation points (typically $N_f \gg N$, here $N_f = 10{,}000$).
Backpropagation computes gradients with respect to both $\theta$ and $\log\alpha$ simultaneously. In PyTorch, both are nn.Parameter objects registered in the same module, so the optimizer sees them as part of the same parameter group. The gradient of the physics loss with respect to $\log\alpha$ carries information about which direction to move $\alpha$ to reduce the PDE residual, this is the signal that drives parameter identification.
def train_step(model, optimizer, x_data, t_data, u_data, x_f, t_f, lam=1.0): optimizer.zero_grad() # --- Data loss --- u_pred = model(x_data, t_data) loss_data = nn.MSELoss()(u_pred, u_data) # --- Physics loss (PDE residual) --- x_f.requires_grad_(True) t_f.requires_grad_(True) u_f = model(x_f, t_f) u_t = torch.autograd.grad(u_f.sum(), t_f, create_graph=True)[0] u_x = torch.autograd.grad(u_f.sum(), x_f, create_graph=True)[0] u_xx = torch.autograd.grad(u_x.sum(), x_f, create_graph=True)[0] residual = u_t - model.alpha * u_xx # model.alpha = exp(log_alpha) loss_physics = (residual**2).mean() loss = loss_data + lam * loss_physics loss.backward() optimizer.step() return loss_data.item(), loss_physics.item(), model.alpha.item()
A natural question: why not simply fit a neural network to the 40 data points without any physics term, and then estimate $\alpha$ from the fitted function by numerical differentiation? The answer is that 40 points in a continuous domain are radically insufficient to uniquely determine a smooth function, let alone its second spatial derivative. An unconstrained network will interpolate the data and produce a $u(x,t)$ that may be arbitrarily wiggly, making $\partial^2 u / \partial x^2$ meaningless.
Classical regularization (Tikhonov, $L^2$ penalty on weights or on derivatives of $u$) enforces smoothness, it prevents the fit from being too wiggly. But smoothness is a geometric property unrelated to the physics. A smooth but physically inadmissible field can still yield a wildly wrong estimate of $\alpha$.
The physics loss enforces a structural constraint: the fitted field must belong to the solution manifold of the heat equation for some value of $\alpha$. This is a far stronger requirement than smoothness. The solution manifold of the heat equation is a much smaller subset of all smooth functions, it is parametrized by $\alpha$ (and the initial/boundary conditions), so fitting onto it simultaneously constrains both the field and the parameter.
The training log from notebook 07_inverse_parameter_estimation.ipynb shows the convergence trajectory. The optimizer is L-BFGS for the final phase (after an Adam warm-up for 3,000 epochs), which is standard for PINN inverse problems because L-BFGS exploits second-order information and converges much faster once the network is in a good basin.
The classical framework for assessing identifiability is the Fisher Information Matrix (FIM). For a parameter vector $\boldsymbol{\theta}$ given noisy observations $\mathbf{y} = \mathbf{f}(\boldsymbol{\theta}) + \boldsymbol{\epsilon}$, $\boldsymbol{\epsilon} \sim \mathcal{N}(0, \sigma^2 I)$, the FIM is
where $\mathbf{J}_{ij} = \partial f_i / \partial \theta_j$ is the Jacobian of the model outputs with respect to parameters. A parameter is identifiable if and only if $\mathcal{I}$ is positive definite. For the single parameter $\alpha$, this reduces to asking whether $\sum_i (\partial u / \partial \alpha)^2 > 0$, i.e., whether the model output is sensitive to the parameter at the observed locations. The PINN computes this sensitivity automatically via backpropagation.
The inverse of the FIM gives the Cramér–Rao lower bound on parameter estimation variance: $\text{Var}(\hat{\alpha}) \geq 1/\mathcal{I}(\alpha)$. This bound is achievable by the maximum likelihood estimator asymptotically, and the PINN inverse problem can be interpreted as an approximate MLE when the noise is Gaussian and the physics loss is used as the likelihood rather than a penalty.
The inverse PINN framework extends naturally to ODE systems. Consider the consecutive reaction A → B → C with mass-action kinetics:
with initial conditions $[A](0) = 1$, $[B](0) = [C](0) = 0$. Suppose the true constants are $k_1^\star = 0.3$, $k_2^\star = 0.7$, and you have only 15 noisy concentration measurements for each species at scattered time points. The PINN inverse setup:
class KineticsInversePINN(nn.Module): def __init__(self, k1_init=1.0, k2_init=1.0): super().__init__() # Three networks: one per species concentration self.net_A = make_mlp(1, 64, 4, 1) self.net_B = make_mlp(1, 64, 4, 1) self.net_C = make_mlp(1, 64, 4, 1) # Both rate constants learnable, stored in log-space self.log_k1 = nn.Parameter(torch.tensor(np.log(k1_init))) self.log_k2 = nn.Parameter(torch.tensor(np.log(k2_init))) def physics_residuals(self, t): t = t.requires_grad_(True) A, B, C = self.net_A(t), self.net_B(t), self.net_C(t) k1, k2 = torch.exp(self.log_k1), torch.exp(self.log_k2) dA = grad(A, t)[0]; dB = grad(B, t)[0]; dC = grad(C, t)[0] r1 = dA + k1 * A # dA/dt = -k1*A r2 = dB - k1 * A + k2 * B # dB/dt = k1*A - k2*B r3 = dC - k2 * B # dC/dt = k2*B return r1, r2, r3
The total loss is $\mathcal{L} = \mathcal{L}_{data,A} + \mathcal{L}_{data,B} + \mathcal{L}_{data,C} + \lambda(\mathcal{L}_{phys,A} + \mathcal{L}_{phys,B} + \mathcal{L}_{phys,C})$. With this setup and 4,000 L-BFGS steps, the notebook recovers $\hat{k}_1 = 0.298 \pm 0.004$ and $\hat{k}_2 = 0.703 \pm 0.006$ from 45 noisy data points, relative errors of 0.7% and 0.4% respectively.
The data loss and physics loss typically have very different magnitudes at initialization. If $\mathcal{L}_{physics} \gg \mathcal{L}_{data}$, the optimizer essentially ignores the data, and the network fits the PDE for the wrong parameter. A simple but effective heuristic: compute both losses at initialization and set $\lambda = \mathcal{L}_{data}^{(0)} / \mathcal{L}_{physics}^{(0)}$. This normalizes both terms to the same order of magnitude at the start of training, preventing one from dominating.
For inverse problems, the choice of $\alpha_{init}$ matters more than in the forward case. If $\alpha_{init}$ is too far from the truth in log-space (say, $|\log\alpha_{init} - \log\alpha^\star| > 3$), the physics loss landscape near the true parameter can be lost in a flat region, and the optimizer converges to a local minimum. A practical strategy: run a coarse grid search over 5–10 values of $\alpha_{init}$ spanning 2–3 decades, identify which gives the lowest combined loss after 500 epochs, and use that as the initialization for the full run. This is cheap because 500 epochs is fast, and it dramatically improves robustness.
Before committing to a full training run, estimate the sensitivity $\partial u / \partial \alpha$ at the proposed sensor locations using the analytic solution (if available) or a forward PINN with a perturbed $\alpha$. If the sensitivity is near zero everywhere, the parameter is not recoverable from those sensor positions regardless of how long you train. This pre-check takes minutes and prevents wasted computation.
A point estimate $\hat{\alpha}$ is insufficient for engineering decisions. Three methods are in common use:
| Method | Description | Cost | Best for |
|---|---|---|---|
| Laplace approximation | Fit a Gaussian to the loss landscape at the MAP estimate using the Hessian | Low (1 Hessian computation) | Near-Gaussian posteriors |
| Ensemble PINNs | Train $M$ independent networks with different random seeds; use spread of $\hat{\alpha}_m$ as uncertainty | $M\times$ training cost | Non-Gaussian, multimodal |
| Bayesian PINN (B-PINN) | Treat all weights and $\alpha$ as random variables; approximate posterior via HMC or variational inference | Very high | Full posterior needed |
| MC Dropout | Enable dropout at inference; average $T$ stochastic forward passes | Low (one network) | Quick confidence estimates |
The FDA's Process Analytical Technology (PAT) guidance (2004) [7] defines PAT as "a system for designing, analyzing, and controlling manufacturing through timely measurements of critical quality and performance attributes of raw and in-process materials and processes, with the goal of ensuring final product quality." This is precisely the problem class that PINN inverse methods are designed to address.
Modern PAT installations include inline NIR and Raman spectrometers, focused beam reflectance measurement (FBRM) for particle size, and inline HPLC. These instruments generate high-frequency concentration-time profiles with a precision of $\pm 1$–$3\%$ relative error. The PINN inverse framework can exploit this data to recover kinetic parameters in near-real-time, closing the loop between measurement and process understanding.
The pharmaceutical development context under ICH Q8(R2) [8] requires a Design Space to be established for critical process parameters (CPPs). Traditionally, this requires a formal DoE, a resource-intensive campaign that may require 16–50 batch runs at varied process conditions. Using PINN-based parameter estimation, a process engineer can recover the kinetic model parameters from a single batch run (with sufficient PAT coverage) and then use the calibrated mechanistic model to compute the Design Space analytically. This is a fundamentally different paradigm: model-enabled rather than empirically enumerated Design Space.