Physics-Informed Neural Networks with Trust-Region Sequential Quadratic Programming

## Abstract

Abstract Physics-Informed Neural Networks (PINNs) represent a significant advancement in Scientific Machine Learning (SciML), which integrate physical domain knowledge into an empirical loss function as soft constraints and apply existing machine learning methods to train the model.
However, recent research has noted that PINNs may fail to learn relatively complex Partial Differential Equations (PDEs). This paper addresses the failure modes of PINNs by introducing a novel, hard-constrained deep learning method — trust-region Sequential Quadratic Programming (trSQP-PINN).
In contrast to directly training the penalized soft-constrained loss as in PINNs, our method performs a linear-quadratic approximation of the hard-constrained loss, while leveraging the soft-constrained loss to adaptively adjust the trust-region radius.
We only trust our model approximations and make updates within the trust region, and such an updating manner can overcome the ill-conditioning issue of PINNs.
We also address the computational bottleneck of second-order SQP methods by employing quasi-Newton updates for second-order information, and importantly, we introduce a simple pretraining step to further enhance training efficiency of our method. We demonstrate the effectiveness of trSQP-PINN through extensive experiments. Compared to existing hard-constrained methods for PINNs, such as penalty methods and augmented Lagrangian methods, trSQP-PINN significantly improves the accuracy of the learned PDE solutions, achieving up to 1-3 orders of magnitude lower errors.
Additionally, our pretraining step is generally effective for other hard-constrained methods, and experiments have shown the robustness of our method against both problem-specific parameters and algorithm tuning parameters.

## 1 Introduction

Partial Differential Equations (PDEs) are essential mathematical tools used to model a variety of physical phenomena, including heat transfer, fluid dynamics, general relativity, and quantum mechanics. These equations are often derived from fundamental principles, such as the conservation of mass and energy, and describe how physical quantities change over space and time.
However, deriving analytical solutions of PDEs in real-world applications is usually infeasible. Numerical methods like finite difference methods (Thomas, ), finite element methods (Zienkiewicz et al., ), and pseudo-spectral methods (Fornberg, ) are commonly employed to approximate solutions. These methods, though effective, require finer meshing of the domain and can be computationally intensive. With the advent of machine learning (ML), there has been a growing interest recently in applying ML methods to solve PDEs. Deep neural networks, known for their large capacity and expressivity, are particularly useful in learning nonlinear mappings between inputs (e.g., spatial/temporal coordinates) and outputs (e.g., measurements of physical quantities), thus leading the way in such applications.
This emerging field, called Scientific Machine Learning (SciML), combines data-driven approaches (mostly neural network training) with traditional scientific modeling (i.e., PDEs/ODEs) to address the computational challenges of classical numerical methods and enhance solution accuracy (von Rueden et al., ; Willard et al., ; Huang et al., ).

The Physics-Informed Neural Network (PINN) is a recent development that integrates domain-driven principles with data-driven ML methods for learning PDE solutions (Cuomo et al., ; Hao et al., ; Raissi et al., ). The idea behind PINN is simple: it incorporates PDEs into an empirical loss function as soft constraints, and then optimizes the soft-constrained loss towards zero using classical (unconstrained) ML methods (e.g., Adam and L-BFGS). PINN is generally applied within a supervised-learning framework, utilizing paired labeled datasets to define the loss. Unlike classical PDE solvers that require meshing, PINN is mesh-free, making it well-suited for large-scale problems with irregular domains. PINN has been successfully employed to various SciML applications, including nano-optics materials (Chen et al., ), biomedicine (Sahli Costabal et al., ), fluid mechanics (Raissi et al., ), optimal control (Mowlavi & Nabi, ; Barry-Straume et al., ), and inverse design (Lu et al., ; Wiecha et al., ).

However, recent research has observed some negative results of PINN: it can fail to learn relatively complex PDEs despite employing deep neural network architectures with ample expressivity power (Krishnapriyan et al., ).
This issue is largely attributed to optimization challenges, including ill-conditioning issue and complex loss landscape, stemming from the presence of differential operators in soft penalty terms within the loss function of PINN (Krishnapriyan et al., ; Rathore et al., ; Basir, ). An imbalance between the empirical loss and the soft penalties can skew training, favoring certain constraints, and thus reducing the accuracy and generalizability of PINN. The nonconvex loss landscape also makes PINN difficult for gradient-based methods to find the optimal solution. Additionally, the complexity of neural network models and their sensitivity to variations in PDEs, coupled with the reliance on fine tuning of hyperparameters, can lead to poor convergence and physically infeasible solutions when applying PINN in practice.

To address the failure modes of PINNs, multiple strategies have recently been proposed, the majority of which attempt to alleviate the ill-conditioning issue and/or improve the loss landscape. For example, Krishnapriyan et al. () proposed two resolutions. In one approach, the authors suggested using curriculum regularization, where the penalty term starts from a simple PDE regularization and becomes progressively more complex as neural network gets trained. Another approach is to pose the problem as a sequence-to-sequence learning task (via spatio-temporal decomposition), rather than learning to predict the entire domain at once. Liu et al. () incorporated a preconditioner into the loss function to mitigate the condition number. Subramanian et al. () introduced an adaptive collection scheme that progressively collects more data points in the regions of high errors. Wang et al. (; ; ) reformulated the loss function by reweighting or resampling to balance different regularization components in the loss. Rathore et al. () developed a hybrid optimization method by integrating Newton-CG with Adam and L-BFGS to better adapt to PINN’s landscape.
While the aforementioned literature has demonstrated the effectiveness of their strategies in learning various PDEs under favorable conditions, they still primarily revolve around regularized, soft-constrained loss functions, which may not yield feasible solutions regardless of how we assign weight to the penalty term. There exist simple examples that, for any finite weight, the soft-constrained loss cannot recover the solution to the hard-constrained problem(^1^11Consider $\min x,\;\text{s.t. }x=0$. The soft penalized loss $x+\mu x^{2}$ leads to the solution $x^{\star}=-0.5/\mu$, which cannot equal 0 for any finite $\mu>0$.).

Motivated by this concern, another series of literature aims to impose PDEs as hard constraints and apply constrained optimization methods to train the model. For example, Nandwani et al. () converted soft-constrained PINN problems into alternating min-max problems and applied gradient descent ascent methods to solve them. Dener et al. () applied (stochastic) augmented Lagrangian methods to approximate the Fokker-Planck-Landau collision operator, using stochastic gradient descent to solve inner unconstrained minimization subproblems.
Subsequently, Lu et al. () proposed using penalty methods and augmented Lagrangian methods within the PINN framework for inverse design problems.
In contrast to dealing with the PINN loss directly, the above literature reverts the soft-constrained regularizations back to hard constraints (aligning with our original goals of solving PDEs), and primarily applies either penalty methods(^2^22Penalty methods differ from vanilla PINN in that they gradually increase the penalty coefficient as needed.) or augmented Lagrangian methods to solve the resulting hard-constrained problems.

Although existing hard-constrained methods based on penalty methods and augmented Lagrangian methods have shown significant improvements over PINNs, these two methods have largely been supplanted by Sequential Quadratic Programming (SQP) methods in the numerical optimization literature, which typically exhibit superior convergence properties (Gill et al., ; Nocedal & Wright, ). Compared to penalty methods, SQP does not regularize the loss in the search direction computation, so it preserves the problem structure (especially important for PDE and control problems) and does not suffer from ill-conditioning issues.
Compared to augmented Lagrangian methods, SQP is more robust to dual initialization, exhibits much faster convergence (quadratic/superlinear for SQP while linear for augmented Lagrangian), and produces solutions with fewer objective function and gradient evaluations (Curtis et al., ; Dener et al., ; Hong et al., ). These promising advantages motivate our study in this paper — design an SQP-based method to leverage recent developments of PINNs for learning PDEs.

In particular, we develop a trust-region SQP method called trSQP-PINN. Following the literature Nandwani et al. (); Dener et al. (); Lu et al. (), we first formulate the PDE problem as a hard-constrained problem, where we minimize the empirical loss subject to nonlinear PDE constraints.
Then, we apply a trust-region SQP method to impose hard constraints. At each step, we perform a quadratic approximation to the empirical loss and a linear approximation to the constraints. Inspired by the local natural of the approximation, we additionally introduce a trust-region constraint to the parameters. We only trust the approximation and make updates within the trust region.
Another critical element of the method is the merit function, which is a scalar-valued function that indicates whether a new iterate is better or worse than the current iterate, in the sense of reducing optimality and feasibility errors. We employ the soft-constrained loss as our merit function.
The trust-region radius selection and the step rejection mechanism are both based on the ratio between the actual reduction and the predicted reduction in the soft-constrained loss. If the ratio is large, indicating that our approximation is reliable, we accept the step and enlarge the radius; otherwise, we reject the step and reduce the radius.

A major difference of trSQP-PINN compared to existing penalty and augmented Lagrangian methods is that the trust-region, linear-quadratic subproblem does not rely on any penalty coefficient, effectively avoiding the ill-conditioning issue. The soft-penalized loss is only used to determine whether or not accept the new update.
Two potential computational bottlenecks when applying our method are (i) solving the trust-region subproblem, and (ii) obtaining the Hessian information of neural network models. For (i), we do not solve the trust-region subproblem exactly (which is still much easier than the ones in penalty and augmented Lagrangian methods), but instead obtain a point that is better than gradient descent (i.e. satisfying the Cauchy reduction). For (ii), we employ quasi-Newton updates such as SR1 and (damped) BFGS to approximate the Hessian information.

To further enhance training efficiency and reduce data intensity, we introduce a simple pretraining step focused solely on minimizing the feasibility error. Throughout extensive experiments, we demonstrate that trSQP-PINN significantly outperforms existing hard-constrained methods, improving the accuracy of learned PDE solutions by 1-3 orders of magnitude. Moreover, the improvement of our method is robust against both problem-specific parameters and algorithm tuning parameters, while our pretraining step is broadly effective for other methods as well.

## 2 From PINN Methods to Hard-Constrained Methods

In this section, we first introduce the PDE problem setup and the soft-constrained PINN method. Then, we provide a brief overview of two hard-constrained methods: penalty method and augmented Lagrangian method, and illustrate the motivation of designing a trust-region SQP method.

Problem setup and PINN method.
We consider an abstraction of the PDE problem with boundary conditions (BCs) and initial conditions (ICs) as follows:

$$ $\displaystyle\mathcal{F}(u(x,t))$ $\displaystyle=0,\quad(x,t)\in\Omega\times\mathcal{T}\subseteq\mathbb{R}^{d} \times\mathbb{R}^{+},$ (1a) $\displaystyle\mathcal{B}(u(x,t))$ $\displaystyle=0,\quad(x,t)\in\partial\Omega\times\mathcal{T},$ (1b) $\displaystyle\mathcal{I}(u(x,0))$ $\displaystyle=0,\quad x\in\Omega.$ (1c) $$

Here, $(x,t)$ denotes spatial-temporal coordinates with domain $\Omega\times\mathcal{T}$, $\mathcal{F}$ denotes a differential operator that can include multiple PDEs $\{\mathcal{F}_{1},\dots,\mathcal{F}_{n}\}$, $\mathcal{B}$ is a general form of a boundary-condition operator, and $\mathcal{I}$ is an initial-condition operator. In the context of PDEs, $\mathcal{F}$ may be classified into parabolic, hyperbolic, or elliptic differential operator, and its solution $u(x,t)$ models the change in physical quantities over space and time. For many problems, the analytical solution is not accessible, and solving (1) relies on classical numerical solvers such as finite differences or finite elements.

The recent data-driven approach, PINN, has become popular due to its mesh-free nature and powerful automatic differentiation techniques. In PINN, we apply neural networks to parameterize the solution

$$ $u_{\theta}(x,t)=u(x,t),$ $$

where $(x,t)$ are the NN inputs and $u_{\theta}(x,t)$ is the output, with $\theta\in\mathbb{R}^{p}$ representing the NN parameters. Let us denote the labeled observations as $\{(x_{i},t_{i},u_{i})\}_{i=1}^{N}$, so that we can define the empirical loss:

$$ $\ell(\theta)=\frac{1}{N}\sum_{i=1}^{N}\left(u_{i}-u_{\theta}(x_{i},t_{i}) \right)^{2}.$ (2) $$

To enforce constraints (1), we also randomly sample three sets of unlabeled points in the spatiotemporal domain $\Omega\times\mathcal{T}$: $\{(x_{i}^{\text{pde}},t_{i}^{\text{pde}})\}_{i=1}^{M_{\text{pde}}}$ for PDE constraints (1a), $\{(x_{i}^{\text{BC}},t_{i}^{\text{BC}})\}_{i=1}^{M_{\text{BC}}}$ for BC constraints (1b), and $\{(x_{i}^{\text{IC}},0)\}_{i=1}^{M_{\text{IC}}}$ for IC constraints (1c). Then, we define the constraint function as

$$ $\hskip-5.40608ptc(\theta)=\left(\{\mathcal{F}(u_{\theta}(x_{i}^{\text{pde}},t_ {i}^{\text{pde}}))\}_{i};\;\{\mathcal{B}(u_{\theta}(x_{i}^{\text{BC}},t_{i}^{ \text{BC}}))\}_{i};\;\{\mathcal{I}(u_{\theta}(x_{i}^{\text{IC}},0))\}_{i} \right)\in\mathbb{R}^{M_{\text{pde}}+M_{\text{BC}}+M_{\text{IC}}\eqqcolon M}.$ (3) $$

With the above setup, PINN solves Problem (1) by incorporating (3) into (2) as soft constraints:

$$ $\min_{\theta}:\;\ell(\theta)+\mu\|c(\theta)\|^{2},$ (4) $$

and applies (unconstrained) gradient-based methods to obtain the solution. Here, $\mu>0$ is the penalty coefficient that balances between the empirical loss and the constraints. When $\mu$ is large, we penalize the constraint violations severely, thereby ensuring the feasibility of the solution. However, $\mu\rightarrow\infty$ also makes the optimization problem (4) ill-conditioned and difficult to converge to a minimum (Krishnapriyan et al., ).
On the other hand, if $\mu$ is small, then the obtained solution will not satisfy the considered PDEs, making it invalid and less useful in general. Our findings suggest that PINN struggles with PDEs with high coefficients. See Appendix C for experimental results.
It is worth mentioning that a more flexible formulation is to use different penalty coefficients for PDEs $\mathcal{F}$, BCs $\mathcal{B}$, and ICs $\mathcal{I}$, while for simplicity, we unify them as $\mu$ in this paper.

Hard-constrained methods.
As introduced in Section 1, researchers have developed different hard-constrained methods to address the failure modes of PINNs, including penalty methods and augmented Lagrangian methods (Nandwani et al., ; Dener et al., ; Lu et al., ). For this type of methods, we formulate Problem (1) as a constrained nonlinear optimization problem

$$ $\min_{\theta}\;\ell(\theta)\quad\;\text{s.t. }\;c(\theta)=0,$ (5) $$

and have the following updating schemes (The pseudocodes are provided in Appendix A).

$\bullet$ Penalty methods.
Given $(\theta_{k},\mu_{k})$, we solve $\theta_{k+1}=\arg\min_{\theta}\ell(\theta)+\mu_{k}\|c(\theta)\|^{2}$ with warm initialization at $\theta_{k}$, and then update $\mu_{k+1}=\rho\mu_{k}$ with a scalar $\rho>1$. Unlike the approach of PINN, penalty methods increase the penalty coefficient $\mu$ gradually and solve the subproblem with warm initializations. However, similar to PINNs, a large $\mu$ leads to an ill-conditioned subproblem and slow convergence, and the gradient-based methods may get stuck at poor local minima. This phenomenon is observed in our experiments in Section 4 as well as in Lu et al. ().

$\bullet$ Augmented Lagrangian methods.
Given a triple $(\theta_{k},\lambda_{k},\mu_{k})$, we solve $\theta_{k+1}=\arg\min_{\theta}\ell(\theta)+\lambda_{k}^{T}c(\theta)+\mu_{k}\|c
(\theta)\|^{2}$ with warm initialization at $\theta_{k}$, and then update $\lambda_{k+1}=\lambda_{k}+\mu_{k}c(\theta_{k})$ and $\mu_{k+1}\geq\mu_{k}$. Compared to penalty methods, $\mu_{k}$ here is not necessarily increased to infinity to converge to a feasible solution. We can determine whether to increase $\mu_{k}$ based on the feasibility error $\|c(\theta_{k})\|$. For example, we let $\mu_{k+1}=\rho\mu_{k}$ if the reduction in $\|c(\theta_{k})\|$ is insufficient.

$\bullet$ Motivation of Sequential Quadratic Programming (SQP) methods.
The above hard-constrained methods, though effective for some cases, have a few limitations. First, they have complex nonlinear subproblems involving differential operators, and their quadratic term $\|c(\theta)\|^{2}$ can destroy the problem structure if any. Second, their convergence is slow. Augmented Lagrangian methods are faster than penalty methods, but only exhibit local linear convergence (Rockafellar, ). In addition, it is observed that augmented Lagrangian methods require a significant number of steps to get into the local region and are sensitive to dual initialization (Curtis et al., ).
In contrast, SQP is one of the most effective methods for both small and large constrained problems. The subproblem of SQP is an inequality-constrained linear-quadratic program, which can be efficiently solved by numerous solvers (notably, we only have to solve it approximately).
The convergence rate of SQP matches that of Newton method, and is superlinear if the Hessian is approximated by quasi-Newton update and quadratic if the Hessian is exact (Nocedal & Wright, , Chapters 17 and 18).
Most importantly, recent ML techniques, such as randomization, sketching, and subsampling, have been introduced into SQP schemes, enabling further reduction of the computational cost of this method (Hong et al., ; Na et al., ; Berahas et al., ), and making it attractive to investigate in modern SciML problems.

## 3 Trust-Region Sequential Quadratic Programming–PINN

In this section, we introduce trSQP-PINN, the first SQP-based method designed for solving PDE problems in machine learning. The method has three steps: (i) obtain a linear-quadratic approximation, (ii) incorporate a trust-region constraint, and (iii) update the iterate and trust-region radius. We also introduce a pretraining step to further enhance the training efficiency that is also effective for other methods.

Step 1: linear-quadratic approximation of Problem (5).
Let us define the Lagrangian function of (5) as $\mathcal{L}(\theta,\lambda)=\ell(\theta)+\lambda^{T}c(\theta)$. Given a primal-dual pair $(\theta_{k},\lambda_{k})$, we denote $\nabla\mathcal{L}_{k}=\nabla\mathcal{L}(\theta_{k},\lambda_{k})$ and $c_{k}=c(\theta_{k})$ (similar for $\nabla c_{k}$, $\nabla\ell_{k}$, etc.). Then, we construct a linear-quadratic approximation of the nonlinear problem (5) at $(\theta_{k},\lambda_{k})$ as:

$$ $\displaystyle\min_{\Delta\theta_{k}\in\mathbb{R}^{p}}\;\;$ $\displaystyle\nabla\ell_{k}^{T}\Delta\theta_{k}+\frac{1}{2}\Delta\theta_{k}^{T }H_{k}\Delta\theta_{k},$ (6a) s.t. $\displaystyle c_{k}+\nabla c_{k}\Delta\theta_{k}=0,$ (6b) $$

where (6a) is a quadratic approximation of the empirical loss $\ell(\theta)$ and (6b) is a linear approximation of PDE constraints $c(\theta)$. By using second-order methods, it is not surprising that the Hessian information in $H_{k}$ enables more significant per-iteration progress compared to first-order methods and better escape from stationary points. Furthermore, the computation (6) does not involve any penalty coefficient, thereby avoiding ill-conditioning issues.

However, when applied to PDE problems, the major concern is obtaining the Hessian matrix $H_{k}$. In SQP methods, $H_{k}$ corresponds to the Lagrangian Hessian $\nabla^{2}_{\theta}\mathcal{L}_{k}$, instead of the objective Hessian $\nabla^{2}\ell_{k}$, to leverage curvature information of constraints. This even complicates the computation of $H_{k}$ due to the differential operators in PDE constraints. To get rid of second-order quantities and boil down to first-order information, we perform quasi-Newton updates for $H_{k}$. We consider two schemes.

$\bullet$ Damped BFGS.
Let $s_{k}\coloneqq\theta_{k}-\theta_{k-1}$ and $y_{k}\coloneqq\nabla_{\theta}\mathcal{L}_{k}-\nabla_{\theta}\mathcal{L}_{k-1}$. Both vanilla and limited-memory BFGS methods require a critical curvature condition $s_{k}^{T}y_{k}>0$, which cannot hold in our case due to the saddle structure of the Lagrangian function. Thus, we consider a damped version of BFGS. In particular, for a scalar $\delta\in(0,1)$, we define

$$ $r_{k}=\gamma_{k}y_{k}+(1-\gamma_{k})H_{k-1}s_{k}\quad\quad\text{ with }\quad \quad\gamma_{k}=\begin{cases*}1&\text{if } $s_{k}^{T}y_{k}\geq\delta s_{k}^{T} H_{k-1}s_{k}$,\\ \frac{(1-\delta)s_{k}^{T}H_{k-1}s_{k}}{s_{k}^{T}H_{k-1}s_{k}-s_{k}^{T}y_{k}}& \text{otherwise}.\end{cases*}$ $$

Then, we compute $H_{k}$ as

$$ $H_{k}=H_{k-1}-\frac{H_{k-1}s_{k}s_{k}^{T}H_{k-1}}{s_{k}^{T}H_{k-1}s_{k}}+\frac {r_{k}r_{k}^{T}}{s_{k}^{T}r_{k}}.$ $$

When $\gamma_{k}=1$, the above updating scheme reduces to vanilla BFGS. When $\gamma_{k}=0$, $H_{k}=H_{k-1}$. We note that even if the curvature condition $s_{k}^{T}y_{k}>0$ may not hold, we ensure $s_{k}^{T}r_{k}\geq\delta s_{k}^{T}H_{k-1}s_{k}>0$, provided $H_{k-1}$ is positive definite (usually, $H_{0}=I$). Common choices of $\delta$ in numerical optimization include $\delta=0.2$ (Powell, ; Goldfarb et al., ) and $\delta=0.25$ (Wang et al., ).

$\bullet$ SR1. Despite ensuring a modified curvature condition $s_{k}^{T}r_{k}>0$, the fatal drawback of Damped BFGS is that it essentially employs a positive definite Hessian approximation ($H_{k-1}\succ 0\Rightarrow H_{k}\succ 0$). However, the exact Hessian $\nabla_{\theta}^{2}\mathcal{L}_{k}$ cannot be positive definite when training deep neural networks, leading to inaccuracies in approximation. Thus, we also perform the SR1 method, which is preferred when $\nabla_{\theta}^{2}\mathcal{L}_{k}$ lacks positive definiteness, as validated by our experiments in Appendix F. We have

$$ $H_{k}=H_{k-1}+\frac{(y_{k}-H_{k-1}s_{k})(y_{k}-H_{k-1}s_{k})^{T}}{(y_{k}-H_{k- 1}s_{k})^{T}s_{k}}.$ $$

The above scheme does not lead to $H_{k}\succ 0$ when $s_{k}^{T}y_{k}<s_{k}^{T}H_{k-1}s_{k}$.

Step 2: trust-region constraint.
The linear-quadratic approximation (6) is intended to be of only local interest. This motivates us to restrict $\theta_{k}+\Delta\theta_{k}$ to remain within a trust region centered at $\theta_{k}$; that is, we incorporate the constraint (6b) along with an additional trust-region constraint ($\Delta_{k}$ is radius)

$$ $\|\Delta\theta_{k}\|\leq\Delta_{k}.$ (7) $$

Constraint (7) offers several advantages. First, without (7), the subproblem (6) is well-defined only if $H_{k}$ is positive definite in the null space $\text{Null}(\nabla c_{k})$. This requires further regularizations of the quasi-

Newton updates that cannot preserve the curvature information.
Such regularizations are suppressed by (7).
Second, (7) serves as a hard regularization on parameters to prevent $\theta_{k}$ from moving to a worse point. As seen in Step 3, we may even skip updating $\theta_{k}$ if the new point is worse. Additionally, in PDE problems, the Jacobian $\nabla c_{k}$ tends to be singular since different (especially close) data points carry overlapping local information of PDE solutions. (7) confines the step within a small, reliable area, preventing excessive deviations from local information and mitigating singularity in Jacobians.

However, constraint (7) may conflict with constraint (6b) and lead to

$$ $\{\Delta\theta_{k}:c_{k}+\nabla c_{k}\Delta\theta_{k}=0\}\cap\{\Delta\theta_{k }:\|\Delta\theta_{k}\|\leq\Delta_{k}\}=\emptyset.$ $$

To resolve this conflict, we employ the relaxation technique originally proposed in Omojokun (). In particular, for a scalar $\nu\in(0,1)$, we adjust (6) in two steps:

$$ $\begin{aligned} \min_{\widetilde{\Delta\theta}_{k}\in\mathbb{R}^{p}}\;\;&\|c_{ k}+\nabla c_{k}\widetilde{\Delta\theta}_{k}\|^{2},\\ \text{s.t. }\;\;\;\;&\|\widetilde{\Delta\theta}_{k}\|\leq\nu\Delta_{k},\end{ aligned}\quad\quad\quad\text{ then }\quad\quad\quad\begin{aligned} \min_{ \Delta\theta_{k}\in\mathbb{R}^{p}}\;\;&\nabla\ell_{k}^{T}\Delta\theta_{k}+ \frac{1}{2}\Delta\theta_{k}^{T}H_{k}\Delta\theta_{k},\[-1.0pt] \text{s.t. }\;\;&\nabla c_{k}\Delta\theta_{k}=\nabla c_{k}\widetilde{\Delta \theta}_{k},\[2.2pt] &\|\Delta\theta_{k}\|\leq\Delta_{k}.\end{aligned}$ $$

Intuitively, the solution of the left subproblem $\widetilde{\Delta\theta}_{k}\in\text{Span}(\nabla c_{k}^{T})$ aims to satisfy the linearized constraint (6b) as much as possible within a shrunk radius of $\nu\Delta_{k}$, while the remaining radius is used to further reduce the quadratic loss (6a) by solving the right subproblem. Although closed-form solutions of these subproblems may be computable via matrix decomposition (cf. Omojokun (), (2.1.9)), we only need approximate solutions that are not worse than the Cauchy step, i.e., the steepest descent step satisfying constraint (7) (cf. Nocedal & Wright (), Chapter 4). To achieve this, we employ the dogleg method and projected conjugate gradient method to approximately solve the two subproblems.

Step 3: iterate and radius update.
With $\Delta\theta_{k}$, we then decide whether to accept the new iterate $\theta_{k}+\Delta\theta_{k}$, depending on how accurate the model problem (6) approximates the original problem (5) and how much progress the new iterate has made towards a (local) solution of (5). To do this, we leverage a non-smooth, soft-constrained loss, called the merit function:

$$ $\phi_{\mu}(\theta)=\ell(\theta)+\mu\|c(\theta)\|.$ (8) $$

Note that the empirical loss $\ell(\theta)$ alone is not suitable to justify the approximation quality since a step that decreases $\ell$ may severely violate the constraint $c$. We define the local model of $\phi_{\mu}(\theta)$ at $\theta_{k}$ as

$$ $q_{\mu}^{k}(\Delta\theta)=\ell_{k}+\nabla\ell_{k}^{T}\Delta\theta+\frac{1}{2} \Delta\theta^{T}H_{k}\Delta\theta+\mu\|c_{k}+\nabla c_{k}\Delta\theta\|.$ $$

Then, we compute the predicted reduction and actual reduction as follows:

$$ $\displaystyle\text{Pred}_{k}\coloneqq q_{\mu}^{k}(0)-q_{\mu}^{k}(\Delta\theta_ {k})\quad\quad\text{and}\quad\quad\text{Ared}_{k}\coloneqq\phi_{\mu}(\theta_{k })-\phi_{\mu}(\theta_{k}+\Delta\theta_{k}).$ $$

Finally, we calculate the ratio $\eta_{k}\coloneqq\frac{\text{Ared}_{k}}{\text{Pred}_{k}}$. For two thresholds $0<\eta_{\text{low}}<\eta_{\text{upp}}<1$, if $\eta_{k}\geq\eta_{\text{upp}}$, then it suggests that our approximation is very accurate; we perform $\theta_{k+1}=\theta_{k}+\Delta\theta_{k}$, $\Delta_{k+1}=\rho\Delta_{k}$ for $\rho>1$, and $\lambda_{k+1}=\arg\min_{\lambda}\|\nabla\ell_{k+1}+\nabla c_{k+1}^{T}\lambda\|
^{2}$. If $\eta_{\text{low}}\leq\eta_{k}<\eta_{\text{upp}}$, then our approximation is moderately accurate; we update $(\theta_{k},\lambda_{k})$ but preserve the radius $\Delta_{k+1}=\Delta_{k}$. If $\eta_{k}<\eta_{\text{low}}$, we reject the update and let $\theta_{k+1}=\theta_{k}$, $\lambda_{k+1}=\lambda_{k}$, and $\Delta_{k+1}=\Delta_{k}/\rho$.

###### Remark 1 .

We note that the soft-constrained loss (8) is only used in Step 3 to decide whether to accept the update, while it does not affect the step computation in Step 2; thus, effectively overcomes the ill-conditioning issues in penalty and augmented Lagrangian methods. In fact, the coefficient $\mu$ can be selected adaptively, instead of fine tuning it as a parameter. See Conn et al. () for more details.

Step 0: a pretraining technique.
Unlike typical ML problems, we restrict our NN parameters to satisfy PDE constraints. Random initializations often fall away from the constraint manifold, complicating training and the search for the optimal solution of the hard-constrained problem.
To address this, we develop a simple but effective pretraining step to first train the network using only the physical domain information, allowing it to initialize closer to the feasible region. We apply L-BFGS to solve

$$ $\theta_{\text{init}}\coloneqq\mathrm{argmin}_{\theta}\;\|c(\theta)\|^{2}.$ (9) $$

The pretraining step addresses the scalability concern of SQP. As observed in our experiments, we can use a different set of unlabeled data in the pretraining step (i.e., construct a different constraint function $\tilde{c}$). Then, $\theta_{\text{init}}$ will significantly reduce the number of PDE constraints required in the trSQP-PINN training. Using (9) for initialization is also beneficial for other hard-constrained methods.

We provide the pseudocode for trSQP-PINN in Appendix A.

###### Remark 1 .

## 4 Experimental Results

We implement trSQP-PINN on three PDE systems: transport equation, reaction equation, and reaction-diffusion equation. These systems may have simple analytical solutions or be solved by Fast Fourier Transform method. However, when increasing system coefficients, solving these systems becomes numerically challenging (see Krishnapriyan et al. ()). We compare our trSQP method with two popular hard-constrained methods: penalty method and augmented Lagrangian method, to demonstrate the superiority of trSQP. We also vary algorithm tuning parameters to illustrate the robustness of the trSQP improvement. Detailed experiment setups are provided in Appendix B. The results of penalty and augmented Lagrangian methods presented here are the ones coupled with the pretraining step, which significantly improve the ones without pretraining as presented in Appendix D.

### 4.1 Transport equation

We first consider a linear PDE called transport equation. This equation is commonly used to model phenomena in fluid dynamics and wave mechanics, and is applicable to scenarios like pollutant dispersion in rivers or air:

$$ $\frac{\partial u}{\partial t}+\beta\frac{\partial u}{\partial x}=0,\quad u(0,t )=u(2\pi,t),\quad u(x,0)=\sin(x),\quad(x,t)\in\Omega\times\mathcal{T}.$ $$

Here, $\beta\in\mathbb{R}$ is the transport coefficient. The second equation is a periodic boundary condition, and the third equation is an initial condition. The analytical solution is given by $u(x,t)=\sin(x-\beta t)$.

We vary the coefficient $\beta$ across a wide range of values in $\{\pm 1,\pm 10,\pm 20,\pm 30,\pm 40,\pm 50\}$. While previous experiments in Krishnapriyan et al. () have varied small $\beta$ ranging from $10^{-4}$ to $10^{-1}$, we are particularly interested in large $\beta$ where PINN fails. This exploration allows us to distinguish the performance of different hard-constrained methods in learning numerically challenging PDEs and resolving failure modes of PINN. All three methods (as well as PINN) can precisely learn the solution for $|\beta|\leq 0.1$ (cf. Krishnapriyan et al. (), Figure 1a).

We plot the absolute errors and relative errors of the three hard-constrained methods in Figure 1(a). From the figure, we observe that although penalty and augmented Lagrangian methods achieve low enough errors when $|\beta|$ is small, trSQP-PINN still decreases the errors by 0.5 to one orders of magnitude. Furthermore, trSQP-PINN distinctly outperforms the penalty and augmented Lagrangian methods when $|\beta|\geq 20$. The performance gap becomes even more pronounced with the increasing problem difficulty, where the latter two methods tend to have high errors while trSQP-PINN improves their errors by as much as one to two orders of magnitude for $|\beta|\geq 30$.

Let us examine $\beta=30$ more closely. We present the solution heatmaps for the three methods in Figure 2(a). From the figure, we clearly see that both the penalty and augmented Lagrangian methods fail to capture the complete solution features, while trSQP-PINN successfully learns the periodic dynamics. More precisely, trSQP-PINN achieves an absolute error of 0.72%, which is almost two orders of magnitude lower than the augmented Lagrangian method at 48.21% and penalty method at 51.57%.

The above results of trSQP-PINN highlight its effectiveness in navigating through complex loss landscapes to reach local optima and in mitigating notorious ill-conditioning issues. Its unique feature, the incorporation of trust-region constraints and second-order approximation, significantly enhances the performance compared to other hard-constrained methods.

### 4.2 Reaction equation

We now consider a semi-linear PDE called reaction equation. This equation is used to describe the temporal dynamics of chemical reaction concentrations and is applicable to pharmacokinetics, such as modeling the concentration of drugs in the bloodstream over time:

$$ $\frac{\partial u}{\partial t}-\alpha u(1-u)=0,\quad u(0,t)=u(2\pi,t),\quad u(x ,0)=e^{-\zeta(x-\pi)^{2}},\quad(x,t)\in\Omega\times\mathcal{T}.$ $$

Here, $\alpha\in\mathbb{R}$ is the reaction coefficient and $\zeta\in\mathbb{R}$ is the initial condition coefficient. The analytical solution is highly nonlinear and given by $u(x,t)=u(x,0)e^{\alpha t}/\left\{u(x,0)e^{\alpha t}+1-u(x,0)\right\}$.

Following the transport equation setup, we fix $\zeta=2$ and vary the reaction coefficient $\alpha$ in the set $\{\pm 1,\\
\pm 10,\pm 20,\pm 30,\pm 40,\pm 50\}$. We plot the absolute and relative errors of the three hard-constrained methods in Figure 1(b).
While all three methods maintain low errors at small values of $|\alpha|$, as the problem difficulty increases with $|\alpha|\geq 20$, the advantage of trSQP-PINN grows significantly, exhibiting errors that are one to two orders of magnitude lower compared to the other methods.
We should mention that the error gaps between trSQP-PINN and the other methods at $\alpha=-50$ to $-20$ are narrower than those at $\alpha=20$ to $50$. Nevertheless, as displayed in Appendix E, the learned trSQP-PINN solutions are still remarkably superior to the others, and the small error differences are largely attributed to the solution’s unique sharp corner point pattern — there is an extreme value at and only at $x\approx\pi$ and $t\approx 0$ (meaning most areas in the domain region cannot distinguish between different methods, and only near the point $(x\approx\pi,t\approx 0)$ can the methods be truly tested for their effectiveness).

Figure: (a) Transport equation
Refer to caption: extracted/5859249/Figures/Errors/transport.absolute_error.png

We take $\alpha=30$ as an illustrative example and draw the solution heatmaps for the three methods in Figure 2(b). From the figure, we see that trSQP-PINN excels in solving the reaction equation with high reaction coefficients. The method precisely captures the initial condition features, periodic boundary features, as well as the sharp transition patterns. In terms of errors, trSQP-PINN achieves a relative error of only 2.52%, which is nearly two orders of magnitude lower than that of the penalty method and augmented Lagrangian method, which have relative errors of 62.62% and 58.32%, respectively.

### 4.3 Reaction-diffusion equation

Finally, we consider a parabolic semi-linear PDE called reaction-diffusion equation, which adds a second-order diffusion term to the reaction equation:

$$ $\frac{\partial u}{\partial t}-\tau\frac{\partial^{2}u}{\partial x^{2}}-\alpha u (1-u)=0,\quad u(0,t)=u(2\pi,t),\quad u(x,0)=e^{-\zeta(x-\pi)^{2}},\quad(x,t) \in\Omega\times\mathcal{T}.$ $$

Here, $\alpha\in\mathbb{R}$ and $\tau>0$ are the reaction and diffusion coefficients, and $\zeta\in\mathbb{R}$ is the initial condition coefficient. The exact solution is solved by Fast Fourier Transform method (Duhamel & Vetterli, ).

We vary the reaction coefficient $\alpha$ in $\{\pm 1,\pm 10,\pm 20,\pm 30,\pm 40,\pm 50\}$ while fixing $\tau=2$ and $\zeta=2$, and vary the diffusion coefficient $\tau$ in $\{2,4,6,8,10\}$ while fixing $\alpha=20$ and $\zeta=2$. The absolute and relative errors of the three methods are plotted in Figure 1(c).
While all three methods maintain low errors at small values of $|\alpha|$, trSQP-PINN still achieves lower errors than the other two methods due to its ability to handle the high nonlinearity of the reaction-diffusion equation. Furthermore, as $|\alpha|$ increases, trSQP-PINN significantly outperforms the other two methods. Specifically, despite the increased problem complexity and ill-conditioning issues introduced by the second-order diffusion term, trSQP-PINN achieves an improvement of one to two orders of magnitude in both absolute and relative errors for $\alpha\geq 20$. For $\alpha\leq-20$, trSQP-PINN even reduces the relative error by two to three orders of magnitude.
We also observe from Figure 1(c) that the improvement of trSQP-PINN is robust across different diffusion coefficients $\tau$. For example, even when $\tau$ is as large as $10$, only trSQP-PINN achieves both absolute and relative errors below $10^{-2}$.

Figure: (a) Transport equation ($\beta=30$)
Refer to caption: x1.jpg

We plot the solution heatmaps for $\alpha=20$ and $\tau=2$ in Figure 2(c). We see that only trSQP-PINN captures the initial condition and periodic boundary features, as well as sharp transition areas. In particular, trSQP-PINN achieves an absolute error of 0.66%, which is one order of magnitude lower than the penalty method at 6.84% and the augmented Lagrangian method at 5.54%.

### 4.4 Sensitivity to Tuning Parameters

We test the sensitivity of trSQP-PINN’s performance to various tuning parameters, including the depth and width of the neural networks, the number of pretraining and training data points, and the choice of quasi-Newton updating schemes (damped BFGS v.s. SR1 v.s. identity Hessians). We also implement
penalty and augmented Lagrangian methods for comparison. The results of testing the depth and width of the neural networks and the number of pretraining data points are deferred to Appendix F, while the results of testing the number of training data points and the choice of quasi-Newton updates are summarized in Figure 3 and Table 1.

In particular, we vary the number of training data points $N$ in the set $\{100,200,500,1000\}$. From Figure 3, we observe that trSQP-PINN consistently outperforms the penalty and augmented Lagrangian methods by achieving one to two orders of magnitude lower errors. This suggests the robustness of trSQP-PINN’s superiority across different sizes of training datasets and highlights its ability to extract more PDE system information from limited data.

We employ either damped BFGS, SR1, or simply the identity matrix to approximate the Lagrangian Hessian in trSQP-PINN, with results presented in Table 1. This experiment focuses exclusively on trSQP-PINN as other methods do not involve estimating the Hessian of the Lagrangian function. As expected in Section 3, SR1 achieves lower errors when learning more complex systems, such as reaction and reaction-diffusion equations, where the exact Hessian lacks positive definiteness. Moreover, when using identity Hessians, trSQP-PINN performs worse but still better than the penalty and augmented Lagrangian methods. This indicates that (i) the use of the trust-region technique can overcome the ill-conditioning issue and better learn PDE solutions, and (ii) quasi-Newton updates to approximate second-order information can significantly boost the performance of trSQP-PINN.

Overall, we observe that the superior performance of trSQP-PINN over other methods remains robust across all tuning parameters. As seen in Appendix F, trSQP-PINN performs reasonably well even with limited pretraining data (e.g., $M=30$), indicating its ability to avoid suboptimal solutions resulting from poor initialization. Furthermore, even when reducing the depth and/or width of the neural network, thereby restricting the expressive power of the network, trSQP-PINN still performs robustly. For example, using a single-layer network with 50 neurons or a 4-layer network with 10 neurons per layer, only trSQP-PINN successfully learns the considered PDE systems.

**Table 1: Different Lagrangian Hessian approximation methods for trSQP-PINN. For each problem, the smaller error between different updating schemes is bold. TrSQP-PINN with SR1 generally achieves lower absolute and relative errors than using damped BFGS. Using identity Hessians would restricts the performance of trSQP-PINN.**
| Hessian Estimation Method | Error $(10^{-1})$ | Transport | Reaction | Reaction-diffusion |
| --- | --- | --- | --- | --- |
| Damped BFGS | Abs_err | 0.043 | 0.132 | 0.056 |
| Rel_err | 0.148 | 0.257 | 0.138 |  |
| SR1 | Abs_err | 0.072 | 0.121 | 0.026 |
| Rel_err | 0.321 | 0.252 | 0.066 |  |
| Identity | Abs_err | 2.312 | 0.172 | 0.086 |
| Rel_err | 5.663 | 0.323 | 0.206 |  |

Figure: Figure 3: Absolute and relative errors of three hard-constrained methods with varying number of training data points. TrSQP-PINN consistently outperforms the other methods and maintains low errors even if it is trained on a small dataset.
Refer to caption: extracted/5859249/Figures/robustness_experiment/Number_of_Training_Data.png

## 5 Conclusion and Future Work

We designed a direct, hard-constrained deep learning method — trust-region Sequential Quadratic Programming (trSQP-PINN) method — to address the failure modes of PINNs for solving PDE problems. At each step, the method performs a linear-quadratic approximation to the original nonlinear PDE problem, and incorporates an additional trust-region constraint to respect the local nature of the approximation. The method employs quasi-Newton updates to approximate second-order information and leverages a soft-constrained loss to adaptively update the trust-region radius. Compared to popular hard-constrained methods such as penalty methods and augmented Lagrangian methods, our extensive experiments have shown that trSQP-PINN robustly exhibits superior performance over a wide range of PDE coefficients, achieving 1-2 orders of magnitude lower errors.
The enhanced performance stems from leveraging trust-region technique and second-order information, helping trSQP-PINN to mitigate ill-conditioning issues and navigate through complex loss landscapes efficiently.

Future work includes exploring more numerical techniques for solving SciML problems and applying SQP-type methods to inverse problems and operators learning problems. For example, we can utilize Randomized Numerical Linear Algebra techniques (e.g., sketching or subsampling) to approximately solve trust-region subproblems to further reduce the computational cost of trSQP. Lu et al. () has demonstrated the effectiveness of penalty and augmented Lagrangian methods for solving inverse problems. Designing a suitable SQP-type method for such cases is also an important open direction.

## Appendix A Pseudocodes for Penalty, Augmented Lagrangian, and trSQP Methods

We define $\phi_{\mu}^{\text{P}}(\theta)=\ell(\theta)+\mu\|c(\theta)\|^{2}$ to be the soft-constrained loss for penalty method and $\phi_{\mu}^{\text{AL}}(\theta,\lambda)=\ell(\theta)+\lambda^{T}c(\theta)+\mu\|
c(\theta)\|^{2}$ to be the augmented Lagrangian function. Our stopping criterion for solving unconstrained subproblems in penalty method, augmented Lagrangian method, and the pretraining step is set as follows (we replace $\phi_{\mu}^{\text{P}}$ by $\phi_{\mu}^{\text{AL}}$ and $c$ for the latter two subproblems):

$$ $\|\nabla\phi_{\mu}^{\text{P}}(\theta^{l})\|_{\infty}\leq\text{g}_\text{tol} \quad\quad\textbf{OR}\quad\quad\|\theta^{l+1}-\theta^{l}\|\leq\text{f}_\text{ tol}\quad\quad\textbf{OR}\quad\quad l\geq l_{\max},$ (A.1) $$

where $\text{g}_\text{tol}$ and $\text{f}_\text{tol}$ are the convergence tolerances, $l_{\max}$ is the maximum number of iterations, and $l$ is the inner-loop iteration index for solving the subproblem.

Algorithms 1, 2, and 3 present the pseudocodes for penalty method, augmented Lagrangian method, and trSQP method. In the case of penalty method and augmented Lagrangian method, we adhere to the approach outlined in the literature Lu et al. (), where we terminate the algorithm when the penalty coefficient $\mu$ exceeds the threshold $\mu_{\max}$ to alleviate their ill-conditioning issues during training.

Figure: Algorithm 1 Penalty Method

Figure: Algorithm 2 Augmented Lagrangian Method

Figure: Algorithm 3 Trust-Region Sequential Quadratic Programming Method

## Appendix B Detailed Experiment Setup

For all three PDE systems (transport, reaction, reaction-diffusion), we define the spatial domain $\Omega=[0,2\pi]$ and the temporal domain $\mathcal{T}=[0,1]$. For learning transport and reaction equations, our training dataset consists of $N$ labeled points and $M=M_{\text{pde}}+M_{\text{BC}}+M_{\text{IC}}$ unlabeled points uniformly sampled from $\Omega\times\mathcal{T}$. However, due to the use of the Fast Fourier Transform in obtaining exact solutions of reaction-diffusion equation, it is time-consuming to generate solutions for each pair $(x_{i},t_{i})\in\Omega\times\mathcal{T}$.
Thus, we let $N_{\text{xgrid}}$ and $N_{\text{tgrid}}$ denote the number of evenly distributed points in $\Omega$ and $\mathcal{T}$, respectively, forming an $N_{\text{xgrid}}\times N_{\text{tgrid}}$ grid of points.
We then obtain the training dataset for learning reaction-diffusion equation by uniformly sampling from the grid points while still sampling unlabeled points from $\Omega\times\mathcal{T}$.
The $N_{\text{xgrid}}\times N_{\text{tgrid}}$ grid is also used for evaluating all the methods. We use $M^{\text{pretrain}}$ to denote the number of unlabeled points for the pretraining phase, which is evenly distributed among $M_{\text{pde}}$, $M_{\text{BC}}$, and $M_{\text{IC}}$. Analogously, we use $M^{\text{train}}$ points for the training phase.

Given a spatial-temporal pair $(x_{i},t_{i})$ on the grid, the observation (i.e., label) is constructed as $u_{i}=u(x_{i},t_{i})+\epsilon_{i}$, where $u(\cdot,\cdot)$ is either an analytical PDE solution (transport, reaction) or a solution approximation derived by Fast Fourier Transform (reaction-diffusion), and $\epsilon_{i}$ is a random noise term added to enhance generalization.
The prediction accuracy is assessed using absolute errors and relative errors, defined as

$$ Abs_err $\displaystyle\coloneqq\frac{1}{N_{\text{xgrid}}\cdot N_{\text{tgrid}}}\sum_{i= 1}^{N_{\text{xgrid}}}\sum_{j=1}^{N_{\text{tgrid}}}\|u_{\theta}(x_{i},t_{j})-u( x_{i},t_{j})\|_{2},$ Rel_err $\displaystyle\coloneqq\frac{1}{N_{\text{xgrid}}\cdot N_{\text{tgrid}}}\sum_{i= 1}^{N_{\text{xgrid}}}\sum_{j=1}^{N_{\text{tgrid}}}\frac{\|u_{\theta}(x_{i},t_{ j})-u(x_{i},t_{j})\|_{2}}{\|u(x_{i},t_{j})\|_{2}}.$ $$

We employ a 4-layer, fully-connected neural network with $50$ neurons in each layer. The activation function is $\tanh$. The subproblems in penalty and augmented Lagrangian methods are solved by L-BFGS coupled with backtracking line search. We utilize Flax to build the neural network and Jax to perform automatic differentiation.

For both penalty method and augmented Lagrangian method, the initial penalty coefficient $\mu_{0}$ is set to 1 with an increasing factor $\rho=1.1$. The initial dual vector is set to $\lambda_{0}=0$. For trSQP-PINN, the initial penalty coefficient $\mu_{-1}$ is set to 1 and adaptively selected at each step. The initial trust-region radius $\Delta_{0}$ is set to 1. We apply SR1 for quasi-Newton update in the method.

We employ a fine grid for evaluation, with $N_{\text{xgrid}}=2560$, $N_{\text{tgrid}}=1000$, and $M^{\text{pretrain}}=150$. We set $M^{\text{train}}=12$ for transport equation and $M^{\text{train}}=7$ for both reaction and reaction-diffusion equations.
The use of different numbers of unlabeled points for learning different PDEs is largely due to varying properties of the PDEs; we prefer fewer unlabeled data points when the PDEs become more complex and make the loss landscape more tortuous.
For all three PDE systems, we set $N=1000$. We have tested the robustness of our method against $N$ and $M^{\text{pretrain}}$ in Appendix F.

## Appendix C Failure Modes of PINN Method

In this section, we implement vanilla PINN method outlined in Section 2 (see (4)) with varying penalty coefficient $\mu$. We illustrate that PINN fails to recover the solutions of all three systems with large system coefficients.

We set $\beta=30$ for transport equation, $\alpha=30$ and $\zeta=2$ for reaction equation, and $\alpha=20$, $\tau=2$, and $\zeta=2$ for reaction-diffusion equation. Our results are summarized in Table 2. Comparing this table with Figure 1, we observe that PINN performs much worse than the three hard-constrained methods. Moreover, as the coefficient $\mu$ increases, the error explodes due to the worsening ill-conditioning issue. To visualize the error table, we take the transport equation as an example and draw the heatmap of the PINN solution in Figure 4. The heatmap of the exact solution is shown in Figure 2(a). From these two figures, we clearly see that PINN with varying coefficient $\mu$ fails to learn the solution of the transport equation. We have similar observations for the reaction and reaction-diffusion equations. This observation is consistent with the recent study (Krishnapriyan et al., ).

**Table 2: Absolute and relative errors of PINN method with varying penalty coefficients. PINN exhibits high prediction errors, which increase as the penalty coefficient $\mu$ increases, indicating worsening ill-conditioning issues. The errors in predicting the solutions of the transport equation peak at $\mu=1000$, while PINN consistently performs poorly in predicting the solutions of the reaction and reaction-diffusion equations due to severe ill-conditioning already occurring at $\mu=1$.**
| Coefficient | Error $(10^{-1})$ | Transport | Reaction | Reaction-diffusion |
| --- | --- | --- | --- | --- |
| $\mu=1$ | Abs_err | 5.746 | 7.806 | 8.954 |
| Rel_err | 9.403 | 9.999 | 9.999 |  |
| $\mu=10$ | Abs_err | 5.912 | 7.806 | 8.954 |
| Rel_err | 9.675 | 10.000 | 10.000 |  |
| $\mu=100$ | Abs_err | 7.655 | 7.806 | 8.954 |
| Rel_err | 13.808 | 10.000 | 10.000 |  |
| $\mu=1000$ | Abs_err | 11.329 | 7.806 | 8.953 |
| Rel_err | 19.176 | 9.999 | 9.999 |  |

Figure: (a) PINN ($\mu=1$)
Refer to caption: x13.jpg

## Appendix D Pretraining v.s. No Pretraining

In this section, we demonstrate that our pretraining step is generally beneficial for both penalty method and augmented Lagrangian method. To this end, we compare the results of these two methods with and without pretraining. For the methods without pretraining, we employ random initializations. The coefficients of the three PDE problems are set as described in Appendix C.

Our results are summarized in Table 3. From the table, it is evident that both penalty and augmented Lagrangian methods without pretraining perform significantly worse than their pretrained counterparts when solving reaction and reaction-diffusion equations. While the unpretrained penalty method may exhibit similar performance to its pretrained counterpart for solving transport equation, we find that the unpretrained augmented Lagrangian method is markedly less effective than the pretrained one. Overall, our pretraining step proves to be effective for hard-constrained methods, improving the prediction error by as much as an order of magnitude.

In Figure 5, we visualize the error table by drawing the solution heatmaps of penalty and augmented Lagrangian methods for learning reaction and reaction-diffusion equations. From the figure, we observe that unpretrained methods fail to capture any solution features, while pretrained methods perform reasonably better (especially for augmented Lagrangian). Thus, we conclude that our pretraining phase sets a robust baseline for further optimization and enhances the efficacy of hard-constrained methods in finding better solutions.

**Table 3: Absolute and relative errors of penalty and augmented Lagrangian methods with/without pretraining. For each method, the smaller error between pretraining and unpretraining is bold. The unpretrained versions of penalty and augmented Lagrangian methods generally perform much worse than their pretrained counterparts in learning three PDEs. One exception is about the penalty method for solving transport equation. The pretrained and unpretrained methods show similar results.**
| PDE | Error $(10^{-1})$ | Penalty | Augmented Lagrangian |  |  |
| --- | --- | --- | --- | --- | --- |
|  |  | Unpretrained | Pretrained | Unpretrained | Pretrained |
| Transport | Abs_err | 5.264 | 5.157 | 9.860 | 4.821 |
| Rel_err | 8.986 | 8.746 | 17.313 | 8.479 |  |
| Reaction | Abs_err | 8.953 | 4.339 | 7.808 | 3.766 |
| Rel_err | 9.999 | 6.262 | 10.003 | 5.832 |  |
| Reaction-diffusion | Abs_err | 8.953 | 0.364 | 9.754 | 0.300 |
| Rel_err | 9.999 | 0.684 | 12.495 | 0.554 |  |

Figure: (1a) Unpretrained v.s. Pretrained penalty
Refer to caption: x17.jpg

Figure: (1a) Penalty method
Refer to caption: x25.jpg

## Appendix E Additional Results for Reaction Equation

In this section, we present additional solution heatmaps for learning the reaction equation with $\alpha\in\{-20,-30,-40,-50\}$ in Figure 6. As explained in Section 4.2, although the errors among the three methods are close, Figure 6 reveals significant differences in their learned solutions. Learning the reaction equation with a negative coefficient is extremely difficult near the point $(x\approx\pi,t\approx 0)$, and only trSQP-PINN can reasonably recover the unique sharp corner point pattern for $\alpha\leq-20$. The small differences in errors are also attributed to such an unique pattern in the solutions; most areas in the domain region cannot distinguish between different methods, and only near the point $(x\approx\pi,t\approx 0)$ can the methods be truly tested for their effectiveness

## Appendix F Additional Results for Sensitivity Experiment

In conjunction with the sensitivity experiments described in Section 4.4, we present in this section additional sensitivity experiments to test the robustness of trSQP-PINN’s performance against some other tuning parameters, including the depth and width of the neural networks and the number of pretraining data points. We also include penalty and augmented Lagrangian methods for comparison. For each experiment, only one parameter is altered from the default settings detailed in Appendix B. The problem coefficients are provided in Appendix C; that is, $\beta=30$ for transport equation, $\alpha=30$ and $\zeta=2$ for reaction equation, and $\alpha=20,\tau=2,$ and $\zeta=2$ for reaction-diffusion equation.

Figure 7 shows the sensitivity results for varying the number of pretraining data points. We vary $M^{\text{pretrain}}$ in the sets $\{30,45,150,300\}$. As discussed in Section 3, having fewer pretraining data results in initializations far from the feasible region. The figure demonstrates that trSQP-PINN is significantly more robust to such initializations compared to the other hard-constrained methods.

Figures 8 and 9 show the sensitivity results for varying the depth and width of the neural networks.
We vary the depth in the set $\{1,2,3,4\}$ (while fixing the width $50$) and vary the width in the set $\{10,\\
20,30,40,50\}$ (while fixing the depth $4$). Our results suggest that only trSQP-PINN can consistently learn PDE solutions with relatively small networks. From Figure 8, we see that trSQP-PINN achieves lowest errors when the depth decreases to one layer across all three problems; however, we also see that the errors double compared to using 4-layer networks. This increase in errors can be attributed to the limited expressive power of these shallow networks.
From Figure 9, trSQP-PINN also consistently outperforms the other methods across all neural network widths. However, when solving the transport and reaction equations, all methods exhibit high errors with only 10 neurons per layer, suggesting low expressivity of the neural networks. When solving the reaction-diffusion equation, trSQP-PINN is able to maintain low errors across all neural network widths, while the penalty and augmented Lagrangian methods are highly unstable, showing very high errors when the neural network width is 20 or 30.

Figure: Figure 7: Absolute and relative errors of three hard-constrained methods with varying number of pretraining data points. TrSQP-PINN consistently outperforms the other methods and maintains low errors even if it is initialized far from the feasible region.
Refer to caption: extracted/5859249/Figures/robustness_experiment/Number_of_Pretrain_Data.png

Figure: Figure 8: Absolute and relative errors of three hard-constrained methods with varying neural networks depth. TrSQP-PINN consistently outperforms the other methods and maintains low errors even when trained on shallower neural networks. However, while trSQP-PINN achieves the lowest errors on shallow networks, the errors double compared to those on 4-layer networks. This suggests that sufficient network expressivity is necessary for better learning of PDEs.
Refer to caption: extracted/5859249/Figures/robustness_experiment/NN_Depth.png

Figure: Figure 9: Absolute and relative errors of three hard-constrained methods with varying neural networks width. TrSQP-PINN consistently outperforms the other methods and maintains low errors even when trained on narrower neural networks. All methods exhibit high errors when learning transport and reaction equations with only 10 neurons per layer, suggesting that the network may lack sufficient expressive power to capture the PDE solution in this case. When learning the reaction-diffusion equation, the penalty and augmented Lagrangian methods perform comparably to trSQP-PINN at neural network widths of 10 or 40. However, these methods are highly unstable, exhibiting very high errors at widths of 20 or 30, whereas trSQP-PINN maintains low errors across all widths.
Refer to caption: extracted/5859249/Figures/robustness_experiment/NN_Width.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: Listed in weakness
soundness: 2
strengths: 1. The idea of using constraint optimization to solve PINN is interesting.
2. The authors gives a detailed description of their method, making it easy to follow.
summary: This paper focuses on the optimization problems in PINN methods. The authors introduce the trust-region sequential quadratic programming to solve the constraint in PINN, e.g., initial condition and PDE constraint. They describe the method with mathematical formulas and pseudo codes. The authors conduct experiments on the well-known failure cases of vanilla PINN. The results show that the proposed method can solve these failure cases of PINN.
weaknesses: 1. The primary concern from the reviewer is the problem setting. According to Eq.(2), (3), the authors regard the empirical loss as loss function and regard initial/boundary condition as the constraint, so they includes several labeled observations in the training to compute the loss function. However, this setting is significantly different from the failure cases in Ref[1], in which the boundary condition and initial condition are the only available observations. Given that the observations in the vanilla PINN are viewed in the same way, i.e., using L2 loss to minimize the difference, it's hard to believe that one can observe failure cases in Ref[1] within this setting. The reviewer wonders why the authors can observe the same failure case as in Ref[1]. Also, the reviewer wonders if these failure cases can be solved by simply reweighting the loss term. If so, the experiments in current manuscript may not be suitable to demonstrate the effectiveness of the proposed method. 
2. Another concern from the reviewer is the performance of the baseline penalty method. According to the footnote from the authors, the only difference  between penalty method and vanilla PINN is the magnitude of penalty. In this paper, the authors also introduce a 'pretraining method',  which tries to minimize the L2 loss of constraint term. This pretraining method, which can also be viewed as a standard PINN training without bulk observations, is shown effective even for the penalty method. However, according to the footnote from the authors, this pretraining method would be very similar to the penalty method and should not lead to significant performance difference. Why does the author can observe a large difference with/without the pretraining method in the penalty method?
3. The final concern from the reviewer is about the significance of the proposed method. The authors proposed a new optimization method to deal with the constraint in PINN. However, in Ref[2], the results from an Adam+L-BFGS optimizer is comparable with the results listed in this paper. Considering about the large number of hyperparameters introduced by the proposed method, simply using adam+L-BFGS may be a better choice. The reviewer suggest that the authors should give more solid evidence that the proposed method outperforms original methods, e.g., Adam+L-BFGS.

[1] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems, 34:26548–26560, 2021.

[2] Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in training pinns: A loss landscape perspective. arXiv preprint arXiv:2402.01868, abs/2402.01868, 2024.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. Problem formulation
This paper considered the forward problem of partial differential equations (PDEs), which is well-posed when given appropriate initial and boundary conditions. However, the authors have formulated the problem as:
\begin{equation}\min_{\theta} \mathcal{l}(\theta) + \mu \Vert c(\theta)\Vert^2, \end{equation}
where $\mathcal{l}$ represents an "empirical" loss term, used to denote additional data, and $c(\theta)$ includes the governing equations, initial conditions, and boundary conditions. Typically, the loss function $\mathcal{l}(\theta)$ is used in inverse problems, where information is missing within the PDE framework. In the most relevant reference, [1], $\mathcal{l}(\theta)$ serves as the loss function specifically for initial and boundary conditions. Consequently, I believe that in this paper's forward problem context, initial and boundary conditions would be more appropriate choices for the hard constraints rather than $\mathcal{l}(\theta)$.

2. Experiments
The authors provide a detailed examination of PINNs’ failure modes, drawing inspiration from prior works such as [1], particularly regarding the transport equation. However, the experiments are limited to comparisons with the penalty method and the augmented Lagrangian method. To better assess the proposed method, a broader comparison with recently proposed methods, including [1,2], would be valuable. Additionally, the experiments are confined to simple, one-dimensional problems. Extending the numerical validation to more complex cases would further support the robustness of the findings.

3. Theoretical contributions
Both baseline algorithms—the penalty and augmented Lagrangian methods—come with convergence guarantees. For example, see [3] for convergence proof on the augmented Lagrangian method. Could similar theoretical contributions be established for the proposed method?

4. References
Key baseline algorithms are not cited in the manuscript (such as [3, 4] and several others proposing the augmented Lagrangian method for PINNs). Moreover, [5] did not actually consider PINNs, contrary to what the authors mentioned. It would be beneficial for the authors to carefully review their references.

[1] A. S. Krishnapriyan, A. Gholami, S. Zhe, R. M. Kirby, and M. W. Mahoney, Characterizing possible failure modes in physics-informed neural networks. NeurIPS. 2021.
[2] S. Wang, S. Sankaran, and P. Perdikaris, Respecting causality for training physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering. 2024.
[3] H. Son, S. W. Cho, and H. J. Hwang, Enhanced physics-informed neural networks with augmented Lagrangian relaxation method (AL-PINNs). Neurocomputing. 2023.
[4] J. Huang, H. Wang, and Tao Zhou, An Augmented Lagrangian Deep Learning Method for Variational Problems with Essential Boundary Conditions. Communications in Computational Physics. 2022.
[5] Y. Nandwani, A. Pathak, Mausam and P. Singla, A Primal-Dual Formulation for Deep Learning with Constraints. NeurIPS. 2019.
soundness: 2
strengths: - Successfully implemented hard constraints within PINNs.
- Mitigated the ill-conditioning issue by alleviating the penalty parameter.
- Numerical experiments show that the proposed method outperforms traditional approaches.
summary: This paper proposed a novel training method for PINNs that aims to address the well-known failure modes of PINNs. They utilized the trust-region sequential quadratic programming to approximate the hard constraints in PINN formulation. The proposed method is assessed through numerical examples compared with the penalty method and the augmented Lagrangian method.
weaknesses: - The problem formulation for PINNs appears unconventional.
- The experiments are insufficient to fully support the authors' claims.
- The paper lacks substantial theoretical contributions.
- There are errors and omissions in the references.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 4
questions: -	Experiments: 
   o	How long have you trained your method? What are the inference time? 
   o	How doest this method behaves/scales on 2d problems? 
   o	Could the author explain with more details Figure 1? More specifically, why the method performs sometimes better on negative parameters compared to positives? Intuitively, with equivalent magnitude, one would expect the performances to be as bad in both directions (if the effect of the coefficient is indeed symmetric, however, no comment is made about this aspect in the experiments). 
-	Why one should use hard constrained method in the PINNs framework? I’ve noted the toy example in the footnote which is well-informative, however is this case being reached with PINNs? 
-	Why does only computing an approximation of the problem is sufficient? Are the theoretical guarantees kept?

Other remarks
-	Preconditionning: PINNs are known to be ill-conditioned even with small $\mu$ coefficients, ie the differential operator in itself is ill-conditioned. This has been shown in several works such as [2, 4] and many others. Weight schedules are proposed for example in [5]. 

[4] Tim De Ryck, Florent Bonnet, Siddhartha Mishra, and Emmanuel de Bézenac. An operator precondi-
tioning perspective on training in physics-informed machine learning, 2023.

[5] Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent
kernel perspective. Journal of Computational Physics
soundness: 2
strengths: - The paper is well written despite the numerous steps of the algorithm: several choices, approximation or constraints have been considered in the paper. 
- The algorithm is well detailed and justified. 
- Moreover, in the experiments, several ablations are conducted to show the robustness of the proposed method.
summary: The proposed paper study the use of SQP method in the context of PINN. In a first section, the authors introduce the algorithm and detail the choices they’ve made to reach a functional algorithm for PINNs. Finally, they show the performances of their algorithm on several benchmarks.
weaknesses: -	The proposed contributions are very specific to the PINNs case. In my opinion, it is difficult to understand in which context this method could be used in practice. 
-	Misunderstanding of PINNs: PINNs designed commonly the work proposed by Raissi [1] which is completely data-free. In the paper, the author stated that PINNs are data driven which might be confusing for the reader, whether they designate [1] or Physics-informed method in general. Physics-Informed Neural Network do not rely on any data and optimize the solution neural network through the residual of the PDE and the boundary loss only (eq. 3 only in the paper). 
-	Missing baselines: Since this work targets PINNs, I think more baselines from this literature should be added, first one being [1]. Other example could include soft-constrained methods, in order to highlight the interest of hard-constrained method w.r.t soft-constrained. This is unclear to me why one should take this entire algorithm if easier to implement method could be used. Moreover, since in this method data are used, data-based method should be used as a comparison or at least hybrid method. As an example, the method of [3] could be used as a base model to evaluate the performances of the proposed method, since the section 4.1 take this work as a starting point. 
-	Missing experiments: The author claimed an improved training time w.r.t other method, however no computational time comparison nor number of steps to convergence is shown in the experiments. 
-	Pretraining: The pretraining of this method is directly a PINNs optimization with L-BFGS, which has been proposed in [2]. 

1] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning
framework for solving forward and inverse problems involving nonlinear partial differential
equations. Journal of Computational Physics, 378:686–707, 2019. ISSN 0021-9991. doi: https://
doi.org/10.1016/j.jcp.2018.10.045

[2] Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in
training pinns: A loss landscape perspective. In Forty-first International Conference on Ma-
chine Learning, ICML 2024

[3] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Char-
acterizing possible failure modes in physics-informed neural networks. Advances in Neural
Information Processing Systems, 34:26548–26560, 2021.