Correcting model misspecification in physics-informed neural networks (PINNs)

Consider the following nonlinear ODE/PDE describing a physical system:

where $\Omega$ is the domain, $\lambda$ is the model parameter, $u$ is the sought solution, $f$ is the source term, $b$ is the boundary term, $\mathcal{F}_{\lambda}$ and $\mathcal{B}_{\lambda}$ are operators parameterized by $\lambda$, defining the equation and the boundary condition. The PINN method is capable of solving ODEs/PDEs given initial/boundary conditions as well as learning governing equations and identifying model parameters from data yang2021b, zou2022neuraluq, lu2021deepxde, cai2021physics, pang2019fpinns, kharazmi2021hp, jagtap2021extended, jagtap2020conservative, zou2023hydra, linka2022bayesian, chen2023leveraging, lin2022multi, leung2022nh, tang2023pinns, tang2023adversarial, meng2021multi, sarabian2022physics, guo2022monte, gao2023failure, YIN2023105424. A schematic view of the PINN method is presented in Fig. 1, in which a DNN with $\sigma$ as the activation function and $\theta$ as the parameter is built to approximate $u$ with $u_{\theta}$ and $f,b$ with $f_{\theta},b_{\theta}$, respectively, via automatic differentiation raissi2019physics.

In the present study, we focus on using PINNs to learn governing equations from data. The dataset for learning governing equations with PINNs is denoted by $\mathcal{D}$, which can be expressed as $\mathcal{D}=\mathcal{D}_{u}\cup\mathcal{D}_{f}\cup\mathcal{D}_{b}$, where $\mathcal{D}_{u}=\{x^{u}_{i},u_{i}\}_{i=1}^{N_{u}}$ is the set of data for $u$, and $\mathcal{D}_{f}=\{x^{f}_{i},f_{i}\}_{i=1}^{N_{f}}$ and $\mathcal{D}_{b}=\{x^{b}_{i},b_{i}\}_{i=1}^{N_{b}}$ are sets of data for the physics. The model parameter $\lambda$, which defines the equation, is then obtained by optimizing the following PINN loss function with the gradient-descent method and its variants ruder2016overview:

where

$w_{u},w_{f},w_{b}$ are belief weights for balancing different terms, and $||\cdot||_{2}$ is the $\ell^{2}$-norm for finite-dimensional vector. We note that in certain cases we may not have data corresponding to the boundary/initial conditions, i.e., the last term in $\mathcal{L}_{PDE}(\theta)$, and we can directly drop this term in the computations.

The computational accuracy of PINNs in learning governing equations strongly depends on the encoded physical models, i.e. $\mathcal{F}_{\lambda}$ and/or $\mathcal{B}_{\lambda}$ in Eq. (1). However, the correct specification of the physics is not always guaranteed, leading to significant discrepancy between the observational data and the assumed governing equations. An example is demonstrated in Fig. 2 where the flow is non-Newtonian but the physical model is misspecified as the one for Newtonian flows (see Sec. 3.3.1 for details). As a result, the physics and the data cannot be fitted well at the same time: the training does not land on the region “IV”, where both the losses for the data as well as the residual of the equations are close to zero, regardless of the choice of belief weights in the loss function Eq. (2), i.e. $w_{f}$ for the physics and $w_{u}$ for the data. Three results from PINNs are sampled and displayed in Fig. 2(b): when $u_{\theta}$ fits the data (which come from a non-Newtonian flow), it does not satisfy the governing equations for Newtonian flows; when $u_{\theta}$ is forced to follow the Newtonian flow by setting large $w_{PDE}$, it presents Newtonian behavior and thus cannot fit the data from the non-Newtonian flow well.

The discrepancy between observational data and the assumed governing equations is caused by the physics/model misspecification. In this regard, we propose a general approach developed upon the original PINN framework to quantify the discrepancy as well as to correct the model misspecification. In this work, we focus on the misspecification of the differential operator $\mathcal{F}_{\lambda}$, but the proposed approach can be easily generalized to the boundary/initial operator. Here, we denote the misspecified differential operator by $\tilde{\mathcal{F}}_{\lambda}$ and then the discrepancy can be quantified as follows:

Inspired by the PINN method, in which the sought solution $u$ is modeled with a DNN parameterized by $\theta$, we employ another DNN parameterized by $\psi$ to model $s(x)$, denoted as $s_{\psi}(x)$. The source term $f(x)$ is approximated by $\tilde{\mathcal{F}}_{\lambda}[u_{\theta}](x)+s_{\psi}(x)$, and hence the term penalizing the violation of the physics in the loss function (2) is reformulated as:

We note that $s_{\psi}$ plays two important roles in our approach handling model misspecification: (1) it approximates the function quantifying the model discrepancy, i.e. $s$, and (2) it allows both the data term $\mathcal{L}_{\mathcal{D}_{u}}(\theta)$ and the PDE term $\mathcal{L}_{PDE}(\theta,\psi)$ to be minimized to very small values simultaneously in the training of PINNs such that data are fitted by $u_{\theta}$ and the (corrected) physics are satisfied.

As mentioned in Sec. 2.2, a key step in the present work is to learn the correction term, i.e., $s_{\psi}$, that is used to fix the issue of model misspecifications in PINNs. Generally, the measurements in real-world applications are noisy and can also be incomplete, leading to the uncertainties in predictions from PINNs psaros2023uncertainty, zou2022neuraluq. We refer to the uncertainties in $s_{\psi}$ as model uncertainty. Here, two typical approaches are employed for quantifying the uncertainties in $s_{\psi}$, i.e., ensemble PINNs psaros2023uncertainty, zou2022neuraluq and B-PINNs yang2021b. The former is computationally efficient but has difficulties preventing overfitting when the data are noisy, while the latter is able to handle the noisy data well but is computationally more expensive. We therefore utilize the ensemble PINNs and B-PINNs for cases with clean and noisy data, respectively.

We now briefly review the ensemble PINNs and B-PINNs. The ensemble PINN method follows the deep ensemble method, proposed in lakshminarayanan2017simple, and employs multiple standard training of PINNs independently. That is, PINNs with the PDE loss term in Eq. (4) are trained multiple times with random initialization of the DNN parameter. Each random initialization corresponds to one initial guess for minimizing the PINN loss function (2) and each training is finding one global minimum. As a result, ensemble PINNs can be considered as identifying different maximum a posteriori (MAP) estimates of the PINN loss function psaros2023uncertainty. On the other hand, the B-PINN method focuses on estimating the posteriors for $\theta$ and $\psi$ based on the Bayes’ rule as follows:

Here $p(\theta|\mathcal{D})$ is the density function for the posterior distribution, $p(\mathcal{D}|\theta)$ for the likelihood, and $p(\theta)$ for the prior. Given independent and identically distributed (i.i.d.) data, the density function for the likelihood distribution can be written as follows:

If we further assume that the data for $u$, $f$ and $b$ follow Gaussian distribution with mean zero and standard deviation, $\sigma_{u}$, $\sigma_{f}$ and $\sigma_{b}$, respectively, then correcting it follows that:

We note that $\sigma_{u}$, $\sigma_{f}$ and $\sigma_{b}$ are often treated as scales of the additive Gaussian noise in measuring data of $u$, $f$and $b$, respectively. Various methods, e.g. Markov Chain Monte Carlo (MCMC) and variational inference, have been proposed to tackle the posterior distribution (5). In this work, we choose to use Hamiltonian Monte Carlo (HMC) method to estimate the posterior distribution by sampling $\theta$ and $\psi$ from it. We also note that here we only consider the homogeneous noise, while more complicated noise models, e.g. the heteroscedastic noise psaros2023uncertainty, are also compatible with the framework.

In both the ensemble PINN and B-PINN methods, we obtain samples of the NN parameters, from which statistics of the prediction are estimated. In general, the mean and the standard deviation are often chosen to represent the prediction and the uncertainty psaros2023uncertainty, respectively. Specifically in this work, we obtain samples of $\theta$ and $\psi$ from ensemble PINNs and/or B-PINNs, denoted by $\{\theta_{i}\}_{i=1}^{M}$ and $\{\psi_{i}\}_{i=1}^{M}$, respectively, where $M$ denotes the number of samples. Then, the predicted means of the sought solution $u$ and the discrepancy $s$ are estimated as follows:

where $\hat{u}$ and $\hat{s}$ denote the mean of $u$ and $s$, respectively. Similarly, the predicted uncertainties are:

where $Var[u]$ and $Var[s]$ denote the variance of $u$ and $s$, respectively. Interested readers are directed to psaros2023uncertainty for a review and zou2022neuraluq for a comprehensive Python library termed NeuralUQ.

We note that upon the computation of $s_{\psi}(x)$, it can be used in various way: (1) its value represents the model discrepancy and its uncetainty represents the uncertainty in the physical model which we refer to as model uncertainty; (2) it corrects the misspecified physical models and identifies a system which explains the data; and (3) it could also be used to obtain mathematical expressions for the missing details of certain physical process via the symbolic regression zhang2023discovering, eastman2022pinn, kiyani2023framework. Specifically in the present work, we employ the symbolic regression method developed in cranmer2023interpretable and the associated Python library PySR to regress analytic mathematical expressions for the unknown physical processes (See Sec. 3.1 for a demonstration example).

We first consider the following one-dimensional ODE:

where $\lambda$ is a non-negative constant, $u_{0}$ is the initial condition, and $f$ is the source term. In certain real-world applications, the reaction term, i.e., $\lambda u(1-u)$ in the right hand side of Eq. (10), may not be exactly known and herein it is misspecified as an empirical model. In this regard, we study the effect of different prior knowledge about the reaction model in the discovery of governing equation. In particular, the following cases are considered:

1. 1.

   Case (A): The reaction term is correctly specified as $\lambda u(1-u)$, where $\lambda\geq 0$ is unknown.

2. 2.

   Case (B): The reaction term is misspecified as $\lambda\cos(u)$, where $\lambda\geq 0$ is unknown.

3. 3.

   Case (C): The reaction term is misspecified as $\lambda\cos(u)$ (with known $\lambda=0.2$) but we add $s(t)$ to correct it, i.e., $\lambda\cos(u)+s(t)$.

We denote the reaction term as $\phi$, and the target here is to identify $\phi$ given data on $u/f$ as well as differently specified models. In Case (A) and (B), the problem degenerates to identifying the model parameter $\lambda$, and $\phi$ in these two cases can then be computed by $\tilde{\lambda}u_{\theta}(1-u_{\theta})$ and $\tilde{\lambda}\cos(u_{\theta})$, respectively, where $\tilde{\lambda}$ is the PINN estimate of $\lambda$. In Case(C), $\phi$ is computed by $\tilde{\lambda}\cos(u_{\theta})+s_{\psi}(t)$ where $s_{\psi}$ is the estimation of the discrepancy induced by the model misspecification. We note that in Case (C) we assume $\lambda$ is known and set $\lambda=0.2$ for simplicity. This is because the discrepancy is already modeled by $s(t)$ and hence treating $\lambda$ as unknown does not make any difference in correcting the misspecified model.

To reduce the effect of the noisy/gappy data on the predicted accuracy, we first test the case with clean and sufficient data. The training data are generated by solving Eq. (10) with $f(t)=\sin(2\pi t)$ and $u_{0}=0$ using MATLAB ode45 MATLAB. We assume that we have 101 data for $u$ and $f$, uniformly sampled from $t\in[0,1]$, and we employ ensemble PINNs zou2022neuraluq, psaros2023uncertainty to learn the governing equation and also quantify the uncertainties in predictions. In all the test cases here, we choose a DNN with two hidden layers, each of which has 50 neurons and is equipped with hyperbolic tangent as the activation function, to approximate $u$. The equations are then encoded in PINNs via the automatic differentiation raissi2019physics. Further, we utilize an additional DNN to correct the model misspecification in Case (C), which has the same architecture as the one used for approximating $u$. The Adam optimizer kingma2014adam is applied to train the ensemble PINNs.

As shown in Fig. 4(a), PINNs with correct physical model (Case (A)) enable accurate inference of $\phi$ and the surrogate model $u_{\theta}$ fits the data and satisfies the ODE simultaneously. In contrast, PINNs with misspecified physical model (Case (B)) yield completely wrong estimate of $\phi$ as displayed in Fig. 4(b). In Case (C), we add a DNN to correct the misspecified reaction model used in Case (B). As illustrated in Fig. 4(c) as well as Table 1, the predicted accuracy for $\phi$ improves significantly as compared to the result in Case (B). In addition, the computational accuracy for $u$, $f$ and $\phi$ in Case (C) is quite similar as in the Case (A), where the physical model is correctly specified. The above results demonstrate the effectiveness of the present approach. We note that the predicted uncertainties in all cases are quite small, which is reasonable since we use sufficient data for both $u$ and $f$ to determine the unknown reaction term here. Furthermore, we present further results as well as the related discussion for the case with clean but gappy data in B.1.

We now keep the same setup as in Case (A) - (C) but assume that the training data are both noisy and gappy. In these cases, $50$ data of $f$ and $u$ are randomly sampled from $t\in[0,1]$ and corrupted by independent additive Gaussian noises with known noise scales ($0.05$ for $f$ and $0.01$ for $u$). As discussed in Sec. 2, we will employ the B-PINNs to quantify the uncertainties arsing from the noise as well as the incompleteness of the training data. As shown in Fig. 5, (1) B-PINNs with correct physical model can obtain accurate predictions for both $u$, $f$, and $\phi$, and the computational errors for $\phi$ are generally bounded by the predicted uncertainties (Case (A), Fig. 5(a)); (2) B-PINNs with the misspecified physical model fails to provide accurate inference for $u$, $f$, or $\phi$ (Case (B), Fig. 5(b)); and (3) B-PINNs with a DNN as correction to the misspecified model is able to improve the predicted accuracy for $u$, $f$, and $\phi$, as compared to the results in Case (B) (Case (C), Fig. 5(c)). Similar as in Case (A), the computational errors for $\phi$ are bounded by the predicted uncertainties here.

In addition to the relation between the computational errors and the predicted uncertainties for $\phi$ discussed above, we would like to further discuss that: (1) it is interesting to observe in Case (A) and Case (C) that the predicted uncertainty near the right end is larger than that near the left end. We note that the results here are reasonable since we do not have data for $f$ on the right end; (2) the predicted uncertainties of $\phi$ for $t\in[0,1]$ in Case (C) are larger than those in Case (A) although we use the same training data as well as prior for the BNN that is used to approximate $u$. We remark that the difference between the uncertainties in these two cases is caused by the lack of knowledge on the physical model and reflects the uncertainty in physical models.

We further verify the learned governing equation from the proposed approach. Specifically, we solve the equation ${du}/{dt}=f(t)+\phi(t)$ with the inferred $f$, $\phi$, and the initial condition $u_{0}=0$, using MATLAB ode45 to see if the solution for $u$ agrees with the reference solution. We denote the solution to the identified system by $\tilde{u}$. In particular, 1000 posterior samples of $f$ and $\phi$ from the B-PINNs are used to compute $\tilde{u}$, as displayed in Fig. 6. We further illustrate the relative $L_{2}$ errors between the mean for $\tilde{u}$ and the reference solution in Table  2. As we can see, PINNs with correct physical model leads to the lowest error in $\tilde{u}$. Adding a DNN to correct the misspecified model can significantly improve the accuracy in $\tilde{u}$ as compared to the case with misspecified model, and the computational errors are slightly larger than the ones from the correct physical model. We also present the computational errors for $\tilde{u}$ from the ensemble PINNs in Table  1, which will not be discussed again since the results are similar as the above ones. Finally, we would like to point out that (1) the uncertainties in the reconstructed $u$ (i.e., $\tilde{u}$) are caused by the predicted uncertainties for $f$ and $\phi$ from B-PINNs, and (2) the computational errors for the mean of $\tilde{u}$ in Case (A) and (C) are of the same order as the noise in the measurements, suggesting the reasonableness of the current results.

As discussed in Sec. 2, the discovered physical models or governing equations up to now are expressed using the trained DNNs (i.e., $s_{\psi}$), which are not explicitly presented. We can then combine the proposed method with the symbolic regression as the further step to obtain the explicit governing equations. Here we take the results in Case (C) from both ensemble PINNs and B-PINNs to demonstrate how to get the analytical expressions for the trained NNs, i.e., $s_{\psi}$. Recall that with the present methodology the reaction term is now expressed as $\phi(t)=\tilde{\lambda}\cos(u_{\theta}(t))+s_{\psi}(t)$. Assume that we now have the knowledge that the reaction term is autonomous with respect to the solution, i.e. $\phi=s(u)$. We can then adopt the symbolic regression developed in cranmer2023interpretable to find the analytical expression of $s(u)$ from $u_{\theta}$ and $s_{\psi}$. The data for the symbolic regression are generated as $\{u_{\theta}(t_{i}),s_{\psi}(t_{i})+\tilde{\lambda}\cos(u_{\theta}(t_{i}))\}_{i=1}^{N}$, where $t_{i},i=1,...,N$ are $101$ uniform grids points on $[0,1]$, and $u_{\theta}$ and $s_{\psi}$ are the predicted means from the trained ensemble PINNs/B-PINNs. As displayed in Table  3, the recovered reaction model via symbolic regression agrees well with the reference solution as we have sufficient clean data. However, for the case with noisy and gappy data, the expression for the reaction model from symbolic regression does not match the reference solution. The result is expected since we can see in Fig. 5(c) that the predicted mean for $\phi$ does not agree well with the reference solution well and the predicted uncertainties are generally large. Generally, increasing the number of training data is able to reduce the predicted uncertainty in B-PINNs psaros2023uncertainty, zou2022neuraluq. Here we also show that the learned reaction model agrees better with the reference solution as we increase the training data from 50 to 101 in B-PINNs. A recommendation here is to use the predicted uncertainty as guidance, i.e., we use the symbolic regression for discovering equation when the predicted uncertainty is small to achieve better accuracy. More results for the symbolic regression can be found in B.2.

We now consider the following time-dependent PDE that describes a reaction-diffusion system:

where $D=0.01$ denotes the diffusion coefficient, $g(u)=u(1-u)$ is the reaction model, $\lambda$ is the reaction rate, and $f(x,t)=0$ is the source term. Similar as in Sec. 3.1, we test the following three cases:

1. 1.

   Case (A): The reaction model is correctly specified and the reaction rate $\lambda\geq 0$ is unknown.

2. 2.

   Case (B): The reaction model is misspecified as $u^{2}$ and the reaction rate $\lambda\geq 0$ is unknown.

3. 3.

   Case (C): The reaction model is misspecified as $u^{2}$ and we correct it with a DNN, i.e., $\lambda u^{2}+s$ (with known $\lambda=1$).

We denote the reaction term as $\phi$ and we would like to identify $\phi$ from data of $u$, $f$ and different physical models. In Case (A)/(B)/(C), the estimate of $\phi$ is $\tilde{\lambda}u_{\theta}(1-u_{\theta})$/$\tilde{\lambda}u_{\theta}^{2}$/$\tilde{\lambda}u_{\theta}^{2}+s_{\psi}$, where $\tilde{\lambda}$ is the estimate of $\lambda$. We note that in Case (C) we assume $\lambda=1$ is known because the model misspecification induced by wrong $\lambda$ can be fully covered by $s$.

Here we only consider a general scenario, i.e., the training data are noisy and gappy. In this case, we employ B-PINNs to identify $\phi$ with uncertainties. We assume that we have 121 measurements on $u$ and $195$ measurements of $f$, which are evenly distributed in the temporal-spatial domain $(x,t)\in[0,1]\times[0,1]$. In addition, the measurement errors for both $u$ and $f$ are assumed to be Gaussian, and the noise scales are $0.02$ and $0.05$ for $u$ and $f$, respectively. Further, the reference solution is obtained by solving Eq. (11) with $g(u)=u(1-u)$, $\lambda=2$, and $f(x,t)=0$. Similar as in Sec. 3.1, we use one DNN to model the sought solution $u$ and another DNN to approximate the discrepancy $s$ in Case (C). Details of the B-PINN method, e.g. hyperparameters of NNs and HMC, can be found in A.

As show in Fig. 7, we observe that: (1) the predicted means for $u$ and $\phi$ agree the best with the reference solutions as the reaction model is correctly specified (Case (A), second column in Fig. 7); (2) the predictions for $u$ and $\phi$ show significant discrepancy as compared to the reference solution when the reaction model is misspecified (Case (B), third column in Fig. 7); (3) the predicted mean for $u$ in Case (C) is similar as that in Case (A), in which the reaction model is correctly specified; and (4) the predicted mean for $\phi$ in Case (C) is not as good as the one in Case (A), but it is much better than the result in Case (B) where the reaction model is misspecified. All the results demonstrate that the added DNNs are able to correct the misspecified model which in turn enhance the computational accuracy in PINNs.

We further illustrate the predicted uncertainties for both $u$ and $\phi$ in Fig. 8. In particular, we only depict the uncertainties at one representative time for simplicity. As we can see, (1) the computational errors for both $u$ and $\phi$ are generally bounded by the predicted uncertainties in Case (A) (first column in Fig. 8) and Case (C) (third column in Fig. 8), while they cannot be bounded by the predicted uncertainties in the case with misspecified reaction model i.e. Case (B) (second column in Fig. 8), and (2) the predicted uncertainties for $u$ and $\phi$ are generally larger in Case (C) (third column in Fig. 8) than those in Case (A) (first column in Fig. 8), which are similar as the results in Sec. 3.1. Again, we attribute the difference between the predicted uncertainties for Case (A) and (C) to the model uncertainty.

In this section, we test the capability of the present approach in non-Newtonian flows where we misspecify the constitutive relation as the Newtonian one, and then we employ DNNs to correct the misspecified models. In particular, we consider the steady flows in a 2D channel and a 2D cavity, and the employed non-Newtonian fluids are power-law ones. The constitutive relation in a power-law fluid is expressed as $\mu=\mu_{0}|\mathbf{S}|^{n-1}$, where $\mu$ is the dynamic viscosity of a fluid, $\mu_{0}$ is a constant, $\mathbf{S}=\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}$ denotes the shear rate ($|\mathbf{S}|=\sqrt{(\mathbf{S}:\mathbf{S})/2}$), $\bm{u}$ is the fluid velocify, and $n$ is the power law index. Without loss of generality, the fluids are shear-thinning (pseudoplastic, $n<1$) and shear-thickening (dilatant, $n>1$) in the channel and cavity flows, respectively.