DIRECTIONAL DIFFERENTIABILITY FOR SHAPE OPTIMIZATION WITH VARIATIONAL INEQUALITIES AS CONSTRAINTS

. For equilibrium constrained optimization problems subject to nonlinear state equations, the property of directional diﬀerentiability with respect to a parameter is studied. An abstract class of parameter dependent shape optimization problems is investigated with penalty constraints linked to variational inequalities. Based on the Lagrange multiplier approach, on smooth penalties due to Lavrentiev regularization, and on adjoint operators, a shape derivative is obtained. The explicit formula provides a descent direction for the gradient algorithm identifying the shape of the breaking-line from a boundary measurement. A numerical example is presented for a nonlinear Poisson problem modeling Barenblatt’s surface energies and non-penetrating cracks

The classical theory of state-constrained optimization problems deals with linear equations, typically, given by partial differential equations [29,38]. In our consideration we study state constraints, given by variational inequalities and their penalization. The challenge consists in the fact that the latter are not Fréchet differentiable (see [30,35]). As a consequence, the directional derivative of Lagrangians and related shape differentiability fails. For the concept of the conical differential of a solution of the Signorini variational inequality, see [36]. Sensitivity estimates in shape optimization problems for a class of semi-linear elliptic variational inequalities based on material derivatives were investigated in [9]. Shape sensitivity analysis for an inverse obstacle problem and its regularization via penalization was performed in [13] utilizing geometric properties of active and biactive sets. We suggest a novel approximation for the shape derivative along specifically linearized directions and based on generalized adjoint techniques (see [31,32]). In particular cases, our approach is closely related to the method of averaged adjoints developed by [26]. We compute a directional derivative with respect to specific shape perturbations, which give descent directions for the shape optimization problem. It coincides with the usual shape derivative in case of linear state equations.
Our research addresses the following features: (i) Optimization subject to nonlinear and nonsmooth equilibrium constraints. Within the Lagrange multiplier approach (see [14]), in Section 2 we consider a convex objective function J with a nonlinear equation as constraint. The linearized Lagrangian L is well-posed when using the associated adjoint operator provided in Lemma 2.2.
(ii) Penalty optimization linked to variational inequalities (VI). In order to treat VI, in Section 3 we extend L to a penalized Lagrangian L ε for ε > 0 (see Lem. 3.2), thus reducing the variational inequalities to case (i). Using adjoints we derive necessary and sufficient optimality conditions for a saddle-point problem providing the optimal value l(ε).
(iv) Limit as ε → 0 + . Taking the limit as the penalty parameter ε → 0 + , in Theorem 3.4 we derive the reference variational inequality, its adjoint equation, as well as primal and adjoint variables that are the Lagrange multipliers for inequality constraints. However, a limit directional derivative fails since the VIs are not differentiable.
(v) Shape derivative. In Section 4 we introduce states depending on a family of geometries Ω t parameterized by t. The diffeomorphic perturbations Ω t+s are defined by a kinematic velocity Λ (see e.g. [10,23]). They are used to characterize the shape derivative using the bijection property of the function spaces V (Ω t+s ) → V (Ω t ).
(vi) Application to non-penetrating Barenblatt's cracks. We apply shape perturbations to the nonlinear crack problem (4.11) under non-penetration [16,17] in the anti-plane setting (see [11]). Beyond the classic Griffith's brittle fracture, Barenblatt's cohesion (see [2,18]) allows crack faces to close smoothly and determines a-priori unknown cracks by those points where opening occurs along a breaking line Σ t .
(vii) Hadamard formula and descent directions. In Theorem 4.3 we specify the shape derivative for the nonlinear Poisson problem described in (vi), and express it by a Hadamard formula over the moving boundary in Theorem 4.4. This formula provides kinematic velocities Λ for a descent direction ∂ + j(ε, 0) < 0 within a gradient method.
(viii) Identification of breaking lines. Finally, in Section 5 we present a numerical simulation of the gradient descent algorithm for the inverse problem of identification of the breaking line Σ t , which minimizes the objective J of least-square misfit from a boundary observation. We report that the faces need to be open for identification within DPA.

Directional differentiability of Lagrangians for equilibrium constraints
In separable Banach spaces V and X, let a linear operator M : V → X map the space of states u ∈ V to observations z ∈ X. We consider an abstract objective function dependent on a positive parameter s ∈ I := [0, s 0 ), s 0 > 0: (2.1) Next we introduce our state constraint. Let the continuous function E(s, u) : I × V → R be the energy functional. For every fixed s we assume that its is differentiable, i.e.
Here and in what follows the brackets · , · stands for the duality pairing between V and its dual space V .
We define the reference state as a solution u 0 ∈ V at s = 0 to the equilibrium equation expressed in the variational form: The variational equation (2.2) constitutes the optimality condition for the minimum Lemma 2.1. Let assumption (E1) and the following hold: (E2) E at s = 0 is coercive: there exist a > 0 and f ∈ V such that Then there exists a solution u 0 ∈ V to (2.2).
Proof. We introduce a Galerkin approximation of (2.2) by nonlinear equations in subspaces V n ⊂ V of finite dimension n ∈ N as follows Since the strong and weak convergences coincide in finite-dimensional spaces, under the coercivity and continuity assumptions (E2) and (E3) solutions u n ∈ V n exist according to the Brouwer fixed point theorem, see e.g. [5]. The solutions are uniformly bounded in V due to (E2). Hence there exists a weakly convergent subsequence u n k and an accumulation point u 0 ∈ V . Taking the limit as n k → ∞ due to (E3) the assertion of the theorem follows.
Induced by the state equation (2.2), we have the optimal value Our aim is to extend the state-constrained optimization (2.4) to a well-posed optimal value function j : I ⊂ R → R in such a way that it has a directional derivative at s = 0: Further we linearize the mapping u → E around the reference solution u 0 to (2.2) and use a Lagrange multiplier method [14]. For this task we employ for fixed (s, u 0 ) ∈ I × V an 'associated to adjoint' operator (E ) (s, u 0 ) ∈ L (V, V ), which is defined by means of the Lagrange identity (see [32], Chap. 1): (2.6) Lemma 2.2. If the following assumption holds: for s ∈ I, and r → E (s, ru 0 ) is continuous for r ∈ [0, 1], then an associated to adjoint operator in (2.6) is given by Proof. From the Newton-Leibniz axiom we have Inserting w = u 0 into (2.7) and using (2.8) implies (2.6).
Based on Lemma 2.2, a linearized Lagrange function L : For the Lagrangian L we consider the saddle-point (minimax) problem: for all (u, v) ∈ V 2 . Following [3], we introduce the optimal values: which determine a multi-valued function [s ⇒ K s × K s ] : I ⇒ V 2 . Later we shall prove that these sets are not empty.
Lemma 2.3. Let (E1)-(E3), (E 1) and the following assumptions hold where · , · X ,X is the duality pairing between X and its dual space X ; (J2) the objective functional is convex: (E 2) the associated to adjoint operator is symmetric: Then for every s ∈ I there exists a state u s ∈ V solving the equation:

12)
and an adjoint state v s ∈ V satisfying the adjoint equation: The pair (u s , v s ) ∈ K s × K s is a saddle point satisfying l(s) := l s = L(s, u 0 , u s , v s ) = l s , s ∈ I. (2.14) If f = 0 in (E 3), then the saddle-point is unique.
The proof of Lemma 2.3 is given in Appendix A. We note that u 0 from Lemma (2.1) is a solution to the s-dependent equation (2.12) at s = 0. The latter also coincides with the reference equation (2.2) due to the Lagrange identity (2.6).
The next lemma establishes a sequential semi-continuity property for the solution set K s × K s as s → 0 + . (J4) s → J (s, M u s ) on solutions u s ∈ K s is continuous as s → 0 + ; Then there exist s k → 0 + , a subsequence of saddle points The proof of Lemma 2.4 is technical and is presented in Appendix B. In the last lemma of this section directional differentiability of Lagrangians [3,4] is recalled.
Lemma 2.5. Let the set of saddle points (u s , v s ) ∈ K s × K s satisfying (2.14) be nonempty for each s ∈ I; assume that a subsequence (u s k , v s k ) ∈ K s k × K s k and an accumulation point (u 0 , v 0 ) ∈ K 0 × K 0 exist satisfying strong convergence (2.15) as s k → 0 + . If the following holds: (L1) There exists a partial derivative ∂L/∂s : I × V 3 → R of the Lagrangian L with respect to the first argument at r ∈ I such that lim inf then for the optimal value function j : I → R defined as j(s) := J (s, M u s ) for u s ∈ V solving (2.12), (2.16) the directional derivative ∂ + j(0) in (2.5) exists. It is equal to a directional derivative ∂ + l(0) for the optimal value function l : I → R from (2.14) and is expressed by the partial derivative ∂L/∂s as follows The proof of Lemma 2.5 follows [4], Chapter 10, Theorem 5.1. Based on Lemmas 2.1-2.5 we state the main theorem of this section. 3 it follows that the set of saddle points (u s , v s ) ∈ K s × K s satisfying (2.14) is nonempty. Together with assumptions (E4), (E5), (J3), (J4), (E 4), (E 5) Lemma 2.4 guarantees the existence of a subsequence (u s k , v s k ) ∈ K s k × K s k and an accumulation point (u 0 , v 0 ) ∈ K 0 × K 0 satisfying the strong convergence (2.15) as s k → 0 + . Utilizing (L1) Lemma 2.5 implies the assertion of the theorem.
In the following section we extend the directional differentiability result of Theorem 2.6 to a penaltyconstrained optimization motivated by variational inequalities.

Directional differentiability of Lagrangians due to penalty constraints
Let H be another Banach space with an order relation denoted by '≥'. We introduce a parameter-dependent family of linear operators B(s) ∈ L (V, H) with s ∈ I, and the associated inequality constraints As a canonical example we may consider a trace operator. Using the decomposition ζ = Compared to (2.3), the constrained problem at s = 0: leads to the variational inequality: In order to bring (3.4) in equality form akin (2.2), we regularize it by a penalty approximation. For a small penalization parameter ε ∈ (0, ε 0 ), ε 0 > 0, we define the penalty as a map β ε (s, ζ) : I × H → H into the dual space H . For the constraint [ζ] − = 0 according to (3.2), the standard penalty function β ε (0, ζ) = −[ζ] − /ε ≤ 0 forces the compliance condition β ε (0, ζ)ζ = ([ζ] − ) 2 /ε. However, the min-based penalty function is not differentiable (see As. (L3)). Therefore, we introduce a Lavrentiev relaxation [27] satisfying: where · , · H ,H denotes the duality pairing between H and H .
For example, a smooth ε-mollification of the minimum function is depicted in Figure 1 together with its derivative. It satisfies (B1) with β = −β ε (0, 0) = exp(−2) and This leads to the penalized problem: find u ε 0 ∈ V such that then there exists a solution u ε 0 ∈ V to (3.7). Proof. The operator of problem (3.7) is coercive due to assumption (E2) and the lower bound in (3.5). It is weakly continuous due to (E3) and (B2). The proof of Lemma 2.1 can be adapted to guarantee existence of a solution.
Following Lemma 2.2 we assume that for ζ, η ∈ H, and the mapping r → β ε (s, rB(0)u ε 0 ) is continuous for r ∈ [0, 1], where s ∈ I. Then the adjoint β ε (s, B(0)u ε 0 ) ∈ L (H, H ) exists, it is given by and satisfies the Lagrange identity for ζ ∈ H, s ∈ I: Using (3.7) and (3.9) we modify (2.9) with a penalized Lagrange function L ε : The penalized saddle-point problem reads: The optimal values and solution sets in (2.11) are We establish results for (3.12) analogous to those of Lemmas 2.3 and 2.4.
replacing u 0 , and the following assumptions hold true: is continuous for s ∈ I. Then for every s ∈ I there exist a state u ε s ∈ V solving the equation: and an adjoint state v ε s ∈ V satisfying the adjoint equation: 14) The The proof of Lemma 3.2 is technical and presented in Appendix C.
Below we state a theorem on differentiability of L ε .
, and the two following assumptions hold: The directional derivative of the optimal value function j : (3.17) and the associated Lagrangian function l : Here the partial derivative is given by solves the penalty problem (3.7) and the adjoint equation (3.14) at s = 0: Proof. The differentiability assumptions (L1)-(L3) together with the continuity in (B4), (B 5) imply the existence of the partial derivative of L ε in (3.19) with respect to s ∈ I and its semi-continuity properties: Therefore, utilizing Lemma 3.2 and proceeding as in Lemma 2.5 we obtain formula (3.18) for the directional derivative. Taking the limit s → 0 + in (3.13) and using (3.9) we arrived at (3.7). The adjoint equation (3.20) follows from (3.14). The proof is complete.
Next we analyze the limit as ε → 0 + . For this task we employ the Lagrangian L from (2.9) at s = 0.
) and the following assumptions hold: Then there exists a quadruple (u 0 , λ 0 , v 0 , µ 0 ) ∈ (V ×H ) 2 , whereH is the dual space toH from (B8) with the duality pairing · , · H ,H , which satisfies the primal problem: the adjoint problem: the complementarity relations: and the compatibility condition Moreover, u 0 satisfies [B(0)u 0 ] − = 0 and the variational inequality (3.4). Together with the Lagrange multiplier λ 0 it solves (3.25) The adjoint v 0 solves the variational equation for all u ∈ V : for µ 0 obtained as an accumulation point in the following limit:  Figure 2. An example geometry Ω t in 2D.
The proof of Theorem 3.4 is technical and it is presented in Appendix D.
It is worth noting that we cannot pass to the limit as ε → 0 + in the derivative β ε of the penalty, since it is unbounded in general, see Figure 1. This would be needed for β ε which enters into ∂L ε /∂s in (3.18).

Shape derivative for breaking-line identification
Now we turn to a model problem for a nonlinear Poisson equation. We derive a shape derivative suitable for shape optimization in the problem of breaking-line identification from a boundary measurement. Let be a parameter dependent family of domains contained in the hold-all domain D. For some fixed t ∈ (t 0 , t 1 ) we refer to Ω t as the reference domain. We assume that Ω t = Ω + t ∪ Ω − t ∪ Σ t is split into two variable subdomains Ω ± t with Lipschitz boundaries ∂Ω ± t and outward normal vectors n ± t . The sub-domains are separated by a one-dimensional breaking line with the chosen normal direction ν t = n − t = −n + t (see Fig. 2). Let the outer boundary be split into two variable parts without intersection ∂Ω t = Γ D t ∪ Γ N t , and the outward normal vector n t be such that n ± t = n t at ∂Ω t . The condition Γ D t ∩ ∂Ω ± t = ∅ on the Dirichlet boundary is assumed to guarantee the Poincaré inequality in Ω ± t . A part of the Neumann boundary Γ O t ⊂ Γ N t builds the observation boundary. Further we introduce We adopt the formalism from Sections 2 and 3 to the geometry-dependent spaces of functions The observation operator M : V (Ω t ) → L 2 (Γ O t ) maps to the boundary traces on Γ O t . The restriction operator B : V (Ω t ) → L 2 (Σ t ) is independent of s and describes a jump across the breaking line Σ t subject to the non-penetration condition (see motivation in [11]): Here we take into account the dissipative interaction phenomenon of cohesion (see [2,18]) described by a surface energy density α(s, ζ). The following conditions are imposed: 6) and the existence of K α1 > 0, K α2 > 0 such that: For example, a mollification of the function (K c /κ) min(κ, |ζ|) as where 0 < δ < κ, κ > 0, and K c > 0 is the fracture toughness parameter. The function from (4.8) is depicted in Figure 3. Let the Lame parameter µ L > 0 and the traction force g ∈ H 1 (D N ), ensuring that g ∈ L 2 (Γ N t ) on Lipschitz curves Γ N t ⊂ D N , be given. The bulk and the surface energies together constitute the total potential energy E(0) : V (Ω t ) → R: We calculate the Gateaux derivative E (0) : V (Ω t ) → V (Ω t ) at u: where denotes the transpose. The constrained optimization (3.3) leads to the variational inequality (3.4), which takes the form: There exists a solution to the variational inequality (4.11). It satisfies the linear complementarity problem: The solution is unique for convex α (hence, monotone α ).
Proof. For u ∈ V (Ω t ) we recall the Poincaré inequality: (4.13) and the trace inequality: both uniform in t ∈ (t 0 , t 1 ). Using the bound K α1 > 0 in (4.7) and (4.13), (4.14) we can estimate E (0, u; Ω t ), u in (4.10) from below and conclude the coercivity property (E2). The weak-to-weak continuity (E3) for E (0, u) holds due to the continuity of α assumed in (4.6). Therefore, by Lemma 3.1 there exists a solution u ε t ∈ V (Ω t ) to the penalty equation (see (3.7)) in the form: for all v ∈ V (Ω t ). It satisfies the mixed boundary value problem: By the compactness argument used in the proof of Theorem 3.4 we get an accumulation point such that u ε k t u t weakly in V (Ω t ) as ε k → 0, which solves the variational inequality (4.11). The derivation of relations (4.12) is standard, see e.g [16], Chapter 1.
. We aim at the shape optimization problem for identification of an unknown breaking line from the observation: find Σ * as the solution to min Σt⊂DΣ j(0, 0) = J(0, u t ; Ω t ) : where J represents J from (2.4), and ρ ≥ 0 stands for the reason of perimeter regularization.
If ρ = 0, then there exists a solution to the shape optimization problem (4.17). In general, the solution is non-unique.
The diffeomorphism (4.21) preserves the bijectivity between the function spaces in (4.4): With the help of (4.26) we transform the perturbed objective J(0,ũ; Ω t+s ) from (4.19) forũ ∈ V (Ω t+s ) such that where ω s will be defined later. Based on the second derivative in the identity (see (2.6)): we linearize at the solution u ε t the perturbed state operator in (4.10): where the terms are In (4.27) and (4.30) we use the chain rule (4.31) and the Jacobian in the domain and at the boundary: for more details, see e.g. [19,20,25]. Similarly, using the following identity analogous to (4.28): we perturb the penalty term in (4.15) linearized at u ε t such that Next we present a formula for the shape derivative.
Theorem 4.3. Let the bound K α2 > 0 in (4.7) be sufficiently small such that a := µ L K P − 2K α2 K 2 tr > 0, (4.36) where K P and K tr are the constants from the Poincaré and the trace estimates (4.13) and (4.14). Then the directional derivative of L ε exists and is given by the formula where the tangential divergence is defined as The saddle point (u ε t , v ε t ) ∈ V (Ω t ) 2 solves the penalty equation (4.15) and the adjoint equation: for which the mixed boundary value formulation is given by: For the proof of Theorem 4.3 one checks the conditions of Theorem 3.3. It is given in Appendix E.
In the following we decompose the velocity into the normal and tangential vectors at the boundary: where τ ± t is the tangential vector positively oriented to n ± t . Theorem 4.4. Let the solution of (4.15), (4.39) be smooth such that (u ε t , v ε t ) ∈ H 2 (Ω ± t ) 2 . Then the shape derivative in Theorem 4.3 satisfies an equivalent Hadamard's representation by the boundary integrals: The terms in (4.42) are given by where κ ± t := div τt n ± t denotes the curvature at ∂Ω ± t , and we utilize the notation at Σ t :
The proof of Theorem 4.4 is based on integration by parts and is presented in Appendix F. The expression (4.42) is important for gradient-based iterative techniques.

Numerical simulation
We set a piecewise-linear breaking line Σ * ⊂ D Σ to be identified: which breaks the rectangle Ω = (0, 1) × (0, 0.5) into two parts Ω ± * . Let the boundary ∂Ω be split into fixed Dirichlet and Neumann parts: see the illustration of the geometry in Figure 2. We choose for the Young's modulus E Y = 73000 (mPa) and Poisson's ratio ν P = 0.34, the Lamé parameter µ L = E Y /(2(1 + ν P )) ≈ 27239, and the linear traction force Then there exists a solution z ∈ H 1 (Ω ± * ) such that z = 0 on Γ D * , [[z]] − = 0 on Σ * , which satisfies the variational equation (4.18) according to Lemma 4.1. Let the observation boundary be Γ O * = Γ N * . Now we discretize the problem. For Σ t ⊂ D Σ breaking Ω into Ω ± t , let Ω ± t,h be a triangulation of mesh size h > 0 of Ω ± t , which is compatible at the interface Σ t,h := Σ t ∩ ∂Ω 1 t,h = Σ t ∩ ∂Ω 2 t,h . At Σ t,h the cohesion function α(0, ζ) is set as in (4.8) with K c = 10 −3 (mPa·m), κ = 10 −2 (m). For small δ and h we rely on the discretization α h (0, ζ) such that After piecewise-linear FE discretization of the problem on a grid of mesh size h = 10 −2 according to (5.1)-(5.4), we solve the variational equation (4.18) by a primal-dual active set (PDAS) iterative algorithm developed in [12]. The numerical solution z h obtained after 3 iterations with zero residual is plotted in Figure 4   According to the proof given in Lemma 4.1 we approximate the variational inequality (4.11) by the penalty equation (4.15). For small ε and h, the penalty operator from (3.6) is discretized as The discrete penalty equation (4.15) determines u ε t,h ∈ V t,h (Ω t,h ) such that and ignoring the singularity of α h the discrete adjoint equation (4.39) reads: find v ε t,h ∈ V t,h (Ω t,h ) such that After solving problems (5.6) and (5.7), according to Theorem 4.4 we calculate D ε 3 at the moving boundary Σ t,h , and D 1 at Σ t,h ∩ Γ D * , where ρ = 1/µ L is set. By the virtue of (5.4), (5.5) here q ε,h = 0 and Since Γ D * and Γ N * = Γ O * are fixed in the identification problem, the normal velocity n t Λ = 0 at ∂Ω when k 2 = k 6 = k 7 = 0 in (4.45). The tangential velocity is set τ t Λ = 0 at Σ t by means of k 4 = k 5 = k 8 = k 9 = 0. Therefore, we get a descent direction when Λ 1,H = 0 and The scaling k 3 = 0.1h/ Λ 2,H C(Σ t,h ) is chosen, and the weight k 1 = k 3 / √ h at Γ D * was found empirically in [6]. We point out that the discrete velocity Λ H at the interface Σ t is defined on a coarser grid of size H > 0, compared to the mesh size h of the problem.
We summarize the optimization algorithm for breaking line identification.  Figure 5. The penalty parameter ε = 10 −10 was taken. In plot (a) the selected iterations n = 0, 10, 20, 40, 100, 200 of Σ (n) according to (5.10) are drawn in Ω in comparison with the true interface Σ * (the thick solid line). In plot (b) of Figure 5 we plot the ratio J (n) /J (0) of the objective optimal values recalled here to be 11) and the shape ratio ψ (n) − ψ * C([0,1]) / ψ (0) − ψ * C([0,1]) . The computed misfit ratios attain as minimum 12% and 88%, respectively. From the simulation we conclude the following feature. In Figure 5 (a) it can be observed that the left part of curve Σ * , where no contact occurs (see Fig. 4 (a)), is recovered well by the identification Algorithm 5.1,   whereas the right part of interface, where contact occurs, the initialization Σ (0) is almost unchanged during the iterations.
To remedy the hidden part of Σ * , we apply to the same configuration a traction force which is more stretching than that in (5.3): As the result, the whole Σ * is open without contact, however, the cohesion occurs at the interface as shown in Figure 6. In this case, the result of Algorithm 5.1 for n ∈ {0, . . . , 400} is depicted in Figure 7. The objective ratio attains the minimum 0, 4%, and the shape error ratio 25%. We observe in Figure 7 (a) that the whole curve Σ * is recovered well compared to the previous case of contacting faces.
On the basis of our numerical simulation, we conclude that the breaking line identification algorithm is consistent with the setup of destructive physical analysis (DPA).

Appendix A. Proof of Lemma 2.3
Let us define the quadratic functional E : It is weakly lower semi-continuous and coercive due to (E 3), Gateaux-differentiable by (E 2), and (E ) (s, u 0 ) = (E ) (s, u 0 ). Adding to E in (A.1) the linear term E (s, 0), v , the above properties provide an argument u s ∈ V of the minimum: with an optimality condition in the form of the variational equation ( After cancelling J (s, M u s ) and testing with v = v s ± w we obtain the variational equation (2.12). Conversely, (2.12) satisfies the first inequality of (2.10) as equality.
On the other side, the second inequality of (2.10) after cancelling the term − E (s, 0), v s reads Substituting here u = u s ± rw, dividing the results with r and passing r → 0, by the virtue of differentiability of J assumed in (J1), this leads to the variational equation (2.13). Conversely, by the convexity assumption (J2) the necessary optimality condition (2.13) is sufficient for the minimum in the second inequality of (2.10) provided by u s . This proves that (u s , v s ) ∈ V 2 is a saddle point to problem (2.10). The definition (2.11) of solution sets K s , K s implies that (u s , v s ) ∈ K s × K s and satisfies the equality (2.14). This completes the proof of Lemma 2.3.

Appendix B. Proof of Lemma 2.4
We test the primal equation (2.12) with v = u s , apply the Cauchy-Schwarz inequality, the coercivity (E 3) with u = u s , and the boundedness assumption (E4) to derive the upper bound Testing the adjoint equation (2.13) with u = v s , from (E 3) with u = v s and (J3) it follows similarly that We combine (B.1) and (B.2) together in the uniform in s ∈ I estimate Then there exist s k → 0 + , a subsequence of saddle points (u s k , v s k ) ∈ K s k × K s k and an accumulation point where (2.12) was tested with v = u s k − u 0 , and (2.6) and property (E 2) were used. Inserting u = v s k − v 0 into (E 3) and using (2.13) Taking the limit as k → ∞ in (B.5) and (B.6), we get (2.15) with the help of the weak convergence in (B.4) and the boundedness properties (E4), (J3), (E 4) of E , J , (E ) . Finally, taking the limits in the primal (2.12) and adjoint (2.13) equations and using the strong convergence (2.15) and the continuity assumptions (E5), (J4) and (E 5), this guarantees that the pair (u 0 , v 0 ) solves (2.2) (due to identity (2.6)) and (2.13) at s = 0. Therefore, (u 0 , v 0 ) ∈ K 0 × K 0 which ends the proof of Lemma 2.4.

Appendix C. Proof of Lemma 3.2
The modified quadratic functional E ε : is weakly lower semi-continuous and coercive due to (E 3), (B3), and (B 3). Using (E 2) and (B 2) its Gateaux derivative is given by Consequently, the variational equation (3.13) is an optimality condition for the minimizer u ε s ∈ V of the following problem: and the adjoint equation (3.14) provides an argument v ε s ∈ V for The uniqueness assertion is similar to Lemma 2.3 and done by coercivity.
The left-hand side of the saddle-point formulation (3.11) is equivalent to the primal problem: which implies equation (3.13). The right-hand side is equivalent to the adjoint equation (3.14) due to the convexity (J2). Then (u ε s , v ε s ) ∈ K s ε × K ε s satisfies the saddle-point condition (3.15).
Appendix D. Proof of Theorem 3.4 Passing s k → 0 + due to the strong convergence (3.16) we refine the estimates (C.4) as follows. Using the lower bound in (3.5) and (E2), from (3.7) tested with v = u ε 0 we get: which is uniform in ε ∈ (0, ε 0 ). From (C.5) as s k → 0 + it follows that From (D.7) we conclude the existence of a solution v 0 ∈ V to the limit adjoint equation (3.26) and the -weak convergence (3.27). Applying the convergences (D.5) and (3.27) to the identity (3.9) at s = 0, in the limit the compatibility condition (3.24) follows. Equation (3.26) is the necessary and sufficient optimality condition for the adjoint limit problem (3.22). The proof of Theorem 3.4 is complete.