Learning Parameters and Constitutive Relationships with Physics Informed Deep Neural Networks

Introduction

Physical models of many complex natural systems are, at best, “partially” known as conservation laws do not provide a closed system of equations. Accurate theoretical models for closing the system of conservation equations are available for homogeneous systems exhibiting time and length scale separation. Examples of accurate closures include Newtonian stress for homogeneous (Newtonian) fluids, Fick’s law for mass flux in diffusion processes, and the Darcy law for fluid flux in porous media. For more complex systems, including non-homogeneous turbulence, non-Newtonian fluid flow, multiphase flow and transport in porous media, and granular materials, accurate theoretical closures are not available. Instead, phenomenological constitutive relationships are used, which are usually accurate for a narrow range of conditions. Even when sufficiently accurate closed-form partial differential equations (PDE) models are available, (space-dependent) parameters are typically unknown.

Computational approaches for parameter and constitutive law estimation cast this inherently ill-posed problem as an optimization or a statistical inference task, usually requiring repeated evaluation of expensive forward solvers until the parameters that minimize a given error metric are found or until space of parameter configurations satisfying available data is explored. This results in significant and often intractable computational cost. Furthermore, minimization via gradient-based methods requires computing gradients from said expensive forward models, which either requires additional computational cost or careful formulation of the adjoint problems. Additionally, forward modeling requires the knowledge of initial and boundary conditions, which usually are not fully known and, therefore, also must be estimated from data together with unknown parameters and constitutive relations, significantly complicating parameter estimation.

Although significant progress has been made over the last two decades involving high-order schemes for PDEs, automatic differentiation of computer code, PDE-constrained optimization, and optimization under uncertainty, parameter estimation in large-scale problems remains a significant challenge .

While there are a number of established methods for parameter estimation in (closed-form) PDE models, such as Bayesian inference and maximum a posteriori probability (MAP) estimation , existing approaches for learning unknown physics from partially known models and data are few and not fully mature.

Unknown forms of equations at a given scale can either be found by upscaling (coarsening) known equations governing the same process at a smaller scale (e.g., ), or learned from data. In this work, we are interested in the latter.

In the most general case, system dynamics can be described as

where $F$ is a function or differential operator. A recent review of methods for learning $F$ can be found in . These methods include NARMAX , equation-free methods , Laplacian spectral analysis , and neural networks . Recently, data-driven approximations of Koopman operators, including dynamic mode decomposition , diffusion maps , delay coordinates , and neural networks , have gained significant attention. This work concerns with problems where conservation laws and other physical knowledge provide significantly more information about the structure of the operator $F(u)$. Specifically, we are interested in a steady state of the problem

where $K(\mathbf{x},u)$ is an unknown constitutive relationship, which is a function of space and the PDE state. The boundary conditions may or may not be known. Among other physical phenomena, this equation describes flow in porous media . Two special cases of this problem are $K(\mathbf{x},u)=K(\mathbf{x})$, where Eq (1) becomes a linear diffusion equation with heterogeneous diffusion coefficient $K(\mathbf{x})$; and $K(\mathbf{x},u)=K(u)$, where Eq (1) becomes a nonlinear diffusion equation with a state-dependent coefficient $K(u)$. We are interested in learning $K(\mathbf{x},u)=K(\mathbf{x})$ when measurements of both $u$ and $K$ are available, and $K(\mathbf{x},u)=K(u)$ when only $u$ measurements are available.

To achieve this objective, we propose a physics informed DNN method, where PDEs and data are used to train DNN representations of the PDE states and the unknown parameters constitutive relations. With application to Eq (1), this approach consists of defining two DNNs, one for $K(\mathbf{x},u)$ and another for $u(\mathbf{x})$, together with auxiliary DNNs obtained by substituting $K$ and $u$ into Eq (1) and using automatic differentiation to evaluate the right-hand side expression and boundary conditions. These networks are trained simultaneously using available data. Our work extends the physics informed DNN method proposed in for finding unknown constants and solutions to PDEs given states observations. It is important to note that the physics informed DNN method was originally developed for time-dependent PDEs, and the extension of this method proposed herein also can be applied to time-dependent problems.

The work is organized as follows: in Section 2, we introduce the PDE model and the physics-informed DNN approach. In Section 3, we provide a detailed study of the accuracy and performance of the proposed approach for a linear parameter estimation problem with $K(\mathbf{x},u)=K(\mathbf{x})$. In Section 4, we demonstrate the method’s accuracy for learning the non-linear constitutive relationship $K(\mathbf{x},u)=K(u)$. Discussion and conclusions are presented in Section 5.

Physics-informed Deep Neural Network Approach

subject to the Dirichlet and Neumann boundary conditions

We assume that $N_{K}$ measurements of $K$, $N_{u}$ measurements of $u$, $N_{D}$ measurements of $g$, and $N_{N}$ measurements of $q$ are collected at the locations $\{\mathbf{x}^{K}_{i}\}^{N_{K}}_{i=1}$, $\{\mathbf{x}^{u}_{i}\}^{N_{u}}_{i=1}$, $\{\mathbf{x}^{D}_{i}\}^{N_{D}}_{i=1}$, and $\{\mathbf{x}^{N}_{i}\}^{N_{N}}_{i=1}$, respectively. The observations are denoted by $K^{*}_{i}\equiv K(\mathbf{x}^{K}_{i},u(\mathbf{x}^{K}_{i}))$ ($i=1,\dots,N_{K}$), $u^{*}\equiv u(\mathbf{x}^{u}_{i})$ ($i=1,\dots,N_{u}$), $g^{*}_{i}\equiv g(\mathbf{x}^{D}_{i})$ ($i=1,\dots,N_{D}$), and $q^{*}_{i}\equiv q(\mathbf{x}^{N}_{i})$ ($i=1,\dots,N_{N}$).

To learn $K(\mathbf{x},u)$, we define the following DNNs for $u(\mathbf{x})$ and $K(\mathbf{x},u)$:

where $\theta$ and $\gamma$ are the DNN parameters. Substituting these two DNNs into the governing equation (2) and the Neumann boundary condition (4), and evaluating the corresponding spatial derivatives via automatic differentiation, yields two additional “auxiliary” DNNs:

Next, we define the following loss function to train these four networks simultaneously:

The first and second terms in $L$ force the $K$ and $u$ DNNs to match the $K$ and $u$ measurements. The third and fourth terms enforce Dirichlet and Neumann boundary conditions. Finally, the fifth term enforces the PDE at $N_{c}$ “collocation” points $\{\mathbf{x}^{c}_{i}\}^{N_{c}}_{i=1}$ that can be chosen uniformly or non-uniformly over $\Omega$ depending on the problem.

The DNNs are trained, i.e., $\theta$ and $\gamma$ are found, by minimizing the loss function:

The minimization is carried out using the L-BFGS-B method together with Xavier’s normal initialization scheme . We use a quasi-Newton optimizer such as L-BFGS-B instead of stochastic gradient descent (a more common optimizer for DNNs) because of its superior rate of convergence and more favorable computational cost for problems with a relatively small number of observed data.

Note that the proposed method does not make any assumptions about the measurement noise. Also, our method can be easily extended to time-dependent PDEs by defining DNNs in Eq (5) as functions of both $\mathbf{x}$ and $t$ .

In the following two sections, we apply the physics-informed DNNs to learn parameters and constitutive relationships in PDE models of the form (2)–(4).

Parameter estimation in a linear diffusion equation

In this section, we consider a linear diffusion equation with unknown diffusion coefficient $K(\mathbf{x})$,

subject to the Dirichlet boundary conditions

Among other problems, this equation describes saturated flow in heterogeneous porous media with hydraulic conductivity $K(\mathbf{x})$ . We assume that $N_{K}$ measurements of $K(\mathbf{x})$ and $N_{u}$ measurements of $u(\mathbf{x})$ are available: $K^{*}_{i}\equiv K(\mathbf{x}^{K}_{i})$ $(i=1,\dots,N_{K})$ and $u^{*}_{i}\equiv u(\mathbf{x}^{u}_{i})$ $(i=1,\dots,N_{u})$. We define DNNs for $K(\mathbf{x})$ and $u(\mathbf{x})$, $\hat{K}(\mathbf{x};\gamma)=\mathcal{N}\mathcal{N}_{K}(\mathbf{x};\gamma)$, and $\hat{u}(\mathbf{x};\theta)=\mathcal{N}\mathcal{N}_{u}(\mathbf{x};\theta),$ together with two auxiliary DNNs $f(\mathbf{x};\gamma,\theta)=\nabla\cdot[\mathcal{N}\mathcal{N}_{K}(\mathbf{x};\gamma)\nabla\mathcal{N}\mathcal{N}_{u}(\mathbf{x};\theta)]=\mathcal{N}\mathcal{N}_{f}(\mathbf{x};\theta,\gamma)$ and $f_{N}(\mathbf{x};\theta)=\partial\mathcal{N}\mathcal{N}_{u}(\mathbf{x};\theta)/\partial x_{2}=\mathcal{N}\mathcal{N}_{N}(\mathbf{x};\theta)$. For this problem, the loss function takes the form:

The DNNs are trained by minimizing the loss function (11) as described in Section 2. Throughout this work, we use feed-forward networks with two hidden layers and 50 units per layer.

To demonstrate the proposed approach, we generate a reference $\ln K(\mathbf{x})$ field as a realization of the Gaussian process with zero mean and covariance function $C(\mathbf{x},\mathbf{x}^{\prime})=\sigma^{2}\exp(-\|\mathbf{x}-\mathbf{x}^{\prime}\|^{2}/2\lambda^{2})$, with $\sigma=1$ and $\lambda=0.15$. The reference $u$ is generated by solving Eqs (8)–(10) using the finite volume (FV) method with the two-point flux approximation and a cell-centered regular mesh with 1024 cells. Figure 1 presents the reference $K$ and $u$ fields. We randomly choose $N_{K}$ and $N_{u}$ FV cell centroids as measurement locations for $K$ and $u$, respectively. These measurement locations are shown in Figure 2. For evaluating the loss function and training the DNNs, we use $N_{c}=1024$ uniformly distributed collocation points.

We quantify the error between estimated and reference $K$ and $u$ fields in terms of the relative $L_{2}$ errors, defined as

Figure 2 shows the estimated $\hat{K}$ and $\hat{u}$ with $N_{K}=250$, $N_{u}=100$, and $N_{c}=1024$. The relative $L_{2}$ errors are $\varepsilon_{u}\approx 0.5$% and $\varepsilon_{K}\approx 1.7$ %. Figure 2 also depicts the point-wise absolute error in estimates of $K$ and $u$. The point errors in $K$ are concentrated in the upper left corner where no $K$ measurements are available with a maximum point error of approximately $30\%$. The point errors in $u$ are much smaller (maximum error is approximately $1\%$) and more uniformly distributed throughout the domain.

Next, we study the effect of DNN initialization on the estimated $K$ and $u$. For this purpose, we draw multiple initializations of the DNNs employing Xavier’s scheme (see Section 2), and, for each initialization, we train the DNNs. We quantify the effect of initialization in terms of the mean and standard deviation of the relative $L_{2}$ errors of the estimated fields obtained for each initialization, i.e.,

where $\varepsilon_{(\cdot),i}$ is the relative $L_{2}$ error for either $K$ or $u$ for the $i$th initialization, and $N_{s}$ is the number of network initializations. Figure 3 shows the mean and standard deviation of $\varepsilon_{u}$ and $\varepsilon_{K}$ as a function of $N=N_{K}=N_{u}$ obtained from $N_{s}=11$ different network initializations. The $u$ and $K$ measurement locations are the same in these simulations, and $N_{c}=1024$ collocation points are used. As before, we see that uncertainty in $\hat{u}$ is much smaller than in $\hat{K}$, i.e., $\sigma_{\varepsilon_{K}}>>\sigma_{\varepsilon_{u}}$. For both $\hat{u}$ and $\hat{K}$, the standard deviation associated with the initialization is approximately 10 times smaller than the mean value (the coefficient of variation is $\approx 0.1$), which indicates that the initialization of the DNNs does not have a significant effect on DNN predictions.

For the DNNs training, the governing PDEs are enforced at $N_{c}$ collocation points. To study the effect of the number and location of collocation points, we compute the relative errors $\varepsilon_{u}$ and $\varepsilon_{K}$ as a function of the number and location of the collocation points. Figures 4(a) and (b) show the mean and standard deviation of the relative errors versus $N_{c}$ for $N=N_{u}=N_{K}=20$. For a given $N_{c}$, the DNNs are trained $N_{s}=11$ times for different locations chosen via Latin hypercube sampling to compute the mean and variance of the relative errors. The error in $\hat{K}$ is about eight times larger than in $\hat{u}$. As expected, the mean and standard deviation of $\varepsilon_{u}$ and $\varepsilon_{K}$ decrease with increasing $N_{c}$ until they reach asymptotic values at approximately $N_{c}=300$, which is approximately 33% of the number of grid points in these simulations. The location of collocation points has a notable effect on the errors, especially for a relatively small $N_{c}$, which is evident from relatively large coefficients of variation $\sigma_{\varepsilon_{K}}/\overline{\varepsilon}_{K}$ and $\sigma_{\varepsilon_{u}}/\overline{\varepsilon}_{u}$.

Figures 4(c) and (d) show that $\overline{\varepsilon}_{u}$ and $\overline{\varepsilon}_{K}$ asymptotically decrease for four considered values of $N$. The asymptotic values of $\overline{\varepsilon}_{u}$ and $\overline{\varepsilon}_{K}$ also decrease with increasing $N$. For all considered $N$, imposing PDE constraints reduces the mean error in $K$ and $u$ by close to 50%. Here, as in Figures 4(a) and (b), $\varepsilon_{u}$ is significantly smaller than $\varepsilon_{K}$.

In Figures 5 and 6, we study how the number of $K$ observations versus the number of $u$ observations affects $\varepsilon_{u}$ and $\varepsilon_{K}$. Here, the number of collocation points is $N_{c}=1024$. Figures 5(a) and (b) depict the effect of $N_{u}$ on the mean and standard deviation of $\varepsilon_{u}$ and $\varepsilon_{K}$ for $N_{K}=20$. For each $N_{u}$ value, we estimate $\hat{K}$ and $\hat{u}$ with $N_{s}=11$ different distributions of the $u$ measurement locations generated with the Latin hypercube sampler. Then, we compute $\overline{\varepsilon}_{u}$, $\overline{\varepsilon}_{K}$, $\sigma_{\varepsilon_{u}}$, and $\sigma_{\varepsilon_{K}}$. We see that $N_{u}$ does not significantly affect $\varepsilon_{K}$ with $\overline{\varepsilon}_{K}\approx 24$% and $\sigma_{\varepsilon_{K}}\approx 6$%. For $\hat{u}$, $\overline{\varepsilon}_{u}$ deceases from more than 3% for $N_{u}=20$ to $\approx 2$% for $N_{u}>100$. The standard deviation $\sigma_{\varepsilon_{u}}$ is 1.5% for $N_{u}=20$ and less than 1% for $N_{u}>50$.

Figures 5 (c) and (d) show $\overline{\varepsilon}_{K}$ and $\overline{\varepsilon}_{u}$ as a function of $N_{u}$ for different $N_{K}$. For all considered $N_{K}$, $\overline{\varepsilon}_{K}$ is practically independent of $N_{u}$ and decreases with increasing $N_{K}$. On the other hand, $\overline{\varepsilon}_{u}$ decreases with increasing $N_{u}$ and/or $N_{K}$.

Figures 6(a) and (b) reveal that all $\overline{\varepsilon}_{u}$, $\overline{\varepsilon}_{K}$, $\sigma_{\varepsilon_{u}}$, and $\sigma_{\varepsilon_{K}}$ decrease with increasing $N_{K}$ for fixed $N_{u}$. Alternatively, Figures 6(c) and (d) demonstrate that $N_{u}$ has a relatively minor effect on $\overline{\varepsilon}_{u}$ and $\overline{\varepsilon}_{K}$. In all considered cases, $\overline{\varepsilon}_{u}$ is almost an order of magnitude smaller than $\overline{\varepsilon}_{K}$.

The main conclusion to be drawn from Figures 5 and 6 is that $K$ measurements are more important than $u$ measurements for reducing error in $\hat{K}$ and $\hat{u}$. This may be explained partially by the fact that $u(\mathbf{x})$ is much smoother than $K(\mathbf{x})$, and a relatively small number of $u$ measurements are needed to describe the $u$ field. Beyond this number (approximately 50 for this example), additional $u$ measurements do not have a significant effect on $\overline{\varepsilon}_{u}$ and $\overline{\varepsilon}_{K}$.

For the considered problems, we find that the smallest $L_{2}$ error is obtained with $\gamma=10^{-6}$. Figures 7 and 8 compare the two methods. Figure 7 presents $\varepsilon_{K}$ as a function of $N=N_{K}=N_{u}$ for $\hat{K}$ found from the two methods. The physics informed DNNs produce $\hat{K}$ with smaller errors for all considered $N$.

Figure 8 depicts $\hat{K}$ obtained from the two methods for $N=50$. The $\hat{K}$ field estimated from the physics informed DNNs is significantly smoother (and closer to the reference $K$ shown in Figure 1) than one estimated from MAP, even though the $\varepsilon_{K}$ errors for the two methods are relatively similar: $19\%$ for NN versus $22\%$ for MAP. The tent-like character of the MAP prediction in Figure 8 stems from the discrepancy with respect to observations is penalized more than the smoothness of the estimate (result of a relatively small $\gamma$). We find that larger $\gamma$ results in a smoother field, but also in a larger prediction error.

Both the MAP and physics-informed DNN methods involve an objective-function minimization. In addition, MAP estimation via gradient-based optimization algorithms (such as the LM algorithm) requires computing the gradient of the predicted observations, $\mathbf{H}_{u}\mathbf{u}$ and $\mathbf{H}_{K}\ln\mathbf{k}$ with respect to $\mathbf{k}$. For an FV discretization, this is done via the discrete adjoint method. The total cost for each iteration of the LM algorithm is one forward solution of the PDE problem to evaluate the objective function and one adjoint solution to compute the gradient. Therefore, MAP requires careful discretization of the PDE problem and formulation and solution of the corresponding adjoint problem. In contrast, in the physics-informed DNN method, both spatial derivatives and the gradients with respect to DNN parameters are computed via automatic differentiation, and the methodology does not require solving the PDE problem or formulating an adjoint problem.

Finally, significant gains can be achieved in the physics informed DNN method performance by employing GPU accelerators.

Nonlinear diffusion equation

In this section, we consider a nonlinear diffusion equation with unknown state-dependent diffusion coefficient $K(u)$,

This equation describes a two-dimensional horizontal unsaturated flow (flow of water and air) in a homogeneous porous medium, where $u(\mathbf{x})$ is the water pressure and $K(u)$ is the pressure-dependent partial conductivity of the porous medium . In practice, $K(u)$ is difficult to measure directly. Therefore, this work assumes that no measurements of $K(u)$ are available, and only $N_{u}$ measurements of $u$ are given.

We define two DNNs for unknown $K(u)$ and $u(\mathbf{x})$,

and two auxiliary DNNs obtained by substituting the DNNs for $K$ and $u$ into (13), (15), and (16),

where $\mathbf{f}_{N}$ is a vector DNN with two DNN components $f^{(x_{1})}_{N}$ and $f^{(x_{2})}_{N}$. Then, the loss function becomes

where $\mathbf{x}_{i}^{N,x_{1}}$ $(i=1,...,N_{N}^{(x_{1})})$ are the collocation points on the Neumann boundary $(x_{1}=0,x_{2})$ and $\mathbf{x}_{i}^{N,x_{2}}$ $(i=1,...,N_{N}^{(x_{2})})$ are the collocation points on the Neumann boundaries $(x_{1},x_{2}=0)$ and $(x_{1},x_{2}=L_{2})$.

This model is tested with data generated using the Subsurface Transport Over Multiple Phases (STOMP) code with the van Genuchten model for the $K(u)$ function

Here, $K_{s}$ is the saturated hydraulic conductivity, $u_{g}=\frac{P_{g}}{\rho g}$, $P_{g}$ is the air pressure, $\rho$ is the density, $g$ is gravity, and $\alpha$ and $m$ are the van Genuchten parameters. The following parameter values are used in the STOMP simulation: $u_{0}=-10$ m, $\alpha=0.1$, $m=0.469$, $q=8.25\times 10^{-5}$ m/s, $u_{g}=0$, and $K_{s}=8.25\times 10^{-4}$ m/s.

Figure 9 shows the reference $u(\mathbf{x})$ field generated with STOMP and the assumed locations of $u$ measurements. Figure 10 presents the estimated $\hat{K}(u)$ function and the reference $K(u)$ function given by Eqs (18) and (19). It is evident that the physics informed DNN method provides an accurate estimate of unknown $K(u)$ without any direct measurements of $K$ as a function of $u$.

Finally, we examine the robustness of the physics informed DNN in the presence of observation noise. We use the exact same setup as before except we add 1% random noise to the values of observed $u$. Figure 11 shows the difference between predicted and referenced $u(x)$ and $K(u)$. The added noise increases maximum error in the reconstructed $u$ from 0.002 to 0.012, but the accuracy of the reconstructed $K(u)$ practically does not change. Note that the relative $L_{2}$ prediction errors for the “noisy” case are quite small, including $7.4\times 10^{-4}$ for $u$ and $6.4\times 10^{-3}$ for $K$. For comparison, the relative $L_{2}$ prediction errors for the “noiseless” case are $5.8\times 10^{-5}$ for $u$ and $5.9\times 10^{-3}$ for $K$.

Conclusions

In this work, we have presented a physics informed DNN method for estimating parameters and unknown physics (constitutive relationships) in PDE models. The proposed method uses both PDEs and measurements to train DNNs to approximate unknown parameters and constitutive relationships, as well as states (the PDE solution). Physical knowledge increases the accuracy of DNN training with small data sets and affords the ability to train DNNs when no direct measurements of the functions of interest are available.

We have tested this method for estimating an unknown space-dependent diffusion coefficient in a linear diffusion equation and an unknown constitutive relationship in a non-linear diffusion equation. For the parameter estimation problem, we assume that partial measurements of the coefficient and state are available and have demonstrated that the proposed method is more accurate than the state-of-the-art MAP parameter estimation method. For the non-linear diffusion PDE model with unknown constitutive relationship (state-dependent diffusion coefficient), the proposed method has proven able to accurately estimate the non-linear diffusion coefficient without any measurements of the diffusion coefficient and with measurements of the state only. We also have demonstrated that adding physics constraints in training DNNs could increase the accuracy of DNN parameter estimation by as much as 50%.

Parameter estimation is an ill-posed problem, and standard parameter estimation methods, including MAP, rely on regularization. In this work, we have trained DNNs for unknown parameters without regularizing the estimated parameter field or unknown function. In the absence of regularization, we have found that the estimates of the parameter and state mildly depend on the DNN Xavier initialization scheme. For the considered problem, the uncertainty (standard deviation) and mean error decreased with an increasing number of parameter measurements. The coefficient of variation of the relative error (the ratio of the relative error standard deviation to the mean value) was found to be approximately $0.1$. In future research, we will investigate the effect of regularization in the physics informed DNN method on parameter estimation accuracy.

Acknowledgments

This research was partially supported by the U.S. Department of Energy (DOE) Advanced Scientific Computing (ASCR) and Biological Environmental Research (BER) offices and the Pacific Northwest National Laboratory (PNNL) “Deep Learning for Scientific Discovery Agile Investment program. PNNL is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.

References