Solving Elliptic Optimal Control Problems via Neural Networks and
  Optimality System
In this work, we investigate a neural network based solver for optimal
control problems (without / with box constraint) for linear and semilinear
second-order elliptic problems. It utilizes a coupled system derived from the
first-order optimality system of the optimal control problem, and employs deep
neural networks to represent the solutions to the reduced system. We present an
error analysis of the scheme, and provide $L^2(\Omega)$ error bounds on the
state, control and adjoint in terms of neural network parameters (e.g., depth,
width, and parameter bounds) and the numbers of sampling points. The main tools
in the analysis include offset Rademacher complexity and boundedness and
Lipschitz continuity of neural network functions. We present several numerical
examples to illustrate the method and compare it with two existing ones.
Yongcheng Dai, Bangti Jin, Ramesh Sau, Zhi Zhou
# Solving Elliptic Optimal Control Problems using Physics Informed Neural Networks+

###### Abstract

In this work, we present and analyze a numerical solver for optimal control problems (without / with box constraint) for linear and semilinear second-order elliptic problems. The approach is based on a coupled system derived from the first-order optimality system of the optimal control problem, and applies physics informed neural networks (PINNs) to solve the coupled system. We present an error analysis of the numerical scheme, and provide \(L^{2}(\Omega)\) error bounds on the state, control and adjoint state in terms of deep neural network parameters (e.g., depth, width, and parameter bounds) and the number of sampling points in the domain and on the boundary. The main tools in the analysis include offset Rademacher complexity and boundedness and Lipschitz continuity of neural network functions. We present several numerical examples to illustrate the approach and compare it with three existing approaches.

**Key words**: optimal control, physics informed neural network, error estimate, elliptic problem

## 1 Introduction

Optimal control problems subject to partial differential equations (PDE) constraints represent a very important class of problems in practice and have found various important applications in science and engineering, e.g., fluid flow, heat conduction, structural optimization, microelectronics, crystal growth, vascular surgery, and cardiac medicine, to name just a few. The mathematical theory of distributed / boundary optimal control problems subject to linear / nonlinear elliptic / parabolic PDEs is well understood (see, e.g., the monographs [27, 39, 30]). Numerically, one of the most common approaches to solve such problems utilizes the first-order optimality system and introduces either adjoint state or Lagrangian multiplier, and then iteratively updates the control using established optimization algorithms, e.g., (accelerated / stochastic) gradient descent method, semismooth Newton method, and alternating direction method of multipliers. In practical implementation, the resulting continuous formulation is then often discretized by finite element method (FEM) and spectral methods.

In the past few years, deep neural networks (DNNs) have demonstrated remarkable performance across a wide range of applications, e.g., computer vision, imaging, and natural language processing. These successes have permeated into several scientific disciplines, and DNNs have also become very popular tools to approximate solutions to various PDEs. The ease of implementation, and the potential to overcome the notorious curse of dimensionality have sparked intensive interest in revisiting classical computational problems from a deep-learning perspective. In the last few years, a large number of neural solvers based on DNNs have been proposed for solving PDEs, e.g., physics informed neural network (PINN) [34], deep Ritz method [14], deep Galerkin method [36], weak adversarial networks [41], and deep least-squares method [9]. These
methods employ DNNs as ansatz functions to approximate the PDE solution directly. Among these neural solvers, PINN is undoubtedly the most prominent one, which is based on PDE residual minimization, penalty method for enforcing the boundary condition and computing the derivatives with automatic differentiation. The approach has been successfully applied to a wide range of scientific and engineering problems (see the reviews [23, 8]). See Section 2.2 for a more detailed description of PINN.

These impressive developments on neural PDE solvers have also sparked interest in solving optimal control problems governed by PDE constraint with deep learning techniques, and several different methods have been proposed. One first class of methods is based on the classical penalty methods, including standard penalty method (soft constraint) [32], and augmented Lagrangian method [28]. The idea of the penalty based methods is to transform the PDE-constrained optimization problem into an unconstrained one by introducing an augmented loss function that incorporates the PDE residual with a suitable penalty parameter. The idea is simple yet quite effective. The main drawback of penalty-based methods is that the optimization problem becomes increasingly ill-conditioned (and hence more challenging to optimize) as the penalty weights become large [5, 19, 28], but imposing the PDE constraint accurately does require a large penalty. Moreover, in the presence of multiple constraints, one has to tune multiple penalty weights, which is crucial to yield reasonable results yet the tuning can be highly nontrivial, time-consuming and in fact very tedious. The augmented Lagrangian method partly alleviates the necessity of using a large penalty weight [19]. The second class of methods employs the adjoint technique [10], which solves the adjoint state equation and computes the derivative of the reduced objective with respect to the control variable. One method in the class is the direct adjoint looping (DAL) [27, 20], which iteratively updates the control using a gradient descent scheme computed via the adjoint system. Using the idea of DAL, very recently, Yin et al [40] proposed the adjoint-oriented neural network (AONN). It can be viewed as a neural network realization of the adjoint-based gradient descent scheme. It involves solving the state and adjoint and then updates the control sequentially, and thus consists of inner and outer loops. The authors reported excellent empirical performance on a number of challenging settings, e.g., (parametric) Navier-Stokes problem. See Section 3.4 for further descriptions of these methods for elliptic optimal control problems. Demo et al [12] proposed an extended PINN, which augments DNN inputs with parameters so that the Karush-Kuhn-Tucker (KKT) system and neural networks are combined for parametric optimal control problems; see also [2]. Recently, Song et al [38] develop a PINN and ADMM based approach for efficiently solving nonsmooth PDE constrained optimization.

In this work, we investigate a class of numerical solvers based on DNNs for solving elliptic optimal control problems using PINN, which is termed coupled Physics Informed Neural Network (C-PINN in short). The solver is derived as follows. First, we derive a coupled system (i.e., KKT system) for the continuous optimal control problem, which consists of a state equation, an adjoint state equation, and an equality / variational inequality for the unconstrained and constrained cases, respectively. The variational equality allows recovering the optimal control from the optimal adjoint, and the variational inequality allows realizing the optimal control as a min-max projection of the optimal adjoint onto the constrained set. Therefore, the KKT system is equivalent to a (reduced) coupled system for the state and adjoint, which consists of a state equation with a min-max projection operator in the source, and an adjoint state equation. Second, we construct a (empirical) PINN loss for the (reduced) coupled system, which consists of the PDE residuals and boundary penalties. The method is illustrated for the unconstrained case, constrained case, and semilinear case. Third, we give an error analysis of the approach for linear elliptic optimal control problems, and derive an \(L^{2}(\Omega)\) error estimate for the control, state, and adjoint state approximations. The error bounds can be decomposed into two parts: approximation and statistical errors. The approximation error arises from approximating the solution of the control problem using neural networks and is bounded using approximation theory of DNNs with general activation functions [15], and the statistical error arises from using Monte Carlo methods to approximate the (potentially) high-dimensional integrals and is bounded using offset Rademacher complexity [26, 13] and certain technical estimates on the DNN functions. Fourth and last, we present a set of numerical experiments in two- and multi-dimensions, and give a comparative study with existing neural solvers for optimal control problems. The derivation of the new computational technique based on the coupled system and PINN, its error analysis and the numerical evaluation represent the main contributions of the work. To the best of our knowledge, the coupled PINN represents the first DNN based solver with convergence guarantees for elliptic optimal control problems.

The rest of the paper is organized as follows. In Section 2, we describe some preliminaries on neural networks and PINN. In section 3, we describe elliptic optimal control problems via PINN based on the
coupled system reformulation of the optimal control problems. In Section 4, we present an error analysis for both constrained and unconstrained cases, which includes a detailed analysis of the approximation error and statistical errors, which in turn relies on suitable technical estimates on DNNs in Section 5. Finally, in Section 6, we present numerical examples to validate the theoretical results. Throughout the notation \(c\) denotes a generic constant which may change at each occurrence, but it is always independent of the neural network parameters.

## 2 Preliminaries

In this section we review preliminaries on deep neural networks and physics informed neural networks.

### Deep neural networks

Due to their extraordinary expressivity for approximating functions, deep neural networks (DNNs) have been widely used, and play a vital role in the great success of deep learning. In this work, we consider fully connected feedforward neural networks. Fix \(L\in\mathbb{N}\), and integers \(\{n_{\ell}\}_{\ell=0}^{L}\subset\mathbb{N}\), with \(n_{0}=d\) and \(n_{L}=1\). A neural network \(f_{\theta}:\mathbb{R}^{d}\to\mathbb{R}\), with \(\theta\in\mathbb{R}^{N_{\theta}}\) being the DNN parameters, is defined recursively by

\[f^{(0)}(x) =x,\] \[f^{(\ell)}(x) =\rho(A^{(\ell)}f^{(\ell-1)}+b^{(\ell)}),\quad\ell=1,\dots,L-1,\] \[f_{\theta}:=f^{(L)}(x) =A^{(L)}f^{(\ell-1)}+b^{(L)},\]

where \(A^{(\ell)}=[a_{ij}^{(\ell)}]\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}\), and \(b^{(\ell)}=(b_{i}^{(\ell)})\in\mathbb{R}^{n_{\ell}}\) are known as the weight matrix and bias vector at the \(\ell\)th layer, respectively. \(L\) is called the depth of the neural network, and \(W:=\max(n_{\ell},\ell=1,\dots,L)\) is called the width of the neural network. We denote by \(\boldsymbol{n}_{i}\), \(i=1,\dots,L\), as the number of nonzero weights on the first \(i\) layers. Then \(\boldsymbol{n}_{L}\) is the total number of nonzero weights. We use the notation \(\mathcal{N}_{\rho}(L,\boldsymbol{n}_{L},R)\) to refer to the collection of DNN functions generated by \(\rho\) of depth \(L\), total number \(\boldsymbol{n}_{L}\) of nonzero weights and each weight bounded by \(R\). Throughout, we focus on two activation functions \(\rho\): hyperbolic tangent \(\rho(t)=\frac{e^{t}-e^{-t}}{e^{t}+e^{-T}}\) and sigmoid \(\rho(t)=\frac{1}{1+e^{-t}}\), which are popular in the literature on PINN.

**Lemma 2.1**.: _The following estimates hold._

1. _If_ \(\rho(t)=\tanh(t)\)_, then_ \(\|\rho^{(i)}\|_{L^{\infty}(\mathbb{R})}\leq 1\)_,_ \(i=0,1,2\)_, and_ \(\|\rho^{\prime\prime\prime}\|_{L^{\infty}(\mathbb{R})}\leq 2\)_._
2. _If_ \(\rho(t)=\frac{1}{1+e^{-t}}\)_, then_ \(\|\rho^{(i)}\|_{L^{\infty}(\Omega)}\leq 1\)_,_ \(i=0,1,2,3\)_._

Proof.: The lemma follows from direct computation. Indeed, for \(\rho(t)=\tanh(t)\), we have \(\|\rho\|_{L^{\infty}(\mathbb{R})}\leq 1\) and \(\rho^{\prime}(t)=1-\rho(t)^{2}\), \(\rho^{\prime\prime}(t)=-2\rho(t)(1-\rho(t)^{2})\), \(\rho^{\prime\prime\prime}(t)=(6\rho(t)^{2}-2)(1-\rho(t)^{2})\). Similarly, for \(\rho(t)=\frac{1}{1+e^{-t}}\), clearly \(\rho(t)\in(0,1)\), and \(\rho^{\prime}(t)=\frac{e^{-t}}{(1+e^{-t})^{2}}=\rho(t)(1-\rho(t))\in(0,\frac{ 1}{2})\), \(\rho^{\prime\prime}(t)=\rho(t)(1-\rho(t))(1-2\rho(t))\), \(\rho^{\prime\prime\prime}(t)=\rho(t)(1-\rho(t))[1-6\rho(t)+6\rho(t)^{2}]\). Then the desired assertions follow immediately. 

### Physics informed neural networks

Now we describe physics informed neural networks (PINNs) of Rassi et al [34], for solving the elliptic boundary value problem

\[\begin{cases}-\Delta y=f,&\text{in }\Omega,\\ \quad y=0,&\text{on }\partial\Omega,\end{cases} \tag{2.1}\]

where \(\Omega\subset\mathbb{R}^{d}\) is a smooth domain with a boundary \(\partial\Omega\), and \(f\in L^{2}(\Omega)\) is the source. PINN is based on the principle of PDE residual minimization. For problem (2.1), the continuous loss \(\mathcal{L}(u)\) is given by

\[\mathcal{L}(y)=\frac{1}{2}\|\Delta y+f\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2} \|y\|_{L^{2}(\partial\Omega)}^{2}, \tag{2.2}\]
where the penalty term \(\|y\|_{L^{2}(\partial\Omega)}^{2}\) (with weight \(\alpha>0\)) is to approximately enforce the zero Dirichlet boundary condition. We seek an approximation \(y_{\theta}\in\mathcal{Y}=\mathcal{N}_{\rho}(L,\mathbf{n}_{L},R)\), and then discretize relevant integrals using Monte Carlo method. Let \(U(D)\) be the uniform distribution over a set \(D\), and \(|D|\) the Lebesgue measure of \(D\), and \(\mathbb{E}_{\nu}\) taking expectation with respect to a probability distribution \(\nu\). Then the loss \(\mathcal{L}(y_{\theta})\) is given by

\[\mathcal{L}(y_{\theta})=\frac{1}{2}|\Omega|\mathbb{E}_{U(\Omega)}[(\Delta y_{ \theta}(X)+f(X))^{2}]+\frac{\alpha}{2}|\partial\Omega|\mathbb{E}_{U(\partial \Omega)}[(y_{\theta}(Y))^{2}],\]

and hence also known as population loss. Let the sampling points \(\{X_{i}\}_{i=1}^{n_{d}}\) and \(\{Y_{j}\}_{j=1}^{n_{b}}\) be identically and independently distributed (i.i.d.), uniformly on the domain \(\Omega\) and boundary \(\partial\Omega\), respectively, i.e., \(\mathbb{X}=\{X_{i}\}_{i=1}^{n_{d}}\sim U(\Omega)\) and \(\mathbb{Y}=\{Y_{j}\}_{j=1}^{n_{b}}\sim U(\partial\Omega)\). Then the empirical loss \(\widehat{\mathcal{L}}(y_{\theta})\) is given by

\[\widehat{\mathcal{L}}(y_{\theta})=\frac{|\Omega|}{2n_{d}}\sum_{i=1}^{n_{d}}( \Delta y_{\theta}(X_{i})+f(X_{i}))^{2}+\frac{\alpha|\partial\Omega|}{2n_{b}} \sum_{j=1}^{n_{b}}(y_{\theta}(Y_{j}))^{2}.\]

Note that the resulting optimization problem is well-posed due to the box constraint on the DNN parameters \(\theta\), i.e., \(|\theta|_{\ell^{\infty}}\leq R\) for suitable \(R\), which induces a compact set in \(\mathbb{R}^{N_{\theta}}\). Meanwhile the empirical loss \(\widehat{\mathcal{L}}(y_{\theta})\) is continuous in \(\theta\), when the activation function \(\rho\) is smooth. In the absence of the box constraint, the optimization problem might not have a finite minimizer, due to the lack of weak closedness of the DNN function class [33]. The theoretical analysis of the PINNs has been carried out in [21, 35, 31, 29].

## 3 Elliptic optimal control via PINN

In this section, we develop the PINN for solving elliptic optimal control problems via the coupled optimality system, and term the resulting method as C-PINN. We also survey three existing neural solvers in the context of elliptic optimal control problems.

### Optimal control without control constraint

Let \(\Omega\subset(-1,1)^{d}\subset\mathbb{R}^{d}\) (\(d\geq 1\)) be a smooth domain with a boundary \(\partial\Omega\). Consider the following distributed optimal control problem

\[\min\;J(y,u):=\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{\lambda}{2}\|u \|_{L^{2}(\Omega)}^{2}, \tag{3.1}\]

subject to the following elliptic PDE constraint

\[\begin{cases}-\Delta y=f+u,&\text{in }\Omega,\\ y=0,&\text{on }\partial\Omega,\end{cases} \tag{3.2}\]

where \(u\in L^{2}(\Omega)\) is the control, \(f\in L^{2}(\Omega)\) and \(y_{d}\in L^{2}(\Omega)\) are the given source and target function, respectively, and the scalar \(\lambda>0\) balances the two terms in the objective \(J(y,u)\). It is known that the optimal control problem (3.1)-(3.2) has a unique solution \((\bar{y},\bar{u})\)[39, Theorem 2.17]. The main goal is to approximate the solution \((\bar{y},\bar{u})\) by a DNN pair \((\bar{y}_{\theta},\bar{u}_{\kappa})\) via PINN.

It is known that problem (3.1) is equivalent to the following Karush-Kuhn-Tucker (KTT) system for the optimal tuple \((\bar{y},\bar{p},\bar{u})\) (see, e.g., [27, Chapter 2, Theorem 1.4]):

\[\begin{cases}-\Delta y=f+u,&\text{in }\Omega,\\ y=0,&\text{on }\partial\Omega,\\ -\Delta p=y-y_{d},&\text{in }\Omega,\\ p=0,&\text{on }\partial\Omega,\\ u=-\lambda^{-1}p,&\text{in }\Omega.\end{cases} \tag{3.3}\]

The last line is the first-order necessary optimality condition of problem (3.1)-(3.2).
The starting point of the coupled approach is system (3.3). Upon eliminating the first-order necessary condition \(u=\lambda^{-1}p\), problem (3.1)-(3.2) amounts to solving the following reduced system in \(y\) and \(p\):

\[\begin{cases}-\Delta y=f-\lambda^{-1}p,&\text{in }\Omega,\\ y=0&\text{on }\partial\Omega,\\ -\Delta p=y-y_{d}&\text{in }\Omega,\\ p=0&\text{on }\partial\Omega.\end{cases} \tag{3.4}\]

In view of the (reduced) coupled system (3.4), we define the following continuous loss

\[\mathcal{L}(y,p)=\left\|\Delta y+f-\lambda^{-1}p\right\|_{L^{2}(\Omega)}^{2}+ \alpha_{i}\left\|\Delta p+y-y_{d}\right\|_{L^{2}(\Omega)}^{2}+\alpha_{b}^{y} \left\|y\right\|_{L^{2}(\partial\Omega)}^{2}+\alpha_{b}^{p}\left\|p\right\|_ {L^{2}(\partial\Omega)}^{2}, \tag{3.5}\]

where the positive scalars \(\alpha_{i},\alpha_{b}^{y}\) and \(\alpha_{b}^{p}\) are penalty parameters to enforce the PDE residual (of \(p\)) and zero Dirichlet boundary conditions of \(y\) and \(p\). Below we denote by \(\mathbf{\alpha}=(\alpha_{i},\alpha_{b}^{y},\alpha_{b}^{p})\). Let \(U(\Omega)\) and \(U(\partial\Omega)\) be the uniform distributions on \(\Omega\) and \(\partial\Omega\), respectively. Then the loss \(\mathcal{L}(y,p)\) can be rewritten as

\[\mathcal{L}(y,p)= |\Omega|\mathbb{E}_{X\sim U(\Omega)}\left[\left(\Delta y(X)+f(X) -\lambda^{-1}p(X)\right)^{2}\right]+\alpha_{i}|\Omega|\mathbb{E}_{X\sim U( \Omega)}\left[\left(\Delta p(X)+y(X)-y_{d}(X)\right)^{2}\right]\] \[+\alpha_{b}^{y}|\partial\Omega|\mathbb{E}_{Y\sim U(\partial \Omega)}\left[\left(y(Y)\right)^{2}\right]+\alpha_{b}^{p}|\partial\Omega| \mathbb{E}_{Y\sim U(\partial\Omega)}\left[\left(p(Y)\right)^{2}\right]. \tag{3.6}\]

Following the standard PINN paradigm in Section 2.2, we employ the DNN spaces \(\mathcal{Y}\) and \(\mathcal{P}\) to approximate the state \(y\) and adjoint state \(p\), respectively, and then approximate the expectations via Monte Carlo methods. Let \((y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}\) be the DNN approximation, and \(\mathbb{X}=\{X_{i}\}_{i=1}^{n_{d}}\) and \(\mathbb{Y}=\{Y_{j}\}_{j=1}^{n_{b}}\) be i.i.d. samples drawn from \(U(\Omega)\) and \(U(\partial\Omega)\), respectively. Then the empirical PINN loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) is given by

\[\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})= \frac{|\Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\left(\Delta y_{\theta}(X _{i})+f(X_{i})-\lambda^{-1}p_{\sigma}(X_{i})\right)^{2}+\alpha_{i}\frac{| \Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\left(\Delta p_{\sigma}(X_{i})+y_{\theta}(X _{i})-y_{d}(X_{i})\right)^{2}\] \[+\alpha_{b}^{y}\frac{|\partial\Omega|}{n_{b}}\sum_{j=1}^{n_{b}} \left(y_{\theta}(Y_{j})\right)^{2}+\alpha_{b}^{p}\frac{|\partial\Omega|}{n_{b }}\sum_{j=1}^{n_{b}}\left(p_{\sigma}(Y_{j})\right)^{2}. \tag{3.7}\]

Let \((y_{\theta^{*}},p_{\sigma^{*}})\in\mathcal{Y}\times\mathcal{P}\) be a minimizer of \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) over \(\mathcal{Y}\times\mathcal{P}\), or equivalently with respect to the DNN parameter vectors \(\theta\) and \(\sigma\). In practice, this is commonly achieved by off-shelf optimizers, e.g., limited memory BFGS [7] and Adam [24], which are implemented in many public software frameworks, e.g., PyTorch or TensorFlow. To evaluate the loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\), one requires computing the derivative of the network outputs \(y_{\theta}\) and \(p_{\sigma}\) with respect to the input \(x\) (i.e., spatial derivatives), and to apply a first-order optimizer, one requires computing the derivative of the loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) with respect to the DNN parameters \(\theta\) and \(\sigma\). Both can be realized using automatic differentiation efficiently [3], e.g., torch.autograd in PyTorch. Given a (local) minimizer \((\theta^{*},\sigma^{*})\), the corresponding DNN approximations of the state \(\bar{y}\in\mathcal{Y}\) and adjoint \(\bar{p}\in\mathcal{P}\) are given by \(y_{\theta^{*}}\) and \(p_{\sigma^{*}}\), respectively. In view of the first-order necessary condition, the DNN approximation \(u^{*}\) of the optimal control \(\bar{u}\) is given by \(-\lambda^{-1}p_{\sigma^{*}}\).

**Remark 3.1**.: _Demo et al [12] proposed an extended PINN, based on the first-order KKT system for parametric optimal control problems, where the DNNs also take the parameters as inputs besides the spatial coordinates. The approach applies also to the deterministic case. In the standard elliptic case, it is based on the full optimality system (3.3), and the corresponding population loss is given by_

\[\mathcal{L}(y,p,u)= \left\|\Delta y+f+u\right\|_{L^{2}(\Omega)}^{2}+\alpha_{i}\left\| \Delta p+y-y_{d}\right\|_{L^{2}(\Omega)}^{2}\] \[+\alpha_{b}^{y}\left\|y\right\|_{L^{2}(\partial\Omega)}^{2}+ \alpha_{b}^{p}\left\|p\right\|_{L^{2}(\partial\Omega)}^{2}+\alpha_{i}^{u}\|u+ \lambda^{-1}p\|_{L^{2}(\Omega)}^{2}.\]

_Compared with C-PINN in (3.5), the loss \(\mathcal{L}(y,p,u)\) requires also one neural network \(u_{\kappa}\in\mathcal{U}\) for approximating the control \(u\), and tuning the weight \(\alpha_{i}^{u}\), which might be nontrivial. Also, the work [12] does not discuss the constrained case, where the control \(u\) is inherently nonsmooth, and thus numerically challenging to approximate using PINN [40]. A similar approach to [12] was suggested in [2]._
### Optimal control with box constraint

Now we turn to the constrained optimal control problem, i.e., the control \(u\) is from a box constrain set \(U=\{u\in L^{2}(\Omega):u_{a}\leq u\leq u_{b}\text{ a.e. in }\Omega\}\), with \(u_{a}<u_{b}\). It is known that the optimal tuple \((\bar{y},\bar{p},\bar{u})\) satisfies the following first-order necessary optimality system [39, Theorem 2.25, equation (2.52), p. 67]

\[\left\{\begin{aligned} -\Delta y&=f+u,& \text{in }\Omega,\\ y&=0,&\text{on }\partial\Omega,\\ -\Delta p&=y-y_{d},&\text{in }\Omega,\\ p&=0,&\text{on }\partial\Omega,\\ (u+\lambda^{-1}p,v-u)&\geq 0,&\forall v \in U.\end{aligned}\right. \tag{3.8}\]

Now we introduce the pointwise projection operator \(P_{U}\) into the constraint set \(U\), defined by \(P_{U}(v)(x)=\min(\max(v(x),u_{a}),u_{b})\). Then the variational inequality in (3.8) can be rewritten as

\[u=P_{U}(-\lambda^{-1}p),\]

and accordingly, the (reduced) optimality system is given by

\[\left\{\begin{aligned} -\Delta y&=f+P_{U}(-\lambda^{-1}p),& \text{in }\Omega,\\ y&=0,&\text{on }\partial\Omega,\\ -\Delta p&=y-y_{d},&\text{in }\Omega,\\ p&=0,&\text{on }\partial\Omega.\end{aligned}\right. \tag{3.9}\]

The reduced optimality system (3.9) naturally motivates the following continuous PINN loss:

\[\mathcal{L}(y_{\theta},p_{\sigma})= |\Omega|\mathbb{E}_{X\sim U(\Omega)}\left[\left(\Delta y_{\theta} (X)+f(X)+P_{U}(-\lambda^{-1}p_{\sigma}(X))\right)^{2}\right]+\alpha_{b}^{y}| \partial\Omega|\mathbb{E}_{Y\sim U(\partial\Omega)}\left[\left(y_{\theta}(Y) \right)^{2}\right]\] \[+\alpha_{i}|\Omega|\mathbb{E}_{X\sim U(\Omega)}\left[\left(\Delta p _{\sigma}(X)+y_{\theta}(X)-y_{d}(X)\right)^{2}\right]+\alpha_{b}^{p}|\partial \Omega|\mathbb{E}_{Y\sim U(\partial\Omega)}\left[\left(p_{\sigma}(Y)\right)^ {2}\right].\]

Upon approximating the expectations with Monte Carlo methods using i.i.d. samples \(\mathbb{X}=\{X_{i}\}_{i=1}^{n_{d}}\sim U(\Omega)\) and \(\mathbb{Y}=\{Y_{j}\}_{j=1}^{n_{b}}\sim U(\partial\Omega)\), we obtain the following empirical PINN loss

\[\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})= \frac{|\Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\Bigl{(}\Delta y_{\theta} (X_{i})+f(X_{i})+P_{U}(-\lambda^{-1})p_{\sigma}(X_{i})\Bigr{)}^{2}+\alpha_{b}^ {y}\frac{|\partial\Omega|}{n_{b}}\sum_{j=1}^{n_{b}}\left(y_{\theta}(Y_{j}) \right)^{2}\] \[+\alpha_{i}\frac{|\Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\left(\Delta p _{\sigma}(X_{i})+y_{\theta}(X_{i})-y_{d}(X_{i})\right)^{2}+\alpha_{b}^{p}\frac {|\partial\Omega|}{n_{b}}\sum_{j=1}^{n_{b}}\left(p_{\sigma}(Y_{j})\right)^{2}.\]

The empirical loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) is then minimized with respect to the DNN parameters \((\theta,\sigma)\), with a minimizer \((\theta^{*},\sigma^{*})\). The approximation to the optimal control \(\bar{u}\) is given by \(P_{U}(-\lambda^{-1}p_{\sigma^{*}})\). Note that the change from the unconstrained case is rather minimal in terms of formulation and implementation.

### Optimal control for semilinear problems

C-PINN applies also to more complex optimal control problems. We illustrate this with semilinear elliptic problems, one popular model class in PDE optimal control. Consider the following optimal control problem

\[\min\;J(y,u):=\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{\lambda}{2}\|u \|_{L^{2}(\Omega)}^{2}, \tag{3.10}\]

subject to the following semilinear elliptic PDE constraint

\[\left\{\begin{aligned} -\Delta y&=f(y)+u,& \text{in }\Omega,\\ y&=0,&\text{on }\partial\Omega,\end{aligned}\right. \tag{3.11}\]
where \(u\in L^{2}(\Omega)\) is the control. The semilinear term \(f:\mathbb{R}\to\mathbb{R}\) should satisfy suitable conditions (e.g., monotone) so that the state problem (3.11) is well-posed [39, Theorem 4.4]. Then the optimal control problem has a locally unique solution \((\bar{y},\bar{u})\)[39, Theorem 4.15, p. 208]. The KKT system for the optimal control problem (3.10)-(3.11) is given by [39, Theorem 4.20, p. 216]:

\[\begin{cases}-\Delta y=f(y)+u,&\text{in }\Omega,\\ y=0,&\text{on }\partial\Omega,\\ -\Delta p=f^{\prime}(y)p+y-y_{d},&\text{in }\Omega,\\ p=0,&\text{on }\partial\Omega,\\ u=-\lambda^{-1}p,&\text{in }\Omega.\end{cases}\]

Following the preceding derivations for linear elliptic PDEs, the empirical PINN loss \(\widehat{L}(y_{\theta},p_{\sigma})\) for problem (3.10)-(3.11) is given by

\[\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})= \frac{|\Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\Bigl{(}\Delta y_{\theta} (X_{i})+f(y_{\theta}(X_{i}))-\lambda^{-1}p_{\sigma}(X_{i})\Bigr{)}^{2}+ \alpha_{b}^{y}\frac{|\partial\Omega|}{n_{b}}\sum_{j=1}^{n_{b}}\big{(}y_{ \theta}(Y_{j})\big{)}^{2}\] \[+\alpha_{i}\frac{|\Omega|}{n_{d}}\sum_{i=1}^{n_{d}}\big{(}\Delta p _{\sigma}(X_{i})+f^{\prime}(y_{\theta}(X_{i}))p_{\sigma}(X_{i})+y_{\theta}(X_{ i})-y_{d}(X_{i}))^{2}+\alpha_{b}^{p}\frac{|\partial\Omega|}{n_{b}}\sum_{j=1}^{n_{b}} \big{(}p_{\sigma}(Y_{j})\big{)}^{2}.\]

The adaptation of the approach to semilinear problems with box constraint on the control \(U\) is direct.

### Existing DNN based approaches to optimal control

The use of DNNs to solve PDE optimal control problems has just started to receive attention recently. Three different methods have been developed, i.e., augmented Lagrangian method (ALM) [28], penalty method (PM) [32], and adjoint oriented neural network (AONN) [40], which we describe below for unconstrained linear elliptic control problems. We mention also necessary changes in the constrained case.

PM transforms the constrained problem (3.1)-(3.2) into an unconstrained one by incorporating the PINN loss \(\mathcal{L}_{\text{pinn}}\) as

\[\mathcal{L}_{\text{pmn}}(y_{\theta},u_{\kappa}):= J(y_{\theta},u_{\kappa})+\mu\mathcal{L}_{\text{pinn}}(y_{ \theta},u_{\kappa})\] \[= \Bigl{(}\frac{1}{2}\|y_{\theta}-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac {\lambda}{2}\|u_{\kappa}\|_{L^{2}(\Omega)}^{2}\Bigr{)}+\mu\Bigl{(}\frac{1}{2} \|F(y_{\theta},u_{\kappa})\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2}\|y_{\theta} \|_{L^{2}(\partial\Omega)}^{2}\Bigr{)},\]

where \(\mu>0\) is the penalty parameter for the PINN loss \(\mathcal{L}_{\text{pinn}}\), and \(F(y_{\theta},u_{\kappa})=\Delta y_{\theta}+f+u_{\kappa}\) is the PDE residual for the DNN approximation \(y_{\theta}\in\mathcal{Y}\) and \(u_{\kappa}\in\mathcal{U}\). Like in PINN, the integrals are approximated using Monte Carlo methods, leading to an empirical loss \(\widehat{\mathcal{L}}_{\text{pmn}}(p_{\theta},u_{\kappa})\) to be minimized in the DNN parameters \(\theta\) and \(\kappa\). The penalty factor \(\mu\) can be fixed _a priori_. Then PM involves solving only one optimization problem but enforces only the PDE constraint approximately. To strictly enforce the constraint, one should send \(\mu\to\infty\) and solve a sequence of problems, commonly via a path-following strategy [28, 18]. Specifically, fix \(\mu_{0}>0\), and set \(\mu_{k+1}=\beta\mu_{k}\), with the increasing factor \(\beta>1\), and then one solves a sequence of problems for \(\mu_{k}\), with the solution \((\theta_{k-1}^{*},\kappa_{k-1}^{*})\) of the optimization problem with weight \(\mu_{k-1}\) as the initial guess of the current one with \(\mu_{k}\). One drawback of PM is that it becomes increasingly more ill-conditioned as \(\mu\to\infty\) (and thus more challenging to resolve) [5, 28]. Further, the presence of additional constraints requires more penalty terms, for which the tuning of the penalty weights can be very tedious [25, 40]. Mowlavi and Nabi [32] presented guidelines for tuning parameters in PM for PDE constrained optimal control.

ALM transforms a PDE constrained optimal control problem into an unconstrained one using Lagrangian multipliers. The use of Lagrangian multiplier can mitigate the ill-conditioning of PM [19]. Let \((\cdot,\cdot)\) and \((\cdot,\cdot)_{L^{2}(\partial\Omega)}\) denote the \(L^{2}(\Omega)\) and \(L^{2}(\partial\Omega)\) inner products. Then the loss \(\mathcal{L}_{\text{alm}}\) reads

\[\mathcal{L}_{\text{alm}}(y_{\theta},u_{\kappa},\eta^{d},\eta^{b}):=J(y_{\theta}, u_{\kappa})+\mu\mathcal{L}_{\text{pinn}}(y_{\theta},u_{\kappa})+(\eta^{d},F(y_{ \theta},u_{\kappa}))+(\eta^{b},y_{\theta})_{L^{2}(\partial\Omega)},\]

where the functions \(\eta^{d}\in L^{2}(\Omega),\eta^{b}\in L^{2}(\partial\Omega)\) are Lagrangian multipliers for the governing equation \(F(y,u)=0\) in \(\Omega\) and boundary condition \(y=0\) on \(\partial\Omega\), respectively, which may also be realized by DNNs. This loss is
often minimized alternatingly: for fixed \((\eta_{k}^{d},\eta_{k}^{b})\), we minimize \(\mathcal{L}_{\mathrm{alm}}(y_{\theta},u_{\kappa},\eta^{d},\eta^{b})\) in \((y_{\theta},u_{\kappa})\):

\[(y_{\theta^{k}},u_{\kappa^{k}})=\arg\min_{(y_{\theta},u_{\kappa})\in\mathcal{Y }\times\mathcal{P}}\mathcal{L}_{\mathrm{alm}}(y_{\theta},u_{\kappa},\eta_{k}^{ d},\eta_{k}^{b}),\]

and then update the Lagrangian multipliers \((\eta_{k+1}^{d},\eta_{k+1}^{b})\) by an Uzawa type algorithm:

\[\eta_{k+1}^{d}=\eta_{k}^{d}+\mu F(y_{\theta^{k}},u_{\kappa^{k}})\quad\text{and }\quad\eta_{k+1}^{b}=\eta_{k}^{b}+\mu\alpha y_{\theta^{k}}.\]

In practice, the integrals are approximated by evaluations on fixed collocation points. Thus it suffices to discretize the multipliers pointwise point-wise [28].

AONN is also based on the optimality system (3.3), but utilizes the adjoint \(p\) to compute the gradient of the reduced cost \(J(u)=J(y(u),u)\) in \(u\) only, so as to combine the advantages of DAL and deep learning techniques. The total derivative \(\mathrm{d}_{u}J(y,u)=\alpha u+p\) is used to update \(u\) as

\[u_{\kappa^{k+1}}=\arg\min_{u_{\kappa}\in\mathcal{U}}\|u_{\kappa}-(u_{\kappa^{ k}}-s^{k}\mathrm{d}_{u}J(y_{\theta^{k}},u_{\kappa^{k}}))\|_{L^{2}(\Omega)}^{2}, \tag{3.12}\]

where the step size \(s^{k}>0\) should be suitably chosen, and the minimization is to project the update into the set \(\mathcal{U}\). In AONN, the state \(y_{\theta}\) and adjoint state \(p_{\sigma}\) are constructed sequentially, and a third DNN \(u_{\kappa}\in\mathcal{U}\) is employed to approximate \(u\). Similar to PM, the initial guess for each subproblem is given by the last minimizer. Unlike PM and ALM, AONN gives an all-at-once approximation of the control, state, and adjoint. Numerically, AONN outperforms penalty based methods over a range of problems [40].

The presence of box constraint on the control \(u\in U\) requires a slightly different treatment. In PM, one can add an additional penalty term to enforce the constraint:

\[\mathcal{L}_{\mathrm{pm-c}}(y_{\theta},u_{\kappa})=J(y_{\theta},u_{\kappa})+ \mu\mathcal{L}_{\mathrm{pinn}}(y_{\theta},u_{\kappa})+\frac{\mu^{\prime}}{2} \|u_{\kappa}-P_{U}(u_{\kappa})\|_{L^{2}(\Omega)}^{2}, \tag{3.13}\]

where \(\mu^{\prime}\) is a penalty weight, and \(u_{\kappa}\in\mathcal{U}\) is a DNN approximation of the control variable \(u\). The drawback is that the obtained approximation \(u_{\kappa}\) may be infeasible. Alternatively, one can also enforce the constraint by applying the point-wise projection \(P_{U}\) on \(u_{\kappa}\) directly

\[\mathcal{L}_{\mathrm{pm}}(y_{\theta},u_{\kappa})=J(y_{\theta},P_{U}(u_{\kappa }))+\mu\mathcal{L}_{\mathrm{pinn}}(y_{\theta},P_{U}(u_{\kappa})).\]

The projection \(P_{U}\) leads to vanishing gradient when the constraint is active, and the optimizer may get stuck in a local minimum and hence requires a good initial guess. Indeed, box constraints also lead to non-smoothness of the loss, causing numerical issues [16, 28]. AONN treats the box constraint \(u\in U\) by

\[u_{\kappa^{k+1}}=\arg\min_{u_{\kappa}\in\mathcal{U}}\|u_{\kappa}-P_{U}(u_{ \kappa^{k}}-s^{k}d_{u}J(y_{\theta^{k}},u_{\kappa^{k}}))\|_{L^{2}(\Omega)}^{2},\]

where the initial guess for the DNN parameters \(\kappa^{k+1}\) is \(\kappa^{k}\). For the coupled PINN, the constraint is enforced by substituting \(P_{U}(p_{\sigma})\) into the PDE system. Numerically, after projection, the control becomes a piece-wise function, which is not fully parameterized by the DNN parameters. We choose the initial DNN parameters with a proper scale so that most part of the function satisfies the constraint. Alternatively one may also apply smoothed version of the projection operator \(P_{U}\) like sigmoid to find a proper initial guess.

## 4 Error analysis

Now we present analyze C-PINN for linear elliptic optimal control problems.

### Fundamental estimates

First, we provide two fundamental results bounding the error in terms of the population loss (i.e., weak coercivity results). We begin by discussing optimal control problems without constraints.
**Lemma 4.1**.: _Let \((\bar{y},\bar{p})\) be the solution tuple to the system (3.4) with \(\bar{u}=-\lambda^{-1}\bar{p}\). Then for any \((y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}\), with \(u_{\sigma}=-\lambda^{-1}p_{\sigma}\), the following estimate holds_

\[\|\bar{y}-y_{\theta}\|_{L^{2}(\Omega)}^{2}+\|\bar{p}-p_{\sigma}\|_{L^{2}( \Omega)}^{2}+\|\bar{u}-u_{\sigma}\|_{L^{2}(\Omega)}^{2}\leq c(\mathbf{\alpha}, \lambda)\mathcal{L}(y_{\theta},p_{\sigma}).\]

Proof.: To deal with the nonzero Dirichlet boundary conditions with \(y_{\theta}\) and \(p_{\sigma}\), we define

\[\begin{cases}-\Delta\zeta^{y}=0,&\text{in}\;\Omega,\\ \quad\quad\zeta^{y}=y_{\theta},&\text{on}\;\partial\Omega,\end{cases}\quad \text{and}\quad\begin{cases}-\Delta\zeta^{p}=0,&\text{in}\;\Omega,\\ \quad\quad\zeta^{p}=p_{\sigma},&\text{on}\;\partial\Omega.\end{cases} \tag{4.1}\]

That is, \(\zeta^{y}\) and \(\zeta^{p}\) are harmonic extensions of \(y_{\theta}\) and \(p_{\sigma}\), respectively. By the standard elliptic regularity theory [4, Theorem 4.2, p. 870], the following stability estimates hold

\[\|\zeta^{y}\|_{L^{2}(\Omega)}\leq c\|y_{\theta}\|_{L^{2}(\partial\Omega)} \quad\text{and}\quad\|\zeta^{p}\|_{L^{2}(\Omega)}\leq c\|p_{\sigma}\|_{L^{2}( \partial\Omega)}. \tag{4.2}\]

Let \(e_{y}=\bar{y}-y_{\theta}\) and \(e_{p}=\bar{p}-p_{\sigma}\), also \(\tilde{e}_{y}=e_{y}+\zeta^{y}\) and \(\tilde{e}_{p}=e_{p}+\zeta^{p}\). Then \(\tilde{e}_{y}\) and \(\tilde{e}_{p}\) satisfy

\[\begin{cases}-\Delta\tilde{e}_{y}+\lambda^{-1}\tilde{e}_{p}=f- \lambda^{-1}p_{\sigma}+\Delta y_{\theta}+\lambda^{-1}\zeta^{p},&\text{in}\; \Omega,\\ \quad\quad\tilde{e}_{y}=0,&\text{on}\;\partial\Omega,\end{cases} \tag{4.3a}\] \[\begin{cases}-\Delta\tilde{e}_{p}-\tilde{e}_{y}=y_{\theta}-y_{d} +\Delta p_{\sigma}-\zeta^{y},&\text{in}\;\Omega,\\ \quad\quad\tilde{e}_{p}=0,&\text{on}\;\partial\Omega.\end{cases} \tag{4.3b}\]

Multiplying (4.3a) by \(\tilde{e}_{y}\) and (4.3b) by \(\lambda^{-1}\tilde{e}_{p}\), integrating over \(\Omega\) and adding the resulting identities yield

\[\|\nabla\tilde{e}_{y}\|_{L^{2}(\Omega)}^{2}+\lambda^{-1}\|\nabla \tilde{e}_{p}\|_{L^{2}(\Omega)}^{2}\leq \|f-\lambda^{-1}p_{\sigma}+\Delta y_{\theta}+\lambda^{-1}\zeta^{ p}\|_{L^{2}(\Omega)}\|\tilde{e}_{y}\|_{L^{2}(\Omega)}\] \[+\lambda^{-1}\|y_{\theta}-y_{d}+\Delta p_{\sigma}-\zeta^{y}\|_{L^ {2}(\Omega)}\|\tilde{e}_{p}\|_{L^{2}(\Omega)}.\]

Then by the triangle inequality and the stability estimates in (4.2), we deduce

\[\|\nabla\tilde{e}_{y}\|_{L^{2}(\Omega)}^{2}+\lambda^{-1}\|\nabla \tilde{e}_{p}\|_{L^{2}(\Omega)}^{2}\leq \|f-\lambda^{-1}p_{\sigma}+\Delta y_{\theta}\|_{L^{2}(\Omega)}\| \tilde{e}_{y}\|_{L^{2}(\Omega)}+c\lambda^{-1}\|p_{\sigma}\|_{L^{2}(\partial \Omega)}\|\tilde{e}_{y}\|_{L^{2}(\Omega)}\] \[+\lambda^{-1}\|y_{\theta}-y_{d}+\Delta p_{\sigma}\|_{L^{2}(\Omega )}\|\tilde{e}_{p}\|_{L^{2}(\Omega)}+c\lambda^{-1}\|y_{\theta}\|_{L^{2}(\partial \Omega)}\|\tilde{e}_{p}\|_{L^{2}(\Omega)}.\]

Now by Poincare inequality, trace inequality, and Young's inequality, we obtain

\[\|\tilde{e}_{y}\|_{L^{2}(\Omega)}^{2}+\lambda^{-1}\|\tilde{e}_{p} \|_{L^{2}(\Omega)}^{2}\leq c\big{(}\|f-\lambda^{-1}p_{\sigma}+\Delta y_{\theta}\|_{L^{2}( \Omega)}^{2}+\lambda^{-1}\|y_{\theta}-y_{d}+\Delta p_{\sigma}\|_{L^{2}(\Omega )}^{2}\] \[+\lambda^{-1}\|y_{\theta}\|_{L^{2}(\partial\Omega)}^{2}+\lambda^{ -1}\|p_{\sigma}\|_{L^{2}(\partial\Omega)}^{2}\big{)}. \tag{4.4}\]

Meanwhile, by the triangle inequality, we have

\[\|e_{y}\|_{L^{2}(\Omega)}^{2}\leq c(\|\tilde{e}_{y}\|_{L^{2}(\Omega)}^{2}+\|\zeta^{y}\|_{L^{2}( \Omega)}^{2})\leq c(\|\tilde{e}_{y}\|_{L^{2}(\Omega)}^{2}+\|y_{\theta}\|_{L^{2}( \partial\Omega)}^{2})\leq c(\mathbf{\alpha},\lambda)\mathcal{L}(y_{\theta},p_{\sigma}),\] \[\|e_{p}\|_{L^{2}(\Omega)}^{2}\leq c(\|\tilde{e}_{p}\|_{L^{2}(\Omega)}^{2}+\|\zeta^{p}\|_{L^{2}( \Omega)}^{2})\leq c(\|\tilde{e}_{p}\|_{L^{2}(\Omega)}^{2}+\|p_{\sigma}\|_{L^{2}( \partial\Omega)}^{2})\leq c(\mathbf{\alpha},\lambda)\mathcal{L}(y_{\theta},p_{\sigma}).\]

Moreover, for the error \(e_{u}=\bar{u}-u_{\sigma}\), from the identity \(\|e_{u}\|_{L^{2}(\Omega)}=\lambda^{-1}\|e_{p}\|_{L^{2}(\Omega)}\), we deduce \(\|e_{u}\|_{L^{2}(\Omega)}^{2}\leq c(\mathbf{\alpha},\lambda)\mathcal{L}(y_{\theta},p_ {\sigma})\). Combining the preceding three estimates completes the proof of the lemma. 

**Remark 4.1**.: _It follows from the estimate (4.4) that one should choose \(\alpha_{i}\) as \(\lambda^{-1}\). The boundary penalty parameters \(\alpha_{b}^{y}\) and \(\alpha_{b}^{p}\) should be scaled accordingly._

Next, we discuss the constrained case, which is more involved, due to the presence of the nonlinear projection operator \(P_{U}\). Thus, we shall employ a different proof technique, based on the idea of reconstruction.

**Lemma 4.2**.: _Let \((\bar{y},\bar{p})\) be the solution tuple of the system (3.9), with \(\bar{u}=P_{U}(-\lambda^{-1}\bar{p})\). Then for any \((y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}\), and \(u_{\sigma}=P_{U}(-\lambda^{-1}p_{\sigma})\), the following estimate holds_

\[\|\bar{y}-y_{\theta}\|_{L^{2}(\Omega)}^{2}+\|\bar{p}-p_{\sigma}\|_{L^{2}(\Omega) }^{2}+\|\bar{u}-u_{\sigma}\|_{L^{2}(\Omega)}^{2}\leq c(\mathbf{\alpha},\lambda) \mathcal{L}(y_{\theta},p_{\sigma}).\]
Proof.: The proof employs crucially the following auxiliary functions \(R\bar{y},R\bar{p}\in H^{1}_{0}(\Omega)\) satisfying

\[\begin{cases}-\Delta R\bar{y}=f+u_{\sigma},&\text{in }\Omega,\\ R\bar{y}=0,&\text{on }\partial\Omega,\end{cases}\quad\text{and}\quad\begin{cases}- \Delta R\bar{p}=y_{\theta}-y_{d},&\text{in }\Omega,\\ R\bar{p}=0,&\text{on }\partial\Omega.\end{cases} \tag{4.5}\]

Then subtracting (4.5) from (3.4) yields

\[\begin{cases}-\Delta(\bar{y}-R\bar{y})=\bar{u}-u_{\sigma},&\text{in }\Omega,\\ \bar{y}-R\bar{y}=0,&\text{on }\partial\Omega,\end{cases} \tag{4.6}\]

and likewise

\[\begin{cases}-\Delta(\bar{p}-R\bar{p})=\bar{y}-y_{\theta},&\text{in }\Omega,\\ \bar{p}-R\bar{p}=0,&\text{on }\partial\Omega.\end{cases} \tag{4.7}\]

Now multiplying (4.6) by \((\bar{p}-R\bar{p})\) and (4.7) by \((\bar{y}-R\bar{y})\), integrating over the domain \(\Omega\), and equating the resulting identities yield

\[(\bar{p}-R\bar{p},\bar{u}-u_{\sigma})=(\bar{y}-R\bar{y},\bar{y}-y_{\theta}). \tag{4.8}\]

It follows from the first-order necessary optimality condition in (3.8) that

\[(\bar{p}+\lambda\bar{u},v-\bar{u})\geq 0,\quad\forall v\in U. \tag{4.9}\]

Also, since \(u_{\sigma}=P_{U}(\lambda^{-1}p_{\sigma})\in U\), we have

\[(p_{\sigma}+\lambda u_{\sigma},v-u_{\sigma})\geq 0,\quad\forall v\in U. \tag{4.10}\]

Putting \(v=u_{\sigma}\in U\) in (4.9), \(v=\bar{u}\in U\) in (4.10) and adding the resulting inequalities give

\[(\bar{p}-p_{\sigma},u_{\sigma}-\bar{u})-\lambda\left\|\bar{u}-u_{\sigma} \right\|_{L^{2}(\Omega)}^{2}\geq 0. \tag{4.11}\]

Meanwhile, it follows from the identity (4.8) that

\[(\bar{p}-p_{\sigma},u_{\sigma}-\bar{u}) =(\bar{p}-R\bar{p},u_{\sigma}-\bar{u})+(R\bar{p}-p_{\sigma},u_{ \sigma}-\bar{u})\] \[=-(\bar{y}-R\bar{y},\bar{y}-y_{\theta})+(R\bar{p}-p_{\sigma},u_{ \sigma}-\bar{u})\] \[=-\left\|\bar{y}-R\bar{y}\right\|_{L^{2}(\Omega)}^{2}-(\bar{y}-R \bar{y},R\bar{y}-y_{\theta})+(R\bar{p}-p_{\sigma},u_{\sigma}-\bar{u}).\]

Together with the inequality (4.11) and Cauchy-Schwarz inequality, we obtain

\[\left\|\bar{y}-R\bar{y}\right\|_{L^{2}(\Omega)}^{2}+\lambda\left\| \bar{u}-u_{\sigma}\right\|_{L^{2}(\Omega)}^{2} \leq(\bar{y}-R\bar{y},y_{\theta}-R\bar{y})+(R\bar{p}-p_{\sigma},u _{\sigma}-\bar{u})\] \[\leq\|\bar{y}-R\bar{y}\|_{L^{2}(\Omega)}\|y_{\theta}-R\bar{y}\|_ {L^{2}(\Omega)}+\|R\bar{p}-p_{\sigma}\|_{L^{2}(\Omega)}\|u_{\sigma}-\bar{u}\|_ {L^{2}(\Omega)}.\]

Then an application of Young's inequality leads to

\[\left\|\bar{y}-R\bar{y}\right\|_{L^{2}(\Omega)}^{2}+\lambda\left\|\bar{u}-u_{ \sigma}\right\|_{L^{2}(\Omega)}^{2}\leq\|y_{\theta}-R\bar{y}\|_{L^{2}(\Omega) }^{2}+\lambda^{-1}\|R\bar{p}-p_{\sigma}\|_{L^{2}(\Omega)}^{2}.\]

By the triangle inequality, we have

\[\left\|\bar{y}-y_{\theta}\right\|_{L^{2}(\Omega)}\leq c(\lambda)\big{(}\|R \bar{y}-y_{\theta}\|_{L^{2}(\Omega)}+\|R\bar{p}-p_{\sigma}\|_{L^{2}(\Omega)} \big{)}. \tag{4.12}\]

By the elliptic regularity theory for problem (4.7), cf. (4.2), we have \(\left\|\bar{p}-R\bar{p}\right\|_{L^{2}(\Omega)}\leq c\left\|\bar{y}-y_{\theta} \right\|_{L^{2}(\Omega)}\). Then using (4.12), we obtain the following estimate

\[\left\|\bar{p}-p_{\sigma}\right\|_{L^{2}(\Omega)} \leq\left\|\bar{p}-R\bar{p}\right\|_{L^{2}(\Omega)}+\|R\bar{p}-p_{ \sigma}\|_{L^{2}(\Omega)}\] \[\leq c\left\|\bar{y}-y_{\theta}\right\|_{L^{2}(\Omega)}+\|R\bar{p}- p_{\sigma}\|_{L^{2}(\Omega)}\] \[\leq c(\lambda)\big{(}\|R\bar{y}-y_{\theta}\|_{L^{2}(\Omega)}+\|R \bar{p}-p_{\sigma}\|_{L^{2}(\Omega)}\big{)}.\]
Next, let \(e_{y}^{R}=R\bar{y}-y_{\theta}\) and \(\tilde{e}_{y}^{R}=e_{y}^{R}+\zeta^{y}\), and similarly, let \(e_{p}^{R}=R\bar{p}-p_{\sigma}\) and \(\tilde{e}_{p}^{R}=e_{p}^{R}+\zeta^{p}\), where \(\zeta^{y}\) and \(\zeta^{p}\) are harmonic extensions defined in (4.1). Clearly, \(\tilde{e}_{y}^{R}\) and \(\tilde{e}_{p}^{R}\) satisfy

\[\begin{cases}-\Delta\tilde{e}_{y}^{R}=f+u_{\sigma}+\Delta y_{\theta},&\text{in } \Omega,\\ \tilde{e}_{y}^{R}=0,&\text{on }\partial\Omega,\end{cases}\quad\text{and}\quad \begin{cases}-\Delta\tilde{e}_{p}^{R}=y_{\theta}-y_{d}+\Delta p_{\sigma},& \text{in }\Omega,\\ \tilde{e}_{p}^{R}=0,&\text{on }\partial\Omega.\end{cases} \tag{4.13}\]

By the standard elliptic regularity estimate, we have

\[\left\|\tilde{e}_{y}^{R}\right\|_{L^{2}(\Omega)}\leq c\left\|f+u_{\sigma}+ \Delta y_{\theta}\right\|_{L^{2}(\Omega)}\quad\text{and}\quad\left\|\tilde{e }_{p}^{R}\right\|_{L^{2}(\Omega)}\leq c\left\|y_{\theta}-y_{d}+\Delta p_{ \sigma}\right\|_{L^{2}(\Omega)}.\]

Hence, we have

\[\left\|R\bar{y}-y_{\theta}\right\|_{L^{2}(\Omega)} =\left\|e_{y}^{R}\right\|_{L^{2}(\Omega)}=\left\|\tilde{e}_{y}^{ R}-\zeta^{y}\right\|_{L^{2}(\Omega)}\leq\left\|\tilde{e}_{y}^{R}\right\|_{L^{2}( \Omega)}+\left\|\zeta^{y}\right\|_{L^{2}(\Omega)}\] \[\leq c\left\|f+u_{\sigma}+\Delta y_{\theta}\right\|_{L^{2}(\Omega )}+c\left\|y_{\theta}\right\|_{L^{2}(\partial\Omega)},\] \[\left\|R\bar{p}-p_{\sigma}\right\|_{L^{2}(\Omega)} =\left\|e_{p}^{R}\right\|_{L^{2}(\Omega)}=\left\|\tilde{e}_{p}^{ R}-\zeta^{p}\right\|_{L^{2}(\Omega)}\leq\left\|\tilde{e}_{p}^{R}\right\|_{L^{2}( \Omega)}+\left\|\zeta^{p}\right\|_{L^{2}(\Omega)}\] \[\leq c\left\|y_{\theta}-y_{d}+\Delta p_{\sigma}\right\|_{L^{2}( \Omega)}+c\left\|p_{\sigma}\right\|_{L^{2}(\partial\Omega)}.\]

Therefore, we obtain the desired estimates and complete the proof of the lemma. 

### Approximation and statistical errors

Now we bound the generalization error via offset Rademacher complexity, first introduced by Liang et al [26] to prove sharp bounds for a two-step star estimator, and recently further developed in [13] for a broad class of problem settings. It is a penalized version of Rademacher complexity which localizes \(\mathcal{F}\) adaptively according to the magnitude of \(f^{2}\).

**Definition 4.1**.: _Let \(\mathcal{F}\) be a class of measurable functions from \(\mathcal{X}\) to \(\mathbb{R}\), \(\mu_{X}\) a probability distribution of \(X\) and \(\mathbb{X}:=\{X_{i}\}_{i=1}^{n}\) i.i.d. samples draw from \(\mu_{X}\). Let \(\{\tau_{i}\}_{i=1}^{n}\) be i.i.d. Rademacher random variables, i.e., \(\mathbb{P}(\tau_{i}=1)=\mathbb{P}(\tau_{i}=-1)=\frac{1}{2}\). Then the empirical offset Rademacher complexity of \(\mathcal{F}\) is defined as_

\[\mathcal{R}_{n}^{\rm off}(\mathcal{F},\beta|\mathbb{X}):=\mathbb{E}_{\tau} \Big{[}\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}f(X_{i})-\beta f (X_{i})^{2}|\mathbb{X}\Big{]},\]

_for some \(\beta>0\), and the offset Rademacher complexity of \(\mathcal{F}\) is defined by_

\[\mathcal{R}_{n}^{\rm off}(\mathcal{F},\beta):=\mathbb{E}_{\mathbb{X}}\mathcal{ R}_{n}^{\rm off}(\mathcal{F},\beta|\mathbb{X})=\mathbb{E}_{\mathbb{X},\tau} \Big{[}\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}f(X_{i})-\beta f (X_{i})^{2}\Big{]}.\]

First, we give a fundamental estimate of the generalization error in terms of the approximation error and statistical error (using offset Rademacher complexity). We define two sets \(\mathcal{G}_{d}\) and \(\mathcal{G}_{b}\) of mappings from \(X\in\Omega\mapsto\mathbb{R}\) and \(Y\in\partial\Omega\mapsto\mathbb{R}\), respectively, by

\[\mathcal{G}_{d} =\{|\Omega|((\Delta y+f-\lambda^{-1}p)(X))^{2}+\alpha_{i}|\Omega|( (\Delta p+y-y_{d})(X))^{2}:y\in\mathcal{Y},p\in\mathcal{P}\},\] \[\mathcal{G}_{b} =\{\alpha_{b}^{p}|\partial\Omega|y(Y)^{2}+\alpha_{b}^{p}|\partial \Omega|p(Y)^{2}:y\in\mathcal{Y},p\in\mathcal{P}\}.\]

Below we write \(g_{d}(y,p,X)=|\Omega|((\big{(}\Delta y+f-\lambda^{-1}p\big{)}(X))^{2}+\alpha_{i} ((\Delta p+y-y_{d})(X))^{2})\in\mathcal{G}_{d}\), and \(g_{b}(y,p,Y)=|\partial\Omega|(\alpha_{b}^{y}y(Y)^{2}+\alpha_{b}^{p}p(Y)^{2})\in \mathcal{G}_{b}\), and let \(b_{d}=\sup_{g\in\mathcal{G}_{d}}\left\|g\right\|_{L^{\infty}(\Omega)}\) and \(b_{b}=\sup_{g\in\mathcal{G}_{b}}\left\|g\right\|_{L^{\infty}(\partial\Omega)}\).

**Theorem 4.2**.: _For any minimizer \((\hat{y},\hat{p})\in\mathcal{Y}\times\mathcal{P}\) of the empirical loss \(\widehat{\mathcal{L}}(y,p)\), the following estimate holds_

\[\mathbb{E}_{\mathbb{X},\mathbb{Y}}[\mathcal{L}(\hat{y},\hat{p})]\leq 4\mathcal{R}_{n_{ d}}^{\rm off}\big{(}\mathcal{G}_{d},(2b_{d})^{-1}\big{)}+4\mathcal{R}_{n_{b}}^{\rm off }\big{(}\mathcal{G}_{b},(2b_{b})^{-1}\big{)}+3\inf_{(y,p)\in\mathcal{Y}\times \mathcal{P}}\mathcal{L}(y,p).\]
Proof.: Clearly, we have the following decomposition for any \((y,p)\in\mathcal{Y}\times\mathcal{P}\),

\[\mathbb{E}_{\mathbb{X},\mathbb{Y}}[\mathcal{L}(\hat{y},\hat{p})]= \mathbb{E}_{\mathbb{X},\mathbb{Y}}\left[\mathbb{E}_{X,Y}\left[g_{d}(\hat{y},\hat{p},X)+g_{b}(\hat{y},\hat{p},Y)\right]\right]\] \[= \mathbb{E}_{\mathbb{X},\mathbb{Y}}\Big{[}\mathbb{E}_{X,Y}\left[g_ {d}(\hat{y},\hat{p},X)+g_{b}(\hat{y},\hat{p},Y)\right]-\frac{3}{n_{d}}\sum_{i= 1}^{n_{d}}g_{d}(\hat{y},\hat{p},X_{i})-\frac{3}{n_{b}}\sum_{i=1}^{n_{b}}g_{b}( \hat{y},\hat{p},Y_{i})\Big{]}\] \[+\mathbb{E}_{\mathbb{X},\mathbb{Y}}\Big{[}\frac{3}{n_{d}}\sum_{i =1}^{n_{d}}g_{d}(\hat{y},\hat{p},X_{i})+\frac{3}{n_{b}}\sum_{i=1}^{n_{d}}g_{b} (\hat{y},\hat{p},Y_{i})\Big{]}:=\mathrm{I}+\mathrm{II}.\]

Since \((\hat{y},\hat{p})\in\mathcal{Y}\times\mathcal{P}\) is an empirical risk minimizer, i.e., \(\widehat{\mathcal{L}}(\hat{y},\hat{p})\leq\widehat{\mathcal{L}}(y,p)\), for any \((y,p)\in\mathcal{Y}\times\mathcal{P}\), we deduce

Meanwhile, the convexity of supremum and Jensen's inequality imply

\[\mathrm{I}\leq \sup_{(y,p)\in\mathcal{Y}\times\mathcal{P}}\mathbb{E}_{\mathbb{X}, \mathbb{Y}}\Big{[}\mathbb{E}_{X,Y}[g_{d}(y,p,X)+g_{b}(y,p,Y)]-\frac{3}{n_{d}} \sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})-\frac{3}{n_{b}}\sum_{i=1}^{n_{b}}g_{b}(y,p, Y_{i})\Big{]}\] \[\leq \mathbb{E}_{\mathbb{X},\mathbb{Y}}\sup_{(y,p)\in\mathcal{Y}\times \mathcal{P}}\Big{[}\mathbb{E}_{X}[g_{d}(y,p,X)]-\frac{3}{n_{d}}\sum_{i=1}^{n_{ d}}g_{d}(y,p,X_{i})+\mathbb{E}_{Y}[g_{b}(y,p,Y)]-\frac{3}{n_{b}}\sum_{i=1}^{n_{b}}g_{b} (y,p,Y_{i})\Big{]}\] \[\leq \mathbb{E}_{\mathbb{X}}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P }}\Big{[}\mathbb{E}_{X}[g_{d}(y,p,X)]-\frac{3}{n_{d}}\sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})\Big{]}+\mathbb{E}_{\mathbb{Y}}\sup_{(y,p)\in\mathcal{Y}\times \mathcal{P}}\Big{[}\mathbb{E}_{Y}[g_{b}(y,p,Y)]-\frac{3}{n_{b}}\sum_{i=1}^{n_{ b}}g_{b}(y,p,Y_{i})\Big{]}.\]

Next we bound the two terms, denoted by \(\mathrm{I}_{d}\) and \(\mathrm{I}_{b}\). Note that \(0\leq g_{d}(\cdot,\cdot,X)\leq b_{d}\) for all \(X\in\Omega\), which implies \(g_{d}{}^{2}(\cdot,\cdot,X)\leq b_{d}g_{d}(\cdot,\cdot,X)\). Consequently,

\[\mathrm{I}_{d}\leq \mathbb{E}_{\mathbb{X}}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P }}\Big{[}2\mathbb{E}_{X}g_{d}(y,p,X)-\frac{1}{b_{d}}\mathbb{E}_{X}g_{d}(y,p,X) ^{2}-\frac{2}{n_{d}}\sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})-\frac{1}{b_{d}n_{d}} \sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})^{2}\Big{]}.\]

Now we introduce independent copies of \(\mathbb{X}\) and \(\mathbb{Y}\), i.e., \(\mathbb{X}^{\prime}:=\{\mathbb{X}^{\prime}_{i}\}_{i=1}^{n_{d}}\) and \(\mathbb{Y}^{\prime}=\{Y^{\prime}_{i}\}_{i=1}^{n_{b}}\), and let \(\tau=\{\tau_{i}\}_{i=1}^{n_{d}}\) be a sequence of i.i.d. Rademacher random variables independent of \(\mathbb{X}\) and \(\mathbb{X}^{\prime}\). By the technique of symmetrization, the convexity of supremum, and Jensen's inequality, we obtain

\[\mathrm{I}_{d}\leq \mathbb{E}_{\mathbb{X}}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P }}\Big{[}\mathbb{E}_{\mathbb{X}^{\prime}}\Big{[}\frac{2}{n_{d}}\sum_{i=1}^{n_{d} }g_{d}(y,p,X^{\prime}_{i})-\frac{1}{b_{d}n_{d}}\sum_{i=1}^{n_{d}}g_{d}(y,p,X^{ \prime}_{i})^{2}\Big{]}-\frac{2}{n_{d}}\sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})-\frac{ 1}{b_{d}n_{d}}\sum_{i=1}^{n_{d}}g_{d}(y,p,X_{i})^{2}\Big{]}\] \[\leq \mathbb{E}_{\mathbb{X}}\mathbb{E}_{\mathbb{X}^{\prime}}\sup_{(y,p )\in\mathcal{Y}\times\mathcal{P}}\Big{[}\frac{2}{n_{d}}\sum_{i=1}^{n_{d}}(g_{d} (y,p,X^{\prime}_{i})-g_{d}(y,p,X_{i}))-\frac{1}{b_{d}n_{d}}\sum_{i=1}^{n_{d}}(g_ {d}(y,p,X^{\prime}_{i})^{2}+g_{d}(y,p,X_{i})^{2})\Big{]}\] \[= \mathbb{E}_{\mathbb{X}}\mathbb{E}_{\mathbb{X}^{\prime}}\mathbb{E}_{ \tau}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P}}\Big{[}\frac{2}{n_{d}}\sum_{i=1}^{n _{d}}\tau_{i}(g_{d}(y,p,X^{\prime}_{i})-g_{d}(y,p,X_{i}))-\frac{1}{b_{d}n_{d}} \sum_{i=1}^{n_{d}}(g_{d}(y,p,X^{\prime}_{i})^{2}+g_{d}(y,p,X_{i})^{2})\Big{]}\] \[= 2\mathbb{E}_{\mathbb{X}^{\prime}}\mathbb{E}_{\tau}\sup_{(y,p)\in \mathcal{Y}\times\mathcal{P}}\frac{1}{n_{d}}\sum_{i=1}^{n_{d}}\Big{(}\tau_{i}g_{ d}(y,p,X^{\prime}_{i})-\frac{1}{2b_{d}}g_{d}(y,p,X^{\prime}_{i})^{2}\Big{)}\] \[= 2\mathbb{E}_{\mathbb{X}^{\prime}}\mathcal{R}^{\mathrm{off}}_{n} \big{(}\mathcal{G}_{d},(2b_{d})^{-1}|\mathbb{X}^{\prime}\big{)}+2\mathbb{E}_{ \mathbb{X}}\mathcal{R}^{\mathrm{off}}_{n_{d}}\big{(}\mathcal{G}_{d},(2b_{d})^{-1} |\mathbb{X}\big{)}=4\mathcal{R}^{\mathrm{off}}_{n_{d}}\big{(}\mathcal{G}_{d},(2b_ {d})^{-1}\big{)}.\]

Similarly, we can deduce \(\mathrm{I}_{b}\leq 4\mathcal{R}^{\mathrm{off}}_{n_{b}}\big{(}\mathcal{G}_{b},(2b_{b})^{-1} \big{)}\). This completes the proof of the theorem. 

Theorem 4.2 indicates that the total error can be decomposed into the approximation error \(\mathcal{E}_{app}=3\inf_{(y,p)\in\mathcal{Y}\times\mathcal{P}}\mathcal{L}(y,p)\) and statistical error \(\mathcal{E}_{stat}:=4\mathcal{R}^{\mathrm{off}}_{n_{d}}(\mathcal{G}_{d},(2b_{d}) ^{-1})+4\mathcal{R}^{\mathrm{off}}_{n_{b}}(\mathcal{G}_{b},(2b_{b})^{-1})\) characterized by
the offset Rademacher complexity. The former arises from restricting the trial spaces of \((y,p)\) from the whole space \(H^{1}_{0}(\Omega)\cap H^{2}(\Omega)\) to \(\mathcal{Y}\times\mathcal{P}\), whereas the latter is due to approximating the integrals with Monte Carlo methods. We have ignored the optimization error: the optimizer only finds an approximate local minimizer of the loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\), due to its highly nonconvex landscape. The analysis of the optimization error is largely open, and hence it is not studied below. Next we bound the approximation error \(\mathcal{E}_{app}\) and statistical error \(\mathcal{E}_{stat}\).

#### 4.2.1 Approximation Error

We recall the following result on the approximation error [15, Proposition 4.8].

**Proposition 4.1**.: _Let \(p\geq 1,s,k,d\in\mathbb{N}\cup\{0\},s\geq k+1\), \(\rho\) be the logistic or \(\tanh\) function, and fix \(\mu>0\) small. Then for any \(\epsilon>0\), \(f\in W^{s,p}([0,1]^{d})\) with \(\|f\|_{W^{s,p}([0,1]^{d})}\leq 1\), there exists a DNN \(f_{\theta}\in\mathcal{N}_{\rho}(c\log(d+s),\,c(d,s,p,k)\epsilon^{-\frac{d}{s- k-\mu}},\,c(d,s,p,k)\epsilon^{-2-\frac{2(d/p+d+k+\mu)+d/p+d}{s-k-\mu}})\) such that \(\|f-f_{\theta}\|_{W^{k,p}([0,1]^{d})}\leq\epsilon\)._

Then we have the following bound on the approximation error \(\mathcal{E}_{app}\).

**Lemma 4.3**.: _Fix a tolerance \(\epsilon\in(0,1)\), and \(\mu>0\) arbitrarily small. If the optimal state \(\bar{y}\in H^{s}(\Omega)\cap H^{1}_{0}(\Omega)\) and adjoint \(\bar{p}\in H^{s}(\Omega)\cap H^{1}_{0}(\Omega)\), \(s\geq 3\), then there exist DNNs \((y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}\) with \(\mathcal{Y}=\mathcal{P}=\mathcal{N}_{\rho}(c\log(d+s),c(d,s)\epsilon^{-\frac{ d}{s-2-\mu}},c(d,s)\epsilon^{-\frac{d+s}{2(s-2-\mu)}})\) and \(c>0\) dependent on \(\alpha\), \(\lambda\), \(s\), \(\|\bar{y}\|_{H^{s}(\Omega)}\) and \(\|\bar{p}\|_{H^{s}(\Omega)}\) such that the approximation error \(\mathcal{E}_{app}\leq c\epsilon^{2}\)._

Proof.: For any \((y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}\), since \(\mathcal{L}(\bar{y},\bar{p})=0\), we deduce

\[\mathcal{L}(y_{\theta},p_{\sigma})-\mathcal{L}(\bar{y},\bar{p})= \mathcal{L}(y_{\theta},p_{\sigma})\] \[= \|\Delta y_{\theta}+f-\lambda^{-1}p_{\sigma}\|_{L^{2}(\Omega)} ^{2}+\alpha_{i}\|\Delta p_{\sigma}+y_{\theta}-y_{d}\|_{L^{2}(\Omega)}^{2}+ \alpha_{b}^{y}\|y_{\theta}\|_{L^{2}(\partial\Omega)}^{2}+\alpha_{b}^{p}\|p_{ \sigma}\|_{L^{2}(\partial\Omega)}^{2}\] \[= \|\Delta(y_{\theta}-\bar{y})+\lambda^{-1}(\bar{p}-p_{\sigma})\|_ {L^{2}(\Omega)}^{2}+\alpha_{i}\|\Delta(p_{\sigma}-\bar{p})+y_{\theta}-\bar{y} \|_{L^{2}(\Omega)}^{2}+\alpha_{b}^{y}\|y_{\theta}-\bar{y}\|_{L^{2}(\partial \Omega)}^{2}+\alpha_{b}^{p}\|p_{\sigma}-\bar{p}\|_{L^{2}(\partial\Omega)}^{2}\] \[\leq c(\boldsymbol{\alpha},\lambda)\big{(}\|\bar{y}-y_{\theta}\|_{H^ {2}(\Omega)}^{2}+\|\bar{p}-p_{\sigma}\|_{H^{2}(\Omega)}^{2}\big{)},\]

by the trace inequality \(\|y_{\theta}-\bar{y}\|_{L^{2}(\partial\Omega)}\leq c\|y_{\theta}-\bar{y}\|_{H^ {2}(\Omega)}\). With \(c_{y}=\|\bar{y}\|_{H^{s}(\Omega)}\), we have

\[\inf_{y_{\theta}\in\mathcal{Y}}\|\bar{y}-y_{\theta}\|_{H^{2}(\Omega)}^{2}=c_{ y}^{2}\inf_{y_{\theta}\in\mathcal{Y}}\Big{\|}\frac{\bar{y}}{c_{y}}-\frac{y_{ \theta}}{c_{y}}\Big{\|}_{H^{2}(\Omega)}^{2}=c_{y}^{2}\inf_{y_{\theta}\in \mathcal{Y}}\Big{\|}\frac{\bar{y}}{c_{y}}-y_{\theta}\Big{\|}_{H^{2}(\Omega)}^{2}.\]

By Proposition 4.1, there exists a DNN \(y_{\theta}\in\mathcal{Y}=\mathcal{N}_{\rho}(c\log(d+s),c(d,s)\epsilon^{- \frac{d}{s-2-\mu}},c(d,s)\epsilon^{-\frac{0d+4\epsilon}{2(s-2-\mu)}})\) such that \(\|\frac{\bar{y}}{c_{y}}-y_{\theta}\|_{H^{2}(\Omega)}\leq\epsilon\). Similarly, there exists a DNN \(p_{\sigma}\in\mathcal{P}\) with the requisite property. Hence,

\[\mathcal{E}_{app}=3\inf_{(y_{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P} }\mathcal{L}(y_{\theta},p_{\sigma})\leq c(\boldsymbol{\alpha},\lambda)\inf_{(y _{\theta},p_{\sigma})\in\mathcal{Y}\times\mathcal{P}}(\|\bar{y}-y_{\theta}\|_{H ^{2}(\Omega)}^{2}+\|\bar{p}-p_{\sigma}\|_{H^{2}(\Omega)}^{2})\leq c\epsilon^{2},\]

where \(c\) depends on \(\boldsymbol{\alpha}\), \(\lambda\), \(s\) and the \(H^{s}(\Omega)\) norms of \(\bar{y}\) and \(\bar{p}\). Then the desired assertion follows. 

**Remark 4.3**.: _The approximation error \(\mathcal{E}_{app}\) depends on the regularity of \(\bar{y}\) and \(\bar{u}\). If \(f,y_{d}\in H^{s-2}(\Omega)\cap L^{\infty}(\Omega)\), by standard elliptic regularity theory, we have \(\bar{y},\bar{p}\in H^{s}(\Omega)\), which satisfies the requirements of Lemma 4.3. In the unconstrained case, the optimal control \(\bar{u}\) has the same regularity as the adjoint state \(\bar{p}\), but in the constrained case, we have \(\bar{u}\in H^{1}(\Omega)\) only, due to the presence of the projection operator \(P_{U}\)._

#### 4.2.2 Statistical error

To bound the statistical error via the offset Rademacher complexities, we employ the covering number [1].

**Definition 4.2**.: _Let \((\mathcal{M},d)\) be a metric space, and the set \(\mathcal{F}\subset\mathcal{M}\). A set \(\mathcal{F}_{\delta}\subset\mathcal{M}\) is called a \(\delta\)-cover of \(\mathcal{F}\) if for each \(f\in\mathcal{F}\), there exists an \(f_{\delta}\in\mathcal{F}_{\delta}\) such that \(d(f,f_{\delta})\leq\delta\). Moreover, \(N(\delta,\mathcal{F},d):=\inf\{|\mathcal{F}_{\delta}|:\mathcal{F}_{\delta}\text{ is a $\delta$-cover of $\mathcal{F}$}\}\) is called the \(\delta\)-covering number of \(\mathcal{F}\)._

We also need Hoeffding's inequality [17].
**Lemma 4.4**.: _Let \(X_{1}\),..., \(X_{n}\) be independent random variables such that \(a_{i}\leq X_{i}\leq b_{i}\) almost surely. Consider the sum of these random variables, \(S_{n}=X_{1}+\cdots+X_{n}\). Then for all \(t>0\)_

\[\mathbb{P}\left(S_{n}-\mathbb{E}[S_{n}]\geq t\right)\leq\exp\left(-\frac{2t^{2} }{\sum_{i=1}^{n}(b_{i}-a_{i})^{2}}\right).\]

Next, we bound the offset Rademacher complexity in terms of the covering number.

**Theorem 4.4**.: _Let \(\kappa_{d}\) and \(\kappa_{b}\) be the Lipschitz constants of \(g_{d}\in\mathcal{G}_{d}\) and \(g_{b}\in\mathcal{G}_{b}\) in \(y\) and \(p\) in the \(W^{2,\infty}(\Omega)\) and \(L^{\infty}(\partial\Omega)\) norms, respectively, i.e.,_

\[|g_{d}(y,p,X)-g_{d}(\tilde{y},\tilde{p},X)| \leq\kappa_{d}(\|y-\tilde{y}\|_{W^{2,\infty}(\Omega)}+\|p-\tilde {p}\|_{W^{2,\infty}(\Omega)}),\quad\forall(y,p),(\tilde{y},\tilde{p})\in \mathcal{Y}\times\mathcal{P},X\in\Omega,\] \[|g_{b}(y,p,Y)-g_{b}(\tilde{y},\tilde{p},Y)| \leq\kappa_{b}(\|y-\tilde{y}\|_{L^{\infty}(\partial\Omega)}+\|p- \tilde{p}\|_{L^{\infty}(\partial\Omega)}),\quad\forall(y,p),(\tilde{y},\tilde{ p})\in\mathcal{Y}\times\mathcal{P},Y\in\partial\Omega.\]

_Then for any \(\delta>0\), the offset Rademacher complexities \(\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{d},\beta)\) and \(\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{b},\beta)\) are bounded by_

\[\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{d},\beta)\leq \frac{1+\log N(\delta,\mathcal{F}\times\mathcal{P},W^{2,\infty}( \Omega))}{2n\beta}+2(1+2b_{d}\beta)\kappa_{d}\delta,\] \[\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{b},\beta)\leq \frac{1+\log N(\delta,\mathcal{F}\times\mathcal{P},L^{\infty}( \partial\Omega))}{2n\beta}+2(1+2b_{b}\beta)\kappa_{b}\delta.\]

Proof.: We give the proof only for \(\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{d},\beta)\), since that for \(\mathcal{R}^{\rm off}_{n}(\mathcal{G}_{b},\beta)\) is identical. Since \(\tau=\{\tau_{i}\}_{i=1}^{n}\) is a sequence of i.i.d. random variables independent of \(\mathbb{X}\), then conditionally on \(\mathbb{X}\), we have

\[\mathbb{E}_{\tau|\mathbb{X}}\Big{[}\sup_{(y,p)\in\mathcal{Y}\times \mathcal{P}}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y,p,X_{i})-\frac{\beta}{n} \sum_{i=1}^{n}g_{d}(y,p,X_{i})^{2}|\mathbb{X}\Big{]}\] \[= \mathbb{E}_{\tau}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P}}\Big{[} \frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y,p,X_{i})-\frac{\beta}{n}\sum_{i=1}^{n} g_{d}(y,p,X_{i})^{2}\Big{]}.\]

Let \(\delta>0\) and let \(\mathcal{Y}_{\delta}\) and \(\mathcal{P}_{\delta}\) be a minimal \(W^{2,\infty}(\Omega)\)\(\delta\)-cover of \(\mathcal{Y}\) and \(\mathcal{P}\), respectively. For any fixed \((y,p)\in\mathcal{Y}\times\mathcal{P}\), there exists a tuple \((y_{\delta},p_{\delta})\in\mathcal{Y}_{\delta}\times\mathcal{P}_{\delta}\) such that \(\|y-y_{\delta}\|_{W^{2,\infty}(\Omega)}\leq\delta\) and \(\|p-p_{\delta}\|_{W^{2,\infty}(\Omega)}\leq\delta\). Therefore, by Lipschitz continuity of \(g_{d}\) with respect to \((y,p)\) in the \(W^{2,\infty}(\Omega)\) norm,

\[\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y,p,X_{i}) \leq\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y_{\delta},p_{\delta}, X_{i})+\frac{1}{n}\sum_{i=1}^{n}|\tau_{i}||(g_{d}(y,p,X_{i})-g_{d}(y_{\delta},q_{ \delta},X_{i})|\] \[\leq\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y_{\delta},p_{\delta}, X_{i})+2\kappa_{d}\delta.\]

Since \(|g_{d}(y,p,X)|\leq b_{d}\) and \(|g_{d}(y_{\delta},p_{\delta},X)|\leq b_{d}\), and by Lipschitz continuity of \(g_{d}(y,p,X)\), we obtain

\[-\frac{1}{n}\sum_{i=1}^{n}g_{d}(y,p,X_{i})^{2} =-\frac{1}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{\delta},X_{i})^{2} +\frac{1}{n}\sum_{i=1}^{n}(g_{d}(y_{\delta},p_{\delta},X_{i})^{2}-g_{d}(y,p,X_ {i})^{2})\] \[\leq-\frac{1}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{\delta},X_{i})^ {2}+\frac{2b_{d}}{n}\sum_{i=1}^{n}|g_{d}(y_{\delta},p_{\delta},X_{i})-g_{d}(y,p,X _{i})|\] \[\leq-\frac{1}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{\delta},X_{i})^ {2}+4b_{d}\kappa_{d}\delta.\]

Hence, it follows that

\[\mathbb{E}_{\tau}\sup_{(y,p)\in\mathcal{Y}\times\mathcal{P}}\Big{[}\frac{1}{n} \sum_{i=1}^{n}\tau_{i}g_{d}(y,p,X_{i})-\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y,p,X_ {i})^{2}\Big{]}\]
\[\leq \mathbb{E}_{\tau}\max_{(y_{\delta},p_{\delta})\in\mathcal{F}_{\delta} \times\mathcal{P}_{\delta}}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y_{ \delta},p_{\delta},X_{i})-\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{ \delta},X_{i})^{2}\Big{]}+2(1+2b_{d}\beta)\kappa_{d}\delta. \tag{4.14}\]

Since \(\{\tau_{i}g_{d}(y_{\delta},p_{\delta};X_{i})\}_{i=1}^{n}\) are independent random variables conditioning on \(\mathbb{X}\),

\[\mathbb{E}_{\tau}[\tau_{i}g_{d}(y_{\delta},p_{\delta},X_{i})] =0,\quad i=1,\ldots,n,\] \[-g_{d}(y_{\delta},p_{\delta},X_{i})\leq\tau_{i}g_{d}(y_{\delta},p _{\delta},X_{i}) \leq g_{d}(y_{\delta},p_{\delta},X_{i}),\quad i=1,\ldots,n.\]

By Hoeffding's inequality from Lemma 4.4, we deduce that for any \((y_{\delta},p_{\delta})\in\mathcal{Y}_{\delta}\times\mathcal{P}_{\delta}\) and \(\xi>0\),

\[\mathbb{P}_{\tau}\Big{\{}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y _{\delta},p_{\delta},X_{i})>\xi+\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta}, p_{\delta},X_{i})^{2}\Big{\}}\] \[\leq \exp\Big{(}-\frac{(n\xi+\beta\sum_{i=1}^{n}g_{d}(y_{\delta},p_{ \delta},X_{i})^{2})^{2}}{2\sum_{i=1}^{n}g_{d}(y_{\delta},p_{\delta},X_{i})^{2} }\Big{)}\leq\exp(-2\beta n\xi),\]

where the last step follows from the elementary inequality \(\frac{(a+y)^{2}}{y}\geq\frac{(a+a)^{2}}{a}=4a\), for any \(y\in\mathbb{R}_{+}\). Therefore, we may bound the tail probability by

\[\mathbb{P}_{\tau}\left\{\max_{(y_{\delta},p_{\delta})\in\mathcal{ Y}_{\delta}\times\mathcal{P}_{\delta}}\Big{(}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d} (y_{\delta},p_{\delta},X_{i})-\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p _{\delta},X_{i})^{2}\Big{)}>\xi\right\}\] \[\leq N(\delta,\mathcal{F}\times\mathcal{P},W^{2,\infty}(\Omega))\max_{ (y_{\delta},p_{\delta})\in\mathcal{Y}_{\delta}\times\mathcal{P}_{\delta}} \mathbb{P}_{\tau}\Big{\{}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y_{\delta},p _{\delta},X_{i})>\xi+\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{\delta}, X_{i})^{2}\Big{\}}\] \[\leq N(\delta,\mathcal{F}\times\mathcal{P},W^{2,\infty}(\Omega))\exp( -2\beta n\xi).\]

Hence, for any \(a>0\), we have

\[\mathbb{E}_{\tau}\Big{[}\max_{(y_{\delta},p_{\delta})\in\mathcal{ Y}_{\delta}\times\mathcal{P}_{\delta}}\Big{(}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d} (y_{\delta},p_{\delta},X_{i})-\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{ \delta},X_{i})^{2}\Big{)}\Big{]}\] \[\leq a+\int_{a}^{\infty}N(\delta,\mathcal{F}\times\mathcal{P},W^{2, \infty}(\Omega))\exp(-2\beta n\xi)\mathrm{d}\xi\leq a+\frac{N(\delta,\mathcal{F }\times\mathcal{P},W^{2,\infty}(\Omega))}{2\beta n}\exp(-2\beta na).\]

Setting \(a=\frac{\log N(\delta,\mathcal{F}\times\mathcal{P},W^{2,\infty}(\Omega))}{2\beta n}\) leads to

\[\mathbb{E}_{\tau}\Big{[}\max_{(y_{\delta},p_{\delta})\in\mathcal{Y}_{\delta} \times\mathcal{P}_{\delta}}\Big{(}\frac{1}{n}\sum_{i=1}^{n}\tau_{i}g_{d}(y_{ \delta},p_{\delta},X_{i})-\frac{\beta}{n}\sum_{i=1}^{n}g_{d}(y_{\delta},p_{ \delta},X_{i})^{2}\Big{]}\leq\frac{1+\log N(\delta,\mathcal{F}\times\mathcal{P},W^{2,\infty}(\Omega))}{2\beta n}.\]

Combining this inequality with the estimate (4.14) yields the desired assertion. 

The next result gives the statistical error, by exploiting the bound on the DNNs in Section 5.

**Theorem 4.5**.: _Let \(\mathcal{Y}=\mathcal{P}=\mathcal{N}_{\rho}(L,\mathbf{n}_{L},R)\), with depth \(L\), \(\mathbf{n}_{L}\) nonzero DNN parameters and maximum bound \(R\). Then the statistical error \(\mathcal{E}_{stat}\) is bounded by_

\[\mathcal{E}_{stat}\leq\frac{cL\mathbf{n}_{L}^{4L+1}R^{4L}(L\log R+L\log\mathbf{n}_{L}+ \log n_{d})}{n_{d}}+\frac{c\mathbf{n}_{L}^{3}R^{2}(L\log R+L\log\mathbf{n}_{L}+\log n_ {b})}{n_{b}},\]

_where the constant \(c\) depends on \(\mathbf{\alpha}\), \(\lambda\), \(d\), \(\|f\|_{L^{\infty}(\Omega)}\) and \(\|y_{d}\|_{L^{\infty}(\Omega)}\) at most polynomially._
Proof.: By Lemmas 5.1 and 5.2 below, for any \(v\in\mathcal{Y}\), \(\|v\|_{C(\overline{\Omega})}\leq\mathbf{n}_{L}R\) and \(\|\Delta v\|_{L^{\infty}(\Omega)}\leq dL\mathbf{n}_{L}^{2L}R^{2L}\), which directly give \(b_{d}=cL\mathbf{n}_{L}^{4L}R^{4L}\) and \(b_{b}=c(\mathbf{n}_{L}R)^{2}\). Next, for any \((y,p),(\tilde{y},\tilde{p})\in\mathcal{Y}\times\mathcal{P}\), we have

\[|g_{d}(y,p,X)-g_{d}(\tilde{y},\tilde{p},X)|\leq \kappa_{d}(\|y-\tilde{y}\|_{W^{2,\infty}(\Omega)}+\|p-\tilde{p}\|_{ W^{2,\infty}(\Omega)}),\quad\text{with }\kappa_{d}=c\mathbf{n}_{L}^{2L}R^{2L},\] \[|g_{b}(y,p,Y)-g_{b}(\tilde{y},\tilde{p},Y)|\leq \kappa_{b}(\|y-\tilde{y}\|_{L^{\infty}(\partial\Omega)}+\|p- \tilde{p}\|_{L^{\infty}(\partial\Omega)}),\quad\text{with }\kappa_{b}=c\mathbf{n}_{L}R.\]

Next we bound \(N(\delta,\mathcal{Y},W^{2,\infty}(\Omega))\). By Lemmas 5.1 and 5.3, for any \(v_{\theta},v_{\tilde{\theta}}\in\mathcal{Y}\),

\[\|v_{\theta}-v_{\tilde{\theta}}\|_{C(\overline{\Omega})} \leq\mathbf{n}_{L}^{L}R^{L-1}\|\theta-\tilde{\theta}\|_{\ell^{2}},\] \[\|\partial_{x_{i}x_{i}}^{2}v_{\theta}-\partial_{x_{i}x_{i}}^{2}v _{\tilde{\theta}}\|_{C(\overline{\Omega})} \leq 4L^{2}\mathbf{n}_{L}^{3L-2}R^{3L-3}\|\theta-\tilde{\theta}\|_{\ell^{ 2}},\quad i=1,\ldots,d.\]

Hence, with \(\Lambda_{d}=cL^{2}\mathbf{n}_{L}^{3L-2}R^{3L-3}\) and \(\Lambda_{b}=\mathbf{n}_{L}^{L}R^{L}\),

\[\|v_{\theta}-v_{\tilde{\theta}}\|_{W^{2,\infty}(\Omega)}\leq \Lambda_{d}\|\theta-\tilde{\theta}\|_{\ell^{2}}\quad\text{and}\quad\|v_{ \theta}-v_{\tilde{\theta}}\|_{L^{\infty}(\partial\Omega)}\leq\Lambda_{b}\| \theta-\tilde{\theta}\|_{\ell^{2}},\quad\forall v_{\theta},v_{\tilde{\theta}} \in\mathcal{Y}. \tag{4.15}\]

Note that for any \(m\in\mathbb{N}\), \(r\in[1,\infty)\), \(\epsilon\in(0,1)\), and \(B_{r}:=\{x\in\mathbb{R}^{m}:\;\|x\|_{\ell^{2}}\leq r\}\), then by a simple counting argument (see, e.g., [11, Proposition 5] or [22, Lemma 5.5]), we have \(\log N(\epsilon,B_{r},\|\cdot\|_{\ell^{2}})\leq m\log(2r\sqrt{m}\epsilon^{-1})\). This, Lipschitz continuity of NN functions in (4.15) and the estimate \(\|\theta\|_{\ell^{2}}\leq\sqrt{\mathbf{n}_{L}}\|\theta\|_{\ell^{\infty}}\leq\sqrt {\mathbf{n}_{L}}R\) imply

\[\log N(\delta,\mathcal{Y},W^{2,\infty}(\Omega)) \leq\log N(\Lambda_{d}^{-1}\delta,\mathcal{N}_{Y},\|\cdot\|_{\ell ^{2}})\leq\mathbf{n}_{L}\log(2R\mathbf{n}_{L}\Lambda_{d}\delta^{-1}),\] \[\log N(\delta,\mathcal{Y},L^{\infty}(\partial\Omega)) \leq\log N(\Lambda_{b}^{-1}\delta,\mathcal{N}_{Y},\|\cdot\|_{\ell ^{2}})\leq\mathbf{n}_{L}\log(2R\mathbf{n}_{L}\Lambda_{b}\delta^{-1}),\]

where \(\mathcal{N}_{Y}\) denotes the parameter space for \(\mathcal{Y}\). These estimates and Theorem 4.4 with \(\beta=(2b_{d})^{-1}\) yield

\[\mathcal{R}_{n}^{\text{off}}\big{(}\mathcal{G}_{d},(2b_{d})^{-1} \big{)} \leq\frac{b_{d}(1+\log N(\delta,\mathcal{F}\times\mathcal{P},W^{ 2,\infty}(\Omega)))}{n}+4\kappa_{d}\delta\] \[\leq\frac{cL\mathbf{n}_{L}^{4L+1}R^{4L}\log(2R\mathbf{n}_{L}\Lambda_{d} \delta^{-1})}{n}+c\mathbf{n}_{L}^{2L}R^{2L}\delta.\]

An analogous bound on \(\mathcal{R}_{n}^{\text{off}}\big{(}\mathcal{G}_{b},(2b_{b})^{-1}\big{)}\) holds. Setting \(\delta=1/n\), substituting \(\Lambda_{d}\) and \(\Lambda_{b}\) and simplifying the resulting expressions yield the desired estimate. 

#### 4.2.3 Final error estimate

Now we state the error of the C-PINN approximation \((y_{\theta^{*}},p_{\sigma^{*}},u_{\sigma^{*}})\). Thus, with the parameters in the loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) chosen suitably, C-PINN can yield an accurate approximation.

**Theorem 4.6**.: _Let the tuple \((\bar{y},\bar{p},\bar{u})\) solve the optimality system (3.3) / (3.8) such that \(\bar{y}\in H^{s}(\Omega)\cap H^{1}_{0}(\Omega)\) and \(\bar{p}\in H^{s}(\Omega)\cap H^{1}_{0}(\Omega)\) with \(s\geq 3\), and \((y_{\theta^{*}},p_{\sigma^{*}},u_{\sigma^{*}})\) be the C-PINN approximation. Fix the tolerance \(\epsilon>0\), and take \(\mathcal{Y}=\mathcal{P}=\mathcal{N}_{\rho}(c\log(d+s),c(d,s)\epsilon^{-\frac{ \theta d+s}{2-\mu}},c(d,s)\epsilon^{-\frac{\theta d+s}{2(\sigma-2\mu)}})\). Then by choosing \(n_{d}=c(d,s)\epsilon^{-\frac{c(d+s)\ln(d+s)}{s-2-\mu}}\) and \(n_{b}=c(d,s)\epsilon^{-\frac{\theta(3d+s)}{s-2-\mu}-2}\) sampling points in \(\Omega\) and on \(\partial\Omega\), there holds_

\[\mathbb{E}_{\mathbb{X},\mathbb{Y}}\left[\|\bar{y}-y_{\theta^{*}}\|_{L^{2}(\Omega)} ^{2}+\|\bar{p}-p_{\sigma^{*}}\|_{L^{2}(\Omega)}^{2}+\|\bar{u}-u_{\sigma^{*}}\|_ {L^{2}(\Omega)}^{2}\right]\leq c\epsilon^{2},\]

_where the constant \(c\) depends on \(\mathbf{\alpha}\), \(\lambda\), \(d\), \(s\), \(\|f\|_{L^{\infty}(\Omega)}\), \(\|y_{d}\|_{L^{\infty}(\Omega)}\), \(\|\bar{y}\|_{H^{s}(\Omega)}\), and \(\|\bar{p}\|_{H^{s}(\Omega)}\)._

Proof.: The fundamental estimates in Lemmas 4.1 and 4.2 imply

\[\mathbb{E}_{\mathbb{X},\mathbb{Y}}[\|\bar{y}-y_{\theta^{*}}\|_{L^{2}(\Omega)}^{2}+ \|\bar{p}-p_{\sigma^{*}}\|_{L^{2}(\Omega)}^{2}+\|\bar{u}-u_{\sigma^{*}}\|_{L^{2} (\Omega)}^{2}]\leq c(\mathbf{\alpha},\gamma)\mathbb{E}_{\mathbb{X},\mathbb{Y}}[ \mathcal{L}(y_{\theta^{*}},p_{\sigma^{*}})].\]

By the error decomposition in Theorem 4.2 and Lemma 4.3 (with the given choice of \(\mathcal{Y}\) and \(\mathcal{P}\)),

\[\mathbb{E}_{\mathbb{X},\mathbb{Y}}[\mathcal{L}(y_{\theta^{*}},p_{\sigma^{*}})] \leq c\epsilon^{2}+\frac{cL\mathbf{n}_{L}^{4L+1}R^{4L}(L\log R+L\log\mathbf{n}_{L}+ \log n_{d})}{n_{d}}+\frac{c\mathbf{n}_{L}^{3}R^{2}(L\log R+L\log\mathbf{n}_{L}+ \log n_{b})}{n_{b}}.\]

Since \(L=c\log(d+s)\), \(\mathbf{n}_{L}=c(d,s)\epsilon^{-\frac{d}{2-\mu}}\) and \(R=c(d,s)\epsilon^{-\frac{\theta d+s}{2(s-2-\mu)}}\), the choices of the numbers \(n_{d}\) and \(n_{b}\) imply that the statistical error \(\mathcal{E}_{stat}\) is also bounded by \(c\epsilon^{2}\). This completes the proof of the theorem.
## 5 Technical estimates on \(f_{\theta}\)

In this part, we derive several technical estimates, especially the bound and Lipschitz continuity of \(\partial_{x_{p}}^{2}f_{\theta}\) for \(f_{\theta}\in\mathcal{N}_{\rho}(L,\mathbf{n}_{L},R)\) in terms of the neural network parameters \(\theta\). First, we recall several bounds on \(f_{\theta}\) and \(\partial_{x_{p}}f_{\theta}\), which will be used extensively below. Throughout, we denote \(\pi_{i}=\prod_{j=1}^{i}n_{j}\), and let \(\mathcal{P}=\mathcal{N}_{\rho}(L,\mathbf{n}_{L},R)\).

**Lemma 5.1**.: _For any \(f\equiv f_{\theta},\tilde{f}\equiv f_{\tilde{\theta}}\in\mathcal{P}\), the following estimates hold: for any \(p\in[d],q\in[n_{\ell}]\),_

\[\|f\|_{L^{\infty}(\Omega)} \leq n_{L-1}R, \tag{5.1}\] \[\|f_{q}^{(\ell)}-\tilde{f}_{q}^{(\ell)}\|_{L^{\infty}(\Omega)} \leq\begin{cases}\quad\pi_{\ell-1}R^{\ell-1}\sum_{j=1}^{\mathbf{n}_{ \ell}}|\theta_{j}-\tilde{\theta}_{j}|,&\ell=1,\dots,L-1,\\ \sqrt{\mathbf{n}_{L}}\pi_{L-1}R^{L-1}\|\theta-\tilde{\theta}\|_{\ell^{2}},&\ell= L,\end{cases}\] (5.2) \[\|\partial_{x_{p}}f^{(\ell)}\|_{L^{\infty}(\Omega)} \leq\pi_{\ell-1}R^{\ell},\quad\ell=1,\dots,L,\] (5.3) \[\|\partial_{x_{p}}f_{q}^{(\ell)}-\partial_{x_{p}}\tilde{f}_{q}^{ (\ell)}\|_{L^{\infty}(\Omega)} \leq(\ell+1)\pi_{\ell-1}^{2}R^{2\ell-1}\sum_{j=1}^{\mathbf{n}_{\ell} }|\theta_{j}-\tilde{\theta}_{j}|,\quad\ell=1,\dots,L-1. \tag{5.4}\]

Proof.: These estimates are contained in [22, Lemmas 5.9-5.11]. Note that the estimate (5.4) improves that in [22, Lemma 5.11] by a factor of \(R\), by slightly improving the bound on page 17, line 7 of [22]. 

We also need the following uniform bound on \(\partial_{x_{p}}^{2}f_{\theta}\).

**Lemma 5.2**.: _Let \(\mathcal{P}=\mathcal{N}_{\rho}(L,\mathbf{n}_{L},R)\). Then for any \(p\in[d]\), \(\left|\partial_{x_{p}}^{2}f_{q}^{(\ell)}\right|\leq\ell\pi_{\ell-1}^{2}R^{2 \ell}\), \(\ell=1,2,\cdots,L\)._

Proof.: It follows from direct computation that

\[\partial_{x_{p}}^{2}f_{q}^{(\ell)}= \Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}\partial_{x_{p}}f_{j }^{(\ell-1)}\Big{)}^{2}\rho^{\prime\prime}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^ {(\ell)}f_{j}^{(\ell-1)}+b_{q}^{(\ell)}\Big{)}\] \[+\rho^{\prime}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}f_{j}^ {(\ell-1)}+b_{q}^{(\ell)}\Big{)}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}\partial _{x_{p}}^{2}f_{j}^{(\ell-1)}. \tag{5.5}\]

When \(\ell=1\), this identity and the uniform bound on \(\rho^{\prime\prime}\) from Lemma 2.1 imply

\[\|\partial_{x_{p}}^{2}f_{q}^{(1)}\|_{L^{\infty}(\Omega)}=\Big{(} \sum_{j=1}^{n_{0}}a_{qj}^{(1)}\Big{)}^{2}\|\rho^{\prime\prime}\Big{(}\sum_{j=1 }^{n_{0}}a_{qj}^{(1)}x+b_{q}^{(1)}\Big{)}\|_{L^{\infty}(\mathbb{R})}\leq n_{0} ^{2}R^{2}. \tag{5.6}\]

Next, we treat the case \(\ell>1\). The bounds on \(\rho^{\prime}\) and \(\rho^{\prime\prime}\) in Lemma 2.1 and the a priori bound on \(a_{qj}^{(\ell)}\) and the estimate (5.3) imply

\[\left|\partial_{x_{p}}^{2}f_{q}^{(\ell)}\right|\leq\Big{(}\sum_{j=1 }^{n_{\ell-1}}\left|a_{qj}^{(\ell)}\right|\Big{)}^{2}\cdot(\pi_{\ell-2}R^{\ell- 1})^{2}+R\sum_{j=1}^{n_{\ell-1}}\left|\partial_{x_{p}}^{2}f_{j}^{(\ell-1)}(x) \right|\leq\pi_{\ell-1}^{2}R^{2\ell}+R\sum_{j=1}^{n_{\ell-1}}\left|\partial_{x _{p}}^{2}f_{j}^{(\ell-1)}\right|.\]

Then applying the recursion repeatedly, the estimate (5.6) and mathematical induction yields

\[\left|\partial_{x_{p}}^{2}f_{q}^{(\ell)}\right|\leq\ell\pi_{\ell-1}^{2}R^{2 \ell},\quad\ell=1,\dots,L-1.\]

The case \(\ell=L\) follows also directly from the definition of \(f\) and the preceding estimate. 

The next result represents one of the main tools in establishing Rademacher complexity bound.

**Lemma 5.3**.: _Let \(f_{\theta},f_{\tilde{\theta}}\in\mathcal{P}\), and define \(\eta=1\) for sigmoid and \(\eta=2\) for hyperbolic tangent. Then,_

\[\left|\partial_{x_{p}}^{2}f_{\theta}(x)-\partial_{x_{p}}^{2}f_{\tilde{\theta}}(x) \right|\leq 2(L-1)L\eta\sqrt{\mathbf{n}_{L}}\pi_{L-1}^{3}R^{3L-3}\|\theta-\tilde{ \theta}\|_{\ell^{2}},\quad\forall p\in[d].\]
Proof.: The proof is based on mathematical induction, and it is divided into three steps.

(i) **Prove the bound for the case \(\ell=1\)**. By the identity (5.5) and triangle inequality, for \(\ell=1\),

\[\big{|}\partial_{x_{p}}^{2}f_{q}^{(1)}-\partial_{x_{p}}^{2}\tilde{f }_{q}^{(1)}\big{|}=\Big{|}(a_{qp}^{(1)})^{2}\rho^{\prime\prime}\Big{(}\sum_{j=1} ^{n_{0}}a_{qj}^{(1)}x_{j}+b_{q}^{(1)}\Big{)}-(\tilde{a}_{qp}^{1})^{2}\rho^{ \prime\prime}\Big{(}\sum_{j=1}^{n_{0}}\tilde{a}_{qj}^{(1)}x_{j}+\tilde{b}_{q}^ {(1)}\Big{)}\Big{|}\] \[\leq\] \[\leq 2R\big{|}a_{qp}^{(1)}-\tilde{a}_{qp}^{(1)}\big{|}+\eta R^{2}\sum _{j=1}^{n_{0}}\big{|}a_{qj}^{(1)}-\tilde{a}_{qj}^{(1)}\big{|}+\eta R^{2}\big{|} b_{q}^{(1)}-\tilde{b}_{q}^{(1)}\big{|}\leq 3\eta R^{2}\sum_{k=1}^{\mathbf{n}_{1}}\big{|} \theta_{k}-\tilde{\theta}_{k}\big{|}, \tag{5.7}\]

since by definition, \(\eta\geq 1\).

(ii) **Derive the recursion relation.** For \(\ell=2,\ldots,L-1\), in view of the identity (5.5), we have

\[\big{|}\partial_{x_{p}}^{2}f_{q}^{(\ell)}-\partial_{x_{p}}^{2} \tilde{f}_{q}^{(\ell)}\big{|}\leq \Big{|}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}\partial_{x_{p} }f_{j}^{(\ell-1)}\Big{)}^{2}-\Big{(}\sum_{j=1}^{n_{\ell-1}}\tilde{a}_{qj}^{( \ell)}\partial_{x_{p}}\tilde{f}_{j}^{(\ell-1)}\Big{)}^{2}\Big{|}\Big{|}\rho^{ \prime\prime}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}f_{j}^{(\ell-1)}+b_{ q}^{(\ell)}\Big{)}\Big{|}\] \[+\Big{|}\rho^{\prime}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell) }f_{j}^{(\ell-1)}+b_{q}^{(\ell)}\Big{)}-\rho^{\prime}\Big{(}\sum_{j=1}^{n_{ \ell-1}}\tilde{a}_{qj}^{(\ell)}\tilde{f}_{j}^{(\ell-1)}+\tilde{b}_{q}^{(\ell)} \Big{)}\Big{|}\Big{|}\Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}\partial_{x_ {p}}^{2}f_{j}^{(\ell-1)}\Big{)}\Big{|}\] \[+\Big{|}\rho^{\prime}\Big{(}\sum_{j=1}^{n_{\ell-1}}\tilde{a}_{qj} ^{(\ell)}\tilde{f}_{j}^{(\ell-1)}+\tilde{b}_{q}^{(\ell)}\Big{)}\Big{|}\Big{|} \Big{(}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}\partial_{x_{p}}^{2}f_{j}^{(\ell- 1)}-\sum_{j=1}^{n_{\ell-1}}\tilde{a}_{qj}^{(\ell)}\partial_{x_{p}}^{2}\tilde{f }_{j}^{(\ell-1)}\Big{|}=\sum_{m=1}^{4}\text{I}_{m}.\]

Below we bound the four terms separately. To bound the term \(\text{I}_{1}\), by the estimate (5.3),

\[\sum_{j=1}^{n_{\ell-1}}\Big{(}\big{|}a_{qj}^{(\ell)}\big{|}\big{|} \partial_{x_{p}}f_{j}^{(\ell-1)}\big{|}+\big{|}\tilde{a}_{qj}^{(\ell)}\big{|} \big{|}\partial_{x_{p}}\tilde{f}_{j}^{(\ell-1)}\big{|}\Big{)}\leq 2R\sum_{j=1}^{n_{\ell-1}} \pi_{\ell-2}R^{\ell-1}=2\pi_{\ell-1}R^{\ell}.\]

Likewise, by the estimates (5.3) and (5.4),

\[\sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}\partial_{x_{p}}f_{j} ^{(\ell-1)}-\tilde{a}_{qj}^{(\ell)}\partial_{x_{p}}\tilde{f}_{j}^{(\ell-1)} \big{|}\] \[\leq \sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{( \ell)}\big{|}\big{|}\partial_{x_{p}}f_{j}^{(\ell-1)}\big{|}+\sum_{j=1}^{n_{ \ell-1}}\big{|}\tilde{a}_{qj}^{(\ell)}\big{|}\big{|}\partial_{x_{p}}f_{j}^{( \ell-1)}-\partial_{x_{p}}\tilde{f}_{j}^{(\ell-1)}\big{|}\] \[\leq \pi_{\ell-2}R^{\ell-1}\sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{( \ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+\sum_{j=1}^{n_{\ell-1}}R\cdot\Big{(}\ell \pi_{\ell-2}^{2}R^{2\ell-2}\sum_{k=1}^{\mathbf{n}_{\ell-1}}\big{|}\theta_{k}- \tilde{\theta}_{k}\big{|}\Big{)}\] \[\leq \pi_{\ell-2}R^{\ell-1}\sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{( \ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+\ell\pi_{\ell-2}\pi_{\ell-1}R^{2\ell-1} \sum_{j=1}^{\mathbf{n}_{\ell-1}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}.\]

Thus we can bound the term \(\text{I}_{1}\) by

\[\text{I}_{1}\leq 2\pi_{\ell-2}\pi_{\ell-1}R^{2\ell-1}\sum_{j=1}^{n_{\ell-1}} \big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+2\ell\pi_{\ell-2}\pi_{ \ell-1}^{2}R^{3\ell-1}\sum_{j=1}^{\mathbf{n}_{\ell-1}}\big{|}\theta_{j}-\tilde{ \theta}_{j}\big{|}.\]

Next, it follows from the estimate (5.2) that

\[\Big{|}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{(\ell)}f_{j}^{(\ell-1)}-\tilde{a}_{qj}^{( \ell)}\tilde{f}_{j}^{(\ell-1)}+b_{q}^{(\ell)}-\tilde{b}_{q}^{(\ell)}\Big{|}\leq \sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+ \big{|}b_{q}^{(\ell)}-\tilde{b}_{q}^{(\ell)}\big{|}+\pi_{\ell-1}R^{\ell-1}\sum_ {j=1}^{\mathbf{n}_{\ell-1}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}.\]
Then for the term \(\mathrm{I}_{2}\), from the estimate (5.3), we deduce

\[\mathrm{I}_{2} \leq\eta\Big{(}\sum_{j=1}^{n_{\ell-1}}\big{|}\tilde{a}_{qj}^{(\ell )}\big{|}^{2}\Big{)}\Big{(}\sum_{j=1}^{n_{\ell-1}}\big{|}\partial_{x_{p}}\tilde{ f}_{j}^{(\ell-1)}\big{|}^{2}\Big{)}\Big{(}\Big{|}\sum_{j=1}^{n_{\ell-1}}a_{qj}^{( \ell)}f_{j}^{(\ell-1)}-\tilde{a}_{qj}^{(\ell)}\tilde{f}_{j}^{(\ell-1)}+b_{q}^{ (\ell)}-\tilde{b}_{q}^{(\ell)}\Big{|}\Big{)}\] \[\leq\eta\pi_{\ell-1}^{2}R^{2\ell}\Big{(}\sum_{j=1}^{n_{\ell-1}} \big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+\big{|}b_{q}^{(\ell)}- \tilde{b}_{q}^{(\ell)}\big{|}+\pi_{\ell-1}R^{\ell-1}\sum_{j=1}^{\boldsymbol{n }_{\ell-1}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}\Big{)}.\]

Similarly, Lemma 5.2 implies

\[\sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}\partial_{x_{p}}^{2}f_{j}^{(\ell- 1)}\Big{|}\leq R\cdot n_{\ell-1}\cdot(\ell-1)\pi_{\ell-2}^{2}R^{2(\ell-1)}=( \ell-1)\pi_{\ell-2}\pi_{\ell-1}R^{2\ell-1},\]

and thus we can bound the term \(\mathrm{I}_{3}\) by

\[\mathrm{I}_{3}\leq(\ell-1)\pi_{\ell-2}\pi_{\ell-1}R^{2\ell-1}\Big{(}\sum_{j=1 }^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+\big{|}b_{ q}^{(\ell)}-\tilde{b}_{q}^{(\ell)}\big{|}+\pi_{\ell-1}R^{\ell}\sum_{j=1}^{ \boldsymbol{n}_{\ell-1}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}\Big{)}.\]

For the last term \(\mathrm{I}_{4}\), using Lemma 5.2 again, we have

\[\mathrm{I}_{4} \leq\sum_{j=1}^{n_{\ell-1}}\big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^ {(\ell)}\big{|}\big{|}\partial_{x_{p}}^{2}f_{j}^{(\ell-1)}\big{|}+R\sum_{j=1}^ {n_{\ell-1}}\big{|}\partial_{\partial x_{p}}^{2}f_{j}^{(\ell-1)}-\partial_{x_{ p}}^{2}\tilde{f}_{j}^{(\ell-1)}\big{|}\] \[\leq(\ell-1)\pi_{\ell-2}\pi_{\ell-1}R^{2(\ell-1)}\sum_{j=1}^{n_{ \ell-1}}\big{|}a_{qj}^{(\ell)}-\tilde{a}_{qj}^{(\ell)}\big{|}+R\sum_{j=1}^{n_{ \ell-1}}\big{|}\partial_{x_{p}}^{2}f_{j}^{(\ell-1)}-\partial_{x_{p}}^{2}\tilde {f}_{j}^{(\ell-1)}\big{|}.\]

Combining the last four estimates gives the crucial recursion

\[\big{|}\partial_{x_{p}}^{2}f_{q}^{(\ell)}-\partial_{x_{p}}^{2}\tilde{f}_{q}^{( \ell)}\big{|}\leq R\sum_{j=1}^{n_{\ell-1}}\big{|}\partial_{x_{p}}^{2}f_{j}^{( \ell-1)}-\partial_{x_{p}}^{2}\tilde{f}_{j}^{(\ell-1)}\big{|}+4\ell\eta\pi_{ \ell-1}^{3}R^{3\ell-1}\sum_{j=1}^{\boldsymbol{n}_{\ell}}\big{|}\theta_{j}- \tilde{\theta}_{j}\big{|}. \tag{5.8}\]

(iii) **Prove the intermediate case by mathematical induction.** Using the recursion (5.8), we claim that for \(\ell=1,2,\ldots,L-1\), there holds

\[\big{|}\partial_{x_{p}}^{2}f_{q}^{(\ell)}(x)-\partial_{x_{p}}^{2}\tilde{f}_{q}^ {(\ell)}(x)\big{|}\leq 2\ell(\ell+1)\eta\pi_{\ell-1}^{3}R^{3\ell-1}\sum_{j=1}^{ \boldsymbol{n}_{\ell}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}. \tag{5.9}\]

This is trivially true for \(\ell=1\), by the estimate (5.7). Indeed, for \(\ell=2\), the recursion (5.8) and the estimate for \(\ell=1\) in (5.7) imply

\[\big{|}\partial_{x_{p}}^{2}f_{q}^{(2)}(x)-\partial_{x_{p}}^{2} \tilde{f}_{q}^{(2)}(x)\big{|} \leq 3\eta n_{1}R^{3}\sum_{j=1}^{\boldsymbol{n}_{1}}\big{|} \theta_{k}-\tilde{\theta}_{k}\big{|}+8\eta R^{5}\pi_{1}^{3}\sum_{j=1}^{ \boldsymbol{n}_{2}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}\] \[\leq 11\eta\pi_{1}^{3}R^{5}\sum_{j=1}^{\boldsymbol{n}_{2}}\big{|} \theta_{j}-\tilde{\theta}_{j}\big{|}\leq 2\cdot 2\cdot 3\eta\pi_{1}^{3}R^{5}\sum_{j=1}^{ \boldsymbol{n}_{2}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|},\]

and hence the claim holds for \(\ell=2\). Now suppose it holds for some \(2\leq\ell<L-1\). Then for the case \(\ell+1\),

\[\Big{|}\partial_{x_{p}}^{2}f_{q}^{(\ell+1)}(x)-\partial_{x_{p}}^{2 }\tilde{f}_{q}^{(\ell+1)}(x)\Big{|} \leq R\sum_{j=1}^{n_{\ell}}\Big{|}\partial_{x_{p}}^{2}f_{j}^{( \ell)}-\partial_{x_{p}}^{2}\tilde{f}_{j}^{(\ell)}\Big{|}+4(\ell+1)\eta\pi_{\ell}^ {3}R^{3\ell+2}\sum_{j=1}^{\boldsymbol{n}_{\ell+1}}\Big{|}\theta_{j}-\tilde{ \theta}_{j}\Big{|}\] \[\leq n_{\ell}\cdot 2\ell(\ell+1)\eta\pi_{\ell-1}^{3}R^{3\ell-1}\sum_{j=1 }^{\boldsymbol{n}_{\ell}}\big{|}\theta_{j}-\tilde{\theta}_{j}\big{|}+4(\ell+1) \eta\pi_{\ell}^{3}R^{3\ell+2}\sum_{j=1}^{\boldsymbol{n}_{\ell+1}}\big{|}\theta_{j} -\tilde{\theta}_{j}\big{|}\]
\[\leq 2(\ell+1)(\ell+2)\pi_{\ell}^{3}R^{3\ell+2}\sum_{j=1}^{\mathbf{n}_{\ell+1} }|\theta_{j}-\tilde{\theta}_{j}|,\]

which by mathematical induction implies the desired claim (5.9).

(iv) **Obtained the final estimate**. Last, for \(\ell=L\), direct computation shows

\[\partial_{x_{p}}^{2}f-\partial_{x_{p}}^{2}\tilde{f}=\sum_{j=1}^{n_{L-1}}a_{1j} ^{(L)}\partial_{x_{p}}^{2}f_{j}^{(L-1)}-\sum_{j=1}^{n_{L-1}}\tilde{a}_{1j}^{(L )}\partial_{x_{p}}^{2}\tilde{f}_{j}^{(L-1)}\]

Thus, by the estimate (5.9) and Lemma 5.2,

\[\left|\partial_{x_{p}}^{2}f-\partial_{x_{p}}^{2}\tilde{f}\right| \leq\sum_{j=1}^{n_{L-1}}|a_{1j}^{(L)}-\tilde{a}_{1j}^{(L)}|| \partial_{x_{p}}^{2}f_{j}^{(L-1)}|+\sum_{j=1}^{n_{L-1}}|\tilde{a}_{1j}^{(L)}| \partial_{x_{p}}^{2}f_{j}^{(L-1)}-\partial_{x_{p}}^{2}\tilde{f}_{j}^{(L-1)}|\] \[\leq(L-1)\pi_{L-2}R^{2L-2}\sum_{j=1}^{n_{L-1}}|a_{1j}^{(L)}- \tilde{a}_{1j}^{(L)}|+\sum_{j=1}^{n_{L-1}}R\cdot 2(L-1)L\eta\pi_{L-2}^{3}R^{3L-4} \sum_{j=1}^{\mathbf{n}_{L-1}}|\theta_{j}-\tilde{\theta}_{j}|\] \[\leq 2(L-1)L\eta\pi_{L-1}^{3}R^{3L-3}\sum_{j=1}^{\mathbf{n}_{L}} \left|\theta_{j}-\tilde{\theta}_{j}\right|.\]

This completes the proof of the lemma. 

## 6 Numerical experiments and discussions

Now we present numerical experiments to illustrate C-PINN, and evaluate its performance against existing techniques, including penalty method (PM), augmented Lagrangian method (ALM), and adjoint oriented neural network (AONN) described in Section 3.4. Throughout, we fix the weight \(\alpha_{i}=\lambda^{-1}\) in C-PINN, and use the same boundary weights for state and adjoint PDEs, i.e. \(\alpha_{b}^{y}=\alpha_{b}^{p}/\alpha_{i}\). Thus we only need to tune one single hyper-parameter \(\alpha_{b}^{y}\) in C-PINN. ANOO PM and ALM involve more hyper-parameters. The boundary weight \(\alpha_{b}\) of the PINN residual is determined by grid search [32]: we apply PINN to solve a direct problem with known ground truth, with a range of values of \(\alpha_{b}\), and then select the value \(\alpha_{b}^{*}\) attaining the smallest validation error after a fixed number of L-BFGS iterations as the optimal one, which is then used in the solver for the state and adjoint problems. Other hyper-parameters, e.g., initial penalty factor and increasing factor \(\beta\) in PM and ALM, and step size and the number \(K\) of gradient descent steps in AONN are all determined manually; see Table 1 for the hyper-parameters for the examples. For AONN, the tuple \(N_{\text{aonn}}\) denotes the iteration number for PDE solvers (both PDE and adjoint PDE) and for the gradient descent step in (3.12); and for PM and ALM, the tuple \(N_{\text{pm}}\) and \(N_{\text{alm}}\) denote the iteration for the first and remaining sub-problems. The multipliers in ALM are initialized to zero. The function \(\rho\) is taken to be \(\tanh\). The DNN parameters are initialized by the Pytorch default Xavier scheme (and zero bias). To measure the accuracy of the approximate state \(y^{*}\) and control \(u^{*}\), we employ the relative \(L^{p}(\Omega)\) error, \(p\in\{2,\infty\}\), defined by \(e_{p}(y^{*})=\|y^{*}-\bar{y}\|_{L^{p}(\Omega)}/\|\bar{y}\|_{L^{p}(\Omega)}\), where \(\bar{y}\) denotes the exact state, and similarly the relative error \(e_{p}(u^{*})\) for the control \(u^{*}\). Examples 6.1, 6.2 and 6.4 are conducted on Macbook air with Apple M2 silicon, 16 GB RAM and mac OS Monterey using double precision. Example 6.3 is conducted on an NVIDIA RTX A6000 GPU in single precision. The Python code for reproducing the results will be made publicly available.

First, we consider the linear elliptic control problem:

\[\min_{y,u}\Big{\{}J(y,u)=\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{ \lambda}{2}\|u\|_{L^{2}(\Omega)}^{2}\Big{\}},\quad\text{subject to}\ \begin{cases}-\Delta y=f+u,&\text{in}\ \Omega,\\ y=g,&\text{on}\ \partial\Omega,\end{cases}\]

where \(f\) and \(g\) are the known problem data and \(\lambda\) is the penalty weight.

The first example is about an unconstrained elliptic optimal control problem.

**Example 6.1**.: _The domain \(\Omega=\{(r,\theta):1\leq r\leq 3,\theta\in[0,2\pi]\}\), and the penalty weight \(\lambda=0.01\), \(y_{d}(r,\theta)=r^{2}+3\lambda(1-\frac{1}{r^{2}})\sin(\theta)\), \(f(r,\theta)=(r-1)(r-3)\sin(\theta)-4\), and Dirichlet data \(g(r,\theta)=r^{2}\). The exact state \(\bar{y}\) and control \(\bar{u}\) are given by \(\bar{y}(r,\theta)=r^{2}\) and \(\bar{u}(r,\theta)=(r-1)(r-3)\sin(\theta)\)._
For both state and adjoint state, we use a DNN architecture with 4 hidden layers and 30 neurons in each layer, which is also used to approximate the control \(u\) in the benchmark methods. When formulating the empirical loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\), we take \(n_{d}=10000\) points i.i.d. from \(U(\Omega)\) for the PDE residual, and \(n_{b}=3000\) points i.i.d. from \(U(\partial\Omega)\). For all four methods, we minimize the loss \(\widehat{\mathcal{L}}(y_{\theta},p_{\sigma})\) using L-BFGS provided by the PyTorch library torch.optim (version 2.0.0). The relative errors of the approximate state \(y^{*}\) and control \(u^{*}\) are presented in Table 2. The computing time is also given, which is determined by choosing a near-minimum number of iterations that still achieves the reported accuracy. All the methods can achieve a relative error about \(10^{-4}\). In view of the well-known accuracy saturation of PINN-type solvers for direct problems [34, 25], the attained accuracy for the optimal control problem based on PINN is expected to be nearly the best possible, at least empirically. In Fig. 2, with a fixed step size \(s=10\), AONN requires about 20 outer iterations to reach convergence. Penalty-based methods can capture the state accurately, but the control approximations are less accurate. In contrast, C-PINN gives the most accurate state approximation, lowest objective value within a shorter computing time.

In Fig. 2 we show the training dynamics of the four methods. It is observed for C-PINN, the optimizer converges rapidly in terms of the objective value, and the PINN loss for the state and adjoint. In contrast, for AONN, it does minimize the objective steadily, but the PINN losses \(\widehat{\mathcal{L}}_{\mathrm{pinn}}(y)\) and \(\widehat{\mathcal{L}}_{\mathrm{pinn}}(p)\) exhibit alternating convergence behavior due to the alternating nature of the algorithm. PM and ALM show very similar convergence behavior in the early stage, and the objective value enjoys fast initial decay but the convergence slows down greatly as the iteration proceeds, and also the PINN loss \(\widehat{\mathcal{L}}_{\mathrm{pinn}}(y)\) also decreases

\begin{table}
\begin{tabular}{c|c c c c|c c c c} \hline \hline Example & \(\alpha_{b}\) & \((s,K)_{\mathrm{aonn}}\) & \((\mu_{0},\beta,K)_{\mathrm{pm}}\) & \((\mu,K)_{\mathrm{alm}}\) & \(N_{\mathrm{cpinn}}\) & \(N_{\mathrm{p}m}\) & \(N_{\mathrm{alm}}\) \\ \hline
6.1 & 5 & (10, 30) & (0.1, 2, 8) & (0.1, 8) & 1.5e4 & (1e3, 5e2) & (6e3, 3e3) & (6e3, 3e3) \\ \hline
6.2 & 5 & (10, 30) & (0.1, 2, 8) & (0.1, 8) & 1.5e4 & (1e3, 5e2) & (1.2e4, 1e4) & (6e3, 3e3) \\ \hline
6.3 & 100 & (100, 15) & (0.5, 2, 5) & (2, 8) & 6e4 & (6e3, 1.5e3) & (2e4, 2e4) & (2e4, 2e4) \\ \hline
6.4 & 100 & (100, 20) & (1, 2, 8) & (16, 8) & 2e4 & (2e3, 5e2) & (1e4, 8e3) & (1e4, 8e3) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Hyper-parameter settings for the Examples. \(\alpha_{b}\) denotes boundary penalty weight, \(s\) denotes the step size in AONN, \(\mu_{0}\) the initial penalty weight and \(\beta\) the increasing factor, \(\mu\) penalty weight in ALM, and \(K\) denotes the update / path-following steps, and \(N\) the maximum number of iterations for the optimizer.

Figure 1: The approximate control \(u^{*}\) (top), its pointwise error \(|u^{*}-\bar{u}|\) (middle) and approximate state \(y^{*}\) (bottom) obtained by four DNN-based methods for Example 6.1.
rapidly throughout. There is a drastic increase of the PINN error in the later stage of ALM, indicating training instability. This can be alleviated by choosing a larger penalty weight, as usually recommended for convex optimization [5, 37, 6], or safeguarding the multipliers [6]. We choose \(\mu=0.1\) in the experiment since it leads to the smallest cost objective with an acceptable PINN loss in the grid search process.

All methods involve some hyper-parameters, and the ease of their tuning is of great importance. Table 3 shows the variation of the performance metric with respect to one hyper-parameter. C-PINN enjoys good accuracy for a wide range of parameters. AONN is effective if the step size \(s\) is carefully tuned: a too large \(s\) can lead to divergence, and a too small \(s\) requires many outer iterations to yield accurate approximations. PM and ALM also stably yield accurate results, but in PM, a larger \(\mu_{0}\) can drastically increase the error, and more iterations are required to make the sub-problem converge, due to their ill-condition.

Next, we consider a constrained problem.

**Example 6.2**.: _The domain \(\Omega\), penalty weight \(\gamma\), and problem data \(y_{d}\) and \(g\) are identical with that for Example 6.1, but with box constraint on the control, i.e., \(u\in U:=\{u\in L^{2}(\Omega):u(x)\in[-0.5,0.7]\) a.e. in \(\Omega\}\). The source \(f(r,\theta)=-4-P_{[-0.5,0.7]}((1-r)(r-3)\sin(\theta))\), and the exact state \(\bar{y}\) and control \(\bar{u}\) are given by \(\bar{y}(r,\theta)=r^{2}\) and \(\bar{u}(r,\theta)=P_{[-0.5,0.7]}((1-r)(r-3)\sin(\theta))\), respectively._

In addition to the hyper-parameters in Example 6.1, in PM and ALM, we impose the constraint by including a penalty term with a weight \(\mu^{\prime}=1\) and increasing factor \(\beta=2\), cf. (3.13). Meanwhile, since the

\begin{table}

\end{table}
Table 2: The numerical results for Example 6.1.

\begin{table}

\end{table}
Table 3: Hyper-parameter sensitivity of the four methods for Example 6.1.

Figure 2: The training dynamics of the four methods, for Example 6.1.
loss is less smooth due to the presence of the projection operator \(P_{U}\), we increase the maximum iteration numbers for PM. The results are reported in Table 4 and Fig. 3. Since the optimal control \(\bar{u}\) is non-smooth near the boundary of the constraint active area (and more challenging to approximate with DNNs), all methods except C-PINN suffer from pronounced errors therein. In C-PINN, this error is greatly alleviated due to the use of \(P_{U}\) a posteriori on the adjoint \(p^{*}\). Both PM and ALM become less efficient with large relative \(L^{2}(\Omega)\) errors in the approximation \(u^{*}\), and numerically, the choice of the penalty parameter \(\mu\) becomes even more crucial for obtaining acceptable DNN approximations.

The next example is about the medium-dimensional case.

**Example 6.3**.: _The domain \(\Omega=(0,1)^{4}\), and penalty weight \(\lambda=0.01\). The data \(f\equiv 0\), \(y_{d}(x)=(1-16\lambda\pi^{4})\prod_{i=1}^{4}\sin(\pi x_{i})\), and \(g=0\). The exact state \(\bar{y}\) and control \(\bar{u}\) are given by \(\bar{y}(x)=\prod_{i=1}^{4}\sin(\pi x_{i})\), and \(\bar{u}(x)=4\pi^{2}\prod_{i=1}^{4}\sin(\pi x_{i})\), respectively._

The experiment setup is similar to Example 6.1, but with the following changes. The DNN has 4 hidden layers, each having 80 neurons. To form the empirical loss, we employ 60000 points in the domain \(\Omega\) and 5000 points on the boundary \(\partial\Omega\), drawn i.i.d. from \(U(\Omega)\) and \(U(\partial\Omega)\), respectively. Due to the larger data size, we conduct all the training in single precision. See Table 1 for the hyper-parameters. Note that for PM, we terminate the iteration at a relatively small weight \(\mu_{4}=8\) because the sub-problem hardly converges within the specified optimizer iterations when the weight continues to grow. Due to the increased complexity of the problem, we employ 60000 Adam iterations for C-PINN with a learning rate \(10^{-3}\), divide it by 10 after the 20000th and 40000th iterations; and 6000 Adam iterations for each PDE solve in AONN with

\begin{table}
\begin{tabular}{c|c|c|c|c|c|c} \hline  & \(e_{2}(y^{*})\) & \(e_{\infty}(y^{*})\) & \(e_{2}(u^{*})\) & \(e_{\infty}(u^{*})\) & \(J(u^{*},y^{*})\) & time (s) \\ \hline C-PINN & 3.49e-5 & 7.91e-5 & 6.16e-3 & 2.39e-2 & 1.026e-3 & 1848.68 \\ AONN & 6.86e-5 & 1.91e-4 & 3.53e-2 & 1.65e-1 & 1.029e-3 & 4676.55 \\ PM & 2.85e-4 & 3.98e-4 & 6.23e-2 & 2.72e-1 & 1.032e-3 & 6492.62 \\ ALM & 3.00e-4 & 7.89e-4 & 6.87e-2 & 2.94e-1 & 1.040e-3 & 2543.71 \\ \hline \end{tabular}
\end{table}
Table 4: The numerical results for Example 6.2.

Figure 3: The approximate control \(u^{*}\) (top), its pointwise error \(|u^{*}-\bar{u}|\) (mid) and approximate state \(y^{*}\) (bottom) obtained by four DNN-based methods for Example 6.2.
15 gradient descent steps in total. In PM and ALM, we employ 20000 iterations for each sub-problem. In AONN, PM, and ALM, we take a learning rate 1e-3 in the first sub-problem, and then 1e-4 for the remaining sub-problems. The numerical results are reported in Table 5. C-PINN and AONN achieve the relative error at about \(10^{-2}\) for both \(y^{*}\) and \(u^{*}\). C-PINN yields a more accurate control approximation \(u^{*}\) in the \(L^{2}(\Omega)\) error, and the achieved objective value is also lower than that by AONN. PM and ALM both report good accuracy for \(y^{*}\), but much larger errors for \(u^{*}\).

The last example is about optimal control with semilinear elliptic PDEs.

**Example 6.4**.: _Consider the following semilinear PDE:_

\[\begin{cases}-\Delta y+y+q(x,y)=f+u,&\text{in }\Omega,\\ y=g,&\text{on }\partial\Omega,\end{cases}\]

_with \(q(x,y)=k(x)y(x)^{3}\). The domain \(\Omega=(0,1)^{2}\) is divided into two parts \(\omega=[0.25,0.75]^{2}\) and \(\Omega\setminus\omega\), and the function \(k\equiv 1\) in \(\omega\), and \(k\equiv 3\) in \(\Omega\setminus\omega\). We take a penalty weight \(\lambda=0.01\). The problem data \(y_{d}\), \(g\) and \(f\) are taken such that the exact state \(\bar{y}\) and control \(\bar{u}\) are given by \(\bar{y}(x)=e^{x_{1}(1-x_{1})}\sin(\pi x_{2})+e^{x_{2}(1-x_{2})}\sin(\pi x_{1})\), and \(\bar{u}(x)=-\lambda^{-1}x_{1}x_{2}(1+\cos(\pi x_{1})(1+\cos(\pi x_{2})))\), respectively._

We employ the DNN architecture for Example 6.1. Due to the nonlinearity of the state problem, each forward solve requires more L-BFGS iterations to reach convergence, cf. Table 1 for the precise setting. The results are reported in Table 6 and Fig. 5. Since exact state \(\bar{y}\) and control \(\bar{u}\) have very good regularity,

\begin{table}
\begin{tabular}{c|c|c|c|c|c|c} \hline \hline  & \(e_{2}(y^{*})\) & \(e_{\infty}(y^{*})\) & \(e_{2}(u^{*})\) & \(e_{\infty}(u^{*})\) & \(J(u^{*},y^{*})\) & time (sec) \\ \hline C-PINN & 2.44e-2 & 8.25e-2 & 9.29e-3 & 4.37e-2 & 8.020 & 6738.26 \\ AONN & 6.83e-3 & 1.34e-2 & 1.03e-2 & 2.38e-2 & 8.031 & 9579.84 \\ PM & 3.27e-2 & 2.10e-2 & 9.49e-2 & 1.05e-1 & 8.024 & 7342.11 \\ ALM & 1.99e-2 & 1.87e-2 & 7.25e-2 & 8.83e-2 & 8.033 & 12485.55 \\ \hline \hline \end{tabular}
\end{table}
Table 5: The numerical results for Example 6.3.

Figure 4: The approximate optimal control \(u^{*}\) (top), its pointwise error \(|u^{*}-\bar{u}|\) (mid), and approximate state \(y^{*}\) (bottom) obtained by four DNN-based methods for Example 6.3, cross section at \(x_{3}=x_{4}=0.5\).
both PM and ALM can yield accurate approximations, and yield smaller cost objectives at the expense of less accuracy in the state solution. In contrast, both C-PINN and AONN can approximate the state and control with smaller \(L^{2}(\Omega)\) and \(L^{\infty}(\Omega)\) errors, but C-PINN approximates the optimal control more accurately within a shorter time. Note that the efficiency of AONN depends heavily on the step size \(s\), but line search is not effective in the context since it involves multiple expensive forward solves. The step size \(s=100\) was determined in a trial-and-error manner in order to achieve good efficiency.

## 7 Conclusion

In this work, we have developed a new approach to solve distributed optimal control problems associated with elliptic PDEs, with / without box constraint on the control variable. The approach is derived from the first-order necessary optimality system, and then applies physics informed neural networks to the reduced optimality system. We provide an error analysis of the approach, using the approximation theory of neural networks and offset Rademacher complexity. The \(L^{2}(\Omega)\) error bound on the control, state, and adjoint is explicit in terms of the DNN parameters (depth, width), and the number of Monte Carlo sampling points in the domain and on the boundary. Several examples confirm the feasibility of the approach for linear and semilinear elliptic problems, indicating its potential for solving a class of optimal control problems.

\begin{table}
\begin{tabular}{c|c|c|c|c|c} \hline \hline  & \(e_{2}(y^{*})\) & \(e_{\infty}(y^{*})\) & \(e_{2}(u^{*})\) & \(e_{\infty}(u^{*})\) & \(J(u^{*},y^{*})\) & time (sec) \\ \hline C-PINN & 6.49e-5 & 3.36e-4 & 2.68e-4 & 1.66e-3 & 16.063 & 2584.97 \\ AONN & 4.70e-5 & 3.48e-4 & 1.17e-3 & 4.37e-3 & 16.063 & 3335.25 \\ PM & 1.20e-3 & 4.04e-3 & 4.16e-2 & 9.32e-2 & 16.060 & 5262,92 \\ ALM & 9.60e-4 & 2.01e-3 & 3.09e-2 & 7.93e-2 & 16.062 & 5077.21 \\ \hline \hline \end{tabular}
\end{table}
Table 6: The numerical results for Example 6.4.

Figure 5: The approximate optimal control \(u^{*}\) (top), its pointwise error \(|u^{*}-\bar{u}|\) (middle), and the state \(y^{*}\) (bottom) obtained by four different methods for Example 6.4.