A posteriori error estimates for a bang-bang optimal control problem
Francisco Fuica





Facultad de Matemáticas, Pontificia Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, Santiago, Chile.

francisco.fuica@uc.cl


[2010]Primary 
49M25,		   
65N15,         
65N30.         








§ ABSTRACT 
We propose and analyze a posteriori error estimates for a control-constrained optimal control problem with bang-bang solutions. 
We consider a solution strategy based on the variational approach, where the control variable is not discretized; no Tikhonov regularization is made. 
We design, for the proposed scheme, a residual-type a posteriori error estimator that can be decomposed as the sum of two individual contributions related to the discretization of the state and adjoint equations.
We explore reliability and efficiency properties of the aforementioned error estimator.
We illustrate the theory with numerical examples.


A posteriori error estimates for a bang-bang optimal control problem
    Francisco Fuica
    
====================================================================














§ INTRODUCTION

The main purpose of this work is the design and analysis of a posteriori error estimates for a bang-bang optimal control problem governed by an elliptic partial differential equation (PDE) as state equation. 
This PDE-constrained optimization problem entails the minimization of a cost functional in which the cost of the control is negligible, so no Tikhonov regularization term is considered.
To make matters precise, let Ω⊂ℝ^d (d∈{2, 3}) be an open, bounded and convex polygonal/polyhedral domain with boundary ∂Ω. 
Given a desired state y_Ω∈ L^2(Ω), we define the cost functional

    J(u):=1/2y_u - y_Ω_Ω^2.
 
We will consider the following optimal control problem:

    min J(u)   subject to   -Δ y_u  =  u   in Ω,   
    y  =  0   on ∂Ω,    and   u∈𝕌_ad,
 
where 𝕌_ad:={v ∈ L^2(Ω):  a ≤ v(x) ≤ b  a.e. x ∈Ω} with a < b.

Problem (<ref>) has been considered in several works.
One of the main sources of difficulty for deriving error estimates within this type of problems is that the cost function J does not incorporate the standard Tikhonov regularization term αu_Ω^2 with α > 0.
Note that without this term we cannot directly derive nor bound the standard error term αu̅ - u̅_h_Ω^2, where u̅ denotes the optimal control and u̅_h denotes a suitable approximation of u̅.
To the best of our knowledge, the first work that provides an a priori error analysis for problem (<ref>) is <cit.>.
In such a work, the authors used the so-called variational approach, introduced in <cit.>, in order to discretize problem (<ref>).
In addition, they proved estimates for the approximation error associated to the optimal state and adjoint state without assuming that the control is of bang-bang type; see <cit.>.
For the case when the optimal control is of bang-bang type, the authors proved an error estimate for all the optimal variables <cit.>.
A suitable Tikhonov regularization of problem (<ref>) and its convergence were studied in <cit.>, under an additional assumption on the structure of the optimal control.
The parabolic counterpart of (<ref>) was studied in <cit.>, where the authors proved, using Petrov-Galerkin schemes in time and conforming finite elements in space, a priori estimates for the error between a discretized regularized problem and the limit problem.
In the particular case of bang-bang controls, the estimates were further improved; see <cit.>. 
We also mention the work <cit.>, in which the authors considered the Tikhonov regularization of problem (<ref>) with a semilinear elliptic PDE and derived a priori regularization error estimates for the control; a suitable extension to sparse optimal control problems was derived as well.
Finally, for second-order analysis and approximation results for bang-bang optimal control problems, we refer the reader to <cit.>.

Among the different numerical methods that exist in the literature to approximate solutions to PDE-constrained optimization problems (and PDEs in general), a particular class stands out for its competitive performance, improving the quality of discrete approximations of the corresponding problem within a prescribed tolerance using a minimal amount of work. 
These are the adaptive finite element methods (AFEMs).
The main tools present in these iterative methods are a posteriori error estimates, which provide global and local information on the error of discrete solutions and that can be easily computed from the given numerical solution and problem data.
Regarding the use of these methods in the context of control–constrained linear–quadratic optimal control problems, several advances has been made in recent years.
For a discussion on this matter, we refer the interested reader to the following non-exhaustive list of works: <cit.>.
As opposed to these advances, the analysis of AFEMs for bang-bang optimal control problems is rather scarce. 
To the best of our knowledge, the work <cit.> appears to be the first one that provides a posteriori error estimates for problem (<ref>). 
In this article, the author investigated Tikhonov regularization and discretization of bang-bang control problems, developing a parameter choice rule that adaptively selects the Tikhonov regularization parameter depending on a posteriori computable quantities. 
However, the error estimates were not robust with respect to α.
We also mention the work <cit.>, where robust global reliability estimates were provided.
We note that no efficiency estimates were provided in <cit.>.

In the present manuscript, we consider the variational discretization <cit.> to approximate the optimal control problem (<ref>).
In particular, we use piecewise linear functions to approximate the solution of the state and adjoint equations whereas the admissible control set is not discretized.
To perform the analysis, we follow <cit.> and do not consider a Tikhonov regularization.
This approach allows us to circumvent the necessity of choosing  a suitable regularization parameter  for each mesh, cf. <cit.>.
 Within this framework, we devise a residual–based a posteriori error estimator that is formed by only two contributions that are related to the discretization of the state and adjoint equations.
In two- and three-dimensional convex polygonal/polyhedral domains, we prove efficiency estimates; reliability properties of the a posteriori error estimator are studied as well.
More precisely, we prove that the corresponding local error indicators associated to the discretization of the state and adjoint equations are locally efficient; see Theorem <ref>.
We recall that reliability estimates are sufficient to obtain a numerical solution with an accuracy below a prescribed tolerance, whereas efficiency estimates are of importance since they ensure that the mesh is correctly refined, so that one obtains a numerical solution with a prescribed tolerance using a (nearly) minimal number of degrees of freedom <cit.>.
Based on the proposed a posteriori error estimator, we design a simple adaptive loop that delivers optimal experimental rates of convergence for all the involved individual contributions of the corresponding error. 
In particular, the aforementioned loop delivers quadratic rates of convergence for the approximation error associated to all the optimal variables. 
In addition, and in contrast to <cit.>, the error indicator that we consider for the adjoint variable in the L^∞(Ω)–norm allows for unbounded forcing terms. 
This is of importance since, as it can be observed from (<ref>), the discrete adjoint equation has y̅_ℓ(𝔲̅_ℓ) - y_Ω as a forcing term and, in general, the latter term does not necessarily belong to L^∞(Ω). 


The rest of the manuscript is organized as follows. 
The remaining of this section is devoted to introduce the notation that we will use throughout the manuscript.
In section <ref> we present a weak formulation for the optimal control problem under consideration and introduce a finite element discretization scheme.
The main part of the paper is section <ref>, where we design an a posteriori error estimator for the proposed approximation scheme and analyze reliability and efficiency properties. 
Finally, in section <ref>, we present a series of two-dimensional numerical examples that illustrate the theory and reveal a competitive performance of the devised AFEMs.






















 §.§ Notation


Throughout this work, we use standard notation for Lebesgue and Sobolev spaces and their corresponding norms. 
Given an open and bounded domain G, we denote by (·,·)_G and ·_G the inner product and norm of L^2(G), respectively.
If 𝒳 and 𝒴 are Banach function spaces, we write 𝒳↪𝒴 to denote that 𝒳 is continuously embedded in 𝒴. 
The relation 𝔞≲𝔟 indicates that 𝔞≤ C 𝔟, with a positive constant that depends neither on 𝔞, 𝔟 nor on the discretization parameters. 
The value of C might change at each occurrence.












§ THE OPTIMAL CONTROL PROBLEM



In this section, we briefly present a weak formulation for problem (<ref>) and recall first-order optimality conditions. 
We also introduce a finite element discretization scheme.



 §.§ Weak formulation


We consider the following weak version of problem (<ref>): Find

    min{ J(u): u ∈𝕌_ad}

subject to 

    y_u∈ H_0^1(Ω)  :  (∇ y_u,∇ v)_Ω=(u,v)_Ω  ∀ v ∈ H_0^1(Ω).

Problem (<ref>)–(<ref>) admits a unique optimal solution u̅∈𝕌_ad <cit.> and the optimal control u̅ satisfies the first-order optimality condition <cit.> 

    (p̅, u - u̅)_Ω≥ 0   ∀ u ∈𝕌_ad,

where the optimal adjoint state p̅∈ H_0^1(Ω) solves the following adjoint equation:

    (∇ v,∇p̅)_Ω=(y̅ - y_Ω,v)_Ω  ∀ v ∈ H_0^1(Ω),

with y̅ := y_u̅.
Moreover, inequality (<ref>) implies that, for a.e. x∈Ω, we have 

    u̅(x) = a   if  p̅(x) > 0,   u̅(x) ∈ [a,b]   if  p̅(x) = 0,   u̅(x) = b    if  p̅(x) < 0;

see <cit.>.







 §.§ Finite element approximation


Let us introduce some ingredients of standard finite element approximations <cit.>. 
Let 𝒯 = {T} be a conforming partition of Ω into simplices T with size h_T := diam(T). 
Let us denote by 𝕋 the collection of conforming and shape regular meshes that are refinements of 𝒯_0, where 𝒯_0 represents an initial mesh. 
Given a mesh 𝒯_ℓ∈𝕋 with ℓ∈ℕ_0, we denote by ℰ_ℓ the set of internal (d-1)-dimensional interelement boundaries e of 𝒯_ℓ. 
For T ∈𝒯_ℓ, we let ℰ_T denote the subset of ℰ_ℓ which contains the sides of the element T. 
We denote by 𝒩_e ⊂𝒯_ℓ the subset that contains the two elements that have e as a side, namely, 𝒩_e={T^+,T^-}, where T^+, T^- ∈𝒯_ℓ are such that e = T^+ ∩ T^-. 
For T ∈𝒯_ℓ, we define the star associated with the element T as

    𝒩_T:= { T^'∈𝒯_ℓ: ℰ_T∩ℰ_T^'≠∅}.

In an abuse of notation, below we denote by 𝒩_T either the set itself or the union of its elements. 

Given a mesh 𝒯_ℓ∈𝕋 with ℓ∈ℕ_0, we define the finite element space of continuous piecewise polynomials of degree one that vanish on the boundary as

    𝕍_ℓ = {v_ℓ∈ C(Ω): v_ℓ|_T∈ℙ_1(T)  ∀ T∈𝒯_ℓ}∩ H_0^1(Ω).


Given a discrete function v_ℓ∈𝕍_ℓ, we define, for any internal side e ∈ℰ_ℓ, the jump or interelement residual ∇ v_ℓ·𝐧 on e by

    ∇ v_ℓ·𝐧|_e:= 𝐧^+·∇ v_ℓ|^_T^+ + 𝐧^-·∇ v_ℓ|^_T^-,

where 𝐧^± denote the unit exterior normal to the element T^±. Here, T^+, T^-∈_ℓ are such that T^+≠ T^- and ∂ T^+∩∂ T^- = e.





We consider the following semi-discrete version of the optimal control problem (<ref>)–(<ref>): Find 

    min{ J(𝔲): 𝔲∈𝕌_ad}

subject to the discrete state equation

    (∇ y_ℓ(𝔲),∇ v_ℓ)_Ω
    =
    (𝔲, v_ℓ)_Ω    ∀ v_ℓ∈𝕍_ℓ.

Existence of an optimal control 𝔲̅ for (<ref>)–(<ref>) can be proved by standard arguments. However, uniqueness of 𝔲̅ is not guaranteed <cit.>. In spite of this fact, we can characterize optimal solutions: every optimal control 𝔲̅ for (<ref>)–(<ref>) satisfies 

    (p̅_ℓ, u - 𝔲̅)_Ω≥ 0   ∀ u ∈𝕌_ad,

where p̅_ℓ∈𝕍_ℓ solves the discrete adjoint equation

    (∇ v_ℓ,∇p̅_ℓ)_Ω=(y̅_ℓ(𝔲̅) - y_Ω,v_ℓ)_Ω  ∀ v_ℓ∈𝕍_ℓ.

Here, y̅_ℓ(𝔲̅) solves problem (<ref>) with 𝔲 = 𝔲̅. We immediately notice that y̅_ℓ(𝔲̅) and p̅_ℓ are unique, even if 𝔲̅ is not unique; see <cit.>.

As in the continuous case, from (<ref>) it stems the following characterization for optimal controls 𝔲̅, for a.e. x∈Ω: 

    𝔲̅(x) = a   if  p̅_ℓ(x) > 0,   𝔲̅(x) ∈ [a,b]   if  p̅_ℓ(x) = 0,   𝔲̅(x) = b    if  p̅_ℓ(x) < 0.

From (<ref>) it follows that, if p̅_ℓ admits a zero level set of measure 0, then 𝔲̅(x)=a or 𝔲̅(x)=b for a.e. x∈Ω and thus 𝔲̅ is both unique and of bang-bang type.

Finally, since 𝔲̅ implicitly depends on 𝒯_ℓ, in what follows we shall adopt the notation 𝔲̅_ℓ.












§ A POSTERIORI ERROR ESTIMATES


In the present section we design an error estimator for optimal control problem (<ref>)–(<ref>).
We explore its reliability properties in section <ref>, whereas its local efficiency is proved in section <ref>.








 §.§ Reliability

The upcoming analysis mainly hinges on approximations of the error between a solution to the discrete optimal control problem and suitable auxiliary variables that we shall define in what follows.








  §.§.§ Auxiliary upper bounds


Let 𝔲̅_ℓ be a solution of the semi-discrete optimal control problem associated to a mesh 𝒯_ℓ∈𝕋. We introduce the auxiliary variable y_𝔲̅_ℓ∈ H_0^1(Ω), defined as the unique solution to 

    (∇ y_𝔲̅_ℓ,∇ v)_Ω
    =
    (𝔲̅_ℓ, v)_Ω  ∀ v ∈ H_0^1(Ω).

We note that the discrete optimal state y̅_ℓ(𝔲̅_ℓ) corresponds to the finite element approximation of y_𝔲̅_ℓ in 𝕍_ℓ. 
Hence, since we have assumed that Ω is convex, we can use the results from <cit.> to obtain that

    y_𝔲̅_ℓ - y̅_ℓ(𝔲̅_ℓ)_Ω≲η_st,2,

where the error estimator η_st,2 and its local error indicators are defined by

    η_st,2^2:=∑_T∈_ℓE_st,T^2, 
        
    E_st,T^2:=h_T^4𝔲̅_ℓ_T^2 + ∑_e∈ℰ_Th_T^3∇y̅_ℓ(𝔲̅_ℓ) ·𝐧|_e_e^2.


We define p_y̅_ℓ∈ H_0^1(Ω) as the unique solution to 

    (∇ v, ∇ p_y̅_ℓ)_Ω
    =
    (y̅_ℓ(𝔲̅_ℓ) - y_Ω, v)_Ω  ∀ v ∈ H_0^1(Ω),

and immediately note that p̅_ℓ corresponds to the finite element approximation of p_y̅_ℓ in 𝕍_ℓ. Using again that Ω is convex, we invoke <cit.> to obtain that 

    p_y̅_ℓ - p̅_ℓ_Ω≲η_adj,2,

where the error estimator η_adj,2 and its corresponding local indicators are given by 

    η_adj,2^2:=∑_T∈_ℓ E_adj,T^2,     E_adj,T^2:=h_T^4y̅_ℓ(𝔲̅_ℓ) - y_Ω_T^2 + ∑_e∈ℰ_Th_T^3∇p̅_ℓ·𝐧|_e_e^2.


We introduce, for each T∈_ℓ, the following a posteriori local error indicators 

    E_adj,∞,T:= h_T^2-d/2y̅_ℓ(𝔲̅_ℓ) - y_Ω_T + h_Tmax_e∈ℰ_T∇p̅_ℓ·𝐧|_e_L^∞(e),

and the error estimator η_adj,∞:=max_T∈𝒯_ℓE_adj,∞,T. The proof of the following reliability estimate can be found in <cit.>:

    p_y̅_ℓ - p̅_ℓ_L^∞(Ω)≲ι_ℓη_adj,∞,

where the term ι_ℓ is defined by

    ι_ℓ:= |log(max_T∈𝒯_ℓ1/h_T)|.

An important feature of the local error indicator in (<ref>) is that it incorporates the L^2(T)-norm of the element residual instead of the L^∞(T)-norm, which is the common consideration in the literature; see, e.g., <cit.>.
 In particular, the error estimator η_adj,∞ allows for a pointwise a posteriori error analysis with unbounded right-hand sides, cf. <cit.>. 
 This is of importance since, as it can be observed in (<ref>), its indicators contain the term  y̅_ℓ(𝔲̅_ℓ) - y_Ω which is not necessarily bounded.










  §.§.§ Reliability estimates


Let us define the set S:={x∈Ω : p̅(x) ≠ 0}⊂Ω, i.e., the support of p̅. 
In what follows we shall assume that (cf. <cit.>)

    ∃β∈ (0,1],  ∃ℭ > 0,  ∀ε > 0   such that   |{x∈ S : |p̅(x)|≤ε}|≤ℭε^β,

where, given a set A⊂Ω, |A| denotes the Lebesgue measure of it.



Let u̅∈𝕌_ad be the unique solution to problem (<ref>)–(<ref>) with y̅ and p̅ being the corresponding state and adjoint state variables, respectively. Let 𝔲̅_ℓ∈𝕌_ad be a solution to the semi-discrete problem with y̅_ℓ(𝔲̅_ℓ) and p̅_ℓ being the corresponding discrete state and discrete adjoint state variables, respectively. If assumption (<ref>) holds, then

    y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2 + p̅ - p̅_ℓ_L^∞(Ω)^2≲η_st,2^2 + η_adj,2 + ι_ℓ^β + 1η_adj,∞^β + 1 + ι_ℓ^2/2 - βη_adj,∞^2/2 - β ,

and

    u̅ - 𝔲̅_ℓ_L^1(S)^2≲η_st,2^2β + η_adj,2^β + ι_ℓ^β(β + 1)η_adj,∞^β(β + 1) + ι_ℓ^2β/2 - βη_adj,∞^2β/2 - β .

The hidden constants are independent of the continuous and discrete optimal variables, the size of the elements in the mesh 𝒯_ℓ, and #𝒯_ℓ.


We invoke assumption (<ref>) and the same arguments that lead to inequality (2.13) in <cit.> to obtain 

    u̅ - 𝔲̅_ℓ_L^1(S)≲p̅ - p̅_ℓ_L^∞(Ω)^β.


We now proceed on the basis of three steps.

Step 1. (estimation of y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω) We consider u=𝔲̅_ℓ in (<ref>) and u=u̅ in (<ref>). 
Adding the resulting inequalities and using the auxiliary term p_y̅_ℓ, defined as the solution to (<ref>), we arrive at 

    0 ≤ (p̅ - p̅_ℓ, 𝔲̅_ℓ - u̅)_Ω = (p_y̅_ℓ - p̅_ℓ, 𝔲̅_ℓ - u̅)_Ω + (p̅ - p_y̅_ℓ, 𝔲̅_ℓ - u̅)_Ω =: 𝐈 + 𝐈𝐈.

We now bound the terms 𝐈 and 𝐈𝐈 in (<ref>). 
To estimate 𝐈, we use (<ref>) and the a posteriori error estimates (<ref>) and (<ref>):

    𝐈   = (p_y̅_ℓ - p̅_ℓ,𝔲̅_ℓ - u̅)_Ω∖ S + (p_y̅_ℓ - p̅_ℓ,𝔲̅_ℓ - u̅)_S
       ≲p_y̅_ℓ - p̅_ℓ_Ω + p_y̅_ℓ - p̅_ℓ_L^∞(Ω)𝔲̅_ℓ - u̅_L^1(S)
       ≲η_adj,2 + ι_ℓη_adj,∞p̅ - p̅_ℓ_L^∞(Ω)^β.

To control the term p̅ - p̅_ℓ_L^∞(Ω)^β in (<ref>), we use the triangle inequality and the a posteriori error estimate (<ref>) to arrive at

    p̅ - p̅_ℓ_L^∞(Ω)^β≲    p̅ - p_y̅_ℓ_L^∞(Ω)^β + p_y̅_ℓ - p̅_ℓ_L^∞(Ω)^β
    ≲    p̅ - p_y̅_ℓ_L^∞(Ω)^β + ι_ℓ^βη_adj,∞^β.

Now, we note that p̅ - p_y̅_ℓ∈ H_0^1(Ω) solves

    (∇ v, ∇ (p̅ - p_y̅_ℓ))_Ω  = (y̅ - y̅_ℓ(𝔲̅_ℓ), v)_Ω  ∀ v∈ H_0^1(Ω).

Hence, the convexity of Ω, the continuous embedding H^2(Ω) ↪ C(Ω), and the stability of problem (<ref>) yield the bound p̅ - p_y̅_ℓ_L^∞(Ω)^β≲p̅ - p_y̅_ℓ_H^2(Ω)^β≲y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^β. Using the latter in (<ref>) results in

    p̅ - p̅_ℓ_L^∞(Ω)^β≲y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^β + ι_ℓ^βη_adj,∞^β.

Thus, the use (<ref>) and Young's inequality in (<ref>) with exponents 𝗉=2β and 𝗊=22-β, yield

    𝐈≤ C(η_adj,2 + ι_ℓ^β + 1η_adj,∞^β + 1 + ι_ℓ^2/2 - βη_adj,∞^2/2 - β) + 1/4y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2,    C > 0.

We notice that if ι_ℓ>1 and η_adj,∞<1, then 

    ι_ℓ^β + 1η_adj,∞^β + 1 + ι_ℓ^2/2 - βη_adj,∞^2/2 - β≤ 2 ι_ℓ^β + 1η_adj,∞^2/2-β≤ 2ι_ℓ^2η_adj,∞^2/2-β,

due to the fact that 2/2 - β≤β + 1 ≤  2 when β∈(0,1].

To estimate 𝐈𝐈, we first note that y_𝔲̅_ℓ - y̅∈ H_0^1(Ω) solves

    (∇ (y_𝔲̅_ℓ - y̅), ∇ v)_Ω  = (𝔲̅_ℓ - u̅, v)_Ω  ∀ v∈ H_0^1(Ω).

Then, replacing v = p̅ - p_y̅_ℓ in (<ref>) and v = y_𝔲̅_ℓ - y̅ in (<ref>) we obtain the identity (𝔲̅_ℓ - u̅, p̅ - p_y̅_ℓ)_Ω =  (y̅ - y̅_ℓ(𝔲̅_ℓ),  y_𝔲̅_ℓ - y̅)_Ω, which, in turns, yields 

    𝐈𝐈 = (y̅ - y̅_ℓ(𝔲̅_ℓ),  y_𝔲̅_ℓ - y̅_ℓ(𝔲̅_ℓ))_Ω - y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2.

This identity, in light of Young's inequality and estimate (<ref>), allows us to obtain 

    𝐈𝐈≤ Cη_st,2^2 - 1/2y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2,      C>0.


We have thus proved, in view of (<ref>), (<ref>), and (<ref>), the  error bound

    y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2≲η_st,2^2 + η_adj,2 + ι_ℓ^β + 1η_adj,∞^β + 1 + ι_ℓ^2/2 - βη_adj,∞^2/2 - β.



Step 2. (estimation of p̅ - p̅_ℓ_L^∞(Ω)) The use of triangle inequality, a posteriori error estimate (<ref>), and the bound p̅ - p_y̅_ℓ_L^∞(Ω)≲y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω results in 

    p̅ - p̅_ℓ_L^∞(Ω)^2≲p̅ - p_y̅_ℓ_L^∞(Ω)^2 + p_y̅_ℓ - p̅_ℓ_L^∞(Ω)^2≲y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2 + ι_ℓ^2η_adj,∞^2.

Consequently, using (<ref>) in the latter estimate we conclude that

    p̅ - p̅_ℓ_L^∞(Ω)^2≲η_st,2^2 +  η_adj,2 + ι_ℓ^β + 1η_adj,∞^β + 1 + ι_ℓ^2/2 - βη_adj,∞^2/2 - β.



Step 3. (estimation of u̅ - 𝔲̅_ℓ_L^1(S)) To obtain the estimate (<ref>) it suffices to use the bound (<ref>) in (<ref>).



When β=1 in (<ref>), we obtain the following simplified upper bound for the total approximation error:

    u̅ - 𝔲̅_ℓ_L^1(S)^2 + y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2 + p̅ - p̅_ℓ_L^∞(Ω)^2≲η_st,2^2 +  ι_ℓ^2η_adj,∞^2 + η_adj,2.

A sufficient condition to ensure that assumption (<ref>) is fulfilled with β=1 was given in <cit.>.









 §.§ Efficiency

In this section, we study efficiency properties for the local a posteriori error estimators η_st,2, η_adj,2, and η_adj,∞, defined in section <ref>. 
Before proceeding with the analysis, we introduce some notation: for an edge, triangle or tetrahedron G, let 𝒱(G) be the set of vertices of G. We also introduce, given T ∈_ℓ and e ∈ℰ_T, the standard interior and edge bubble functions φ_T and φ_e, respectively; see, e.g., <cit.>.

Let T ∈_ℓ and e ∈ℰ_T. Recall that 𝒩_e denotes the patch composed by the two elements T^+ and T^- sharing e. We introduce the following edge/face bubble function 

    ψ_e|_𝒩_e=d^4d(∏_∈𝒱(e)ϕ_^T^+ϕ_^T^-)^2,

where, for ∈𝒱(e), ϕ_^T^± denotes the barycentric coordinates of T^±, respectively, which are understood as functions over 𝒩_e. 
The bubble function ψ_e has the following properties: ψ_e∈ℙ_4d(𝒩_e), ψ_e∈ C^2(𝒩_e), and ψ_e = 0 on ∂𝒩_e. 
In addition, it satisfies

    ∇ψ_e = 0  on ∂𝒩_e,   ∇ψ_e·𝐧= 0  on  e.


Given T∈𝒯_ℓ, let Π_T : L^2(T) →ℙ_0(T) be the orthogonal projection operator into constant functions over T. In particular, for v∈ L^2(Ω), we have Π_Tv:=1|T|∫_Tv and Π_T v_T≤v_T.


With all these ingredients at hand, we are ready to prove the local efficiency of the aforementioned error estimators. We start with η_st,2, defined in (<ref>).


Let u̅∈𝕌_ad be the unique solution to problem (<ref>)–(<ref>) with y̅ being the corresponding optimal state. Let 𝔲̅_ℓ∈𝕌_ad be a solution to the semi-discrete problem with y̅_ℓ(𝔲̅_ℓ) being the corresponding discrete state variable. 
Then, for T∈𝒯_ℓ, the local error indicator E_st,T, defined as in (<ref>), satisfies

    E_st,T^2
    ≲y̅ - y̅_ℓ(𝔲̅_ℓ)_𝒩_T^2
    +
    h_T^4 - du̅-𝔲̅_ℓ_L^1(𝒩_T)^2
    +
    ∑_T'∈𝒩_Th_T^4𝔲̅_ℓ - Π_T𝔲̅_ℓ_T'^2,

where 𝒩_T is defined as in (<ref>) and the hidden constant is independent of continuous and discrete optimal variables, the size of the elements in the mesh 𝒯_ℓ, and #𝒯_ℓ.

 Let v ∈ H_0^1(Ω) be such that v|_T∈ C^2(T) for all T∈𝒯_ℓ. Use v as a test function in (<ref>) and apply elementwise integration by parts to obtain

    ( ∇ (y̅ - y̅_ℓ(𝔲̅_ℓ)), ∇ v)_Ω -  (u̅ - 𝔲̅_ℓ,v)_Ω 
    =
    ∑_T∈_ℓ ( 𝔲̅_ℓ, v )_T  - ∑_e∈ℰ_ℓ(∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e,v)_e.

At the same time, we integrate by parts, again, to arrive at

    (∇ (y̅ - y̅_ℓ(𝔲̅_ℓ)), ∇ v)_Ω = ∑_e∈ℰ_ℓ (∇ v·𝐧|_e, y̅-y̅_ℓ(𝔲̅_ℓ))_e - ∑_T∈_ℓ (y̅-y̅_ℓ(𝔲̅_ℓ), Δ v)_T.

Hence, combining the two previous identities we obtain, for all v ∈ H_0^1(Ω)  such that v|_T∈ C^2(T) for all T∈𝒯_ℓ, the equality

    ∑_e∈ℰ_ℓ (∇ v·𝐧|_e, y̅-y̅_ℓ(𝔲̅_ℓ))_e  - ∑_T∈_ℓ (y̅-y̅_ℓ(𝔲̅_ℓ), Δ v)_T  -  (u̅ - 𝔲̅_ℓ,v)_Ω
        =
    ∑_T∈_ℓ[ (Π_T𝔲̅_ℓ, v )_T + (𝔲̅_ℓ - Π_T𝔲̅_ℓ, v )_T]  - ∑_e∈ℰ_ℓ ( ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e, v)_e.


We now proceed on the basis of two steps.

Step 1. (estimation of h_T^4𝔲̅_ℓ_T^2) Let T∈_ℓ. 
An application of the triangle inequality gives

    h_T^4𝔲̅_ℓ_T^2≲
    h_T^4Π_T𝔲̅_ℓ_T^2 + h_T^4𝔲̅_ℓ - Π_T𝔲̅_ℓ_T^2.

To estimate h_T^4Π_T𝔲̅_ℓ_T^2, we evaluate v = φ_T^2Π_T𝔲̅_ℓ in (<ref>). 
Then, we use that ∇ (φ^2Π_T𝔲̅_ℓ)=0 on ∂ T, standard properties of the bubble function φ_T, and the inverse estimate φ_T^2Π_T𝔲̅_ℓ_L^∞(T)≲ h_T^-d/2φ_T^2Π_T𝔲̅_ℓ_T. 
These arguments yield

    Π_T𝔲̅_ℓ_T^2≲    y̅-y̅_ℓ(𝔲̅_ℓ)_TΔ (φ_T^2Π_T𝔲̅_ℓ)_T  
        + (h_T^-d/2u̅ - 𝔲̅_ℓ_L^1(T) + 𝔲̅_ℓ - Π_T𝔲̅_ℓ_T)Π_T𝔲̅_ℓ_T.

Using that Δ (φ_T^2Π_T𝔲̅_ℓ) = Π_T(𝔲̅_ℓ)Δφ_T^2 in combination with properties of φ_T it follows that Δ (φ_T^2Π_T𝔲̅_ℓ)_T≲ h_T^-2Π_T𝔲̅_ℓ_T. 
Thus, we obtain 

    h_T^4Π_T(𝔲̅_ℓ)_T^2≲y̅-y̅_ℓ(𝔲̅_ℓ)_T^2 + h_T^4-du̅ - 𝔲̅_ℓ_L^1(T)^2 + h_T^4𝔲̅_ℓ - Π_T𝔲̅_ℓ_T^2.

A combination of (<ref>) and (<ref>) concludes the estimation.


Step 2. (estimation of h_T^3∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e_e^2) 
Let T∈_ℓ and e∈ℰ_T.  
In what follows, we extend the jump term ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e, defined on e, to the patch 𝒩_e. 
Such extension can be done, for example, by using a continuation operator as in <cit.>.
Hereinafter  we make no distinction between the jump term and its extension.

We invoke the bubble function ψ_e, evaluate v = ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_eψ_e in (<ref>)  and use that ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e∈ℝ, ψ_e ∈ H^2_0(𝒩_e), and ∇ψ_e·𝐧|_e = 0. 
These, the fact Δ v =  ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_eΔψ_e since ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e∈ℝ, and basic inequalities imply that

    ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_eψ_e^1/2_e^2≲   ∑_T' ∈𝒩_e( h_T'^-2y̅-y̅_ℓ(𝔲̅_ℓ)_T'  +𝔲̅_ℓ_T'. 
       . +    h_T'^-d/2u̅ - 𝔲̅_ℓ_L^1(T')) h_T^1/2∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e_e,

upon using the bound ∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_eψ_e_T'≲ h_T^1/2∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e_e.

 With these estimates at hand, we thus use standard bubble functions arguments, estimate (<ref>), and the shape regularity property of the family {𝒯_ℓ} to arrive at

    h_T^3/2∇y̅_ℓ(𝔲̅_ℓ)·𝐧|_e_e≲∑_T'∈𝒩_e  
    (
    y̅-y̅_ℓ(𝔲̅_ℓ)_T' + h_T^2𝔲̅_ℓ_T' + h_T^2-d/2u̅ - 𝔲̅_ℓ_L^1(T')
    ),

which concludes the proof.


We now continue with the study of the local efficiency properties of the estimator η_adj,2 defined in (<ref>).


Let u̅∈𝕌_ad be the unique solution to problem (<ref>)–(<ref>) with y̅ and p̅ being the corresponding optimal state and adjoint state, respectively. 
Let 𝔲̅_ℓ∈𝕌_ad be a solution to the semi-discrete problem with y̅_ℓ(𝔲̅_ℓ) and p̅_ℓ being the corresponding discrete state and adjoint state variables, respectively.  
Then, for T∈𝒯_ℓ, the local error indicator E_adj,T, defined as in (<ref>), satisfies that

    E_adj,T^2
    ≲
    h_T^dp̅-p̅_ℓ_L^∞(𝒩_T)^2 + h_T^4y̅-y̅_ℓ(𝔲̅_ℓ)_𝒩_T^2 + ∑_T'∈𝒩_Th_T^4y_Ω - Π_T(y_Ω)_T'^2,

where 𝒩_T is defined as in (<ref>) and the hidden constant is independent of continuous and discrete optimal variables, the size of the elements in the mesh 𝒯_ℓ, and #𝒯_ℓ.


Similar arguments to the ones that lead to (<ref>) yield, for every v ∈ H_0^1(Ω)  such that v|_T∈ C^2(T) (T∈𝒯_ℓ), the identity

    ∑_e∈ℰ_ℓ (∇ v·𝐧|_e, p̅-p̅_ℓ)_e - ∑_T∈_ℓ (p̅-p̅_ℓ, Δ v)_T  -  (y̅ - y̅_ℓ(𝔲̅_ℓ),v)_Ω
       =
    ∑_T∈_ℓ( Π_Ty_Ω - y_Ω,v)_T +
    ∑_T∈_ℓ ( y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω, v )_T  - ∑_e∈ℰ_ℓ ( ∇p̅_ℓ·𝐧|_e, v)_e.

We now estimate the terms in (<ref>) on the basis of two steps.

Step 1. (estimation of h_T^4y̅_ℓ(𝔲̅_ℓ) - y_Ω_T^2) Let T∈_ℓ. The use of the triangle inequality yields

    h_T^4y̅_ℓ(𝔲̅_ℓ) - y_Ω_T^2≲
    h_T^4y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2  
    +
    h_T^4Π_Ty_Ω - y_Ω_T^2.

We recall that Π_T denotes the orthogonal projection operator into constant functions over T. 
To control h_T^4y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2 in the above inequality, we invoke the interior bubble function φ_T and  choose v=φ_T^2(y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω) in (<ref>). Hence, using standard properties of φ_T and the inequality  Δ(φ_T^2(y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω))_L^1(T)≲ 
h_T^d/2-2y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T, which follows from the inverse estimate <cit.>, we obtain

    y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T≲Π_Ty_Ω - y_Ω_T + y̅ - y̅_ℓ(𝔲̅_ℓ)_T + h_T^d/2-2p̅ - p̅_ℓ_L^∞(T),

from which we conclude that

    h_T^4y̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2≲    
    h_T^4Π_Ty_Ω - y_Ω_T^2 
        + h_T^4y̅ - y̅_ℓ(𝔲̅_ℓ)_T^2 + h_T^dp̅ - p̅_ℓ_L^∞(T)^2.

The use of the latter in (<ref>) results in the desired bound.

Step 2. (estimation of h_T^3∇p̅_ℓ·𝐧|_e_e^2) 
Let T∈_ℓ and e∈ℰ_T. 
As in Lemma <ref>, we extend the jump term ∇p̅_ℓ·𝐧|_e to the patch 𝒩_e, making no distinction between the jump term and its extension.
We choose v = ∇p̅_ℓ·𝐧|_eψ_e in (<ref>)  and use that ∇p̅_ℓ·𝐧|_e∈ℝ, ψ_e ∈ H^2_0(𝒩_e), and ∇ψ_e·𝐧|_e = 0. 
Consequently, basic inequalities imply that

    ∇p̅_ℓ·𝐧|_e_e^2≲    ∑_T'∈𝒩_e(h_T^d/2-2p̅ - p̅_ℓ_L^∞(T')
        + y̅ - y̅_ℓ(𝔲̅_ℓ)_T' + y̅_ℓ(𝔲̅_ℓ) - y_Ω_T')ψ_e∇p̅_ℓ·𝐧|_e_T'.

Therefore, using that ψ_e∇p̅_ℓ·𝐧|_e_T'≲ h_T^1/2∇p̅_ℓ·𝐧|_e_e in combination with estimate (<ref>), we conclude that 

    h_T^3∇p̅_ℓ·𝐧|_e_e^2≲   ∑_T'∈𝒩_e(h_T^dp̅ - p̅_ℓ_L^∞(T')^2
         +  h_T^4y̅ - y̅_ℓ(𝔲̅_ℓ)_T'^2 + h_T^4Π_Ty_Ω - y_Ω_T'^2),

which ends the proof.


We now continue with the study of the local efficiency properties of the estimator η_adj,2 defined in (<ref>).


In the framework of Lemma <ref> we have that, for T∈𝒯_ℓ, the local error indicator E_adj,∞,T, defined as in (<ref>), satisfies 

    E_adj,T,∞^2
    ≲p̅-p̅_ℓ_L^∞(𝒩_T)^2 + h_T^4 - dy̅-y̅_ℓ(𝔲̅_ℓ)_𝒩_T^2 + ∑_T'∈𝒩_Th_T^4-dy_Ω - Π_Ty_Ω_T'^2,

where 𝒩_T is defined as in (<ref>) and the hidden constant is independent of continuous and discrete optimal variables, the size of the elements in the mesh 𝒯_ℓ, and #𝒯_ℓ.


We proceed on the basis of two steps.

Step 1. (estimation of h_T^4-dy̅_ℓ(𝔲̅_ℓ) - y_Ω_T^2) Let T∈_ℓ. To estimate this term in (<ref>), we first use a triangle inequality and obtain

    h_T^4-dy̅_ℓ(𝔲̅_ℓ) - y_Ω_T^2≤ 
    h_T^4-dy̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2
    +
    h_T^4-dΠ_Ty_Ω - y_Ω_T^2.

The estimate for h_T^4-dy̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2 stems from (<ref>) after multiplying by h^-d_T:

    h_T^4-dy̅_ℓ(𝔲̅_ℓ) - Π_Ty_Ω_T^2≲p̅ - p̅_ℓ_L^∞(T)^2 + h_T^4-dy̅ - y̅_ℓ(𝔲̅_ℓ)_T^2 + h_T^4-dΠ_Ty_Ω - y_Ω_T^2.


Step 2. (estimation of h_T^2∇p̅_ℓ·𝐧|_e_L^∞(e)^2)  
Let T∈_ℓ and e∈ℰ_T. 
Since ∇p̅_ℓ·𝐧|_e∈ℝ, we deduce that 

    ∇p̅_ℓ·𝐧|_e_L^∞(e)^2 
    =
    |∇p̅_ℓ·𝐧|_e|^2
    =
    |e|^-1∇p̅_ℓ·𝐧|_e_e^2,

where |e| denotes the measure of e. 
In view of the shape regularity of the mesh _ℓ we have that |e|≈ h_T^d-1 and consequently h_T^2+d∇p̅_ℓ·𝐧|_e_L^∞(e)^2≈
h_T^3 ∇p̅_ℓ·𝐧|_e_e^2. 
Hence, using the latter in estimate (<ref>) and multiplying by h_T^-d, we arrive at

    h_T^2∇p̅_ℓ·𝐧|_e_L^∞(e)^2≲    ∑_T'∈𝒩_e(p̅ - p̅_ℓ_L^∞(T')^2
        +  h_T^4-dy̅ - y̅_ℓ(𝔲̅_ℓ)_T'^2 + h_T^4-dΠ_Ty_Ω - y_Ω_T'^2),

 which concludes the proof.


The next result is a direct consequence of Lemmas <ref>, <ref> and <ref>.


In the framework of Lemma <ref> we have, for T∈𝒯_ℓ, that

    E_st,T^2 + E_adj,T^2 + 
    E_adj,T,∞^2 
    ≲ 
    (1 + h_T^4 + h_T^4-d)y̅ - y̅_ℓ(𝔲̅_ℓ)_𝒩_T^2
    + (1 + h_T^d)p̅-p̅_ℓ_L^∞(𝒩_T)^2 
    
     + h_T^4 - du̅-𝔲̅_ℓ_L^1(𝒩_T)^2  + 
    ∑_T'∈𝒩_Th_T^4𝔲̅_ℓ - Π_T𝔲̅_ℓ_T'^2  + ∑_T'∈𝒩_T(h_T^4 + h_T^4-d)y_Ω - Π_T(y_Ω)_T'^2,

where 𝒩_T is defined as in (<ref>) and the hidden constant is independent of continuous and discrete optimal variables, the size of the elements in the mesh 𝒯_ℓ, and #𝒯_ℓ.










§ NUMERICAL EXAMPLES

In this section we conduct a series of numerical examples in 2D that support our theoretical findings and illustrate the performance of the error estimator 

    E^2
    :=
    η_st,2^2 +  η_adj,∞^2 + η_adj,2^2

that we proposed and analyzed in section <ref>.
We immediately note that we have considered the instance when β = 1 (cf. Remark <ref>) in the total error estimator E. 
We shall consider this configuration for E in all the following numerical examples.
 In order to simplify the construction of exact optimal solutions, we have incorporated an extra forcing term f∈ L^2(Ω) in the state equation. 
 The right-hand side of the state equation reads now as follows: (f  + u̅, v)_Ω.
Moreover, in sections <ref> and <ref> below, we go beyond the presented theory and perform numerical experiments with a non-convex domain. 
The considered numerical examples have been carried out with the help of a code that we implemented using . 
When assembling all system matrices we have used exact integration whereas right-hand sides, approximation errors, and error indicators are computed by a quadrature formula which is exact for polynomials of degree 19.


For a given partition _ℓ, we seek y̅_ℓ(𝔲̅_ℓ) ∈𝕍_ℓ,  p̅_ℓ∈𝕍_ℓ, and 𝔲̅_ℓ∈ U_ad that solve (<ref>), (<ref>), and (<ref>).
We solve such a nonlinear system of equations by using the solution technique devised in <cit.>. 
Once a discrete solution is obtained, we compute the error indicator 

    E_T := (E_st,T^2 + E_adj,T^2 + E_adj,T,∞^2)^1/2

to drive the adaptive procedure described in Algorithm <ref>.
For the numerical results, we define the total number of degrees of freedom Ndofs = 2 dim𝕍_ℓ.
Finally, we define the effectivity index 

    ℐ_eff:= E/(u̅ - 𝔲̅_ℓ_L^1(Ω)^2 + y̅ - y̅_ℓ(𝔲̅_ℓ)_Ω^2 + p̅ - p̅_ℓ_L^∞(Ω)^2)^1/2.








 §.§ Exact solution on convex domain

Following <cit.> (see also <cit.>) we set Ω:=(0,1)^2, a=-1, b=1, and take f and y_Ω such that 

    y̅(x_1,x_2) = sin(π x_1)sin(π x_2), 
       p̅(x_1,x_2) = -sin(2π x_1)sin(2π x_2)/(8π^2),   u̅(x_1,x_2) = -sign(p̅(x_1,x_2))

for  (x_1,x_2)∈Ω.





In Figures <ref> and <ref> we display the results obtained for this example. 
We show, in Fig. <ref>, experimental rates of convergence for each contribution of the total error when uniform and adaptive refinement are considered.
We observe that all the approximation errors obtained for both schemes exhibit optimal rates of convergence.
In Fig. <ref> we show experimental rates of convergence for all the individual contributions of the error estimator E (see (<ref>)) and the effectivity index, when adaptive refinement is considered.
When the total number of degrees of freedom increases, we observe that the effectivity index stabilizes around the values 3 and 5.



 §.§ Exact solution on non-convex domain

We set Ω=(-1,1)^2∖[0,1)×(-1,0], a=-1, b=1, and take f and y_Ω such that the exact optimal state and adjoint state are given, in polar coordinates (ρ,ω) with ω∈[0,3π/2], by

    y̅(ρ,ω)   =sin(π(ρsin(ω)+1)/2)sin(π(ρcos(ω)+1)/2)ρ^2/3sin(2ω/3), 
    
    p̅(ρ,ω)   = (0.5-ρ)y̅(ρ,ω).

The purpose of this example is to investigate the performance of the devised a posteriori error estimator when we violate the convexity assumption considered on the domain.  







We present  the results obtained for this example in Figures <ref>, <ref> and <ref>.
In Fig. <ref>, we display experimental rates of convergence for each contribution of the total error when uniform and adaptive refinement are considered.
We observe that the designed adaptive procedure outperforms uniform refinement.
In particular, it exhibits optimal rates of convergence for each contribution of the total error.
In Fig. <ref> we show experimental rates of convergence for all the individual contributions of the error estimator E and the effectivity index, when adaptive refinement is considered.
We observe that the effectivity seems to stabilize around the values 3 and 4 when the total number of degrees of freedom increases.
Finally, in Fig. <ref>, we present the finite element approximations y̅_ℓ(𝔲̅_ℓ) and p̅_ℓ, and adaptively refined meshes. 
It can be observed that the refinement is being concentrated on the re-entrant corner (0, 0).




 §.§ Unknown solution on non-convex domain


We set Ω=(-1,1)^2∖[0,1)×[0,1), a=-0.5, b=0.5, and data

    y_Ω(x_1,x_2) = 1/√(x_1^2 + x_2^2) - 10sin(x_1x_2),      (x_1,x_2)∈Ω.

We note that y_Ω∈ L^2(Ω) ∖ L^∞(Ω). 



In Fig. <ref> we show the results obtained for this example. 
Similar conclusions to the ones presented for the example from section <ref> can be derived. 
Particularly, we observe optimal experimental rates of convergence for all the individual contributions of the error estimator E within the adaptive loop.









siam

















