Mesh Adaptivity in Optimal Controlof Elliptic Variational Inequalities with Point-Tracking of the State

An adaptive finite element method is developed for a class of optimal control problems with elliptic variational inequality constraints and objective functionals defined on the space of continuous functions, necessitated by a point-tracking requirement with respect to the state variable. A suitable first order stationarity concept is derived for the problem class via a penalty technique. The dual-weighted residual approach for goal-oriented adaptive finite elements is applied and relies on the stationarity system. It yields primal residuals weighted by approxi- mate dual quantities and vice versa as well as complementarity mismatch errors. A report on numerical tests, including the critical case of biactivity, completes this work.

2010 Mathematics Subject Classification: [2010]49K20, 65K15, 90C33

Keywords: adaptive finite element methods, optimal control of variational inequalities in func- tion space, point-tracking, C-stationarity, goal-oriented error estimation.

1

Mesh adaptivity in optimal control of elliptic variational

inequalities with point-tracking of the state

Charles Brett

MASDOC CDT, Mathematics Institute, University of Warwick, Coventry, CV4 7AL, UK

c.brett@warwick.ac.uk Charles M. Elliott Mathematics Institute,

University of Warwick, Coventry CV4 7AL, UK C.M.Elliott@warwick.ac.uk

Michael Hinterm¨uller

Department of Mathematics, Humboldt-Universit¨at zu Berlin, Unter den Linden 6, D-10099 Berlin, Germany

hint@math.hu-berlin.de Caroline L¨obhard

Department of Mathematics, Humboldt-Universit¨at zu Berlin, Unter den Linden 6, D-10099 Berlin, Germany

loebhard@math.hu-berlin.de

Abstract

The final publication is available at European Mathematical Society Publishing House:

DOI 10.4171/IFB/332

© European Mathematical Society 2015

1 Introduction

In this paper we study a class of optimal control problems with variational inequality and pointwise control constraints, as well as an objective containing the tracking of the state variable at specific points. Such models are of interest, for example, in inverse problems where a parameter has to be identified from measurements at predefined locations in the underlying domain; see [1] for an application in mathematical finance.

Since variational inequalities often model an equilibrium condition, the problems treated here are called mathematical programs with equilibrium constraints (MPECs) in the literature. This problem class has been studied since the 1970s in finite dimensional spaces and later also in function space;

see [38, 43] and [6, 41] for extensive accounts of the progress in the respective subjects. Despite the research efforts ever since, the problem class considered in this paper has not yet been analyzed.

The degeneracy of the constraint set which is typical for MPECs makes standard methods for the derivation of necessary optimality conditions inapplicable. As an alternative, penalty methods (see e.g. [6] and the references therein) and methods relying on generalized differentiability concepts (see e.g. [40]) have been used. Concerning the state system of the present work, useful analytic properties (such as stability and directional differentiability) of solution operators of variational inequalities are discussed in [15, 34, 45]. These properties were employed successfully to derive alternative forms of stationarity. For such concepts in function space we refer to [28, 42], and to [47]

for finite dimensional problems.

Besides their value for the design of mesh independent solution algorithms, function space sta- tionarity conditions can be employed to derive goal oriented a posteriori error indicators with respect to the objective functional in the context of finite element discretizations. In optimization with partial differential equation constraints, this concept was pioneered in [10, 4], and an associated adaptive method for an MPEC is given in [27]. The method relies on local error indicators for the imple- mentation of an automatic adaption of the discrete space. For different approaches to adaptive finite element methods (AFEMs) we refer to the monographs [49, 4, 44, 3] and the references therein. The typical dual-weighted residual error estimators have been successfully applied to optimal control problems with partial differential constraints e.g. in [37, 50, 24, 23, 11, 25, 26, 46]. Alternatively, functional error estimation for optimal control problems is analyzed in [20], and residual based es- timators are studied in [37, 36, 31, 30, 33, 32, 19]. In particular, a residual based estimator for the optimal control of an obstacle problem is suggested in [19].

In the present paper, we analyze the new problem class and derive suitable stationarity condi- tions with a smooth penalty approach and an averaging technique for the function evaluations in the objective. Some of the analytic complications associated with the underlying problem class are due to the point tracking which requires a state space which embeds into continuous functions, and thus leads to reduced regularity of adjoint states. Based on the resulting so-called-almost C-stationarity system, we develop a goal-oriented indicator for the error in the objective functional of the optimal control problem which contains dual weighted residuals and terms covering the mismatch in com-

plementarity of primal variables and the respective multipliers. These are similar to the ones in [27], but with extra terms concerning the supplementary control constraints.

The paper is structured as follows: We provide a collection of basic notation in the rest of this section. In Section 2 we state the considered MPEC and prove a regularity result for the solution of the variational inequality. In addition, the penalization scheme with a new averaging technique for the objective is defined. Section 3 shows the limiting process that takes the first order optimality conditions of the auxiliary problem to the-almost C-stationarity conditions for the MPEC. These conditions are employed for the derivation of an abstract error representation formula and the a pos- teriori error indicators in Section 4. We finish in Section 5 with a brief description of the solution method and document numerical results.

Notation

For a Banach spaceY and its dualY^∗, the dual pairing of elements y^∗ ∈ Y^∗ and y ∈Y is written as⟨y^∗, y⟩Y^∗ =y^∗(y). In a reflexive Banach space Y with the identificationi ∶ Y →Y^∗∗, we write

⟨y, y^∗⟩^Y ∶= ⟨i(y), y^∗⟩^Y∗∗. If X is a Hilbert space, we write (x₁, x₂)^X for the scalar product of x₁, x₂ ∈ X. For sequences(y_k)k∈N ⊂ Y and (y_k^∗)k∈N ⊂ Y^∗, we denote the strong convergence to a limity∈Y (i.e.,∥y_k−y∥^Y →0fork → ∞) byy_k →y, whereas weak and weak^∗convergence are denoted byy_k⇀y(i.e.,∀z^∗∈Y^∗,⟨z^∗, y_k−y⟩^Y∗→0) andy_k^∗⇀^∗ y^∗(i.e.,∀z∈Y,⟨y^∗_k−y^∗, z⟩^Y∗ →0).

The standard Lebesgue an Sobolev spaces on an open, bounded Lipschitz domainΩare denoted in the usual way and we use the norms ∥v∥H₀¹(Ω) ∶= ∥∇v∥L²(Ω) and ∥v∥W₀^1,q(Ω) ∶= ∥∇v∥L^q(Ω). Fur- thermore we writeC⁰(Ω¯)for continuous functions on the closure of Ωand identify its dual space (C⁰(Ω¯))^∗with the space of bounded Borel measuresMb(Ω¯)by the Riesz-Alexandroff representa- tion theorem.

For a functionψ ∈L²(Ω), we define(ψ)+=max(0, ψ)in a pointwise almost everywhere (a.e.) sense. The characteristic function of a subsetM ⊂Ωis denoted byχ_M ∶Ω→ {0,1}, while∣M∣is the Lebesgue measure ofM providedM is measurable. We writeR^>0= {t∈R∣t>0}and, for a subset M of a Banach spaceY,M⁺= {y^∗ ∈Y^∗∣ ∀y∈M ∶ y^∗(y) ≥0} ⊂Y^∗. The number of elements in a finite setI is denoted bycard(I).

Throughout the text,C >0is a generic constant which depends only on the problem input data, but not on solutions, or on the discretization.

2 Model problem and the penalization scheme

2.1 Statement of the problem

For an open bounded domain Ω ⊂ R², we consider a_ij ∈ L^∞(Ω) (i, j ∈ {1,2}) such that for all ζ= (ζ1, ζ2) ∈R²and almost allx∈Ω,

∑2 i,j=1

a_ij(x)ζ_iζ_j ≥α∣ζ∣²,

whereα>0, and we define the uniformly coercive operatorA∶H₀¹(Ω) →H⁻¹(Ω)by A= − ∑²

i,j=1

∂

∂x_ia_ij ∂

∂x_j = −div((a_ij)∇). (2.1)

Here and below, the coefficient matrix in A is denoted by (a_ij) ∈ L^∞(Ω)^2×2. Let V = H₀¹(Ω), A ∶ V → V^∗ and K = {z ∈ V ∣z ≥ 0a.e. inΩ}. By the Lions-Stampacchia theorem (cf. [45, Thm.4:3.1]), the variational inequality problem

Findy∈K such that ⟨Ay, z−y⟩H⁻¹(Ω)≥ ⟨f+u, z−y⟩H⁻¹(Ω) ∀z∈K, (2.2) has a unique solution for all u, f ∈ V^∗. This solution can be equivalently characterized by the following complementarity system with slack variableξ,

ξ=Ay−f−u∈V^∗, y≥0inV, ξ≥0inV^∗, ⟨ξ, y⟩H⁻¹(Ω)=0. (2.3) In a slight misuse of notation, we denote the solution operator of (2.2) byy(⋅), so this mapsu∈V^∗ to y = y(u). It is well-known from [13] (see also [45, Thm.5:2.5 (i)]) that for ∂Ω ∈ C¹, a_ij ∈ L^∞(Ω) ∀i, j∈ {1,2}andf+u∈W^−1,q(Ω), there existsqˆ>2such that for all2<q<qˆthe solution of (2.2) satisfiesy(u) ∈W₀^1,q(Ω). The following proposition shows that in fact, less regularity of the boundary∂Ωis necessary to prove that the solution operator of (2.2) mapsL²(Ω)intoW₀^1,q(Ω). Proposition 2.1. Assume thatΩis of classC^0,1(i.e. has a Lipschitz boundary, which includes non- convex polygonal domains),a_ij ∈L^∞(Ω) ∀i, j∈ {1,2}andf+u∈L²(Ω). Then there existsQ>2 such that for all2<q<Qthe solution of (2.2) satisfiesy(u) ∈W₀^1,q(Ω)and the following estimate holds,

∥y(u)∥W₀^1,q(Ω)≤C∥f +u∥^W−1,q(Ω). (2.4) Proof. In the same way as in [35, IV, Theorem 2.3], we defineϑ_τ ∶R→ [0,1]forτ >0as

ϑτ(t) =⎧⎪⎪⎪

⎨⎪⎪⎪⎩

1, t≤0, 1−τ^t, 0<t<τ,

0, t≥τ

and approximate the variational inequality by the problem

Find y∈H₀¹(Ω) such that Ay=max(0,−f−u)ϑ_τ(y) +f+u. (2.5) Lemma 2.2 in [35, IV] yields a unique solution y_τ ∈ H₀¹(Ω). The Lipschitz domain Ωis regular in the sense of [22, Def. 2], and thus by [22, Thm. 3] there existsqˆ> 2such that the duality map

−∆∶H₀¹(Ω) →H⁻¹(Ω)mapsW₀^1,ˆq(Ω)ontoW^−1,ˆq(Ω)and∥φ∥W₀^1,q(Ω)≤C∥ −∆φ∥^W−1,q(Ω)for all φ∈W₀^1,q(Ω). So from Meyer’s estimate (see [39, Thm.1]) we obtainQ∈ (2,qˆ], depending only on Ω, the ellipticity constant and the norm of (a_ij), such that AmapsW₀^1,q(Ω)ontoW^−1,q(Ω)for all 2<q<Q, and we have the respective norm estimate. The right hand side of (2.5) is inW^−1,q(Ω)so the solution satisfiesy_τ ∈W₀^1,q(Ω)and

∥y_τ∥W₀^1,q(Ω)≤C∥Ay_τ∥W^−1,q(Ω)≤C∥max(0,−f−u)ϑ_τ(y_τ) +f+u∥W^−1,q(Ω)

≤C∥f+u∥^W−1,q(Ω).

Thus for a sequence τ_k → 0, (y_τk)^k∈N is uniformly bounded in W₀^1,q(Ω) and therefore contains a weakly convergent subsequence with limity¯∈W₀^1,q(Ω). Following the arguments in [35], we find that the limit satisfies the variational inequality, i.e.y¯=y(u). The uniqueness of this solution implies the weak convergence of the full sequence. The weak lower semi-continuity of the norm implies

∥y(u)∥W₀^1,q(Ω)≤lim inf{∥y_τk∥W₀^1,q(Ω)∣k∈N} ≤C∥f+u∥^W−1,q(Ω).

In the rest of this paper we assume thatΩis a Lipschitz domain inR² and thatq∈ (2, Q), with Q>2from Proposition 2.1. In this way,W₀^1,q(Ω)embeds intoC(Ω¯), and thus an optimal control problem with objectiveJ ∶C(Ω¯) ×L²(Ω) →Rcan be treated. SinceL²(Ω) ≅L²(Ω)^∗⊂W^−1,q(Ω), we can consider the solution operator as a mappingy∶L²(Ω) →W₀^1,q(Ω).

We analyze optimal control problems with point-tracking of the state variableyin the objective, and with the constraint (2.2) as well as optionalL²(Ω)-box-constraints on the control variable u, i.e.,

Minimize J(y, u) = 1

2∑

w∈I

(y(w) −y_w)²+ν

2∥u∥²L²(Ω) (2.6a)

over (y, u) ∈W₀^1,q(Ω) ×L²(Ω), (2.6b)

subject to (s.t.) y=y(u)solves (2.2), (2.6c)

u∈U_ad= {u∈L²(Ω) ∣a≤u≤ba.e. inΩ}, (2.6d) wherewdenotes a tracking point in a finite set I ⊂Ω, y_w ∈R is a desired (or measured) value of the statey at w, ν > 0 is the cost of the control, f ∈ L²(Ω), and the bounds a, b ∈ L²(Ω) satisfy a < b almost everywhere in Ω, or a = −∞ and b = ∞. The function evaluation of y ∈ W₀^1,q(Ω) inw is also denoted by ⟨δw, y⟩W^−1,q′(Ω), where δw ∈ W^−1,q′(Ω) = (W₀^1,q(Ω))^∗ with q^′ = _q−1^q . The

objective functional J ∶ W₀^1,q(Ω) ×L²(Ω) → Ris continuous and convex, and thus weakly lower semi-continuous. In its reduced form, the control problem reads

MinimizeJˆ(u) = 1 2∑

w∈I

∣y(u)(w) −y_w∣²+ν

2∥u∥²L²(Ω) over u∈U_ad. (2.7) Remark 2.2. The differential operatorAis a topological isomorphism fromW₀^1,q(Ω)intoW^−1,q(Ω) according to the proof of Proposition 2.1. We denote the adjoint ofAbyA^∗∶W₀^1,q′(Ω) →W^−1,q′(Ω), and a ∶ W₀^1,q(Ω) ×W₀^1,q′(Ω) → R is the bilinear form induced by A. Note that ˆa ∶ W₀^1,q(Ω) → R,ˆa(v) ∶=a(v, v)is convex and continuous, and thus weakly lower semi-continuous.

Remark 2.3. In the case of increased regularity, a, b, f ∈ L^p(Ω) for p ≥ 2, [15, Thm.II.1] yields thatAy∈L^p(Ω)and∥Ay∥^Lp(Ω)can be bounded by∥f+u∥^Lp(Ω)(see [15, Rem.II.2]). In particular, the slack variable ξ satisfies ξ = Ay−u−f ∈ L^p(Ω). Furthermore, ifu+f ∈ L^∞(Ω), then [48, Rem.2.6] yields H¨older continuity of y ∶ L²(Ω) → W₀^1,q(Ω) (or as a mapping to a more regular space, depending on the domainΩ) with exponent ¹₂.

The following lemma serves as a tool to prove solvability of the model problem (2.6).

Lemma 2.4. Letf ∈L²(Ω). Ifu_k ⇀u˜inL²(Ω), then(y(u_k))^k∈Nhas a subsequence converging weakly inW₀^1,q(Ω)toy(u˜).

Proof. The sequence (u_k)^k∈N is bounded owing to its weak convergence. We use the estimate in Proposition 2.1 to obtain the uniform boundedness of(y(u_k))^k∈N in W₀^1,q(Ω). Next we note that (y_k)k∈Nhas a subsequence which converges weakly inW₀^1,q(Ω)to a limity˜∈K∩W₀^1,q(Ω)due to the weak closedness ofK. To complete the proof, insert the elements of this subsequence into the variational inequality (2.2) with controlukto show that the limit satisfiesy˜=y(u˜).

Lemma 2.4 and a standard infimizing sequence (Weierstraß) argument now prove the following result.

Proposition 2.5. Problem (2.7) has a solution.

Foru∈U_adof (2.7) we definey=y(u)andξ=Ay−u−f ∈L²(Ω). The set{x∈Ω∣y(x) =0}is calledactive set, its complement{x∈Ω∣y(x) >0}is calledinactive set, and the set{x∈Ω∣ξ(x) = 0a.e. andy(x) =0}is calledbiactive set.

For the derivation of a stationarity system for (2.6) we next study a penalized version of the problem.

2.2 Penalized variational inequality

The variational inequality (2.2) can be approximated by the semi-linear partial differential equation Ay−γmax(0,−y) =u+f inH⁻¹(Ω), (2.8)

which is the first order optimality condition of the problem Minimize 1

2⟨Ay, y⟩H⁻¹(Ω)− (u+f, y)L²(Ω)+γ

2∥max(0,−y)∥²L²(Ω) overy∈H₀¹(Ω). (2.9) Above,γ >0is given andmax(0,⋅)is understood in the pointwise a.e. sense. Equation (2.8) has a unique solutiony_γ^max(u) ∈H₀¹(Ω). According to [21, Thm. 3.1] it holds for allu∈L²(Ω)that

y^max_γ (u) →y(u) inH₀¹(Ω)asγ→ ∞. (2.10) In a smoothed version of (2.8), themax(0,⋅)-operator is approximated by aC¹-functionmax(0,⋅), which depends on a parameter>0, such that

max(0,⋅) →max(0,⋅)pointwise a.e., as→0.

For now we consider the local variant from [29, Eq. (2.4)], i.e.

max(0, t) =⎧⎪⎪⎪

⎨⎪⎪⎪⎩

t ift≥,

t²

4 +₂^t+₄ ift∈ (−, ), 0 ift≤ −.

(2.11)

The smoothed penalized constraint then reads

Ay−γmax(0,−y) =u+f inH⁻¹(Ω). (2.12) Proposition 2.6(Solution operator of the semi-linear equation (2.12)). Forγ, >0,ϕ∈L²(Ω), the linearH₀¹(Ω)-elliptic partial differential operatorA as defined in (2.1) andmax from (2.11), the equation

Ay−γmax(0,−y) =ϕ (2.13)

admits a unique solution y˜ ∈ W₀^1,q(Ω). Furthermore there exists a constant C > 0 such that the following estimate holds true:

max{γ∥max(0,−y˜)∥^L2(Ω),∥y˜∥W₀^1,q(Ω)} ≤C(∥ϕ∥^L2(Ω)+γ), from some constantC=C(Ω) >0which is independent ofγand.

Proof. To analyze the solvability of the semi-linear equation (2.12), we defineψ∶H₀¹(Ω) →Rvia ψ(y) = ∫_Ω∫₀^y(x)max(0,−s)ds dx.

This mapping is monotone and bounded from below, and its derivative reads d

dyψ(y) ⋅φ= − ∫_Ωmax(0,−y)φ dx.

Now equation (2.12) can be interpreted as the optimality system of the problem MinimizeF(y) ∶=1

2⟨Ay, y⟩^H−1(Ω)− ⟨ϕ, y⟩^H−1(Ω)+γψ(y)overy∈H₀¹(Ω), (2.14) which has a unique solutiony, because˜ F is strictly convex and coercive in H₀¹(Ω). For all ≥ 0, it holds thatmax(0,−y˜) ∈ H¹(Ω). Therefore the functiony˜can be interpreted as a solution to the linear elliptic equationAy =f˜with the right hand side f˜=γmax(0,−y˜) +ϕ∈L²(Ω). According to the proof of Proposition 2.1 this yields thaty˜∈W₀^1,q(Ω)and

∥y˜∥W₀^1,q(Ω)=∥A⁻¹(γmax(0,−y˜) +ϕ)∥W₀^1,q(Ω)≤C(∥γmax(0,−y˜)∥L²(Ω)+ ∥ϕ∥L²(Ω)), (2.15) whereA⁻¹ denotes the solution operator associated with Av= f , v˜ ∈ W₀^1,q(Ω). Next, we establish the bound on ∥γmax(0,−y˜)∥L²(Ω). The semi-linear partial differential equation (2.13) holds in W^−1,q(Ω)and, owing tof˜=γmax(0,−y˜) +ϕ∈L²(Ω), even inL²(Ω).

Sincemax(0,−y˜) ∈H₀¹(Ω)and applying Green’s theorem we get from (2.13) γ∥max(0,−y˜)∥²L²(Ω) ≤ γ∫_Ωmax(0,−y˜)max(0,−y˜)dx

= − ∫_Ωg(x)∇y˜(x)^T(aij(x))∇y˜(x)dx− ∫_Ωϕmax(0,−y˜)dx, where we also use that0≤max(0,−y˜) ≤max(0,−y˜)a.e. inΩ. Here,0≤g(x) ∈∂max(0,−y˜(x)) for a.e. x ∈ Ω, where∂ denotes the subdifferential from convex analysis, cf. [14, Prop.6.45]. The positive definiteness of(a_ij)a.e. inΩand H¨older’s inequality then yield

γ∥max(0,−y˜)∥L²(Ω)≤ ∥ϕ∥L²(Ω). (2.16) Next we observe that

∥max(0,−y˜)∥L²(Ω) ≤ ∥max(0,−y˜)∥L²(Ω)+ ∥max(0,−y˜) −max(0,−y˜)∥L²(Ω)

≤ γ⁻¹∥ϕ∥L²(Ω)+ 4∣Ω∣^1/2. From (2.15) we consequently obtain

∥y˜∥W₀^1,q(Ω)≤C(∥ϕ∥L²(Ω)+γ) from some constantC=C(Ω) >0which is independent ofγand.

As in [29] we invoke the following assumption onγ and, and we denote the solution operator of (2.12) with=(γ)byy_γ∶L²(Ω) →H₀¹(Ω).

Assumption 2.7. For eachγ>0, let(γ) >0be given such thatlimγ→∞γ(γ) =0.

Theorem 2.8. Let (uk)^k∈Nbe a sequence in Uad converging weakly in L²(Ω)to u˜ ∈Uad. Assume furthermore that (γ_k)^k∈N ⊂ R^>0 withγ_k → ∞ and let _k =(γ_k) satisfy Assumption 2.7. By y_k = y_γk(u_k) and y˜ = y(u˜) we denote the solution of the smoothed penalized equation (2.12) and the variational inequality (2.2), respectively. For eachk ∈Nwe define ξ_k ∶= γ_kmax_k(0,−y_k) and set ξ˜∶= A˜y−u˜−f. Then there exists a subsequence of (ξ_k)^k∈N and (y_k)^k∈N (denoted the same) such that fork→ ∞,

ξ_k⇀ξ˜ inL²(Ω), y_γk(u_k) →y(u˜) inW₀^1,q(Ω).

Proof. Theorem 2.3 in [29] provides the strong convergencey_k→y˜inH₀¹(Ω). The uniform bound- edness of the weakly convergent sequence (∥uk∥L²(Ω))^k∈N, together with Proposition 2.6 and the boundedness of the productγ_kk, yields a uniform bound on∥ξ_k∥L²(Ω), and thus a subsequence con- verging weakly to the limitξ˜=Ay˜−u˜−f inL²(Ω). Owing to the compact embedding ofL²(Ω) intoW^−1,q(Ω), this convergence and also the convergence of(u_k)k∈Nis strong in W^−1,q(Ω). The continuity ofA⁻¹ ∶W^−1,q(Ω) →W₀^1,q(Ω)(see Remark 2.2) then implies strong convergence of the subsequence(yk)^k∈NinW₀^1,q(Ω).

2.3 Penalized optimal control problem with smoothed objective

We define the smoothed penalized optimal control problem with smoothed objective as follows:

MinimizeJ_r(y, u) = ∑

w∈I

1

2∣B_r(w)∣ ∥y−y_w∥²L²(Br(w))+ν

2∥u∥²L²(Ω) (2.17a)

over(y, u) ∈H₀¹(Ω) ×U_ad (2.17b)

s.t.Ay−γmax(0,−y) =u+f inH⁻¹(Ω), (2.17c) wherer > 0is sufficiently small and B_r(w) = {x∈ Ω∣ ∣x−w∣ < r}. For each r > 0, J_r is weakly lower semi-continuous. The smoothing of the objective in the sense of (2.17a) is used below in order to establish a suitable stationarity system for (2.6). In this section, we show that the averaged smooth-penalty scheme (2.17) is consistent with the optimal control problem (2.6) in the sense of Theorem 2.11 below.

Proposition 2.9. For allγ >0,>0andr>0, problem (2.17) has a solution.

Proof. The functionalJ_r is bounded from below. Therefore the set

{J_r(y, u) ∣ (y, u) ∈H₀¹(Ω) ×U_adsolving (2.17c)}

has an infimum denoted by j, and we can choose an infimizing sequence (y_k, u_k)k∈N with limit lim_k→∞J_r(y_k, u_k) =j. The sequence(u_k)^k∈Nis uniformly bounded inL²(Ω)and therefore contains a subsequence, which we also denote by(uk), with weak limitu∈Uad. Together with(∥uk∥L²(Ω))^k∈N, by Proposition 2.6, the sequence (∥y_k∥W₀^1,q(Ω))^k∈N is bounded. By the continuous embedding of

W₀^1,q(Ω)intoH₀¹(Ω), (∥yk∥H₀¹(Ω))^k∈Nis also bounded and thus contains a weakly convergent sub- sequence with limity ∈ H₀¹(Ω). The limiting pair(y,˜ u˜) is feasible for problem (2.17). The weak lower semi-continuity of the objectiveJ_rfinally implies

j=lim inf{J_r(y_k, u_k) ∣k∈N} ≥J_r(y,˜ u˜).

Lemma 2.10. Let the sequence (r_k)^k∈N ⊂R^>0 converge to zero and(G_k)^k∈N ⊂C⁰(Ω¯)converge to GinC⁰(Ω¯). Then for everyw∈Ω,

1

∣Br_k(w)∣ ∫^Brk(w)G_k(x)dx→G(w). Proof. For everyk∈N, we define the mapping

F_k∶C⁰(Ω¯) →R, g↦ 1

∣B_rk(w)∣ ∫^Brk(w)g(x)dx,

which is linear and bounded and thus an element of the dual spaceC⁰(Ω¯)^∗. Every pointw∈Ωis a Lebesgue point of the continuous functiong∈C⁰(Ω¯), i.e. fork→ ∞

⟨F_k, g⟩C⁰(Ω¯)^∗ = 1

∣B_rk(w)∣ ∫^Brk(w)g(x)dx→g(w) = ⟨δ_w, g⟩C⁰(Ω¯)^∗.

From this we concludeF_k⇀^∗δ_w inC⁰(Ω¯)^∗. Together with the strong convergence ofG_kinC⁰(Ω¯), this yields the convergence of the product

⟨F_k, G_k⟩C⁰(Ω¯)^∗ → ⟨δ_w, G⟩C⁰(Ω¯)^∗=G(w).

In order to formulate the central theorem on consistency in this section, we denote a solution of the smoothed penalized problem (2.17) with parameters(γ, , r) = (γ_k, _k, r_k)by(y_k, u_k).

Theorem 2.11. Assume that (γ_k)k∈N,(r_k)k∈N ⊂ R^>0 converge to zero and let _k = (γ_k) satisfy Assumption 2.7. Then there exists a subsequence of(y_k, u_k)^k∈N(denoted the same) and a solution (y,˜ u˜)of the original problem (2.7) such that

y_k→y˜ inW₀^1,q(Ω), u_k⇀u˜ inL²(Ω).

Proof. We choose u¯∈U_ad and denote the solution of the semi-linear equation (2.17c) withγ =γ_k and=_kbyy_γk(u¯)for allk∈N. Proposition 2.6, together with Assumption 2.7, yields a bound on (∥y_γk(u¯)∥W₀^1,q(Ω))^k∈Nthat does not depend onk∈N, and we can estimate

α

2∥u_k∥²L²(Ω)≤J_rk(y_k, u_k) ≤J_rk(y_γk(u¯),u¯) ≤C.

The uniform boundedness of(u_k)^k∈NinL²(Ω)thus yields a subsequence, still denoted by(u_k)^k∈N, with weak limit u˜ ∈ Uad. Theorem 2.8 then yields the strong convergence of a (sub-)subsequence

(yk)^k∈Nto the solutiony˜=y(u˜)of (2.2) inW₀^1,q(Ω). To prove optimality of the limiting pair(y,˜ u˜), let(y^⋆, u^⋆)denote a solution of problem (2.7). Using Lemma 2.10 for the first term inJ_rk(y_γk(u^⋆), u^⋆), asy_γk(u^⋆) →y(u^⋆)inC⁰(Ω¯)by Theorem 2.8 and the embedding ofW₀^1,q(Ω)intoC⁰(Ω¯), we ob- serve that

lim inf

k→∞ J_rk(y_γk(u^⋆), u^⋆) = lim

k→∞J_rk(y_γk(u^⋆), u^⋆) =J(y(u^⋆), u^⋆) =J(y^⋆, u^⋆).

Furthermore, with the weak lower semi-continuity of the second term inJ_rk(y_k, u_k), and the same arguments for the first term, we note that

lim inf

k→∞ Jr_k(yk, uk) ≥J(y,˜ u˜).

Now we exploit the optimality of(y_k, u_k)with respect to the objectiveJ_rk as well as the feasibility of(y_γk(u^⋆), u^⋆) for the smoothed penalized problem to deriveJ_rk(y_k, u_k) ≤J_rk(y_γk(u^⋆), u^⋆) for allk∈N. Altogether it holds that

J(y,˜ u˜) ≤lim inf

k→∞ J_rk(y_k, u_k) ≤lim inf

k→∞ J_rk(y_γk(u^⋆), u^⋆) =J(y^⋆, u^⋆).

3 Derivation of stationarity conditions

3.1 Optimality system for the smoothed penalized problem

For the derivation of a first order optimality system for (2.17) we apply [51, Thm. 3.1]. For this purpose we have to check the associated constraint qualification, i.e. regularity of a solution.

Lemma 3.1. Every feasible point(y_k, u_k) ∈H₀¹(Ω) ×L²(Ω)of (2.17) with parametersγ_k, _k∈R^>0 is regular in the sense of [51].

Proof. For(y, u) ∈ X ∶=H₀¹(Ω)×L²(Ω), defineg(y, u) =Ay−γ_kmax_k(0,−y)−u−fwith Fr´echet derivative

g^′(y_k, u_k)(y, u) =Ay+γ_kmax^′_k(0,−y_k)y−u. (3.1) We formulate (2.17c) asg(y, u) ∈ {0} ⊂H⁻¹(Ω) and, according to [51, Equ. (1.4)], have to show that

g^′(yk, uk)C(yk, uk) =H⁻¹(Ω),

whereC(yk, uk) ∶= {t((y, u)−(yk, uk)) ∣t≥0,(y, u) ∈H₀¹(Ω)×Uad}. Note that the definitions ofA andmax_k in (2.1) and (2.11) yield immediately thatg^′(y_k, u_k)C(y_k, u_k) ⊂ H⁻¹(Ω). Furthermore, the operatorg^′(y_k, u_k) ∶ C(y_k, u_k) →H⁻¹(Ω)is surjective, because for any given φ∈H⁻¹(Ω), we can choose u = u_k ∈ U_ad and y = y˜+y_k ∈ H₀¹(Ω), where y˜ is the solution of the linear elliptic partial differential equation Ay+γ_kmax^′

k(0,−y_k)y = φ. Then with t = 1, (y,˜ 0) ∈ C(y_k, u_k) and g^′(yk, uk)(y,˜ 0) =φ.

As a consequence of [51, Thm. 3.1] and Lemma 3.1 we obtain stationarity conditions for the smoothed penalized problem (2.17) as stated in the next proposition.

Proposition 3.2. If(y_k, u_k) ∈ H₀¹(Ω) ×U_ad withU_ad ≠L²(Ω)is optimal for problem (2.17) with parameters(γ_k, _k, r_k) ∈ (R^>0)³, then there existsp_k ∈ H₀¹(Ω)such that the following first order conditions hold:

u_k= 1

νp_k+ (a−1 νp_k)

+− (1

νp_k−b)

+ inL²(Ω), (3.2a) A^∗pk+γkmax^′_k(0,−yk)pk= − ∑

w∈I

1

∣B_rk(w)∣(yk−yw)χ_Brk(w) inH⁻¹(Ω). (3.2b) IfU_ad=L²(Ω), then (3.2a) is replaced byu_k=_ν¹p_k.

Remark 3.3. We introduce the variable σ_k = νu_k−p_k ∈ L²(Ω), which can be decomposed as σ_k =σ_ak −σ_bk withσ_ak =max(0, σ_k)andσ_bk = −min(0, σ_k)pointwise almost everywhere. Then the following conditions are equivalent to (3.2a):

(νu_k, u−u_k)L²(Ω)− (p_k, u−u_k)L²(Ω)≥0 ∀u∈U_ad, (3.3) uk=ΠU_ad(1

νpk), (3.4)

u_k∈U_ad, σ_ak ≥0, σ_bk ≥0, σ_ak(a−u_k) =σ_bk(b−u_k) =0. (3.5) HereΠ_Uad denotes theL²(Ω)-projection into the closed convex setU_ad.

3.2 Stationarity system for the optimal control problem

Borrowing terminology from [28], next we define the stationarity concept which is relevant in our context.

Definition 3.4(Limiting ε-almost C-stationarity). We call(y, u, ξ) ∈ W₀^1,q(Ω) ×U_ad×W^−1,q(Ω), U_ad ≠ L²(Ω) limiting ε-almost C-stationary for problem (2.6) if y = y(u) solves the variational inequality (2.2), ξ is defined as ξ = Ay−u−f, and if there exist p ∈ W₀^1,q′(Ω), λ ∈ M^b(Ω¯) and sequences (p_k)^k∈N ⊂ H₀¹(Ω) with p_k ⇀ p in W₀^1,q′(Ω), and (λ_k)^k∈N ⊂ H⁻¹(Ω) with λ_k ⇀^∗ λ in Mb(Ω¯)such that the following conditions are satisfied,

u−1

νp− (a−1 νp)

++ (1 νp−b)

+=0 inL²(Ω), (3.6a)

∀ψ∈W₀^1,q(Ω) ∶ ⟨ψ, A^∗p⟩W₀^1,q(Ω)− ⟨ψ, λ⟩W₀^1,q(Ω)+ ∑

w∈I

(y(w) −y_w)ψ(w) =0, (3.6b)

⟨λ, y⟩W^−1,q′(Ω)=0, (3.6c)

∀τ >0 ∃E_τ ⊂ {y>0}such that∣{y>0} ∖E_τ∣ <τ and

∀ϕ∈C⁰(Ω¯), ϕ∣^Ω∖Eτ =0, ⟨λ, ϕ⟩Mb(Ω¯)=0, (3.6d) lim sup{⟨λ_k, p_k⟩^H−1(Ω)∣k∈N} ≤0, (3.6e)

⟨ξ, p⟩^W−1,q(Ω)=0. (3.6f)

IfUad=L²(Ω), then (3.6a) is again replaced byu= ¹_νp∈W₀^1,q′(Ω).

Theorem 3.5. For eachk ∈N, let γk, k =(γk), rk >0be penalization and smoothing parameters which satisfy Assumption 2.7, whereγ_k→ ∞andr_k→0.

Furthermore let(y_k, u_k, p_k)be stationary for problem (2.17) in the sense that the tuple is feasible and (3.2a), (3.2b) hold, and assume that(∥u_k∥L²(Ω))^k∈Nis bounded.

Then there exist(y,˜ u,˜ ξ,˜p,˜ λ˜) ∈W₀^1,q(Ω)×L²(Ω)×L²(Ω)×H₀¹(Ω)×Mb(Ω¯)and a subsequence (also denoted by indexk) such that

yk→y˜ inW₀^1,q(Ω), (3.7a)

u_k⇀u˜ inL²(Ω), (3.7b)

ξ_k∶=γ_kmax_k(0,−y_k) ⇀ξ˜ inL²(Ω), (3.7c)

p_k⇀p˜ inH₀¹(Ω), (3.7d)

λ_k∶= −γ_kmax^′

k(0,−y_k)p_k⇀^∗ ˜λ inC⁰(Ω¯)^∗, (3.7e) and(y,˜ u,˜ ξ˜)is a limitingε-almost C-stationary point for problem (2.6) with multipliersp,˜λ.˜ Proof. Convergence and feasibility of yk, uk, ξk andpk. Theorem 2.8 yields the assertions for yk

andξ_k. Testing the adjoint equation withp_k∈H₀¹(Ω), we obtain a uniform bound on(∥p_k∥H₀¹(Ω))^k∈N

as follows,

C∥p_k∥²H₀¹(Ω)≤⟨p_k, A^∗p_k⟩H⁻¹(Ω)

= −γ_k∫_Ωmax^′

k(0,−y_k)p²_kdx− ∑^w∈I 1

∣B_rk(w)∣∫^Brk(w)(y_k−y_w)p_kdx

≤C(card(I) ∥y_k∥L^∞(Ω)+sup{∣y_w∣ ∣w∈I}) ∥p_k∥L²(Ω).

The weak convergence of a subsequence of(u_k)^k∈NinL²(Ω)follows from the boundedness of the operator (⋅)+ ∶ L²(Ω) → L²(Ω). Note that if a, b ∈ H¹(Ω), and a∣^∂Ω < 0 < b∣^∂Ω, then we have convergence along a subsequence inH₀¹(Ω).

Convergence ofλ_k. Forδ>0we define the function

ρ_δ∶R→R, ρ_δ(p) =⎧⎪⎪⎪

⎨⎪⎪⎪⎩

−1 forp< −δ,

p

δ forp∈ [−δ, δ], 1 forp>δ,

and note that lim_δ→0+ρ_δ(p) = sign(p) for allp ∈ R. Furthermore, ρ_δ(p_k) ∈ H₀¹(Ω) ∩L^∞(Ω) is a feasible test function for the adjoint equation (3.2b), yielding

⟨A^∗p_k, ρ_δ(p_k)⟩^H−1(Ω)+γ_k∫_Ωmax^′_k(0,−y_k)p_kρ_δ(p_k)dx

= − ∑

w∈I

1

∣Br_k(w)∣ ∫^Brk(w)(y_k−y_w)ρ_δ(p_k)dx≤card(I)∥y_k∥^L∞(Ω)+sup{∣y_w∣ ∣w∈I} ≤C.

The first term on the left hand side is non-negative as

⟨A^∗p_k, ρ_δ(p_k)⟩H⁻¹(Ω)≥C∫_Ωρ^′(p_k)∣∇p_k∣²dx≥0.

This yields the uniform boundedness of γ_k∫^Ωmax^′_k(0,−y_k)p_kρ_δ(p_k)dx ≥ 0. Sending δ → 0, we obtain the bound

∥λ_k∥^L1(Ω)= ∫_Ω∣λ_k∣dx=lim

δ→0∫_Ωγ_kmax^′

k(0,−y_k)p_kρ_δ(p_k)dx≤C. (3.8) The embedding ofL¹(Ω)into M^b(Ω¯) = C⁰(Ω¯)^∗ and the Banach-Alaoglu theorem then provides the existence of a subsequence converging weakly^∗ toλ˜∈ Mb(Ω¯). Note that because of the embed- ding ofW₀^1,q(Ω)intoC⁰(Ω¯), the subsequence also converges weakly inW^−1,q′(Ω)and(λ_k)^k∈N is bounded in this space.

Adjoint equation(3.6b) for y,˜ p,˜ λ.˜ Since p_k ⇀ p˜in W₀^1,q′(Ω) ⊃ H₀¹(Ω), so by the properties ofAwe haveA^∗p_k ⇀A^∗p˜inW^−1,q′(Ω). Ifϕ∈C⁰(Ω¯), then∑^w∈I(y_k−y_w)ϕ→ ∑^w∈I(y˜−y_w)ϕin C⁰(Ω¯)and by Lemma 2.10,

⟨∑w∈I

1

∣B_rk(w)∣(y_k−y_w)χ_Brk(w), ϕ⟩C⁰(Ω¯)^∗ = ∑

w∈I

1

∣B_rk(w)∣ ∫^Brk(w)(y_k−y_w)ϕ dx

→ ∑

w∈I

(y˜(w) −y_w)ϕ(w). The embedding ofC⁰(Ω¯)^∗intoW^−1,q′(Ω)then gives

∑

w∈I

1

∣B_rk(w)∣(y_k−y_w)χ_Brk(w)⇀ ∑

w∈I

(y˜−y_w)δ_w inW^−1,q′(Ω). Together with the weak convergence ofλkinW^−1,q′(Ω), this yields equation (3.6b).

Complementarity ofy˜ andλ˜ (3.6c). From the convergence y_k → y˜in W₀^1,q(Ω) we infer that (−y_k)+→ (−y˜)+=0inW₀^1,q(Ω), so

⟨λ_k,(−y_k)+⟩W^−1,q′(Ω)→ ⟨λ,˜ 0⟩W^−1,q′(Ω)=0.

Furthermore we observe that

∣⟨λk,(yk)+⟩W^−1,q′(Ω)∣ =γk∫_Ωmax^′_k(0,−yk)∣pk∣(yk)+dx

=γ_k∫_{0≤y

k≤_k}max^′

k(0,−y_k)∣p_k∣y_kdx+γ_k∫{yk>k}0dx≤γ_kk∥p_k∥L¹(Ω),

which converges to zero due to the boundedness of (pk)^k∈N inW₀^1,q′(Ω)and Assumption 2.7. So decomposingy_k= (y_k)+− (−y_k)+we get

⟨λ_k, y_k⟩W^−1,q′(Ω)= ⟨λ_k,(y_k)+⟩W^−1,q′(Ω)− ⟨λ_k,(−y_k)+⟩W^−1,q′(Ω)→0.

We also have the convergence⟨λk, yk⟩W^−1,q′(Ω)→ ⟨˜λ,y˜⟩W^−1,q′(Ω), which yields the assertion (3.6c).