Optimization Methods for Dirichlet Control Problems
Mariano Mateos









Optimization Methods for Dirichlet Control Problems


    
Mariano Mateos
Dpto. de Matemáticas,
Universidad de Oviedo, Campus de Gijón, 33203 Gijón, Spain 
Email: mmateos@uniovi.es.
    January 2017
===============================================================================================================================






We discuss several optimization procedures to solve finite element approximations of linear-quadratic Dirichlet optimal control problems governed by an elliptic partial differential equation posed on a 2D or 3D Lipschitz domain. The control is discretized explicitely using continuous piecewise linear approximations. Unconstrained, control-constrained, state-constrained and control-and-state constrained problems are analyzed. A preconditioned conjugate method for a reduced problem in the control variable is proposed to solve the unconstrained problem, whereas semismooth Newton methods are discussed for the solution of constrained problems. State constraints are treated via a Moreau-Yosida penalization. Convergence is studied for both the continuous problems and the finite dimensional approximations. In the finite dimensional case, we are able to show convergence of the optimization procedures even in the absence of Tikhonov regularization parameter. Computational aspects are also treated and several numerical examples are included to illustrate the theoretical results.



Dirichlet optimal control, discretization, constrained optimization, preconditionate conjugate gradient, semismooth Newton methods




 49J20, 
49M05, 
49M15,  
65N30 



myheadings
plain
M. Mateos Optimization for Dirichlet Control Problems




§ INTRODUCTION


Only in the last ten years it has been possible to develop a systematic study of Dirichlet optimal control problems governed by elliptic equations, which started with the seminal paper <cit.>. In that work a 2D control-constrained problem governed by a semilinear equation posed on a convex polygonal domain is studied. Several other works have been published about numerical error estimates; see <cit.> or <cit.>  for a variational approach to control-constrained problems posed on smooth 2D or 3D domains, <cit.> for a superconvergence result on the state approximation for unconstrained 2D problems, <cit.> for control-and-state constrained problems posed on convex polygonal domains. In <cit.> the authors study the control-constrained problem in a 2D smooth convex domain taking into account the problems derived by the boundary approximation and in <cit.> an apparent paradox between <cit.> and <cit.> is explained. The regularity of the solution in possibly nonconvex polygonal plane domains is studied in <cit.>; see also the introduction of that paper for further references about related problems. In the recent publication <cit.>, error estimates for a Dirichlet control problem governed by a parabolic equation are obtained. In this work, the spatial discretization of the control is studied in both the cases of continuous piecewise linear finite elements and variational approach.



Just before that, several works dealing with efficient optimization methods for control or state constrained problems had appeared; in <cit.> the semismooth Newton method is thoroughly studied. Since all these papers about optimization are previous to the papers about the numerical analysis of Dirichlet control problems, very little or no reference is made in them to their applicability to the problems we are going to deal with. Only in <cit.> two different Dirichlet control problems with more regular solutions than the ones we treat here are studied in the infinite-dimensional case. Let us also mention that in <cit.> the Dirichlet boundary condition is replaced by a Robin penalization. For this kind of penalization the methods developed in the aforementioned references are directly applicable. In <cit.>, the semismooth Newton method is studied in the context of the variational approach of the control. Although the authors only exemplify their results through distributed and Robin boundary control, a combination of their results and some of the results we present in Section <ref> can be applied to Dirichlet control. See Remark <ref> below.

In this work we describe optimization methods for Dirichlet control problems in both the infinite and finite dimensional cases. Convergence proofs, examples and practical implementation details are discussed for all the algorithms through the paper.


In Section <ref> we state the problem and prove that it is well posed in Lipschitz domains; see Lemma <ref>. So far, only smooth, convex, polygonal or polyhedral domains had been studied.

Next, we discretize the problem. As is usual in control problems, we have three choices to discretize the control: the first option is not to discretize the control, using a variational discretization as introduced in <cit.> for distributed problems; as a second option we may use piecewise constant functions; and finally we can use continuous piecewise linear functions.

 The choice is not trivial because first order optimality conditions for Dirichlet control problems involve the normal derivative of the adjoint state.

 If we discretize directly the optimality system using e.g. continuous piecewise linear finite elements, we would obtain two different approximations of the control: the trace of the discrete state would be continuous piecewise linear and the normal derivative of the discrete adjoint state would be piecewise constant.
In <cit.> and in <cit.> this difficulty is solved using a mixed formulation of the state equation which is discretized with the lowest order Raviart-Thomas element. In <cit.> a convergence proof for this kind of discretization is given. In this way, both the trace of the discrete state and the normal derivative of the discrete adjoint state are piecewise constant functions, so the identification of both with the discrete control is meaningful. The discretization of the control in the presence of control constraints is carried out in these papers using a variational discretization. This kind of approximation may be convenient when the gradient of the state is the variable of interest, since this quantity is treated as one of the variables of the problem.


Another approach, in finite dimension, is to use the variational discrete normal derivative of the discrete adjoint state introduced in <cit.>, which is a continuous piecewise linear function. Doing so, both the trace of the discrete state and the normal derivative of the discrete adjoint state are continuous piecewise linear functions, so the identification of both with the discrete control is meaningful. Following this idea, in <cit.> the control is not discretized for the control-constrained case, but a variational approach is followed.

In this work we investigate optimization methods for the case of discretizing the control explicitely using continuous piecewise linear functions.

The end of Section <ref> is devoted to describe with some detail some computational aspects that will be important in the rest of the work.

We next provide in Section <ref> an efficient method to solve the unconstrained problem.  We propose in Subsection <ref> a preconditioned conjugate gradient (pcg in the rest of the work) method for a reduced problem in the spirit of <cit.>. We are able to prove convergence of the conjugate gradient even for the case where the Tikhonov regularization parameter ν vanishes.

In sections <ref>, <ref> and <ref> we study the convergence of the semismooth Newton method for the constrained problems and write practical algorithms for the solution of the finite dimensional approximations.
In Section <ref> we deal with the control-constrained problem. In many of the aforementioned references about the semismooth Newton method for control-constrained problems the authors study convergence of an abstract problem or an infinite dimensional problem. For distributed and Neumann control problems the results are immediately applicable to the finite dimensional approximations because the controls can be discretized using piecewise constant functions, and therefore the variational inequality arising from first order necessary optimality conditions can be written in an element-wise form. The same idea applies when we are dealing with a variational approach as proposed in <cit.> or a mixed formulation, as studied in <cit.>. 
When we use continuous piecewise linear elements,
 the variational inequality cannot be written in a point-wise or elementwise form; see  (<ref>) and Remark <ref>. We include the analysis of Newton methods for both the original problem and the discretized one.

In Section <ref> we study the state-constrained problem using a Moreau-Yosida penalization. Since there are no control constraints, the analysis of the semismooth Newton method for the infinite dimensional problem is applicable to the finite dimensional one, so we do not need to repeat it. We prove that the finite-element approximation of the penalized problems converge to the solution of the penalized problem. This result cannot be deduced straightforward from the ones in the literature since the penalized functional is not of class C^2. A continuation strategy as proposed in <cit.> is developed. Finally, in Section <ref> we discuss the problem with both control and state constraints.

It is well known that the main difficulty with Dirichlet control problems is the low regularity of the solutions. This regularity and related error estimates are, so far, well established in 2D polygonal domains <cit.> and 2D or 3D smooth domains <cit.> but there is not, up to our best knowledge, a general study in 3D polyhedral domains. Although the main focus of this work is on optimization methods, we also study the regularity of the solution and error estimates of the approximations in one example case in a polyhedron; see Example <ref>.




§ STATEMENT OF THE PROBLEM AND FIRST PROPERTIES


Let Ω⊂ℝ^d, d=2 or d=3, be a bounded domain with  Lipschitz boundary Γ  and in this domain consider a target state y_Ω∈ L^2(Ω). Consider also the continuous linear operator S:L^2(Γ)→ L^2(Ω) such that y=Su if and only if y is the solution in the transposition sense (see Definition <ref> below)  of

    -Δ y = 0Ω, y=uΓ.

Let ω⊂Ω be a domain such that ω̅⊂Ω and define two pairs of functions α, β∈ C(Γ) and a, b∈ C(ω̅) such that α(x)<β(x) for all x∈Γ and a(x)<b(x) for all x∈ω̅. For some fixed regularization parameter ν≥0, define

    J(u) =1/2Su-y_Ω^2_L^2(Ω)+ν/2u^2_L^2(Γ)

and consider the sets

    U_α,β={u∈ L^2(Γ): α(x)≤ u(x)≤β(x)x∈Γ},

and

    K_a,b={y∈ L^2(Ω)∩ C(ω̅): a(x)≤ y(x)≤β(x)x∈ω̅}.


In this work we will study optimization procedures for

the following four Dirichlet control problems:

    (P^U)min_u∈ L^2(Γ) J(u),   (P^C)min_u∈ U_α,βJ(u),   (P^S)min_Su∈ K_a,bJ(u),   (P^CS)min_[  u∈ U_α,β; Su∈ K_a,b ]J(u),

namely the unconstrained, control-constrained, state-constrained and control-and-state constrained problems.
Almost all the literature related to these problems is written using the state equation (<ref>). There would be no problem in taking into account an equation of the form

    Ay :=-∑_i=1^d∂_i(a_i,j∂_j y) + a_0 y = FΩ,   y = u + GΓ

with regular enough F, G,  a_0≥ 0 and a_i,j=a_j,i  satysfing an uniform ellipticity condition.



Let us state precisely what we mean by solution in the transposition sense. Consider the space

    Φ = {ϕ: ϕ∈ H^1_0(Ω)Δϕ∈ L^2(Ω)}.

This is a Banach space with the graph norm ϕ_Φ = ϕ_H^1(Ω)+ Δϕ_L^2(Ω).
Further, the functions in this space satisfy ∂_n ϕ∈ L^2(Γ). This is known to be true for smooth domains; convex domains, see <cit.>; plane polygonal domains, see <cit.>; or polyhedral domains, see <cit.> and the usual trace theorem. We have not been able to find a proof of this fact for general Lipschitz domains.  In <cit.> the regularity ϕ∈ H^3/2(Ω) is proved. Nevertheless, as the authors notice in page 165 of this reference, the trace theorem is not valid neither in H^3/2(Ω) nor in H^1/2(Ω), so we cannot deduce immediately that ∂_nϕ∈ L^2(Γ). The results in <cit.> imply that the usual trace result can be extended to harmonic functions in the limit cases. We show  next how to take advantage of this  to prove that ∂_nϕ∈ L^2(Γ) in Lipschitz domains. Regarding the analysis of semismooth Newton methods –see Lemma <ref> below– we also prove L^q(Γ) regularity for some q>2.
Let Ω⊂ℝ^d, d=2 or d=3, be a bounded domain with  Lipschitz boundary Γ and consider ϕ∈Φ. Then, there exists q_0>2 depending on the domain such that, for 2≤ q<q_0, we have
∂_n ϕ∈ L^q(Γ) and

    ∂_n ϕ_L^q(Γ)≤ C Δϕ_L^2(Ω).


If, further, Ω is smooth or convex or polygonal or polyhedral, then there exists t>0 such that ∂_n ϕ∈ H^t(Γ) and

    ∂_n ϕ_H^t(Γ)≤ C Δϕ_L^2(Ω).




Denote z=-Δϕ, extend z by zero to ℝ^d, consider the Newtonian potential centered at the origin N(x) and define w = z * N, the convolution product of z and N. Then w∈ H^2(Ω) and ∇ w∈ H^1(Ω)^d, so it is clear that Tr (∇ w)∈ H^1/2(Γ)^d↪ L^q(Γ)^d for all q<∞ if d=2 and all q≤ 4 if d=3, since the dimension of Γ is d-1. This implies that:

(a) ∂_n w = ∇ w· n∈ L^q(Γ) because the boundary Γ is Lipschitz and therefore it has a unit normal vector defined almost everywhere; and there exists C>0 such that

    ∂_n w_L^q(Γ)≤ CΔϕ_L^2(Ω);


(b) g = Tr(w)∈ W^1,q(Γ) due precisely to the definition of W^1,q(Γ), and there exists C>0 such that

    g_W^1,q(Γ)≤ C Δϕ_L^2(Ω).

Define now v∈ H^1(Ω) the unique solution of -Δ v = 0 in Ω, v = g on Γ.
Using <cit.>, we have that there exists q_0>2 such that if 2≤ q <q_0, then the nontangential maximal function of the gradient of v satisfies M(∇ v)∈ L^q(Γ) and there exists C>0 such that

    M(∇ v)_L^q(Γ)≤ Cg_W^1,q(Γ).

As is pointed out in <cit.>, this implies that ∇ v has nontangential limit a.e. on Γ and we can define the normal derivative of v at a point s∈Γ  as the nontangential limit as x→ s of ∇ v(x)· n(s). For a precise definition of nontangential limit and nontangential maximal function see, e.g., the introduction of the work <cit.>. This, together with inequalities (<ref>) and (<ref>) imply that ∂_n v ∈ L^q(Ω) and

    ∂_n v_L^q(Γ)≤ C Δϕ_L^2(Ω).


So we have that ϕ = v-w and we can define in a natural way ∂_n ϕ  = ∂_n v -∂_n w∈ L^q(Γ). The estimate (<ref>) follows from (<ref>) and (<ref>).

For smooth, convex, polygonal or polyhedral domains, the second result follows from the regularity ϕ∈ H^3/2+t(Ω) and the usual trace theorem for ∇ϕ. See <cit.> for convex domains,  <cit.> for plane polygonal domains and <cit.> for polyhedral domains.



We will say that y∈ L^2(Ω) is the solution in the transposition sense of (<ref>) if

    (y,-Δϕ)_Ω = -(u,∂_n ϕ)_Γϕ∈Φ.


Here and in the rest of the work (·,·)_X  stands for the standard inner product in L^2(X).


The adjoint operator of S is S^*:L^2(Ω)→ L^2(Γ) defined by S^*z = -∂_nϕ, where ϕ is the unique weak solution of

    -Δϕ = zΩ, ϕ = 0Γ.

We can write now

    J(u)    =   1/2(Su-y_Ω,Su-y_Ω)_Ω+ν/2(u,u)_Γ =1/2(S^*Su+ν u,u)_Γ - (S^*y_Ω,u)_Γ + c_Ω

where c_Ω=0.5 y_Ω_L^2(Ω)^2 is a constant.

From this expression, we can easily compute the derivative of J at a point u∈ L^2(Γ) in the direction v∈ L^2(Γ):

    J'(u)v=(S^*Su+ν u,v)_Γ -(S^*y_Ω,v)_Γ.

For later use, we will define now for every u∈ L^2(Γ), y_u=Su∈ H^1/2(Γ) and φ_u∈ H^1_0( Ω) the unique solution of

    -Δφ_u = y_u-y_ΩΩ, φ_u=0Γ.








 §.§ Discretization.



To discretize the problems. consider  {𝒯_h}_h a regular family of triangulations of Ω̅. To simplify the notation, we will suppose that Γ is polygonal or polyhedral. Related to the mesh, we will call N the number of nodes and {x_j}_j=1^N the nodes of the mesh. We define the sets of interior indexes, boundary indexes and indexes in ω̅ as I = {j:x_j∈Ω}, B={j:x_j∈Γ} and J = {j:x_j∈ω̅}.
For the discretization of the state and the adjoint state we use the space of linear finite elements Y_h⊂ H^1(Ω),

    Y_h={y_h∈ C(Ω̅) y_h∈ P^1(T) ∀ T∈𝒯_h}.

As usual, we will abbreviate Y_h0=Y_h∩ H^1_0(Ω). For the control we use the space U_h of continuous piecewise linear functions on Γ

    U_h = {u_h∈ C(Γ) u_h∈ P^1(T∩Γ) ∀ T∈𝒯_h}.

Notice that the elements of U_h are the traces of the
elements of Y_h.
We will denote I_h:C(Ω̅)→ Y_h or I_h:C(Γ)→ U_h the interpolation operator related to these spaces and Π_h:L^2(Ω)→ Y_h or Π_h:L^2(Γ)→ U_h the projection onto this spaces in the L^2 sense. For all y∈ L^2(Ω) and all u∈ L^2(Γ):

    (Π_h y,y_h)_Ω = (y,y_h)_Ω ∀ y_h∈ Y_h
    (Π_h u,u_h)_Γ = (u,u_h)_Γ ∀ u_h∈ U_h.



We discretize the state equation following the work by Berggren <cit.>: define S_h:L^2(Γ)→ L^2(Ω) such that for u∈ L^2(Γ), y_h=S_hu if and only if y_h∈ Y_h is  the unique solution of

    a(y_h, z_h)=0 ∀ z_h∈ Y_h0, (y_h,v_h)_Γ = (u,v_h)_Γ ∀ v_h∈ U_h,

where a(·,·) is the bilinear form associated to the operator in the PDE. In the case of the Laplace operator a(y,z)=∫_Ω∇^T y∇ z dx.
The discrete functional is thus defined as

    J_h(u)=1/2(S_hu-y_Ω,S_hu-y_Ω)_Ω+ν/2(u,u)_Γ.


We define now

    U^h_α,β = {u_h∈ U_h: α(x_j)≤ u_h(x_j)≤β(x_j) ∀ j∈B},

 and

    K^h_a,b = {y_h∈ Y_h: a(x_j)≤ y_h(x_j)≤ b(x_j) ∀ j∈J}.

The discrete problems reads thus as

    (P_h^U)min_u_h∈ U_h J_h(u_h),   (P_h^C)min_u∈ U^h_α,βJ_h(u_h),   (P_h^S)min_S_hu_h∈ K^h_a,bJ_h(u_h),   (P_h^CS)min_[    u_h∈ U^h_α,β; S_hu_h∈ K^h_a,b ]J_h(u_h).




The adjoint operator of S_h is given by the discrete variational normal derivative. See <cit.>. S_h^*:L^2(Ω)→ L^2(Γ) and w_h=S_h^*y if w_h∈ U_h satisfies

    (w_h,z_h)_Γ =  (y,z_h)_Ω - a(z_h,ϕ_h)z_h∈ Y_h,

where ϕ_h∈ Y_h0 is the unique solution of

    a(z_h, ϕ_h) = (y,z_h)_Ω z_h∈ Y_h0.

It is customary to write S_h^*y=-∂_n^hϕ_h. For u∈ L^2(Γ) and y∈ L^2(Ω), we have

    (S_hu,y)_Ω = (u,S^*_h y)_Γ.

We can then write

    J_h(u) =1/2(S^*_hS_hu+ν u,u)_Γ -(S_h^* y_Ω,u)_Γ+c_Ω.

The computation of the derivative of J_h at a point u∈ L^2(Γ) in the direction v∈ L^2(Γ) is then obtained with

    J'_h(u)v =(S_h^*S_hu+ν u,v)_Γ - (S^*_hy_Ω,v)_Γ.

Since our final objective is to optimize in U_h, let us see in more detail how to make the computations in this space.

Consider {e_j}_1^N the nodal basis in Y_h, where N is the dimension of Y_h and satisfies N=N_I+N_B, the latter being respectively the number of interior and boundary nodes. With an abuse of notation, we will also denote e_j the restriction of e_j to Γ. Define the usual finite element stress, mass and boundary mass matrices by

    K_i,j=a( e_i,e_j), M_i,j=(e_i,e_j)_Ω, B_i,j =(e_i,e_j)_Γ1≤ i,j≤ N.



We will also use I∈ℝ^N× N for the identity matrix, O∈ℝ^N× N for the zero matrix and usually refer to submatrices as K_I,I or K_I,: defined by taking the rows or columns designed by the sets of indexes in the subscripts, the semicolon meaning “all the indexes”. For instace, the usual boundary mass matrix is B_B,B.

Given u_h=∑_j∈Bu_j e_j, we denote u∈ℝ^N_B× 1 the vector (u_1,…,u_N_B)^T and for y_h=∑_j=1^N y_je_j we denote y∈ℝ^N× 1 the vector (y_1,…,y_N)^T. Using (<ref>), we have that y_h=S_hu_h iff

    [[ K_I,I K_I,B; O_B,I B_B,B ]]
    [
    [ y_I; y_B ]]
    =
    [
    [      0; B_B,Bu ]].

Since B_B,B is nonsingular, we can write this as

    [ K_I,Iy_I        = -K_I,Bu,;      y_B        =       u. ]

If we define S∈ℝ^N_B× N as

    S=[[ K_I,I K_I,B; O_B,I I_B,B ]]^-1I_:,B

we have that y_h=S_hu_h if and only if y=Su.

Given y_h∈ Y_h, let ϕ_h=∑_j∈Iϕ_j e_j be the solution of (<ref>) for y=y_h. Denoting ϕ∈ℝ^N_I× 1 the corresponding vector, it can be computed as the solution of

    K_I,Iϕ = M_I,:y.

Then we could compute w∈ℝ^N_B× 1, the vector whose components are the coefficients of the (minus) discrete normal derivative  -∂_n^hϕ_h = S_h^*y_h=w_h=∑_j∈Bw_je_j solving the system

    B_B,Bw = M_B,:y -
    K_B,Iϕ.

Formally, we can also write that
w = B_B,B^-1S^TMy.
To finish this section we also define the matrix A∈ℝ^N_B,×  N_B and the vector f∈ℝ^N_B× 1 by

    A_i,j= (S_h^*S_h e_i+ν e_i,e_j)_Γ f_i= (S^*_hy_Ω,e_i)_Γ.

We have that

    J_h(u_h)=1/2u^TAu - f^Tu+c_Ω.

We also note that

    A = S^TMS+νB_BB.

To compute the vector f,  we consider the projection of y_Ω onto Y_h in the L^2(Ω) sense, y_Ω,h = Π_h y_Ω and denote y_Ω∈ℝ^N× 1, the vector whose j-th component is y_Ω,h(x_j), and 
    f=S^TMy_Ω.


So for u_h,v_h∈ U_h, the latter represented by the vector v, we can compute

    J'(u_h)v_h=v^TB_B,Bw +νv^T B_B,Bu-v^Tf.

Notice that applying the equality (<ref>), the explicit computation of w is not necessary, and we can write

    J'(u_h)v_h=v^T(M_B,:y-
    K_B,Iϕ +νB_B,Bu-f).





§ UNCONSTRAINED PROBLEM


Problem (P^U) has a unique solution u̅∈ L^2(Γ) that satisfies J'(u̅)=0.
For every 0<h<h_0, problem (P^U_h) has also a unique solution u̅_h that satisfies J'_h(u̅_h)=0. The problems being convex, these conditions are also sufficient. Moreover it is known that u̅_h→u̅ strongly in L^2(Γ) and also error estimates are available in terms of the mesh size h if ν>0 and the domain is either 2D and polygonal or smooth. See <cit.>.

Let us comment two different approaches for the solution of the FEM approximation of (P^U).











 §.§ Solve the optimality system for the state, the adjoint state and the control.

We write first order optimality conditions for (P^U_h). There exists some h_0>0 such that for every 0<h<h_0, there exist unique u̅_h∈ U_h, y̅_h∈ Y_h and φ̅_h∈ Y_h0 satisfying the optimality system:


    a(y̅_h,z_h)   =   0 z_h∈ Y_h0,
    y̅_h   ≡   u̅_hΓ,
    
    a( z_h, φ̅_h)     =   (y̅_h-y_Ω,z_h)_Ωz_h∈ Y_h0,
    νu̅_h    ≡   ∂_n^hφ̅_hΓ.


Taking into account the definition of discrete variational normal derivative and relations (<ref>)–(<ref>), we can write this optimality system as

    [[  K_I,I  O_I,B  K_I,B  O_I,I;  O_B,I  I_B,B -I_B,B  O_B,I; -M_I,I -M_I,B  O_I,B  K_I,I;  M_B,I  M_B,B νB_B,B -K_B,I ]]
    [[ y_I; y_B;   u; φ_I ]] =
    [[       0_I;       0_B; -M_I,:y_Ω;  M_B,:y_Ω ]].

We may eliminate u with the boundary condition, and reorder the equations  to get the linear system

    [[   M+νB -K_:,I; -K_I,:  O_I,I ]]
    [[   y; φ_I ]] =
    [[ My_Ω;    0 ]].

Notice that system (<ref>) is solvable even for ν = 0.
Solving this equation would completely solve the problem for the unconstrained case.
When the discretization is very fine, the number of unknowns may make the solution of the system by direct methods too difficult.
The preconditioned conjugate gradient method for this kind of linear systems is studied by Schöberl and Zulehner in <cit.> and  Herzog and Sachs in <cit.>. A preconditioner can be built using matrices representing scalars products in Y_h and Y_h0. Following Algorithm 1 in the last-mentioned reference, at each iteration three linear systems must be solved: two of size N and one of size N_I. In <cit.>, the systems are solved using a multigrid method.

Reduced problems in the adjoint state variable and related pcg have also been studied in <cit.>. Nevertheless, the structure of the matrix we have in (<ref>) is different to the one treated in that reference and application of their results is not straightforward.



 §.§ Use an iterative method to solve a reduced problem in the control variable.

Let us see another way of solving (P^U_h). Using (<ref>), first order optimality conditions can also be written as

    Au = f.

The matrix A is symmetric and positive definite for all ν≥ 0 and there exists C>0 independent of h and ν such that its condition number is bounded by

    κ( A ) ≤
    C λ_N_B(B_B,B)/λ_1(M)λ_N(M) +νλ_N_B(B_B,B)/λ_1(M) +νλ_1(B_B,B),

where 0<λ_1(B_B,B)< λ_N_B(B_B,B) and 0<λ_1(M)<λ_N(M) are respectively the smallest and greatest eigenvalues of the matrices B_B,B and M.

It is clear from (<ref>) that A is symmetric. The mass matrices M and B_B,B are symmetric and positive definite and therefore λ_1(M)>0 and λ_1(B_B,B)>0.  Since the boundary components of Su are exactly u, we have that Su_ℝ^N≥u_ℝ^N_B and hence

    u^TAu = u^T(S^TMS+νB_B,B) u =
    (Su)^TM(Su) +νu^T B_B,Bu)≥ (λ_1(M)+ νλ_1(B_B,B)) u_ℝ^N_B^2

so A is positive definite.

From <cit.>, we know that there exists some C>0 such that S_h u_h_L^2(Ω)≤√(C)u_h_L^2(Γ). Since

    S_hu_h_L^2(Ω)^2 = (Su)^TM(Su)≥λ_1(M)Su_ℝ^N^2

and

    u_h_L^2(Γ)^2 = uB_B,Bu≤λ_N_B(B_B,B)u_ℝ^N_B^2,

we obtain

    u^TAu =
    (Su)^TM(Su) +νu^T B_B,Bu)
    ≤ (C λ_N_B(B_B,B)/λ_1(M)λ_N(M)+ νλ_N_B(B_B,B)) u_ℝ^N_B^2.

The bound for the condition number follows directly from the last inequality and (<ref>).

For a graded family of mesh with grading parameter 0<μ≤ 1 (see <cit.>), we usually define h=N^-1/d, N being the number of nodes of the mesh and d the dimension of Ω. The case μ=1 corresponds to that of a quasi-uniform mesh family.
If 𝒯_h is a family of graded meshes with grading parameter 0<μ≤ 1, then there exists C>0 independent of h and ν such that

    κ(A)≤ C h^(2d-1)(1-1/μ)1 + ν h^(1-1/μ)d/h^1/μ+ν.

In particular, for a quasi-uniform family, there exists C>0 such that

    κ(A)≤ C1+ν/ν+h.



 from <cit.> and <cit.>, we have that there exists constants 0<C_1<C_2

    λ_N_B(B_B,B)≤ C_2 h^d-1,λ_1(B_B,B)≥ C_1 h^d-1/μ,λ_N(M)≤ C_2 h^d,λ_1(M) ≥ C_1 h^d/μ

Estimate (<ref>) follows then from (<ref>). Estimate (<ref>) follows from (<ref>) for μ=1.


The explicit computation of the matrix A is out of the question, since it requires the solution of 2N_B systems of size N_I. Much better performance can be achieved using an iterative method which only requires the computation of d=Au. This is shown in Algorithm <ref>.



Preconditioned conjugate gradient methods can be carried out at the cost of one evaluation of Au per iteration. The computation of f can be done in a similar way, but of course it only must be computed once; see Algorithm <ref>.



To finish this section, let us say a word about a preconditioner to solve (<ref>). In practice, matrices like B_B,B or  P=M_B,B+νB_B,B make acceptable preconditioners. This is specially important when using graded meshes; see Example <ref> below.
Notice that at each pcg iterate, we only need to solve two linear systems, each of size N_I, to compute Au, plus another one of size N_B  to compute P^-1r.





All the code for all the examples has been made by the author using Matlab R2015b and has been run on a desktop PC with Windows 10 Pro 64bits with 16GB RAM DIMM 1333Mhz on an Intel Core i7 CP 870@2.93Ghz.


In this example we compare the different execution times, t_D and t_P, that we  obtain when we solve the optimality system using a direct solver for equation (<ref>) or using an iterative method, the preconditioned conjugate gradient method in our case, to solve the reduced problem (<ref>).
We have used respectively Matlab builtin commands  and .

We use the example in <cit.>, where Ω the interior of the convex hull of the points (-0.5,-0.5), (0.5,-0.5), (0.5,0), (0,0.5), (-0.5,0.5); y_Ω≡ 1 and ν=1.

A rough initial mesh is built with Matlab PDEToolbox (version 2.1) commands  and  set to 0.2 and  subsequent nested meshes are obtained by regular diadic refinement using . We use our own code to assemble the matrices K, M and B.

For the pcg we use as initial guess the null vector, a relative tolerance of 1E-10 and the preconditioner  P = M_B,B+νB_B,B. At each iteration we have to solve the linear systems (<ref>) and (<ref>). To do this, we first obtain the Cholesky factors of K_I,I so that at each iteration we only have to solve triangular systems. The interior mesh nodes were ordered using a symmetric approximate minimum degree permutation in order to minimize the number of nonzero elements in the Cholesky factors of K_I,I using Matlab command . An analogous reordering of the boundary nodes is done to optimize the sparsity pattern of the Cholesky factors of the preconditioner. This reorderings and factorizations take more or less the same cputime  than each single iteration of the pcg. This time is included in t_P.

For reference and possible double-check, we also report on the optimal value of the functional. We have marked with ∞ the experiments in which we have run out of RAM.

From the data in Table <ref> it is clear that the proposed iterative method for the reduced problem (<ref>) has an overall better performance for Dirichlet control problems.








We have also coded Algorithm 1 in <cit.>, though the practical comparison with the other two methods is less clear.


The main difficulty is the choice of appropriate scalar products. For instance, we have tried (using the notation in  <cit.>) X=K+B, X=K+M or  X=K+M+B for the scalar product in Y_h and Q=K_I,I  for the scalar product in Y_h0, but we have observed that the number of iterations needed to achieve the prescribed accuracy grows as the number of unknowns N increases. For the first three meshes in the 2D example described in Table <ref> and ν=1 we obtain 126, 173 and 236 iterations respectively.



In this example we show the effect of the use of P=M_B,B+νB_B,B as a preconditioner to solve (<ref>). In this case ν = 0.01, Ω = {x=re^iθ∈ℝ^2: r<1, 0<θ<11π/12}, and we use a mesh family graded at the origin with parameter μ = 13/33≈ 0.4; see <cit.>.
In Table <ref> we compare the number of iterations used to reach the prescribed accuracy of 1E-10 without and with preconditioner. This example shows very clearly the effectiveness of this strategy of preconditioning.





In this example we show that the number of pcg iterations for each fixed h is independent of the value of ν≥ 0 for all ν<ν_0. We take the same domain and meshes as in the 2D problem in Example <ref> and set y_Ω(x)=|x|^2, so that J(u̅)>0 even for ν=0. The results are summarized in Table <ref>.



A 3D example in a polyhedron.
Up to our best knowledge, there is not a general theory about regularity of the solutions or approximation results for 3D Dirichlet control problems posed on polyhedral domains. In <cit.>, the authors study a semi-discrete or variational approach to 3D problems on regular domains. Although the semi-discrete approach coincides with the full approach for unconstrained problems, we cannot profit their results since the regularity of the solutions in a smooth domain is (see <cit.>) much higher than the one we may obtain in polyhedral domains.

For this example we will take Ω the unit cube (-0.5,0.5)^3, y_Ω≡ 1 and ν = 1. First, we obtain a regularity result for the solution of our problem.

Let Ω⊂ℝ^3 be the interior of a rectangular parallelepiped and y_Ω∈ L^p(Ω) for some 3<p<+∞. Consider u̅ the solution of problem (P^U). Then, u̅∈ W^1-1/p(Γ), y̅∈ W^1,p(Ω) and φ̅∈ W^2,p(Ω). Moreover, u̅≡ 0 on the edges and corners of Γ.

Since u̅∈ L^2(Γ), using transposition and interpolation it is clear that y∈ H^1/2(Ω). Classical Sobolev embedding in 3D leads to y∈ L^3(Ω). Since Ω is a parallelepiped and y_Ω is regular enough, <cit.> states that φ̅∈ W^2,3(Ω), and hence ∂_nφ̅∈Π_i W^2/3,3(Γ_i), where Γ_i, i=1:6, are the faces of Γ. This does not imply immediately that u̅ belongs to W^2/3,3(Γ) because 2/3· 3=2, which is the topological dimension of Γ, and some integral compatibility condition should be imposed on the edges; but for all q<3, it is true that u̅∈ W^1-1/q,q(Γ), and therefore y̅∈ W^1,q(Ω). Using again Sobolev embeddings, we have that y̅∈ L^s(Ω) for all s<+∞ and hence y̅-y_Ω∈ L^p(Ω). Applying once again <cit.>, we have that φ̅∈ W^2,p(Ω).

Now we have that ∂_nφ̅∈Π_i W^1-1/p,p(Γ_i) and we can prove that if we define ∂_nφ̅=0 on the edges of Ω, then we obtain a continuous function.
To do this,
we use a similar argument to the one used in <cit.> for 2D problems. Since p>3, φ̅∈ C^1(Ω̅). Consider two faces A and B with a common edge AB and let τ^1_X, τ^2_X be two linearly independent vectors tangent to face X (X=A or X=B) such that τ^1_A=τ^1_B is also tangent to the edge AB. Since φ̅= 0 on Γ, we have that for every x∈ A, ∇φ̅(x)·τ^1_A = ∇φ̅(x)·τ^2_A = 0 and for every x∈ B,
∇φ̅(x)·τ^1_B = ∇φ̅(x)·τ^2_B = 0. So for x∈ AB, we have that ∇φ̅(x)=0  and therefore ∇φ̅(x)· n can be seen as a continuous function  if we set its value to 0 on the edges, despite the jump discontinuities of the normal vector n.

So u̅ is continuous and hence u̅∈ W^1-1/p,p(Γ). The regularity of the optimal state follows from the trace theorem; see e.g. <cit.>.

Using this regularity, an error estimate can be proved for the example problem.
Let Ω be a rectangular parallelepiped and y_Ω∈ L^p(Ω) for some 3<p<+∞. Consider u̅ be the solution of problem (P^U) and u̅_h be the solution of (P^U_h). Then there exists some h_0>0 and C>0 such that for all 0<h<h_0

    u̅-u̅_h_L^2(Γ)≤ C h^1-1/p.


The proof follows the same guidelines as those of <cit.> or <cit.> and thus will be omitted. An interesting remark is that for 2D problems, we can deduce uniform convergence for the controls using the L^2(Γ) estimate and an inverse inequality, since the boundary is 1D. This does not work for 3D problems with the error estimate at hand, since now Γ is 2D.


To solve the problem we have used a regular mesh of identical cubes of size h, each of them divided into 6 tetrahedra according to the Kuhn triangulation of the cube (see e.g. <cit.>). Up to our best knowledge, the current version of the Matlab PDEToolbox (version 2.1) computes an approximation of B using the barycenter formula. Although this does not affect the order of convergence for the FEM, this matrix plays a central role in Dirichlet control problems (for instance, B_B,B is not singular, but the barycenter approximation may be singular), so we have computed it in an exact way (using the mid-sides formula). Mesh data, computation times and optimal values are displayed in Table <ref>.






§ CONTROL CONSTRAINED PROBLEMS


Problem (P^C) has a unique solution u̅∈ L^2(Γ) that satisfies J'(u̅)(u-u̅)≥ 0 for all u∈ U_α,β.
For every 0<h<h_0, problem (P^C_h) has also a unique solution u̅_h that satisfies J'_h(u̅_h)(u_h-u̅_h)≥ 0 for all u_h∈ U^h_α,β. The problems being convex, these conditions are also sufficient. Moreover it is known that u̅_h→u̅ strongly in L^2(Γ) and also error estimates are available in terms of the mesh size h when ν>0 and the domain is convex in 2D. See <cit.>. In <cit.> the  smooth (2D and 3D) case is treated using variational approach, whose optimization process is different from the one we are presenting in the work at hand: in <cit.>, the problem is solved using a fixed point method, convergent for ν large enough; see <cit.> and Remark <ref> below for the convergence of the semismooth Newton method.


 §.§ Continuous problem

We can formulate first order optimality conditions as: there exist unique u̅∈ L^2(Γ), y̅∈ H^1/2(Ω) and φ̅∈ H^1_0(Ω) such that


    -Δy̅   =    0xxxxΩ, y̅=u̅Γ,
    
    -Δφ̅   =   y̅-y_ΩΩ, φ̅=0Γ,
    
    (-∂_nφ̅+ νu̅, u-u̅)   ≥    0xxxx u∈ U_α,β.


We first describe a semismooth Newton method to solve this optimality system and prove a convergence result for it (see Theorem <ref>). Next, we will reformulate the optimality system in terms of the Lagrange multipliers related to the constraints. We will see that this approach is better suited to the discrete problem, but has the drawback that the Newton method related to it is not semismooth; nonetheless, we will prove a convergence result for it; cf. Theorem <ref>.
To facilitate the notation, we will skip the lower bound α, and will work only with the constraint u≤β.

The variational inequality (<ref>) is a projection in L^2(Γ). In this case, it is equivalent to a pointwise projection:

    νu̅(x)=min{νβ(x),∂_nφ̅(x)}Γ.


In order to analyze the semismooth Newton method to solve the optimality system, we define G:L^2(Γ)→ L^2(Γ) by

    G(u) = ν u-min{νβ,-S^*(S u-y_Ω)}.

Thus, solving the optimality system is equivalent to solving

    G(u)=0.

Given u∈ L^2(Γ), we define the sets of active and free points related to u as

    Γ_A(u)={x∈Γ: -S^*(S u-y_Ω)>νβ}, Γ_F(u) = Γ∖Γ_A(u).

Abusing notation and when this does not lead to confusion, we will often drop the u and only write Γ_A and Γ_F.  (Although it is customary to use the word “inactive”, we have preferred to use “free” since we have already used the letter I for the identity matrix I and the interior nodes I.) χ_B will denote the characteristic function of a set B.
The relation G:L^2(Γ)→ L^2(Γ)
is slantly differentiable in the sense stated in <cit.> and  ∂^CLG semismooth in the sense stated in <cit.>, where ∂^CLG is Clarke's generalized differential <cit.>.

A slanting functional M(u)∈ℒ(L^2(Γ),L^2(Γ)) is given by

    M(u)v = ν v+χ_Γ_F(u)S^*Sv

for all v∈ L^2(Γ).

Finally, if ν>0, M(u) has an inverse uniformly bounded in ℒ(L^2(Γ),L^2(Γ)) for all u∈ L^2(Γ).

Using  (<ref>) in Lemma <ref>, we have  that there exists some q>2 such that S^*S∈ℒ(L^2(Γ),L^q(Γ)). The slant differentibility then follows directly from <cit.> and the semismoothness from <cit.>;  see also <cit.>.

The expression for the slanting functional follows from the slant derivative of the function min(0,·) and the chain rule <cit.>.

Finally, given z∈ L^2(Γ)

    M(u)v = z{[         ν v = z         Γ_A(u),; ν v =  z -S^*Sv         Γ_F(u). ].

This equations can be read as

    -Δ y  =0Ω,       y =vΓ, 
    
    -Δφ  = yΩ,      φ =0Γ, 
    ν v = z  Γ_A(u),      ν v = z + ∂_nφΓ_F(u).

Let us define w=χ_Γ_F(u)v. Taking into account the definition of S and S^*, we have that
ν v = ν w +  zχ_Γ_A(u) and

    -Δ y_w  =0Ω,       y_w =wΓ,   
    
    -Δφ_w  = y_wΩ,       φ_w =0Γ, 
    
    w = 0  Γ_A(u),       ν w = z + ∂_nφ_w -  S^*S z χ_Γ_A(u)Γ_F(u).

On one hand, using the continuity of S^* (cf. (<ref>)) and of S, we have

    S^*S z χ_Γ_A(u)_L^2(Γ)≤ C  z_L^2(Γ_A(u))≤ C z_L^2(Γ).

On the other hand, using the definition of solution in the transposition sense, we have that

    (y_w,y_w)_Ω = -(w,∂_nφ_w)_Γ = - ν (w,w)_Γ +  (w,z)_Γ -  (w,S^*Sz χ_Γ_A(u))_Γ.

So

    νw^2_L^2(Γ)   =     (w,z)_Γ - (w,S^*S z χ_Γ_A(u))_Γ- y_w^2_L^2(Ω)
       ≤   w_L^2(Γ)z_L^2(Γ) + w_L^2(Γ)S^*Sz χ_Γ_A(u)_L^2(Γ)
       ≤    Cw_L^2(Γ)z_L^2(Γ).

And we get w_L^2(Γ)≤ C/νz_L^2(Γ). Taking into account the definition of w and the condition for v on the active set, we get that v_L^2(Γ)≤ C/νz_L^2(Γ), where C is independent of u, and hence M(u) has a uniformly bounded inverse for each ν>0.

Notice that in the infinite dimensional case, if ν=0 then M(u)v = χ_Γ_F(u)S^*Sv. In practical cases, it is known (cf. (<ref>) or <cit.>) that there exists t>0 such that S^*Sv∈ H^t(Γ) which is compactly embedded in L^2(Γ), and hence M(u) does not have a bounded inverse.


Given a current iterate u∈ L^2(Γ), we may compute next iterate u^+ using a semismooth Newton method as follows:

    M(u)(u^+-u) = -G(u).

Writing this in detail leads to

    ν u^+-ν u +χ_Γ_F(u)S^*S(u^+-u) = -ν u+min{νβ,-S^*(Su-y_Ω)},

which means that, if ν>0,

    ν u^+-ν u = -ν u+νβ⇒ u^+ = βΓ_A(u),

and

    u^+ = -1/νS^*(S u^+-y_Ω)Γ_F(u).

This equations can be read as


    -Δ y^+ = 0Ω,        y^+ = u^+Γ,
    
    -Δφ^+ = y^+-y_ΩΩ,       φ^+ = 0Γ,
    
    u^+ = βΓ_A(u),       
    u^+ = 1/ν∂_nφ^+Γ_F(u).


With all these considerations, we can write a semismooth Newton method to solve the optimality system in the infinite dimensional case.




The semismooth Newton method described in Algorithm <ref> converges q-superlinealy to u̅ provided u^0 is close enough to u̅ in the sense of L^2(Γ).

Once we have proved Lemma <ref>, this result is a direct consequence of <cit.>; see also <cit.> or <cit.>.
 

The discrete version of the variational inequality (<ref>) –see (<ref>) below– does not have a pointwise version analog to (<ref>). A more convenient approach to the continuous problem, from the point of view of the discretized problem, is obtained using Lagrange multipliers associated to the bound control constraints.  The con of this approach is that we do not obtain a semismooth Newton method, in the sense that the involved functions are known not to be semismooth.

 It must noticed that, when the variational discretization of the control is used, <cit.>, the discrete optimal control is obtained as the pointwise projection of the discrete optimal adjoint state as in (<ref>), so convergence of the semismooth Newton method for such a discretization would follow from  Theorem 4.2  and Lemma <ref> with S replaced by S_h.




To simplify the notation, we will reduce the exposition to the case of having only an upper bound u≤β on Γ.


Condition (<ref>) can be replaced by the following pair of equations: there exists also λ̅∈ L^2(Γ) such that equation (<ref>) can be written as


    νu̅(x)   =   ∂_nφ̅(x)-λ̅(x)Γ.
    u̅(x)≤β(x), λ̅(x)≥ 0, λ̅(x)(u̅(x)-β(x))   =   0x∈Γ.


Condition (<ref>) can be written as

    λ(x)-max{0,λ(x)+c(u(x)-β(x))} =0x∈Γ,

which is true for any c>0.
Since u∈ L^2(Γ) appears inside the max operation, it is known that the operator described in (<ref>) is not semismooth (following <cit.>) or has not a slant derivative (as in <cit.>), so it is not clear that  Newton's method applied to (<ref>), (<ref>), (<ref>) converges.
Define

    F_c(u,λ) = ([ ν u+S^*Su-S^*y_Ω+λ;   λ-max{0,λ+cu-cβ} ])

and for each pair control-multiplier (u,λ)∈ L^2(Γ)× L^2(Γ), set

    Γ_A(u,λ) = {x∈Γ: 1/cλ+u>β}, Γ_F(u,λ) = Γ∖Γ_A(u,λ).

We have that

    ([        ν I+S^*S               I; - cχ_Γ_A(u,λ) I    χ_Γ_F(u,λ) I ])∈∂^CL F_c(u,λ).

So, given a current iterate (u,λ), a Newton-like iterate to obtain (u^+,λ^+) reads as


    ν u^+ = -S^*(Su^+-y_Ω) - λ^+,
    
    u^+ = βΓ_A(u,λ),
    λ^+ =0 Γ_F(u,λ).


It is remarkable that the parameter c only appears in this equations hidden in the definition of the active set.



If c=ν>0, the sequence u^k generated by the Newton-like method described in Algorithm <ref> converges q-superlinearly to u̅ in L^2(Γ) provided u^0 is close enough to u̅ in the sense of L^2(Γ) and λ^0=S^*y_Ω-(S^*Su^0 + ν u^0).

In this case we cannot apply directly the results in <cit.>; see also <cit.> or <cit.> since u∈ L^2(Γ) appears inside the max operation, and it is known that the operator described in (<ref>) is not semismooth (following <cit.>) or has not a slant derivative (as in <cit.>).

We will follow instead the method of proof of <cit.> and we will show that the sequence u^k generated by Algorithm <ref> is exactly the same as the one generated by Algorithm <ref>, which we will call ũ^k.

Writing in detail (<ref>) we obtain


    -Δ y^+ = 0 Ω,        y^+ = u^+Γ,
    
    -Δφ^+ = y^+-y_ΩΩ,       φ^+ = 0Γ,
          ν u^+ = ∂_nφ^+ -λ^+Γ.


Using now (<ref>) and (<ref>) we obtain


    -Δ y^+ = 0Ω,        y^+ = u^+Γ,
    
    -Δφ^+ = y^+-y_ΩΩ,       φ^+ = 0Γ,
    
    u^+ = βΓ_A(u,λ),       
    u^+ = 1/ν∂_nφ^+Γ_F(u,λ),
    λ^+ = ∂_nφ^+-νβΓ_A(u,λ),       λ^+=0Γ_F(u,λ).


Notice that  (<ref>), (<ref>) and (<ref>) is exactly as (<ref>), (<ref>) and (<ref>), provided that Γ_A(u,λ)= Γ_A(u).

To finish, let us prove by induction that u^k=ũ^k and Γ_A(u^k,λ^k)= Γ_A(ũ^k) for all k∈ℕ∪{0}. For k=0 it is clear from the definition of the active sets and the choice of λ^0 made in the assumption of this theorem. Suppose now u^k=ũ^k and Γ_A(u^k,λ^k)= Γ_A(ũ^k). It is clear that u^k+1=ũ^k+1. On the other hand, from (<ref>) and (<ref>) and the choice c=ν we have that c u^k+1+λ^k+1= - S^*S(u^k+1-y_Ω), and therefore Γ_A(u^k+1,λ^k+1)= Γ_A(ũ^k+1).

The result, hence, follows from Theorem <ref>.

We want to remark here that writting the Newton step as in (<ref>), (<ref>), (<ref>) and (<ref>) we obtain exactly the same kind of algorithm as the one described in <cit.>



 §.§ Finite dimensional approximation

Let us focus now on the finite dimensional approximation.
For every 0<h<h_0, there exist a unique solution u̅_h∈ U_h of (P_h^C) and unique, y̅_h∈ Y_h and φ̅_h∈ Y_h0 such that


    a(y̅_h, z_h)   =   0 z_h∈ Y_h0,
    y̅_h   ≡   u̅_hΓ,
    
    a( z_h,φ̅_h)     =   (y̅_h-y_Ω,z_h)_Ωz_h∈ Y_h0,
    
    (-∂_n^hφ̅_h+νu̅_h,u_h-u̅_h)    ≥   0u_h∈ U^h_α,β.


This way of writting the optimality system is useful to obtain error estimates (cf. <cit.>). Nevertheless, we cannot deduce a pointwise projection formula for the optimal control from the variational inequality.
Think of this naive example. Let Γ=[-1,1] and take a mesh with nodes {-1,0,1}. For every u_h,v_h∈ U_h, we have that
(u_h,v_h)_Γ = u^T Bv, where B is the mass matrix

    B=1/6([ 2 1 0; 1 4 1; 0 1 2 ]).

Consider α = -∞, β = 0, and take , e.g., w =(-2,1,1)^T. Then we have that

    (-w_h+u̅_h,u_h-u̅_h)≥ 0 ∀ u_h∈ U^h_α,βu̅_h = -1.5e_1 + 0 e_2 + 0 e_3u = (-1.5,0,0)^T,

but

    min(w,0)= (-2,0,0)^T,

and we get different results with the L^2(Γ) projection and the pointwise projection.

To circumvent this difficulty, we rewrite conditions (<ref>)–(<ref>) in order to use standard semismooth Newton method. Taking into account (<ref>), we can write (P^C_h) as

    [ min1/2u^T Au -f^Tu;             α≤u≤β, ]

where A and f are defined in (<ref>) and (<ref>) and α,β∈ℝ^N_B× 1 are the vectors whose j-component are respectively α(x_j) and β(x_j) and the inequalities are understood componentwise. In order to simplify the notation, we will restrict ourselves again to the case without lower bound.
The optimality system can thus be written in the following form: if u is the solution of problem (<ref>), then there exists a unique λ∈ℝ^N_B× 1 such that

    [             Au+λ -f                   =                  0,; λ≥ 0, u≤β, λ^T(u-β)                   =                  0. ]

Notice that although we may name λ̅_h(x) = ∑_j∈Bλ_j e_j(x) for x∈Γ to obtain a Lagrange multiplier λ̅_h∈ U_h, the complementarity condition with respect to the dot product in ℝ^N_B, does not imply that
(λ̅_h,u̅_h-β)_Γ =0,
which shows again the convenience of using (<ref>) instead of (<ref>)–(<ref>). To continue, we rewrite again the second condition in (<ref>) to obtain

    [           Au+λ -f                 =                0,; λ-max(0,λ+c(u-β))                 =                0. ]

Notice that (<ref>) looks like a discrete version of the optimality system formed by (<ref>), (<ref>) and (<ref>).
Acting in an analogous way as we did for the continuous problem, we define

    F_h,c(u,λ) = ([           Au+λ -f; λ-max(0,λ+c(u-β)) ]).

and for every pair (u,λ), we define the
sets of active  and free
  indexes as

    A(u,λ) = {j:λ_j+c(u_j-β(x_j))> 0},F(u,λ) = {j:λ_j+c(u_j-β(x_j))≤ 0}.

Abusing notation and when this does not lead to confusion, we will often just write A and F. Notice that A∪F=B for every possible pair control-multiplier.
The function F_h,c is slantly differentiable. A slanting function for F_h,c is

    M_h,c(u,λ) = (
    [             A         I_B,B; -c I_B,AI_A,B    I_B,FI_F,B ]).

Finally, for all ν≥ 0, the inverse of M_h,c(u,λ) is uniformly bounded w.r.t (u,λ).


In finite dimension, the function max{0,·} is semismooth due to Rademacher's theorem (see e.g. <cit.> or <cit.>).

A forward computation shows that M_h,c is a slanting function.


Let us prove the uniform boundedness of the inverse. Given z∈ℝ^N_B and η∈ℝ^N_B, we have that

    M_h,c(u,λ) [[ v; μ ]] =
    [[ z; η ]]{[  Av+μ     =    z,; -cv_A     =  η_A,;   μ_F     =  η_F. ].

We write the last equalities as

    Av+
    [
    [ μ_A; 0_F ]]
        =    z-
    [
    [ 0_A; η_F ]] =:ζ, 
    
    -cv_A    =    η_A.

This is the optimality system of the  equality constrained optimization problem

    v = min_w∈ℝ^N_B1/2w^TAw - ζ^Tw,
          subject to w_A = -η_A/c.

Writing w = I_ℬ,ℱw_F + I_ℬ,𝒜w_A and taking into account the equality constraint, we have that v_F is the solution of the following unconstrained optimization problem


    v_F = min_w_F1/2w_F^T
    (I_F,BAI_B,F)
    w_F -
    (I_F,B( ζ  + AI_B,Aη_A/c ))^T
    w_F

and therefore v_F is the solution of the following linear system

    I_F,BAI_B,Fv_F = I_F,B( ζ +AI_B,Aη_A /c)

Since A is symmetric and positive definite, so is I_F,BAI_B,F, and its smallest eigenvalue is bounded from below by 0<λ_1(M)+νλ_1(B_B,B); see  (<ref>). Therefore, the previous system is solvable and there exists a constant C>0, that may depend on h, ν and c, but is independent of u,  such that v_F≤ C (z+η). From this it is straight to deduce that  v+ μ≤ C (z+η) and hence M_h,c has a uniformly bounded inverse.


With these considerations, given a current iterate (u,λ) with active and free index sets A=A(u,λ) and F=F(u,λ), we can compute the next step of Newton's method
(u^+,λ^+) solving

    [ u^+ = min_u∈ℝ^N_B× 11/2u^T Au -f^Tu;                       xxxxu_A= β_A.;                     λ^+ =  f -Au^+. ]

At each iteration, this is equivalent to solving the following unconstrained optimization problem in the lower-dimensional space
ℝ^N_F× 1:

    [ u^+_F = min_u_F∈ℝ^N_F× 11/2u_F^T
    (I_F,BAI_B,F)
    u_F -
    (I_F,B( f -AI_B,Aβ_A ))^T
    u_F;                                                                         u^+_A= β_A;                                                                     λ^+ =  f -Au^+ ]

Again the preconditioned gradient method works fine to solve this problem. A good preconditioner in practice is P = M_F,F+νB_F,F.

Alternatively, taking into account the definition of A
 and f (see also Algorithms <ref> and <ref> respectively), we may write one step of the semismooth Newton algorithm as

    [ [[   M+νB -K_:,I  I_:,A; -K_I,:  O_I,I  O_I,A;  I_A,:  O_A,I  O_A,A ]]
    [[   y^+; φ^+_I; λ^+_A ]]                                                                                                =                                                                     [[   My_Ω;    0_I; I_A,Bβ ]];                                                                                              u^+                                                                                                =                                                                                            y^+_B;                                                                                            λ^+_F                                                                                                =                                                                                              0_F ]

As for the unconstrained problem, direct methods for small size problems or preconditioned conjugate gradient techniques described in <cit.> can be applied to solve this system at each Newton step.



Notice that the only information from iteration k used to compute iteration k+1 is the set of active indexes. Therefore, when A(u_k+1,λ_k+1)=A(u_k,λ_k) we have reached an stationary point. This is usually the criterion used to stop the semismooth Newton method. Another consequence of this is that to initiate the algorithm, in principle we do not not need an initial guess u^0 and λ^0, but only an initial guess for the active set. We include nevertheless an initial guess for the control variable in Algorithm <ref> because we use u_k as the initial guess for the pcg to obtain u_k+1.





The following convergence result is a direct consequence of Lemma <ref> and <cit.>; see also <cit.> or <cit.>.


The sequence u_k generated by Algorithm <ref> converges q-superlinearly to u, the solution of (<ref>).

We resume the 2D problem described in Example <ref> and the 3D problem described in Example <ref>, with the upper constraint  β≡0.16. We test Algorithm <ref>.

Following the tip of Theorem <ref>, we have taken the parameter c=ν and λ_0=f-Au_0. Nevertheless, we have not been able to observe any problem for different values of c.  The seed is set to u_0=0.

To solve the optimality system of the unconstrained optimization problem at each Newton iterate, we have used the preconditioned conjugate gradient method for (<ref>) with initial guess I_F,Bu_k. The use of the reduced problem is even more advisable in this case than it was in Examples <ref> and <ref>, because the size of the system (<ref>) is N+N_I+N_A, which in any case is greater or equal than the size of the system (<ref>), which is N+N_I.

In Tables <ref> and <ref> and we report on the number of Newton iterations for each mesh size as well as the total number of conjugate gradient iterations. For reference, we also report on the number of active nodes and the optimal solution of the discrete problem being approximated.









Choosing the initial guess for Algorithm <ref>. If possible, it is a good idea to select an initial point for the Newton method close to the solution. If we are dealing with the problem for some mesh size h, a good candidate for the initial iteration is I_hu̅_h_-1, the solution in a coarser mesh with h_-1≥ h. This idea can be iterated with a mesh family with parameters h_-M,…,h_0=h. The computation time should decrease for fine meshes provided that the interpolation can be carried out in an effective way. For instance, using nested meshes.
Solving the last 2D problem in Table <ref> takes a cputime of 620 seconds. Using a nested iteration, we solve the problem with 5120 boundary nodes in 436  seconds, with just 2 Newton iterates and 7 conjugate gradient iterates at the finest level. The times include mesh generations and matrix assembly. To generate the meshes and interpolate the solution, we have used Matlab PDEToolbox command .

In the 3D case (see Table <ref>), solving the problem for h=2^-6 takes 681 seconds. Using a nested mesh strategy this time is reduced to 404 seconds with 2 Newton iterates and 9 pcg iterations at the finest level. In our example we have been able to make an efficient interpolation of the solution at the previous level just using  Matlab's 
.

Absence of Tikhonov parameter In the same 2D pentagonal domain, we take again β≡ 0.16, but now α≡ -1.2, ν = 0, y_Ω≡ 1 if x_1>0.25, and y_Ω≡ -1 if x_1<0.25. We are able to solve the finite dimensional problem, but up to our best knowledge, there is no proof available about the convergence of the discrete optimal solutions to the continuous optimal solution.

To solve the problem we follow a nested mesh strategy. The number of pcg iterations per Newton iterate grow as h tends to zero, as is to be expected from Corollary <ref>; see Table <ref>.

The constraints have been chosen in such way that we seemingly find both a bang-bang part and a singular arc. We have sketched a plot of the discrete optimal control for h=0.2× 2^-8 in Figure <ref>. The boundary has been stretched on a 1D line and its corners have been marked with circles on the graph of u̅_h.









§ STATE CONSTRAINTS


In the rest of the work, we will suppose that d=2, Ω is convex and Γ is polygonal. According to <cit.>, if problem (P^S) admits a feasible Slater point, then it has a unique solution u̅∈ H^1/2(Γ) and there exist y̅∈ H^1(Ω)∩ C(ω̅), φ̅∈ W^1,t_0(Ω) for all t<2 and two nonnegative measures μ̅^+, μ̅^-∈ℳ(ω̅) such that


    -Δy̅ = 0Ω, y̅ = u̅Γ,
    
    -Δφ̅= y̅-y_Ω+μ̅^+-μ̅^-Ω, φ̅= 0Γ,
    u̅(x) = 1/ν∂_nφ̅(x)Γ,
    ⟨μ̅^+-μ̅^-,y-y̅⟩≤ 0 ∀ y∈ K_a,b,


and μ̅^+⊂{y̅ = b}, μ̅^-⊂{y̅ = a}. In this case the adjoint state equation (<ref>) must be understood in the transposition sense, <cit.>, and ⟨·,·⟩ denotes the duality product between ℳ(ω̅) and C(ω̅).

As is pointed out in <cit.>, or <cit.>, the PDAS or semismooth Newton methods described in the previous section are not applicable to this problem, since the multiplier is a measure. In <cit.> a Moreau-Yosida regularization is proposed. This strategy is further investigated in <cit.>. Let us briefly describe this method. Again, to simplify the notation, we will do the exposition for the unilateral constraint y≤ b in ω̅.

Given a shift function μ^*∈ L^q(ω) for some q>2 and a parameter γ>0, we will solve the unconstrained problem

    (Q^γ) min_u∈ L^2(Γ) J(u) +1/2γ∫_ω̅max{0,μ^*+γ(Su-b)}^2 dx.

This problem has a unique solution u^γ∈ H^1/2(Γ) for every γ >0. Since the functional is of class C^1, first order optimality conditions can be derived directly from the work <cit.>: there exist unique y^γ∈ H^1(Ω) and φ^γ∈ H^1_0(Ω) such that


    -Δ y^γ   = 0Ω, y^γ = u^γΓ,
    
    -Δφ^γ   = y^γ-y_Ω+χ_ω̅max{0,μ^*+γ(y^γ-b)}Ω, φ^γ = 0Γ,
    ν u^γ   = ∂_nφ^γΓ.



The semismooth Newton method to solve (<ref>), (<ref>), (<ref>) converges locally q-superlinearly.


We may write the system (<ref>), (<ref>), (<ref>) as the equation G(u)=0, where G:L^2(Γ)→ L^2(Γ) is given by

    G(u) = ν u+S^*S u-S^*y_Ω +S^*( χ_ω̅max{0,μ^*+γ(Su-b)}).

Using the regularity results in <cit.>, together with the election of μ^* in L^q(ω)  we have that μ^*+γ(Su-b)∈ L^q(Ω) for some q>2, and hence G is semismooth in the sense stated in <cit.>; see Theorems 2.13 and 2.10(c) in the aforementioned reference.

Define now ω_A(u)={x∈ω̅: μ^*+γ(Su-b) > 0}. A slant differential of G(u) is given by the expression


    M(u) v =  ν v + S^*Sv +γ S^*χ_ω_A(u) Sv.

Let us see that the inverse of M(u) is uniformly bounded in ℒ(L^2(Γ),L^2(Γ)) for all u∈ L^2(Γ). For any z∈ L^2(Γ), the equation M(u)v=z can be read as

    -Δ y = 0Ω,       y=vΓ,
    
    -Δφ = y+γχ_ω_A(u) yΩ,      φ=0Γ,
          ν v  = ∂_nφ + zΓ.

Using the definition of solution in the transposition sense and the last equation of this system, we have that

    (y,y+γχ_ω_A(u) y)_Ω = -(v,∂_nφ) = -ν(v,v)_Γ+ν(v,z)_Γ

and hence

    νv_L^2(Γ)^2 = ν(v,z)_Γ-∫_Ω y^2 -γ∫_ω_A(u)y^2≤νv_L^2(Γ)z_L^2(Γ).

So we have that v_L^2(Γ)≤z_L^2(Γ) and therefore M(u) is uniformly bounded in ℒ(L^2(Γ),L^2(Γ)) for all u∈ L^2(Γ).

 When γ→+∞, u^γ→u̅; see <cit.> or <cit.>.



Let us turn now to the discrete problem (P^S_h). Consider the space ℳ_h⊂ℳ(ω̅) which is spanned by the Dirac measures corresponding to the nodes {x_j}_j∈J, where J={j x_j∈ω̅}. Following <cit.>, if the continuous problem (P^S) has a regular feasible Slater point, then (P^S_h) has a unique solution u̅_h. Moreover  there exist unique y̅_h∈ K_a,b^h, φ̅_h∈ Y_h,0 and two nonnegative measures μ̅_h^+, μ̅_h^-∈ℳ_h such that


    a( y̅_h,z_h)   =0 ∀ z_h∈ Y_h0, y̅_h = u̅_hΓ,
    
    a( z_h,φ̅_h)   =(y̅_h-y_Ω,z_h)+
    ⟨μ̅_h^+-μ̅_h^-,z_h⟩ ∀ z_h∈ Y_h0,
    ⟨μ̅_h^+-μ̅_h^-,y_h-y̅_h⟩   ≤ 0 ∀ y_h∈ K_a,b^h,
    νu̅_h   =∂_n^hφ̅_hΓ,


and μ̅_h^+⊂{x_j: y̅(x_j)=b(x_j)}, μ̅_h^-⊂{x_j: y̅(x_j)=a(x_j)}.

Since we are dealing with a finite dimensional problem and the max function is known to be semismooth in finite dimension, we could think about applying directly a semismooth Newton method as described in <cit.>. Nevertheless, the other assumption fails here: the slant derivative may not have an inverse.
Let us show this. Consider

    F(u,μ) = ([           Au-f+S^Tμ; μ -max{ 0,μ+γ(Su-b) ]).

It is clear that the optimality system (<ref>)–(<ref>) is equivalent to the equation F(u,μ)=0. We define the sets of active and free nodes related to a pair control-multiplier as

    A(u_h,μ) = {j∈J:μ_j+γ(S_h(u_h)(x_j)-b_j) >0}, F = J∖A.

The slant derivative of F(u,μ) is given by

    M(u,μ) = ([             A           S^T; -γI_J,AI_A,:S    I_J,FI_F,J ]).

The inverse of M(u,μ) may not exist for some (u,μ).


For any (z,δ)∈ℝ^N_B×ℝ^N_J, we have that


    M(u,μ)([ v; η ]) =
    ([ z; δ ])
    {[                    Av + S^Tη                            =                           z,; -γI_J,AI_A,:Sv +
    I_J,FI_F,Jη                            =                           δ. ].

From the second equation we have that

    γI_A,:Sv = - I_A,Jδ.

If N_A > N_B we have a linear system with more equations than variables, and hence it will not be consistent for at least one value of δ and M will not have an inverse.





Instead, we use a Moreau-Yosida penalization of (P^S_h). Again we will write only the case of unilateral upper constraint. For some shift function μ^*∈ Y_h such that μ^*(x_j)=0 if j∉J, and a parameter γ>0, a direct discretization of (Q^γ) could be

    (Q̃^γ_h)min_u_h∈ U_hJ_h(u_h)+1/2γ∫_ω̅max{0,μ^*+γ(S_hu_h-b)}^2dx.

Problem (Q̃^γ_h) has a unique solution u^γ_h∈ U_h that converges to u^γ as h→ 0. Since the functional is not of class C^2, the problem does not fit exactly into the framework of <cit.>, so we will give a sketch of the proof of this fact.
Problem (Q̃^γ_h) has a unique solution u^γ_h∈ U_h. Moreover, there exist unique y_h^γ∈ Y_h and φ_h^γ∈ Y_h0 such that


    a( y^γ_h,z_h)   =0 ∀ z_h∈ Y_h0,  y^γ_h =  u^γ_hΓ,
    
    a(z_h,φ^γ_h)   =( y^γ_h-y_Ω+
    max{0,μ^*+γ(y_h^γ-b)},z_h)_Ω ∀ z_h∈ Y_h0,
    ν u^γ_h   =∂_n^hφ^γ_hΓ.


Finally, u_h^γ→ u^γ in L^2(Γ) and there exists some λ>0 such that

    u^γ-u^γ_h_L^2(Γ)≤ C h^λ.


Existence and uniqueness of solution follows immediately from the coerciveness of the discrete penalized functional. Since this discrete functional is also of class C^1, first order optimality conditions are standard.

Consider two controls u_1,u_2∈ L^2(Γ) and its related states y_i and adjoint states φ_i w.r.t. problem (Q^γ). Using the definition of solution in the transposition sense for the states and both the state and adjoint state equations (<ref>) and (<ref>) we have (shorting z_i = μ^*+γ(y_i-b), i=1,2)

    (-∂_nφ_1 + ∂_nφ_2,u_1-u_2)_Γ   =   
    (y_1-y_2,-Δφ_1+Δφ_2)_Ω
        =     (y_1-y_2, y_1+χ_ωmax{0,z_1}-y_2-χ_ωmax{0,z_2})_Ω
       =   y_1-y_2^2_L^2(Ω)+1/γ(z_1 - z_2,max{0,z_1}-max{0,z_2})_ω
       ≥    0.

The rest is quite standard:

    νu^γ- u^γ_h_L^2(Γ)^2   ≤   
    (-∂_nφ^γ+ν u^γ+∂_nφ_u_h^γ- ν u_h^γ, u^γ-u_h^γ)_Γ
        =     (-∂_nφ^γ+∂_n^hφ_h^γ+ν u^γ-ν u_h^γ, u^γ-  u_h^γ)_Γ
           +(-∂_n^hφ_h^γ+∂_nφ_u_h^γ,u^γ- u_h^γ)_Γ.

For fixed γ, u_h^γ is uniformly bounded in L^2(Γ) because u_h^γ_L^2(Γ)^2≤ 1/νy_Ω_L^2(Ω)^2 + 1/γ/νmax{0,μ^*-γ b}_L^2(Ω)^2, and therefore using <cit.> we have that -∂_n^hφ_h^γ+∂_nφ_u_h^γ≤ C h^λ for some λ>0. For the first addend we obtain directly 0 using

first order optimality conditions (<ref>) and (<ref>).









A practical way of computing the integral in the penalty term is to use the lumped mass matrix L∈ℝ^N× N. Since we are going to use it only for integrals in ω̅, we define it as

    L_i,j = 0i≠ j, L_j,j = 0j∉J, L_j,j= ∑_k=1^N M_k,j.j∈J.

Therefore, we will be solving

    (Q^γ_h)min_u_h∈ U_hJ_h(u_h)+1/2γ∑_j∈JL_j,jmax{0,μ^*(x_j)+γ(S_hu_h(x_j)-b(x_j))}^2.

 First order optimality conditions for this problem read like


    a( y^γ_h,z_h)   =0 ∀ z_h∈ Y_h0,  y^γ_h =  u^γ_hΓ,
    
    a( z_h,φ^γ_h)   =( y^γ_h-y_Ω,z_h)_Ω
          +∑_j∈JL_j,jmax{0,μ^*(x_j)+γ(y_h^γ(x_j)-b(x_j))}z_h(x_j)
     ∀ z_h∈ Y_h0,
    ν u^γ_h   =∂_n^hφ^γ_hΓ.


This nonlinear system can be solved using a semismooth Newton method.
To this end, we will use the nonlinear operator G_h:ℝ^N_B→ℝ^N_B defined by

    G_h(u) =
    Au-f+
    S^TLmax{0,μ^*+γ(Su-b)}.

Thus, solving (<ref>)–(<ref>) is equivalent to solving G_h(u)=0.

For every u_h∈ U_h
we define sets of active and free nodes as

    A_ω(u_h,γ,h) = {j∈J: μ^*_j+γ(S_hu_h(x_j)-b_j) > 0}, F_ω(u_h,γ,h) =I∖A_ω(u_h,γ,h).

Abusing notation we will drop some or all of the arguments or will use the vector notation when this does not lead to confusion, e.g. A_ω, A_ω(u).
Notice that if μ^*≡ 0, then A_ω(u_h,γ,h) is independent of γ.
For the sake of notation, it is also convenient to define for every set A_ω⊆J the diagonal matrix H(A_ω)∈ℝ^N× N such that

    H_i,j=δ_i,j{[      0  j∉A_ω;  L_j,j j∈A_ω. ].

Abusing notation, we will often write H(u)=H(A_ω(u)) or even we will write H when this does not lead to confusion. The proof of the following result is as the corresponding one in infinite dimension.
G_h(u) is slantly differentiable, a slant differential is given by

    M(u)v = Av+γS^TH(u) Sv

and it has a uniformly bounded inverse w.r.t. u.



Notice that using this H notation, we can write

    G_h(u) =
    Au-f+
    S^TH(u)(μ^*+γ(Su-b)),

and therefore, for a given u, and denoting H=H(u), one Newton step reads like

    Au^+ +γS^THSu^+ = f+
    S^TH(γb-μ^*).

Let us comment that for the computation of w = S^THy for some y∈ℝ^N, first we solve

    K_I,Iϕ_s = H_I,:y,

and, next, we have

    w = H_B,:y-K_B,Iϕ_s = -K_B,Iϕ_s

because H is diagonal and its nonzero components  correspond to nodes that lie in ω̅, and are hence interior to Ω.

Again a preconditioned conjugate gradient method can be used to solve this system, provided an efficient way of computing d=(A + γS^THS) v and the second member of the system. Notice that in each of the algorithms <ref> and <ref> the computation of ϕ can be done with just one system solve.










Alternatively, the solution to (<ref>) can be obtained solving


    [ [[ M+νB+γH  -K_:,I;  -K_I,:   O_I,I ]]
    [[   y^+; φ_I^+ ]]                                                         =                    [[ My_Ω+H(γb-μ^*));               0 ]];                                                       u^+                                                         =                                                    y^+_B. ]









Finally, the semismooth Newton algorithm to solve our problem reads as:



As we noticed in Remark <ref>,  the only information from iteration k used to compute iteration k+1 is the set of active indexes. Therefore, when A_ω(u_k+1)=A_ω(u_k) we have reached an stationary point. This is usually the criterion used to stop the semismooth Newton method.

Nevertheless, as we will see below
 in the context of algorithms <ref> and <ref>, the main use for the solution of (Q_h^γ) will be to provide an initial guess for the next step in those procedures, so it does not seem necessary to solve exactly (Q_h^γ) in general. Following <cit.>, we may implement the following stopping criterion for Algorithm <ref>. After step <ref>, we compute an approximation of the multiplier μ_k+1∈ℝ^N
as

    μ_k+1,j = 0j∉A_ω,k, μ_k+1,j = μ^*_j+γ(y_k+1,j-b_j)j∈A_ω,k

and we may stop the algorithm if

    μ_k+1-max(μ^*+γ(y_k+1-b),0_L^2(Ω)< ε_λ.

In practice, this quantity gives a good measure of the change in the active set between iterates as well as the unfeasibility combined with the penalization parameter γ.


We state the convergence result for the Newton method. It follows directly from Theorem <ref>
The semismooth Newton method to solve (Q^γ_h) described in Algorithm <ref> converges locally q-superlinearly.

In all the examples below we have taken the shift μ^*=0.

In general it is not a good idea to solve directly (Q_h^γ) for some  γ big enough. The resulting intermediate problems are usually very ill conditioned.

We resume the 2D problem taken from <cit.> and described in Example <ref>. As in <cit.>, we take ω the ball centered at (-0.1,-0.1) with radius 0.2 and b≡ 0.15. We use Algorithm <ref> taking u_0=0. The results are summarized in Table <ref>. We measure the unfeasibility of the state as the maximum constraint violation mcv(y_h) = max(y_h-b,0)_L^∞(ω̅).


The first thing that can be observed is that γ cannot be too big w.r.t. the problem size. Although the algorithm should converge in finite time, a value too big for γ will make the active sets of the intermediate steps fluctuate, and in practice it may not stop. This is what happens in this case with meshes coarser than the first one exposed in the table. On the other hand, it can also be noticed that the amount of Newton iterations looks mesh independent for h small enough. This was to be expected because the infinite dimensional version of the method is convergent (cf. Theorem <ref>).

Nevertheless, the computational effort can be better measured by the total number of conjugate gradient iterations made when we solve (<ref>) in step <ref> of Algorithm <ref>. Solving the solution in the finest mesh takes 8246 seconds.




In <cit.> a continuation strategy is proposed. This reduces considerably the computational effort. We will use subscripts for the iteration number of the semismooth Newton method described in Algorithm <ref> and superscripts for the numbering of the iterates of the continuation strategy provided in Algorithm <ref>. For some sequence {τ^n}_n≥ 1 such that τ^n>1 for all n≥ 1, we have



There are still three important details to be explained about Algorithm <ref>: first and most important, the choice of the initial guess u^0; next, the appropriate values for τ^n; and finally, a suitable stopping criterion.

To obtain a good initial guess, we will solve the problem in a coarser mesh with a smaller value of γ. Since we are going to deal with meshes of different sizes, we will use the notation A_ω(u_h,γ,h) when needed. We also recall that I_h is the pointwise interpolation operator from C(Γ̅) onto U_h. Given (h,γ^0), we fix M>0 and pick a finite sequence (h_n,γ_n)_n=0^M such that h_n is decreasing, h_M = h, γ_n is increasing and γ_M≤γ^0.


Although Algorithm <ref> is meaningful for non-nested meshes, the computational effort needed to make the interpolation in step <ref> in non-nested meshes can be considerable, specially for 3D problems.

In <cit.> some criteria are given to choose τ^n+1>1. In practice, if τ^n+1  is very small the algorithm would not advance; if it is very big, we would lose the advantage given by the continuation strategy.


In <cit.> the authors stop the continuation strategy if residuals related to the state, the adjoint state and the multiplier are smaller than a certain tolerance of order O(h^2). In our case, since the problem is linear quadratic, the residuals related to the state and the adjoint state are zero (at least up to roundoff error). The residual related to the multiplier can be computed as

    r_d = ∑_j∈A_ω^n+1L_j,j (y^n+1_j-b_j).

As an alternative, a tolerance for the maximum constraint violation e_∞ can be used stop if mcv(y_h^γ_n)≤ e_∞.

With all these considerations, we propose the following algorithm. Fix a mesh sequence such that {h_j}_j=0^M is decreasing, γ_0>0, C>0, u_h_0∈ U_h_0, n_max∈ℕ and a sequence {τ^n}_n=1^n_max+1 such that τ^n>1 for all n.



Now we apply Algorithm <ref> for γ_0=1, τ^n=10 and a family of 9 nested meshes of sizes h_j=0.2*2^-j, j=0:8, the coarsest for h_0 with 52 nodes and the finest for h_8 with  2689537 nodes.

A summary of the results can be found in Table <ref>. The total computation time was 1570 seconds, which is a significant improvement compared with the 8246 seconds used by Algorithm <ref> (see Example <ref>).



Determination of (the support of) the Lagrange multipliers of (P^S). Comparison of the adjoint state equations for (P_h^S) and (Q_h^γ) leads to the approximation formula

    μ̅_h^+ ≈∑_j∈A_ω(u_h^γ)L_j,j (μ^*_j +γ(y^γ_h(x_j)-b_j))δ_x_j.


Since ω̅⊂Ω, it is very common that the Lagrange multipliers of the original problem are finite sums of Dirac measures centered at points on the boundary of ω, say μ̅^+ = ∑_k=1^n μ_kδ_X_k for some n>0 and X_k∈∂ω.

Increasing γ will have as a result a better approximation of the Lagrange multiplier and its support. Notice, nevertheless, this will done with great effort and the approximation of both the control, the state and the functional optimal value will not improve in a significant way.

In Example 5.7 the number of active nodes for the last mesh is 748. All of them are on ∂ω but they are not isolated. We can see in Table <ref> how the number of active nodes decreases as γ increases, with little variation of the functional or the state.





§ CONTROL AND STATE CONSTRAINTS


Again according to <cit.>, if problem (P^CS) admits a feasible Slater point, then it has a unique solution u̅∈ H^1/2(Γ) and there exist y̅∈ H^1(Ω)∩ C(ω̅), φ̅∈ W^1,t_0(Ω) for all t<2 and two nonnegative measures μ̅^+, μ̅^-∈ℳ(ω̅) such that


    -Δy̅ = 0Ω, y̅ = u̅Γ,
    
    -Δφ̅= y̅-y_Ω+μ̅^+-μ̅^-Ω, φ̅= 0Γ,
    u̅(x) = min{β(x),max{α(x),1/ν∂_nφ̅(x)}}Γ,
    ⟨μ̅^+-μ̅^-,y-y̅⟩≤ 0 ∀ y∈ K_a,b


and μ̅^+⊂{y̅ = b}, μ̅^-⊂{y̅ = a}. As we said in the previous section, a semismooth Newton strategy for this problem is meaningless, so instead we are going to deal with a Moreau-Yosida approximation.
As we did in the previous sections, we will consider only unilateral constraints u≤β on Γ and y≤ b in ω̅ to simplify the notation.
For a shift function μ^*∈ L^q(ω) for some q>2 and a parameter γ>0, we consider the problem

    (Q^C,γ) min_u∈ U_-∞,βJ(u)+1/2γ∫_ω̅max{0,μ^*+γ(Su-b)}^2dx.

This problem has a unique solution u^γ∈ H^1/2(Γ). Moreover, there exist y^γ∈ H^1(Ω), φ^γ∈ H^s(Ω), s>3/2, such that


    -Δ y^γ =    0Ω,  y^γ =  u^γΓ,
    
    -Δφ^γ =    y^γ-y_Ω+max{0,μ^*+γ(y^γ-b)}Ω, φ^γ = 0Γ,
    
     (-∂_nφ^γ+ν u^γ,u-u^γ)   ≥ 0u∈ U_-∞,β.


Define

    G^γ(u) = ν u-min{νβ,-S^*Su+S^*y_Ω-S^*χ_ω̅max{0,μ^*+γ(Su-b)}}.

It is clear that u^γ is the unique solution of (Q^c,γ) if and only if G^γ(u^γ)=0.
For some fixed shift function μ^*∈ L^q(ω), q>2, consider the active sets

    Γ_A(u,γ) = {x∈Γ: -S^*Su+S^*y_Ω-S^*χ_ωmax{0,μ^*+γ(Su-b)>νβ}

and

    ω_A(u,γ) = {x∈ω̅: μ^*+γ(Su-b)>0}.

A slant differential of G^γ(u) is given by

    M^γ(u)v = ν v+χ_Γ_A(u,γ)S^*(1+ γχ_ω_A(u,γ))Sv.



G^γ is slantly differentiable, M^γ(u) is a slant differential of G^γ and for every fixed ν>0, M^γ(u) has an inverse in ℒ(L^2(Γ),L^2(Γ))  uniformly bounded for all u∈ L^2(Γ).
The semismooth Newton method M^γ(u^+-u) = -G^γ(u) converges q-superlinearly.

The proof follows the same lines as those of Lemma <ref> and Theorems <ref> and <ref>

As we did for the pure control-constrained case, to deal with a problem better suited to the finite dimensional case, we write the optimality condition (<ref>) with the help of a Lagrange multiplier.
There exists λ^γ∈ L^2(Γ) such that, for any c>0,


    ν u^γ =    ∂_nφ^γ-λ^γΓ,
    λ^γ=   max{0,λ^γ+c(u^γ-β)}Γ.


Define now

    F_c^γ(u,λ) = (
    [ ν u + S^*Su-S^*y_Ω + S^*χ_ω̅max{0,μ^*+γ(Su-b)}+λ;                            λ - max{0,λ + cu -cβ} ])

We have that

    ([ ν I+S^*(1+ γχ_ω_A(u,γ))S                        I;             cχ_Γ_A(u,λ)I              χ_Γ_F(u,λ)I ])∈∂^CLF_c^γ(u,λ)

and hence a Newton-like method to solve F_c^γ(u,λ)=0  is given by

    ν u^+    =    -S^*(1+ γχ_ω_A(u,γ))Su^+-y_Ω) -λ^+,
    
    u^+    =   βΓ_A(u,λ),  
    λ^+    =    0Γ_F(u,λ).

Although F_c^γ is known not to be slantly differentiable, it can be proved as in Theorem <ref> that the sequence generated by the above described Newton-like method to solve F_c^γ(u,λ)=0 is the same than the one generated by the semismooth Newton method to solve G^γ(u)=0 provided we take the same initial guess u_0, c=ν and λ_0=S^*(y_Ω-(1+ γχ_ω_A(u_0,γ))Su_0)-ν u_0.


Let us turn now to the finite dimensional problem.
We can write the approximation of problem (Q^C,γ_h) as a constrained optimization problem in ℝ^N_B,1:

    (Q^C,γ_h){[ min1/2u^TAu-f^Tu+
    1/2γmax{0,μ^*+γ(Su-b)}^TLmax{0,μ^*+γ(Su-b)};                                                          u≤β. ].

Existence and uniqueness of solution of this problem, as well as error estimates for the difference u^γ-u^γ_h_L^2(Γ) can be proved as we did for (Q̃^γ) in Lemma <ref>. Since we have control constraints, instead of the point-wise interpolation, the Casas-Raymond interpolate <cit.> should be used in the last step.


First order optimality conditions read as

    [ Au+S^TLmax{0,μ^*+γ(Su-b)}+λ                           =                          f,;           λ-max(0,λ+c(u-β))                           =                          0. ]

Using the definitions (<ref>) and (<ref>) for the active sets of indexes A(u,λ) and A_ω(u,γ,h) and the matrix H related to A_ω, we have that one step of Newton's method can be written as

    [                                                                        u^+_F                                                                            =                                 min_u_F1/2u_F^T
    (I_F,B(A + γS^THS)I_B,F)
    u_F;                                                                             ;                                                                              xxxxxx -
    ( I_F,B
    ( f
    - AI_B,Aβ_A
    - S^TH ( μ^* +γ (SI_B,Aβ_A -b) )
     )
    )^T
    u_F;                                                                        u^+_A                                                                            =                                                                          β_A;                                                                          λ^+                                                                            =                                         f - Au^+ - S^TH ( μ^* +γ (Su^+ -
    b)) ]

or, alternatively, as

    [ [[ M+νB+γH  -K_:,I   I_:,A;  -K_I,:   O_I,I   O_I,A;   I_B,:   O_B,I   O_B,A;         ]]
    [[   y^+; φ_I^+;   λ^+ ]]                                                                                                                  =                                                               [[ My_Ω+H(γb-μ^*);              0;         I_A,Bβ ]];                                                                                                                u^+                                                                                                                  =                                                                                                              y^+_B;                                                                                                              λ^+_F                                                                                                                  =                                                                                                                 0. ]



An adaptation of algorithms <ref> and <ref> to solve (Q^C,γ_h) is straightforward, and so is an adaptation of Algorithm <ref> to use a continuation strategy together with a nested mesh strategy.
We repeat Example <ref> adding the control constraint u≤ 0.16. We obtain the results summarized in Table <ref>. The total computation time was 2894 seconds. Although this may seem a lot of time, we have to remember that we are solving an optimal control problem with pointwise constraints in both the control and the state using a mesh with almost 2.7E+06 nodes and 5.4E+06 elements to test the algorithm.





§ FUNDING

The author was partially
supported by the Spanish Ministerio de Ciencia e Innovación under Project MTM2014-57531-P.




10



Casas-Raymond2006
Casas E, Raymond JP. Error estimates for the numerical approximation of
  Dirichlet boundary control for semilinear elliptic equations. SIAM J
  Control Optim. 2006;45(5):1586–1611.
  <http://dx.doi.org/10.1137/050626600>.

Deckelnick-Gunther-Hinze2009
Deckelnick K, Günther A, Hinze M. Finite element approximation of
  Dirichlet boundary control for elliptic PDEs on two- and
  three-dimensional curved domains. SIAM J Control Optim.
  2009;48(4):2798–2819.
  <http://dx.doi.org/10.1137/080735369>.

Gunther2010
Günther A. Structure exploiting Galerkin schemes for optimal control of
  PDEs with constraints on the involved variables [dissertation].
  Universität Hamburg; 2010.
  <http://ediss.sub.uni-hamburg.de/volltexte/2010/4926/pdf/PhDThesis_Guenther.pdf>.

MayRannacherVexler2013
May S, Rannacher R, Vexler B. Error analysis for a finite element approximation
  of elliptic Dirichlet boundary control problems. SIAM J Control Optim.
  2013;51(3):2585–2611.
  <http://dx.doi.org/10.1137/080735734>.

Mateos-Neitzel2015
Mateos M, Neitzel I. Dirichlet control of elliptic state constrained problems.
  Comput Optim Appl. 2016;63(3):825–853.
  <http://dx.doi.org/10.1007/s10589-015-9784-y>.

Casas-Sokolowski2010
Casas E, Sokolowski J. Approximation of boundary control problems on curved
  domains. SIAM J Control Optim. 2010;48(6):3746–3780.
  <http://dx.doi.org/10.1137/090761550>.

Casas-Gunther-Mateos2011
Casas E, Günther A, Mateos M. A paradox in the approximation of Dirichlet
  control problems in curved domains. SIAM J Control Optim.
  2011;49(5):1998–2007.
  <http://dx.doi.org/10.1137/100794882>.

AMPR2015
Apel T, Mateos M, Pfefferer J, et al. On the regularity of the solutions of
  Dirichlet optimal control problems in polygonal domains. SIAM J Control
  Optim. 2015;53(6):3620–3641.
  <http://dx.doi.org/10.1137/140994186>.

Gong-Hinze-Zhou2016
Gong W, Hinze M, Zhou Z. Finite element method and a priori error estimates for
  Dirichlet boundary control problems governed by parabolic PDEs. J Sci
  Comput. 2016;66(3):941–967.
  <http://dx.doi.org/10.1007/s10915-015-0051-2>.

Chen-Nashed-Qi2000
Chen X, Nashed Z, Qi L. Smoothing methods and semismooth methods for
  nondifferentiable operator equations. SIAM J Numer Anal.
  2000;38(4):1200–1216 (electronic).
  <http://dx.doi.org/10.1137/S0036142999356719>.

Bergounioux-Ito-Kunisch1999
Bergounioux M, Ito K, Kunisch K. Primal-dual strategy for constrained optimal
  control problems. SIAM J Control Optim. 1999;37(4):1176–1194
  (electronic). <http://dx.doi.org/10.1137/S0363012997328609>.

Bergounioux-Kunisch2002
Bergounioux Mt, Kunisch K. Primal-dual strategy for state-constrained optimal
  control problems. Comput Optim Appl. 2002;22(2):193–224.
  <http://dx.doi.org/10.1023/A:1015489608037>.

Kunisch-Rosch2002
Kunisch K, Rösch A. Primal-dual active set strategy for a general class of
  constrained optimal control problems. SIAM J Optim.
  2002;13(2):321–334.
  <http://dx.doi.org/10.1137/S1052623499358008>.

Ulbrich2003
Ulbrich M. Semismooth Newton methods for operator equations in function
  spaces. SIAM J Optim. 2002;13(3):805–842 (2003).
  <http://dx.doi.org/10.1137/S1052623400371569>.

Hintermuller-Ito-Kunisch2003
Hintermüller M, Ito K, Kunisch K. The primal-dual active set strategy as a
  semismooth Newton method. SIAM J Optim. 2003;13(3):865–888
  (2003). <http://dx.doi.org/10.1137/S1052623401383558>.

Hintermuller-Kunisch2006
Hintermüller M, Kunisch K. Feasible and noninterior path-following in
  constrained minimization with low multiplier regularity. SIAM J Control
  Optim. 2006;45(4):1198–1221.
  <http://dx.doi.org/10.1137/050637480>.

Casas-Mateos-Raymond-2009
Casas E, Mateos M, Raymond JP. Penalization of Dirichlet optimal control
  problems. ESAIM Control Optim Calc Var. 2009;15(4):782–809.
  <http://dx.doi.org/10.1051/cocv:2008049>.

Hinze-Vierling2012
Hinze M, Vierling M. The semi-smooth Newton method for variationally
  discretized control constrained elliptic optimal control problems;
  implementation, convergence and globalization. Optim Methods Softw.
  2012;27(6):933–950.
  <http://dx.doi.org/10.1080/10556788.2012.676046>.

Hinze2005
Hinze M. A variational discretization concept in control constrained
  optimization: the linear-quadratic case. Comput Optim Appl.
  2005;30(1):45–61.
  <http://dx.doi.org/10.1007/s10589-005-4559-5>.

Hinze_Pinnau_Ulbrich_Ulbrich2009
Hinze M, Pinnau R, Ulbrich M, et al. Optimization with PDE constraints.
  (Mathematical Modelling: Theory and Applications; Vol. 23). New York:
  Springer; 2009.

Gong-Yan2011
Gong W, Yan N. Mixed finite element method for Dirichlet boundary control
  problem governed by elliptic PDEs. SIAM J Control Optim.
  2011;49(3):984–1014.
  <http://dx.doi.org/10.1137/100795632>.

Casas-Mateos-Troeltz2005
Casas E, Mateos M, Tröltzsch F. Error estimates for the numerical
  approximation of boundary semilinear elliptic control problems. Comput Optim
  Appl. 2005;31(2):193–219.
  <http://dx.doi.org/10.1007/s10589-005-2180-2>.

Grisvard1992
Grisvard P. Singularities in boundary value problems. (Recherches en
  Mathématiques Appliquées [Research in Applied Mathematics]; Vol. 22).
  Masson, Paris; Springer-Verlag, Berlin; 1992.

Jerison-Kenig1995
Jerison D, Kenig C. The inhomogeneous Dirichlet problem in Lipschitz
  domains. J Funct Anal. 1995;130(1):161–219.
  <http://dx.doi.org/10.1006/jfan.1995.1067>.

Jerison-Kenig1981N
Jerison DS, Kenig CE. The Neumann problem on Lipschitz domains. Bull Amer
  Math Soc (NS). 1981;4(2):203–207.
  <http://dx.doi.org/10.1090/S0273-0979-1981-14884-9>.

Dahlberg-Kenig1987
Dahlberg BEJ, Kenig CE. Hardy spaces and the Neumann problem in L^p for
  Laplace's equation in Lipschitz domains. Ann of Math (2).
  1987;125(3):437–465.
  <http://dx.doi.org/10.2307/1971407>.

Berggren2004
Berggren M. Approximations of very weak solutions to boundary-value problems.
  SIAM J Numer Anal. 2004;42(2):860–877 (electronic).
  <http://dx.doi.org/10.1137/S0036142903382048>.

AMPR2017
Apel T, Mateos M, Pfefferer J, et al. Error estimates for Dirichlet control
  problems in polygonal domains. 2017;In preparation.

Schoberl-Zulehner2007
Schöberl J, Zulehner W. Symmetric indefinite preconditioners for saddle
  point problems with applications to PDE-constrained optimization problems.
  SIAM J Matrix Anal Appl. 2007;29(3):752–773 (electronic).
  <http://dx.doi.org/10.1137/060660977>.

Herzog-Sachs2010
Herzog R, Sachs E. Preconditioned conjugate gradient method for optimal control
  problems with control and state constraints. SIAM J Matrix Anal Appl.
  2010;31(5):2291–2317.
  <http://dx.doi.org/10.1137/090779127>.

Schiela-Ulbrich2014
Schiela A, Ulbrich S. Operator preconditioning for a class of inequality
  constrained optimal control problems. SIAM J Optim.
  2014;24(1):435–466.
  <http://dx.doi.org/10.1137/120877532>.

ASW1996
Apel T, Sändig AM, Whiteman JR. Graded mesh refinement and error estimates
  for finite element solutions of elliptic boundary value problems in
  non-smooth domains. Math Methods Appl Sci. 1996;19(1):63–85.

Kamenski-Huang-Xu2014
Kamenski L, Huang W, Xu H. Conditioning of finite element equations with
  arbitrary anisotropic meshes. Math Comp. 2014;83(289):2187–2211.
  <http://dx.doi.org/10.1090/S0025-5718-2014-02822-6>.

Dauge1989
Dauge M. Problèmes mixtes pour le laplacien dans des domaines polyédraux
  courbes. C R Acad Sci Paris Sér I Math. 1989;309(8):553–558.

nevcas2012direct
Nečas J. Direct methods in the theory of elliptic equations. Corrected
  2nd ed. Springer Berlin Heidelberg; 2012. Monographs and Studies in
  Mathematics.

Plaza2007
Plaza A. The eight-tetrahedra longest-edge partition and Kuhn triangulations.
  Comput Math Appl. 2007;54(3):427–433.
  <http://dx.doi.org/10.1016/j.camwa.2007.01.023>.

Clarke1983
Clarke F. Optimization and nonsmooth analysis. Toronto: John Wiley & Sons;
  1983.

Ito-Kunisch2003
Ito K, Kunisch K. Semi-smooth Newton methods for state-constrained optimal
  control problems. Systems Control Lett. 2003;50(3):221–228.
  <http://dx.doi.org/10.1016/S0167-6911(03)00156-7>.



