Optimization Methods for Dirichlet Control Problems

TL;DR

This paper explores various optimization algorithms for finite element approximations of linear-quadratic Dirichlet optimal control problems governed by elliptic PDEs, analyzing convergence and computational efficiency.

Contribution

It introduces a preconditioned conjugate method and semismooth Newton methods for constrained and unconstrained problems, with convergence analysis and numerical validation.

Findings

Convergence of optimization procedures without Tikhonov regularization.

Effective solution methods for control and state constrained problems.

Numerical examples demonstrating theoretical results.

Abstract

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.

Equations

$$-\Delta y=0\text{mbox}in\Omega ,y=u\text{mbox}on\Gamma .$$

$$J(u)=\frac{1}{2}\Vert Su-y_{\Omega }{\Vert }_{L^2(\Omega )}^2+\frac{\nu }{2}\Vert u{\Vert }_{L^2(\Gamma )}^2$$

$$U_{\alpha ,\beta }=\{u\in L^2(\Gamma ):\alpha (x)\le u(x)\le \beta (x)\text{mbox}fora.e.x\in \Gamma \},$$

$$K_{a,b}=\{y\in L^2(\Omega )\cap C(\bar{\omega }):a(x)\le y(x)\le \beta (x)\text{mbox}fora.e.x\in \bar{\omega }\}.$$

$$(P^{U})\min_{u\in L^{2}(\Gamma)}J(u),\ \ \ \ (P^{C})\min_{u\in U_{\alpha,\beta}}J(u),\ \ \ \ (P^{S})\min_{Su\in K_{a,b}}J(u),\ \ \ \ (P^{CS})\min_{\begin{array}[]{c}\scriptstyle{u\in U_{\alpha,\beta}}\\ \scriptstyle{Su\in K_{a,b}}\end{array}}J(u),$$

$$Ay:=-\sum _{i=1}^d{\partial }_i(a_{i,j}{\partial }_{j}y)+a_{0}y=F\text{mbox}in\Omega ,y=u+G\text{mbox}on\Gamma$$

$$\Phi =\{\varphi :\varphi \in H_0^1(\Omega )\text{mbox}and\Delta \varphi \in L^2(\Omega )\}.$$

$$\Vert {\partial }_n\varphi {\Vert }_{L^q(\Gamma )}\le C\Vert \Delta \varphi {\Vert }_{L^2(\Omega )}.$$

$$\Vert {\partial }_n\varphi {\Vert }_{H^t(\Gamma )}\le C\Vert \Delta \varphi {\Vert }_{L^2(\Omega )}.$$

$$\Vert {\partial }_{n}w{\Vert }_{L^q(\Gamma )}\le C\Vert \Delta \varphi {\Vert }_{L^2(\Omega )};$$

$$\Vert g{\Vert }_{W^{1,q}(\Gamma )}\le C\Vert \Delta \varphi {\Vert }_{L^2(\Omega )}.$$

$$\Vert \mathrm{M}(\nabla v){\Vert }_{L^q(\Gamma )}\le C\Vert g{\Vert }_{W^{1,q}(\Gamma )}.$$

$$\Vert {\partial }_{n}v{\Vert }_{L^q(\Gamma )}\le C\Vert \Delta \varphi {\Vert }_{L^2(\Omega )}.$$

$$(y,-\Delta \varphi {)}_{\Omega }=-(u,{\partial }_n\varphi {)}_{\Gamma }\text{mbox}forall\varphi \in \Phi .$$

$$-\Delta \varphi =z\text{mbox}in\Omega ,\varphi =0\text{mbox}on\Gamma .$$

$$J(u)$$

$$J^{\prime }(u)v=(S^{\ast }Su+\nu u,v{)}_{\Gamma }-(S^{\ast }{y}_{\Omega },v{)}_{\Gamma }.$$

$$-\Delta {\phi }_u=y_u-y_{\Omega }\text{mbox}in\Omega ,{\phi }_u=0\text{mbox}on\Gamma .$$

$$Y_h=\{y_h\in C(\bar{\Omega }):y_h\in P^1(T)\forall T\in {\mathcal{T}}_h\}.$$

$$U_h=\{u_h\in C(\Gamma ):u_h\in P^1(T\cap \Gamma )\forall T\in {\mathcal{T}}_h\}.$$

$$({\Pi }_{h}y,y_h{)}_{\Omega }=(y,y_h{)}_{\Omega }\forall y_h\in Y_h\text{mbox}and({\Pi }_{h}u,u_h{)}_{\Gamma }=(u,u_h{)}_{\Gamma }\forall u_h\in U_h.$$

$$a(y_h,z_h)=0\forall z_h\in Y_{h0},(y_h,v_h{)}_{\Gamma }=(u,v_h{)}_{\Gamma }\forall v_h\in U_h,$$

$$J_h(u)=\frac{1}{2}(S_{h}u-y_{\Omega },S_{h}u-y_{\Omega }{)}_{\Omega }+\frac{\nu }{2}(u,u{)}_{\Gamma }.$$

$$U_{\alpha ,\beta }^h=\{u_h\in U_h:\alpha (x_j)\le u_h(x_j)\le \beta (x_j)\forall j\in \mathbb{B}\},$$

$$K_{a,b}^h=\{y_h\in Y_h:a(x_j)\le y_h(x_j)\le b(x_j)\forall j\in \mathbb{J}\}.$$

$$(P_{h}^{U})\min_{u_{h}\in U_{h}}J_{h}(u_{h}),\ \ \ \ (P_{h}^{C})\min_{u\in U^{h}_{\alpha,\beta}}J_{h}(u_{h}),\ \ \ \ (P_{h}^{S})\min_{S_{h}u_{h}\in K^{h}_{a,b}}J_{h}(u_{h}),\ \ \ \ (P_{h}^{CS})\min_{\begin{array}[]{c}\scriptstyle{u_{h}\in U^{h}_{\alpha,\beta}}\\ \scriptstyle{S_{h}u_{h}\in K^{h}_{a,b}}\end{array}}J_{h}(u_{h}).$$

$$(w_h,z_h{)}_{\Gamma }=(y,z_h{)}_{\Omega }-a(z_h,{\varphi }_h)\text{mbox}forall{z}_h\in Y_h,$$

$$a(z_h,{\varphi }_h)=(y,z_h{)}_{\Omega }\text{mbox}forall{z}_h\in Y_{h0}.$$

$$(S_{h}u,y{)}_{\Omega }=(u,S_h^{\ast }y{)}_{\Gamma }.$$

$$J_h(u)=\frac{1}{2}(S_h^{\ast }{S}_{h}u+\nu u,u{)}_{\Gamma }-(S_h^{\ast }{y}_{\Omega },u{)}_{\Gamma }+c_{\Omega }.$$

Full text

{amscode}

49J20, 49M05, 49M15, 65N30

1 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 [1]. 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 [2] or [3, Ch 2.1] for a variational approach to control-constrained problems posed on smooth 2D or 3D domains, [4] for a superconvergence result on the state approximation for unconstrained 2D problems, [5] for control-and-state constrained problems posed on convex polygonal domains. In [6] 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 [7] an apparent paradox between [2] and [6] is explained. The regularity of the solution in possibly nonconvex polygonal plane domains is studied in [8]; see also the introduction of that paper for further references about related problems. In the recent publication [9], 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 [10, 11, 12, 13, 14, 15, 16] 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 [16] 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 [17] 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 [18], 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 4 can be applied to Dirichlet control. See Remark 4 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 2 we state the problem and prove that it is well posed in Lipschitz domains; see Lemma 2.1. 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 [19] 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 [20, Ch. 3.2.7.3] and in [3, Example 2.1.11] this difficulty is solved using a mixed formulation of the state equation which is discretized with the lowest order Raviart-Thomas element. In [21] 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 [1], 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 [2] 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 2 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 3 an efficient method to solve the unconstrained problem. We propose in Subsection 3.2 a preconditioned conjugate gradient (pcg in the rest of the work) method for a reduced problem in the spirit of [22, Section 5]. We are able to prove convergence of the conjugate gradient even for the case where the Tikhonov regularization parameter $\nu$ vanishes.

In sections 4, 5 and 6 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 4 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 [2, 3, 19] or a mixed formulation, as studied in [21, 20, 3]. When we use continuous piecewise linear elements,the variational inequality cannot be written in a point-wise or elementwise form; see (36d) and Remark 5. We include the analysis of Newton methods for both the original problem and the discretized one.

In Section 5 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 [16] is developed. Finally, in Section 6 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 [1, 5, 8] and 2D or 3D smooth domains [2] 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 3.6.

2 Statement of the problem and first properties

Let $\Omega\subset\mathbb{R}^{d}$, $d=2$ or $d=3$, be a bounded domain with Lipschitz boundary $\Gamma$ and in this domain consider a target state $y_{\Omega}\in L^{2}(\Omega)$. Consider also the continuous linear operator $S:L^{2}(\Gamma)\to L^{2}(\Omega)$ such that $y=Su$ if and only if $y$ is the solution in the transposition sense (see Definition 2.2 below) of

[TABLE]

Let $\omega\subset\Omega$ be a domain such that $\bar{\omega}\subset\Omega$ and define two pairs of functions $\alpha,\ \beta\in C(\Gamma)$ and $a,\ b\in C(\bar{\omega})$ such that $\alpha(x)<\beta(x)$ for all $x\in\Gamma$ and $a(x)<b(x)$ for all $x\in\bar{\omega}$. For some fixed regularization parameter $\nu\geq 0$, define

[TABLE]

and consider the sets

[TABLE]

and

[TABLE]

In this work we will study optimization procedures for the following four Dirichlet control problems:

[TABLE]

namely the unconstrained, control-constrained, state-constrained and control-and-state constrained problems.

Remark 1*.*

Almost all the literature related to these problems is written using the state equation (1). There would be no problem in taking into account an equation of the form

[TABLE]

with regular enough $F$, $G$, $a_{0}\geq 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

[TABLE]

This is a Banach space with the graph norm $\|\phi\|_{\Phi}=\|\phi\|_{H^{1}(\Omega)}+\|\Delta\phi\|_{L^{2}(\Omega)}$. Further, the functions in this space satisfy $\partial_{n}\phi\in L^{2}(\Gamma)$. This is known to be true for smooth domains; convex domains, see [17, Lemma A.1]; plane polygonal domains, see [8, Corollary 2.3]; or polyhedral domains, see [23, Theorem 2.6.7] and the usual trace theorem. We have not been able to find a proof of this fact for general Lipschitz domains. In [24, Theorem B.2] the regularity $\phi\in H^{3/2}(\Omega)$ is proved. Nevertheless, as the authors notice in page 165 of this reference, the trace theorem is not valid neither in $H^{3/2}(\Omega)$ nor in $H^{1/2}(\Omega)$, so we cannot deduce immediately that $\partial_{n}\phi\in L^{2}(\Gamma)$. The results in [25, 24] 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 $\partial_{n}\phi\in L^{2}(\Gamma)$ in Lipschitz domains. Regarding the analysis of semismooth Newton methods –see Lemma 4.1 below– we also prove $L^{q}(\Gamma)$ regularity for some $q>2$.

Lemma 2.1**.**

Let $\Omega\subset\mathbb{R}^{d}$, $d=2$ or $d=3$, be a bounded domain with Lipschitz boundary $\Gamma$ and consider $\phi\in\Phi$. Then, there exists $q_{0}>2$ depending on the domain such that, for $2\leq q<q_{0}$, we have $\partial_{n}\phi\in L^{q}(\Gamma)$ and

[TABLE]

If, further, $\Omega$ is smooth or convex or polygonal or polyhedral, then there exists $t>0$ such that $\partial_{n}\phi\in H^{t}(\Gamma)$ and

[TABLE]

Proof.

Denote $z=-\Delta\phi$, extend $z$ by zero to $\mathbb{R}^{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\in H^{2}(\Omega)$ and $\nabla w\in H^{1}(\Omega)^{d}$, so it is clear that $\mathrm{Tr}(\nabla w)\in H^{1/2}(\Gamma)^{d}\hookrightarrow L^{q}(\Gamma)^{d}$ for all $q<\infty$ if $d=2$ and all $q\leq 4$ if $d=3$, since the dimension of $\Gamma$ is $d-1$. This implies that:

(a) $\partial_{n}w=\nabla w\cdot n\in L^{q}(\Gamma)$ because the boundary $\Gamma$ is Lipschitz and therefore it has a unit normal vector defined almost everywhere; and there exists $C>0$ such that

[TABLE]

(b) $g=\mathrm{Tr}(w)\in W^{1,q}(\Gamma)$ due precisely to the definition of $W^{1,q}(\Gamma)$, and there exists $C>0$ such that

[TABLE]

Define now $v\in H^{1}(\Omega)$ the unique solution of $-\Delta v=0$ in $\Omega$, $v=g$ on $\Gamma$. Using [24, Theorem 5.6], we have that there exists $q_{0}>2$ such that if $2\leq q<q_{0}$, then the nontangential maximal function of the gradient of $v$ satisfies $\mathrm{M}(\nabla v)\in L^{q}(\Gamma)$ and there exists $C>0$ such that

[TABLE]

As is pointed out in [26, p. 438], this implies that $\nabla v$ has nontangential limit a.e. on $\Gamma$ and we can define the normal derivative of $v$ at a point $s\in\Gamma$ as the nontangential limit as $x\to s$ of $\nabla v(x)\cdot n(s)$. For a precise definition of nontangential limit and nontangential maximal function see, e.g., the introduction of the work [26]. This, together with inequalities (6) and (5) imply that $\partial_{n}v\in L^{q}(\Omega)$ and

[TABLE]

So we have that $\phi=v-w$ and we can define in a natural way $\partial_{n}\phi=\partial_{n}v-\partial_{n}w\in L^{q}(\Gamma)$. The estimate (2) follows from (4) and (7).

For smooth, convex, polygonal or polyhedral domains, the second result follows from the regularity $\phi\in H^{3/2+t}(\Omega)$ and the usual trace theorem for $\nabla\phi$. See [17, Lemma A.1] for convex domains, [8, Corollary 2.3] for plane polygonal domains and [23, Theorem 2.6.7] for polyhedral domains. ∎

Definition 2.2**.**

We will say that $y\in L^{2}(\Omega)$ is the solution in the transposition sense of (1) if

[TABLE]

Here and in the rest of the work $(\cdot,\cdot)_{X}$ stands for the standard inner product in $L^{2}(X)$.

The adjoint operator of $S$ is $S^{*}:L^{2}(\Omega)\to L^{2}(\Gamma)$ defined by $S^{*}z=-\partial_{n}\phi$, where $\phi$ is the unique weak solution of

[TABLE]

We can write now

[TABLE]

where $c_{\Omega}=0.5\|y_{\Omega}\|_{L^{2}(\Omega)}^{2}$ is a constant.

From this expression, we can easily compute the derivative of $J$ at a point $u\in L^{2}(\Gamma)$ in the direction $v\in L^{2}(\Gamma)$:

[TABLE]

For later use, we will define now for every $u\in L^{2}(\Gamma)$, $y_{u}=Su\in H^{1/2}(\Gamma)$ and $\varphi_{u}\in H^{1}_{0}(\Omega)$ the unique solution of

[TABLE]

Discretization.

To discretize the problems. consider $\{\mathcal{T}_{h}\}_{h}$ a regular family of triangulations of $\bar{\Omega}$. To simplify the notation, we will suppose that $\Gamma$ 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 $\bar{\omega}$ as $\mathbb{I}=\{j:x_{j}\in\Omega\}$, $\mathbb{B}=\{j:x_{j}\in\Gamma\}$ and $\mathbb{J}=\{j:x_{j}\in\bar{\omega}\}$. For the discretization of the state and the adjoint state we use the space of linear finite elements $Y_{h}\subset H^{1}(\Omega)$,

[TABLE]

As usual, we will abbreviate $Y_{h0}=Y_{h}\cap H^{1}_{0}(\Omega)$. For the control we use the space $U_{h}$ of continuous piecewise linear functions on $\Gamma$

[TABLE]

Notice that the elements of $U_{h}$ are the traces of the elements of $Y_{h}$. We will denote $I_{h}:C(\bar{\Omega})\to Y_{h}$ or $I_{h}:C(\Gamma)\to U_{h}$ the interpolation operator related to these spaces and $\Pi_{h}:L^{2}(\Omega)\to Y_{h}$ or $\Pi_{h}:L^{2}(\Gamma)\to U_{h}$ the projection onto this spaces in the $L^{2}$ sense. For all $y\in L^{2}(\Omega)$ and all $u\in L^{2}(\Gamma)$:

[TABLE]

We discretize the state equation following the work by Berggren [27]: define $S_{h}:L^{2}(\Gamma)\to L^{2}(\Omega)$ such that for $u\in L^{2}(\Gamma)$, $y_{h}=S_{h}u$ if and only if $y_{h}\in Y_{h}$ is the unique solution of

[TABLE]

where $a(\cdot,\cdot)$ is the bilinear form associated to the operator in the PDE. In the case of the Laplace operator $a(y,z)=\int_{\Omega}\nabla^{T}y\nabla zdx$. The discrete functional is thus defined as

[TABLE]

We define now

[TABLE]

and

[TABLE]

The discrete problems reads thus as

[TABLE]

The adjoint operator of $S_{h}$ is given by the discrete variational normal derivative. See [1]. $S_{h}^{*}:L^{2}(\Omega)\to L^{2}(\Gamma)$ and $w_{h}=S_{h}^{*}y$ if $w_{h}\in U_{h}$ satisfies

[TABLE]

where $\phi_{h}\in Y_{h0}$ is the unique solution of

[TABLE]

It is customary to write $S_{h}^{*}y=-\partial_{n}^{h}\phi_{h}$. For $u\in L^{2}(\Gamma)$ and $y\in L^{2}(\Omega)$, we have

[TABLE]

We can then write

[TABLE]

The computation of the derivative of $J_{h}$ at a point $u\in L^{2}(\Gamma)$ in the direction $v\in L^{2}(\Gamma)$ is then obtained with

[TABLE]

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_{\mathbb{I}}+N_{\mathbb{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 $\Gamma$. Define the usual finite element stress, mass and boundary mass matrices by

[TABLE]

We will also use $\mathcal{I}\in\mathbb{R}^{N\times N}$ for the identity matrix, $\mathcal{O}\in\mathbb{R}^{N\times N}$ for the zero matrix and usually refer to submatrices as $\mathcal{K}_{\mathbb{I},\mathbb{I}}$ or $\mathcal{K}_{\mathbb{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 $\mathcal{B}_{\mathbb{B},\mathbb{B}}$.

Given $u_{h}=\sum_{j\in\mathbb{B}}u_{j}e_{j}$, we denote $\bm{u}\in\mathbb{R}^{N_{\mathbb{B}}\times 1}$ the vector $(u_{1},\dots,u_{N_{\mathbb{B}}})^{T}$ and for $y_{h}=\sum_{j=1}^{N}y_{j}e_{j}$ we denote $\bm{y}\in\mathbb{R}^{N\times 1}$ the vector $(y_{1},\dots,y_{N})^{T}$. Using (10), we have that $y_{h}=S_{h}u_{h}$ iff

[TABLE]

Since $\mathcal{B}_{\mathbb{B},\mathbb{B}}$ is nonsingular, we can write this as

[TABLE]

If we define $\mathcal{S}\in\mathbb{R}^{N_{\mathbb{B}}\times N}$ as

[TABLE]

we have that $y_{h}=S_{h}u_{h}$ if and only if $\bm{y}=\mathcal{S}\bm{u}$.

Given $y_{h}\in Y_{h}$, let $\phi_{h}=\sum_{j\in\mathbb{I}}\phi_{j}e_{j}$ be the solution of (12) for $y=y_{h}$. Denoting $\bm{\phi}\in\mathbb{R}^{N_{\mathbb{I}}\times 1}$ the corresponding vector, it can be computed as the solution of

[TABLE]

Then we could compute $\bm{w}\in\mathbb{R}^{N_{\mathbb{B}}\times 1}$, the vector whose components are the coefficients of the (minus) discrete normal derivative $-\partial_{n}^{h}\phi_{h}=S_{h}^{*}y_{h}=w_{h}=\sum_{j\in\mathbb{B}}w_{j}e_{j}$ solving the system

[TABLE]

Formally, we can also write that $\bm{w}=\mathcal{B}_{\mathbb{B},\mathbb{B}}^{-1}\mathcal{S}^{T}\mathcal{M}\bm{y}.$ To finish this section we also define the matrix $\mathcal{A}\in\mathbb{R}^{N_{\mathbb{B}},\times N_{\mathbb{B}}}$ and the vector $\bm{f}\in\mathbb{R}^{N_{\mathbb{B}}\times 1}$ by

[TABLE]

We have that

[TABLE]

We also note that

[TABLE]

To compute the vector $\bm{f}$, we consider the projection of $y_{\Omega}$ onto $Y_{h}$ in the $L^{2}(\Omega)$ sense, $y_{\Omega,h}=\Pi_{h}y_{\Omega}$ and denote $\bm{y}_{\Omega}\in\mathbb{R}^{N\times 1}$, the vector whose $j$-th component is $y_{\Omega,h}(x_{j})$, and

[TABLE]

So for $u_{h},v_{h}\in U_{h}$, the latter represented by the vector $\bm{v}$, we can compute

[TABLE]

Notice that applying the equality (16), the explicit computation of $\bm{w}$ is not necessary, and we can write

[TABLE]

3 Unconstrained problem

Problem $(P^{U})$ has a unique solution $\bar{u}\in L^{2}(\Gamma)$ that satisfies $J^{\prime}(\bar{u})=0$. For every $0<h<h_{0}$, problem $(P^{U}_{h})$ has also a unique solution $\bar{u}_{h}$ that satisfies $J^{\prime}_{h}(\bar{u}_{h})=0$. The problems being convex, these conditions are also sufficient. Moreover it is known that $\bar{u}_{h}\to\bar{u}$ strongly in $L^{2}(\Gamma)$ and also error estimates are available in terms of the mesh size $h$ if $\nu>0$ and the domain is either 2D and polygonal or smooth. See [1, 2, 6, 4, 5, 28].

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

3.1 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 $\bar{u}_{h}\in U_{h}$, $\bar{y}_{h}\in Y_{h}$ and $\bar{\varphi}_{h}\in Y_{h0}$ satisfying the optimality system:

[TABLE]

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

[TABLE]

We may eliminate $\bm{u}$ with the boundary condition, and reorder the equations to get the linear system

[TABLE]

Notice that system (21) is solvable even for $\nu=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 [29] and Herzog and Sachs in [30]. 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_{\mathbb{I}}$. In [30, Algorithm 3], the systems are solved using a multigrid method.

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

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

Let us see another way of solving $(P^{U}_{h})$. Using (17), first order optimality conditions can also be written as

[TABLE]

Lemma 3.1**.**

The matrix $\mathcal{A}$ is symmetric and positive definite for all $\nu\geq 0$ and there exists $C>0$ independent of $h$ and $\nu$ such that its condition number is bounded by

[TABLE]

where $0<\lambda_{1}(\mathcal{B}_{\mathbb{B},\mathbb{B}})<\lambda_{N_{\mathbb{B}}}(\mathcal{B}_{\mathbb{B},\mathbb{B}})$ and $0<\lambda_{1}(\mathcal{M})<\lambda_{N}(\mathcal{M})$ are respectively the smallest and greatest eigenvalues of the matrices $\mathcal{B}_{\mathbb{B},\mathbb{B}}$ and $\mathcal{M}$.

Proof.

It is clear from (18) that $\mathcal{A}$ is symmetric. The mass matrices $\mathcal{M}$ and $\mathcal{B}_{\mathbb{B},\mathbb{B}}$ are symmetric and positive definite and therefore $\lambda_{1}(\mathcal{M})>0$ and $\lambda_{1}(\mathcal{B}_{\mathbb{B},\mathbb{B}})>0$. Since the boundary components of $S\bm{u}$ are exactly $\bm{u}$, we have that $\|\mathcal{S}\bm{u}\|_{\mathbb{R}^{N}}\geq\|\bm{u}\|_{\mathbb{R}^{N_{\mathbb{B}}}}$ and hence

[TABLE]

so $\mathcal{A}$ is positive definite.

From [4, Equation (3.8)], we know that there exists some $C>0$ such that $\|S_{h}u_{h}\|_{L^{2}(\Omega)}\leq\sqrt{C}\|u_{h}\|_{L^{2}(\Gamma)}$. Since

[TABLE]

and

[TABLE]

we obtain

[TABLE]

The bound for the condition number follows directly from the last inequality and (24). ∎

For a graded family of mesh with grading parameter $0<\mu\leq 1$ (see [32]), we usually define $h=N^{-1/d}$, $N$ being the number of nodes of the mesh and $d$ the dimension of $\Omega$. The case $\mu=1$ corresponds to that of a quasi-uniform mesh family.

Corollary 3.2**.**

If $\mathcal{T}_{h}$ is a family of graded meshes with grading parameter $0<\mu\leq 1$, then there exists $C>0$ independent of $h$ and $\nu$ such that

[TABLE]

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

[TABLE]

Proof.

from [33] and [32], we have that there exists constants $0<C_{1}<C_{2}$

[TABLE]

Estimate (25) follows then from (23). Estimate (26) follows from (25) for $\mu=1$. ∎

The explicit computation of the matrix $\mathcal{A}$ is out of the question, since it requires the solution of $2N_{\mathbb{B}}$ systems of size $N_{\mathbb{I}}$. Much better performance can be achieved using an iterative method which only requires the computation of $\bm{d}=\mathcal{A}\bm{u}$. This is shown in Algorithm 1.

Preconditioned conjugate gradient methods can be carried out at the cost of one evaluation of $\mathcal{A}\bm{u}$ per iteration. The computation of $\bm{f}$ can be done in a similar way, but of course it only must be computed once; see Algorithm 2.

To finish this section, let us say a word about a preconditioner to solve (22). In practice, matrices like $\mathcal{B}_{\mathbb{B},\mathbb{B}}$ or $\mathcal{P}=\mathcal{M}_{\mathbb{B},\mathbb{B}}+\nu\mathcal{B}_{\mathbb{B},\mathbb{B}}$ make acceptable preconditioners. This is specially important when using graded meshes; see Example 3.4 below. Notice that at each pcg iterate, we only need to solve two linear systems, each of size $N_{\mathbb{I}}$, to compute $\mathcal{A}\bm{u}$, plus another one of size $N_{\mathbb{B}}$ to compute $\mathcal{P}^{-1}\bm{r}$.

Example 3.3**.**

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 (21) or using an iterative method, the preconditioned conjugate gradient method in our case, to solve the reduced problem (22). We have used respectively Matlab builtin commands mldivide and pcg.

We use the example in [5], where $\Omega$ 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_{\Omega}\equiv 1$ and $\nu=1$.

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

For the pcg we use as initial guess the null vector, a relative tolerance of 1E-10 and the preconditioner $\mathcal{P}=\mathcal{M}_{\mathbb{B},\mathbb{B}}+\nu\mathcal{B}_{\mathbb{B},\mathbb{B}}$. At each iteration we have to solve the linear systems (14) and (15). To do this, we first obtain the Cholesky factors of $\mathcal{K}_{\mathbb{I},\mathbb{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 $\mathcal{K}_{\mathbb{I},\mathbb{I}}$ using Matlab command symamd. 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 $\infty$ the experiments in which we have run out of RAM.

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

Remark 2*.*

We have also coded Algorithm 1 in [30], 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 [30]) $\mathcal{X}=\mathcal{K}+\mathcal{B}$, $\mathcal{X}=\mathcal{K}+\mathcal{M}$ or $\mathcal{X}=\mathcal{K}+\mathcal{M}+\mathcal{B}$ for the scalar product in $Y_{h}$ and $\mathcal{Q}=\mathcal{K}_{\mathbb{I},\mathbb{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 3.3 and $\nu=1$ we obtain 126, 173 and 236 iterations respectively.

Example 3.4**.**

In this example we show the effect of the use of $\mathcal{P}=\mathcal{M}_{\mathbb{B},\mathbb{B}}+\nu\mathcal{B}_{\mathbb{B},\mathbb{B}}$ as a preconditioner to solve (22). In this case $\nu=0.01$, $\Omega=\{x=re^{i\theta}\in\mathbb{R}^{2}:\ r<1,\ 0<\theta<11\pi/12\}$, and we use a mesh family graded at the origin with parameter $\mu=13/33\approx 0.4$; see [32]. In Table 3.4 we compare the number of iterations used to reach the prescribed accuracy of $1$E$-10$ without and with preconditioner. This example shows very clearly the effectiveness of this strategy of preconditioning.

Example 3.5**.**

In this example we show that the number of pcg iterations for each fixed $h$ is independent of the value of $\nu\geq 0$ for all $\nu<\nu_{0}$. We take the same domain and meshes as in the 2D problem in Example 3.3 and set $y_{\Omega}(x)=|x|^{2}$, so that $J(\bar{u})>0$ even for $\nu=0$. The results are summarized in Table 3.5.

Example 3.6**.**

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 [2], 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 [2, Lemma 2.1]) much higher than the one we may obtain in polyhedral domains.

For this example we will take $\Omega$ the unit cube $(-0.5,0.5)^{3}$, $y_{\Omega}\equiv 1$ and $\nu=1$. First, we obtain a regularity result for the solution of our problem.

Proposition 3.7**.**

Let $\Omega\subset\mathbb{R}^{3}$ be the interior of a rectangular parallelepiped and $y_{\Omega}\in L^{p}(\Omega)$ for some $3<p<+\infty$. Consider $\bar{u}$ the solution of problem $(P^{U})$. Then, $\bar{u}\in W^{1-1/p}(\Gamma)$, $\bar{y}\in W^{1,p}(\Omega)$ and $\bar{\varphi}\in W^{2,p}(\Omega)$. Moreover, $\bar{u}\equiv 0$ on the edges and corners of $\Gamma$.

Proof.

Since $\bar{u}\in L^{2}(\Gamma)$, using transposition and interpolation it is clear that $y\in H^{1/2}(\Omega)$. Classical Sobolev embedding in 3D leads to $y\in L^{3}(\Omega)$. Since $\Omega$ is a parallelepiped and $y_{\Omega}$ is regular enough, [34, Theorem 1] states that $\bar{\varphi}\in W^{2,3}(\Omega)$, and hence $\partial_{n}\bar{\varphi}\in\Pi_{i}W^{2/3,3}(\Gamma_{i})$, where $\Gamma_{i}$, $i=1:6$, are the faces of $\Gamma$. This does not imply immediately that $\bar{u}$ belongs to $W^{2/3,3}(\Gamma)$ because $2/3\cdot 3=2$, which is the topological dimension of $\Gamma$, and some integral compatibility condition should be imposed on the edges; but for all $q<3$, it is true that $\bar{u}\in W^{1-1/q,q}(\Gamma)$, and therefore $\bar{y}\in W^{1,q}(\Omega)$. Using again Sobolev embeddings, we have that $\bar{y}\in L^{s}(\Omega)$ for all $s<+\infty$ and hence $\bar{y}-y_{\Omega}\in L^{p}(\Omega)$. Applying once again [34, Theorem 1], we have that $\bar{\varphi}\in W^{2,p}(\Omega)$.

Now we have that $\partial_{n}\bar{\varphi}\in\Pi_{i}W^{1-1/p,p}(\Gamma_{i})$ and we can prove that if we define $\partial_{n}\bar{\varphi}=0$ on the edges of $\Omega$, then we obtain a continuous function. To do this, we use a similar argument to the one used in [7, Section 4] for 2D problems. Since $p>3$, $\bar{\varphi}\in C^{1}(\bar{\Omega})$. Consider two faces $A$ and $B$ with a common edge $AB$ and let $\tau^{1}_{X}$, $\tau^{2}_{X}$ be two linearly independent vectors tangent to face $X$ ($X=A$ or $X=B$) such that $\tau^{1}_{A}=\tau^{1}_{B}$ is also tangent to the edge $AB$. Since $\bar{\varphi}=0$ on $\Gamma$, we have that for every $x\in A$, $\nabla\bar{\varphi}(x)\cdot\tau^{1}_{A}=\nabla\bar{\varphi}(x)\cdot\tau^{2}_{A}=0$ and for every $x\in B$, $\nabla\bar{\varphi}(x)\cdot\tau^{1}_{B}=\nabla\bar{\varphi}(x)\cdot\tau^{2}_{B}=0$. So for $x\in AB$, we have that $\nabla\bar{\varphi}(x)=0$ and therefore $\nabla\bar{\varphi}(x)\cdot n$ can be seen as a continuous function if we set its value to [math] on the edges, despite the jump discontinuities of the normal vector $n$.

So $\bar{u}$ is continuous and hence $\bar{u}\in W^{1-1/p,p}(\Gamma)$. The regularity of the optimal state follows from the trace theorem; see e.g. [35]. ∎

Using this regularity, an error estimate can be proved for the example problem.

Proposition 3.8**.**

Let $\Omega$ be a rectangular parallelepiped and $y_{\Omega}\in L^{p}(\Omega)$ for some $3<p<+\infty$. Consider $\bar{u}$ be the solution of problem $(P^{U})$ and $\bar{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}$

[TABLE]

The proof follows the same guidelines as those of [1] or [5, Section 6.1] 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}(\Gamma)$ 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 $\Gamma$ 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. [36]). Up to our best knowledge, the current version of the Matlab PDEToolbox (version 2.1) computes an approximation of $\mathcal{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, $\mathcal{B}_{\mathbb{B},\mathbb{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 3.6.

4 Control constrained problems

Problem $(P^{C})$ has a unique solution $\bar{u}\in L^{2}(\Gamma)$ that satisfies $J^{\prime}(\bar{u})(u-\bar{u})\geq 0$ for all $u\in U_{\alpha,\beta}$. For every $0<h<h_{0}$, problem $(P^{C}_{h})$ has also a unique solution $\bar{u}_{h}$ that satisfies $J^{\prime}_{h}(\bar{u}_{h})(u_{h}-\bar{u}_{h})\geq 0$ for all $u_{h}\in U^{h}_{\alpha,\beta}$. The problems being convex, these conditions are also sufficient. Moreover it is known that $\bar{u}_{h}\to\bar{u}$ strongly in $L^{2}(\Gamma)$ and also error estimates are available in terms of the mesh size $h$ when $\nu>0$ and the domain is convex in 2D. See [1, 6]. In [2] 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 [3, Ch. 2.1], the problem is solved using a fixed point method, convergent for $\nu$ large enough; see [18] and Remark 4 below for the convergence of the semismooth Newton method.

4.1 Continuous problem

We can formulate first order optimality conditions as: there exist unique $\bar{u}\in L^{2}(\Gamma)$, $\bar{y}\in H^{1/2}(\Omega)$ and $\bar{\varphi}\in H^{1}_{0}(\Omega)$ such that

[TABLE]

We first describe a semismooth Newton method to solve this optimality system and prove a convergence result for it (see Theorem 4.2). 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 4.3. To facilitate the notation, we will skip the lower bound $\alpha$, and will work only with the constraint $u\leq\beta.$

The variational inequality (27c) is a projection in $L^{2}(\Gamma)$. In this case, it is equivalent to a pointwise projection:

[TABLE]

In order to analyze the semismooth Newton method to solve the optimality system, we define $G:L^{2}(\Gamma)\to L^{2}(\Gamma)$ by

[TABLE]

Thus, solving the optimality system is equivalent to solving

[TABLE]

Given $u\in L^{2}(\Gamma)$, we define the sets of active and free points related to $u$ as

[TABLE]

Abusing notation and when this does not lead to confusion, we will often drop the $u$ and only write $\Gamma_{A}$ and $\Gamma_{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 $\mathcal{I}$ and the interior nodes $\mathbb{I}$.) $\chi_{B}$ will denote the characteristic function of a set $B$.

Lemma 4.1**.**

The relation $G:L^{2}(\Gamma)\to L^{2}(\Gamma)$ is slantly differentiable in the sense stated in [15] and $\partial^{CL}G$ semismooth in the sense stated in [20, Ch. 2], where $\partial^{CL}G$ is Clarke’s generalized differential [37].

A slanting functional $M(u)\in\mathcal{L}(L^{2}(\Gamma),L^{2}(\Gamma))$ is given by

[TABLE]

for all $v\in L^{2}(\Gamma)$.

Finally, if $\nu>0$, $M(u)$ has an inverse uniformly bounded in $\mathcal{L}(L^{2}(\Gamma),L^{2}(\Gamma))$ for all $u\in L^{2}(\Gamma)$.

Proof.

Using (2) in Lemma 2.1, we have that there exists some $q>2$ such that $S^{*}S\in\mathcal{L}(L^{2}(\Gamma),L^{q}(\Gamma))$. The slant differentibility then follows directly from [15, Pr 4.1 b)] and the semismoothness from [20, Theorem 2.13]; see also [14].

The expression for the slanting functional follows from the slant derivative of the function $\min(0,\cdot)$ and the chain rule [20, Theorem 2.10(c)].

Finally, given $z\in L^{2}(\Gamma)$

[TABLE]

This equations can be read as

[TABLE]

Let us define $w=\chi_{\Gamma_{F}(u)}v$. Taking into account the definition of $S$ and $S^{*}$, we have that $\nu v=\nu w+z\chi_{\Gamma_{A}(u)}$ and

[TABLE]

On one hand, using the continuity of $S^{*}$ (cf. (2)) and of $S$, we have

[TABLE]

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

[TABLE]

So

[TABLE]

And we get $\|w\|_{L^{2}(\Gamma)}\leq C/\nu\|z\|_{L^{2}(\Gamma)}$. Taking into account the definition of $w$ and the condition for $v$ on the active set, we get that $\|v\|_{L^{2}(\Gamma)}\leq C/\nu\|z\|_{L^{2}(\Gamma)}$, where $C$ is independent of $u$, and hence $M(u)$ has a uniformly bounded inverse for each $\nu>0$. ∎

Remark 3*.*

Notice that in the infinite dimensional case, if $\nu=0$ then $M(u)v=\chi_{\Gamma_{F}(u)}S^{*}Sv$. In practical cases, it is known (cf. (3) or [8]) that there exists $t>0$ such that $S^{*}Sv\in H^{t}(\Gamma)$ which is compactly embedded in $L^{2}(\Gamma)$, and hence $M(u)$ does not have a bounded inverse.

Given a current iterate $u\in L^{2}(\Gamma)$, we may compute next iterate $u^{+}$ using a semismooth Newton method as follows:

[TABLE]

Writing this in detail leads to

[TABLE]

which means that, if $\nu>0$,

[TABLE]

and

[TABLE]

This equations can be read as

[TABLE]

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

Theorem 4.2**.**

The semismooth Newton method described in Algorithm 3 converges q-superlinealy to $\bar{u}$ provided $u^{0}$ is close enough to $\bar{u}$ in the sense of $L^{2}(\Gamma)$.

Proof.

Once we have proved Lemma 4.1, this result is a direct consequence of [10, Theorem 3.4]; see also [15, Theorem 1.1] or [20, Theorem 2.12]. ∎

The discrete version of the variational inequality (27c) –see (36d) below– does not have a pointwise version analog to (27d). 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.

Remark 4*.*

It must noticed that, when the variational discretization of the control is used, [19, 2, 3, 18], the discrete optimal control is obtained as the pointwise projection of the discrete optimal adjoint state as in (27d), so convergence of the semismooth Newton method for such a discretization would follow from Theorem 4.2 and Lemma 4.1 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\leq\beta$ on $\Gamma$.

Condition (27c) can be replaced by the following pair of equations: there exists also $\bar{\lambda}\in L^{2}(\Gamma)$ such that equation (27d) can be written as

[TABLE]

Condition (31b) can be written as

[TABLE]

which is true for any $c>0$. Since $u\in L^{2}(\Gamma)$ appears inside the $\max$ operation, it is known that the operator described in (32) is not semismooth (following [20, Lemma 2.7]) or has not a slant derivative (as in [15, Pr 4.1]), so it is not clear that Newton’s method applied to (27a), (27b), (32) converges. Define

[TABLE]

and for each pair control-multiplier $(u,\lambda)\in L^{2}(\Gamma)\times L^{2}(\Gamma)$, set

[TABLE]

We have that

[TABLE]

So, given a current iterate $(u,\lambda)$, a Newton-like iterate to obtain $(u^{+},\lambda^{+})$ reads as

[TABLE]

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

Theorem 4.3**.**

If $c=\nu>0$, the sequence $u^{k}$ generated by the Newton-like method described in Algorithm 4 converges q-superlinearly to $\bar{u}$ in $L^{2}(\Gamma)$ provided $u^{0}$ is close enough to $\bar{u}$ in the sense of $L^{2}(\Gamma)$ and $\lambda^{0}=S^{*}y_{\Omega}-(S^{*}Su^{0}+\nu u^{0})$.

Proof.

In this case we cannot apply directly the results in [10, Theorem 3.4]; see also [15, Theorem 1.1] or [20, Theorem 2.12] since $u\in L^{2}(\Gamma)$ appears inside the $\max$ operation, and it is known that the operator described in (32) is not semismooth (following [20, Lemma 2.7]) or has not a slant derivative (as in [15, Pr 4.1]).

We will follow instead the method of proof of [15, Theorem 4.1] and we will show that the sequence $u^{k}$ generated by Algorithm 4 is exactly the same as the one generated by Algorithm 3, which we will call $\tilde{u}^{k}$.

Writing in detail (33a) we obtain

[TABLE]

Using now (33b) and (33c) we obtain

[TABLE]

Notice that (35a), (35b) and (35c) is exactly as (30a), (30b) and (30c), provided that $\Gamma_{A}(u,\lambda)=\Gamma_{A}(u)$.

To finish, let us prove by induction that $u^{k}=\tilde{u}^{k}$ and $\Gamma_{A}(u^{k},\lambda^{k})=\Gamma_{A}(\tilde{u}^{k})$ for all $k\in\mathbb{N}\cup\{0\}$. For $k=0$ it is clear from the definition of the active sets and the choice of $\lambda^{0}$ made in the assumption of this theorem. Suppose now $u^{k}=\tilde{u}^{k}$ and $\Gamma_{A}(u^{k},\lambda^{k})=\Gamma_{A}(\tilde{u}^{k})$. It is clear that $u^{k+1}=\tilde{u}^{k+1}$. On the other hand, from (35c) and (35d) and the choice $c=\nu$ we have that $cu^{k+1}+\lambda^{k+1}=-S^{*}S(u^{k+1}-y_{\Omega})$, and therefore $\Gamma_{A}(u^{k+1},\lambda^{k+1})=\Gamma_{A}(\tilde{u}^{k+1})$.

The result, hence, follows from Theorem 4.2. ∎

We want to remark here that writting the Newton step as in (35a), (35b), (35c) and (35d) we obtain exactly the same kind of algorithm as the one described in [13]

4.2 Finite dimensional approximation

Let us focus now on the finite dimensional approximation. For every $0<h<h_{0}$, there exist a unique solution $\bar{u}_{h}\in U_{h}$ of $(P_{h}^{C})$ and unique, $\bar{y}_{h}\in Y_{h}$ and $\bar{\varphi}_{h}\in Y_{h0}$ such that

[TABLE]

This way of writting the optimality system is useful to obtain error estimates (cf. [1, 2, 5]). Nevertheless, we cannot deduce a pointwise projection formula for the optimal control from the variational inequality.

Remark 5*.*

Think of this naive example. Let $\Gamma=[-1,1]$ and take a mesh with nodes $\{-1,0,1\}$. For every $u_{h},v_{h}\in U_{h}$, we have that $(u_{h},v_{h})_{\Gamma}=\bm{u}^{T}\mathcal{B}\bm{v}$, where $\mathcal{B}$ is the mass matrix

[TABLE]

Consider $\alpha=-\infty$, $\beta=0$, and take , e.g., $\bm{w}=(-2,1,1)^{T}$. Then we have that

[TABLE]

but

[TABLE]

and we get different results with the $L^{2}(\Gamma)$ projection and the pointwise projection.

To circumvent this difficulty, we rewrite conditions (36a)–(36d) in order to use standard semismooth Newton method. Taking into account (17), we can write $(P^{C}_{h})$ as

[TABLE]

where $\mathcal{A}$ and $\bm{f}$ are defined in (18) and (19) and $\bm{\alpha},\bm{\beta}\in\mathbb{R}^{N_{\mathbb{B}\times 1}}$ are the vectors whose $j$-component are respectively $\alpha(x_{j})$ and $\beta(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 $\bm{u}$ is the solution of problem (37), then there exists a unique $\bm{\lambda}\in\mathbb{R}^{N_{\mathbb{B}}\times 1}$ such that

[TABLE]

Notice that although we may name $\bar{\lambda}_{h}(x)=\sum_{j\in\mathbb{B}}\lambda_{j}e_{j}(x)$ for $x\in\Gamma$ to obtain a Lagrange multiplier $\bar{\lambda}_{h}\in U_{h}$, the complementarity condition with respect to the dot product in $\mathbb{R}^{N_{\mathbb{B}}}$, does not imply that $(\bar{\lambda}_{h},\bar{u}_{h}-\beta)_{\Gamma}=0,$ which shows again the convenience of using (38) instead of (36a)–(36d). To continue, we rewrite again the second condition in (38) to obtain

[TABLE]

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

[TABLE]

and for every pair $(\bm{u},\bm{\lambda})$, we define the sets of active and free indexes as

[TABLE]

Abusing notation and when this does not lead to confusion, we will often just write $\mathbb{A}$ and $\mathbb{F}$. Notice that $\mathbb{A}\cup\mathbb{F}=\mathbb{B}$ for every possible pair control-multiplier.

Lemma 4.4**.**

The function $F_{h,c}$ is slantly differentiable. A slanting function for $F_{h,c}$ is

[TABLE]

Finally, for all $\nu\geq 0$, the inverse of $M_{h,c}(\bm{u},\bm{\lambda})$ is uniformly bounded w.r.t $(\bm{u},\bm{\lambda})$.

Proof.

In finite dimension, the function $\max\{\bm{0},\cdot\}$ is semismooth due to Rademacher’s theorem (see e.g. [20, Example 2.4] or [15, Lemma 3.1]).

A forward computation shows that $M_{h,c}$ is a slanting function.

Let us prove the uniform boundedness of the inverse. Given $\bm{z}\in\mathbb{R}^{N_{\mathbb{B}}}$ and $\bm{\eta}\in\mathbb{R}^{N_{\mathbb{B}}}$, we have that

[TABLE]

We write the last equalities as

[TABLE]

This is the optimality system of the equality constrained optimization problem

[TABLE]

Writing $\bm{w}=\mathcal{I}_{\mathcal{B},\mathcal{F}}\bm{w}_{\mathbb{F}}+\mathcal{I}_{\mathcal{B},\mathcal{A}}\bm{w}_{\mathbb{A}}$ and taking into account the equality constraint, we have that $\bm{v}_{\mathbb{F}}$ is the solution of the following unconstrained optimization problem

[TABLE]

and therefore $\bm{v}_{\mathbb{F}}$ is the solution of the following linear system

[TABLE]

Since $\mathcal{A}$ is symmetric and positive definite, so is $\mathcal{I}_{\mathbb{F},\mathbb{B}}\mathcal{A}\mathcal{I}_{\mathbb{B},\mathbb{F}}$, and its smallest eigenvalue is bounded from below by $0<\lambda_{1}(\mathcal{M})+\nu\lambda_{1}(\mathcal{B}_{\mathbb{B},\mathbb{B}})$; see (24). Therefore, the previous system is solvable and there exists a constant $C>0$, that may depend on $h$, $\nu$ and $c$, but is independent of $\bm{u}$, such that $\|\bm{v}_{\mathbb{F}}\|\leq C(\|\bm{z}\|+\|\bm{\eta}\|)$. From this it is straight to deduce that $\|\bm{v}\|+\|\bm{\mu}\|\leq C(\|\bm{z}\|+\|\bm{\eta}\|)$ and hence $M_{h,c}$ has a uniformly bounded inverse. ∎

With these considerations, given a current iterate $(\bm{u},\bm{\lambda})$ with active and free index sets $\mathbb{A}=\mathbb{A}(\bm{u},\bm{\lambda})$ and $\mathbb{F}=\mathbb{F}(\bm{u},\bm{\lambda})$, we can compute the next step of Newton’s method $(\bm{u}^{+},\bm{\lambda}^{+})$ solving

[TABLE]

At each iteration, this is equivalent to solving the following unconstrained optimization problem in the lower-dimensional space $\mathbb{R}^{N_{\mathbb{F}}\times 1}$:

[TABLE]

Again the preconditioned gradient method works fine to solve this problem. A good preconditioner in practice is $\mathcal{P}=\mathcal{M}_{\mathbb{F},\mathbb{F}}+\nu\mathcal{B}_{\mathbb{F},\mathbb{F}}$.

Alternatively, taking into account the definition of $\mathcal{A}$ and $\bm{f}$ (see also Algorithms 1 and 2 respectively), we may write one step of the semismooth Newton algorithm as

[TABLE]

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

Remark 6*.*

Notice that the only information from iteration $k$ used to compute iteration $k+1$ is the set of active indexes. Therefore, when $\mathbb{A}(\bm{u}_{k+1},\bm{\lambda}_{k+1})=\mathbb{A}(\bm{u}_{k},\bm{\lambda}_{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 $\bm{u}^{0}$ and $\bm{\lambda}^{0}$, but only an initial guess for the active set. We include nevertheless an initial guess for the control variable in Algorithm 5 because we use $\bm{u}_{k}$ as the initial guess for the pcg to obtain $\bm{u}_{k+1}$.

The following convergence result is a direct consequence of Lemma 4.1 and [10, Theorem 3.4]; see also [15, Theorem 1.1] or [20, Theorem 2.12].

Theorem 4.5**.**

The sequence $\bm{u}_{k}$ generated by Algorithm 5 converges q-superlinearly to $\bm{u}$, the solution of (37).

Example 4.6**.**

We resume the 2D problem described in Example 3.3 and the 3D problem described in Example 3.6, with the upper constraint $\beta\equiv 0.16$. We test Algorithm 5.

Following the tip of Theorem 4.3, we have taken the parameter $c=\nu$ and $\bm{\lambda}_{0}=\bm{f}-\mathcal{A}\bm{u}_{0}$. Nevertheless, we have not been able to observe any problem for different values of $c$. The seed is set to $\bm{u}_{0}=\bm{0}$.

To solve the optimality system of the unconstrained optimization problem at each Newton iterate, we have used the preconditioned conjugate gradient method for (42) with initial guess $\mathcal{I}_{\mathbb{F},\mathbb{B}}\bm{u}_{k}$. The use of the reduced problem is even more advisable in this case than it was in Examples 3.3 and 3.6, because the size of the system (43) is $N+N_{\mathbb{I}}+N_{\mathbb{A}}$, which in any case is greater or equal than the size of the system (21), which is $N+N_{\mathbb{I}}$.

In Tables 4.6 and 4.6 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 5. 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_{h}\bar{u}_{h_{-1}}$, the solution in a coarser mesh with $h_{-1}\geq h$. This idea can be iterated with a mesh family with parameters $h_{-M},\dots,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.

Example 4.7**.**

Solving the last 2D problem in Table 4.6 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 refinemesh.

In the 3D case (see Table 4.6), 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 interp3.

Example 4.8**.**

Absence of Tikhonov parameter In the same 2D pentagonal domain, we take again $\beta\equiv 0.16$, but now $\alpha\equiv-1.2$, $\nu=0$, $y_{\Omega}\equiv 1$ if $x_{1}>0.25$, and $y_{\Omega}\equiv-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 3.2; see Table 4.8.

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\times 2^{-8}$ in Figure 1. The boundary has been stretched on a 1D line and its corners have been marked with circles on the graph of $\bar{u}_{h}$.

5 State constraints

In the rest of the work, we will suppose that $d=2$, $\Omega$ is convex and $\Gamma$ is polygonal. According to [5], if problem $(P^{S})$ admits a feasible Slater point, then it has a unique solution $\bar{u}\in H^{1/2}(\Gamma)$ and there exist $\bar{y}\in H^{1}(\Omega)\cap C(\bar{\omega})$, $\bar{\varphi}\in W^{1,t}_{0}(\Omega)$ for all $t<2$ and two nonnegative measures $\bar{\mu}^{+},\ \bar{\mu}^{-}\in\mathcal{M}(\bar{\omega})$ such that

[TABLE]

and $\operatorname{supp}\bar{\mu}^{+}\subset\{\bar{y}=b\},$ $\operatorname{supp}\bar{\mu}^{-}\subset\{\bar{y}=a\}$. In this case the adjoint state equation (44b) must be understood in the transposition sense, [5, Equation (4)], and $\langle\cdot,\cdot\rangle$ denotes the duality product between $\mathcal{M}(\bar{\omega})$ and $C(\bar{\omega})$.

As is pointed out in [12, page 197], or [15, page 878], the PDAS or semismooth Newton methods described in the previous section are not applicable to this problem, since the multiplier is a measure. In [38] a Moreau-Yosida regularization is proposed. This strategy is further investigated in [16]. Let us briefly describe this method. Again, to simplify the notation, we will do the exposition for the unilateral constraint $y\leq b$ in $\bar{\omega}$.

Given a shift function $\mu^{*}\in L^{q}(\omega)$ for some $q>2$ and a parameter $\gamma>0$, we will solve the unconstrained problem

[TABLE]

This problem has a unique solution $u^{\gamma}\in H^{1/2}(\Gamma)$ for every $\gamma>0$. Since the functional is of class $C^{1}$, first order optimality conditions can be derived directly from the work [1]: there exist unique $y^{\gamma}\in H^{1}(\Omega)$ and $\varphi^{\gamma}\in H^{1}_{0}(\Omega)$ such that

[TABLE]

Theorem 5.1**.**

The semismooth Newton method to solve (45a), (45b), (45c) converges locally $q$-superlinearly.

Proof.

We may write the system (45a), (45b), (45c) as the equation $G(u)=0$, where $G:L^{2}(\Gamma)\to L^{2}(\Gamma)$ is given by

[TABLE]

Using the regularity results in [8, 5], together with the election of $\mu^{*}$ in $L^{q}(\omega)$ we have that $\mu^{*}+\gamma(Su-b)\in L^{q}(\Omega)$ for some $q>2$, and hence $G$ is semismooth in the sense stated in [20]; see Theorems 2.13 and 2.10(c) in the aforementioned reference.

Define now ${\omega_{A}}(u)=\{x\in\bar{\omega}:\ \mu^{*}+\gamma(Su-b)>0\}$. A slant differential of $G(u)$ is given by the expression

[TABLE]

Let us see that the inverse of $M(u)$ is uniformly bounded in $\mathcal{L}(L^{2}(\Gamma),L^{2}(\Gamma))$ for all $u\in L^{2}(\Gamma)$. For any $z\in L^{2}(\Gamma)$, the equation $M(u)v=z$ can be read as

[TABLE]

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

[TABLE]

and hence

[TABLE]

So we have that $\|v\|_{L^{2}(\Gamma)}\leq\|z\|_{L^{2}(\Gamma)}$ and therefore $M(u)$ is uniformly bounded in $\mathcal{L}(L^{2}(\Gamma),L^{2}(\Gamma))$ for all $u\in L^{2}(\Gamma)$. ∎

When $\gamma\to+\infty$, $u^{\gamma}\to\bar{u}$; see [38, Theorem 3.1] or [16, Pr 2.1].

Let us turn now to the discrete problem $(P^{S}_{h})$. Consider the space $\mathcal{M}_{h}\subset\mathcal{M}(\bar{\omega})$ which is spanned by the Dirac measures corresponding to the nodes $\{x_{j}\}_{j\in\mathbb{J}}$, where $\mathbb{J}=\{j\colon x_{j}\in\bar{\omega}\}$. Following [5, Corollary 2], if the continuous problem $(P^{S})$ has a regular feasible Slater point, then $(P^{S}_{h})$ has a unique solution $\bar{u}_{h}$. Moreover there exist unique $\bar{y}_{h}\in K_{a,b}^{h}$, $\bar{\varphi}_{h}\in Y_{h,0}$ and two nonnegative measures $\bar{\mu}_{h}^{+},\bar{\mu}_{h}^{-}\in\mathcal{M}_{h}$ such that

[TABLE]

and $\operatorname{supp}\bar{\mu}_{h}^{+}\subset\{x_{j}:\,\bar{y}(x_{j})=b(x_{j})\}$, $\operatorname{supp}\bar{\mu}_{h}^{-}\subset\{x_{j}:\,\bar{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 [12]. Nevertheless, the other assumption fails here: the slant derivative may not have an inverse. Let us show this. Consider

[TABLE]

It is clear that the optimality system (46a)–(46d) is equivalent to the equation $F(\bm{u},\bm{\mu})=\bm{0}$. We define the sets of active and free nodes related to a pair control-multiplier as

[TABLE]

The slant derivative of $F(\bm{u},\bm{\mu})$ is given by

[TABLE]

Lemma 5.2**.**

The inverse of $M(\bm{u},\bm{\mu})$ may not exist for some $(\bm{u},\bm{\mu})$.

Proof.

For any $(\bm{z},\bm{\delta})\in\mathbb{R}^{N_{\mathbb{B}}}\times\mathbb{R}^{N_{\mathbb{J}}}$, we have that

[TABLE]

From the second equation we have that

[TABLE]

If $N_{\mathbb{A}}>N_{\mathbb{B}}$ we have a linear system with more equations than variables, and hence it will not be consistent for at least one value of $\bm{\delta}$ 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 $\mu^{*}\in Y_{h}$ such that $\mu^{*}(x_{j})=0$ if $j\not\in\mathbb{J}$, and a parameter $\gamma>0$, a direct discretization of $(Q^{\gamma})$ could be

[TABLE]

Problem $(\tilde{Q}^{\gamma}_{h})$ has a unique solution $u^{\gamma}_{h}\in U_{h}$ that converges to $u^{\gamma}$ as $h\to 0$. Since the functional is not of class $C^{2}$, the problem does not fit exactly into the framework of [1, 2, 5], so we will give a sketch of the proof of this fact.

Lemma 5.3**.**

Problem $(\tilde{Q}^{\gamma}_{h})$ has a unique solution $u^{\gamma}_{h}\in U_{h}$. Moreover, there exist unique $y_{h}^{\gamma}\in Y_{h}$ and $\varphi_{h}^{\gamma}\in Y_{h0}$ such that

[TABLE]

Finally, $u_{h}^{\gamma}\to u^{\gamma}$ in $L^{2}(\Gamma)$ and there exists some $\lambda>0$ such that

[TABLE]

Proof.

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}\in L^{2}(\Gamma)$ and its related states $y_{i}$ and adjoint states $\varphi_{i}$ w.r.t. problem $(Q^{\gamma})$. Using the definition of solution in the transposition sense for the states and both the state and adjoint state equations (45a) and (45b) we have (shorting $z_{i}=\mu^{*}+\gamma(y_{i}-b)$, $i=1,2$)

[TABLE]

The rest is quite standard:

[TABLE]

For fixed $\gamma$, $u_{h}^{\gamma}$ is uniformly bounded in $L^{2}(\Gamma)$ because $\|u_{h}^{\gamma}\|_{L^{2}(\Gamma)}^{2}\leq 1/\nu\|y_{\Omega}\|_{L^{2}(\Omega)}^{2}+1/\gamma/\nu\|\max\{0,\mu^{*}-\gamma b\}\|_{L^{2}(\Omega)}^{2}$, and therefore using [1, Theorem 5.7] we have that $\|-\partial_{n}^{h}\varphi_{h}^{\gamma}+\partial_{n}\varphi_{u_{h}^{\gamma}}\|\leq Ch^{\lambda}$ for some $\lambda>0$. For the first addend we obtain directly [math] using first order optimality conditions (45c) and (47c). ∎

A practical way of computing the integral in the penalty term is to use the lumped mass matrix $\mathcal{L}\in\mathbb{R}^{N\times N}$. Since we are going to use it only for integrals in $\bar{\omega}$, we define it as

[TABLE]

Therefore, we will be solving

[TABLE]

First order optimality conditions for this problem read like

[TABLE]

This nonlinear system can be solved using a semismooth Newton method. To this end, we will use the nonlinear operator $G_{h}:\mathbb{R}^{N_{\mathbb{B}}}\to\mathbb{R}^{N_{\mathbb{B}}}$ defined by

[TABLE]

Thus, solving (48a)–(48c) is equivalent to solving $G_{h}(\bm{u})=\bm{0}$.

For every $u_{h}\in U_{h}$ we define sets of active and free nodes as

[TABLE]

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. $\mathbb{A}_{\omega}$, $\mathbb{A}_{\omega}(\bm{u})$. Notice that if $\mu^{*}\equiv 0$, then $\mathbb{A}_{\omega}(u_{h},\gamma,h)$ is independent of $\gamma$. For the sake of notation, it is also convenient to define for every set $\mathbb{A}_{\omega}\subseteq\mathbb{J}$ the diagonal matrix $\mathcal{H}(\mathbb{A}_{\omega})\in\mathbb{R}^{N\times N}$ such that

[TABLE]

Abusing notation, we will often write $\mathcal{H}(\bm{u})=\mathcal{H}(\mathbb{A}_{\omega}(\bm{u}))$ or even we will write $\mathcal{H}$ when this does not lead to confusion. The proof of the following result is as the corresponding one in infinite dimension.

Theorem 5.4**.**

$G_{h}(\bm{u})$* is slantly differentiable, a slant differential is given by*

[TABLE]

and it has a uniformly bounded inverse w.r.t. $\bm{u}$.

Notice that using this $\mathcal{H}$ notation, we can write

[TABLE]

and therefore, for a given $\bm{u}$, and denoting $\mathcal{H}=\mathcal{H}(\bm{u})$, one Newton step reads like

[TABLE]

Let us comment that for the computation of $\bm{w}=\mathcal{S}^{T}\mathcal{H}\bm{y}$ for some $\bm{y}\in\mathbb{R}^{N}$, first we solve

[TABLE]

and, next, we have

[TABLE]

because $\mathcal{H}$ is diagonal and its nonzero components correspond to nodes that lie in $\bar{\omega}$, and are hence interior to $\Omega$.

Again a preconditioned conjugate gradient method can be used to solve this system, provided an efficient way of computing $\bm{d}=(\mathcal{A}+\gamma\mathcal{S}^{T}\mathcal{H}\mathcal{S})\bm{v}$ and the second member of the system. Notice that in each of the algorithms 6 and 7 the computation of $\bm{\phi}$ can be done with just one system solve.

Alternatively, the solution to (50) can be obtained solving

[TABLE]

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

Remark 7*.*

As we noticed in Remark 6, the only information from iteration $k$ used to compute iteration $k+1$ is the set of active indexes. Therefore, when $\mathbb{A}_{\omega}(\bm{u}_{k+1})=\mathbb{A}_{\omega}(\bm{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 9 and 10, the main use for the solution of $(Q_{h}^{\gamma})$ 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}^{\gamma})$ in general. Following [16], we may implement the following stopping criterion for Algorithm 8. After step 8, we compute an approximation of the multiplier $\bm{\mu}_{k+1}\in\mathbb{R}^{N}$ as

[TABLE]

and we may stop the algorithm if

[TABLE]

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 $\gamma$.

We state the convergence result for the Newton method. It follows directly from Theorem 5.4

Theorem 5.5**.**

The semismooth Newton method to solve $(Q^{\gamma}_{h})$ described in Algorithm 8 converges locally q-superlinearly.

In all the examples below we have taken the shift $\mu^{*}=0$.

Example 5.6**.**

In general it is not a good idea to solve directly $(Q_{h}^{\gamma})$ for some $\gamma$ big enough. The resulting intermediate problems are usually very ill conditioned.

We resume the 2D problem taken from [5] and described in Example 3.3. As in [5], we take $\omega$ the ball centered at $(-0.1,-0.1)$ with radius $0.2$ and $b\equiv 0.15$. We use Algorithm 8 taking $\bm{u}_{0}=\bm{0}$. The results are summarized in Table 5.6. We measure the unfeasibility of the state as the maximum constraint violation $\mathrm{mcv}(y_{h})=\|\max(y_{h}-b,0)\|_{L^{\infty}(\bar{\omega})}$.

The first thing that can be observed is that $\gamma$ cannot be too big w.r.t. the problem size. Although the algorithm should converge in finite time, a value too big for $\gamma$ 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 5.1).

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

In [16] 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 8 and superscripts for the numbering of the iterates of the continuation strategy provided in Algorithm 9. For some sequence $\{\tau^{n}\}_{n\geq 1}$ such that $\tau^{n}>1$ for all $n\geq 1$, we have

There are still three important details to be explained about Algorithm 9: first and most important, the choice of the initial guess $\bm{u}^{0}$; next, the appropriate values for $\tau^{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 $\gamma$. Since we are going to deal with meshes of different sizes, we will use the notation $\mathbb{A}_{\omega}({u}_{h},\gamma,h)$ when needed. We also recall that $I_{h}$ is the pointwise interpolation operator from $C(\bar{\Gamma})$ onto $U_{h}$. Given $(h,\gamma^{0})$, we fix $M>0$ and pick a finite sequence $(h_{n},\gamma_{n})_{n=0}^{M}$ such that $h_{n}$ is decreasing, $h_{M}=h$, $\gamma_{n}$ is increasing and $\gamma_{M}\leq\gamma^{0}$.

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

In [16] some criteria are given to choose $\tau^{n+1}>1$. In practice, if $\tau^{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 [16] 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

[TABLE]

As an alternative, a tolerance for the maximum constraint violation $e_{\infty}$ can be used stop if $\mathrm{mcv}(y_{h}^{\gamma_{n}})\leq e_{\infty}$.

With all these considerations, we propose the following algorithm. Fix a mesh sequence such that $\{h_{j}\}_{j=0}^{M}$ is decreasing, $\gamma_{0}>0$, $C>0$, $u_{h_{0}}\in U_{h_{0}}$, $n_{\max}\in\mathbb{N}$ and a sequence $\{\tau^{n}\}_{n=1}^{n_{\max}+1}$ such that $\tau^{n}>1$ for all $n$.

Example 5.7**.**

Now we apply Algorithm 11 for $\gamma_{0}=1$, $\tau^{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 5.7. The total computation time was 1570 seconds, which is a significant improvement compared with the 8246 seconds used by Algorithm 8 (see Example 5.6).

Determination of (the support of) the Lagrange multipliers of $(P^{S})$. Comparison of the adjoint state equations for $(P_{h}^{S})$ and $(Q_{h}^{\gamma})$ leads to the approximation formula

[TABLE]

Since $\bar{\omega}\subset\Omega$, 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 $\omega$, say $\bar{\mu}^{+}=\sum_{k=1}^{n}\mu_{k}\delta_{X_{k}}$ for some $n>0$ and $X_{k}\in\partial\omega$.

Increasing $\gamma$ 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.

Example 5.8**.**

In Example 5.7 the number of active nodes for the last mesh is 748. All of them are on $\partial\omega$ but they are not isolated. We can see in Table 5.8 how the number of active nodes decreases as $\gamma$ increases, with little variation of the functional or the state.

6 Control and state constraints

Again according to [5], if problem $(P^{CS})$ admits a feasible Slater point, then it has a unique solution $\bar{u}\in H^{1/2}(\Gamma)$ and there exist $\bar{y}\in H^{1}(\Omega)\cap C(\bar{\omega})$, $\bar{\varphi}\in W^{1,t}_{0}(\Omega)$ for all $t<2$ and two nonnegative measures $\bar{\mu}^{+},\ \bar{\mu}^{-}\in\mathcal{M}(\bar{\omega})$ such that

[TABLE]

and $\operatorname{supp}\bar{\mu}^{+}\subset\{\bar{y}=b\},$ $\operatorname{supp}\bar{\mu}^{-}\subset\{\bar{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\leq\beta$ on $\Gamma$ and $y\leq b$ in $\bar{\omega}$ to simplify the notation. For a shift function $\mu^{*}\in L^{q}(\omega)$ for some $q>2$ and a parameter $\gamma>0$, we consider the problem

[TABLE]

This problem has a unique solution $u^{\gamma}\in H^{1/2}(\Gamma)$. Moreover, there exist $y^{\gamma}\in H^{1}(\Omega)$, $\varphi^{\gamma}\in H^{s}(\Omega)$, $s>3/2$, such that

[TABLE]

Define

[TABLE]

It is clear that $u^{\gamma}$ is the unique solution of $(Q^{c,\gamma})$ if and only if $G^{\gamma}(u^{\gamma})=0$. For some fixed shift function $\mu^{*}\in L^{q}(\omega)$, $q>2$, consider the active sets

[TABLE]

and

[TABLE]

A slant differential of $G^{\gamma}(u)$ is given by

[TABLE]

Theorem 6.1**.**

$G^{\gamma}$* is slantly differentiable, $M^{\gamma}(u)$ is a slant differential of $G^{\gamma}$ and for every fixed $\nu>0$, $M^{\gamma}(u)$ has an inverse in $\mathcal{L}(L^{2}(\Gamma),L^{2}(\Gamma))$ uniformly bounded for all $u\in L^{2}(\Gamma)$. The semismooth Newton method $M^{\gamma}(u^{+}-u)=-G^{\gamma}(u)$ converges q-superlinearly.*

Proof.

The proof follows the same lines as those of Lemma 4.1 and Theorems 4.2 and 5.1∎

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 (54c) with the help of a Lagrange multiplier. There exists $\lambda^{\gamma}\in L^{2}(\Gamma)$ such that, for any $c>0$,

[TABLE]

Define now

[TABLE]

We have that

[TABLE]

and hence a Newton-like method to solve $F_{c}^{\gamma}(u,\lambda)=0$ is given by

[TABLE]

Although $F_{c}^{\gamma}$ is known not to be slantly differentiable, it can be proved as in Theorem 4.3 that the sequence generated by the above described Newton-like method to solve $F_{c}^{\gamma}(u,\lambda)=0$ is the same than the one generated by the semismooth Newton method to solve $G^{\gamma}(u)=0$ provided we take the same initial guess $u_{0}$, $c=\nu$ and $\lambda_{0}=S^{*}(y_{\Omega}-(1+\gamma\chi_{\omega_{\mathbb{A}(u_{0},\gamma)}})Su_{0})-\nu u_{0}$.

Let us turn now to the finite dimensional problem. We can write the approximation of problem $(Q^{C,\gamma}_{h})$ as a constrained optimization problem in $\mathbb{R}^{N_{\mathbb{B}},1}$:

[TABLE]

Existence and uniqueness of solution of this problem, as well as error estimates for the difference $\|u^{\gamma}-u^{\gamma}_{h}\|_{L^{2}(\Gamma)}$ can be proved as we did for $(\tilde{Q}^{\gamma})$ in Lemma 5.3. Since we have control constraints, instead of the point-wise interpolation, the Casas-Raymond interpolate [1, Equation (7.9)] should be used in the last step.

First order optimality conditions read as

[TABLE]

Using the definitions (40) and (49) for the active sets of indexes $\mathbb{A}(\bm{u},\bm{\lambda})$ and $\mathbb{A}_{\omega}(\bm{u},\gamma,h)$ and the matrix $\mathcal{H}$ related to $\mathbb{A}_{\omega}$, we have that one step of Newton’s method can be written as

[TABLE]

or, alternatively, as

[TABLE]

An adaptation of algorithms 5 and 8 to solve $(Q^{C,\gamma}_{h})$ is straightforward, and so is an adaptation of Algorithm 11 to use a continuation strategy together with a nested mesh strategy.

Example 6.2**.**

We repeat Example 5.7 adding the control constraint $u\leq 0.16$. We obtain the results summarized in Table 6.2. 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.