Finite Element Method for Constrained Optimal Control Problems Governed by Nonlinear Elliptic Pdes

In this paper, we study the finite element method for constrained optimal control problems governed by nonlinear elliptic PDEs. Instead of the standard error estimates under L 2-or H 1-norm, we apply the goal-oriented error estimates in order to avoid the difficulties which are generated by the nonsmoothness of the problem. We derive the a priori error estimates of the goal function, and the error bound is O(h 2), which is the same as one for some well known quadratic optimal control problems governed by linear elliptic PDEs. Moreover, two kinds of practical algorithms are introduced to solve the underlying problem. Numerical experiments are provided to confirm our theoretical results.

Note that the state equation (2) is nonlinear, the theoretical analysis and numerical methods of the problem (1)- (2) are difficult that lead to only few results in the literatures. As far as we know there are some papers dealing with the numerical approximation of the problem by Chen et al. [6,7]. In these papers, first order conditions of the optimal control problem are discussed, and the smoothing method and semismooth method are introduced to solve the underlying problem.
Instead of finite difference methods used by Chen, in this paper we consider finite element methods to approximate the optimal control problem (1)- (2), and discuss the related error estimates. As mentioned in [6], since our problem is nonsmooth and globally nonconvex, the solutions are not unique in many cases and the first order condition is neither sufficient nor necessary for the original problem. Therefore the standard error analysis of optimal control problems (see, e.g., [9], [11], [16], [18]) can not be applied in this kind of problem. Moreover, since the solutions of the problem are not unique in many cases, the standard error estimates under H 1 −norm or L 2 −norm are not proper because that we can not point out to which specific solution (y, u) is approximated by finite element solution (y h , u h ). To overcome these difficulties, we resort to goal-oriented error estimates with the form |J(y, u) − J h (y h , u h )|. Similar idea has been used for a posteriori error estimates (see, e.g., [2]). In our paper, we prove that J h (y h , u h ) is convergent to J(y, u), and the convergence order is O(h 2 ), which is the same as one for some well known quadratic optimal control problems governed by linear elliptic PDEs. It indicates that the finite element solution tends to an exact solution in the sense of goal function.
The outline of the paper is as follows: In the next section, we will introduce the weak form of the problem and provide its finite element scheme. In section 3, we discuss the a priori goal-oriented error estimates of the finite element approximation for the problem. In section 4, we introduce smoothing SQP algorithm ( [7]) and projection gradient method for solving the underlying problem. Moreover, some numerical examples are presented in this section to confirm our theoretical results.
Let us consider the finite element approximation to the problem (4)- (5). We consider only n-simplex elements and conforming finite element spaces. Let Ω h be a polygonal approximation to Ω with a boundary ∂Ω h . To be simple, we assume that Ω is a polygon such that Ω h = Ω. Let T h be a regular partitioning of Ω into disjoint n-simplices τ , so thatΩ = τ ∈T hτ ,τ andτ have either only one common vertex or a whole edge or face if τ, τ ∈ T h andτ ∩τ = φ. Let h τ denote the diameter of τ , h = max τ ∈T h h τ .
Associated with T h are finite dimensional subspaces V h of C(Ω) and U h of L 2 (Ω), such that for all τ ∈ T h , χ| τ are polynomials of k-order (k ≥ 1), ∀χ ∈ V h , and ψ| τ are polynomials of m-order ( For simplicity, we will only consider the conforming piecewise linear finite element space, i.e., k = 1 for V h , and the piecewise constant finite element space, i.e., m = 0 for U h . Then the standard FEM approximation to (4)-(5) is 3. Goal-oriented error estimate. In this section, we will focus on the goaloriented error estimates of the finite element approximation, i.e., |J(y, u)−J h (y h , u h )|, where the functionals J(·, ·) and J h (·, ·) are defined in (4) and (6), while (y, u) and (y h , u h ) are one of the solutions of the problems (4)-(5) and (6)-(7), respectively. In order to discuss the goal-oriented error estimates, we firstly introduce some results of the nonlinear elliptic equation (5), which have been proved in [5].
where f ∈ L 2 (Ω) and Ω is an open, bounded convex domain with a Lipschitz boundary. Let w h be the finite element approximation to w, i.e., Then the solutions of (8) and (9) are existent and unique, w ∈ H 2 (Ω) ∩ H 1 0 (Ω), and where the constant c is independent of h.
Using Lemma 3.1, we can obtain the following result: Let (y, u) and (y h , u h ) be one of the solutions of the problems (4)- (5) and (6)-(7), respectively. Then we have Proof. We introduce an auxiliary problem It is apparent that y h is the standard finite element approximation of y(u h ). Note that u h ∈ L 2 (Ω) and Ω is convex. Thus y(u h ) ∈ H 2 (Ω). It can be deduced from Lemma 3.1 that Note that (y, u) is one of the solutions of the problem (4)- (5), and y(u h ) satisfies the equation (11). Thus it is easy to see that Then it follows from (12) and (13) that Here the constant c in (12) and (14) comes from Lemma 3.1, and has been proven to be independent of h in [5]. Then the proof of this Lemma is completed.
Next we estimate the goal-oriented error on the other side of inequality.
Lemma 3.3. Let (y, u) and (y h , u h ) be one of the solutions of the problems (4)- (5) and (6)-(7), respectively. Assume that u ∈ H 1 (Ω), then we have It is well-known that if u ∈ H 1 (Ω), then (see e.g., [16]) We resort to two auxiliary equations: Setting (17) and (18), we have Note that on any point in the domain Ω, the function Thus, inequalities (19) and (20) imply that It is easy to see that where Q h v is the L 2 −projection of v defined as (15). Thus, it follows from (21) and (22) that (23) On the other hand, note that y h (u) is the finite element approximation of y. Using Lemma 3.1, we have that (24) Resorted to some basic properties of projection operator, it follows from (16), (23) and (24) that Since (y h , u h ) is one of the solutions of the problem (6)-(7), then Therefore, it follows from (25) and (26) that which completes the proof.
Summarizing the results of Lemmas 3.2 and 3.3, we have the following goaloriented error estimate.
Theorem 3.4. Let (y, u) and (y h , u h ) be one of the solutions of the problems (4)- (5) and (6)-(7), respectively. Assume that u ∈ H 1 (Ω), then Remark 1. In Theorem 3.4, a goal-oriented error estimate is obtained under the condition u ∈ H 1 (Ω). This condition can be satisfied in many practical cases. If only the weaker condition can be satisfied, e.g., u ∈ H γ (Ω), 0 < γ < 1, the project error estimate (16) and (22) can be replaced by Then it can be deduced similarly as the proof in Lemmas 3.2 and 3.3 that Consider the worst case, i.e., u only belongs to L 2 (Ω). Define Π h : Then it is easy to see that Π h 0,Ω is uniformly bounded. On the other hand, lim h→0 Π h u = 0 for any u ∈ H 1 (Ω). Note that H 1 (Ω) is dense in L 2 (Ω). We have that for all u ∈ L 2 (Ω), Then it can be deduced as Lemmas 3.2 and 3.3 that This means that the finite element approximation is always convergent to the exact solution in the sense of goal functional, this result is valid even in the worst case: u ∈ L 2 (Ω).

Remark 2.
The results of Theorem 3.4 and Remark 1 also can be extended to more general nonsmooth constrained optimal control problems. For example, J(y, u) can be extended to J(y, u) = g(y) + j(u) satisfying the following conditions: ,Ω . The state equation (2) can be extended to general nonlinear partial differential equations Ay(u) = f + u satisfying the following conditions: ,Ω , where y h (u) is the finite element approximation of y(u).

Remark 3.
If the functional J further satisfies that there is δ > 0 such that where u is the solution of (4)-(5) and B(u; δ) = {v ∈ K : v − u 0,Ω ≤ δ}, which is a weaker condition than locally convex condition, then from Theorem 3.4 and (12) we can deduce that where y(u h ) is the solution of (13). Thus we can obtain the error estimate: which is same with the standard error estimate for the constrained quadratic optimal control problems governed by linear elliptic PDEs. Next, we present an example which satisfies the condition (27). Let q(x) < 0 in Ω, where q(x) is defined in (3). For any v ∈ U ad , we have v(x) < 0 in Ω. Since (y(v), v) is the solution of (5), we have According to the principle of the maximum, the maximum of y(v) can be reached only on the boundary ∂Ω. Since y(v) = 0 on ∂Ω, we have y(v) ≤ 0 in Ω. Therefore the nonlinear term of (5) vanishes, and the control problem (4)-(5) can be treated as a linear constraint problem. Then we have the first optimal condition Such conclusion can be found in [15](see, Chapt. 2, Sec. 1.3). Then This means that the condition (27) is satisfied.

Numerical algorithms.
In the first place, we provide a practical algorithm called smoothing SQP method(SSQP), which has been introduced by Chen in [7].
Here we use the function Φ (y h ) = 1 2 (y h + y 2 h + 2 ) to approximate the nonlinear function max(0, y h ), choosing small enough. Moreover, in the (k + 1)-th iteration, is the solution of the k-th iteration. Therefore, the problem (6)- (7) in the (k + 1)-th is transformed into This is a linearized problem. We use a preconditioned projection algorithm proposed in [16](see, Chapt. 8, Sec. 8.2) to solve the discretized optimization problem (29)-(30), which is shown to be efficient and converge within a few iteration steps by many numerical experiments for problems with parameter α > 1.0e − 3. For small α, e.g., α < 1.0e − 3, we resort to the primal-dual active set method proposed in [3] to solve the problem (29)-(30). Then we can set the solution of (29)-(30) as (y k+1 h , u k+1 h ) and go to the next iteration. We stop the iteration when where σ is a given small positive number.
Now we restate the algorithm as follows.
Smoothing SQP method 1. Initialization: choose initial value u 0 h ∈ K h and calculate y 0 h by state equation, set iteration termination tolerance σ, smoothing parameter and set k = 0.
2. Solve the following sub-problem by the projection gradient method h ) be the solution of the above problem.
In [7], Chen has pointed out that this algorithm is convergent and effective when certain conditions are satisfied.
Next, we provide another algorithm without smoothing parameter. We introduce function φ(y h ) = 1 2 (sgn(y h )+1) as one of subgradients of max(0, y h ) (see [8]). Here, we still resort to the projection gradient method in which the adjoint state equation is needed in order to obtain the descent direction of u h (see [11] and [16]). Note that if we linearize the state equation as where p h ∈ V h 0 . Therefore in the k-th iteration, we solve the y k h directly by the equation Chang has provided two methods to solve the above equation in [5]. And the adjoint state p k h can be solved by We restate the algorithm as follows.
Projection gradient method with subgradient 1. Initialization: choose initial value u 0 h ∈ K h , set iteration termination tolerance σ and set k = 0.

Solve the state equation
, and then the adjoint state equation 4. Else, update k = k + 1 and goto 2.
In this algorithm, we use the linearized form only to obtain the adjoint state equation. The state equation itself is not linearized in computation. Since the smoothing parameter is not used here, we can treat this algorithm as the first algorithm with approaching to zero. Numerical tests in the following show that both algorithms are convergent to the solution well.
Remark 4. It is necessary to point out that although Chen in [7] has provided a convergent condition of the SSQP method, it is not a condition easy to check. There are more methods to solve the nonsmooth optimal problem. They are almost related with the concept subgradient (see [8]), such as Proximal Bundle Method (see in [10], [17]) and Gradient Sampling algorithm (see in [4]). However, to our best knowledge, they are focus on some specific problems and none of them can provide a convergent condition easy to be checked.
In the meantime, since the nonlinear part is not a dominated part of the partial differential equation in our problem, the simple linearized method, such as SSQP method and projection gradient method are effective. In fact, we have chosen different subgradients in the second algorithm which is the main idea of the modified methods such as PBM or GS, the numerical results are almost the same due to the weak nonsmoothness of the problem.

Numerical examples.
In this subsection we will present some numerical results obtained for the optimal control problem (1)-(2) by using the discretization scheme described in section 2. The purpose is to clarify the convergence rates to be expected and the efficiency of the algorithms. For the computation the software package AFEPack ( [14]) has been used.
In the following Tables, we report the functional error |J(y, u) − J h (y h , u h )| for different values of "Node" and , where "Node" represents the number of nodes of the mesh, and "order" represents the error orders. The "reference solutions" are calculated on a fine mesh with more than 8 × 10 5 nodes and = 10 −6 to replace the exact solutions, because that we can not have the exact solutions in these numerical examples. Since we don't have the well known standard to choose the smoothing parameter in SSQP method we provide the results with different values to see how it effects the final result.
Firstly, we use the SSQP method. It is shown from Table 4.1 and Table 4.2 that the order of goal-oriented error |J(y, u) − J h (y h , u h )| is O(h 2 ) when is small enough, which is consistent with our theoretical results. The tables show that the order of convergence is not good enough as expected by theoretical analysis when = 10 −2 in both examples. Generally speaking, the choice of has little effect on the result. It is because that only effects the nonlinear term λ max(0, y), but this term is not the main term compared with the term −∆y in the state equation.  Secondly, we show the results computed by the project gradient method in Table  4.3 and Table 4.4. It is easy to see that the convergence order for goal-oriented error |J(y, u) − J h (y h , u h )| is also O(h 2 ) and the results are almost the same with the SSQP method when = 10 −6 .