{'fraction_non_alphanumeric': 0.029983898341964403, 'fraction_numerical': 0.00983506680012185, 'mean_word_length': 4.92713438225432, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physics and Equality Constrained Artificial Neural Networks: Application to Partial Differential Equations

Physics-informed neural networks (PINNs) have been proposed to learn the solution of partial differential equations (PDE). In PINNs, the residual form of the PDE of interest and its boundary conditions are lumped into a composite objective function as an unconstrained optimization problem, which is then used to train a deep feed-forward neural network. Here, we show that this specific way of formulating the objective function is the source of severe limitations in the PINN approach when applied to different kinds of PDEs. To address these limitations, we propose a versatile framework that can tackle both inverse and forward problems. The framework is adept at multi-fidelity data fusion and can seamlessly constrain the governing physics equations with proper initial and boundary conditions. The backbone of the proposed framework is a nonlinear, equality-constrained optimization problem formulation aimed at minimizing a loss functional, where an augmented Lagrangian method (ALM) is used to formally convert a constrained-optimization problem into an unconstrained-optimization problem. We implement the ALM within a stochastic, gradient-descent type training algorithm in a way that scrupulously focuses on meeting the constraints without sacrificing other loss terms. Additionally, as a modification of the original residual layers, we propose lean residual layers in our neural network architecture to address the so-called vanishing-gradient problem. We demonstrate the efficacy and versatility of our physics- and equality-constrained deep-learning framework by applying it to learn the solutions of various multi-dimensional PDEs, including a nonlinear inverse problem from the hydrology field with multi-fidelity data fusion. The results produced with our proposed model match exact solutions very closely for all the cases considered.


Introduction
Machine learning has been highly impactful in a plethora of fields such as pattern recognition [1,2], speech recognition [3], natural language processing [4,5,6] and in the solution of partial differential equations (PDE) for forward and inverse problems. The history of machine learning has been discussed in great depths in [7,8]. It has been argued that the history of machine learning can be traced back to 300 BC, but the first primitive neural network was introduced by Warren McCulloch and Walter Pitts in 1943 [9].Years later, Rosenblatt [10] proposed the Perceptron that converges for linearly separable data. However, it was the seminal work of Hubel and Wiesel on the primary visual cortex of cats that had showed the existence of hierarchical layers of processing for visual stimulus, which has inspired modern neural networks [11]. Motivated by that innovation, the first deep neural network, termed the neocognitron, which incorporated neurophysiological insights, was introduced by Fukushima [12] in 1980, which came to be known as convolutional neural networks (CNN).
The simplest form of a neural network is called a feed-forward neural network or multilayer perceptrons, which are parametric functions that learn to map a set of inputs to a set of outputs [13]. These networks can be specialized for certain tasks, such as CNN for object recognition, recurrent neural networks (RNN) for sequential data processing [14,15,16] and generative adversarial networks, or GANs for short, and autoencoders for unsupervised learning tasks [17,18,19,20,21].
The simplest neural-network model can be a parametric linear model, such as a linear regression model for representing linear functions. The parameters of this linear model can be discovered reliably either in closed form or with convex optimization methods. Convex optimization is a subfield of mathematical optimization dealing with optimization of convex functions with convex constraints where every local minimum is a global minimum [22]. Linear models, however, have limited capacities as they are only suitable to represent linear relations between input and output data. On the other hand, deep neural networks include stacks of parametric nonlinear transformations that can be optimized for a good representation of the target nonlinear function [13]. Consequently, these nonlinearities can result in a nonconvex objective function with respect to the parameters of the network. Hence, closed-form optimization of neural networks' parameters becomes impractical. Thus, iterative gradientbased algorithms are usually preferred for optimizing the parameters of neural networks. For gradient-based learning algorithms, the network parameters are randomly initialized [23,24].
They are then adjusted to reduce the objective function. In neural networks, gradients of the objective function with respect to the parameters of the neural network are calculated using automatic differentiation [25] . Once the gradients are calculated, they can be fed to any optimization algorithm to minimize the objective function [26,27,13].
Development of the universal approximation theorem [28] has been a key milestone in the evolution of neural networks. According to this theorem a feed-forward network with as few as a single layer is sufficient to represent any function [28,29]. This attribute of neural networks has sparked an interest in the scientific community to use neural networks to learn the solution of PDEs and ordinary differential equations (ODE) in the 90s [30,31,32,33,34]. Two of the earliest examples of using neural networks to solve PDEs are the works of Dissanayake and Phan-Thien [31] and van Milligen et al. [32]. We describe the method proposed in those works in more details in Section 2.1 as the overall framework proposed in those works is common in several variations of the modern physics-constrained learning approach. Dissanayake and Phan-Thien [31] used a feed-forward neural network with two hidden layers to solve a two-dimensional inhomogeneous Poisson's equation with a sinusoidal forcing term, where they lumped the residual form of the PDE and its boundary conditions into a composite objective function for training their neural network model. van Milligen et al. [32] used a similar approach to solve a two-dimensional magnetohydrodynamic plasma equilibrium problem. Lagaris et al. [33] used neural networks for the solution of ODEs and PDEs by creating custom trial functions that satisfied the boundary conditions automatically. The authors demonstrated their method by solving two-dimensional Poisson's equations with simple exponential and sinusoidal forcing terms. However, their work is not easily extensible to non-rectangular domains as it requires construction of auxiliary functions, which is only feasible for rectangular domains. Other researchers extended the approach of Dissanayake and Phan-Thien to the solution of nonlinear Schrodinger equation [35], chemical reactor problems [36], self-gravitating systems of N bodies [37], and the solution of one-dimensional Burgers equation [38].
The overall technical approach for using neural networks to solve PDEs and ODEs that was adopted in the aforementioned works, particularly the method described in [31,32], has found a resurgent interest in recent years [39,40,41,42]. Raissi et al. [43] dubbed the term physics-informed neural networks (PINNs). They, too, formulated the loss function as a composite objective function using the residual forms of the PDE and its boundary conditions as mean squared error (MSE). The composite loss/objective function was then minimized using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [44] to train a neural network and discover its parameters. Raissi et al. made use of the modern implementation of automatic differentiation techniques [25] in the popular deep-learning platform TensorFlow [45], whereas early researchers performed automatic differentiation manually partly because their networks were composed of a few hidden layers only. Raissi et al. [46] also demonstrated the potential of the PINN approach for fluid mechanics to learn the velocity and pressure fields in arbitrary domains by encoding the Navier-Stokes equations into their neural networks using passive concentration data at sparse locations.
The physics-informed machine learning has been growing fast in popularity and applied to several unique forward and inverse problems [47,48,49,50,51,52,53,54]. Extensive reviews of the current state in physics-informed machine learning are available in literature [55,56].
One of the motivations for using neural networks in solving differential equations is their potential to break the curse of dimensionality [57,58,59]. Unlike the conventional mesh-based techniques such as finite difference, finite element, and volume methods, the meshless nature of neural networks also make them an attractive approach, albeit their slower performance compared to conventional numerical methods. Neural networks are advantageous for developing data-driven techniques for forward and inverse modeling problems as well [60].
In contrast to conventional numerical methods, once the model is trained, it can produce results at any point in the domain. In what follows, we present the PINN approach and identify a weak link in the construction of the loss function as an unconstrained optimization problem. We then formulate the machine learning of the solution of PDEs as a constrainedoptimization problem in the first place and adopt an augmented Lagrangian method (ALM) [61,62] with a unique implementation to formally convert a constrained optimization problem into an unconstrained optimization problem. Additionally, we will present a modification of the residual layers networks to address the so-called vanishing gradient problem. We present several challenging numerical examples to demonstrate the efficacy and versatility of our approach over the original PINN approach.

Technical Background
In this section, we start with describing the PINN approach and discuss some of the technical issues arising from the unconstrained optimization approach adopted in PINNs.
We then present the augmented-Lagrangian method [62,61], an approach that has enjoyed a renewed interest in constrained-optimization problems [63].

Physics-informed Neural Networks
We focus our presentation on the common elements of the physics-informed learning framework that was presented in the works of Dissanayake and Phan-Thien [31] and van Milligen et al. [32], and, as part of contemporary developments in deep learning methods, in the work of Raissi et al. [43]. Consider a scalar function u(x) : R d → R on the domain Ω ⊂ R d with its boundary ∂Ω satisfying the following partial differential equation where D and B are differential operators on the interior and boundary respectively. F and G are source functions. The above PDE and its boundary conditions are then expressed in residual forms as follows: whereû(x, θ) is an approximate solution in some underlying space V, which satisfies Eqs. (1) and (2), and θ represents its parameters. Dissanayake and Phan-Thien proposed minimizing the sum of squared error (SSE) loss expressed as a composite objective function as follows: where the weights ω Ω and ω ∂Ω were set to unity originally and the integration was approximated by sampling collocation points on Ω and ∂Ω. In van Milligen et al., these weights were tuned manually. It is worth noting that the strong form of the partial differential equation is used in Eq. (5). For this purpose, the only explicit information that is needed to feed to the network is the boundary conditions. Therefore, this approach of learning the solution of PDEs while only feeding the boundary conditions can be viewed as a form of semi-supervised learning.
Finding the proper weights that appear in the composite objective function is of utmost importance as it helps balance the contributions to the composite objective function from Eqs. (3) and (4). Parisi et al. [36] used the above approach to solve an unsteady solid-gas reactor problem composed of three coupled first-order PDEs. Parisi et al. made these weights a function of the number of collocation points where N Ω and N ∂Ω are the number collocation points in the domain and on the boundary, respectively. With the above change to the weights, the composite loss function given in Eq. (5) becomes equivalent to a sum of the individual mean of squared errors (MSE) resulting from the PDE and the boundary and initial conditions. In practice, Eq. (6) works better only because N Ω is typically much larger than N ∂Ω , which essentially gives more weight to the boundary conditions in the loss function. Based on our own experience, we find that manual tuning of these weights is not ideal, because it creates a ripple effect as we then need to tune other hyperparameters, such as the number of collocations points, the learning rate, and the architecture.
Ad-hoc nature of determining these weights undermines the predictive capability of physics-informed learning in general. Alternative approaches to determining these weights have been proposed by other researchers. E and Yu [39] proposed to express the loss function in a variational form, where trial functions were represented by the neural networks and the boundary conditions on the PDE were used as a soft-penalty term in the loss function, but the penalty term included an ad-hoc parameter that was tuned for each problem. E and Yu acknowledged the non-convex nature of the variational formulation and potential difficulties in matching the boundary conditions. This variational form of the loss function was also adopted in Zhu et al. [42]. van der Meer et al. [64] suggested a heuristic approach to approximate the weights in the loss function by casting the optimization problem as a multi-objective optimization problem with separate objectives arising from the solution of the PDE on the interior domain and satisfying the boundary conditions. The multi-objective optimization problem was then reduced to a single objective function through scalarization such that the loss function becomes a scaled version of the original form (i.e. Eq. 5). This normalized weight in the scaled loss function was then learned by a neural network using some form of the relative error. As stated by van der Meer et al., this approach is prone to learning a value from a local minimum, leading to poor predictions.
Wang et al. [65] also stressed the undesirable issues associated with manually tuning the Constrained optimization is subfield of the optimization field in general and there are well-established methods for enforcing constraints of various forms. In the present work, we consider the application of the so-called augmented-Lagrangian method (ALM), which was proposed independently by Powell [62] and Hestenes [61], to physics-constrained learning of PDEs. ALM is also called the method of multipliers and should not be confused with the method of Lagrange multipliers, although they are related. Modern applications of ALM are presented in Birgin and Martínez [63] in great details.

Augmented-Lagrangian Method for Constrained Optimization
Consider the following nonlinear, equality-constrained optimization problem with n decision variables, and m equality constraints where F is a nonlinear function of R n in R, C is a nonlinear function of R n in R m and Θ is a given subset of R n , n-dimensional Euclidean space. There are mainly three classes of methods for solving Eq. (7). The penalty method [68] converts the constrained optimization problem of Eq.(7) into an unconstrained optimization problem by adding the constraints to the objective function as quadratic penalty terms as follows: for a scalar sequence {µ k } such that µ k ≤ µ k+1 for all k and µ k → ∞. However, the main disadvantage of this method is that as µ → ∞ the penalty term prevails which induces ill-conditioning and makes the training unstable [68].
Another class of method for solving the constrained optimization problem of Eq. (7) is based on the sequential minimization of the Lagrangian function [68] defined for every θ in R n and λ in R m as follows: where one minimizes the Lagrangian L(θ, λ) over θ ∈ Θ for a sequence of multiplier vectors {λ k }. This sequence is generated iteratively where θ k is the minimizing point of L(θ, λ) over θ ∈ Θ, and α k is the stepsize (scalar) parameter. The above formula is a steepest ascent iteration aimed at finding an optimal solution of an associated dual. For this reason, it is called a primal-dual method [22]. However, the disadvantage of this method is that the problem must have a locally convex structure and the Lagrangian must be minimized a large number of times since Eq. (10) converges moderately [68]. The final class of method for solving the constrained optimization problem of Eq. (7) is the augmented Lagrangian method (ALM), which was proposed independently by Powell [62] and Hestenes [61] and is also the method of choice in the present work. In ALM, the objective function is minimized possibly by violating the constraints. Subsequently, the feasibility is restored progressively as the iterations proceed [69]. ALM combines the desirable features of the penalty and the Lagrange multipliers methods for solving the constrained optimization problem of (7). In ALM, the penalty term is added not to the objective function F(θ), but rather to its Lagrangian function, thus forming the "augmented" Lagrangian function as follows: where {µ k } is a sequence of positive penalty parameters. Minimization of Eq. (11) can be performed with any variant of gradient descent methods [70,26]. The multiplier sequence is generated from the following update rule An important aspect of the ALM is that convergence of the method may be induced by not just increasing the penalty parameter µ, but also invoking multiplier iterations as in Eq. (12). As a result, convergence in ALM may occur with finite µ k , effectively removing the ill-conditioning issue that arises in the conventional penalty method. Furthermore, the multiplier iteration (12) converges to a Lagrange multiplier vector much faster than the iteration of the primal-dual method (Eq. 10), and the optimization problem does not even have to possess a locally convex structure [68]. These aspects of the ALM make it a suitable choice for neural networks as their objective functions are typically non-convex with respect to the parameters of the network.

Proposed Method: Physics & Equality Constrained Artificial Neural Networks
Thus far, we have discussed the current state in PINNs and explained the overall framework in sufficient details to be able to elaborate on potential areas for improvement. Here, we present our contributions to the current state in physics-informed (constrained) machine learning. Our first contribution is to the way the optimization of the network parameters is formulated in the first place. We propose to cast the minimization of the residual form of the PDE as a constrained optimization problem using the boundary conditions on the PDE. To this end, we use the augmented Lagrangian method [62,61] and implement it in a unique way to make it work for constraining a PDE solution. The principled application of the ALM eliminates the need for manual tuning of weights in the loss function and produces much more improved predictions as we show in later sections. Our second contribution to physics constrained learning is addressing the infamous problem of vanishing gradients by adapting the original residual layers, which are the building blocks for residual neural networks (ResNets) [2]. We find that the straightforward application of the ResNets for learning a PDE is not optimal. Consequently, we modify the original residual layers in a unique way to make them learn PDEs with better accuracy. Together, these two contributions comprise, what we call, the physics and equality constrained artificial neural networks (PECANN, for short).

Objective Function Formulation
Our goal is to constrain the residual form of the PDE with its boundary and initial conditions and, should they exist, with multi-fidelity observed data in a principled fashion.
Therefore, our focus is on enforcing equality constraints in an optimization problem. We should stress that inequality constraints can be handled within ALM, should they exist for a problem of interest [63]. An example of enforcing inequality constraints is presented in Lu et al. [71], where they constrained an inverse design, topology optimization problem with PDE solutions using ALM.
Definition 1. Distance function Q(·) : R n → R n is a quadratic function that is acting element-wise on its input vector t such that where is an element-wise multiplication operator. Before we discuss our neural networks architecture, we construct a general objective function that can admit low-fidelity measurement data, high-fidelity measurement data, boundary conditions, initial conditions, and governing equations for the physics of the problem at hand. For notation clarification, we stress that (·) and(·) represent low fidelity and high fidelity data respectively. Given a set of low-fidelity constrained by its boundary conditions Eq. (4), we can construct the following unconstrained objective function, where µ is the penalty parameter, D denotes an In addition,D represents an N-vector noisy measurement data and the model predictions on its collocation points viewed as a column vector. Also,D represents an M-vector between high fidelity data and the model predictions on its collocation points viewed as a column vector. Finally,λ andλ denote the vectors of multipliers for the clean data and the boundary conditions, respectively. Therefore, any equality constraints, should they arise, can be incorporated in a similar fashion to the objective function L as in Eq.(13a). We can reduce our objective function into the canonical form of ALM by introducing the following parameters, The particular form of the distance function Q deserves a discussion. We use a quadratic function to calibrate the violation in our constraints. We choose this distance function intentionally for multiple reasons. To begin with, we avoid the issue of residual cancellation due to differing signs since the loss is calculated over a set of training points in machine learning applications. Otherwise, the accumulated errors might end up being either minuscule or cancel out to the point where the network has not learned anything but the user might think it has. Also, a quadratic function is differentiable. We note that one could also take the absolute value function, but this function is not differentiable and its gradient is constant with respect to its inputs. In other words, its gradient for a large error and a small error is the same. Another option would be to choose a fourth-order function or quartic function but the gradient of this function near the origin is flat which causes stiff learning. Therefore, a quadratic function is desirable for neural networks in the ALM context. Now that we have the formulation down, it is necessary to apply it under the stochastic training setting of neural networks.

Training Algorithm
There are mainly three variants of gradient descent algorithms for training neural networks [13,26]. Batch gradient descent (BGD) processes the entire training set simultaneously.
However, gradient calculations can be slow and impractical for the training sets that cannot fit in the memory [26]. On the other hand, stochastic gradient descent (SGD) technique updates the network parameters using a single training example at a time. While SGD allows using large training sets and makes frequent updates to the network parameters, the objective function heavily fluctuates as the updates are based on the gradient of a single training example. Most machine learning algorithms are known as mini-batch gradient descent methods that enjoy the merits of both BGD and SGD methods. Mini-batch gradient descents are known to be efficient in utilizing multi-core architectures [13].
In Algorithm 1, we present a training algorithm using the objective function presented in (13a). The input to the algorithm is an initialized set of parameters for the network, a maximum value µ ∞ for safeguarding the penalty term, the number of epochs E, and the training set S. We should note that over-focusing on the constraints might result in a trivial prediction, where the constraints are satisfied, but the solution has not been found. Therefore, we tackle this issue by updating the multipliers when two conditions are met simultaneously: First, the ratio of the penalty loss term from successive iterations has not decreased. Second, the maximum allowable violation on the constraints have not been met. The first condition helps prevent aggressive updating of multipliers that might cause the aforementioned issue.
In the second condition, we basically relax updating the multipliers if a satisfactory precision set by the user has been achieved. This, in turn, enables the network to freely choose to optimize any loss terms in the objective function so as to not sacrifice any loss term.

Neural Networks Architecture
In this section, we introduce our neural networks architecture that we use in both the PINN and our proposed PECANN models. We start by describing a conventional layer that has a weight layer and a nonlinear activation function as presented in Fig. 1(a). For conventional feed-forward neural networks, these nonlinear activations may saturate and cause the notoriously known vanishing-gradient problem, which makes learning significantly stiff. He et al. [2] proposed residual learning to address the performance degradation of deep convolutional networks (CNN) caused by the vanishing gradient problem. CNNs are designed for object recognition tasks. We should also note that the authors used extremely deep networks with up to 1000 layers in their study as such deep networks are not used in our current study.
He et al. used a shortcut connection to explicitly force a stack of nonlinear layers G to learn a residual mapping R := G + x, as shown in Fig.1(b). This shortcut connection perform identity mapping and their output is added to the output of the nonlinear layers, which is sent through another nonlinear σ transformation. In the original introduction of residual layers this non-linearity function was the dying relu. However, for our applications we adopt tangent hyperbolic (tanh) actitvation function because it is symmetric and many times differentiable, as the governing laws incorporated into our objective functions require taking multiple derivatives of the output of the network with respect to the input of the network. We note that sending the output R through another non-linearity after + ○ might still cause the vanishing gradient problem as pointed by the authors in [72]. Therefore, we remove this non-linearity for two purposes: First, we reduce the risk of vanishing gradient to reoccur by not having this nonlinearity function, which consequently allows direct propagation of information not just in the residual layer but throughout the entire network of residual layers. Second,it allows us to make use of the already learned features from the bottom layers and to connect each layer to every other upper layers. Having removed this last nonlinear activation function after the summation junction + ○, we observe that the weight layer before the junction becomes redundant because the output of the current residual layer will be fed to another residual layer that processes its input through a weight layer. In other words, linearly stacking two weight layers can be collapsed into a single weight layer.
Therefore, we eliminate this extra weight layer to have a modified, leaner residual layer that has no extra weight or activation function in comparison to its conventional counterpart. A schematic representation of our proposed modified residual layer is shown in Fig. 1 From Eq.(15), we observe that the gradient of loss with respect to input ∂J ∂x will not vanish even if ∂N ∂x → 0. To gain additional insight into our proposed residual neural network, we explore the behavior of a three-layer residual network with modified residual layers. Consider N 1 , N 2 ,and N 3 representing the first, the second and the third modified residual layers. We can expand the mappings from the input layer x to the output layer as follows: We observe that the third layer N 3 is not only a function N 2 , but also a function of N 1 and P(x). Consequently, we write an output representation of a modified residual network with K residual layers and one output weight layer with parameters θ L as follows: where θ L represents another vector of n−dimensional Euclidean space and T represents transposition. We observe from Eq.(17) that our network elegantly maps increasingly complex functions and aggregate them with proper weights to make a final prediction. We will demonstrate the performance gains and efficiency of our architecture with modified residual layers over another architecture with feed-forward neural networks with conventional layers as well as a network with the original residual layers on a pedagogical example presented in our results section. Fig.1(d) presents a schematic representation of a neural network architecture with K modified residual layers.

Training Details
In this section, we present the hyperparameters for efficiently training our networks.
We choose Adam [70] with its initial learning rate set to 10 −2 as our optimizer because it is adapting the learning rates for each parameter in the network during training. We use ReduceLROnPlateau(patience=100,factor=0.95) that is built-in Pytorch [73] framework to dynamically reduce the learning rate during training. We choose the hyperbolic tangent As mentioned in the preliminary section for ALM, µ ∞ can be a finite number since it does not have to approach to infinity for ALM to converge. We also justified this attribute of ALM with a pedagogical example by investigating the dependency and sensitivity of our model predictions on the limiting value of the penalty parameter µ ∞ . We also provide a heuristic approach to setting µ ∞ . First, we assume that we update our multipliers λ at every mini-batch and every epoch, for a given number of epochs E and a number of mini-batches |S|.
Also, assuming that |C i | → O(1) for all i, we show that the update rule for the multipliers follow a geometric sequence as follows: where µ ∞ corresponds to the value of 2 n in the above equation. The limiting value of our multipliers will approach the sum of the above series. Assuming |λ i | → O(10 5 ) and for all i then it is easy to see that C i → O(10 −8 ). In this work, we set the limiting penalty parameter µ ∞ = 500, unless otherwise stated. All the codes in this work are developed in PyTorch framework [73] and will be publicly available on https://github.com/HiPerSimLab/PECANN.

Numerical Experiments
We provide several PDE problems with increasing complexity to demonstrate the efficacy of our proposed PECANN model and compare its accuracy relative to the PINN approach.
Two additional examples, one for a three-dimensional problem and another for a time- In the forward model examples, we assess the accuracy of our results by providing the L ∞ and the relative L 2,r error separately on the domain and on the boundaries for both the PINN and PECANN models. Given two sets of vectorsû, the learned solution, and u, the exact solution, we define the relative L 2,r ofû as follows, We should note that the relative L 2,r error is a percentage accuracy measure of the model and it is different than the L 2 norm used in mesh-based numerical methods to assess the order of accuracy of a numerical solution. In Table 1, we present a summary of the L ∞ and L 2,r norm of the errors resulting from applying the PINN and PECANN models to all the forward modeling problems considered in this study. The errors are broken into errors generated from the PDE domain and from the boundary. Comparison of the errors resulting from the application of the PINN and the PECANN models shows that there is a significant difference in the accuracy levels between these two approaches in favor of the PECANN model, despite adopting the same neural network architectures (i.e. modified residual layers with same training details). We observe that L 2 and L ∞ error levels from the PECANN model are consistently two to three orders of magnitude lower than the error levels from the PINN model.

Pedagogical Example
The aims of this pedagogical example are twofold: First, we thoroughly demonstrate the implementation intricacies of our proposed method and highlight its advantages over the PINN approach. Second, we demonstrate the significant improvement achieved in the model prediction using our modified residual architecture relative to the original residual networks [2].  Let us consider the following one-dimensional Poisson's equation where Ω = {x | 0 ≤ x ≤ 1} and ∂Ω is its boundary. The exact solution to the above problem is a sinusoidal nonlinear function u(x) = cos(15πx). Considering a neural network solution for the above equation asû(x; θ) parameterized with θ, we write the residual form of this one-dimensional Poisson's equation as follows: Next, we use the above residual form of this differential equation to construct an objective function as proposed earlier in Eq.(13a) where µ is the penalty parameter and D is an N Ω -vector {D(x (1) ,û (1) ), · · · , D(x (N Ω ) ,û (N Ω ) )} viewed as a column vector representing the outputs of N Ω number of collocation points sampled from Ω calculated on the residual form of the differential equation as in Eq. (22).
λ ∈ R 2 is a vector of multipliers for the boundary conditions and superscript T represents transposition. Finally, Q is a distance function defined in Eq. (1). In contrast to our constrained optimization with the ALM, the composite objective function adopted in PINNs (i.e. Eq. 5) yields the following loss function for the current example Having constructed the objective functions using the constrained-optimization method in the present work and the composite approach adopted in PINNs, we design three networks in such a way that they have the same number of neurons and hidden layers to allow a fair comparison.
We   Table 2, for both the PINN and PECANN models. Clearly, the use of the modified residual layers as part of the neural network architecture improves results for both PINN and PECANN models, but the best results, by and large, are obtained when the neural network architecture with modified residual layers is used as part of the PECANN model. Based on these results, we adopt the neural network architecture with modified residual layers (see Fig. 1d) for the rest of the problems reported in the present work.
In Fig. 3, we plot the loss values during training to gain a better understanding of improved results resulting from the use of modified residual layers in connection with the PECANN model. These results highlight that our PECANN model yields 7-10 orders of magnitude lower losses for the boundary conditions while efficiently producing lower losses on the PDE in comparison to the PINN approach. This is not entirely surprising, because   aspect of the PECANN model that we observe from the same plots is that the average loss rapidly reduces at the early stages of the training, which implies that the PECANN model can produce acceptable results with a small number of epochs.
To gain further insight into the intricacy of the PECANN model during training, we plot the evolution of learning rate decay, the multipliers, and the penalty parameter during training in Fig. 4(a)-(c), respectively. From the learning rate evolution (Fig. 4a), we observe that PECANN model has a higher learning rate in the second half of training epochs than PINN's learning rate. This implies that the optimizer has reached a fairly better minimum point than that of the PINN model. In contrast, the learning rate from the PINN model keeps decreasing during the entire training. This shows that the optimizer is still searching for a better local minimum. From the penalty parameter curve in Fig. 4(c), we observe that the specified limiting penalty parameter is reached at an early stage during training. We further observe that while the penalty parameters reaches its limiting value, the multipliers still improves and plateaus much later during training as evident from Fig. 4(b). These observations substantiate the fact that both the multipliers and the penalty parameter induce convergence of ALM.

Two-dimensional Poisson's Equation
We consider a two-dimensional non-homogeneous Poisson's equation as follows: subject to the following boundary conditions where Ω = {(x, y) | − l 2 ≤ x ≤ l 2 , −l 2 ≤ y ≤ l 2 } and l = 2 is the size of the domain. We manufacture a complex solution that oscillates along one side of the domain and dies out as we move away from the boundary. The exact solution is The architecture and hyperparameters are kept the same as in section 3.4. The results obtained with the PINN and PECANN models are presented in Fig. 5. We observe that the

Two-dimensional Helmholtz Equation
Let us now consider the following two-dimensional Helmholtz equation: subject to the following boundary conditions where Ω = {(x, y) | −l 2 ≤ x ≤ l 2 , −l 2 ≤ y ≤ l 2 } and ∂Ω is its boundary and u(x, y, t) is the unknown solution to be determined. Following the equation presented above, we manufacture an oscillatory solution that satisfy the Helmholtz equation Eq. (29) with its boundary conditions Eq. (30) as follows As illustrated in Fig. 6(a)
For this problem, we keep the same hyperparameters and network architecture. The limiting penalty parameter µ ∞ is set as 1000 for this problem. For the present case, the predictions of both models for the entire domain are juxtaposed in Fig. 7.
As presented in Fig. 7(a), predicted solution from the PINN model is visibly different from the exact solution shown in Fig. 7(b). This difference is more pronounced on the boundaries as can be seen from the error distribution in Fig. 7(c). Fig. 7(d) presents the predictions from the PECANN model with much higher prediction accuracy than the PINN model due to constraining the boundary conditions in a principled fashion in the optimization problem.
This key aspect is confirmed by the uniform error distribution of our method in the domain as shown in Fig. 7(f).

Nonlinear Inverse-PDE Problem
Our PECANN model is also suitable to solve inverse-PDE problems using multi-fidelity data. With multi-fidelity, we mean that the observed data may include both data with low accuracy and data with very high accuracy. As part of our objective function formulation, we can constrain the high-fidelity data in a principled fashion and take advantage of the low-fidelity data to regularize our hypothesis space. To demonstrate our framework, we study one of the example problem that was tackled in Meng and Karniadakis [53] with composite neural networks. This particular inverse-PDE problem arises in unsaturated flows as they are central in characterizing contaminant transport [75], soil-atmosphere interaction [76], soil-plant-water interaction [77], ground-subsurface water interaction zone [78] to name a few. Describing processes involving soil-water interactions at a microscopic level is very complex due to the existence of tortuous, irregular, and interconnected pores [79]. Therefore, these flows are generally characterized in terms of their macroscopic characteristics. An important quantity that is essential in describing flows through unsaturated soil is hydraulic conductivity, which is a nonlinear parameter that is highly dependent on the geometry of the porous media [79]. Let us consider the following nonlinear differential equation representing an unsaturated one-dimensional (1D) soil column with variable water content: subject to the following boundary conditions, where is the pressure head (cm) and K(h) is the hydraulic conductivity (cm h −1 ) which is described as follows: where K s is the saturated hydraulic conductivity (cm h −1 ), and S e is the effective saturation expressed as follows [80]: where α is an empirical parameter that is inversely related to the air-entry pressure value (cm −1 ) and m is an empirical parameter related to the pore-size distribution that is hard to measure due to the complex geometry of the porous media. We aim to infer the unknown empirical parameters α, and m from sparse measurements of pressure head h. To generate multi-fidelity synthetic measurements or experimental data, we select the soil type loam for which the empirical parameters are as follows: α = 0.036 and m = 0.36.
We generate high-fidelity pressure data using the exact empirical parameters and lowfidelity data with α = 0.015 and m = 0.31. Using the built-in bvp5c MATLAB function, we solve the governing PDE as given in Eq. Eqs. 35 through Eq. 38 using the selected empirical parameters to generate multi fidelity training data as shown in Fig 8(a). In Fig 8(b) we also depict the corresponding hydraulic conductivity k(h) values for the pressure head data, which shows that low-fidelity hydraulic conductivity has a significant deviation from the exact hydraulic conductivity distribution. To highlight the robustness, efficiency, and accuracy of our framework on an inverse-PDE with multi-fidelity data fusion, we compare our results with the results reported in Meng and Karniadakis [53]. For comparison purposes, we also choose a feed-forward neural network with two hidden layers with 20 neurons per layer as in [53] for their physics-informed neural network trained on high fidelity alone which failed to discover the parameters of interest. However, Meng and Karniadakis [53] constructed customized networks for high fidelity data and low fidelity data separately and then aggregated them together by manually crafted correlations. Therefore, they refer to their approach as a composite neural networks. Unlike Meng and Karniadakis [53], we do not need to make any inductive bias about the data and, therefore, use a single network initialized with Xavier initialization technique [23] that we separately train on high-fidelity and multi-fidelity data. This shows that robustness and efficiency of our approach that we can train the same network on multi-fidelity data without the need to design customized networks to process data differently. We let a single network discover and extract features from multi-fidelity data with the help of known physics. We use Adam with 10 −2 initial learning rate and ReduceLROnPlateau(patience=500,factor=0.95) as our learning rate scheduler. We set the total number of epochs and the limiting penalty parameter µ ∞ to 1500. As for the collocation points, we uniformly generate 400 collocation points from across our domain. As considered in [53], we assume the flux at the inlet q 0 is known, which allows us to use the integral form of Eq. (35) given as follows, Fig. 8(a) and Fig. 8(b) depict the reconstructed pressure head and the corresponding hydraulic conductivity distributions obtained by our PECANN model trained on high-fidelity and multi-fidelity data separately. Comparing with the exact solution, it is seen that the inferred results are highly accurate, which shows the robustness and efficiency of our method.
Furthermore, in Table 3, we report the average and standard deviation of inferred α and m from our model along with the results from Meng and Karniadakis [53].
For this inverse problem, the average CPU training time was approximately 5.3s using a Dell OptiPlex 5050 desktop with an Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz processor and 8 GB of RAM. It is seen that the same neural network architecture trained with our proposed approach on multi-fidelity data has improved averaged values for α and m which shows the importance of low-fidelity data. However, the same model trained on just a highfidelity model has better more stable predictions since the data considered are high fidelity.
However, owing to processing the low-fidelity data, our model trained on multi-fidelity data has slightly increased CPU time compared to the same model trained on just high fidelity data.

Conclusion
Physics Our findings suggest that not only the choice of the neural network architecture, but also the optimization problem formulation is crucial in learning PDEs using artificial neural networks. We conjecture that future progress in physics-constrained (informed) learning of PDEs would come from exploring new approaches in the field of non-convex constrained optimization field. Future endeavours could shed light on challenging questions such as: how does the loss landscape of neural networks change with respect to the optimization problem formulation? What is the optimal neural network architecture for PDE learning? And, is there a physics-based approach in searching for optimal architectures? astrophysics, cosmology, classical mechanics. In this equation, we have different types of constraints that arise from physical boundary conditions and initial conditions. Here, we show that the PECANN model is applicable to equality constraints arising from both the boundary conditions and initial conditions. shows a similar behavior as before. In Table 1 we report L 2,r and L ∞ errors in the domain and on the boundary, which indicates superior predictions of our PECANN model over our PINN model. We reiterate that this improved results are obtained thanks to constraining the boundary/initial conditions using Augmented Lagrange Method and our proposed modified residual neural network architecture.