¹

1 Introduction↩︎

The control theory of nonlinear partial differential equations (PDEs) has been the subject of sustained and intensive investigation over the past several decades; we refer, for instance, to [1]–[3] and the extensive bibliography therein. Despite this substantial body of work, the literature remains comparatively limited on the numerical reconstruction of controls for nonlinear systems. In particular, while many controllability results, including local, global, approximate, and null controllability, have been established under diverse structural assumptions, considerably less is known about how to compute reliable and accurate numerical approximations of the corresponding control functions.

This difficulty is not merely technical. The intrinsic nonlinearity of the dynamics typically precludes the use of superposition principles and may lead to instabilities, loss of compactness, or sensitivity with respect to perturbations. These features complicate both the analytical and numerical reconstruction of controls. Consequently, the design of stable, consistent, and convergent numerical methods for control reconstruction remains a challenging problem, especially in genuinely nonlinear regimes.

This question appears to have been investigated only relatively recently, starting with the initial contribution [4], where a theoretical proof of a numerical scheme was provided without numerical results. Since then, it has begun to develop into an emerging research direction, although it remains far from fully established and is still comparatively unexplored. Constructive and least-squares approaches have been proposed for semilinear wave and heat equations, including distributed and boundary control problems, as well as space–time methods for the approximation of controls; see [5]–[9]. These contributions highlight both the difficulty and the importance of obtaining accurate numerical approximations of controls in nonlinear PDE settings. It is also worth noting that classical numerical methods for PDEs may face increasing computational challenges in high-dimensional problems, a phenomenon commonly referred to as the curse of dimensionality. This provides further motivation for developing mesh-free or sampling-based methods for high-dimensional PDE control problems.

Physics-informed neural networks (PINNs) are neural-network-based methods for solving differential equations by incorporating the governing equations directly into the training loss. Introduced in [10]–[12], PINNs approximate the solution by a neural network trained with data, boundary conditions, and PDE residuals computed by automatic differentiation. They have become a central tool in scientific machine learning, with applications to forward and inverse problems, particularly in data-sparse settings [12]. From a computational perspective, substantial effort has been devoted to improving the robustness and accuracy of PINNs, especially in connection with optimization difficulties, stiffness, and loss imbalance [13]. Important developments include conservative formulations, domain decomposition methods, adaptive activation functions, adaptive weighting strategies, and self-adaptive PINNs; see, for instance, [14]–[16]. On the theoretical side, convergence and generalization properties have begun to be understood for certain classes of PDEs [17], [18]. We also refer to the recent surveys [19], [20] for further developments and perspectives.

1.1 Positioning with respect to existing PINN methodologies↩︎

The PINN literature has grown rapidly since [10]–[12], and many methodological refinements have been proposed. Since the present work introduces an adaptive PINN framework designed specifically for the numerical reconstruction of controls for nonlinear PDEs, it is useful to position our contribution within the broader landscape of PINN variants. We organize the discussion into eight categories. In each case, we explain both the principal mechanism of the existing method and the aspects of the control reconstruction problem that remain unaddressed.

1.1.0.1 (A) Loss-balancing methods at the level of whole loss terms.

A standard PINN loss is typically written as a weighted sum \[\mathcal{L}_{\mathrm{PINN}} = \mathcal{L}_{\mathrm{PDE}} + \lambda_{\mathrm{bd}}\mathcal{L}_{\mathrm{bd}} + \lambda_{\mathrm{ic}}\mathcal{L}_{\mathrm{ic}} +\cdots .\] The different components of this loss may have very different magnitudes and gradient scales during training, leading to severe imbalance. This issue is analyzed through gradient-flow pathologies in [13] and from a neural tangent kernel perspective in [21]. This line of work also proposes dynamic update rules for the scalar coefficients \(\lambda_\bullet\). The inverse Dirichlet weighting method of [22] provides an alternative rebalancing rule based on residual variance. These methods improve the scaling among different loss components, but they operate at the level of whole loss terms. They do not distinguish the individual operator components inside the PDE residual itself, such as the time derivative, the diffusion operator, the nonlinear response, or the control term. Nor do they separately treat the control constraint that defines the controllability problem.

1.1.0.2 (B) Self-adaptive PINNs with pointwise weights.

A second family attaches trainable per-collocation-point weights to the residual, replacing \[\sum_j |R[u_\theta](t_j,x_j)|^2\] by \[\sum_j \lambda_j^2 |R[u_\theta](t_j,x_j)|^2,\] where the weights \(\lambda_j\) are updated to emphasize ill-resolved points [16]. This creates a soft attention mechanism over the collocation points and improves local accuracy where the residual is large. Such adaptivity, however, is attached to points or samples rather than to the differential operator. It does not identify which component of the PDE residual, for example the time derivative, diffusion term, nonlinear term, or control term, is responsible for the local error. It is also not specifically designed for reconstructing controls that enter the equation through a particular structural mechanism.

1.1.0.3 (C) Residual-based adaptive sampling.

A third family modifies neither the residual norm nor the loss weights, but instead changes the distribution of collocation points. Residual-adaptive refinement (RAR) adds collocation points where the residual is large; the family of RAR, RAR-D, and RAD strategies has been systematized in [23] and is implemented in the DeepXDE framework [24]. Related ideas include retain-and-resample (\(R^3\)) sampling [25]. These methods improve the resolution of localized features by changing where the residual is evaluated. However, the residual itself remains a single scalar quantity with no operator-level decomposition, and the method is agnostic to the mechanism through which the control enters the PDE.

1.1.0.4 (D) Weighted residual and variational PINNs.

A fourth family replaces \(R[u_\theta]\) by \(\phi(t,x)R[u_\theta]\), where \(\phi\) is a suitable test or weight function; equivalently, the residual is measured in a weighted norm. Variational PINNs (VPINNs) and their \(hp\)-adaptive versions (hp-VPINNs) test the residual against piecewise polynomial spaces [26], [27]. Gradient-enhanced PINNs (gPINNs) [28] add gradients of the residual to the loss. In these approaches, the weight, test function, or derivative term acts on the residual as a whole. It does not separately weight the components \[\partial_t u_\theta,\qquad \partial_{x_i x_i}u_\theta,\qquad v(u_\theta),\qquad f_{\bar\theta}\mathbf{1}_\omega, \qquad\text{or}\qquad f_{\bar\theta}u_\theta.\] This distinction is essential in the control setting, because the unknown control is identified through the balance among these individual components.

1.1.0.5 (E) Causal and curriculum-style PINNs.

For time-dependent PDEs, several works improve the temporal structure of training. Causal PINNs reweight the residual at later times so that earlier times must be resolved first [29], while time-marching and sequence-to-sequence formulations train the network on growing time windows. These methods address temporal causality and improve training stability in evolutionary problems. However, they do not distinguish among operator components at a fixed point in space-time, and they have not been formulated specifically for the reconstruction of controls.

1.1.0.6 (F) Domain decomposition methods.

Conservative PINNs (cPINNs) [14], extended PINNs (XPINNs) [30], finite-basis PINNs (FBPINNs) [31], and parallel PINN variants split the space-time domain into subdomains, each with its own neural representation. These methods are powerful for multiscale problems, but they primarily address the geometry and decomposition of the training domain rather than the structure of the PDE residual. They do not provide an operator-decomposed adaptive mechanism for control reconstruction.

1.1.0.7 (G) Architecture-level adaptivity and hard constraints.

Locally adaptive activation functions (LAAF) [15] and adaptive activation slopes modify the neural representation itself. Hard-constrained PINNs (hPINNs) [32] embed constraints directly into the network architecture, while augmented-Lagrangian PINNs (AL-PINNs) [33] treat boundary and initial conditions as constraints rather than penalty terms. Physics-constrained neural networks (PCNNs), exact boundary lifting, and related formulations share a similar philosophy. These methods are effective for enforcing constraints, but their primary objective is not the numerical reconstruction of controls, and they do not address the operator-component imbalance that arises when an unknown control is identified through a multiplicative or localized term.

1.1.0.8 (H) Other specialized variants.

Bayesian PINNs (B-PINNs) [34] introduce uncertainty quantification; fractional PINNs (fPINNs) [35] handle fractional derivatives; and PINNs combined with Kolmogorov–Arnold networks (PIKANs) [36] explore alternative neural representations. These contributions are important and complementary to ours, but they do not directly address the control reconstruction problem considered in this paper.

To the best of our knowledge, none of the adaptive PINN methodologies described above has been systematically developed for the numerical reconstruction of controls for nonlinear PDEs, with internal or bilinear controls, treated here.

1.2 The control reconstruction setting and the WeightedPINN formulation↩︎

The control setting presents a different and more delicate difficulty from that of a forward PDE problem. In a forward problem, the primary objective is to approximate the state \(u\). In a control reconstruction problem, by contrast, both the state and the control are unknown, and the control is identified only indirectly through the state equation, the boundary and initial conditions, and the terminal target. For example, in the internal control setting one has \[\partial_t u-\Delta u+v(u)=f\mathbf{1}_\omega,\] whereas in the bilinear control setting one has \[\partial_t u-\Delta u+v(u)=fu.\] Thus, the control enters through a specific structural mechanism. In the bilinear case, the control appears through the product \(fu\), so errors in the state and errors in the control interact multiplicatively. Consequently, a small averaged PDE residual does not necessarily imply that the recovered control is accurate. Likewise, an improved approximation of the state alone may still leave the control poorly reconstructed.

Although PINN-based methodologies have recently been applied to optimal deterministic control problems [37]–[40], as well as to stochastic control settings [41], [42], our objective here is different. We develop a new PINN-based approach, which we call WeightedPINNs; see Subsection 2.3 for an illustrative explanation. The method is tailored to the computation of reliable and accurate numerical approximations for the control reconstruction problem initiated in [4].

The central idea is to introduce adaptive weights directly into the PINN formulation at the level of the controlled PDE itself. Instead of the standard residual \[\partial_t u_\theta -\Delta u_\theta +v(u_\theta) -f_{\bar\theta}\mathbf{1}_\omega,\] we consider weighted residuals of the form \[\phi_{\theta^1}\partial_t u_\theta -\sum_{i=1}^d \phi_{\bar\theta^i}\partial_{x_i x_i}u_\theta +\phi_{\theta^2} v(u_\theta) -\phi_{\theta^3} f_{\bar\theta}\mathbf{1}_\omega .\] For bilinear control, the control component is instead weighted as \[-\phi_{\theta^3} f_{\bar\theta}u_\theta .\] Thus, the time derivative, the directional diffusion terms, the nonlinear response, and the control term are weighted independently. The same principle is applied to the boundary, initial, and terminal constraints. The resulting method is therefore not merely a scalar reweighting of a PINN loss. It is an operator-decomposed and control-aware adaptive residual method.

1.3 Structural features of WeightedPINN↩︎

The position of WeightedPINN relative to existing PINN methodologies has been discussed in Subsection 1.1. We now summarize the structural features that, taken together, define the proposed framework and justify its use for control reconstruction.

(i) Operator-level adaptivity. In a standard PINN formulation, the PDE residual is treated as a single quantity, so that the time derivative, the diffusion operator, the nonlinear term, and the control term are all penalized through the same residual norm. This can lead to imbalance during training, particularly in nonlinear or stiff regimes. WeightedPINN introduces space–time-dependent weights at the level of each operator component, so that the differential, nonlinear, and control contributions can be rebalanced individually during training.

(ii) Control-sensitive residual metric. In control reconstruction, the unknown control is not observed directly. It is identified through the state equation together with the initial, boundary, and terminal constraints. Therefore, a small averaged PDE residual does not necessarily imply an accurate reconstructed control. This issue is especially delicate in bilinear control problems, where the control appears through the product \(f_{\bar\theta}u_\theta\) and where errors in the state and errors in the control interact multiplicatively. WeightedPINN addresses this difficulty by assigning an independent adaptive weight to the control term, making the loss more sensitive to the part of the residual through which the control is identified.

(iii) Directional and space–time weighting. The Laplacian is decomposed into its directional components \(\partial_{x_i x_i}u_\theta\), and each component is weighted independently. This produces an anisotropic residual metric that distinguishes errors associated with different spatial directions. Because the weights are also functions of space and time, the method can emphasize regions where the approximation is less accurate or where the control is more difficult to recover, for instance near the control region, near transition layers, or in portions of the domain where the state is small and control identification becomes ill-conditioned.

(iv) Penalized min–max formulation consistent with the original PINN. The state and control networks minimize the weighted residuals, while the weight networks amplify poorly resolved components of the equation and of the constraints. A quadratic penalty keeps the weights close to \(1\), preventing the learned weighted operator from drifting too far from the original PDE. In particular, when all weights are equal to \(1\), the WeightedPINN loss reduces exactly to the standard PINN loss. Thus, the standard PINN is contained as the limiting case in which the penalty dominates, and the proposed method may be viewed as a controlled adaptive perturbation of the standard PINN formulation.

(v) Structured treatment of control constraints. The weighted formulation is also applied to the boundary, initial, and terminal conditions. This is essential in controllability problems, because the terminal condition is not auxiliary data; it is the target condition that defines the control problem. A neural pair \((u_\theta,f_{\bar\theta})\) that satisfies the PDE residual but fails to reach the prescribed terminal state cannot be regarded as a reliable reconstructed control. By weighting both the dynamical residual and the control constraints in a structured manner, the method enforces the governing equation and the target condition more effectively.

Overall, the main contribution of this work is not simply the introduction of adaptive weights into PINNs. Rather, we introduce an operator-decomposed and control-aware WeightedPINN framework for nonlinear PDE control reconstruction. The method is new at the methodological level because the adaptive weights act on individual differential, nonlinear, control, and constraint components. It is also new in the control setting, because previous adaptive PINN methodologies have not been systematically developed for the numerical reconstruction of controls for nonlinear PDEs of the type considered here.

The WeightedPINN method is introduced in Section 2, and its theoretical analysis is presented in Section 3. The numerical results in Section 4 indicate that WeightedPINN improves over the standard PINN baseline on the high-dimensional control benchmarks considered in this work. In dimension \(d=10\), the method produces smaller residual-based testing errors for the internal and bilinear control problems governed by both heat and wave equations. The improvement is most pronounced in the equation residual, where the adaptive weighting mechanism frequently reduces the PDE error substantially, while keeping the boundary and endpoint errors comparable or smaller.

The lower-dimensional manufactured-solution benchmarks in Section 5 further show that WeightedPINN attains the smallest mean reconstruction errors for the state and the control among the methods tested. The numerical evidence suggests that a principal advantage of WeightedPINN, in the benchmarks considered here, appears in regimes where the residual contains componentwise imbalances, such as singular perturbations, anisotropic or variable-coefficient operators, and control-identification mechanisms in which the unknown control enters through a localized or multiplicative term. This is consistent with the design of the method: different components of the PDE residual are weighted separately, allowing the training procedure to respond to componentwise residual imbalances that are not explicitly treated by scalar loss reweighting, adaptive sampling, or causal weighting.

Acknowledgments↩︎

The authors would like to thank Emmanuel Trélat for fruitful discussions on the topic.

2 Control problems, PINNs and WeightedPINNs↩︎

We first recall the definition of a feedforward neural network mapping \(\mathbb{R}^d\) to \(\mathbb{R}\) with width \(\mathcal{W}\) and depth \(\mathcal{L}\). The width structure of the network is specified by a sequence \((\mathcal{W}_i)_{i=0}^{\mathcal{L}+1}\) such that \(\mathcal{W}_0 = d\), \(\mathcal{W}_{\mathcal{L}+1} = 1\), and \(\mathcal{W} := \max_{1\le i\le\mathcal{L}}\mathcal{W}_i\).

Given parameters \(\theta = (w_0, b_0, \dots, w_{\mathcal{L}}, b_{\mathcal{L}})\), where \(w_i \in \mathbb{R}^{\mathcal{W}_{i+1} \times \mathcal{W}_i}\) and \(b_i \in \mathbb{R}^{\mathcal{W}_{i+1}}\), we define the affine transformations \[\begin{align} \ell_i(x) = w_i x + b_i, \quad i = 0, \dots, \mathcal{L}. \end{align}\] A neural network \(F_\theta\) with activation function \(\sigma : \mathbb{R} \to \mathbb{R}\) is defined by \[\begin{align} F_{\theta}(x) := \ell_{\mathcal{L}} \circ \sigma_{\mathcal{L}} \circ \ell_{\mathcal{L}-1}\circ \cdots \circ \sigma_{1} \circ \ell_0(x), \end{align}\] where, for each \(i = 1, \dots, \mathcal{L}\), the map \(\sigma_i : \mathbb{R}^{\mathcal{W}_i} \to \mathbb{R}^{\mathcal{W}_i}\) denotes the componentwise application of \(\sigma\), i.e., \[\begin{align} \sigma_{i}\big((x_1, \dots, x_{\mathcal{W}_i})^\top\big) = \big(\sigma(x_1), \dots, \sigma(x_{\mathcal{W}_i})\big)^\top. \end{align}\]

Suppose we aim to solve the partial differential equation (PDE) \[\label{eq:general95PDE} \mathcal{D}u(x) = f(x), \quad \text{in } Q, \qquad \mathcal{B}u(x) = g(x), \quad \text{on } \Gamma,\tag{1}\] where \(x \in Q\) denotes the spatial variable, and may also include a temporal component \(t\) when appropriate. Here, \(\mathcal{D}u = f\) represents the governing PDE, and \(\mathcal{B}u = g\) represents the associated boundary and/or initial conditions. The sets \(Q\) and \(\Gamma\) denote the corresponding domain and boundary, respectively.

In a standard PINN approach, we seek a neural network \(u(x;\theta)\) that approximates the solution \(u\) of 1 . The residuals associated with the PDE and the boundary/initial conditions are defined by \[\mathcal{R}_{Q}(u(x;\theta)) := \mathcal{D}u(x;\theta) - f(x), \qquad \mathcal{R}_{\Gamma}(u(x;\theta)) := \mathcal{B}u(x;\theta) - g(x).\]

The PDE can then be solved by seeking parameters \(\theta\) that minimize the residual errors [43]: \[\min_{\theta} \left\|\mathcal{R}_{Q}(u(x;\theta))\right\|_{Q}^{2} + \lambda \left\|\mathcal{R}_{\Gamma}(u(x;\theta))\right\|_{\Gamma}^{2},\] where \(\|\cdot\|\) typically denotes the \(L^{2}\)-norm, and \(\lambda > 0\) is a weighting parameter. For example, this minimization problem can be written as \[\label{StandardNET} \min_{\theta} \mathbb{E}_{x \in Q} \left[ \left|\mathcal{D}u(x;\theta) - f(x)\right|^{2} \right] + \lambda \, \mathbb{E}_{x \in \Gamma} \left[ \left|\mathcal{B}u(x;\theta) - g(x)\right|^{2} \right].\tag{2}\]

In practice, the \(L^2\)-norm is often used for simplicity, although the PDE solution may require higher regularity.

Solving the minimization problem 2 yields the standard Physics-Informed Neural Network (PINN) formulation.

Below, we consider control problems for semilinear heat and wave equations to illustrate our method. First, we formally state the control problems and then apply the PINN method to approximate the associated control. Next, our new WeightedPINN method introduces adaptability to the learning process. We consider both the internal control for the heat and wave equations and the bilinear control for these equations. We introduce the notation \(\mathcal{R}_Q^{\mathcal{I}},\mathcal{R}_\Gamma^{\mathcal{I}}\), with \(\mathcal{I}\in \{1,2,3,4\}\), to distinguish residuals corresponding to different PDE and control settings. Similarly, we denote by \(\bar{\mathcal{R}}_Q^{\mathcal{I}},\bar{\mathcal{R}}_\Gamma^{\mathcal{I}}\) for the residuals in WeightedPINN settings. We refer to Section [Sec:internalsettings] and Section [Sec:bilinearsettings] for the detailed definitions.

2.1 Internal control problems for heat and wave equations in high dimensions↩︎

  

Let \(\Omega \subset \mathbb{R}^d\) be a bounded domain with smooth boundary, let \(\omega \subset \Omega\) be an open subset, and let \(T>0\). We define \[Q_T := (0,T) \times \Omega, \quad q_T := (0,T) \times \omega, \quad \Sigma_T := (0,T) \times \partial\Omega,\] and \[\Gamma := \Sigma_T \cup (\{0,T\} \times \Omega).\] We denote by \(\mathbf{1}_\omega\) the characteristic function of \(\omega\).

For a given globally Lipschitz function \(v:\mathbb{R}\to\mathbb{R}\), we consider the semilinear heat equation \[\begin{align} \label{Heat1} \partial_t u - \Delta u + v(u) &= f \mathbf{1}_\omega \quad \text{in } Q_T, \nonumber\\ u &= 0 \quad \text{on } \Sigma_T, \\ u(0) &= u_0 \quad \text{in } \Omega, \nonumber \end{align}\tag{3}\] where \(u_0 \in L^2(\Omega)\).

Equation 3 is said to be exactly controllable at time \(T\) if, for any \(u_0 \in L^2(\Omega)\) and \(z_0 \in L^2(\Omega)\), there exists a control \(f \in L^2(q_T)\) such that the corresponding solution satisfies \(u(T)=z_0\).

The existence of such controls is not guaranteed in general; see [44], [45]. In this work, we assume that a control exists for 3 for given initial and terminal conditions.

We also consider the semilinear wave equation \[\begin{align} \label{Wave1} \partial_{tt} u - \Delta u + v(u) &= f \mathbf{1}_\omega \quad \text{in } Q_T, \nonumber\\ u &= 0 \quad \text{on } \Sigma_T, \\ (u(0), \partial_t u(0)) &= (u_0,u_1) \quad \text{in } \Omega, \nonumber \end{align}\tag{4}\] where \(v:\mathbb{R} \to \mathbb{R}\) is globally Lipschitz and \((u_0,u_1)\in H^1(\Omega)\times L^2(\Omega)\).

The wave equation 4 is said to be exactly controllable at time \(T\) if, for any initial data \((u_0,u_1)\in H^1(\Omega)\times L^2(\Omega)\) and target data \((z_0,z_1)\in H^1(\Omega)\times L^2(\Omega)\), there exists a control \(f \in L^2(q_T)\) such that \((u(T), \partial_t u(T)) = (z_0,z_1).\)

We assume the existence of such controls for 4 under given initial and terminal conditions. We refer to [46], [47] for relevant conditions. Specifically, we discuss the conditions further in Remark 13.

• First approach: standard PINN method

  Hereafter, let \(u_\theta\) and \(f_{\bar{\theta}}\) denote neural networks with parameters \(\theta\) and \(\bar{\theta}\) that approximate, respectively, the solution and the control for the model equation.

  Heat equation: We consider the control problem 3 with initial condition \(u_0\) and terminal condition \(z_0\). First, we recall the definition of \(H^{1/2}\)-norm on \(\Sigma_T\), which is \[\begin{align} \|\varphi\|^2_{H^{1/2}(\Sigma_T)}:=\|\varphi\|^2_{L^2(\Sigma_T)}+[\varphi]^2_{H^{1/2}(\Sigma_T)}, \end{align}\] where \[\begin{align} [\varphi]^2_{H^{1/2}(\Sigma_T)}:=\int_{\Sigma_T}\int_{\Sigma_T}\frac{(\varphi(x)-\varphi(y))^2}{\|x-y\|^{d+1}}dS_x dS_y. \end{align}\] The PINN method is associated with the minimization formulation \[\begin{align} \inf_{\theta,\bar\theta}\;&\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_t u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} \mathbf{1}_\omega \right)^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|u_\theta\|_{H^{1/2}(\Sigma_T)}^2 +\|u_\theta(0,\cdot)-u_0\|_{L^2(\Omega)}^2 +\|u_\theta(T,\cdot)-z_0\|_{L^2(\Omega)}^2 \Big), \end{align}\] where \(\lambda>0\) is a given constant, \(|Q_T|\) and \(|\Gamma|\) denote the \((d+1)\)- and \(d\)-dimensional measures of \(Q_T\) and \(\Gamma\), respectively.

  Remark 2. From a theoretical perspective, even though Equation 3 has homogeneous boundary conditions, the approximation \(u_\theta\) is non-homogeneous. Thus, \(u_\theta\) must possess sufficient regularity to approximate the solution of the heat equation.

  In [48], when we need a weak solution in \(L^2(0,T;H^1(\Omega))\), the trace space for parabolic equation is at least \[L^2(0,T;H^{1/2}(\partial \Omega))\cap H^{1/4}(0,T;L^2(\partial \Omega)).\]

  For the ease of the formulation, we utilized the stronger \(H^{1/2}\)-norm on the boundary in the PINN formulation, this norm controls the anisotropic norm on the trace space. In practice, however, we adopt only the \(L^2\)-norm for simplicity.

  We define the equation residual by \[\begin{align} \|\mathcal{R}_{Q}^{1}(u_{\theta},f_{\bar\theta})\|_{Q}^2 :=\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_t u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} \mathbf{1}_\omega \right)^2 dx\,dt, \end{align}\] and define the boundary residual by \[\begin{align} \|\mathcal{R}_{\Gamma}^{1}(u_{\theta})\|_{\Gamma}^2 :=\frac{1}{|\Gamma|} \Big( \|u_\theta\|_{H^{1/2}(\Sigma_T)}^2 +\|u_\theta(0,\cdot)-u_0\|_{L^2(\Omega)}^2 +\|u_\theta(T,\cdot)-z_0\|_{L^2(\Omega)}^2 \Big). \end{align}\]

  The standard PINN formulation for the internal control on the heat equation, with \(\lambda > 0\), can be written in the following compact form: \[\begin{align} \inf_{\theta,\bar{\theta}} \|\mathcal{R}_{Q}^{1}(u_{\theta},f_{\bar{\theta}})\|_{Q}^2 +\lambda \|\mathcal{R}_{\Gamma}^{1}(u_{\theta})\|_{\Gamma}^2. \end{align}\]

  Wave equation: We now consider the control problem in 4 with initial condition \((u_0,u_1)\) and terminal condition \((z_0,z_1)\). The definition of \(H^{3/2}\)-norm on \(\Sigma_T\) is as follow \[\begin{align} \|\varphi\|^2_{H^{3/2}(\Sigma_T)}:=\|\varphi\|^2_{L^2(\Sigma_T)}+\|\partial_t \varphi\|_{H^{1/2}(\Sigma_T)}^2+\sum_{i=1}^{d-1}\bigl\|\frac{\partial\varphi}{\partial \sigma_i}\bigr\|_{H^{1/2}(\Sigma_T)}^2, \end{align}\] where \((\sigma_1,\dots,\sigma_{d-1})\) is a smooth local coordinate of \(\partial \Omega\).

  The PINN method in this control problem is associated with the minimization formulation \[\begin{align} \inf_{\theta,\bar\theta}\;&\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_{tt} u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} \mathbf{1}_\omega \right)^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|u_\theta\|_{H^{3/2}(\Sigma_T)}^2 +\|(u_\theta(0,\cdot),\partial_t u_\theta(0,\cdot))-(u_0,u_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \nonumber\\ &\qquad\quad +\|(u_\theta(T,\cdot),\partial_t u_\theta(T,\cdot))-(z_0,z_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \Big), \end{align}\] where \(\lambda>0\) is a given constant.

  Remark 3. In the same spirit as Remark 2, one would expect the boundary conditions has higher regularity than \(L^2\). For wave equation, we choose \(H^{3/2}\)-norm for boundary conditions.

  Furthermore, the initial and terminal data \(u_0\) and \(z_0\) belong to \(H^1(\Omega)\). Therefore, from a theoretical standpoint, the approximation should be carried out in the \(H^1\)-norm.

  In practice, however, we adopt the \(L^2\)-norm for boundary, initial and terminal conditions.

  The equation residual is defined by \[\begin{align} \|\mathcal{R}_{Q}^{2}(u_{\theta},f_{\bar\theta})\|_{Q}^2 :=\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_{tt} u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} \mathbf{1}_\omega \right)^2 dx\,dt, \end{align}\] and the boundary residual is defined by \[\begin{align} \|\mathcal{R}_{\Gamma}^{2}(u_{\theta})\|_{\Gamma}^2 &:=\frac{1}{|\Gamma|} \Big( \|u_\theta\|_{H^{3/2}(\Sigma_T)}^2 +\|(u_\theta(0,\cdot),\partial_t u_\theta(0,\cdot))-(u_0,u_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \nonumber\\ &\qquad\quad +\|(u_\theta(T,\cdot),\partial_t u_\theta(T,\cdot))-(z_0,z_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \Big). \end{align}\]

  The standard PINN formulation for the internal control on the wave equation, with \(\lambda > 0\), can be written in the following compact form: \[\begin{align} \inf_{\theta,\bar\theta} \|\mathcal{R}_{Q}^{2}(u_{\theta},f_{\bar\theta})\|_{Q}^2 +\lambda \|\mathcal{R}_{\Gamma}^{2}(u_{\theta})\|_{\Gamma}^2. \end{align}\]

• Second approach: WeightedPINN method

  We introduce additional neural networks \(\phi_{\theta}:\mathbb{R}^{d+1}\to \mathbb{R}\) that act as adaptive weights. We impose on the weight networks a uniform boundedness assumption, which will be specified later in [con:weight].

  Heat equation: We again consider the control problem 3 with initial condition \(u_0\) and terminal condition \(z_0\). The WeightedPINN method is associated with the following min–max formulation: \[\begin{align} \label{eq:lossheat} \inf_{\theta,\bar\theta}\,\sup_{(\theta^i)_{i=1}^{8},(\bar\theta^i)_{i=1}^d}\; &\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_t u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} \mathbf{1}_\omega \Big]^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|_{L^2(\Omega)}^2\nonumber\\ &\qquad\quad + \|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^{8}}(T,\cdot)z_0\|_{L^2(\Omega)}^2 \Big)\nonumber\\ &-\rho\Big( \sum_{i=1}^{8}\|\phi_{\theta^i}-1\|^2_{Q^i} +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)} \Big)\nonumber, \end{align}\tag{5}\] where \(\lambda,\rho>0\) are given constants, \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8\).

  In the proposed formulation, the additional neural networks serve as adaptive weights, designed to maximize discrepancies in the equation and boundary residuals, thereby enhancing the effectiveness of the minimization. To regulate these weights, we introduce additional \(L^2\)-norm penalization terms. These terms ensure that the weights remain controlled and do not deviate significantly from the constant function \(1\).

  Remark 4. The weights \(\phi_{\theta^5}\) and \(\phi_{\theta^6}\) can be viewed as functions from \(\mathbb{R}^d\) to \(\mathbb{R}\), i.e., without explicit dependence on \(t\) at \(t=0\). For simplicity, we employ the same network architecture as for the other weight functions, namely neural networks defined on \(\mathbb{R}^{d+1}\). These networks are then restricted to \(t=0\).

  The weights \(\phi_{\theta^7}, \phi_{\theta^{8}}\) are treated analogously at \(t=T\).

  As in the standard PINN framework, we adopt a compact notation for the WeightedPINN formulation by introducing weighted residuals. For readability, we omit their explicit dependence on the neural network parameters.

  We define the weighted residual for the heat equation by \[\begin{align} \label{eq:barRQ1} \|\bar{\mathcal{R}}_{Q}^{1}\|_{Q}^2 :=\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_t u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} \mathbf{1}_\omega \Big]^2 dx\,dt. \end{align}\tag{6}\]

  The weighted boundary residual is \[\begin{align} \label{eq:barRG1} \|\bar{\mathcal{R}}_{\Gamma}^{1}\|_{\Gamma}^2 &:=\frac{1}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|_{L^2(\Omega)}^2\\ &\qquad\quad + \|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^{8}}(T,\cdot)z_0\|_{L^2(\Omega)}^2 \Big).\nonumber \end{align}\tag{7}\]

  In WeightedPINN approach, we also define the auxiliary neural residual by \[\begin{align} \label{eq:RNN1} \|\mathcal{R}_{NN}^{1}\|_{NN}^2 :=\sum_{i=1}^{8}\|\phi_{\theta^i}-1\|^2_{Q^i} +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)}. \end{align}\tag{8}\] Here, we recall \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8\).

  The compact form for WeightedPINN formulation for the internal control on the heat equation, with \(\lambda,\rho>0\), is given by \[\begin{align} \inf_{\theta,\bar\theta}\,\sup_{(\theta^i)_{i=1}^8,(\bar\theta^i)_{i=1}^d} \;\|\bar{\mathcal{R}}_{Q}^{1}\|_{Q}^2 +\lambda \|\bar{\mathcal{R}}_{\Gamma}^{1}\|_{\Gamma}^2 -\rho \|\mathcal{R}_{NN}^{1}\|_{NN}^2. \end{align}\]

  Wave equation: We extend the idea of WeightedPINN to control problems of the wave equation. Let us consider 4 with initial condition \((u_0,u_1)\) and \((z_0,z_1)\). The WeightedPINN is associated with the min–max formulation \[\begin{align} \label{eq:losswave} \inf_{\theta,\bar\theta}\,\sup_{(\theta^i)_{i=1}^{12},(\bar\theta^i)_{i=1}^d,(\tilde{\theta}^i)_{i=1}^d,(\hat{\theta}^i)_{i=1}^{2d},(\check\theta^i)_{i=1}^{2d}}\; &\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_{tt}u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} \mathbf{1}_\omega \Big]^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{L^2(\Sigma_T)}^2 +\|\phi_{\tilde{\theta}^1}\partial_t u_\theta\|^2_{H^{1/2}(\Sigma_T)} +\sum_{i=2}^{d}\bigl\|\phi_{\tilde{\theta}^i}\frac{\partial u_\theta}{\partial \sigma_{i-1}}\bigr\|^2_{H^{1/2}(\Sigma_T)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^{9}}(0,\cdot)\partial_t u_{\theta}(0,\cdot)-\phi_{\theta^{10}}(0,\cdot)u_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^{8}}(T,\cdot)z_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad+\|\phi_{\theta^{11}}(T,\cdot)\partial_t u_{\theta}(T,\cdot)-\phi_{\theta^{12}}(T,\cdot)z_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\hat{\theta}^{2i-1}}(0,\cdot)\partial_{x_i}u_{\theta}(0,\cdot)-\phi_{\hat{\theta}^{2i}}(0,\cdot)\partial_{x_i}u_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\check{\theta}^{2i-1}}(T,\cdot)\partial_{x_i}u_{\theta}(T,\cdot)-\phi_{\check{\theta}^{2i}}(T,\cdot)\partial_{x_i}z_0\|^2_{L^2(\Omega)}\Big)\nonumber\\ &-\rho\Big( \sum_{i=1}^{12}\|\phi_{\theta^i}-1\|^2_{Q^i}\\ &\qquad\quad +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)}+\frac{1}{|\Gamma|}\sum_{i=1}^{d}\|\phi_{\tilde{\theta}^i}-1\|^2_{L^2(\Sigma_T)}\nonumber\\ &\qquad\quad +\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\hat{\theta}^i}(0,\cdot)-1\|^2_{L^2(\Omega)} +\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\check{\theta}^i}(T,\cdot)-1\|^2_{L^2(\Omega)}\nonumber \Big), \end{align}\tag{9}\] where \(\lambda,\rho>0\) are given constants, \((\sigma_1,\dots,\sigma_{d-1})\) is a smooth local coordinate of \(\partial\Omega\), \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6,9,10\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8,11,12\).

  Similar to the control problem for heat equation, as in Remark 4, we take \(\phi_{\theta^5},\phi_{\theta^6},\phi_{\theta^{9}},\phi_{\theta^{10}}\) and \(\phi_{\hat{\theta}^i}\) to be neural networks from \(\mathbb{R}^{d+1}\) to \(\mathbb{R}\), restricted to \(t = 0\). We treat \(\phi_{\theta^7},\phi_{\theta^{8}},\phi_{\theta^{11}},\phi_{\theta^{12}}\) and \(\phi_{\check\theta^i}\) analogously, with restriction to \(t = T\).

  For the internal control on the wave equation, the weighted residuals are defined by \[\begin{align} \label{eq:barRQ2} \|\bar{\mathcal{R}}_{Q}^{2}\|_{Q}^2 :=\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_{tt}u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} \mathbf{1}_\omega \Big]^2 dx\,dt, \end{align}\tag{10}\] and \[\begin{align} \label{eq:barRG2} \|\bar{\mathcal{R}}_{\Gamma}^2\|_{\Gamma}^{2} &:=\frac{1}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{L^2(\Sigma_T)}^2 +\|\phi_{\tilde{\theta}^0}\partial_t u_\theta\|^2_{H^{1/2}(\Sigma_T)} +\sum_{i=1}^{d}\|\phi_{\tilde{\theta}^i}\partial_{x_i}u_\theta\|^2_{H^{1/2}(\Sigma_T)}\\ &\qquad\quad +\|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|^2_{L^2(\Omega)} +\|\phi_{\theta^{9}}(0,\cdot)\partial_t u_{\theta}(0,\cdot)-\phi_{\theta^{10}}(0,\cdot)u_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^{8}}(T,\cdot)z_0\|^2_{L^2(\Omega)} +\|\phi_{\theta^{11}}(T,\cdot)\partial_t u_{\theta}(T,\cdot)-\phi_{\theta^{12}}(T,\cdot)z_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\hat{\theta}^{2i-1}}(0,\cdot)\partial_{x_i}u_{\theta}(0,\cdot)-\phi_{\hat{\theta}^{2i}}(0,\cdot)\partial_{x_i}u_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\check{\theta}^{2i-1}}(T,\cdot)\partial_{x_i}u_{\theta}(T,\cdot)-\phi_{\check{\theta}^{2i}}(T,\cdot)\partial_{x_i}z_0\|^2_{L^2(\Omega)}\Big).\nonumber \end{align}\tag{11}\]

  We also define the neural residual by \[\begin{align} \label{eq:RNN2} \|\mathcal{R}_{NN}^{2}\|_{NN}^2 &:=\sum_{i=1}^{12}\|\phi_{\theta^i}-1\|^2_{Q^i} +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)} +\frac{1}{|\Gamma|}\sum_{i=1}^{d}\|\phi_{\tilde{\theta}^i}-1\|^2_{L^2(\Sigma_T)}\\ &\quad+\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\hat{\theta}^i}(0,\cdot)-1\|^2_{L^2(\Omega)} +\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\check{\theta}^i}(T,\cdot)-1\|^2_{L^2(\Omega)},\nonumber \end{align}\tag{12}\] where \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6,9,10\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8,11,12\).

  Finally, for given parameters \(\lambda,\rho>0\), the WeightedPINN approach for the wave equation has the compact form \[\begin{align} \inf_{\theta,\bar{\theta}}\, \sup_{(\theta^i)_{i=1}^{12},(\bar\theta^i)_{i=1}^d,(\tilde{\theta}^i)_{i=1}^{d},(\hat{\theta}^i)_{i=1}^{2d},(\check\theta^i)_{i=1}^{2d}} \;\|\bar{\mathcal{R}}_Q^{2}\|_{Q}^2+\lambda\|\bar{\mathcal{R}}_{\Gamma}^{2}\|_{\Gamma}^2-\rho\|\mathcal{R}_{NN}^{2}\|_{NN}^2. \end{align}\]

2.2 Bilinear control problems for heat and wave equations in high dimensions↩︎

  

We use the same notations for the spatial domains and boundaries as in the internal control setting. For completeness, we recall them below.

Let \(\Omega \subset \mathbb{R}^d\) be a bounded domain with smooth boundary, and let \(T>0\). We define \[Q_T := (0,T) \times \Omega, \quad \Sigma_T := (0,T) \times \partial\Omega, \quad \Gamma := \Sigma_T \cup (\{0,T\} \times \Omega).\]

Let \(v:\mathbb{R}\to\mathbb{R}\) be a globally Lipschitz function. We consider the heat equation with bilinear control: \[\begin{align} \label{Heat1b} \partial_t u - \Delta u + v(u) &= f u \quad \text{in } Q_T, \nonumber\\ u &= 0 \quad \text{on } \Sigma_T, \\ u(0) &= u_0 \quad \text{in } \Omega, \nonumber \end{align}\tag{13}\]

We say that 13 is bilinearly controllable at time \(T\) from \(u_0 \in L^2(\Omega)\) to a target \(z_0 \in L^2(\Omega)\) if there exists a control \(f \in L^\infty(Q_T)\) such that the corresponding solution satisfies \(u(T)=z_0\). Sufficient conditions for such controllability are given in [49].

We next consider the wave equation with the globally Lipschitz function \(v:\mathbb{R} \to \mathbb{R}\) and bilinear control: \[\begin{align} \label{Wave1b} \partial_{tt} u - \Delta u + v(u) &= f u \quad \text{in } Q_T, \nonumber\\ u &= 0 \quad \text{on } \Sigma_T, \\ (u(0), \partial_t u(0)) &= (u_0,u_1) \quad \text{in } \Omega. \nonumber \end{align}\tag{14}\]

We say that 14 is bilinearly controllable at time \(T\) from \((u_0,u_1) \in H^1(\Omega)\times L^2(\Omega)\) to \((z_0,z_1) \in H^1(\Omega)\times L^2(\Omega)\) if there exists a control \(f \in L^\infty(Q_T)\) such that the solution satisfies \((u(T), \partial_t u(T)) = (z_0,z_1).\) Relevant results on bilinear controllability of the wave equation can be found in [50], [51].

• First approach: standard PINN method

  Heat equation: The standard PINN formulation of bilinear control on the heat equation is given by \[\begin{align} \inf_{\theta,\bar\theta}\;&\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_t u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} u_\theta \right)^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|u_\theta\|_{H^{1/2}(\Sigma_T)}^2 +\|u_\theta(0,\cdot)-u_0\|_{L^2(\Omega)}^2 +\|u_\theta(T,\cdot)-z_0\|_{L^2(\Omega)}^2 \Big), \end{align}\] where \(\lambda>0\) is a given constant. We recall that \(|Q_T|\) and \(|\Gamma|\) denote the \((d+1)\)- and \(d\)-dimensional measures of \(Q_T\) and \(\Gamma\), respectively.

  The equation residual is given by \[\begin{align} \|\mathcal{R}_{Q}^{3}(u_{\theta},f_{\bar\theta})\|_{Q}^2 :=\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_t u_{\theta} - \Delta u_{\theta} + v(u_\theta) - f_{\bar\theta} u_{\theta} \right)^2 dx\,dt, \end{align}\] and the boundary residual is \[\begin{align} \|\mathcal{R}_{\Gamma}^{3}(u_{\theta})\|_{\Gamma}^2 :=\frac{1}{|\Gamma|} \Big( \|u_\theta\|_{H^{1/2}(\Sigma_T)}^2 +\|u_\theta(0,\cdot)-u_0\|_{L^2(\Omega)}^2 +\|u_\theta(T,\cdot)-z_0\|_{L^2(\Omega)}^2 \Big). \end{align}\]

  The standard training formulation for bilinear controls on the heat equation, with \(\lambda>0\), is given in compact form by \[\begin{align} \inf_{\theta,\bar\theta} \|\mathcal{R}_{Q}^{3}(u_{\theta},f_{\bar\theta})\|_{Q}^2 +\lambda \|\mathcal{R}_{\Gamma}^{3}(u_{\theta})\|_{\Gamma}^2. \end{align}\]

  Wave equation: Similar to the heat equation, the PINN formulation for bilinear control on the wave equation is given by \[\begin{align} \inf_{\theta,\bar\theta}\;&\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_{tt} u_{\theta} - \Delta u_{\theta} + v(u_{\theta}) - f_{\bar\theta} u_\theta \right)^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|u_\theta\|_{H^{3/2}(\Sigma_T)}^2 +\|(u_\theta(0,\cdot),\partial_t u_\theta(0,\cdot))-(u_0,u_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \nonumber\\ &\qquad\quad +\|(u_\theta(T,\cdot),\partial_t u_\theta(T,\cdot))-(z_0,z_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \Big), \end{align}\] where \(\lambda>0\) is a given constant.

  The residuals are defined as follows: \[\begin{align} \|\mathcal{R}_{Q}^{4}(u_{\theta},f_{\bar\theta})\|_{Q}^2 &:=\frac{1}{|Q_T|}\int_{Q_T} \left( \partial_{tt} u_{\theta} - \Delta u_{\theta} + v(u_\theta) - f_{\bar\theta} u_{\theta} \right)^2 dx\,dt, \end{align}\] and \[\begin{align} \|\mathcal{R}_{\Gamma}^{4}(u_{\theta})\|_{\Gamma}^2 &:=\frac{1}{|\Gamma|} \Big( \|u_\theta\|_{H^{3/2}(\Sigma_T)}^2 +\|(u_\theta(0,\cdot),\partial_t u_\theta(0,\cdot))-(u_0,u_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \\ &\qquad\quad +\|(u_\theta(T,\cdot),\partial_t u_\theta(T,\cdot))-(z_0,z_1)\|_{H^1(\Omega)\times L^2(\Omega)}^2 \Big).\nonumber \end{align}\]

  The compact form of PINN formulation for bilinear control on the wave equation is given by \[\begin{align} \inf_{\theta,\bar\theta} \|\mathcal{R}_{Q}^{4}(u_{\theta},f_{\bar\theta})\|_{Q}^2 +\lambda \|\mathcal{R}_{\Gamma}^{4}(u_{\theta})\|_{\Gamma}^2. \end{align}\]

• Second approach: WeightedPINN method. We further define a weighted formulation of the residual errors for the heat and wave equations in the context of bilinear control problems. Similar to internal control settings, for readability, we omit the dependence of the residuals on the parameters of the related neural networks.

  Heat equation. We consider the min–max formulation that defines the WeightedPINN for bilinear control in 13 : \[\begin{align} \inf_{\theta,\bar\theta}\,\sup_{(\theta^i)_{i=1}^8,(\bar\theta^i)_{i=1}^d}\; &\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_t u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} u_\theta \Big]^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|_{L^2(\Omega)}^2\\ &\qquad\quad + \|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^8}(T,\cdot)z_0\|_{L^2(\Omega)}^2 \Big)\\ &-\rho\Big( \sum_{i=1}^{8}\|\phi_{\theta^i}-1\|^2_{Q^i} +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)} \Big), \end{align}\] where \(\lambda,\rho>0\) are given constants, \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8\).

  To shorten the formulation, we define \[\begin{align} \label{eq:barRQ3} \|\bar{\mathcal{R}}_Q^3\|_Q^2 :=\frac{1}{|Q_T|}\int_{Q_T}\bigg[ \phi_{\theta^1}\partial_{t}u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^{i}}\partial_{x_ix_i} u_{\theta} + \phi_{\theta^2}v(u_\theta) - \phi_{\theta^{3}}f_{\bar\theta} u_{\theta} \bigg]^2\,dx\,dt. \end{align}\tag{15}\]

  Next, we define the weighted boundary residual: \[\begin{align} \label{eq:barRG3} \|\bar{\mathcal{R}}_{\Gamma}^3\|_{\Gamma}^2 &:=\frac{1}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|_{L^2(\Omega)}^2\\ &\quad\quad + \|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^8}(T,\cdot)z_0\|_{L^2(\Omega)}^2 \Big).\nonumber \end{align}\tag{16}\]

  The auxiliary neural residual is \[\begin{align} \label{eq:RNN3} \|\mathcal{R}_{NN}^{3}\|_{NN}^2 :=\sum_{i=1}^{8}\|\phi_{\theta^i}-1\|^2_{Q^i} +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)}, \end{align}\tag{17}\] where \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8\).

  Thus, the compact form of the WeightedPINN in this case is \[\begin{align} \inf_{\theta,\bar{\theta}}\, \sup_{(\theta^i)_{i=1}^8,(\bar\theta^i)_{i=1}^d} \|\bar{\mathcal{R}}_Q^{3}\|_{Q}^2 +\lambda\|\bar{\mathcal{R}}_{\Gamma}^{3}\|_{\Gamma}^2 -\rho\|\mathcal{R}_{NN}^{3}\|_{NN}^2. \end{align}\]

  Wave equation. We consider the min–max formulation that defines the WeightedPINN for bilinear control in 14 : \[\begin{align} \inf_{\theta,\bar\theta}\,\sup_{(\theta^i)_{i=1}^{12},(\bar\theta^i)_{i=1}^d,(\tilde{\theta}^i)_{i=1}^d,(\hat{\theta}^i)_{i=1}^{2d},(\check\theta^i)_{i=1}^{2d}}\; &\frac{1}{|Q_T|}\int_{Q_T} \Big[ \phi_{\theta^1}\partial_{tt}u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^i}\partial_{x_i x_i} u_{\theta} + \phi_{\theta^{2}} v(u_{\theta}) - \phi_{\theta^{3}} f_{\bar\theta} u_\theta \Big]^2 dx\,dt\\ &+\frac{\lambda}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{L^2(\Sigma_T)}^2 + \|\phi_{\tilde{\theta}^1}\partial_t u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \sum_{i=2}^{d} \bigl\|\phi_{\tilde{\theta}^i}\frac{\partial u_\theta}{\partial \sigma_{i-1}}\bigr\|_{H^{1/2}(\Sigma_T)}^2\\ &\qquad\quad +\|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|^2_{L^2(\Omega)}\\ &\qquad\quad +\|\phi_{\theta^{9}}(0,\cdot)\partial_t u_{\theta}(0,\cdot)-\phi_{\theta^{10}}(0,\cdot)u_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^8}(T,\cdot)z_0\|^2_{L^2(\Omega)}\\ &\qquad\quad +\|\phi_{\theta^{11}}(T,\cdot)\partial_t u_{\theta}(T,\cdot)-\phi_{\theta^{12}}(T,\cdot)z_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\hat{\theta}^{2i-1}}(0,\cdot)\partial_{x_i}u_{\theta}(0,\cdot)-\phi_{\hat{\theta}^{2i}}(0,\cdot)\partial_{x_i}u_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\check{\theta}^{2i-1}}(T,\cdot)\partial_{x_i}u_{\theta}(T,\cdot)-\phi_{\check{\theta}^{2i}}(T,\cdot)\partial_{x_i}z_0\|^2_{L^2(\Omega)}\Big)\nonumber\\ &-\rho\Big( \sum_{i=1}^{12}\|\phi_{\theta^i}-1\|^2_{Q^i}\\ &\qquad\quad +\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)} +\frac{1}{|\Gamma|}\sum_{i=1}^{d}\|\phi_{\tilde{\theta}^i}-1\|^2_{L^2(\Sigma_T)}\\ &\qquad\quad +\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\hat{\theta}^i}(0,\cdot)-1\|^2_{L^2(\Omega)} +\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\check{\theta}^i}(T,\cdot)-1\|^2_{L^2(\Omega)}\nonumber \Big), \end{align}\] where \(\lambda,\rho>0\) are given constants, \((\sigma_1,\dots,\sigma_{d-1})\) is a smooth local coordinate of \(\partial\Omega\), \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6,9,10\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8,11,12\).

  To write a compact form for the formulation, we define \[\begin{align} \label{eq:barRQ4} \|\bar{\mathcal{R}}_Q^4\|_Q^2 :=\frac{1}{|Q_T|}\int_{Q_T}\bigg[ \phi_{\theta^1}\partial_{tt}u_{\theta} - \sum_{i=1}^{d}\phi_{\bar{\theta}^{i}}\partial_{x_ix_i} u_{\theta} + \phi_{\theta^2}v(u_\theta) - \phi_{\theta^{3}}f_{\bar\theta} u_{\theta} \bigg]^2\,dx\,dt. \end{align}\tag{18}\]

  The boundary residual for the wave equation is given by \[\begin{align} \label{eq:barRG4} \|\bar{\mathcal{R}}_{\Gamma}^4\|_{\Gamma}^2 &:=\frac{1}{|\Gamma|} \Big( \|\phi_{\theta^4}u_\theta\|_{L^2(\Sigma_T)}^2 + \|\phi_{\tilde{\theta}^0}\partial_t u_\theta\|_{H^{1/2}(\Sigma_T)}^2 + \sum_{i=1}^{d} \|\phi_{\tilde{\theta}^i}\partial_{x_i}u_\theta\|_{H^{1/2}(\Sigma_T)}^2\\ &\qquad\quad +\|\phi_{\theta^5}(0,\cdot)u_{\theta}(0,\cdot)-\phi_{\theta^6}(0,\cdot)u_0\|^2_{L^2(\Omega)} +\|\phi_{\theta^{9}}(0,\cdot)\partial_t u_{\theta}(0,\cdot)-\phi_{\theta^{10}}(0,\cdot)u_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\|\phi_{\theta^7}(T,\cdot)u_{\theta}(T,\cdot)-\phi_{\theta^8}(T,\cdot)z_0\|^2_{L^2(\Omega)} +\|\phi_{\theta^{11}}(T,\cdot)\partial_t u_{\theta}(T,\cdot)-\phi_{\theta^{12}}(T,\cdot)z_1\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\hat{\theta}^{2i-1}}(0,\cdot)\partial_{x_i}u_{\theta}(0,\cdot)-\phi_{\hat{\theta}^{2i}}(0,\cdot)\partial_{x_i}u_0\|^2_{L^2(\Omega)}\nonumber\\ &\qquad\quad +\sum_{i=1}^{d}\|\phi_{\check{\theta}^{2i-1}}(T,\cdot)\partial_{x_i}u_{\theta}(T,\cdot)-\phi_{\check{\theta}^{2i}}(T,\cdot)\partial_{x_i}z_0\|^2_{L^2(\Omega)}\Big).\nonumber \end{align}\tag{19}\]

  For the wave equation, the auxiliary neural residual is given by \[\begin{align} \label{eq:RNN4} \|\mathcal{R}_{NN}^{4}\|_{NN}^2 &:=\sum_{i=1}^{12}\|\phi_{\theta^i}-1\|^2_{Q^i}+\frac{1}{|Q_T|}\sum_{i=1}^{d}\|\phi_{\bar{\theta}^i}-1\|^2_{L^2(Q_T)}+\frac{1}{|\Gamma|}\sum_{i=1}^{d}\|\phi_{\tilde{\theta}^i}-1\|^2_{L^2(\Sigma_T)}\\ &\quad+\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\hat{\theta}^i}(0,\cdot)-1\|^2_{L^2(\Omega)}+\frac{1}{|\Gamma|}\sum_{i=1}^{2d}\|\phi_{\check{\theta}^i}(T,\cdot)-1\|^2_{L^2(\Omega)},\nonumber \end{align}\tag{20}\] where \(\|\cdot\|^2_{Q^i}=\frac{\|\cdot\|^2_{L^2(Q_T)}}{|Q_T|}\) for \(i=1,2,3\), \(\|\cdot\|^2_{Q^4}=\frac{\|\cdot\|^2_{L^2(\Sigma_T)}}{|\Gamma|}\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(0,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=5,6,9,10\), \(\|\phi\|^2_{Q^i}=\frac{\|\phi(T,\cdot)\|^2_{L^2(\Omega)}}{|\Gamma|}\) for \(i=7,8,11,12\).

  Finally, the WeightedPINN formulation for bilinear control on the wave equation takes the compact form

  \[\begin{align} \inf_{\theta,\bar{\theta}}\, \sup_{(\theta^i)_{i=1}^{12},(\bar\theta^i)_{i=1}^d,(\tilde{\theta}^i)_{i=1}^d,(\hat{\theta}^i)_{i=1}^{2d},(\check\theta^i)_{i=1}^{2d}} \|\bar{\mathcal{R}}_Q^{4}\|_{Q}^2 +\lambda\|\bar{\mathcal{R}}_{\Gamma}^{4}\|_{\Gamma}^2 -\rho\|\mathcal{R}_{NN}^{4}\|_{NN}^2. \end{align}\]

2.3 Illustrative examples in two space dimensions↩︎

We present a representative example in two space dimensions to illustrate the behavior of the standard PINN method and the proposed WeightedPINN formulation for bilinear problems. The emphasis is on how the residual structure is modified and how this impacts the reconstruction of the control.

2.3.1 Example: Two-dimensional heat equation with bilinear control↩︎

Let \(Q_T := (0,T)\times \Omega\), where \(\Omega=(0,1)^2\). For simplicity, we take \(v \equiv 0\) in this example. We consider the bilinear control problem associated with the heat equation \[\begin{align} \partial_t u - \Delta u = f(x)\,u, \qquad (t,x)\in Q_T, \end{align}\] supplemented with homogeneous boundary condition \[\begin{align} u = 0 \quad \text{on } \Sigma_T := (0,T)\times \partial\Omega, \end{align}\] and prescribed initial and terminal data \[\begin{align} u(0,x)=u_0(x), \qquad u(T,x)=z_0(x). \end{align}\]

In the standard PINN approach, one introduces neural network approximations \(u_\theta\) and \(f_{\bar\theta}\), and minimizes the loss \[\begin{align} \mathcal{L}_{\mathrm{PINN}}(\theta,\bar\theta) = \frac{1}{|Q_T|}\int_{Q_T} \left( \partial_t u_\theta - \Delta u_\theta - f_{\bar\theta} u_\theta \right)^2\,dx\,dt +\lambda \|\mathcal{R}_\Gamma(u_\theta)\|^2_{\Gamma}, \end{align}\] where \[\begin{align} \|\mathcal{R}_\Gamma(u_\theta)\|^2_{\Gamma} = \frac{1}{|\Gamma|} \Big( \|u_\theta\|^2_{H^{1/2}(\Sigma_T)} +\|u_\theta(0,\cdot)-u_0\|^2_{L^2(\Omega)} +\|u_\theta(T,\cdot)-z_0\|^2_{L^2(\Omega)} \Big). \end{align}\] This corresponds to a uniform least-squares enforcement of the PDE and boundary conditions. In this formulation, the residual is measured globally over the whole space-time domain, so the training process is driven by an averaged error criterion. As a consequence, the approximation may fit the state reasonably well while still leaving noticeable inaccuracies in the recovered control, especially in our control settings where the dependence on the coefficient \(f\) is sensitive.

In the WeightedPINN approach, one introduces additional neural networks \(\phi\) acting as adaptive weights. In two space dimensions, the weighted PDE residual becomes \[\begin{align} \bar{\mathcal{R}}_Q = \phi_{\theta^1}\,\partial_t u_\theta -\phi_{\bar\theta^1}\,\partial_{x_1x_1}u_\theta -\phi_{\bar\theta^2}\,\partial_{x_2x_2}u_\theta -\phi_{\theta^3}\,f_{\bar\theta}u_\theta. \end{align}\]

We emphasize that, in two space dimensions, the Laplacian operator is decomposed into its directional components, namely \[\Delta u_\theta = \partial_{x_1x_1} u_\theta + \partial_{x_2x_2} u_\theta,\] and the WeightedPINN formulation assigns independent weights to each second-order derivative. This leads to the terms \[-\phi_{\bar\theta^1}\,\partial_{x_1x_1}u_\theta \quad \text{and} \quad -\phi_{\bar\theta^2}\,\partial_{x_2x_2}u_\theta,\] which are treated separately in the residual.

This directional weighting provides additional flexibility compared to a single scalar weight in front of the Laplacian. In particular, it allows the training procedure to distinguish anisotropic features of the approximation error, as the contributions of the second-order derivatives along the \(x_1\) and \(x_2\) directions may differ significantly across the domain. As a result, the maximization step can selectively amplify discrepancies associated with specific spatial directions, thereby improving the sensitivity of the loss to localized and direction-dependent errors.

From a functional viewpoint, this corresponds to replacing the standard isotropic \(L^2\) residual norm by an adaptive anisotropic metric, in which each differential component of the operator is weighted independently. This refinement is particularly relevant in control problems, where the identifiability of the coefficient \(f\) may depend differently on variations along different spatial directions.

In the same spirit, the weighted formulation for the boundary, initial, and terminal conditions is given by \[\begin{align} \|\bar{\mathcal{R}}_\Gamma\|_{\Gamma}^2 &= \frac{1}{|\Gamma|}\bigg( \|\phi_{\theta^4} u_\theta\|^2_{H^{1/2}(\Sigma_T)} + \|\phi_{\tilde{\theta}^1} \partial_{x_1} u_\theta\|^2_{L^2(\Sigma_T)} + \|\phi_{\tilde{\theta}^2} \partial_{x_2} u_\theta\|^2_{L^2(\Sigma_T)} \\ &\quad\quad + \|\phi_{\theta^5}(0,\cdot) u_{\theta}(0,\cdot) - \phi_{\theta^6}(0,\cdot) u_{0}\|_{L^2(\Omega)}^2 \\ &\quad\quad + \|\phi_{\theta^7}(T,\cdot) u_{\theta}(T,\cdot) - \phi_{\theta^8}(T,\cdot) z_{0}\|_{L^2(\Omega)}^2 \bigg). \end{align}\]

The corresponding loss is embedded into the min–max problem \[\begin{align} \inf_{\theta,\bar\theta}\sup_{(\theta^i)_{i=1}^8,(\bar\theta^i)_{i=1}^2} \Bigg( &\|\bar{\mathcal{R}}_Q\|_Q^2 +\lambda \|\bar{\mathcal{R}}_\Gamma\|_\Gamma^2 \\ &-\rho \Big( \frac{1}{|Q_T|}\big(\|\phi_{\theta^1}-1\|^2_{L^2(Q_T)} +\|\phi_{\theta^3}-1\|^2_{L^2(Q_T)}\\ &\quad\quad +\|\phi_{\bar\theta^1}-1\|^2_{L^2(Q_T)} +\|\phi_{\bar\theta^2}-1\|^2_{L^2(Q_T)}\big)\\ &\quad\quad +\frac{1}{|\Gamma|}\big(\|\phi_{\theta^4}-1\|^2_{L^2(\Sigma_T)} +\|\phi_{\theta^5}(0,\cdot)-1\|^2_{L^2(\Omega)} +\|\phi_{\theta^6}(0,\cdot)-1\|^2_{L^2(\Omega)}\\ &\quad\quad +\|\phi_{\theta^7}(T,\cdot)-1\|^2_{L^2(\Omega)} +\|\phi_{\theta^8}(T,\cdot)-1\|^2_{L^2(\Omega)} \big) \Big) \Bigg). \end{align}\]

When \(\|\bar{\mathcal{R}}_Q\|_Q=\|\bar{\mathcal{R}}_\Gamma\|_\Gamma=0\) and \[\phi_{\theta^1}=\phi_{\bar\theta^1}=\phi_{\bar\theta^2}=\phi_{\theta^3}=\phi_{\theta^4}=\phi_{\theta^5}=\phi_{\theta^6}=\phi_{\theta^7}=\phi_{\theta^8}=1,\] the formulation reduces to the original PDE. This formulation changes the role of the residual in the training process. Instead of being measured in a fixed norm, it is evaluated in an adaptive metric determined by the weight networks. The maximization step increases the influence of those regions where the equation is less accurately satisfied, while the penalization term prevents the weights from drifting too far from \(1\) when such amplification is unnecessary. In this way, the method allocates more attention to the most informative parts of the domain and produces a more accurate recovery of both the state and the control.

3 Theoretical error estimate↩︎

3.1 PINNs and WeightedPINNs for internal control problems↩︎

The internal control problem for heat and wave equations has been extensively studied in the literature. Representative results can be found in [44]–[47]. The main result of this section establishes a convergence estimate for the WeightedPINN approximation in terms of the network architecture and the residual errors.

We exploit the universal approximation properties of neural networks to construct numerical approximations of the control. To this end, we first introduce the necessary definitions and assumptions for the neural network framework.

We recall the Heaviside step function \[\begin{align} \mathcal{H}(x) = \begin{cases} 1, & x>0,\\ 0, & x\le 0. \end{cases} \end{align}\]

A function \(\sigma \in W^{m,\infty}_{\mathrm{loc}}(\mathbb{R})\) is called \(m\)-th order quasi-decay if there exist constants \(C_1, C_2 > 0\) such that, for every \(x\neq 0\), \[\begin{align} |\sigma(x)-\mathcal{H}(x)| \le \min\left\{\frac{C_1}{|x|},\, C_2\right\}, \quad |\sigma^{(k)}(x)| \le \min\left\{\frac{C_1}{|x|^{k+1}},\, C_2\right\}, \quad 1 \le k \le m. \end{align}\]

Furthermore, we say that \(\sigma\) exhibits nonlinear behavior if there exist \(a \in \mathbb{R}\) and \(\delta > 0\) such that \[\begin{align} \sigma''(a) \neq 0, \quad \text{and} \quad \sigma \in C^\infty( (a-\delta, a+\delta) ). \end{align}\]

We impose the following assumption on the activation function \(\sigma\):

1. \(\sigma \in W^{m,\infty}_{\mathrm{loc}}(\mathbb{R})\) exhibits nonlinear behavior, and either \(\sigma(x)\) or \(\sigma(x)/x\) is \(m\)-th order quasi-decaying.

In [52], several examples satisfying Assumption [con:sigma] are discussed. These include sigmoid, HardSigmoid, tanh, HardTanh, Softsign, Soft-root-sign (SRS), ELU, SELU, CELU, GELU, SiLU, dSiLU, Softplus, and Mish.

The weight neural networks \(\phi_{\theta}\) may have flexible architectures in the sense that we can use many different activation functions to construct \(\phi_{\theta}\). We require them and their derivatives to be uniformly bounded. Specifically, we assume

1. There exists \(M>0\) such that \[\sup_{(t,x)\in\overline{Q_T}} \bigl\{|\phi_\theta(t,x)|+|\partial_t\phi_\theta(t,x)|+\sum_{i=1}^{d}|\partial_{x_i}\phi_\theta(t,x)| \bigr\}\le M\] for any parameter \(\theta\) in a given parameter set \(\Theta\).

As examples, we consider the weight networks with \(3\) hidden layers. On each hidden layer, we take a \(C^1\) activation function, which can be sigmoid, tanh, or \(\text{ReLU}^{m+1},m\ge 1\). The parameter set \(\Theta\) is assumed to be compact subset of \[\mathbb{R}^{1\times \mathcal{W}_3}\times\mathbb{R}^{1}\times \dots\times\mathbb{R}^{\mathcal{W}_1\times d}\times\mathbb{R}^{\mathcal{W}_1},\] which ensures the boundedness condition in [con:weight].

We now extend the previous result to the wave equation.

Theorem 1 and Theorem 2 establish the existence of deep neural networks that approximate both the state and the control given sufficient regularity. In particular, the results guarantee that the total WeightedPINN loss can be made arbitrarily small, provided the network width and depth are chosen sufficiently large.

Now, let us further assume a multiplicative structure on 3 . For the heat equation, we assume

1. There exists \(V\in L^{\infty}(Q_T)\) such that the control problem on the heat equation under consideration is \[\begin{align} \partial_t u - \Delta u + V u &= f \mathbf{1}_{\omega} \quad \text{in }Q_T,\\ u &= 0 \quad \text{on }\Sigma_T,\\ u(0) &= u_0 \quad \text{in } \Omega,\\ u(T) &= z_0 \quad \text{in } \Omega. \end{align}\]

In other words, we simply replace \(v(u)\) in Equation 3 and the Residual 6 by \(V u\).

In addition to Theorem 1, we show in Theorem 3 that, under the multiplicative structure, if we can find the neural networks \(u_\theta\) and \(f_{\bar{\theta}}\) with sufficiently small WeightedPINN loss, then these networks approximate an admissible pair of solution and control of the PDE.

For brevity, we denote the space of admissible pairs of solution and control associated with heat equation with given initial and terminal conditions \(u_0,z_0\in L^2(\Omega)\) by \(\mathcal{A}_{u_0}^{z_0}\).

As discussed in Section [Sec:internalsettings], \(\mathcal{A}_{u_0}^{z_0}\) is not guaranteed to be non-empty. We assume \(\mathcal{A}_{u_0}^{z_0}\) to be non-empty. This also implies that \(z_0\) is reachable with some control \(f\mathbf{1}_\omega\). If \(f_{\bar\theta}\) approximates \(f\), then we can expect to have an approximate control from \(u_0\) to \(z_0\). We only keep the terminal \(z_0\) to be exact. To make it precise, in this work, for \(V\) in [con:Vheat], we assume:

1. For given neural network \(f_{\bar \theta}\) and \(\delta>0\), there exist \(F_\delta\in L^2(q_T),u_\delta \in L^2(\Omega)\) and \[\bar u\in C([0,T];L^2(\Omega))\cap L^2(0,T;H^1_0(\Omega))\] such that \[\|F_\delta\|_{L^2(q_T)}^2 + \|u_\delta\|_{L^2(\Omega)}^2\le \delta,\] and \[\begin{align} \partial_t\bar u-\Delta\bar u+V\bar u &= (f_{\bar\theta}+F_\delta)\mathbf{1}_\omega \quad\text{in }Q_T,\\ \bar u &= 0 \quad\text{on }\Sigma_T,\\ \bar u(0)&=u_0+u_\delta \quad\text{in }\Omega,\\ \bar u(T)&=z_0\quad\text{in }\Omega. \end{align}\]

We note that \(\mathcal{A}_{u_0}^{z_0}\neq\emptyset\) and the fact that \(f_{\bar\theta}\) approximates a control \(f\) do not directly implies [con:heat]. Here, it stays as a structural hypothesis. We do not claim [con:heat] is a necessary condition.

Now, let \(\bar u\in C([0,T];L^2(\Omega))\cap L^2(0,T;H^1(\Omega))\) and \(f\in L^2(q_T)\), the distance of \((\bar u,\bar f)\) to \(\mathcal{A}_{u_0}^{z_0}\) is given by \[\begin{align} {\rm dist}_h( (\bar u,\bar f),\mathcal{A}_{u_0}^{z_0}):= \inf_{(u,f)\in \mathcal{A}_{u_0}^{z_0}}\Big(\max_{t\in [0,T]}\| u(t)-\bar u(t)\|^2_{L^2(\Omega)} +\|u-\bar u\|^2_{L^2(0,T;H^1(\Omega))} +\|f-\bar f\|^2_{L^2(q_T)}\Big)^{1/2}. \end{align}\]

We extend the result to the wave equation 4 . The multiplicative structure on the wave equation is written as follows:

1. There exists \(V\in L^{\infty}(Q_T)\) such that the control problem on the wave equation is \[\begin{align} \partial_{tt} u - \Delta u + V u &= f \mathbf{1}_{\omega} \quad \text{in }Q_T,\\ u &= 0 \quad \text{on }\Sigma_T,\\ (u(0),\partial_t u(0)) &= (u_0,u_1) \quad \text{in } \Omega,\\ (u(T),\partial_t u(T)) &= (z_0,z_1) \quad \text{in } \Omega. \end{align}\]

We replace \(v(u)\) in Equation 4 and the Residual 10 by \(V u\).

The space of admissible pairs of solution and control associated with 4 with given initial and terminal conditions \((u_0,u_1),(z_0,z_1) \in H^1(\Omega)\times L^2(\Omega)\) is denoted by \(\mathcal{A}_{(u_0,u_1)}^{(z_0,z_1)}\).

Let \(\bar u \in C([0,T];H^1(\Omega))\cap C^1([0,T];L^2(\Omega))\) and \(\bar f\in L^2(q_T)\), the distance of \((\bar u,\bar f)\) to \(\mathcal{A}_{(u_0,u_1)}^{(z_0,z_1)}\) is given by \[\begin{align} &{\rm dist}_w( (\bar u,\bar f),\mathcal{A}_{(u_0,u_1)}^{(z_0,z_1)})\\ &\quad\quad:= \inf_{(u,f)\in \mathcal{A}_{(u_0,u_1)}^{(z_0,z_1)}}\Big(\max_{t\in [0,T]}\{\| u(t)-\bar u(t)\|^2_{H^1(\Omega)} +\|\partial_t u(t)-\partial_t \bar u(t)\|^2_{L^2(\Omega)}\} +\|f-\bar f\|^2_{L^2(q_T)}\Big)^{1/2}. \end{align}\]

We consider the following assumptions:

1. The domain \(\omega\subset\Omega\) contains a neighborhood of the boundary \(\partial \Omega\).

2. \(T > \text{diameter} (\Omega)\).

3.2 PINNs and WeightedPINNs for bilinear controls↩︎

  

We extend the results from Theorem 1 and Theorem 2 to bilinear control settings.

We now state the corresponding results for bilinear control of the wave equation.

4 Numerical results: High dimensional Problems↩︎

4.1 Set up↩︎

We fix the dimension \(d = 10\). We consider several configurations of the domain \(\Omega\) and the controlled region \(\omega\).

We set the constant \(\lambda = \rho = 1\) in the formulation \[\|\bar{\mathcal{R}}_Q\|_{Q}^2+\lambda\|\bar{\mathcal{R}}_{\Gamma}\|_{\Gamma}^2-\rho\|\mathcal{R}_{NN}\|_{NN}^2.\]

In the numerical experiments, derivatives and the Laplace operator are approximated using Finite Difference (FD) schemes with step size \(h = 0.001\). We denotes the time derivative, second-order derivative and the Laplace in the FD schemes by \(\bar\partial_t,\bar\partial_{tt},\bar{\Delta}_x\), respectively.

Now, we set \[v(u) := -u + u^2.\] The numerical experiments are conducted in a locally Lipschitz regime. The theoretical assumptions are not directly satisfied by this choice of \(v\), but the experiments remain stable because the learned states remain uniformly bounded.

• Internal control:

  • Heat equation: We set \(T = 1\) and define \[\label{eq:heat95disc} \mathscr{P}_h u := \bar{\partial}_t u - \bar{\Delta}_x u + v(u).\tag{46}\] The discrete equation is \[\mathscr{P}_h u = f \mathbf{1}_\omega.\]

  • Wave equation: We set \(T = 3\) and define \[\label{eq:wave95disc} \mathscr{P}_w u := \bar{\partial}_{tt} u - \bar{\Delta}_x u + v(u).\tag{47}\] The discrete equation is \[\mathscr{P}_w u = f \mathbf{1}_\omega.\]

• Bilinear control:

  • Heat equation: We set \(T = 1\) and use \(\mathscr{P}_h\) as in 46 . The discrete equation is \[\mathscr{P}_h u = f u.\]

  • Wave equation: We set \(T = 3\) and use \(\mathscr{P}_w\) as in 47 . The discrete equation is \[\mathscr{P}_w u = f u.\]

We next describe the architectures of the neural networks.

We recall the ReLU activation function \(\sigma:\mathbb{R}\to[0,\infty)\), defined by \[\begin{align} \label{ReLU} \mathrm{ReLU}(x) := \max\{0,x\}. \end{align}\tag{48}\]

To ensure twice differentiability, we use the activation function \(\mathrm{ReLU}^3\). This function, denoted by \(\sigma_3\), is defined by \[\begin{align} \label{ReLU3} \sigma_3(x) := \mathrm{ReLU}(x)^3. \end{align}\tag{49}\]

We also recall the sigmoid activation function \(\sigma_s:\mathbb{R}\to(0,1)\), defined by \[\begin{align} \label{sigmoid} \sigma_s(x) := \frac{1}{1+e^{-x}}. \end{align}\tag{50}\]

We consider two types of neural networks with different activation functions: one for \(u_\theta\) and \(f_{\bar{\theta}}\), and the other for the weight networks \(\phi_{\theta^i}\).

• Solution and control networks: These networks have width \(\mathcal{W} = 100\) and depth \(\mathcal{L} = 3\), and use \(\sigma_3\) as the activation function.

• Weight networks: These networks have width \(\mathcal{W} = 40\) and depth \(\mathcal{L} = 3\). We use \(\sigma_3\) as the activation function in the first two hidden layers and \(\sigma_s\) in the final hidden layer. This structure allows us to control the range of the weight networks. In the numerical experiments, the output of the sigmoid layer is rescaled affinely to \([0.2,5]\) via \(\phi=0.2+4.8\sigma_s(\cdot)\). The aforementioned structure satisfies the condition mentioned in Remark 6.

Finally, we describe the training and testing procedures.

Training procedures: We consider eight different configurations. For each configuration, we perform \(10000\) training iterations. At each iteration, the following steps are carried out:

1. We randomly generate \(N_1 = 1000\) training points \((t_i,x_i)\) in the interior of the domain, and \(N_1\) training points \((\hat{t}_i,\hat{x}_i)\) on the boundary. We also generate \(N_2 = 100\) training points \((0,\bar{x}_i)\) for the initial condition and \(N_2\) training points \((T,\tilde{x}_i)\) for the terminal condition.

2. These training points are used to evaluate the loss functions in 5 and 9 , which correspond to discrete approximations of \[\|\bar{\mathcal{R}}_Q\|_{Q}^2 + \|\bar{\mathcal{R}}_\Gamma\|_{\Gamma}^2 - \|\mathcal{R}_{NN}\|_{NN}^2,\] where \(\bar{\mathcal{R}}_Q\), \(\bar{\mathcal{R}}_\Gamma\), and \(\mathcal{R}_{NN}\) are defined in Sections [Sec:internalsettings] and [Sec:bilinearsettings]. For simplicity, we use the \(L^2\)-norm formulation of the loss.

3. We then employ the ADAM optimizer implemented in the torch package to solve the resulting min–max optimization problem. We set the learning rate at \(0.001\).

Testing errors: We randomly generate \(N_3 = 10000\) test points in the interior and on the boundary, and \(N_4 = 1000\) test points for the initial and terminal conditions. These test points are used to evaluate the approximation error in the numerical experiments.

Presenting results: We display the training loss as a function of the number of iterations and report the testing error for each numerical experiment. In the high-dimensional tests of Sections 4.2 and 4.3, the learned control and the corresponding state are computed but are not displayed, due to the high-dimensional nature of the problem. The reported testing errors in these sections are therefore residual-based quantities, measuring the PDE residual together with the boundary, initial, and terminal constraint errors. In the lower-dimensional tests of Section 5, we instead evaluate the reconstruction errors directly against known manufactured solutions.

4.2 Internal control problems↩︎

4.2.1 Heat equation↩︎

For the heat equation, the loss 5 is numerically evaluated using training points and the \(H^{1/2}\)-norm is reduced to \(L^2\)-norm.

For the standard PINN formulation, we set \(\phi_{\theta^j},\phi_{\bar\theta^j} \equiv 1\) for all \(j\).

Using the test points, we numerically evaluate the error as \[\begin{align} H_h &= \|\mathscr{P}_h u - f \mathbf{1}_\omega\|_{L^2([0,1]\times\Omega)} + \|u\|_{L^2([0,1]\times\partial\Omega)} + \|u(0,\cdot) - u_0\|_{L^2(\Omega)} + \|u(T,\cdot) - z_0\|_{L^2(\Omega)}. \end{align}\]

For convenience, we refer to \(\|\mathscr{P}_h u - f \mathbf{1}_\omega\|_{L^2([0,1]\times\Omega)}\) as the equation error, and the remaining terms as the boundary error.

1. Situation 1: The domain \(\Omega\) is the unit ball in \(\mathbb{R}^d\), and the control region is given by \[\omega = \Big(-\sqrt{\frac{1}{d}-\frac{1}{2d^2}},\sqrt{\frac{1}{d}-\frac{1}{2d^2}}\Big)^{d}.\]

   We consider the initial condition \(u_0 = 0\) and the terminal condition \[z_0 = (e-1)\sin\!\left(\frac{\pi}{2}(1-|x|)^{2.5}\right).\]

   In Figure 1, we plot the total training loss as a function of the number of iterations for both methods (blue: PINN, red: WeightedPINN).

   

   We observe in Figure 1 that the initial loss is at \(10^3\) and decreases to approximately \(10^{-3}\) by the \(10000\)-th iteration, demonstrating the effectiveness of both the standard PINN and WeightedPINN.

   We further observe that the WeightedPINN achieves a lower total training loss compared to the standard PINN.

   Table 1 reports the test errors at the final iteration.

   Table 1: Test errors at the final iteration in Situation 1

   	Total	Eqn.	Bnd.
   Standard	\(6.082\times10^{-2}\)	\(2.004\times10^{-2}\)	\(4.078\times10^{-2}\)
   WeightedPINN	\(4.411\times10^{-2}\)	\(3.425\times10^{-3}\)	\(4.068\times10^{-2}\)

   From Table 1, both methods achieve relatively small overall errors. However, the WeightedPINN significantly reduces the equation error, while the boundary error remains comparable between the two approaches.

2. Situation 2: The domain is \(\Omega = (-1,1)^d\), and the control region \(\omega\) is the unit ball in \(\mathbb{R}^d\). We consider the initial condition \[u_0 = \prod_{i=1}^{d} (1 - |x_i|),\] and the terminal condition \[z_0 = e \prod_{i=1}^{d} (1 - |x_i|).\]

   Similar to Situation 1, WeightedPINN achieves improved numerical performance in this setting.

   In Figure 2, we plot the total training loss as a function of the number of iterations for both methods (blue: PINN, red: WeightedPINN).

   

   We observe in Figure 2 that the initial loss is at \(10^3\). The loss of the standard PINN method decreases to approximately \(10^{-3}\) by the \(10000\)-th iteration and the loss of the WeightedPINN method reaches approximately \(10^{-4}\) by the \(10000\)-th iteration. This again shows the effectiveness of both methods, especially the WeightedPINN.

   Table 2 reports the test errors at the final iteration.

   Table 2: Test errors at the final iteration in Situation 2

   	Total	Eqn.	Bnd.
   Standard	\(5.114\times10^{-2}\)	\(2.733\times10^{-2}\)	\(2.382\times10^{-2}\)
   WeightedPINN	\(2.299\times10^{-2}\)	\(3.171\times10^{-3}\)	\(1.982\times10^{-2}\)

4.2.2 Wave equation↩︎

For the wave equation, the loss 9 is numerically evaluated using the training points and the \(H^{3/2}\)-norm on \(\Sigma_T\), the \(H^1\)-norms on \(u_0,z_0\) are reduced to \(L^2\)-norms.

Similarly to Section 4.2.1, for the wave equation we numerically evaluate the error as \[\begin{align} H_w &= \|\mathscr{P}_w u - f \mathbf{1}_{\omega}\|_{L^2([0,T]\times\Omega)} + \|u\|_{L^2([0,T]\times\partial\Omega)} \\ &\quad + \|(u(0,\cdot),\bar{\partial}_t u(0,\cdot)) - (u_0,u_1)\|_{L^2(\Omega)\times L^2(\Omega)} + \|(u(T,\cdot),\bar{\partial}_t u(T,\cdot)) - (z_0,z_1)\|_{L^2(\Omega)\times L^2(\Omega)} \end{align}\] We refer to \(\|\mathscr{P}_w u - f \mathbf{1}_{\omega}\|_{L^2([0,T]\times\Omega)}\) as the equation error, and to the remaining terms collectively as the boundary error.

1. Situation 3: The domain \(\Omega\) is the unit ball in \(\mathbb{R}^d\), and the control region \(\omega\) is defined such that \(\Omega \setminus \omega\) is the closed ball centered at \((1/d,\dots,1/d)\) with radius \(3/4\). We choose \(\omega\) so that Assumption [con:wave] is satisfied. For simplicity, we further assume that \(\omega\) contains the entire boundary \(\partial\Omega\).

   We consider the initial and terminal conditions \[\begin{align} u_0 &= u_1 = 0, \\ z_0 &= (e-1)\sin\!\left(\frac{\pi}{2}(1-|x|)^{2.5}\right), \\ z_1 &= 2e\,\sin\!\left(\frac{\pi}{2}(1-|x|)^{2.5}\right). \end{align}\]

   Figure 3 shows the total training loss in Situation 3 (blue: PINN, red: WeightedPINN).

   

   In this situation, the WeightedPINN requires more iterations to stabilize; however, it eventually achieves a lower training loss compared to the standard PINN.

   We observe in Figure 3 that the loss of the standard PINN method can only reach approximately \(10^{-2}\) at \(10000\)-th iteration. On the other hand, the loss of WeightedPINN can reach approximately \(10^{-4}\) at \(10000\)-th iteration.

   Table 3 reports the test errors at the final iteration.

   Table 3: Test errors at the final iteration in Situation 3

   	Total	Eqn.	Bnd.
   Standard	\(1.745\times10^{-1}\)	\(3.611\times10^{-2}\)	\(1.384\times10^{-1}\)
   WeightedPINN	\(4.833\times10^{-2}\)	\(3.505\times10^{-3}\)	\(4.482\times10^{-2}\)

   In Situation 3, the WeightedPINN demonstrates significant improvements in both the equation error and the boundary error.

2. Situation 4: The domain is \(\Omega = (-1,1)^d\), and the control region \(\omega\) is defined such that \[\Omega \setminus \omega = \Big[-\sqrt{\frac{1}{d}-\frac{1}{2d^2}},\sqrt{\frac{1}{d}-\frac{1}{2d^2}}\Big]^d.\]

   We consider the initial and terminal conditions: \[\begin{align} u_0 &= u_1 = \prod_{i=1}^{d}(1 - |x_i|), \\ z_0 &= z_1 = e \prod_{i=1}^{d}(1 - |x_i|). \end{align}\]

   Figure 4 shows the total training loss in Situation 4 (blue: PINN, red: WeightedPINN).

   

   In figure 4, after \(10000\) iterations, the loss of the standard PINN method reaches approximately \(10^{-3}\) and the loss of the WeightedPINN method reaches approximately \(10^{-4}\).

   Table 4 reports the test errors at the final iteration.

   Table 4: Test errors at the final iteration in Situation 4

   	Total	Eqn.	Bnd.
   Standard	\(7.561\times10^{-2}\)	\(3.047\times10^{-2}\)	\(4.514\times10^{-2}\)
   WeightedPINN	\(4.225\times10^{-2}\)	\(2.566\times10^{-3}\)	\(3.968\times10^{-2}\)

   The WeightedPINN consistently achieves better numerical performance in this setting.

4.3 Bilinear control problems↩︎

The loss functions and error metrics for the bilinear control problems are identical to those in Section 4.2, with the only difference being that \(f_{\bar{\theta}}(t_i,x_i)\mathbf{1}_\omega(x_i)\) is replaced by \(f_{\bar{\theta}}(t_i,x_i)u_\theta(t_i,x_i)\) throughout. We test the same four situations with the same domain and initial/terminal conditions. For brevity, we directly report the numerical results.

4.3.1 Heat equation↩︎

1. Situation 5: We consider the same domain \(\Omega\), initial condition \(u_0\), and terminal condition \(z_0\) as in Situation 1 of Section 4.2.1.

   Figure 5 shows the total training loss in Situation 5 (blue: PINN, red: WeightedPINN).

   

   Table 5 reports the test errors at the final iteration.

   Table 5: Test errors at the final iteration in Situation 5

   	Total	Eqn.	Bnd.
   Standard	\(4.801\times10^{-2}\)	\(2.572\times10^{-2}\)	\(2.228\times10^{-2}\)
   WeightedPINN	\(2.263\times10^{-2}\)	\(2.780\times10^{-3}\)	\(1.985\times10^{-2}\)

2. Situation 6: We consider the same domain \(\Omega\), initial condition \(u_0\), and terminal condition \(z_0\) as in Situation 2 of Section 4.2.1.

   Figure 6 shows the total training loss in Situation 6 (blue: PINN, red: WeightedPINN).

   

   Table 6 reports the test errors at the final iteration.

   Table 6: Test errors at the final iteration in Situation 6

   	Total	Eqn.	Bnd.
   Standard	\(8.053\times10^{-2}\)	\(3.106\times10^{-2}\)	\(4.948\times10^{-2}\)
   WeightedPINN	\(4.358\times10^{-2}\)	\(2.950\times10^{-3}\)	\(4.063\times10^{-2}\)

4.3.2 Wave equation↩︎

1. Situation 7: We consider the same domain \(\Omega\), initial condition, and terminal condition as in Situation 3 of Section 4.2.2.

   Figure 7 shows the total training loss in Situation 7 (blue: PINN, red: WeightedPINN).

   

   Table 7 reports the test errors at the final iteration.

   Table 7: Test errors at the final iteration in Situation 7

   	Total	Eqn.	Bnd.
   Standard	\(1.986\times10^{-1}\)	\(2.350\times10^{-2}\)	\(1.751\times10^{-1}\)
   WeightedPINN	\(1.755\times10^{-1}\)	\(1.005\times10^{-2}\)	\(1.655\times10^{-1}\)

2. Situation 8: We consider the same domain \(\Omega\), initial condition, and terminal condition as in Situation 4 of Section 4.2.2.

   Figure 8 shows the total training loss in Situation 8 (blue: PINN, red: WeightedPINN).

   

   Table 8 reports the test errors at the final iteration.

   Table 8: Test errors at the final iteration in Situation 8

   	Total	Eqn.	Bnd.
   Standard	\(7.538\times10^{-2}\)	\(3.024\times10^{-2}\)	\(4.514\times10^{-2}\)
   WeightedPINN	\(4.258\times10^{-2}\)	\(2.824\times10^{-3}\)	\(3.975\times10^{-2}\)

Conclusion: Under the settings of our numerical experiments, both the standard PINN and the WeightedPINN achieve small total losses, demonstrating the effectiveness of the proposed approach for solving the considered control problems.

Moreover, across the high-dimensional situations considered here for both internal and bilinear control problems, WeightedPINN yields lower residual-based testing errors than the standard PINN baseline, with the most pronounced improvement usually appearing in the equation residual.

5 Numerical Results: Comparison with other PINN-Type Methods↩︎

5.1 Comparison with Loss-balanced PINN, RAR-D PINN and Causal PINN Methods↩︎

The high-dimensional experiments of Sections 4.2 and 4.3 indicate that WeightedPINN yields smaller residual-based testing errors than the standard PINN baseline for the high-dimensional control problems considered. To probe the operator-level adaptive mechanism of WeightedPINN in a controlled one-dimensional setting, we now design a benchmark in which the residual operators \(\partial_t u\), \(\partial_{xx} u\), \(v(u)\), and \(f\,u\) have markedly different magnitudes by construction. We then compare WeightedPINN against four reference methods: the standard PINN of [12], the loss-balanced PINN of [13], the residual-adaptive (RAR-D) sampler of [23], and the causal PINN of [29]. The self-adaptive PINN of [16], which acts at the level of individual collocation points rather than at the operator level, is treated separately. In contrast to Sections 4.2 and 4.3, where the learned control and state are not displayed due to the high-dimensional architecture, in the present section both the state and the control are presented, and their relative reconstruction errors are computed.

5.1.1 A singularly perturbed bilinear-control heat problem.↩︎

We consider the bilinear-control heat equation in \((0,T)\times\Omega\) with \(T=1\) and \(\Omega=(0,1)\), in which the diffusion coefficient is a small parameter \(\varepsilon\): \[\label{eq:stiff-pde} \begin{align} \partial_t u - \varepsilon\,\partial_{xx} u + v(u) &= f\,u && \text{in }(0,T)\times\Omega,\\ u(t,0) = u(t,1) &= 0 && \text{for }t\in(0,T),\\ u(0,x) &= u_0(x) && \text{in }\Omega,\\ u(T,x) &= z_0(x) && \text{in }\Omega, \end{align}\tag{51}\] with \(v(u) = -u + u^2\), as in the rest of the paper, and \(\varepsilon = 10^{-2}\). We use the manufactured solution \[u^\star(t,x) = \sin(\pi x)\exp\!\bigl(\sin(2\pi t)+t\bigr),\] which yields the reference control \[f^\star(t,x) = 2\pi\cos(2\pi t) + \varepsilon\,\pi^2 - 1 + u^\star(t,x),\] together with \(u_0(x) = \sin(\pi x)\) and \(z_0(x) = e\sin(\pi x)\).

5.1.2 Operator imbalance.↩︎

The motivation for this choice of \(\varepsilon\) is the following. Along the truth, \[\frac{\partial_t u^\star}{u^\star} = 2\pi\cos(2\pi t)+1,\qquad \frac{\varepsilon\,\partial_{xx}u^\star}{u^\star} = -\varepsilon\,\pi^2,\qquad \frac{v(u^\star)}{u^\star} = -1+u^\star,\qquad \frac{f^\star u^\star}{u^\star} = f^\star,\] so that the magnitude of the diffusion contribution \(|\varepsilon\,\partial_{xx} u^\star|\) is of order \(\varepsilon\,\pi^2 \approx 10^{-1}\), while the magnitudes of the time derivative \(|\partial_t u^\star|\) and the control term \(|f^\star u^\star|\) are bounded by \(2\pi+1\approx 7\) times \(|u^\star|\). In other words, the diffusion operator is roughly seventy times smaller than the time derivative and the control term in the residual. This is precisely the type of operator-magnitude imbalance that WeightedPINN is designed to handle through its operator-level adaptive weighting.

5.1.3 Common architecture and training.↩︎

Each method shares the same backbone: the state network \(u_\theta\) and the control network \(f_{\bar\theta}\) are MLPs of width \(48\) and depth \(3\) with \(\tanh\) activation. WeightedPINN has in addition nine weight networks: one per operator component (\(\phi_{\theta^1}\) on \(\partial_t u\), \(\phi_{\bar\theta^1}\) on \(\partial_{xx}u\), \(\phi_{\theta^2}\) on \(v(u)\), and \(\phi_{\theta^3}\) on \(fu\)), separate weights for the two boundary components \(x=0\) and \(x=1\), one weight for the initial constraint, one weight for the terminal constraint, and one weight for the interior observation term, each of width \(20\) and depth \(3\). The weight networks output values in the range \(\phi\in[0.2,5]\), with initialization at \(\phi=1\). The min–max formulation of WeightedPINN is solved with ADAM, learning rate \(10^{-3}\) on the state and control networks and \(2\times 10^{-4}\) on the weight networks, with penalty parameter \(\rho = 8\). All methods are trained for \(4{,}000\) ADAM iterations.

To anchor the otherwise non-unique bilinear control reconstruction, each method receives the same fixed set of \(40\) interior observations of \(u^\star\), with relative observation weight equal to the PDE residual weight. The four reference methods are implemented as in Subsection 1.1: the loss-balanced PINN updates the constraint-loss multiplier every \(50\) iterations according to the gradient-norm rule of [13]; the RAR-D PINN resamples the interior collocation set every \(500\) iterations according to a probability proportional to the squared residual; the causal PINN partitions \((0,T)\) into \(16\) chunks and applies the temporal causal weighting \(w_k = \exp(-\varepsilon_c\sum_{j<k} L_j)\) with \(\varepsilon_c = 1\). In this Section, we use Automatic Differentiation (AD) instead of Finite Differences (FD).

5.1.4 Reported quantities.↩︎

For each method we report the relative \(L^2\) error of the recovered state and recovered control against the manufactured solution, \[\mathrm{u\_err} =\frac{\|u_\theta-u^\star\|_{L^2(Q_T)}}{\|u^\star\|_{L^2(Q_T)}}, \qquad \mathrm{f\_err} =\frac{\|f_{\bar\theta}-f^\star\|_{L^2(Q_T)}}{\|f^\star\|_{L^2(Q_T)}},\] averaged over \(3\) independent random runs. Table 9 reports the mean and the standard deviation across independent runs for each method.

5.1.5 Discussion.↩︎

WeightedPINN attains the smallest mean state error and the smallest mean control error among the methods considered in this benchmark. Quantitatively, on the state error, it improves over the standard PINN by about \(7\%\), over Causal PINN by about \(11\%\), and over Loss-balanced and RAR-D PINN by about \(20\%\); on the control error, the corresponding improvements are about \(4\%\), \(8\%\), \(14\%\), and \(19\%\). The standard deviation of WeightedPINN across independent runs is also among the smallest of the five methods, which indicates that the operator-level adaptive weighting does not come at the cost of training instability.

A natural explanation for the WeightedPINN advantage in this benchmark is the operator-magnitude imbalance discussed above. In a standard PINN, the squared residual \[\bigl(\partial_t u_\theta - \varepsilon\,\partial_{xx} u_\theta + v(u_\theta) - f_{\bar\theta}\,u_\theta\bigr)^2\] is dominated, in magnitude and in gradient signal, by the \(\partial_t u_\theta\) and \(f_{\bar\theta}\,u_\theta\) contributions, whose values are roughly \(70\) times larger than the \(\varepsilon\,\partial_{xx} u_\theta\) contribution. The optimizer of the standard PINN may therefore primarily reduce the time-derivative and control-term contributions, while a sizeable error in \(\partial_{xx}u_\theta\) produces only a small increase in the squared loss and may consequently remain poorly resolved.

The loss-balanced, residual-adaptive-sampling, and causal-weighting methodologies act on the residual as a scalar quantity, or on a scalar quantity localized in time in the causal case, and therefore do not explicitly rescale the individual operator components. WeightedPINN, in contrast, has independent space–time-dependent weights \(\phi_{\theta^1},\phi_{\bar\theta^1},\phi_{\theta^2},\phi_{\theta^3}\) on the four operators, taking values in \([0.2,5]\). The maximization step in the min–max problem can increase the weight associated with the diffusion component in regions where the diffusion residual is poorly resolved, thereby amplifying its contribution to the loss and exposing the operator-magnitude imbalance to the minimization step. This encourages the minimization step to reduce the diffusion residual along with the other operator components.

This mechanism is also consistent with the behavior observed in the high-dimensional experiments of Sections 4.2 and 4.3, where the spatial dimension is \(d=10\) and the Laplacian decomposes into \(d\) directional contributions whose magnitudes may vary across the domain. The fact that WeightedPINN already improves the reconstruction in the present one-dimensional benchmark, where operator imbalance is built in by construction, is consistent with interpreting the high-dimensional residual-error improvements as being connected to the operator-level adaptive mechanism.

5.1.6 Summary.↩︎

The benchmark of this subsection confirms two complementary findings.

1. WeightedPINN performs better than the standard PINN of [12], the loss-balanced PINN of [13], the RAR-D sampler of [23], and the causal PINN of [29] on a one-dimensional bilinear control problem in a singularly perturbed regime, on both the state and control errors.

2. The results are consistent with the interpretation that the operator-level adaptive weights help expose and compensate for operator-magnitude imbalances in the PDE residual. This effect is explicitly built into the present 1D stiff benchmark and appears to be consistent with the residual-error improvements observed in the high-dimensional experiments of Sections 4.2 and 4.3. By contrast, the comparison methods considered here act on the residual as a scalar quantity and therefore do not explicitly rescale the individual operator components.

5.2 Comparison with Control-Specific PINN Methods↩︎

Section 5.1 compared WeightedPINN against generic adaptive PINN methodologies that were not specifically designed for control problems. To complement that benchmark, we now compare WeightedPINN against PINN methodologies that are explicitly built for the optimal control of PDEs:

1. the direct cost-based PINN of Mowlavi and Nabi [59];

2. the Control PINN of Barry-Straume, Sarshar, Popov, and Sandu [60], which couples state and control through optimality conditions in a one-stage architecture;

3. the OCP-PINN of Alla, Bertaglia, and Calzola [40], which derives optimality conditions from the Lagrangian multipliers and represents the state, the adjoint, and the control by a single PINN;

4. the indirect PINN of Zhang, Liu, Alla, Darbon, and Karniadakis [37], which enforces a first-order Pontryagin-type optimality system.

Similar to Section 5.1, the state and control errors are also computed.

5.2.1 Control reconstruction versus PDE-constrained optimal control↩︎

The four methods listed above were designed for PDE-constrained optimal control problems, in which the unknowns \((u,f)\) minimize a prescribed cost objective \[\begin{align} \label{eq:Jocp} J(u,f) = \tfrac{1}{2}\,\|u-u^\star\|^{2}_{L^{2}(Q_T)} +\tfrac{\alpha}{2}\,\|f\|^{2}_{L^{2}(Q_T)}, \qquad \alpha>0, \end{align}\tag{52}\] subject to the state equation and to the prescribed initial condition \(u(0)=u_0\). The Tikhonov term \(\tfrac{\alpha}{2}\|f\|^{2}\) is essential: it convexifies the problem and guarantees that the control is uniquely determined.

The setting of the present paper, originating in [4], is different. We do not minimize a cost objective. We solve the control reconstruction problem: given \(u_0\) and a target terminal state \(z_0\), find a pair \((u,f)\) such that \(u\) satisfies the controlled state equation together with \(u(0)=u_0\) and \(u(T)=z_0\). The terminal target \(z_0\) is a hard constraint, not the minimizer of a tracking term, and there is no Tikhonov regularization of \(f\). The set of admissible pairs \(\mathcal{A}_{u_0}^{z_0}\) (or \(\mathcal{A}_{(u_0,u_1)}^{(z_0,z_1)}\) in the wave case) is in general not a singleton, so the reconstruction problem is intrinsically non-unique.

Comparing the two settings is therefore non-trivial. Methods \((1)\)–\((4)\) above must first be adapted to the reconstruction setting before they can be evaluated against WeightedPINN on the same benchmark. In what follows we describe each method briefly and indicate how it is adapted. In this Section, we use Automatic Differentiation (AD) instead of Finite Differences (FD).

5.2.2 Adapted formulations↩︎

We retain the manufactured-solution benchmark of Section 5.1: the bilinear-control heat equation in \((0,T)\times(0,1)\) with \(T=1\), \(v(u)=-u+u^2\), and \[u^\star(t,x)=\sin(\pi x)\exp\!\bigl(\sin(2\pi t)+t\bigr), \quad f^\star(t,x)=2\pi\cos(2\pi t)+\pi^2+u^\star(t,x).\] The initial and terminal data are \(u_0(x)=\sin(\pi x)\) and \(z_0(x)=e\sin(\pi x)\), and a fixed set of \(N_{\mathrm{obs}}=30\) interior observations of \(u^\star\) is provided to all methods, as in Subsection 5.1.

For methods \((1)\)–\((4)\), the cost objective is built from the interior observations and from a Tikhonov regularization of the control. We use the same tracking-plus-control regularization form as in the original references, in which the tracking term is defined on the observation set: \[\begin{align} \label{eq:Jdata} \widetilde{J}(u,f) = \tfrac{1}{2N_{\mathrm{obs}}} \sum_{i=1}^{N_{\mathrm{obs}}}\bigl(u(t_i,x_i)-u^\star(t_i,x_i)\bigr)^2 +\tfrac{\alpha}{2}\,\|f\|^{2}_{L^{2}(Q_T)}. \end{align}\tag{53}\] The Tikhonov coefficient is fixed to \(\alpha=10^{-3}\). The boundary, initial, and terminal constraints are enforced by penalty terms identical to those used for all other methods in Section 5.1.

For simplicity, in the adjoint-based adapted methods below we enforce the smooth bulk adjoint equation and the stationarity condition in a penalty form, while the observation mismatch and terminal target are included through the same soft penalty terms used by all methods. Thus these methods should be interpreted as reconstruction-oriented adaptations of the corresponding optimal-control PINNs, rather than as exact implementations of their original first-order optimality systems.

5.2.2.1 (1) Direct cost-based PINN [59].

The state \(u_\theta\) and the control \(f_{\bar\theta}\) are represented by feedforward networks. The training loss is \[\begin{align} \mathcal{L}_{\mathrm{MN}} = w_{\mathrm{PDE}}\,\|\mathcal{R}_Q^{3}\|^2_{Q} +\|\mathcal{R}_\Gamma^{3}\|^2_{\Gamma} +\widetilde{J}(u_\theta,f_{\bar\theta}), \end{align}\] where \(\mathcal{R}_Q^{3}\) and \(\mathcal{R}_\Gamma^{3}\) are the standard PINN residuals defined in Section [Sec:bilinearsettings] and the scalar weight \(w_{\mathrm{PDE}}\) on the PDE residual is selected by a short two-step line search over \(\{0.5, 1, 2, 5\}\) at the beginning of training, in the spirit of the line-search strategy of [59]. No adjoint variable is introduced.

5.2.2.2 (2) Adapted Control PINN [60].

This formulation introduces a third neural network \(p_{\tilde{\theta}}\) for the adjoint variable and trains the state, control, and adjoint networks jointly. Since the original Control PINN was developed for PDE-constrained optimal control, we use here a penalty-based adaptation to the present reconstruction benchmark. The loss is \[\mathcal{L}_{\mathrm{CP}} = \|\mathcal{R}_Q^{3}\|^2_{Q} +\|\mathcal{R}_\Gamma^{3}\|^2_{\Gamma} +\widetilde{J}(u_\theta,f_{\bar\theta}) +\|\mathcal{R}_{\mathrm{adj}}\|^2_{Q} +\|\mathcal{R}_{\mathrm{opt}}\|^2_{Q} +\|\mathcal{R}_{\mathrm{adj},\Gamma}\|^2_{\Gamma}.\] Here \[\mathcal{R}_{\mathrm{adj}} = -\partial_t p_{\tilde{\theta}} -\partial_{xx}p_{\tilde{\theta}} + \bigl(v'(u_\theta)-f_{\bar\theta}\bigr)p_{\tilde{\theta}}\] is the bulk adjoint residual associated with the bilinear heat operator, and \[\mathcal{R}_{\mathrm{opt}} = \alpha f_{\bar\theta} - u_\theta p_{\tilde{\theta}}\] is the stationarity residual, up to the sign convention used in the Lagrangian. The tracking mismatch is included through the cost term \(\widetilde{J}(u_\theta,f_{\bar\theta})\). Thus this should be understood as a penalty-based adaptation of Control PINN to the present reconstruction benchmark, rather than as an exact enforcement of the full distributional first-order optimality system.

5.2.2.3 (3) OCP-PINN [40].

As in [40], a single neural network with three output channels represents the triple \((u_\theta,p_\theta,f_\theta)\). The loss aggregates the state PDE residual, the adjoint PDE residual, the stationarity residual, and the observation-tracking term, together with boundary, initial, terminal, and adjoint-terminal penalties. This formulation differs from \((2)\) primarily in its shared-backbone architecture, which forces \((u,p,f)\) to share representational features and reduces the number of trainable parameters.

5.2.2.4 (4) Indirect PINN [37].

Three separate networks represent the state, the adjoint, and the control. They are trained to satisfy a penalty-based adaptation of the first-order optimality system in the residual sense: the forward state equation, the smooth bulk adjoint equation, and the stationarity condition. As suggested in [37], the stationarity residual is given a stronger penalty weight than in \((2)\). The observation mismatch is included through the tracking term in \(\widetilde{J}\), as in the other cost-based methods.

5.2.3 Setup↩︎

All methods use the same backbone: state, control, and adjoint networks (when present) are MLPs of width \(32\) and depth \(3\) with \(\tanh\) activation. WeightedPINN additionally uses seven weight networks of width \(20\), depth \(3\), with output \(\phi(t,x)=1+\delta\tanh(\cdot)\), \(\delta=0.5\), covering the time-derivative, spatial-derivative, nonlinear-response, and control terms of the PDE residual, and the boundary, initial, and terminal constraints. All methods are trained for \(6000\) ADAM iterations with learning rate \(10^{-3}\) on the state, control, and adjoint networks; WeightedPINN uses an additional learning rate of \(10^{-4}\) for its weight networks, and penalty \(\rho=20\). Each iteration uses \(800\) interior, \(160\) boundary, \(160\) initial, \(160\) terminal collocation points, and the same fixed set of \(N_{\mathrm{obs}}=30\) interior observations of \(u^\star\). The Tikhonov coefficient is \(\alpha=10^{-3}\) for the cost-based methods (1)–(4); the PDE-residual weight \(w_{\mathrm{PDE}}\) in method (1) is chosen by the two-step line search of Subsection 5.2.2.

For each method we report the relative \(L^2\) errors of the recovered state and the recovered control against the manufactured truth, \[\mathrm{u\_err} = \frac{\|u_\theta-u^\star\|_{L^2(Q_T)}}{\|u^\star\|_{L^2(Q_T)}}, \qquad \mathrm{f\_err} = \frac{\|f_{\bar\theta}-f^\star\|_{L^2(Q_T)}}{\|f^\star\|_{L^2(Q_T)}},\] averaged over \(3\) independent random runs. These two quantities measure the quality of the reconstruction of \((u^\star,f^\star)\), which is precisely the objective of the control reconstruction problem considered in the present paper.

5.2.4 Results and Discussions↩︎

Table 10 reports the relative \(L^{2}\) errors of the recovered state and control against the manufactured truth, averaged over \(3\) independent runs.

The results in Table 10 expose a clear separation between the two families of methods on the quantities that define the control reconstruction problem.

(i) WeightedPINN attains the smallest reconstruction errors. WeightedPINN achieves the smallest mean \(\mathrm{u\_err}\) (\(3.10\times 10^{-1}\)) and the smallest mean \(\mathrm{f\_err}\) (\(2.91\times 10^{-1}\)) of all six methods. The standard PINN baseline is the second best on both metrics, with mean errors essentially comparable to those of WeightedPINN. The four control-aware methods (1)–(4) all yield substantially larger mean errors, with \(\mathrm{u\_err}\) in \([0.49,\,0.54]\) and \(\mathrm{f\_err}\) in \([0.35,\,0.39]\). On the metrics that correspond directly to the control reconstruction problem, the PINN-style methods perform better than the control-aware methods, and WeightedPINN slightly improves upon the standard PINN baseline.

(ii) The Tikhonov bias of the control-aware methods. The systematic gap between the two families on the reconstruction errors admits a structural explanation. The four control-aware methods all incorporate a Tikhonov term \(\tfrac{\alpha}{2}\|f\|^{2}\) in their training objective, see 53 . This term shrinks the recovered control toward zero, which prevents \(f_{\bar\theta}\) from reaching the magnitude of the true control \(f^\star\). The PINN-style methods do not include this term, and they are therefore not biased toward small controls. The asymmetric pattern in Table 10 is exactly what one would expect from a well-regularized but biased estimator of \(f^\star\): the control-aware methods exhibit very low across-run variance (the standard deviations of methods (1)–(4) are an order of magnitude smaller than for the PINN-style methods), but the mean errors are systematically larger because the recovered control underestimates \(f^\star\).

(iii) WeightedPINN reduces variance with respect to the standard PINN baseline. Compared to the standard PINN, WeightedPINN reduces the across-run standard deviation by approximately \(40\%\) on the state error (\(6.6\times 10^{-2}\) versus \(1.1\times 10^{-1}\)) and approximately \(30\%\) on the control error (\(5.5\times 10^{-2}\) versus \(7.7\times 10^{-2}\)), while preserving the smallest mean errors. This is consistent with the variance-reduction finding reported in Section 5.1 and indicates that the operator-decomposed adaptive weighting of WeightedPINN also acts as a stabilizer of the training dynamics in the present comparison.

(iv) The two families occupy complementary niches. The pattern in Table 10 can be summarized as follows. Cost-based control-aware methods are the right choice when the goal is to identify an optimally regularized control that trades faithfulness to a target trajectory against a control-energy penalty; the Tikhonov term then convexifies the problem and yields low across-run variance. PINN-style methods, and in particular WeightedPINN, are the right choice when the goal is to reconstruct a specific admissible pair \((u,f)\in\mathcal{A}_{u_0}^{z_0}\), since they do not introduce the Tikhonov bias on \(f\) and therefore allow the recovered control to reach the magnitude of \(f^\star\). The two families are complementary, and the choice between them is determined by whether the application requires an optimally regularized control or a faithful reconstruction of an unbiased control.

5.2.5 Summary.↩︎

The comparison of this section confirms that WeightedPINN occupies a methodological niche that is distinct from the existing control-aware PINN methodologies. The cost-based methods of [37], [40], [59], [60] are designed for PDE-constrained optimal control problems in which a cost objective is to be minimized and the control is regularized, which biases the recovered control toward smaller magnitudes. WeightedPINN is designed for the control reconstruction problem, in which the terminal state is a hard target and the control is identified directly through the structure of the controlled state equation. On the one-dimensional benchmark of this section, WeightedPINN attains the smallest mean relative \(L^2\) errors against the true \((u^\star,f^\star)\) of all six methods. The two families of methods are therefore complementary rather than competing, and the present comparison clarifies the regime in which each is preferable.

5.3 Numerical Results: WeightedPINN versus Self-Adaptive PINN on a Variable-Coefficient Wave Control Problem↩︎

This section presents a comparison in a setting where the structural distinction between WeightedPINN and Self-Adaptive PINN becomes particularly transparent: a variable-coefficient wave equation in which the leading-order spatial operator itself depends on space.

5.3.1 Why a variable coefficient is a natural test for operator-level adaptivity↩︎

The self-adaptive PINN of [16] replaces the PDE residual \(\sum_j |R[u_\theta](t_j,x_j)|^2\) by \(\sum_j \lambda_j^2\,|R[u_\theta](t_j,x_j)|^2\), where the per-collocation-point weights \(\lambda_j\) are trainable. The adaptive mechanism is therefore localized at the level of collocation points: it can identify regions in which the residual is hard to drive to zero and emphasize them. However, at any given point \((t_j,x_j)\), the weight \(\lambda_j\) acts on the entire residual as a single scalar. It cannot distinguish, at that point, between an error coming from the time-derivative term and an error coming from the spatial elliptic term.

WeightedPINN, in contrast, introduces space–time-dependent weights at the level of each operator component of the residual; see Subsection 1.2 and Subsection 1.3 of the present paper. At any given point, the time-derivative, the diffusion operator, the nonlinear response, and the control term carry independent adaptive weights.

In the constant-coefficient setting considered in Sections 5.1–4, the operator components of the residual have spatially uniform strength: \(\partial_{x_i x_i} u_\theta\) appears with a constant multiplier. The benefit of separating the residual into operator components is therefore most apparent when the coefficients of the PDE vary in space, so that different operator components dominate the residual in different regions of the domain. The variable-coefficient wave equation \[\label{eq:varcoef-wave} \partial_{tt} u \;-\;\nabla\cdot\bigl(a(x)\,\nabla u\bigr) \;+\;v(u) \;=\;f\tag{54}\] is the natural test case: in regions where \(a(x)\) is large the spatial diffusion term dominates, in regions where \(a(x)\) is small the time-derivative, nonlinear, and forcing terms relatively dominate, and the cross term \(\nabla a\cdot\nabla u\) only contributes significantly in regions where \(\nabla a\) is significant. A per-point scalar weight can emphasize difficult regions, but it does not explicitly distinguish which operator component is responsible for the local residual imbalance. A per-operator spatially varying weight, on the other hand, is designed to separate these componentwise effects.

5.3.2 Setting↩︎

We consider the controlled variable-coefficient wave equation 54 in \((0,T)\times\Omega\) with \(\Omega=(0,1)^{2}\), \(T=1\), and \(v(u)=-u+u^{2}\). Homogeneous Dirichlet conditions are imposed on the lateral boundary, and initial and terminal data \((u_0,u_1)\) and \((z_0,z_1)\) are prescribed in \(H^{1}(\Omega)\times L^{2}(\Omega)\) as in Section [Sec:internalsettings]. The control \(f\) is internal and supported on all of \(\Omega\).

The variable coefficient is taken as \[\begin{align} \label{eq:varcoef-a} a(x_1,x_2) = 1 + 0.7\,\cos(2\pi x_1)\cos(2\pi x_2), \end{align}\tag{55}\] so that \(a\in[0.3,1.7]\), i.e., the leading-order spatial operator has a \(5.7\)-fold variation across the domain. The cross term \(\nabla a\cdot\nabla u\) has the structure \[\begin{align} \nabla a\cdot\nabla u = -1.4\pi\sin(2\pi x_1)\cos(2\pi x_2)\,\partial_{x_1}u - 1.4\pi\cos(2\pi x_1)\sin(2\pi x_2)\,\partial_{x_2}u, \end{align}\] which is significant in the transition regions of the domain. The four operator components \(\partial_{tt}u\), \(a\,\partial_{x_i x_i}u\), \(\partial_{x_i}a\cdot\partial_{x_i}u\), and \(v(u)\) therefore all carry spatial structure of distinct character, and their relative magnitudes vary significantly across the domain.

We benchmark on a manufactured-solution setup. The reference state is \[\begin{align} u^{\star}(t,x_1,x_2) = \sin(\pi x_1)\sin(\pi x_2)\,(1+2t-t^{2}), \end{align}\] which satisfies the Dirichlet boundary condition. The corresponding reference control is computed from 54 : \[\begin{align} f^{\star}(t,x_1,x_2) = \partial_{tt}u^{\star} - \nabla\cdot\bigl(a(x)\,\nabla u^{\star}\bigr) + v(u^{\star}). \end{align}\] The initial and terminal data \(u_0(x)=u^{\star}(0,x)\), \(u_1(x)=\partial_t u^{\star}(0,x)\), \(z_0(x)=u^{\star}(T,x)\), \(z_1(x)=\partial_t u^{\star}(T,x)\) are read off from \(u^{\star}\), and a fixed set of \(N_{\mathrm{obs}}=60\) interior observations of \(u^{\star}\) is provided to both methods to anchor the reconstruction. We also employ AD instead of FD. We note that the theoretical results established in Section 3 extend naturally to the case of variable coefficients.

5.3.3 Methods↩︎

5.3.3.1 Self-Adaptive PINN [16].

The state \(u_\theta\) and the control \(f_{\bar\theta}\) are represented by feedforward networks. The training loss is a sum of per-collocation-point weighted residuals, \[\begin{align} \mathcal{L}_{\mathrm{SA}} &= \tfrac{1}{N_{\mathrm{int}}} \sum_{j=1}^{N_{\mathrm{int}}} \lambda_j^{2}\, \bigl|R^{\mathrm{wave}}[u_\theta,f_{\bar\theta}](t_j,x_j)\bigr|^{2} \;+\; \tfrac{1}{N_{\mathrm{bd}}} \sum_{j=1}^{N_{\mathrm{bd}}} \lambda^{\mathrm{bd}}_j{}^{2}\, \bigl|u_\theta(\hat{t}_j,\hat{x}_j)\bigr|^{2} \\ &\quad +\;\text{(initial, terminal, and observation terms, with their own per-point trainable weights)}, \end{align}\] where \(R^{\mathrm{wave}}[u_\theta,f_{\bar\theta}] =\partial_{tt}u_\theta-\nabla\cdot(a\nabla u_\theta)+v(u_\theta) -f_{\bar\theta}\), and the collocation points are fixed at the beginning of training. The state, control, and trainable weights are updated jointly: the network parameters \((\theta,\bar\theta)\) by gradient descent, the per-point weights by gradient ascent, as in [16]. At each collocation point the weight acts on the entire residual as a single scalar.

5.3.3.2 WeightedPINN

The state \(u_\theta\) and the control \(f_{\bar\theta}\) are represented by feedforward networks of the same architecture as for SA-PINN. The PDE residual is decomposed at the operator level: \[\begin{align} \label{eq:wpinn-varcoef-residual} \bar{R}_Q^{\mathrm{vc}} &= \phi_{\theta^1}\,\partial_{tt}u_\theta -\sum_{i=1}^{d}\phi_{\bar{\theta}^{i}}\,a\,\partial_{x_i x_i}u_\theta -\sum_{i=1}^{d}\phi_{\tilde{\theta}^{i}}\,(\partial_{x_i}a)\,\partial_{x_i}u_\theta\\ &\quad +\phi_{\theta^2}\,v(u_\theta) -\phi_{\theta^3}\,f_{\bar\theta},\nonumber \end{align}\tag{56}\] where \(d=2\) here and each weight network \(\phi(t,x)\) is a smooth function of space–time. Together with weights on the boundary condition and on the four initial/terminal constraints \(u(0), \partial_t u(0), u(T), \partial_t u(T)\), this gives \(12\) weight networks in total. The training problem is the min–max formulation of Section 2, in which the state and control networks minimize the weighted residual and constraints while the weight networks maximize them, subject to a quadratic penalty that keeps each \(\phi\) close to \(1\). The decomposition 56 is the natural extension of the constant-coefficient WeightedPINN formulation of Section [Sec:internalsettings] to the variable-coefficient operator \(\nabla\cdot(a(x)\nabla u)\). In contrast to Sections 4.2 and 4.3, in the present section both the state and the control are presented, and their relative reconstruction errors are computed.

5.3.4 Setup↩︎

For both methods, the state network \(u_\theta\) and the control network \(f_{\bar\theta}\) are MLPs of width \(48\) and depth \(4\) with \(\tanh\) activation. WeightedPINN additionally uses \(12\) weight networks of width \(20\) and depth \(3\) with output \(\phi(t,x)=1+\delta\tanh(\cdot)\), \(\delta=0.5\). Both methods are trained for \(2500\) ADAM iterations with learning rate \(10^{-3}\) on the state and control networks; WeightedPINN uses an additional learning rate of \(5\times 10^{-4}\) on its weight networks, and penalty \(\rho=3\). Each iteration uses \(600\) interior, \(200\) boundary, \(150\) initial, \(150\) terminal collocation points, and the same fixed set of \(N_{\mathrm{obs}}=60\) interior observations of \(u^{\star}\). Both methods are evaluated on a dense \(30\times 30\times 30\) space–time grid. Derivatives of the network outputs are computed by automatic differentiation. We report the relative \(L^{2}\) errors of the recovered state and the recovered control against the manufactured truth, \[\mathrm{u\_err} = \frac{\|u_\theta-u^{\star}\|_{L^{2}(Q_T)}}{\|u^{\star}\|_{L^{2}(Q_T)}}, \qquad \mathrm{f\_err} = \frac{\|f_{\bar\theta}-f^{\star}\|_{L^{2}(Q_T)}}{\|f^{\star}\|_{L^{2}(Q_T)}},\] averaged over \(3\) independent random runs, together with the across-run standard deviation.

5.3.5 State and control errors↩︎

Table 11 reports the relative \(L^{2}\) errors of the recovered state and control against the manufactured truth, averaged over \(3\) independent runs; Table 12 gives the per-independent run values for the same comparison.

5.3.6 Discussion↩︎

The comparison of WeightedPINN and the Self-Adaptive PINN in Tables 11–12 yields three observations.

(i) WeightedPINN attains substantially smaller mean errors than the Self-Adaptive PINN on both metrics. On the recovered state, WeightedPINN reaches a mean error of \(2.51\times 10^{-1}\), versus \(4.55\times 10^{-1}\) for the Self-Adaptive PINN, a \(45\%\) reduction. On the recovered control, WeightedPINN reaches a mean error of \(6.61\times 10^{-1}\), versus \(8.52\times 10^{-1}\) for the Self-Adaptive PINN, a \(22\%\) reduction. Moreover, as shown in Table 12, this advantage is not an artifact of averaging: WeightedPINN attains a strictly smaller \(\mathrm{u\_err}\) and a strictly smaller \(\mathrm{f\_err}\) than the Self-Adaptive PINN on every independent run.

(ii) WeightedPINN is substantially more stable across random initializations. The across-run standard deviation of the Self-Adaptive PINN is large on both metrics (\(1.7\times 10^{-1}\) on \(\mathrm{u\_err}\) and \(1.6\times 10^{-1}\) on \(\mathrm{f\_err}\)), driven primarily by independent run \(2\) in Table 12 where the Self-Adaptive PINN essentially fails to fit the manufactured target. WeightedPINN reduces the across-run standard deviation by \(71\%\) on \(\mathrm{u\_err}\) and by \(69\%\) on \(\mathrm{f\_err}\). The structured operator-level adaptive weighting of WeightedPINN is therefore much more robust to random initialization than the unstructured per-point adaptive weighting of the Self-Adaptive PINN. This instability of the Self-Adaptive PINN is consistent with the possibility that per-point trainable weights may amplify residual noise in regions where the network is initially poorly resolved, thereby making the training dynamics more sensitive to initialization.

(iii) The performance gap reflects the structural difference between point-level and operator-level adaptivity. Both methods access the same training data, the same manufactured target, and the same network architectures for the state and the control. They differ only in the form of their adaptive weights:

• Self-Adaptive PINN: one trainable scalar \(\lambda_j\) per fixed collocation point \((t_j,x_j)\), multiplying the entire PDE residual at that point.

• WeightedPINN: one trainable space–time function \(\phi(t,x)\) per operator component \(\partial_{tt}u_\theta\), \(a\partial_{x_i x_i}u_\theta\), \(\partial_{x_i}a\,\partial_{x_i}u_\theta\), \(v(u_\theta)\), \(f_{\bar\theta}\), multiplying that specific operator component throughout the domain.

The two adaptive families therefore differ in where their flexibility is spent. Per-point adaptivity is more expressive at the point level but cannot distinguish which operator component is responsible for the local residual error. Per-operator spatially varying adaptivity is less expressive at any single point but is consistent with the fact that, in the variable-coefficient setting, the local balance of operator components in the PDE residual is itself spatially varying: in regions where \(a(x)\) is large the diffusion term \(a\,\partial_{x_i x_i}u\) dominates, in regions where \(\nabla a\) is significant the cross term \(\partial_{x_i}a\,\partial_{x_i}u\) becomes non-negligible, and where \(a\) is small the time-derivative and forcing terms relatively dominate. WeightedPINN’s per-operator weights are designed to match this spatial operator structure, whereas SA-PINN’s per-point weights do not explicitly distinguish among the operator components.

5.3.7 Summary.↩︎

The two-dimensional variable-coefficient wave benchmark of this section compares two adaptive PINN methodologies on the same control reconstruction problem: the Self-Adaptive PINN of [16], with per-collocation-point scalar weights, and the WeightedPINN of the present paper, with per-operator space–time weights. On this benchmark, WeightedPINN attains smaller relative errors than the Self-Adaptive PINN on every independent run, both for the recovered state and for the recovered control, and the across-run standard deviation is reduced by approximately \(70\%\) on both metrics. The gap is consistent with the structural argument of Subsection 5.3.1: in the variable-coefficient setting, the local balance of operator components in the PDE residual varies spatially, and matching this structure favors adaptive weights that are operator-decomposed rather than only point-decomposed.

References↩︎

---

1. M.-N. Phung and M.-B. T are funded in part by NSF CAREER DMS-2303146, and NSF Grants DMS-2204795, DMS-2305523, DMS-2306379.↩︎