Penalty Adversarial Network (PAN): A neural network-based method to solve PDE-constrained optimal control problems
Shilin Ma Yukun Yue

 arabic

Non-Relativistic Holography from AdS_5/CFT_4
    Juan Miguel Nieto García
    September 9, 2024
============================================



§ ABSTRACT 
In this work, we introduce a novel strategy for tackling constrained optimization problems through a modified penalty method. Conventional penalty methods convert constrained problems into unconstrained ones by incorporating constraints into the loss function via a penalty term. However, selecting an optimal penalty parameter remains challenging; an improper choice, whether excessively high or low, can significantly impede the discovery of the true solution. This challenge is particularly evident when training neural networks for constrained optimization, where tuning parameters can become an extensive and laborious task. To overcome these issues, we propose an adversarial approach that redefines the conventional penalty method by simultaneously considering two competing penalty problems—a technique we term the penalty adversarial problem. Within linear settings, our method not only ensures the fulfillment of constraints but also guarantees solvability, leading to more precise solutions compared to traditional approaches. We further reveal that our method effectively performs an automatic adjustment of penalty parameters by leveraging the relationship between the objective and loss functions, thereby obviating the need for manual parameter tuning. Additionally, we extend this adversarial framework to develop a neural network-based solution for optimal control problems governed by linear or nonlinear partial differential equations. We demonstrate the efficacy of this innovative approach through a series of numerical examples. 












§ INTRODUCTION


Optimal control problems are fundamental in various scientific and engineering disciplines. These problems involve finding a control state y that determines the desired state u through governing physical constraints, aiming to minimize or maximize a given performance criterion, typically expressed as an objective functional <cit.>. In many practical scenarios, the system dynamics are governed by partial differential equations (PDEs), leading to PDE-constrained optimal control problems <cit.>. These problems have gained significant attention due to their critical applications in fields such as aerospace engineering <cit.>, environmental marine sciences <cit.>, medical treatment planning for radiation therapy <cit.>, heat transfer <cit.>, fluid dynamics <cit.>, liquid crystals <cit.>, and wave propagation <cit.>. Achieving optimal performance while adhering to physical laws and constraints is crucial in these applications.


Mathematically, we can formulate the problem as 

    min_u ∈U, y ∈Y J(u,y),   subject to  F(u,y) = 0.

Here, J(u,y) represents the performance criterion to be minimized, often referred to as the objective functional <cit.>. The term F(u,y) contains the constraints that u and y must satisfy, including the differential operators in the form of PDEs and the boundary or initial conditions for the PDEs. We denote U and Y as the appropriate spaces in which u and y reside, respectively.

The highly nonlinear nature and multi-scale structure of PDE-constrained optimal control problems <cit.> necessitate using complex numerical methods. Over the years, various approaches have been developed to create robust and accurate numerical algorithms and tools to solve these problems. Among the prevalent methods, adjoint-based techniques are particularly notable for their effectiveness in gradient computation, which is crucial for iterative optimization algorithms <cit.>. These methods are often combined with traditional numerical techniques, such as the finite difference method or finite element methods, to handle spatial and temporal discretizations, allowing for the management of complex geometries and boundary conditions <cit.>.


Recently, with the rapid development of neural networks, machine learning-based numerical methods have been extensively developed to solve PDEs <cit.>. As an important extension of this work, considerable research efforts have focused on applying these methods to solve PDE-constrained optimal control problems. Several successful examples have emerged in this area <cit.>. Among them, the most prevalent approaches can be classified into three main categories:


    
  * Training surrogate models to obtain solvers for PDEs, then using the trained solver to map inputs to solutions of the PDE, enforcing the PDE constraint while minimizing an objective cost functional <cit.>;
    
  * Using the Lagrangian approach to reformulate the constrained optimization problem and solve systems associated with the Karush–Kuhn–Tucker (KKT) conditions <cit.> via classical neural network-based methods <cit.>;
    
  * Adding the cost functional to the standard loss induced by the PDEs and minimizing the total loss simultaneously <cit.>.


A common challenge these approaches encounter is enforcing the constraints during optimization. The most direct approach is to treat the constraints as penalty terms and add them to the original loss function minimized by the neural network, transforming a constrained problem into an unconstrained one. This transformation makes it easier to apply various iterative optimization tools to solve the problem <cit.>, as implemented in the third approach above. Though not directly employing this method when enforcing PDE constraint, the other two approaches implicitly address the same challenge while solving the PDE by ensuring that initial or boundary conditions are satisfied <cit.>. For simplicity, we will explain our main idea concerning the last approach as an example. At the same time, the same analysis can be applied to the other two approaches when penalty terms are introduced.

Specifically, instead of solving problem (<ref>), we consider a penalty problem which can be formulated as:

    min_u∈U,y∈Y(u,y),

where  is defined as

    (u,y) = J(u,y) + λ/2 F(u,y)^2.

Here, λ > 0 is a tunable penalty parameter, and · denotes the standard L^2 norm. It is important to note that problem (<ref>) is not equivalent to the constrained problem (<ref>). However, as λ approaches infinity, the solution of this unconstrained problem converges to the solution of the constrained one <cit.> (We will provide more details on this in the linear case in Section <ref>). Conversely, when the penalty parameter becomes too large, the problem can become ill-posed and difficult to solve <cit.>. If a neural network is used to solve this optimization problem, achieving convergence during training can be challenging. On the other hand, if the penalty parameter is too small, the PDE constraints will not be adequately satisfied, resulting in a solution that deviates significantly from the desired solution of the original constrained problem. 

To address this challenge, <cit.> proposes a two-step line-search approach to determine the optimal penalty parameter that minimizes the cost function while ensuring the PDE constraints are satisfied within an acceptable tolerance. In contrast, <cit.> suggests dynamically adjusting the penalty parameters based on a problem-dependent update rule <cit.>. These approaches commonly face issues such as the extensive effort required for parameter tuning and the absence of general guidelines for selecting an appropriate penalty parameter for various problems. This complicates the implementation of existing methods for solving optimal control-related problems.

In this paper, we propose a novel framework based on penalty methods and present it using a neural network structure that does not require an artificial selection process for the varying value of the penalty parameter. Instead, we construct two neural networks to minimize (u,y) with different fixed values of λ and train them simultaneously. One network competes with the other during training to ensure ease of training and convergence while maintaining the constraints within an acceptable tolerance. This concept is inspired by the well-known generative adversarial network (GAN) <cit.>, which has extensive applications, for example, in image translation <cit.>, video generation <cit.>, and speech synthesis <cit.>. The adversarial network structure incorporates a generator and a discriminator, with the discriminator providing feedback to the generator to help it produce more realistic information to deceive the discriminator. This idea has also benefited research on numerical methods based on adversarial network structures for solving PDEs. We refer readers to <cit.> and the references therein.

In particular, we create a solver network (corresponding to the generator in GAN) and a discriminator network, and we choose two real numbers λ_1, λ_2 > 0 with λ_1 being much larger than λ_2, and λ_2 being relatively small. The discriminator network is set to minimize the objective functional (u,y) as defined in (<ref>). In practice, the smallness of λ_2 ensures the convergence of the discriminator network as long as the optimal control problem (<ref>) is well-posed. Conversely, the solver network aims to minimize the sum of (u,y) and an additional cost based on feedback from the discriminator network. The exact functional form of this extra cost will be detailed in Section <ref>. 

With feedback from the discriminator network, the solver network can automatically adjust its focus during training without manual tuning. If the solver network emphasizes reducing the objective functional at the expense of not maintaining the PDE constraint, the weight of the PDE constraint will increase accordingly. Conversely, if the solver network focuses too much on satisfying the PDE constraint but fails to reduce the objective cost functional, it will adjust its weights on the objective functionals to correct this imbalance.


To this end, the proposed framework utilizes the strengths of both traditional penalty methods and the adversarial training paradigm, offering a novel solution to the challenges inherent in PDE-constrained optimal control problems. Numerical examples, presented in Section <ref>, demonstrate that our approach ensures robust convergence and effective enforcement of PDE constraints without requiring extensive parameter tuning. This dual-network strategy not only simplifies the training process but also enhances the overall stability and performance of the optimization.


The rest of this paper is structured as follows: In Section <ref>, we describe various PDE-constrained optimal control problems in general forms that will be the focus of this paper. Section <ref> serves as a motivation for our methodology, where a linear problem is discussed to illustrate the effectiveness of the proposed penalty adversarial framework in solving problems with penalty formulations. We demonstrate that in the linear case, the solution to the adversarial problem better adheres to the constraints than simply solving the problem with a small penalty parameter under certain conditions. Following this analysis, Section <ref> provides a detailed construction of a neural network-based method utilizing this concept, culminating in a practical algorithm. Finally, in Section <ref>, we conduct numerical experiments on both linear and nonlinear problems in 1D and 2D to validate the effectiveness of the proposed approach.







§ PROBLEM SETUP


This section formally outlines the various types of optimal control problems, including distributed, boundary, and initial value control problems. We consider an open bounded physical domain Ω⊂ℝ^d, where d is a positive integer denoting the spatial dimension, and a time span [0, T]. The problem is governed by the following system:


    
    ℒ[u(x,t),y_f(x,t)]    = 0,   ∀ x ∈Ω,   t ∈ [0,T],  
    ℬ[u(x,t),y_b(x,t)]    = 0,   ∀ x ∈∂Ω,   t ∈ [0,T], 
    ℐ[u(x,0),y_i(x)]    = 0,   ∀ x ∈Ω.


Here, x and t denote the spatial and temporal variables, respectively. ℒ is an operator involving differentials that represents the PDE to be satisfied by u, ℬ denotes the boundary condition, and ℐ represents the initial condition. Common choices for boundary conditions include Dirichlet, Neumann, or Robin conditions, and our framework imposes no specific restrictions on these choices. The spaces U and Y, in which u and y reside, are selected to ensure the well-posedness of the PDE problem. For instance, if ℒ(u,y) = Δ u - y, corresponding to a standard second-order elliptic equation with no initial condition and a homogeneous Dirichlet boundary condition, appropriate spaces to consider are U = H^1_0(Ω) and Y = L^2(Ω) <cit.>.


If we set y = (y_f, y_b, y_i) and define the constraint F(u,y) = (ℒ(u,y), ℬ(u,y), ℐ(u,y)), we recover the constraint given in (<ref>). Specifically, the entire system (<ref>) or any individual equation within it can be viewed as a specific example of the general form of constraints F(u,y) = 0.

By setting different components of y in (<ref>) to be tunable, we obtain various types of control problems. For example, if we consider y_b and y_i to be fixed and take y_f as a tunable control, we obtain a distributed control problem, initially introduced in <cit.>. The objective function to be minimized in this case is:

    J_d(u,y) = 1/2u - û^2 + ρ/2y_f^2,

where û denotes the desired state that we aim for our solution of the PDE system to approximate by tuning the value of y_f. The second term in this functional is a Tikhonov regularization term <cit.>. Generally, the problem can be ill-posed without such a regularization term, and the Tikhonov regularization parameter ρ value is typically chosen in advance <cit.>. 



If we consider y_b to be tunable and y_f and y_i to be fixed, then we obtain a boundary control problem, which minimizes the following objective function:

    J_b(u,y) = 1/2u - û^2 + ρ/2y_b^2,

with the same notation for û and ρ. Similarly, one can define an initial value optimal control problem. 


Additionally, to clarify our terminology: from now on, we refer to functionals like J(u,y) as objective functionals, as they represent the goal that we aim to minimize with (u,y) found by our algorithms. On the other hand, we refer to functionals like (u,y) as cost functionals or loss functionals.



It is important to note that the optimal control problems listed here only encompass some possible applications of our proposed method. General PDE-constrained optimization problems can be adapted to fit within our framework. Our methodology can be viewed as a variant of the penalty method. As long as a problem can be resolved or approximated using the penalty method, it is feasible to implement our approach. This study will focus on distributed and boundary optimal control problems because they are typically classic and representative examples.







§ DISCRETIZED PROBLEM


This section is devoted to discussing our motivation for setting up our method. In PDE-constrained optimization problems, there is ongoing debate on whether to use the discretize-then-optimize or optimize-then-discretize strategy <cit.>. Since our method is based on a neural network, the autodifferentiation technique <cit.> naturally leads to a discretized system to solve. We adopt the discretize-then-optimize approach, beginning with an analysis of a discretized problem.

We consider the following discretized constrained optimization problem: find u ∈ℝ^n and y ∈ℝ^m to minimize

    J(u,y) = 1/2Au - b^2 + ρ/2y^2,

subject to

    Ku = y,

where 0 < m ≤ n, A ∈ℝ^k × n, K ∈ℝ^m × n with k>0 are matrices, b ∈ℝ^k is a given vector, and ρ > 0 is a given Tikhonov regularization parameter. This is a discretized version of problem (<ref>), with the objective functional chosen to match the type of optimal control problem set up in Section <ref>. As a natural choice for discretizing the optimal control problem, we can take A = I, where I is the n × n identity matrix. However, for the generality of our analysis, we only require that the row vectors of A and K span ℝ^n. Under this assumption, we know that for any α > 0, the matrix G_α is invertible, where G_α is defined as

    G_α = A^T A + α K^T K.


We will start from here to discuss how to solve problem (<ref>)-(<ref>).




 §.§ Explicit Solution

We note that problem (<ref>)-(<ref>) has an explicit solution that can be computed using a Lagrangian formulation. We consider the Lagrangian form:

    L(u,y) = J(u,y) + ζ (Ku - y) = 1/2Au - b^2 + ρ/2y^2 + ζ^T (Ku - y),

with ζ∈ℝ^m as an auxiliary Lagrangian parameter vector. By differentiating (<ref>) with respect to u, y, ζ respectively, the first-order optimality conditions <cit.> are given by the following equations:

    {[ A^T A u - A^T b + K^T ζ = 0 ,;                  ρ y - ζ = 0,;                   Ku - y = 0. ].

Solving this system results in:

    {[        û = (A^T A + ρ K^T K)^-1 A^T b,; ŷ = Kû = K (A^T A + ρ K^T K)^-1 A^T b. ].

Thus, we have found (û, ŷ) to be the analytical solution to problem (<ref>)-(<ref>), and this notation will continue to be used throughout this paper. However, computing the inverse matrix can be challenging when dealing with large-scale systems, and its numerical stability may become an issue <cit.>. Therefore, in practice, the penalty method is often preferred for implementation, as it is more suitable for applying iterative methods that do not require direct computation of the inverse <cit.>. As introduced in (<ref>), the problem is formulated as follows: Find u ∈ℝ^n and y ∈ℝ^m to minimize

    (u,y) = J(u,y) + λ/2Ku - y^2 = 1/2Au - b^2 + ρ/2y^2 + λ/2Ku - y^2,

where λ > 0 is a penalty parameter. For simplicity of notation, we will denote the remainder function corresponding to the constraint (<ref>) as R(u,y), defined as

    R(u,y) = Ku - y^2.
 
The first-order optimality conditions to minimize (<ref>) are given by:

    {[ A^T A u - A^T b + λ K^T K u - λ K^T y = 0 ,;                       ρ y - λ Ku + λ y = 0. ].

Solving this system results in:

    {[                         u^λ = (A^T A + ρλ/ρ + λ K^T K)^-1 A^T b ,; y^λ = λ/ρ + λ Ku^λ = λ/ρ + λ K (A^T A + ρλ/ρ + λ K^T K)^-1 A^T b. ].

Comparing (<ref>) with (<ref>), we observe that as λ tends to infinity, u^λ will converge to û because lim_λ→∞ρλ/ρ + λ = ρ, and y^λ will converge to ŷ. Therefore, although problem (<ref>) is distinct from problem (<ref>)-(<ref>), we can consider (<ref>) with a sufficiently large λ as an acceptable approximation to problem (<ref>)-(<ref>).

However, in practice, as λ increases, the difficulty of solving the unconstrained optimization problem associated with (<ref>) also increases. Mathematically, this can be observed by examining the Hessian matrix for (u,y) with respect to u and y, denoted as . Direct computation shows that

    = [ A^TA+λ K^TK      -λ K^T;        -λ K  (ρ +λ)I_m, ]

where I_m denotes an m× m identity matrix. As λ tends to infinity,  tends to become a singular matrix, making it difficult to find the minimizer of (u,y). 

To illustrate this, we present a graphical example. In Figure <ref>, we consider a one-dimensional problem with m=n=1 and aim to find u, y∈ℝ that minimize

    J^*(u,y) = 1/2| u - 2|^2 + 1/2| y|^2,   subject to  2u = y.

The corresponding penalty formulation is to find u, y∈ℝ such that they minimize

    P^λ,*(u,y) = 1/2| u - 2|^2 + 1/2| y|^2 + λ/2| 2u-y|^2.

Consistent with the notation used above, we denote the exact solution as (û, ŷ) and mark it as a red point in the figure. The solution for the minimization problem with λ_1 = 5, denoted as (, ), is marked as a blue point, while the solution for the minimization problem with λ_2 = 0.5, denoted as (, ), is marked as a green point. The range for u and y is chosen to be [-0.5, 2]. We can observe that (, ) is closer to the exact solution than (, ). 

Additionally, Figure <ref>, parts (a) and (b), plot the contour of each level set of P^λ,* for these two values of λ. When λ is relatively larger, the contour is more flattened, indicating that the conditioning of the problem is worse <cit.>. This results in greater difficulty in finding the optimal value, which aligns with our analysis above.





In a neural network setting, a large λ makes the network difficult to train and may not yield the correct solution. Conversely, if λ is not sufficiently large, the solution to the corresponding unconstrained problem may not satisfy the constraint adequately, and u^λ might not approximate u^* well. Therefore, we seek a practical method that ensures the solution satisfies the constraints while effectively approximating the true solution.





 §.§ Penalty Adversarial Problem

To resolve the problem mentioned in the end of last subsection, now, we propose a new unconstrained optimization approach based on the penalty method. Instead of using a single penalty parameter λ, we simultaneously consider two problems with different penalty parameters. Let λ_1 > λ_2 > 0, and let (u^λ_1, y^λ_1) and (u^λ_2, y^λ_2) denote the solutions that minimize the functionals 𝒫^λ_1(u,y) and 𝒫^λ_2(u,y), respectively. We assume that the former problem is hard to solve in practice, while the latter is easier to solve. The corresponding objective functionals for (u^λ_1, y^λ_1) and (u^λ_2, y^λ_2) can be computed using (<ref>) and are denoted as J(u^λ_1, y^λ_1) and J(u^λ_2, y^λ_2). We now consider the following problem:
Find u ∈ℝ^n and y ∈ℝ^m to minimize

    (u,y) = {   J(u,y) + λ_1/2 R(u,y) + ω[J(u, y) - J(u^λ_2, y^λ_2)]^2,      if  (u,y) ∈Ω_1,
       J(u,y) + λ_1/2 R(u,y),      if  (u,y) ∈Ω_2.
            .

Here, ω > 0 is a tunable parameter, and 

    Ω_1 = {(u,y): J(u, y) > J(u^λ_2, y^λ_2)},   Ω_2 = {(u,y): J(u, y) ≤ J(u^λ_2, y^λ_2)}.

Clearly, Ω_1 is an open set and Ω_2 is a closed set due to the continuity of J(u,y). We call this minimization problem the penalty adversarial problem (PAP).

Compared to the standard penalty problem, (u,y) incorporates an additional term to penalize the failure of achieving a sufficiently small objective functional. By designing this problem, we aim to find a solution that adheres closely to the constraint while maintaining feasibility in its solvability. The problem's feasibility is difficult to assert through a mathematical criterion, as it is highly problem-dependent. However, we will demonstrate this advantage through graphical examples in Section <ref> and numerical examples in Section <ref>. On the other hand, a mathematical measure to determine if a solution better adheres to the constraints is straightforward to obtain, which is the value of the remainder function evaluated at the solution. Therefore, mathematically, we aim to show that by choosing an appropriate value of ω, the minimizer of (u,y), denoted as (, ) hereafter, will satisfy the following:

    R(, ) < R(, ).

Here R(u,y) is defined in (<ref>). As long as this holds, it can be seen that the minimizer for (u,y) is more accurate than the minimizer for (u,y) in terms of satisfying the PDE constraints. Proving this argument will be the main focus of our analysis hereafter.

First, consider the relationship between J(u,y) and R(u,y): For a fixed λ∈ℝ^+, (u^λ, y^λ) minimizes (u,y), which is a combination of two parts: one involves the objective functional J(u,y), while the other comprises the remainder function R(u,y). If the value of J(u,y) is fixed, then minimizing (u,y) is equivalent to minimizing R(u,y). Therefore, intuitively, the fact that (u^λ, y^λ) minimizes (u,y) indicates that (u^λ, y^λ) minimizes the remainder function R(u,y) among all (u,y) ∈ℝ^n ×ℝ^m such that J(u,y) = J(u^λ, y^λ). The following lemma justifies this.



    
    For a fixed λ > 0, let (u^λ, y^λ) ∈ℝ^n ×ℝ^m minimize P^λ(u,y) as given in (<ref>). If (u,y) ∈ℝ^n ×ℝ^m, then we have the following assertions:
    
        
  * If J(u,y) ≤ J(u^λ, y^λ), then R(u,y) ≥ R(u^λ, y^λ).
        
  * If R(u,y) ≤ R(u^λ, y^λ), then J(u,y) ≥ J(u^λ, y^λ).
        
  * In each of the previous two assertions, equality can only be attained simultaneously. Specifically, if J(u,y) < J(u^λ, y^λ), then R(u,y) > R(u^λ, y^λ). If R(u,y) < R(u^λ, y^λ), then J(u,y) > J(u^λ, y^λ).
    


    We will only prove the first argument and the proof for the second will follow in the same way. By definition, (u^λ, y^λ) minimizes (u,y). Thus, for any (u,y) ∈ℝ^n ×ℝ^m, we have (u,y) ≥(u^λ, y^λ). This implies that
    
    J(u,y) + λ/2 R(u,y) ≥ J(u^λ, y^λ) + λ/2 R(u^λ, y^λ).

    Since J(u,y) ≤ J(u^λ, y^λ), it follows that
    
    R(u,y) ≥ R(u^λ, y^λ),

    as λ > 0. The third argument also follows immediately from (<ref>) when the inequality relation between J(u,y) and J(u^λ, y^λ) is strict.



This result indeed provides a sufficient condition to restrict (u,y) in domain Ω_1. Namely, using the third argument of Lemma <ref>, as long as R(u,y) < R(, ), then (u,y) ∈Ω_1 and the value of (u,y) will differ from (u,y) at these points. Another important deduction we can make from here is that the minimum of (u,y) in the closed region Ω_2 could always be attained at (, ). We conclude this in the following result:


    
    If (u,y) ∈Ω_2, then 
    
    (u,y) ≥ J(, ) + λ_1/2 R(, ).

    In other words,
    
    min_(u,y) ∈Ω_2(u,y) = J(, ) + λ_1/2 R(, ).



    (u,y) ∈Ω_2 implies that J(u,y) ≤ J(, ). According to the first assertion in Lemma <ref>, we know that
    
    R(u,y) ≥ R(, ).

    On the other hand, since (, ) minimizes (u,y), we have
    
    J(u,y) + λ_2/2 R(u,y) = (u,y) ≥(, ) = J(, ) + λ_2/2 R(, ).

    Combining these two inequalities, we get
    
    J(u,y) + λ_1/2 R(u,y)    = J(u,y) + λ_2/2 R(u,y) + λ_1 - λ_2/2 R(u,y)
       ≥ J(, ) + λ_2/2 R(, ) + λ_1 - λ_2/2 R(, )
       = J(, ) + λ_1/2 R(, ).

    This proves the claim.


Revealed by this, if for any (u,y) ∈Ω_1, (u,y) is always greater than [J(u^λ_2, y^λ_2) + λ_1/2 R(,)], then (u^λ_2, y^λ_2) can be a solution to the minimization problem related to (<ref>). In this case, our claim that (, ), the minimizer of this problem, will always adhere more closely to the constraint than (, ) would no longer hold. Therefore, to make our claim valid, it is necessary for the minimizer (, ) to fall within Ω_1. Enlightened by this, we present the following result:



Assuming ω > 0 is a fixed positive constant parameter, the following two arguments are equivalent:


    
  * (u,y) ∈Ω_1 and 
    
    (u,y) < J(u^λ_2, y^λ_2) + λ_1/2 R(, ),

    
  * (u,y) satisfies R(u,y) < R(u^λ_2, y^λ_2) and
    
    J(u^λ_2, y^λ_2) < J(u,y) < J(u^λ_2, y^λ_2) + λ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ]/1 + √(1 + 2ωλ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ]).





    We first assume the former argument holds and try to deduce the latter one. To start, we define the difference between J(u,y) and J(u^λ_2, y^λ_2) as D^λ_1, λ_2(u,y), namely,
    
    D^λ_1, λ_2(u,y) = J(u,y) - J(u^λ_2, y^λ_2).

    As (u,y) ∈Ω_1, (u,y) satisfies
    
    D^λ_1, λ_2(u,y) > 0.

    Using the definition of (u,y) given in (<ref>), we see that (<ref>) is equivalent to
    
    ω[ D^λ_1, λ_2(u,y) ]^2 + D^λ_1, λ_2(u,y) + λ_1/2[ R(u,y) - R(u^λ_2, y^λ_2) ] < 0.

    Since ω > 0, this inequality will have solutions only when the discriminant of the corresponding quadratic formula is greater than zero, which is equivalent to:
    
    1 - 2ωλ_1 [ R(u,y) - R(u^λ_2, y^λ_2) ] > 0.

    In this case, D^λ_1, λ_2(u,y) should satisfy
    
    -1 - √(1 + 2ωλ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ])/2ω < D^λ_1, λ_2(u,y) < -1 + √(1 + 2ωλ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ])/2ω.

    Since we require that D^λ_1, λ_2(u,y) > 0, it is necessary to have
    
    -1 + √(1 + 2ωλ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ])>0,

    which is equivalent to
    
    R(u,y) < R(u^λ_2, y^λ_2).

    With this, J(u,y) can be estimated as
    
    J(u,y) = J(, ) + D^λ_1, λ_2(u,y) < J(, ) + √(1 + 2ωλ_1 [ R(u^λ_2, y^λ_2) - R(u,y) ]) - 1/2ω,

    which is equivalent to (<ref>).
    This proves the latter argument. 
    
    For the other direction, Lemma <ref> ensures that R(u,y) < R(, ) implies (u,y) ∈Ω_1. In addition, (<ref>) implies that the value of D^λ_1, λ_2(u,y) will ensure (<ref>), which is equivalent to (<ref>), thereby completing the proof.


Here, we have provided an equivalent condition to characterize the case where the minimum of (u,y) falls in Ω_1. As long as this condition is met, the minimum point (, ) will correspond to a smaller value of R(u,y) compared to R(,), which is the desired outcome. The remaining task is to justify this condition. To do so, we need to demonstrate that there exists a point in Ω_1 that satisfies (<ref>), and then we can invoke the continuity of J(u,y) to conclude the existence of a minimum point. This leads to the main result of this section.

    If K≠ 0 and there exists σ > 0 such that
    
    λ_1/2 R(u^λ_2, y^λ_2) > σ + J(û, ŷ) - J(, ),

    where (û, ŷ) given in (<ref>) is the exact solution that minimizes J(u,y) subject to Ku = y, then there exists ω > 0 such that (, ), which minimizes (u,y) defined in (<ref>), satisfies
    
    R(, ) < R(u^λ_2, y^λ_2).




    Recalling the definition for (û, ŷ), we will show that (û, ŷ) satisfies 
    
    (û, ŷ) < J(u^λ_2, y^λ_2) + λ_1/2 R(, ).

    Using Lemma <ref>, it is equivalent to show R(û, ŷ) < R(, ) and (û, ŷ) satisfies (<ref>).
    
    It is clear that R(û, ŷ) = 0. On the other hand,
    
    R(u^λ_2, y^λ_2)    = Ku^λ_2 - y^λ_2^2 = Ku^λ_2 - λ_2/ρ + λ_2 Ku^λ_2^2
            = ρ/ρ + λ_2 Ku^λ_2^2
           
            > 0,

    
    according to our assumption. Thus, the relation R(û, ŷ) < R(, ) holds. 
    
    Meanwhile, using the fact that R(û, ŷ) = 0 again, (<ref>) in this case reduces to 
    
    J(u^λ_2, y^λ_2) < J(û,ŷ) < J(u^λ_2, y^λ_2) + λ_1 R(u^λ_2, y^λ_2)/1 + √(1 + 2ωλ_1 R(u^λ_2, y^λ_2)).

    Since σ > 0 and (<ref>) holds, we can find ω > 0 to satisfy this condition. Therefore, we have shown (<ref>). This indicates that (, ) is not the minimizer of (u,y) for small enough positive ω, as (u,y) reaches a smaller value at (û, ŷ) compared to (, ). It remains to show that any minimizer of (u, y) satisfies R(, ) < R(u^λ_2, y^λ_2). 

    To see this, we consider a subdomain of Ω_1 = {(u,y) : J(u,y) ≥ J(u^λ_2, y^λ_2)}, the closure of Ω_1, defined as
    
    Ω^λ_1, λ_2 = {(u,y) : J(u^λ_2, y^λ_2) ≤ J(u,y) ≤ J(u^λ_2, y^λ_2) + λ_1/2 R(u^λ_2, y^λ_2)}.

    Due to the continuity of J(u,y) with respect to (u,y) and the definition of J(u,y), we see that Ω^λ_1, λ_2 is a closed and bounded domain. Hence, it is compact. When (u,y) ∈Ω_1 ∖Ω^λ_1, λ_2,
    
    (u,y) ≥ J(u,y) > J(u^λ_2, y^λ_2) + λ_1/2 R(u^λ_2, y^λ_2) = (, ),

    and so it cannot be a minimizer of (u,y) in Ω_1. On the other hand, as a continuous function defined on a compact set Ω^λ_1, λ_2, (u,y) attains its minimum at some point in Ω^λ_1, λ_2. For consistent notation, let us denote this point as (, ). Therefore, 
    
    J(, ) + λ_1/2   R(, ) + ω[J(, ) - J(u^λ_2, y^λ_2)]^2
       = (, ) ≤(û, ŷ) < J(u^λ_2, y^λ_2) + λ_1/2 R(u^λ_2, y^λ_2).

    Finally, we will show that R(, ) < R(, ). By Lemma <ref>, it is sufficient to show that J(, ) > J(, ). Therefore, we need to prove that J(, ) ≠ J(, ) since (, )∈Ω^λ_1, λ_2. In fact, if J(, ) = J(, ), then by Lemma <ref>, we know that R(, ) ≥ R(, ). From (<ref>), we have
    
    J(, ) + λ_1/2 R(, )
       ≤ J(, ) + λ_1/2 R(, ) 
       = J(, ) + λ_1/2 R(, ) + ω[J(, ) - J(u^λ_2, y^λ_2)]^2 
       =(, )≤(û,ŷ)< (,)= J(u^λ_2, y^λ_2) + λ_1/2 R(u^λ_2, y^λ_2).

    This is a contradiction, and so we conclude that J(, ) is strictly larger than J(, ) and thus R(, ) is strictly smaller than R(, ). This completes the proof.


To conclude this part, we will make two comments on our strategy of proof for Theorem <ref>.


     The proof for Theorem <ref> presented here focuses on the existence of a minimizer, which is sufficient to support our claim that R(, ) < R(, ). For a deeper study of the properties of (, ), one can consider a quantitative analysis of the functional form of (u,y). 



    In this context, we use (û, ŷ) as a reference point to demonstrate that there exists at least one point satisfying (<ref>). However, this choice is not mandatory. In fact, any point in Ω_1 that can be computed and shown to result in a value smaller than (, ) could replace (û, ŷ) in our proof.





 §.§ Choice of ω


The last part discusses the possibility of designing a penalty adversarial problem to ensure its minimizer adheres to the constraints more closely than (,). Once λ_1 and λ_2 are fixed, the choice of ω will determine the formulation of the corresponding problem. In this section, we will explore this choice's influence through theoretical analysis and computational examples.



  §.§.§ Upper Bound of ω


We note that (<ref>) provides an explicit condition that can be used to determine if the value of λ_1 is large enough to implement this penalty adversarial strategy independent of ω. Then, after fixing the values for λ_1 and λ_2, based on this proof, we can provide an estimate for the upper bound of ω. This conclusion is formalized in the following proposition.


    Assuming (<ref>) holds, then ω > 0 should satisfy the following relation:


    ω < λ_1 R(u^λ_2, ) - 2[J(û, ŷ) - J(, )]/2[J(û, ŷ) - J(, )]^2.


    In the proof of Theorem <ref>, we deduced the following relation:
    
    J(û, ŷ) < J(u^λ_2, y^λ_2) + λ_1 R(u^λ_2, y^λ_2)/1 + √(1 + 2ωλ_1 R(u^λ_2, y^λ_2)).

    Rearranging this inequality results in (<ref>).


We will provide an intuitive understanding of why there exists an upper bound on the choice of ω in Section <ref>.






  §.§.§ Comparison between (u,y) and (u,y)




We have discussed the advantage of minimizing (u,y) over (u,y), specifically its ability to adhere to the constraint more effectively. On the other hand, the main disadvantage of trying to minimize (u,y), as mentioned earlier, is its poor conditioning, which makes it difficult to solve. Therefore, by choosing to minimize (u,y) instead of (u,y), we aim to make the corresponding minimization problem easier to solve.

Here, we will present a graphical example to illustrate this point. Using the same example as in Section <ref>, we consider a one-dimensional problem with m = n = 1 and aim to find u, y ∈ℝ to minimize



    J^*(u,y) = 1/2| u - 2|^2 + 1/2| y|^2,   subject to  2u = y.

The corresponding penalty formulation has been stated in (<ref>), and the corresponding penalty adversarial problem is as follows:
find u ∈ℝ and y ∈ℝ to minimize


    (u,y) = {   1/2| u - 2|^2 + 1/2| y|^2 + λ_1/2| 2u - y|^2 + ω[1/2| u - 2|^2 + 1/2| y|^2 - J(u^λ_2, y^λ_2)]^2,      if  (u,y) ∈Ω_1,
       1/2| u - 2|^2 + 1/2| y|^2 + λ_1/2| 2u - y|^2,      if  (u,y) ∈Ω_2.
            .

where 

    J(, ) = 1/2| - 2|^2 + 1/2||^2

is a fixed number that can be calculated using (<ref>) once λ_2 is fixed. 

We plot the level sets of (u,y) and (u,y) with different parameters in Figure <ref>. To provide more information, we extend the range of each figure to [-2, 3] instead of [-0.5, 2] as in Figure <ref>. As shown in Figure <ref>, the first row displays the contour plots of (u,y) with λ_1 = 5 and λ_2 = 0.5, while the second row displays the contour plots of (u,y) with the same values for λ_1 and λ_2, and three different choices of ω: 0.1, 1, and 10. The points (û, ŷ), (, ), and (, ) are marked as red, blue, and green points, respectively, in each figure. In the second row, a fourth point representing the minimizer of (u,y) for the corresponding value of ω is also marked.






We make the following observations:


    
  * The point (, ) always lies between (, ) and (, ) and is closer to the analytical solution (û, ŷ) than (, ).
    
  * The smaller the value of ω, the closer (, ) is to (, ). In Figure <ref> (c), (, ) almost overlaps with (, ) for ω = 0.1.
    
  * As the value of ω increases, the contour shape changes from elliptical to bowl-shaped, with the direction towards (, ) generally expanding.


The third point illustrates why it is easier to minimize (u,y) compared to (u,y). At points away from the minimizer, for a suitable choice of ω, the condition number is more favorable since the contour of the corresponding level set is not as flattened as it is for (u,y). This behavior enables the effective use of standard iterative methods to find the minimizer.



  §.§.§ Relation to Auto-Tuning Parameter Strategy

To address the central challenge discussed in this paper regarding the penalty method—specifically, the difficulty of selecting a suitable parameter, as small parameters lead to inaccurate solutions and large parameters lead to ill-conditioned problems—a common strategy is to propose an approach with varying penalty parameters. This strategy is also commonly used when employing neural networks to solve PDE-related problems; for example, see <cit.>. 

To illustrate this, consider using λ^k to denote the penalty parameter used in the k-th iteration, and as an example given in <cit.>, one possible choice is to take λ^k+1=β^k λ^k,
where β^k > 0 is a constant parameter selected for the k-th iteration. The main difficulty with this strategy is that the choices of {β^k}_k are problem-dependent <cit.>, often requiring fine-tuning in practice. Such a process is usually challenging or highly technical. It would be favorable if there were an automatic strategy to tune the value of λ^k. Here, we point out that the penalty adversarial strategy proposed here provides such an automatic mechanism. Specifically, when J(u, y) ≥ J(u^λ_2, y^λ_2), we have

    (u,y)    = J(u,y) + λ_1/2 R(u,y) + ω[J(u, y) - J(u^λ_2, y^λ_2)]^2,

and its gradient can be computed as

    ∇(u,y)=(1+2ω(J(u,y)-J(,)))∇ J(u,y)+λ_1/2∇ R(u,y).

On the other hand, for (u,y) with a general λ>0, its gradient can be computed as

    ∇(u,y) = ∇ J(u,y)+λ/2∇ R(u,y).

If one is using an iterative method to find the minimizer and employs an explicit scheme, then minimizing (u,y) at a fixed stage (u,y) is equivalent to minimizing P^λ̃(u,y) in the sense that their gradients have the same direction, with 

    λ̃ = λ_1/1 + 2ω(J(u,y) - J(u^λ_2, y^λ_2)) < λ_1.

Thus, the penalty adversarial problem becomes more accessible to minimize under such conditions, as it is equivalent to using a smaller penalty parameter than λ_1 when J(u,y)>J(,). In addition, as J(u,y) gets closer to J(,), λ̃ will approach λ_1.

To conclude, solving the penalty adversarial problem is a strategy to dynamically tune the penalty parameter according to the closeness of the solution to the real solution. When the solution is far from the correct one, the corresponding problem will be equivalent to minimizing P^λ(u,y) with a small λ, allowing it to converge to the actual solution more quickly. As the solution nears the correct solution, the problem transitions to one with a larger penalty term, thereby adhering more closely to the constraints. Compared to other traditional penalty methods with dynamic penalty terms, this approach is more straightforward to implement since no specific strategy is needed to determine the dynamics of the penalty parameters in advance.

To this end, we use this equivalence to provide an intuitive reasoning for the existence of an upper bound for the choice of ω. As seen from (<ref>), if ω is too large, when J(u,y)>J(,), λ̃ can become very small, even smaller than λ_2. In this case, minimizing (u,y) behaves similarly to minimizing a penalty problem with a penalty parameter smaller than λ_2. Consequently, the assertion that R(,) < R(,) will no longer hold.





 §.§ Alternate Formulations of Penalty Adversarial Problem


As implemented in (<ref>), the additional penalty term associated with ω is chosen as a quadratic form of the difference | J(u,y) - J(, ) |. However, this choice is not unique. Other non-negative functions of | J(u,y) - J(, ) | can also be used to add this penalty term. A natural choice is to use a monomial with a power of k, and the corresponding function can be defined as:


    (u,y) = {   J(u,y) + λ_1/2 R(u,y) + ω| J(u, y) - J(u^λ_2, y^λ_2)|^k,      if  (u,y) ∈Ω_1,
       J(u,y) + λ_1/2 R(u,y),      if  (u,y) ∈Ω_2.
            .


The corresponding penalty adversarial problem will be modified to find the minimizer of (u,y) instead of (u,y) = A^λ_1,λ_2_ω,2(u,y). Using the same example as presented in (<ref>), we plot the contour of level sets for the corresponding functional (u,y) with k ranging from 1 to 9, and fixed parameters: λ_1 = 5, λ_2 = 0.5, ω = 5. As shown in Figure <ref>, the corresponding minimizer, denoted as (, ), moves closer to the exact solution as the penalty power order k increases. Additionally, as k increases, the shape of the contour generally becomes more bowl-like: the curvature of the part near the actual solution remains similar, while the curvature on the other side becomes more flattened, and the curve extends to be wider. This property enhances the solvability of the problem. Moreover, we can observe that (,) always lies between the point (,) and (,), and it remains closer to the actual solution than (,). In conclusion, in the previous subsection, we considered the case of k = 2 for simplicity of analysis, but in practice, other values of k could also be valid choices.








§ PENALTY ADVERSARIAL NETWORK


This section presents the neural network-based algorithm inspired by the penalty adversarial problem discussed in the previous section. Following the standard approach in physics-informed neural networks <cit.>, the input to the neural network consists of the spatial variable x and the temporal variable t, while the output is the neural network solution to the corresponding optimal control problem constrained by PDEs, as introduced in (<ref>). We denote this neural network solution as ((x,t;θ), (x,t;θ)), with parameters θ learned to minimize the corresponding loss function L[x, t; θ]. By defining the loss function in different ways, various methods can be implemented to solve optimal control problems through neural networks.

As introduced above, a standard approach <cit.> is to add the constraints as penalty terms to the objective function and then minimize it. We denote the loss function as L_P[x, t; θ], which is defined as:



    L_P[x, t; θ] 
        = 1/N_J∑_m=1^N_JJ[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)]_Objective loss + λ_p/N_p∑_m=1^N_pℒ[(x_m^p, t_m^p; θ), (x_m^p, t_m^p; θ)]^2_PDE residual loss
           + λ_b/N_b∑_j=1^N_bℬ[(x_j^b, t_j^b; θ), (x_j^b, t_j^b; θ)]^2_Boundary loss + λ_i/N_i∑_n=1^N_iℐ[(x_n^i, t_n^i; θ), (x_n^i, t_n^i; θ)]^2_Initial loss.

The variables in the loss function are defined as follows: x and t represent the vectors of spatial and temporal variables used as input to the neural network, which are collections of sample points {x_m^J}_m=1^N_J, {x_m^p}_m=1^N_p, {x_j^b}_j=1^N_b, {x_n^i}_n=1^N_i, and {t_m^J}_m=1^N_J, {t_m^p}_m=1^N_p, {t_j^b}_j=1^N_b, {t_n^i}_n=1^N_i, respectively. Here, N_J, N_p, N_b, and N_i denote the number of sample points for the objective functional, PDE residual, boundary conditions, and initial conditions, respectively. The objective loss, evaluated at sample points (x_m^J, t_m^J), is denoted by J[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)]. The PDE residual loss, boundary loss, and initial loss, evaluated at their respective sample points, are represented by ℒ[(x_m^p, t_m^p; θ), (x_m^p, t_m^p; θ)], ℬ[(x_j^b, t_j^b; θ), (x_j^b, t_j^b; θ)], and ℐ[(x_n^i, t_n^i; θ), (x_n^i, t_n^i; θ)]. We note that the loss function based on the penalty formulation is slightly different from the general penalty formulation we discussed above, as we choose not to divide the penalty parameters by 2 to keep the form of the loss functions simple. One can, however, use the previous formulations to construct the loss functions here.

In contrast to this penalty approach, we introduce an adversarial network structure enlightened by the penalty adversarial problem. Following the standard terminology used in adversarial network related works <cit.>, we call the two different network solver network and discriminator network, and denote their solution respectively as ((x,t;θ),(x,t;θ)) and ((x,t;θ),(x,t;θ)). The discriminator network is simply minimizing (<ref>) with some small penalty parameters. Namely, it focuses on minimizing  

    L^d[x, t; θ] 
        = 1/N_J∑_m=1^N_J 
       J[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)] + λ_p^d/N_p∑_m=1^N_pℒ[(x_m^p, t_m^p; θ), (x_m^p, t_m^p; θ)]^2 
           + λ_b^d/N_b∑_j=1^N_bℬ[(x_j^b, t_j^b; θ), (x_j^b, t_j^b; θ)]^2 + λ_i^d/N_i∑_n=1^N_iℐ[(x_n^i, t_n^i; θ), (x_n^i, t_n^i; θ)]^2,

where λ_p^d, λ_b^d, λ_i^d>0 are penalty weights for PDE loss, boundary loss and initial loss of discriminator network, respectively. In practice, we usually choose them to be relatively small to guarantee ease of training. On the other hand, for the solver network, in addition to then standard penalty formulation, we consider an extra term measuring the difference between objectives values of the discriminator network and the solver network. Namely, it minimizes the following loss function:


    L^s[x, t; θ] 
        = 1/N_J∑_m=1^N_J 
       J[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)] + λ_p^s/N_p∑_m=1^N_pℒ[(x_m^p, t_m^p; θ), (x_m^p, t_m^p; θ)]^2 
           + λ_b^s/N_b∑_j=1^N_bℬ[(x_j^b, t_j^b; θ), (x_j^b, t_j^b; θ)]^2 + λ_i^s/N_i∑_n=1^N_iℐ[(x_n^i, t_n^i; θ), (x_n^i, t_n^i; θ)]^2
           +ω|1/N_J∑_m=1^N_J J[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)]-1/N_J∑_m=1^N_J J[(x_m^J, t_m^J; θ), (x_m^J, t_m^J; θ)]|^2,

where ω is a tunable hyperparameter which plays the same role as it is in the penalty adversarial problem.

In practice, we initialize both networks and training data, then train these two networks together. Thanks to the fast development of machine learning toolboxes, the training of this problem is relatively standard and can be implemented directly in the TensorFlow framework <cit.> since it supports automatic differentiation to calculate derivatives of the loss functions with respect to the weights. Backpropagation <cit.> is then applied to update the weights in the network. We summarize the training process in Algorithm <ref>.




To conclude this section, we note that Algorithm <ref> is not exactly the neural network version of the penalty adversarial problem discussed in Section <ref>. The penalty adversarial problem uses the exact objective value at (,) to compute the additional penalty term, which is a fixed number. However, in the neural network algorithm, we use (,), which is not the exact minimum but a solution to the discriminator network trained simultaneously with the solver network. To recover a neural network version of the penalty adversarial problem, it is possible to train the discriminator as a surrogate network in advance and use it to construct the penalty term. At this stage, we are not able to assert which approach is better. We present this work with the current choice for its simplicity in implementation, as it provides an all-at-once solution.




§ NUMERICAL RESULTS


We will present and discuss the performance of the penalty adversarial network when it is applied to solve optimal control problems constrained by different types of equations, including both linear and nonlinear problems. To begin, we will provide additional details to complement Algorithm <ref> for practical implementation.




 §.§ Learning Rate Scheduling


Algorithm <ref> outlines a general workflow for simultaneously training the discriminator and solver networks. In practice, to enhance the training process, we employ a learning rate scheduling strategy <cit.>. Specifically, we set a certain number of epochs as the patience parameter P. If the loss has not improved after P epochs and the learning rate is still larger than a preset minimum learning rate, the learning rate is reduced to half of its original value. This approach helps automatically fine-tune the training process by decreasing the learning rate when performance stagnates, allowing for more precise adjustments. Additionally, during training, the network often gets stuck in local minima in the first few epochs, so we typically discard the initial epochs and begin recording the best weights only after this initial phase.



 §.§ Numerical Examples


Here, we will provide a few examples demonstrating the effectiveness of our proposed strategy when applied to various types of optimal control problems constrained by different equations. In each instance, we will present the effect of the direct penalty formulation with a large penalty parameter and then compare the results with the application of our proposed penalty adversarial approach. We will consider three different problems constrained by the 1D Poisson equation, 2D Poisson equation, and 2D Allen-Cahn equation, respectively. 





  §.§.§ Example 1: Boundary Control Problem Constrained by 1D Poisson Equation


In this example, we consider a boundary control problem constrained by a 1D Poisson equation, following the problem setup described in (<ref>). We implement Algorithm <ref> and the learning rate adjustment strategy discussed in Section <ref>. The specific optimal control problem we consider is given by:


    Minimize       J(u) = 1/2∫_0^1 [ u(x) - u_d(x) ]^2   dx + ρ/2[ u(0)^2 + u(1)^2 ], 
       subject to       -d^2 u/dx^2 = A sin(2π x),    x ∈ [0,1],

where the desired state u_d(x) = A/4π^2sin(2π x) + bx + a, with parameters a = -10, b = 65, and A = 8π^2. The control for this problem is the boundary values u(0) and u(1). In fact, if the boundary values are fixed, the solution to the Poisson equation is uniquely determined, which in turn determines the value of J(u). 

This problem has an analytical solution, given by:

    u_analytical(x) = A/4 π^2sin(2 π x) + b^* x + a^*,

where a^* = 1/1 + 2 ρ a + 2 ρ/(1 + 2 ρ)(1 + 6 ρ) b and b^* = b/1 + 6 ρ. In this example, we take ρ = 2, resulting in a^* = 2 and b^* = 5.

Firstly, we consider using a standard penalty formulation, as given in (<ref>), with a large penalty parameter to solve the problem. Namely, the corresponding neural network aims to output (,) to minimize the following loss function:

    L[x; θ] = ρ/2[ (0;θ)^2 + (1;θ)^2 ]   +1/2N∑_m=1^N[ (x_m;θ) - u_d(x) ]^2
        + λ_p/N∑_m=1^N[d^2 /dx^2(x_m;θ) + A sin(2π x_m)]^2.

The training data {x_m}_m are N = 32 points uniformly chosen from [0,1]. The hyperparameters for this example chosen in implementation are as follows. We set the neural network depth to 4, width to 40, the penalty parameter λ_p=5000. The initial learning rate is 0.001, which can be reduced to a minimum learning rate of 0.0001 based on the learning rate scheduling strategy with patience set to be 3000. The training is conducted over a maximum of 200000 epochs. The numerical findings are presented in Figure <ref>.



We observe that after 200000 epochs, the loss function continues to decrease. However, the numerical solution computed via the neural network remains significantly distant from the analytical solution. The derivative of the numerical solution exhibits a nearly constant error compared to the analytical derivative. In contrast, the second derivative of the numerical solution, which corresponds to the constraint in this problem, remains close to its analytical value. This outcome aligns with our expectations, as we enforce a large penalty parameter on the PDE constraint. While this enforcement ensures the solution adheres well to the constraint, it also makes it difficult for the network to be trained effectively and to reach the actual solution.

It is important to note that the training process here is standard without applying specific techniques. While refining the training procedure or increasing the number of training epochs can bring the numerical solution closer to the analytical solution, the results demonstrate that a basic neural network may struggle to find the optimal solution when a large penalty parameter is employed. In comparison, we now focus on solving the same problem using our proposed penalty adversarial strategy.

Instead of using one single network, to utilize penalty adversarial network approach, we create and train the solver and the discriminator network together. Namely, the discriminator network aims to output (,) to minimize the following loss function

    L^d[x; θ] = ρ/2[ (0;θ)^2 + (1;θ)^2 ]   +1/2N∑_m=1^N[ (x_m;θ) - u_d(x) ]^2
        + λ_p^d/N∑_m=1^N[d^2 /dx^2(x_m;θ) + A sin(2π x_m)]^2,

and the solver network aims to output (,) to minimize the following loss function

    L^s[x; θ] =   ρ/2[ (0;θ)^2 + (1;θ)^2 ]+1/2N∑_m=1^N[ (x_m;θ) - u_d(x) ]^2
       + λ_p^s/N∑_m=1^N[d^2 /dx^2(x_m;θ) + A sin(2π x_m)]^2
       +ω {ρ/2[ (0;θ)^2 + (1;θ)^2 ]-ρ/2[ (0;θ)^2+(1;θ)^2 ] 
               +1/2N∑_m=1^N[ (x_m;θ) - u_d(x) ]^2-1/2N∑_m=1^N[ (x_m;θ) - u_d(x) ]^2}^2.

The parameters here are chosen as λ_p^d = 1, λ_p^s= 5000 and ω = 1.

Other hyperparameters for the construction of the networks remain the same as in the previous setting. The neural network's depth and width are set to 4 and 40, respectively. The training data consists of N=32 points uniformly drawn from [0,1]. The initial learning rate is set to 0.001, which can be reduced to a minimum learning rate of 0.0001 based on the learning rate scheduling strategy, with patience set to 3000. The training is conducted over a maximum of 200000 epochs.



The numerical findings are presented in Figure <ref>. The first row compares the neural network approximations from both the solver and discriminator networks with the analytical solution for the value of u, explicitly showing the magnitude of the difference between these two numerical solutions and the analytical solution. Both numerical solutions approximate the analytical solution well. However, the solver network's solution is closer to the real solution, with a maximum error of 0.08, while the maximum error for the discriminator network reaches 0.16. 

A more notable difference is observed in the second row, which compares the second-order derivatives. The accuracy of the second-order derivative reflects the degree of constraint satisfaction. Here, the solver network's solution maintains the error of d^2u/dx^2 under 0.025, which is quite small compared to the magnitude of the second-order derivative. However, the discriminator network exhibits a much larger error in this term, reaching a maximum error of 1.2.

This example demonstrates the advantage of using an adversarial network compared to a standard penalty formulation with both large and small penalty parameters. As shown in Figure <ref>, using a large penalty parameter makes the neural network harder to train, resulting in a numerical solution that deviates from the real solution after a certain number of epochs. The adversarial approach produces a much more accurate approximation while maintaining a large penalty parameter to guarantee satisfying the underlying PDE. Conversely, using a small penalty parameter results in a formulation whose theoretical solution is away from the real analytical solution u_analytical, as given in (<ref>). This fact means that even if the discriminator neural network is easier to train and converges, the obtained solution will not be accurate enough, particularly not satisfying the PDE constraint well, as shown in Figure <ref>. The adversarial approach successfully achieves a solution closer to the real analytical solution without requiring complicated training techniques beyond introducing the adversarial structure. This example highlights the effectiveness of applying the penalty adversarial approach.



  §.§.§ Example 2: Distributed Control Problem Constrained by 2D Poisson Equation


Here we will present that our approach will also work for multi-dimensional problems. We consider a distributed control problem constrained by a 2D Poisson equation. The same problem is investigated in <cit.> as well.

The specific optimal control problem we consider is given by:

    Minimize       J(u, f) = 1/2∫_0^1 ∫_0^1 [ u(x, y) - u_d(x, y) ]^2   dx   dy + ρ/2∫_0^1 ∫_0^1 f(x, y)^2   dx   dy, 
       subject to       -Δ  u(x,y) = f(x, y),    (x, y) ∈ [0,1] × [0,1],

with u(x,y)=0 on the boundary and the desired state 

    u_d(x, y) = Asin(π x) sin(π y),
 with parameter A = 10. The control for this problem is the distributed control f(x, y) over the domain. This problem also has an analytical solution, given by:

    u_analytical(x, y) = A/1 + 4 ρπ^4sin(π x) sin(π y),

and the exact control is:

    f_analytical(x, y) = 2 π^2 A/1 + 4 ρπ^4sin(π x) sin(π y).

In this example, we take ρ = 0.01. We present the desired state u_d and the analytical solution u_analytical in Figure <ref>. Under this choice of ρ, there is an apparent difference between the desired state and the analytical solution to the control problem. We deliberately choose such an example instead of a problem with a much smaller regularization parameter (for example, ρ = 0.0001) to ensure that the neural network genuinely finds the solution to the optimal control problem rather than merely solving an approximation problem that approximates the desired state as if learning a given function.




Firstly, as a naive approach, we aim to solve this problem using a penalty formulation. We consider using a neural network with a depth of 4 layers and width of 60 neurons to minimize the following loss function:

    L[x, y; θ] = 
       1/2N^2∑_m=1^N∑_n=1^N   [ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 + ρ/2N^2∑_m=1^N∑_n=1^N  (x_m, y_n; θ)^2  
        + λ_p/N^2∑_m=1^N∑_n=1^N[ ∂^2 /∂ x^2(x_m, y_n; θ) + ∂^2 /∂ y^2(x_m, y_n; θ) + (x_m, y_n; θ) ]^2
             + λ_b/(N_b)∑_k=1^N_b[ (x_k^b, y_k^b; θ) ]^2,

The weights for the equation and boundary loss terms are set to λ_e = λ_b = 2000. The training data used to compute the objective function and equation loss are selected as {x_m}_m=1^N, {y_n}_n=1^N with N = 16, resulting in a total of 256 interior points. The training data used to compute the boundary loss are selected as {(x_k^b,y_k^b)}_k=1^N_b with N_b = 32. Both the interior points and boundary data are uniformly sampled from their corresponding domains. For the training details, the initial learning rate is set to 0.001, which can be reduced to a minimum learning rate of 0.0001 based on the learning rate scheduling strategy, with patience set to 3000 epochs. The training is conducted over a maximum of 450000 epochs. The numerical findings are presented in Figure <ref>.




The results indicate that the numerical solution does not approximate the analytical solution accurately. The exact value for u ranges from 0 to approximately 2.04, while the numerical value for u remains below 0.6, indicating a significant discrepancy. Although the training loss continues to decrease, it is evident that the model does not converge to the correct solution even after 450000 epochs. While we cannot exclude the possibility that the network might eventually converge to the actual solution with more training epochs, this convergence has not been observed at this point after such a large number of epochs. This phenomenon showcases the limitations of a direct approach using penalty formulation with large penalty parameters to implement the PDE constraints.



In contrast, let us apply the penalty adversarial network to solve the same problem. The loss function that the discriminator network aims to minimize is given by:

    L^d[x, y; θ] = 
       1/2N^2   ∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 + ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2  
        + λ_p^d/N^2∑_m=1^N∑_n=1^N[ ∂^2 /∂ x^2(x_m, y_n; θ) + ∂^2 /∂ y^2(x_m, y_n; θ) + (x_m, y_n; θ) ]^2
             + λ^d_b/(N_b)∑_k=1^N_b[ (x_k^b, y_k^b; θ) ]^2.

On the other hand, the solver network aims to minimize:

    L^s[x, y; θ]    = 
       1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 + ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2  
            + λ_p^s/N^2∑_m=1^N∑_n=1^N[ ∂^2 /∂ x^2(x_m, y_n; θ) + ∂^2 /∂ y^2(x_m, y_n; θ) + (x_m, y_n; θ) ]^2
           + λ^s_b/(N_b)∑_k=1^N_b[ (x_k^b, y_k^b; θ) ]^2
           +ω {ρ/2N∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2-ρ/2N∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2
                    +  1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2  
                            -   1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 }^2.

The parameters are chosen as λ_p^s=λ_b^s = 2000, λ_p^d=λ_b^d = 10 and ω = 100. All other parameters related to the structure and the training data remain the same. 

The numerical results of applying the penalty adversarial network to solve problem (<ref>) are presented in Figure <ref>. We observe that the solver successfully finds the real solution. Subplot (c) in Figure <ref> shows that the difference between the exact solution, shown in Subplot (a), and the predicted result from the solver, shown in Subplot (b), is minimal. However, the discriminator fails to reach the exact solution, exhibiting a significant error in the bottom-right area compared to other regions, as shown in Subplot (f).




Focusing on the ability of the networks to adhere to the constraints, as reflected by the PDE losses, we see that the solver's prediction for u and f generally follows the Poisson equation quite well in most areas, except for some deviation at the top boundary. In comparison, while the discriminator's prediction also generally follows the Poisson equation, it deviates strongly in the bottom-right area, corresponding to the same region where the discriminator makes a significant error in predicting the values of u.

Examining the evolution of the loss function over epochs, we observe an interesting phenomenon. As seen in Subplot (i) of Figure <ref>, the loss for the discriminator decreases slowly and remains almost flat during the first 400000 epochs. However, after 400000 epochs, the loss decreases much more rapidly. At the same time, the loss for the solver increases quickly due to the difference in the objective value between the solver and the discriminator, which causes the additional penalty term introduced by the penalty adversarial approach to increase rapidly.

Though we do not fully understand this behavior at present, we propose one potential explanation: both the solver and the discriminator gradually converge toward the analytical solution at first. However, since the discriminator uses a small penalty parameter, theoretically, the minimum that the discriminator network can converge to deviates a certain distance from the analytical solution u_analytical. This deviation arises from a tradeoff between disobeying the constraints and decreasing the objective function, leading to a situation where the corresponding numerical solution might attempt to decrease the objective function by strongly disobeying the PDE constraints in certain areas. This explanation aligns with the observed phenomenon in the discriminator's result, which significantly deviates from satisfying the Poisson equation in a specific region. Comparatively, although the loss function for the solver increases after 400000 epochs, the minimum that the solver network reaches does not change significantly and remains close to the real analytical solution.






  §.§.§ Example 3: Distributed Control Problem Constrained by 2D Allen-Cahn Equation


As the last example presented here, we will showcase that the penalty adversarial network can also be applied to solve nonlinear problems, albeit more complicated. We consider a distributed control problem constrained by a 2D Allen-Cahn equation, which has received considerable research interest from the computational community <cit.>.

The specific optimal control problem has the same formulation for the objective functional as the previous example (<ref>) but is constrained by a different equation, given by:

    Minimize       J(u, f) = 1/2∫_0^1 ∫_0^1 [ u(x, y) - u_d(x, y) ]^2   dx   dy + ρ/2∫_0^1 ∫_0^1 f(x, y)^2   dx   dy, 
       subject to       -Δ  u(x,y) + 1/ϵ^2[u(x,y)^3-u(x,y)] = f(x, y),    (x, y) ∈ [0,1] × [0,1],

with u(x,y) = 0 on the boundary and a chosen desired state. The parameter ϵ > 0 is related to the PDE constraint itself. Using an adjoint formulation <cit.>, we can find the following pair of functions that serve as a solution to problem (<ref>): We set


    u_analytical(x, y) = αsin(π x) sin(π y) + βsin(2π x) sin(2π y),

and

    f_analytical(x, y) 
       =  2π^2 [αsin(π x) sin(π y) + 4βsin(2π x) sin(2π y)]
           + 1/ϵ^2[(αsin(π x) sin(π y) + βsin(2π x) sin(2π y))^3 - αsin(π x) sin(π y) - βsin(2π x) sin(2π y)],

while the desired state is

    u_d(x, y)    = u_analytical(x, y) + ρ Δ^2 u_analytical(x, y) -  3ρ/ϵ^2u_analytical^2(x, y) Δ u_analytical(x, y) 
          - ρ/ϵ^2[  6   u_analytical(x, y) |∇ u_analytical(x, y)|^2 + Δ u_analytical(x, y) ] 
          - ρ/ϵ^2[  Δ u_analytical(x, y) - 1/ϵ^2( u_analytical^3(x, y) - u_analytical(x, y) ) ] (3 u_analytical^2(x, y) - 1).

Here, α, β∈ℝ are parameters. Note that if ρ = 0, then u_d = u_analytical, which coincides with the expectation that, in this case, the minimization problem (<ref>) is simply equivalent to learning a known function. However, when ρ > 0, u_d can have obvious differences from u_analytical. In the example presented here, we choose ϵ = 0.4, α = 0.45, β = 0.55, and ρ = 0.0001. Under these choices, the corresponding desired state u_d and u_analytical are plotted in Figure <ref>. Obvious differences in values at the same points can be observed from this comparison, preventing the network from simply learning a function.




For the sake of brevity of our presentation, we will skip the experiment of using a naive approach as a penalty formulation with a large penalty parameter. As expected, it will fail to work. We will only present the results of applying a penalty adversarial network. For this approach, we define the corresponding loss function for the discriminator and solver network, respectively, as

    L^d[x, y; θ] = 
       1/2N^2   ∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 + ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2  
        + λ_p^d/N^2∑_m=1^N∑_n=1^N[ ∂^2 /∂ x^2(x_m, y_n; θ) + ∂^2 /∂ y^2(x_m, y_n; θ) 
                        -1/ϵ^2((x_m, y_n; θ)^3-(x_m, y_n; θ)) +(x_m, y_n; θ) ]^2
       + λ^d_b/(N_b)∑_k=1^N_b[ (x_k^b, y_k^b; θ) ]^2,

and

    L^s[x, y; θ]    = 
       1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 + ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2  
            + λ_p^d/N^2∑_m=1^N∑_n=1^N[ ∂^2 /∂ x^2(x_m, y_n; θ) + ∂^2 /∂ y^2(x_m, y_n; θ) 
                        -1/ϵ^2((x_m, y_n; θ)^3-(x_m, y_n; θ)) +(x_m, y_n; θ) ]^2
           + λ^s_b/(N_b)∑_k=1^N_b[ (x_k^b, y_k^b; θ) ]^2
           +ω {ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2-ρ/2N^2∑_m=1^N∑_n=1^N (x_m, y_n; θ)^2
                    +  1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2  
                            -   1/2N^2∑_m=1^N∑_n=1^N[ (x_m, y_n; θ) - u_d(x_m, y_n) ]^2 }^2.


The hyperparameters for this experiment are configured as follows: The training data for both the objective function and equation loss are sampled as {x_m}_m=1^N and {y_n}_n=1^N, with N = 32. The boundary training data consist of {(x_k^b, y_k^b)}_b=1^N_b, with N_b = 32. Both the interior and boundary points are uniformly sampled from their respective domains. For the training specifics, the initial learning rate is set at 0.001, potentially decreasing to a minimum of 0.0001 using a learning rate scheduling strategy with a patience parameter of 10000 epochs. The learning rate adjustment does not begin until after 300000 epochs to avoid early-stage inaccuracies. The total number of epochs is extended to 1.5 million to capture the full training dynamics, as the solver network, as seen in the experiment results, continues to improve throughout. This extended training period also highlights the common challenges of addressing nonlinear equations in PINN applications <cit.>.


The penalty parameters selected for this example are as follows: λ_p^s = λ_b^s = 1000, λ_p^d = λ_b^d = 0.2, and ω = 20000. It is important to note that we opted for a relatively large value for ω due to the complexity of the associated nonlinear problem. A small ω would not significantly accelerate the training process or help achieve convergence, which is the main advantage over using a large penalty parameter alone. Although we have established that ω should have an upper bound, as demonstrated in Proposition <ref>, this example illustrates that using a relatively large value for ω in some instances is practical. The results of this experiment are shown in Figure <ref>.




We observe that both the solver and discriminator networks converge to solutions close to the exact solution, with the solver network being notably more accurate than the discriminator network. This is evident in Subfigures (c) and (f), which explicitly show the differences between the solver's and discriminator's predictions compared to the exact solution, respectively. The inaccuracy of the discriminator network stems from the fact that its penalty parameter, set at λ_p^d = λ_b^d = 0.2, is too small to ensure a solution that closely approximates the exact solution of the constrained optimization problem. On the other hand, since ω is chosen to be large in this example, the expected advantage of the solver network adhering better to the constraint compared to the discriminator network becomes less pronounced, as observed in Subfigures (g) and (h). However, even so, we can still observe that the PDE loss for the solver's network is slightly smaller than that of the discriminator's network, aligning with theoretical expectations.

Another important observation is that the discriminator network converges quickly, as expected, due to its small penalty parameter. As shown in Figure <ref>, the discriminator network finds a solution with similar errors to the final result after just 300000 epochs. As we can see, Subfigures (b) and (c) are very similar, indicating that the subsequent 1200000 epochs yield minimal improvement for the discriminator network. Additionally, the discriminator's prediction shows larger inaccuracies on the boundary compared to the interior points. This phenomenon may be related to the limited number of training points assigned to each boundary, with only 8 points on each side.

In contrast, the solver network converges much more slowly, as evident in Figure <ref>. The solver network continuously reduces the error between its prediction and the exact solution, ultimately resulting in a more accurate solution than the discriminator network. However, it's important to note that while the solver network requires more time to converge than the discriminator, it is still far more efficient than simply using a traditional penalty formulation with a large penalty parameter, which will still be far from the actual solution after 1.5 million epochs. From Figure <ref>, we observe that after 500000 epochs, the prediction still deviates from the exact solution by a certain margin, but the result after 1000000 epochs is already quite good. The maximum error shown in Subfigure (b) is just slightly larger than in Subfigure (c) while the latter one took an additional 500000 epochs to achieve. This observation suggests a practical tradeoff between efficiency and accuracy that should be considered in real-world applications.









§ ACKNOWLEDGMENT



 §.§ Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the author(s) used ChatGPT in order to improve the writing of this paper. After using this tool, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.

abbrv

