Efficient Techniques for Shape Optimization with Variational Inequalities using Adjoints

TL;DR

This paper develops efficient shape optimization techniques for obstacle problems constrained by variational inequalities, establishing adjoint existence, deriving shape derivatives, and demonstrating convergence with numerical validation.

Contribution

It introduces a novel approach to handle shape optimization with variational inequalities using adjoints, including existence proofs, derivative formulas, and an effective algorithm.

Findings

Existence of adjoints for regularized problems

Convergence of solutions to unregularized limits

Numerical validation of the optimization algorithm

Abstract

In general, standard necessary optimality conditions cannot be formulated in a straightforward manner for semi-smooth shape optimization problems. In this paper, we consider shape optimization problems constrained by variational inequalities of the first kind, so-called obstacle-type problems. Under appropriate assumptions, we prove existence of adjoints for regularized problems and convergence to limiting objects of the unregularized problem. Moreover, we derive existence and closed form of shape derivatives for the regularized problem and prove convergence to a limit object. Based on this analysis, an efficient optimization algorithm is devised and tested numerically.

Equations

$$J(y,\Omega ):=\mathcal{J}(y,\Omega )+{\mathcal{J}}_{\text{reg}}(\Omega ):=\frac{1}{2}{\int }_{\Omega }{\left|y-\bar{y}\right|}^{2}dx+\nu {\int }_{{\Gamma }_{\text{int}}}1ds$$

$$\underset{{\Gamma }_{\text{int}}\in \mathcal{X}}{\min }J(y,\Omega )$$

$$a(y,v-y)\ge ⟨f,v-y⟩\forall v\in K:=\{\theta \in H_0^1(\Omega ):\theta (x)\le \phi (x)\text{ in }\Omega \},$$

$$\begin{array}{rl}a:H_0^1(\Omega )\times H_0^1(\Omega )& \to \mathbb{R}\\ (y,v)& \mapsto {\int }_{\Omega }\sum _{i,j}{a}_{i,j}{\partial }_{i}y{\partial }_{j}v+\sum _i{d}_i({\partial }_{i}yv+y{\partial }_{i}v)+byvdx\end{array}$$

$$a(y,v)+(\lambda ,v{)}_{L^2(\Omega )}$$

$$\begin{array}{rl}\lambda & \ge 0\text{in }\Omega \\ y& \le \phi \text{in }\Omega \\ \lambda (y-\phi )& =0\text{in }\Omega \end{array}$$

$$\lambda=\max\big{(}0,\lambda+c(y-\varphi)\big{)}\quad\text{for any }c>0.$$

$$\lambda=\alpha\cdot\max\big{(}0,\lambda+\varrho(y-\varphi)\big{)}\quad\text{for any }\varrho>0$$

$$\lambda=\max\big{(}0,\bar{\lambda}+c(y-\varphi)\big{)}\quad c\in(0,\infty)$$

$$a(y_{c},v)+(\max\big{(}0,\bar{\lambda}+c(y_{c}-\varphi)\big{)},v)_{L^{2}(\Omega)}=(f,v)_{L^{2}(\Omega)}\qquad\forall v\in H^{1}_{0}(\Omega),$$

$$a(y_{\gamma,c},v)+\left(\text{max}_{\gamma}\big{(}\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)},v\right)_{L^{2}(\Omega)}=(f,v)_{L^{2}(\Omega)}\quad\forall v\in H_{0}^{1}(\Omega).$$

$$\begin{array}{rl}& a(p_{\gamma ,c},v)+c\cdot {\left({\text{sign}}_{\gamma }(\bar{\lambda }+c(y_{\gamma ,c}-\phi ))\cdot p_{\gamma ,c},v\right)}_{L^2(\Omega )}\\ & =-(y_{\gamma ,c}-\bar{y},v{)}_{L^2(\Omega )}\forall v\in H_0^1(\Omega )\end{array}$$

$$y_{\gamma ,c}$$

$$y_c$$

$$\Phi_{\gamma}\colon H^{1}(\Omega)\rightarrow L^{2}(\Omega),y\mapsto\text{max}_{\gamma}\big{(}\bar{\lambda}+c\cdot(y-\varphi)\big{)}$$

$${\Psi }_c:H^1(\Omega )\to H^1(\Omega ),y\mapsto \bar{\lambda }+c\cdot (y-\phi )$$

$$\Phi\colon H^{1}(\Omega)\rightarrow L^{2}(\Omega),y\mapsto\text{max}\big{(}0,\bar{\lambda}+c\cdot(y-\varphi)\big{)}.$$

$$0$$

$$\le a(y_{\gamma ,c}-y_c,y_{\gamma ,c}-y_c)$$

$$\displaystyle\hskip 11.38092pt+\Big{(}\text{max}\big{(}0,\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)}-\text{max}\big{(}0,\bar{\lambda}+c(y_{c}-\varphi)\big{)},y_{\gamma,c}-y_{c}\Big{)}_{L^{2}(\Omega)}$$

$$\displaystyle=\Big{(}\text{max}_{\gamma}\big{(}\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)}-\text{max}\big{(}0,\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)},y_{\gamma,c}-y_{c}\Big{)}_{L^{2}(\Omega)}$$

$$\displaystyle\leq\int_{\Omega}|\text{max}_{\gamma}\big{(}\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)}-\text{max}\big{(}0,\bar{\lambda}+c(y_{\gamma,c}-\varphi)\big{)}|\cdot|y_{\gamma,c}-y_{c}|\;dx$$

$$\le g(\gamma )\cdot \text{vol}(\Omega {)}^{\frac{1}{2}}\cdot ||y_{\gamma ,c}-y_c|{|}_{H^1(\Omega )},$$

$$A_c\subseteq A\forall c>0,$$

$$\Vert {\text{sign}}_{\gamma }(\bar{\lambda }+c\cdot (y_{\gamma ,c}-\phi ))-\text{sign}(\bar{\lambda }+c\cdot (y_c-\phi )){\Vert }_{L^1(\Omega )}\to 0\text{for }\gamma \to \infty .$$

$$\begin{array}{rl}a(p_c,v)+c\cdot {\int }_{\Omega }\text{mathbbm}{1}_{A_c}\cdot p_c\cdot vdx& =-{\int }_{\Omega }(y_c-\bar{y})\cdot vdx\forall v\in H_0^1(\Omega ).\end{array}$$

$$p=\{\begin{array}{l}\tilde{p}\text{ in }\Omega \setminus A\\ 0\text{ in }A\end{array},$$

$$\begin{array}{rl}{a}_{\Omega \setminus A}(\tilde{p},v)& =-{\int }_{\Omega \setminus A}(y-\bar{y})vdx\forall v\in H_0^1(\Omega \setminus A)\end{array}$$

$$\begin{array}{rl}{a}_{\Omega \setminus A}:& H_0^1(\Omega \setminus A)\times H_0^1(\Omega \setminus A)\to \mathbb{R}\\ & (\tilde{p},v)\mapsto {\int }_{\Omega \setminus A}\sum _{i,j}{a}_{i,j}{\partial }_i\tilde{p}{\partial }_{j}v+\sum _i{d}_i({\partial }_i\tilde{p}v+\tilde{p}{\partial }_{i}v)+b\tilde{p}vdx\end{array}$$

$$\Vert {\text{sign}}_{\gamma }(\bar{\lambda }+c\cdot (y_{\gamma ,c}-\phi ))-\text{sign}(\bar{\lambda }+c\cdot (y_c-\phi )){\Vert }_{L^p(\Omega )}$$

$$\le$$

$$H_0^1(\Omega )$$

$$\Vert S_{\gamma ,c}-S_c{\Vert }_{\text{op}}$$

Full text

Key words: Semi-smooth optimization, variational inequality, obstacle problem, shape optimization, numerical methods, adjoint methods.

AMS subject classifications: 65K15, 49Q10, 49M29, 35Q93, 35J86, 49J40.

1 Introduction

We consider shape optimization problems constrained by variational inequalities (VI) of the first kind, so-called obstacle-type problems. Applications are manifold and arise, whenever a shape is to be constructed in a way not to violate constraints for the state solutions of partial differential equation depending on a geometry to be optimized. Just think of a heat equation depending on a shape, where the temperature is not allowed to surpass a certain threshold. This example is basically the model problem that we are formulating in section 2. Applications of general VI’s include contact problems in solid state mechanics, viscoplasticity and network equilibrium problems, and thus a wide range of industrial problems (cf. [38, 1, 33, 14]).

Shape optimization problem constraints in the form of VIs are challenging, since classical constraint qualifications for deriving Lagrange multipliers generically fail. Therefore, not only the development of stable numerical solution schemes but also the development of suitable first order optimality conditions is an issue.

By usage of tools of modern analysis, such as monotone operators in Banach spaces, significant results on properties of solution operators of variational inequalities have been achieved since the 1960s (cf. [6, 7, 28]). However, there are only very few approaches in literature to the problem class of VI constrained shape optimization problems so far. In [27], shape optimization of 2D elasto-plastic bodies is studied, where the shape is simplified to a graph such that one dimension can be written as a function of the other. The non-trivial existence of solutions of VI constrained shape optimization problems is discussed in [10, 44]. E.g., in [44, Chap. 4], shape derivatives of elliptic variational inequality problems are presented in the form of solutions to again variational inequalities. In [35], shape optimization for 2D graph-like domains are investigated. Also [29, 30] present existence results for shape optimization problems which can be reformulated as optimal control problems, whereas [12, 16] show existence of solutions in a more general set-up. In [36, 37], level-set methods are proposed and applied to graph-like two-dimensional problems. Moreover, [20] presents a regularization approach to the computation of shape and topological derivatives in the context of elliptic variational inequalities and, thus, circumventing the numerical problems in [44, Chap. 4]. Recently, in [17], a sensitivity analysis is performed for a class of semi-linear variational inequalities and a strong convergence property is shown for the material derivative. Furthermore, state-shape derivatives are established under regularity assumptions.

In this paper, we aim at optimality conditions for VI constrained shape optimization in the flavor of optimality conditions for VI constrained optimal control problems as in [18, 19, 21]. In general, standard necessary optimality conditions cannot be formulated in a straightforward manner for semi-smooth shape optimization problems. Under appropriate assumptions, we prove existence of adjoints and convergence of adjoints resulting from regularized variational inequalities. These analytical results are also verified numerically. Moreover, convergence of shape derivatives related to the smoothed problem is shown and the limit object is identified. Furthermore, we build on the resulting optimality conditions and devise an optimization algorithm giving specific numerical results. This algorithm does no longer depend on smoothing strategies as in [15]. In [15], a shape optimization method based on a regularized variant of the variational inequality has been devised and observed that the performance of this algorithm strongly depends on the tightness of the obstacle. This problem does no longer arise with the strategy developed in the present paper. On the contrary, the algorithms gets even faster, the more degrees of freedom are constrained by the obstacle.

This paper is structured as follows. In section 2, we formulate the VI constrained shape optimization model with general elliptic coefficients on which we focus in this paper. The necessary optimality conditions, including the existence of adjoint variables under certain regularity assumptions to the model problem are formulated in section 3. In section 4, we formulate an algorithm to solve the model problem based on these analytical results and compare numerically this approach with several regularized strategies.

2 Problem class

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded domain equipped with a sufficiently smooth boundary $\partial\Omega$, where $n\in\mathbb{N}$ is the dimension. For typical applications $n=2$ or $n=3$. This domain is assumed to be partitioned in a subdomain $\Omega_{\text{out}}\subset\Omega$ and an interior domain $\Omega_{\text{int}}\subset\Omega$ with boundary $\Gamma_{\text{int}}:=\partial\Omega_{\text{int}}$ such that $\Omega_{\text{out}}\sqcup\Omega_{\text{int}}\sqcup\Gamma_{\text{int}}=\Omega$, where $\sqcup$ denotes the disjoint union. The closure of $\Omega$ is denoted by $\bar{\Omega}$. We consider $\Omega$ depending on $\Gamma_{\text{int}}$, i.e., $\Omega=\Omega(\Gamma_{\text{int}})$. Figure 1 illustrates this situation. In the following, the boundary $\Gamma_{\text{int}}$ of the interior domain is called the interface and an element of an appropriate shape space $\mathcal{X}$ (cf. remark 1). In contrast to the outer boundary $\partial\Omega$, which is assumed to be fixed, the inner boundary $\Gamma_{\text{int}}$ is variable. If $\Gamma_{\text{int}}$ changes, then the subdomains $\Omega_{\text{int}},\Omega_{\text{out}}\subset\Omega$ change in a natural manner.

Let $\nu>0$ be an arbitrary constant. For the objective function

[TABLE]

we consider the following shape optimization problem:

[TABLE]

constrained by the following obstacle type variational inequality:

[TABLE]

where $y\in K$ is the solution of the VI, $f\in L^{2}(\Omega)$ is explicitly dependent on the shape, $\left<\cdot,\cdot\right>$ denotes the duality pairing and $a(\cdot,\cdot)$ is a general strongly elliptic, i.e. coercive, symmetric bilinear form

[TABLE]

defined by coefficient functions $a_{i,j},d_{j},b\in L^{\infty}(\Omega)$, fulfilling the weak maximum principle. However, the results of this paper still remain correct if symmetry of $a(.,.)$ is dropped as an assumption by simple modifications of proofs.

With the tracking-type objective $\mathcal{J}$ the model is fitted to data measurements $\bar{y}\in H^{1}(\Omega)$. The second term $\mathcal{J}_{\text{reg}}$ in the objective function $J$ is a perimeter regularization. A perimeter regularization is frequently used to overcome ill-posedness of inverse problems, e.g., [3] investigates the regularization and numerical solution of geometric inverse problems related to linear elasticity. In eq. 3, $\varphi$ denotes an obstacle which needs to be an element of $L^{1}_{\text{loc}}(\Omega)$ such that the set of admissible functions $K$ is non-empty (cf. [44]). If additionally $\partial\Omega$ is Lipschitzian and $\varphi\in H^{1}(\Omega)$ with $\varphi_{|\partial\Omega}\geq 0$, then there is a unique solution to eq. 3 satisfying $y\in H_{0}^{1}(\Omega)$, given that the assumptions from above hold (cf. [22, 9, 46]). Further, eq. 3 can be equivalently expressed as

[TABLE]

[TABLE]

with $(\cdot,\cdot)_{L^{2}(\Omega)}$ denoting the $L^{2}$-scalar product and $\lambda\in L^{2}(\Omega)$.

It is well-known, e.g., from [9], that under these assumptions there exists a unique solution $y$ to the obstacle type variational inequality (3) and an associated Lagrange multiplier $\lambda$. The existence of solutions of any shape optimization problem is a non-trivial question. Shape optimization problems constrained by VIs are especially challenging because, in general, it is not guaranteed that an adjoint state can be introduced (cf. [44, Example in Chap. 1, Chap. 4]). An essential theoretical tool for the study of the existence of solutions is the derivation of optimality conditions, i.e., in particular, the formulation of an adjoint equation. Therefore, section 3 investigates the model problem analytically, also in view of formulating a numerically applicable algorithm in section 4.

Remark 1**.**

The interface $\Gamma_{\text{int}}$ is an element of an appropriate shape space. Please note that there exists no common shape space suitable for all applications. It should be mentioned, that the existence of shape derivatives and their form is not dependent on the explicit choice of a shape space, hence only requirements noted in the according theorems are necessary. From a computational point of view one has to deal with polygonal shape representations arising in the setting of constrained shape optimization. This is owed to the fact that finite element methods usually discretize the models. In this paper, we use Steklov-Poincaré metric as introduced in [42]. These metrics can be considered e.g. on the space $B_{e}$ (cf. [32]), or more generally on the space of $H^{1/2}$-shapes (cf. [48]). In [43], it is outlined that this is an essential step towards applying efficient FE solvers. Of course, it is possible to choose other shape space models, but this is beyond the scope of this paper.

3 Convergence results for adjoints and shape derivatives

We assume the situation mentioned in section 2, which is also found in [23], giving us $\lambda\in L^{2}(\Omega)$. It can be easily verified that this in turn gives the possibility to summarize the conditions (6) equivalently into a single condition of the form

[TABLE]

The direct handling of general obstacle-type variational inequalities formulated as in (5)-(6), with eq. 6 being equivalently substitutes by eq. 7, poses several challenges. One challenge of the solution of eq. 5 is the occurrence of distributional numerical iterates for $\lambda$ in $H^{-1}(\Omega)$ when an augmented Lagrangian approach is applied to eq. 5 constrained by eq. 7, despite the analytical solution $\lambda$ having $L^{2}(\Omega)$-regularity. For a more detailed discussion of this, see [23, p.2]. In order to circumvent the occurrence of distributions in the solution scheme, the authors of [23] introduce a relaxation for relation eq. 7 with a given regularization parameter $\alpha\in(0,1)$

[TABLE]

which in turn is equivalent to

[TABLE]

if $\bar{\lambda}=0$ and $c=\frac{\varrho\alpha}{1-\alpha}$, where $\bar{\lambda}\in L^{2}(\Omega)$ can be motivated by updates of the augmented Lagrangian. This results in the equation

[TABLE]

which in the following is called regularized state equation or relaxed obstacle problem. Explicit dependence on $\lambda$ is avoided, making the resulting semi-linear elliptic equation tractable, for example by semi-smooth Newton methods, see, e.g, [23]. Moreover, the authors of [23] prove $L^{2}$-convergence of the regularized multiplier \max\big{(}0,\bar{\lambda}+c\cdot(y_{c}-\varphi)\big{)} to the original $\lambda$ for their proposed semi-smooth Newton method.

With problem (10) we are still left to solve a nonlinear, semi-smooth problem, giving rise to problems concerning existence of adjoints for the shape optimization problem. Hence, standard smoothing strategies can be applied to render this problem smooth enough to show existence of adjoints and to apply techniques such as Newton iterations.

In light of [41] and [8], we pose the following assumptions on the smoothed $\max$-function, which from now on is called $\max_{\gamma}\colon\mathbb{R}\rightarrow[0,\infty)$, with $\gamma>0$ being the smoothing parameter:

Assumption 1** (on smoothed $\max$-function).**

1. (i)

$\max_{\gamma}\in C^{1}(\Omega)$* for all $\gamma>0$;* 2. (ii)

there exists a function $g\colon(0,\infty)\rightarrow[0,\infty)$ with $g(\gamma)\rightarrow 0$ as $\gamma\rightarrow\infty$, s.t. $|\max_{\gamma}(x)-\max(0,x)|\leq g(\gamma)$ for all $x\in\mathbb{R}$ and for all $\gamma>0$; 3. (iii)

$\max_{\gamma}^{\prime}(x)\in[0,1]$* and monotonically nondecreasing for all $x\in\mathbb{R}$ and all $\gamma>0$;* 4. (iv)

$\max_{\gamma}^{\prime}$* converges uniformly to [math] on $(-\infty,-\delta)$ and $1$ on $(\delta,\infty)$ for all $\delta>0$ for $\gamma\rightarrow\infty$.*

In the following, let $\text{sign}_{\gamma}$ denote the derivative of $\max_{\gamma}$. An example satisfying these assumptions is given in (41). Applying $\max_{\gamma}$ instead of $\max$ in (10) gives the following equation, which we call fully regularized state equation in the subsequent chapters:

[TABLE]

So linearizing the corresponding Lagrangian with respect to $y_{\gamma,c}$ results in the typical adjoint equation

[TABLE]

(see, e.g., [18] or [41] in the context of optimal control).

Remark 2**.**

As in [41], smoothness of the state equation (11) in $y_{\gamma,c}$ guarantees existence of solutions to the linearized equation (12) for a given $L^{2}(\Omega)$ right-hand side and, thus, existence of adjoints in the case of the considered tracking-type objective functional (1).

3.1 State and adjoint equation

We first show that solutions of eq. 11 converge strongly in $H^{1}$ to solutions of (5)-(6) for $\gamma,c\rightarrow\infty$. This is proven in [41] for stronger assumptions on the smoothed function $\max_{\gamma}$ and under $\gamma=c$ . Since we rely on the general case $\gamma\neq c$ for the proofs in ongoing discussions, we state an according result. The first part of the following theorem is in analogy to [8, Lemma 4.2]. However, the difference is that we consider general elliptic bilinear forms and—more importantly—a modified argument in the maximum function resulting in different regularized state equations. These generalizations are necessary for our further analytical investigations leading to an adjoint equation.

Proposition 1** ($H^{1}$-convergence of the state).**

Let $y_{\gamma,c}$, $y_{c}$ and $y$ be solutions to eq. 11, eq. 10 and eq. 5, respectively. Here, $a(\cdot,\cdot)$ is chosen by an elliptic bilinear form as in (4) on a bounded, open domain $\Omega\subset{\mathbb{R}}^{n}$ with Lipschitz boundary, $f\in L^{2}(\Omega)$ and $\gamma,c>0$. Moreover, assume $\varphi\in H^{1}(\Omega)$, $\bar{\lambda}\in L^{2}(\Omega)$ and let $\max_{\gamma}\colon\mathbb{R}\rightarrow\mathbb{R}$ satisfy 1.

Then eq. 11 and eq. 10 posses unique solutions and

[TABLE]

Proof.

We prove statement (13) of the theorem. For a proof of statement (14), we refer to [23, Theorem 3.1].

We start by ensuring the existence of solutions to eq. 11 and eq. 10. For this, we show that the Nemetskii-operator defined by

[TABLE]

is a monotone operator for all $\gamma,c>0$. Due to 1, it is clear that $\text{max}_{\gamma}\colon\mathbb{R}\rightarrow\mathbb{R}$ is a point-wise monotone function, implying that $\text{max}_{\gamma}\colon H^{1}(\Omega)\rightarrow H^{1}(\Omega),y\mapsto\text{max}_{\gamma}(y)$ is a monotone operator. Since

[TABLE]

is an affine linear operator, and, thus monotone, the composition $\text{max}_{\gamma}\circ\Psi_{c}=\Phi_{\gamma}$ is also monotone. The same argument holds for the non-smoothed operator

[TABLE]

Therefore, applying the Browder-Minty theorem for monotone operators yields the existence of unique solutions to eq. 11 and eq. 10 in $H^{1}(\Omega)$ for all $f\in L^{2}(\Omega)$ if $\Omega$ is bounded and we operate in Hilbert spaces.

Now, we prove the second convergence (13). For fixed $c>0$, let $y_{\gamma,c}$ and $y_{c}$ be solutions to eq. 11 and eq. 10, respectively. Assumption 1 (ii) together with the monotonicity of $\Phi$, the coercivity of $a(\cdot,\cdot)$ with constant $K>0$ and $y_{\gamma,c}-y_{c}\in H^{1}(\Omega)$ acting as a test-function yields

[TABLE]

which gives the desired convergence (13). ∎

The following definition is needed to state the first main result of this paper, the convergence of adjoints.

Definition 1**.**

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded, open domain with Lipschitz boundary. A set $A\subseteq\Omega$ is called regularly decomposable, if there exists an $N\in\mathbb{N}$ and path-connected, bounded and open $A_{i}\subset\Omega$ with Lipschitz boundaries $\partial A_{i}$ such that $A=\overset{N}{\underset{i=1}{\sqcup}}\bar{A_{i}}$ is a disjoint union.

With this definition it is possible to formulate the first main theorem concerning the convergence of adjoints corresponding to the fully regularized problems and characterization of the limit object.

Theorem 1** (Convergence of the adjoints).**

Let $\Omega\subset\mathbb{R}^{n}$ for $n\leq 4$ be a bounded, open domain with Lipschitz boundary. Moreover, let the following assumptions are satisfied:

• (i)

$\varphi\in H^{1}(\Omega)$, $f\in L^{2}(\Omega)$, $\bar{y}\in H^{1}(\Omega)$ and coefficient functions $a_{i,j},d_{j},b\in L^{\infty}(\Omega)$ in (5)-(6);

• (ii)

the active set $A=\{x\in\Omega\;|\;y-\varphi\geq 0\}$ corresponding to (5)-(6) is regularly decomposable;

• (iii)

$A_{c}:=\{x\in\Omega\;|\;\bar{\lambda}+c\cdot(y_{c}-\varphi)\geq 0\}$* is regularly decomposable and*

[TABLE]

where $y_{c}$ solves the regularized state equation (10);

• (iv)

the following convergence holds:

[TABLE]

Then the adjoints $p_{\gamma,c}\rightarrow p_{c}$ in $H^{1}_{0}(\Omega)$ for $\gamma\rightarrow\infty$ for all $c>0$, where $p_{c}$ is the solution to

[TABLE]

Moreover, there exists $p\in H^{-1}(\Omega)$ to (5)-(6) and $p$ is representable as an $H^{1}_{0}$-function given by the extension of $\tilde{p}\in H^{1}_{0}(\Omega\setminus A)$ to $\bar{\Omega}$, i.e.,

[TABLE]

where $\tilde{p}\in H^{1}_{0}(\Omega\setminus A)$ is the solution of the elliptic problem

[TABLE]

with

[TABLE]

being the restriction of bilinear form $a(\cdot,\cdot)$ to $\Omega\setminus A$.

*Further, the solutions $p_{c}$ of eq. 18 converge strongly in $H^{1}_{0}(\Omega)$ to the

$H^{1}_{0}$-representation of $p$.*

Proof.

Let us consider the regularized problem (11) for $\gamma,c>0$. Existence and uniqueness of solutions $y_{\gamma,c},p_{\gamma,c}$ of regularized state and adjoint are guaranteed by application of the Minty-Browder Theorem in analogy to theorem 1 for $y_{\gamma,c}$ and the Lax-Milgram Theorem for $p_{\gamma,c}$, respectively.

This proof consists of two main parts:

Showing the $H^{1}$-convergence of the smoothed to the non-smoothed regularized adjoint $p_{\gamma,c}\rightarrow p_{c}$ for $\gamma\rightarrow\infty$.

• 2.

Analyzing the limit PDE (18) for $c\rightarrow\infty$ and proving that $p_{c}\rightarrow p$ in $H^{1}(\Omega)$ for $c\rightarrow\infty$, where $p$ is defined as in (19).

To 1. We start to show the $H^{1}$-convergence of the smoothed to the non-smoothed regularized adjoint $p_{\gamma,c}\rightarrow p_{c}$ for $\gamma\rightarrow\infty$.

The assumption (17) of $L^{1}$-convergence of $\textup{sign}_{\gamma}(\bar{\lambda}+c\cdot(y_{\gamma,c}-\varphi))$ is equivalent to $L^{p}$-convergence for all $p\in[1,\infty)$ in our setting, since

[TABLE]

by monotony of the integral and 1 (ii)-(iv). Denote by $S_{\gamma,c}\colon H_{0}^{1}(\Omega)\rightarrow H^{-1}(\Omega)$ the linear operator corresponding to the left-hand side of the smoothed adjoint equation (12) and $S_{c}\colon H^{1}(\Omega)\rightarrow H^{-1}(\Omega)$ the one to (18). We establish convergence of $S_{\gamma,c}$ to $S_{c}$ in the operator norm. In the following, we apply Hölder’s inequality and moreover, we use $L^{p}$-convergence of $\textup{sign}_{\gamma}(\bar{\lambda}+c\cdot(y_{\gamma,c}-\varphi))$ for all $p\in[1,\infty)$ as well as boundedness of $\textup{sign}_{\gamma}$ and sign.

Further, since we are in the situation $\Omega\subset{\mathbb{R}}^{n}$ for $n\leq 4$, we have the following embedding with embedding constant $C>0$ (cf. [40, Thm. 4.12 Part I, Case C])

[TABLE]

Combining all this yields

[TABLE]

which gives the desired convergence in the operator norm. Using analyticity of the inversion $\mathcal{I}\colon S\mapsto S^{-1}$ in the domain of invertible, bounded, linear operators given in our setting, convergence of the solution operators $S^{-1}_{\gamma,c}\rightarrow S^{-1}_{c}$ in operator norm is implied immediately, see e.g., [49, page 237]. Combining this with the convergence of $y_{\gamma,c}\rightarrow y_{c}$ in $H^{1}_{0}(\Omega)$ established by theorem 1 yields

[TABLE]

since $\|S^{-1}_{\gamma,c}\|_{\text{op}}$ can be bounded due to convergence.

To 2. Next, we analyze the limit PDE (18) for $c\rightarrow\infty$. We show that $p_{c}\rightarrow p$ in $H^{1}(\Omega)$ for $c\rightarrow\infty$, where $p$ is defined as in (19). For this, we first notice that our assumption concerning regular decomposability of $A=\{x\in\Omega\;|\;y-\varphi\geq 0\}$ ensures that $\partial A=\{x\in\Omega\;|\;y-\varphi=0\}$ forms a $C^{0,1}$-manifold embedded in $\Omega$. This in turn leads to well definedness of the restricted bilinear form $a_{\Omega\setminus A}(.,.)$ and the well-posedness of the variational problem (20) and, thus, of $p\in H^{1}_{0}(\Omega)$. Our next step is to show

[TABLE]

To show this, we artificially constrain problem (18) to $A\subseteq\Omega$. So denote by $a_{A}(\cdot,\cdot)$ the restriction of the bilinear form $a(\cdot,\cdot)$ to $A\subseteq\Omega$, defined in analogy to eq. 21. The corresponding restricted problem becomes

[TABLE]

where the Dirichlet condition $p_{c_{|A}}=p_{c}$ on $\partial A$ is incorporated in the usual way. Dividing by $c>0$ gives an equivalent equation in the sense that a solution $p_{c_{|A}}\in H^{1}(A)$ to eq. 24 also solves the equivalent equation

[TABLE]

The differential operator corresponding to the left-hand side of the equivalent equation (25) is given by :

[TABLE]

Next, we show that the differential operators $S_{A,c}$ converge in the linear operator norm $\|\cdot\|_{\text{op}}$ with the limit operator

[TABLE]

Whence

[TABLE]

due to eq. 22 and since $\text{vol}(A\setminus A_{c})\rightarrow 0\text{ for }c\rightarrow\infty$, which would otherwise contradict $y_{c}\rightarrow y$ in $H^{1}_{0}(\Omega)$. We can now apply a similar argument as in Step 1, namely the analyticity of the inversion operator $\mathcal{I}\colon S\mapsto S^{-1}$, giving us convergence of the solution operators $S_{A,c}^{-1}\rightarrow S_{A}^{-1}$ in $\|\cdot\|_{\text{op}}$. Also notice that we can obtain the sequence of solutions $p_{c_{|A}}$ by solving eq. 25 with the corresponding right hand sides $-\frac{1}{c}(y_{c}-\bar{y})$ instead of the original equation (24) and that the right hand sides converge to [math] in $H^{1}(\Omega)$ as $c\rightarrow\infty$, as $y_{c}$ is convergent by proposition 1. We conclude

[TABLE]

For the proof of convergence it remains to address the convergence of $p_{c}$ on $\Omega\setminus A$. We can artificially restrict eq. 18 to $\Omega\setminus A$ by imposing the Dirichlet boundary $p_{c|A}$ on $\partial A$, since $\partial A$ forms a $C^{0,1}$-submanifold of $\Omega$ as we assumed regular decomposability (cf. definition 1) of the active set $A$. To distinguish the corresponding bilinear forms, we denote the restricted bilinear form by $a_{\Omega\setminus A}$. Since the unrestricted bilinear form $a(\cdot,\cdot)$ is strongly elliptic, coercivity for some constant $K>0$ also holds for $a_{\Omega\setminus A}$. This together with Hölders inequality, assumption $A_{c}\subseteq A$ for all $c>0$ and the fact that $p_{c}-\tilde{p}\in H^{1}_{0}(\Omega\setminus A)$ can act as a testfunction gives

[TABLE]

where $\tilde{p}\in H^{1}(\Omega\setminus A)$ is defined as in (20). This results in

[TABLE]

due to our assumptions and $y_{c}\rightarrow y$ in $H^{1}(\Omega)$ as by proposition 1. Together with eq. 23 this gives the desired convergence $p_{c}\rightarrow p$ in $H^{1}_{0}(\Omega)$. ∎

There are a few non-trivial assumptions in theorem 1: assumption (iii) and (iv). In the following, we formulate two remarks in which we address these assumptions (cf. remark 3 for (iii) and remark 4 for (iv)).

Remark 3**.**

*It is possible to fulfill assumption (16) on inclusion of the active sets $A_{c}\subset A$ by choosing a sufficient $\bar{\lambda}\in L^{2}(\Omega)$. To be more precisely, if we assume $\varphi\in H^{2}(\Omega)$, we can choose $\bar{\lambda}:=\max\{0,f-S\varphi\}$ with $S$ being the differential operator corresponding to the elliptic bilinear form $a(\cdot,\cdot)$ in (5), guaranteeing feasibility $y_{c_{1}}\leq y_{c_{2}}\leq y\leq\varphi$ for all $0<c_{1}\leq c_{2}$. For the proof of this, we refer to [23, Section 3.2]. *

Remark 4**.**

Assumption 17 ensures that convergence of $\textup{sign}_{\gamma}$ is compatible with convergence of $y_{\gamma,c}$ for $\gamma\rightarrow\infty$. For giving a working example, we verify this assumption in the numerical section under eq. 45 for several demonstrative cases.

Remark 5**.**

The limit object $p\in H^{1}_{0}(\Omega)$ of the adjoints $p_{\gamma,c}$ as defined in (19) is the solution of an elliptic problem (20) on a domain $\Omega\setminus A$ with topological dimension greater than [math]. This can be exploited in numerical computations, for instance by a fat boundary method for finite elements on domains with holes as proposed by the authors of [34].

Remark 6**.**

We remind the reader, that $p$ is not necessarily an adjoint to the original problem eq. 1 constrained by eq. 5, eq. 6, but merely solution of a part of the limit of the optimality conditions for the regularized problem. For a discussion of a similar phenomenon in context of optimal control, we refer the interested reader to [8, Section 4.2].

3.2 Shape derivatives

In this section, we apply our convergence results for the regularized state and adjoint equations to derive similar convergence results for the shape derivatives of the shape optimization problem constrained by the fully regularized state equation (11). In general, shape derivatives of the unregularized VI constrained shape optimization problems do not exist (cf., e.g., [44, Chapter 1.1]). Nevertheless, we show existence of shape derivatives for the shape optimization problem constrained by the fully regularized VI eq. 11. Then, a limiting object corresponding to the unregularized equation eq. 5-eq. 6 is derived.

In the following, we split the main results into two theorems, the first one being the shape derivative for the fully regularized equation, the second one being convergence of the former for $\gamma,c\rightarrow\infty$.

The shape derivative of a general shape functional $H$ at $\Omega$ in direction of a sufficiently smooth vector field $V$ is denoted by $DH(\Omega)[V]$. For the definition of shape derivatives or a detailed introduction into shape calculus, we refer to the monographs [11, 44]. In general, we have to deal with so-called material and shape derivatives of generic functions $h\colon\Omega\to{\mathbb{R}}$ in order to derive shape derivatives of objective shape functions. For their definitions and more details we refer to the literature, e.g., [39]. In the following, we denote the material derivative of $h$ by $\dot{h}$ or $D_{m}(h)$ and the shape derivative of $h$ in the direction of a vector field $V$ is denoted by $h^{\prime}$.

Remark 7**.**

In this section, we only consider the shape functional $J$ defined in (1) without regularization term $\mathcal{J}_{\text{\emph{reg}}}$, i.e., we focus only on $\mathcal{J}$. The shape derivative of $J$ is given by the sum of the shape derivative of $\mathcal{J}$ and $\mathcal{J}_{\text{\emph{reg}}}$, where $D\mathcal{J}_{\text{\emph{reg}}}(\Omega)[V]=\nu\int_{\Gamma_{\text{\emph{int}}}}\kappa\left<V,n\right>ds$ with $\kappa:=\text{\emph{div}}_{\Gamma_{\text{\emph{int}}}}(n)$ denoting the mean curvature of $\Gamma_{\text{\emph{int}}}$. Please note that the objective functional and the shape derivative in correlation with the regularized VI (11) depends on the parameters $\gamma$ and $c$. In order to denote this dependency, we use the notation $\mathcal{J}_{\gamma,c}$ and $D\mathcal{J}_{\gamma,c}(\Omega)[V]$ for the objective functional and its shape derivative, respectively.

We state the first theorem, which presents the shape derivative of the objective functional $\mathcal{J}$ defined in (1) constrained by the fully regularized VI (11).

Theorem 2**.**

Assume the setting of the shape optimization problem formulated in section 2. Let the assumptions of theorem 1 hold. Moreover, let $M:=(a_{i,j})_{i,j=1,2,\dots,n}$ be the matrix of coefficient functions to the leading order terms in (4). Assume $y_{\gamma,c},p_{\gamma,c}\in W^{1,4}(\Omega)$, $a_{ij},b_{i},d\in L^{\infty}(\Omega)\cap W^{1,4}(\Omega)$ and $f\in H^{1}(\Omega)$. Furthermore, let $D_{m}(y_{\gamma,c}),D_{m}(p_{\gamma,c})\in H^{1}_{0}(\Omega)$ for all $\gamma,c>0$. Then the shape derivatives of $\mathcal{J}$ defined in (1) constrained by a fully regularized VI (11) in direction of a vector field $V\in H^{1}_{0}(\Omega,{\mathbb{R}}^{n})$ exist and are given by

[TABLE]

Proof.

Let us consider the shape optimization problem with fully regularized state equations with parameters $\gamma,c>0$ as in (11) and fixed shape $\Gamma_{\text{int}}$ to derive corresponding shape derivative. The first part of the proof, consisting of the existence of shape derivatives $D\mathcal{J}_{\gamma,c}$ for all $\gamma,c>0$, is found in appendix A.

As the second part of the proof, we derive the shape derivative expression. Note that eq. 29 and following integrals are well defined due to our assumptions on integrability combined with eq. 22 and dimension $n\leq 4$. By applying standard shape calculus techniques (cf. [4, 47]) to the target functional part of the Lagrangian we get

[TABLE]

since the target $\bar{y}\in L^{2}(\Omega)$ does not depend on the shape. Next, as similarly found in, e.g., [47], we calculate the shape derivative of the bilinear form $a(\cdot,\cdot)$. For avoiding confusion with the active sets $A$ and $A_{c}$, we call the coefficient matrix $(a_{i,j})_{i,j}$ of the leading order parts of the bilinear form $M$. As before we have

[TABLE]

We use linearity, chain rules, product rules and gradient identities for the material derivative $D_{m}(\cdot)$, as found in [4], to reformulate D_{m}\big{(}a(y_{\gamma,c},p_{\gamma,c})\big{)}. For readability, we analyze each term individually. We start with the leading order terms:

[TABLE]

For the first order terms of $a(\cdot,\cdot)$ we only compute $y_{\gamma,c}d^{T}\nabla p_{\gamma,c}$, since calculations are analogous for the second term by switching the roles of $y_{\gamma,c}$ and $p_{\gamma,c}$. We get

[TABLE]

where we again use shape independence of the coefficient functions of $a(\cdot,\cdot)$. For the term of order zero we apply the product rule for material derivatives and shape independence of coefficient functions:

[TABLE]

Combining these formulas, plugging them into eq. 30 and collecting all material derivatives of $y_{\gamma,c}$ and $p_{\gamma,c}$ result in the shape derivative of the bilinear form $a(\cdot,\cdot)$:

[TABLE]

The shape derivative of the term including $\max_{\gamma}$ is calculated by chain rule, which is applicable since we assume sufficient smoothness of $\max_{\gamma}$:

[TABLE]

[TABLE]

due to $D_{m}(\bar{\lambda})=0$ and $D_{m}(\varphi)=\nabla\varphi^{T}V$, as $\varphi$ is invariant under perturbation of the domain by problem definition. The shape derivative of the last term in the Lagrangian (50) is given by a simple product rule

[TABLE]

We now use the assumptions $D_{m}(y_{\gamma,c}),D_{m}(p_{\gamma,c})\in H^{1}_{0}(\Omega)$. If we rearrange the terms with $D_{m}(y_{\gamma,c})$ and $D_{m}(p_{\gamma,c})$ acting as test functions and applying the saddle point conditions, which means that the state equation (11) and adjoint equation (12) are fulfilled, the terms consisting $D_{m}(y_{\gamma,c})$ and $D_{m}(p_{\gamma,c})$ cancel. By adding all terms of eq. 50, the shape derivative $D\mathcal{J}_{\gamma,c}(\Omega)[V]$ as in (29) is established. ∎

Remark 8**.**

One can fulfill the assumptions of the averaged adjoint theorem and, thus, guarantee existence of the shape derivative without a computation of the material derivatives $D_{m}(y_{\gamma,c}),D_{m}(p_{\gamma,c})$. The assumption $D_{m}(y_{\gamma,c}),D_{m}(p_{\gamma,c})\in H^{1}_{0}(\Omega)$ for all $\gamma,c>0$ in theorem 2 is only needed in order to calculate the shape derivative expression (29).

Remark 9**.**

The assumption $D_{m}(y_{\gamma,c})\in H^{1}_{0}(\Omega)$ in theorem 2 is needed to apply the saddle point conditions and get the closed form of the shape derivative. For this, it is sufficient that $y_{\gamma,c}\in H^{2}_{0}(\Omega)$. For example, this regularity can be ensured by additionally assuming $a_{i,j}\in C^{1}(\bar{\Omega})$, $d\equiv 0$, $b\equiv 0$ for the coefficients of the strongly elliptic bilinear form $a(\cdot,\cdot)$ together with $\varphi>0$ and choosing $\bar{\lambda}\in L^{\infty}(\Omega)$. The latter two assumptions imply that the maximal monotone Nemetskii-operator (15) is equal to [math] for $y_{\gamma,c}=0$ and sufficiently large $\gamma,c>0$. In combination with the former assumptions, [5, Theorem A.1.] can be applied to get $y_{\gamma,c}\in H^{2}_{0}(\Omega)$ for all sufficiently large $\gamma,c>0$.

Remark 10**.**

The assumption $D_{m}(p_{\gamma,c})\in H^{1}_{0}(\Omega)$ in theorem 2 can be fullfilled, e.g., by assuming additional regularity $a_{i,j}\in C^{1}(\bar{\Omega})$ of the leading coefficients of the bilinear form $a(.,.)$ and $C^{2}$-regularity of the boundary $\partial\Omega$. This together with the fact that $c\cdot\textup{sign}_{\gamma}(\bar{\lambda}+c\cdot(y_{\gamma,c}-\varphi))\in L^{\infty}(\Omega)$ acts as part of the coefficient function of the zero order terms permits application of a regularity theorem for linear elliptic problems (cf. [13, p. 317, theorem 4]) giving $p_{\gamma,c}\in H^{2}(\Omega)$. This in turn guarantees $D_{m}(p_{\gamma,c})\in H^{1}_{0}(\Omega)$.

Next, we formulate the second main theorem of this section, which states the convergence of the shape derivatives of the fully regularized problem.

Theorem 3**.**

Assume the setting of the shape optimization problem formulated in (2). Let the assumptions of theorem 1 hold and $\varphi\in H^{2}(\Omega)$. Moreover, let $M:=(a_{i,j})_{i,j=1,2}$ be the matrix of coefficient functions to the leading order terms in (4). Then, for all $V\in H_{0}^{1}(\Omega,{\mathbb{R}}^{n})$, the shape derivatives $D\mathcal{J}_{\gamma,c}(\Omega)[V]$ in (29) converge to $D\mathcal{J}(\Omega)[V]$ for $\gamma,c\rightarrow\infty$, where

[TABLE]

Proof.

We see that (29) already resembles (32) except for the two terms

[TABLE]

We proceed in two steps: First, we show convergence for $T_{1}$ and $T_{2}$ as restricted operators on $C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$. Second, we show that the limiting operators can be continuously extended to $H^{1}_{0}(\Omega,\mathbb{R}^{n})$.

Let $V\in C^{\infty}_{0}(\Omega,{\mathbb{R}}^{n})$. By this, we have $\textup{div}(V)\cdot p_{\gamma,c},\nabla\varphi^{T}V\in H^{1}_{0}(\Omega)$ for all $\gamma,c>0$, which enables to use these functions as test functions for the state and adjoint equations. This leads to

[TABLE]

and

[TABLE]

due to theorem 1 proposition 1 and our assumption $\varphi\in H^{2}(\Omega)$.

Next, we lift the convergence from $V\in C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$ to $H^{1}_{0}(\Omega,\mathbb{R}^{n})$ by continuous extension. Since $C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$ is a dense subspace of $H^{1}_{0}(\Omega,\mathbb{R}^{n})$ and the latter being the completion of the former by the $\|\cdot\|_{H^{1}_{0}(\Omega,\mathbb{R}^{n})}$ norm, it is sufficient to show that the limits of $T_{1}(V_{n})$ and $T_{2}(V_{n})$ form a Cauchy sequence for a given Cauchy sequence $(V_{n})_{n\in\mathbb{N}}\subset C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$ under the $\|\cdot\|_{H^{1}_{0}(\Omega,\mathbb{R}^{n})}$-norm. So let $(V_{n})_{n\in\mathbb{N}}\subset C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$ with $\|V_{n}-V_{m}\|_{H^{1}_{0}(\Omega,\mathbb{R}^{n})}\rightarrow 0$ for $m,n\rightarrow\infty$. For the limit of $T_{1}$ we have

[TABLE]

[TABLE]

Here, we use integration by parts, Gauss’ Theorem, $p_{|\partial\Omega}=0$, $(V_{n}-V_{m})_{|\partial\Omega}=0$ and $L^{2}(\Omega)\hookrightarrow L^{1}(\Omega)$ with constant $C>0$ as in eq. 55. Thus, $(\tilde{T}_{1}(V_{n}))_{n\in\mathbb{N}}$ forms a Cauchy sequence and, therefore, gives a value for the continuous extension of $\tilde{T}_{1}$ for the limit of $V_{n}$ in $H^{1}_{0}(\Omega,\mathbb{R}^{n})$. For $T_{2}$ we use the same techniques, leading to

[TABLE]

With these convergences $T_{1},T_{2}$ converge to the continuously extended limit objects, which we from now on denote by $T_{1},T_{2}$, for all $V\in H^{1}_{0}(\Omega,\mathbb{R}^{n})$. Next, we simplify the sum of these two limiting objects. Let $V\in C^{\infty}_{0}(\Omega,\mathbb{R}^{n})$. Then

[TABLE]

where we use the definition of $p$, complementary slackness of $\lambda\in L^{2}(\Omega)$, test function properties of $\nabla\varphi^{T}V$ and $p\cdot\textup{div}(V)$, the state and adjoint equations. We apply again a continuity argument to gain this identity for all $V\in H^{1}_{0}(\Omega,\mathbb{R}^{n})$. We see that the limit object in (35) is exactly the missing term in the limit of the shape derivatives $D\mathcal{J}_{\gamma,c}(\Omega)[V]$ (cf. (32)). ∎

Remark 11**.**

Theorem 2 and theorem 3 are also valid when $f\in H^{1}(\Omega)$ or $\varphi\in H^{2}(\Omega)$ depend explicitly on the shape $\Omega$ with shape derivatives $f^{\prime},\varphi^{\prime}\in H^{1}_{0}(\Omega)$. Then the shape derivatives need to be modified accordingly by replacing terms including $\nabla f^{T}V$ and $\nabla\varphi^{T}V$ by $\nabla f^{T}V+f^{\prime}$ and $\nabla\varphi^{T}V+\varphi^{\prime}$. Further, Theorem 2 and theorem 3 remain valid for piecewise constant $f\in L^{\infty}(\Omega)$ depending on the shape $\Omega$ by adjusting the proofs applying integral splitting techniques as found in [47, Remark 4.21, Thrm. 4.23].

Remark 12**.**

It is common knowledge that by pushing the obstacle $\varphi$ to infinity, i.e., $\varphi(x)\uparrow\infty$ for all $x\in\Omega$, the state equation representing the variational inequality (5) becomes a regular elliptic PDE in weak formulation

[TABLE]

due to (6). This means that we encounter shape optimization problems with elliptic PDE constraints. Formula (32) remains valid by applying $A=\varnothing$, giving a shape derivative for a general elliptic problem.

Remark 13**.**

The limiting object eq. 32 is in general not the shape derivative of the unregularized problem. It can be regarded as part of the limit system arising during convergence of the optimality systems of the fully regularized problem. Finding a framework in shape optimization to describe the type of this limiting object will be part of further research and is beyond the scope of this article. For readers interested in a treatise on limiting systems of optimality systems in non-smooth optimal control we recommend [8].

Remark 14**.**

*The limiting objects of the convergence results for adjoint variables (cf. theorem 1) and shape derivatives (cf. theorem 2) can be put into relation by conditions resembling C-stationarity, e.g., as found in [19, Definition 4.1.].

Using our terminology, it is necessary for C-stationarity conditions to hold that a $\xi\in H^{-1}(\Omega)$ exists such that the adjoint equation can be formulated in the form*

[TABLE]

We can define such a $\xi\in H^{-1}(\Omega)$ by emulating the definition of $p$ in (19), including enforcement of the Dirichlet condition $p=0$ on $\partial A$ with Nitsche’s method using boundary terms (cf. [24]). The state equation, corresponding complementarity conditions, and the design equation, which in our setting can be viewed as the shape derivative identity (29), hold in analogy to the cited definition of C-stationarity. The remaining conditions

[TABLE]

by the definitions of $\xi$ and $p$, are satisfied as well. It is worth mentioning that—to knowledge of the authors—no type of C-stationarity-like conditions for optimality of VI constrained shape optimization problems have been investigated or defined before. By defining C-stationarity in this context, as outlined above, we can sum up the theorems by stating that the solutions of the regularized equations converge to a C-stationary system.

4 Algorithmic aspects and numerical investigations

In this section, we put the theoretical treatise highlighted in the previous section into numerical practice on domains in ${\mathbb{R}}^{2}$. We employ a steepest descent algorithm with backtracking linesearch in order to perform the optimization procedures with various regularized as well as unregularized versions of the specialized variational inequality (see (39)). Also, we propose a way to incorporate the unregularized approach in an algorithm and compare it to the different regularizations.

For convenience, we specialize the more general constraint (5) to a Laplacian version:

[TABLE]

We use $\nu=10^{-5}$ for all computations in this section. As the right-hand side of the state equation in eq. 38 we choose following piecewise constant function $f\in L^{\infty}(\Omega)$ defined by

[TABLE]

For calculations of the smoothed state and adjoint we have to specify $\max_{\gamma}$ satisfying 1. For demonstrative purpose, we choose a similar smoothing procedure as in [23, Section 2]:

[TABLE]

A different, more regular smoothing is, e.g., given in [41, (1.10)]. Both smoothing techniques mentioned satisfy 1. For the sake of completeness, we also give the first derivative formula

[TABLE]

In this setting, the shape derivative (32) simplifies to

[TABLE]

and analogously the shape derivative for the fully regularized equation in eq. 29. Notice that the shape derivative of the perimeter regularization is also included in our computations (cf. remark 7).

In the following numerical experiments, we consider two different obstacles:

[TABLE]

The calculations are performed with Python using the finite element package FEniCS. For detailed informations on FEniCS, we refer to [2] and [31]. As initial shape we choose a centered circle with radius $0.15$, illustrated in fig. 7. The computational grid of the initial shape, which is embedded in the hold-all-domain $(0,1)^{2}\subset\mathbb{R}^{2}$, consists of $2\,184$ vertices with $4\,206$ cells, having a maximum cell diameter of 0.0359 and a minimum cell diameter of 0.018. The algorithm employed for the shape optimization is summarized in algorithm 2. In the following, we describe the algorithm and the chosen parameters in detail.

The target data $\bar{y}\in L^{2}(\Omega)$ is computed by using the mesh of the target interface to calculate a corresponding state solution of eq. 39 by the semi-smooth Newton method proposed in [23]. These are visualized in fig. 2 for both obstacles $\varphi_{1}$ and $\varphi_{2}$. We apply the same method for calculating state variables $y$ in the unregularized optimization approach.

For the regularized and smoothed states $y_{\gamma,c}$ and $y_{c}$ we use a Newton- and semi-smooth Newton method provided by the FEniCS package in order to solve the linear systems assembled by using first order polynomials on the computational grids. All state calculations in our routines are performed with a stopping criterion of $\varepsilon_{\text{state}}=3.e-4$ for the error norms. In light of remark 3 we choose $\bar{\lambda}=\textup{max}\{0,f+\Delta\varphi\}$, which is possible due to sufficient regularity of $\varphi_{1}$ and $\varphi_{2}$.

To ensure assumptions of theorem 1, theorem 2 and theorem 3, it is necessary to fulfill

[TABLE]

We calculate the corresponding norm using various $c>0$ and both, $\varphi_{1}$ and $\varphi_{2}$, on refined meshes having $212\,642$ vertices, $423\,682$ cells and maximum and minimum cell diameter of $0.0038$ and $0.0015$, respectively. An example convergence plot can be found in fig. 3.

We want to point out that as $\gamma\rightarrow\infty$, the norm in (45) converges to an $\varepsilon>0$ which is close to [math]. This is due to numerical errors resulting from the state equation, since their solution determines the active set, which is needed to calculate the values of sign and $\textup{sign}_{\gamma}$. The functions, whose $L^{1}$-norms are of interest, are illustrated in fig. 4 on a refined mesh. Furthermore, we observe that these functions, and hence the errors go to [math] for ever finer grid widths and more precisely calculated states $y_{\gamma,c},y$. This is supported by a study successively evaluating the mentioned $L^{1}$-norms in the circular start shape, as depicted in fig. 7, on meshes generated by adaptive refinement at the boundary of the active set for large $\gamma,c>0$. The results are seen in fig. 5.

The adjoints $p_{\gamma,c}$ and $p_{c}$ are calculated by solving eq. 12 and eq. 18 with first order elements by using the FEniCS standard linear algebra back end solver PETSc.

Calculating the limit $p$ of the adjoints $p_{\gamma,c}$ as in eq. 19 and eq. 20 are performed in several steps. First, a linear system corresponding to

[TABLE]

is assembled without incorporation of information from the active set $A$. Afterwards, the vertex indices corresponding to the points in the active set $A=\{x\in\Omega\;|\;y-\varphi\geq 0\}$ are collected by checking the condition

[TABLE]

for some error bound $\varepsilon_{\text{adj}}>0$. The error bound $\varepsilon_{\text{adj}}$ is incorporated since $y$ is feasibly approximated by $y_{i}$ with the semi-smooth Newton method from [23], i.e., $y_{i}\leq\varphi$ for all $i\in\mathbb{N}$. After this, the collected vertex indices are used to incorporate the Dirichlet boundary conditions $p=0\quad\text{in }A$ into the linear system corresponding to (46). To solve the resulting system, we use the same procedures as to solve for $p_{\gamma,c},p_{c}$, i.e.m the standard PETSc back end conjugate gradient solver. An exemplary solution $p$ of the unregularized adjoint equation is illustrated in fig. 6. We want to point out that the active set and consequently the zero level set resulting from the Dirichlet conditions can be observed in fig. 6.

To calculate gradients $U\in H^{1}_{0}(\Omega,\mathbb{R}^{2})$ used in the steepest descent method for solving eq. 38, we use a Steklov-Poincaré metric induced by the linear elasticity equation, as proposed in [42]. In particular, we assemble the shape derivatives given in theorem 2 and theorem 3 as the right-hand side of the linear elasticity equation

[TABLE]

with the so called Lamé-parameters $\lambda_{\text{elas}}$ and $\mu_{\text{elas}}$. Here, we choose $\lambda_{\text{elas}}=0$ and $\mu_{\text{elas}}$ as the solution of the Poisson problem

[TABLE]

for $\mu_{\text{max}},\mu_{\text{min}}>0$. As a physical interpretation, this enables to control stiffness of the grid by choosing $\mu_{\text{max}}$ and $\mu_{\text{min}}$ in order to influence $\mu_{\text{elas}}$ acting as a coefficient function in the linear elasticity equation (48). Thus, larger values of $\mu_{\text{max}}$ lead to more stiffness at the interface $\Gamma$ and larger values of $\mu_{\text{min}}$ to more stiffness at the boundary $\partial\Omega$ of the hold-all domain $\Omega$. For our calculations, we choose $\mu_{\text{min}}=0$ and $\mu_{\text{max}}=25$ for $\varphi_{1}$ and $\mu_{\text{max}}=55$ for $\varphi_{2}$. It is important to notice that we set all right-hand side values of (48) which do not have a neighboring vertex on the interface to [math]. For a more detailed discussion of this we refer to [42].

To complete the description of our optimization we shortly explain the linesearch we will employ in our numerical calculations. We use a simple backtracking linesearch with sufficient descent criterion, where $U_{k}$ denotes the shape derivative calculated at the corresponding interface in $\Omega_{k}$ in step number $k$, $\mathcal{T}_{\tilde{U}}(\Omega_{k}):=\{y\in\mathbb{R}^{2}:y=x+\tilde{U}(x)\text{ for some }x\in\Omega_{k}\}$ the linearized vector transport by $\tilde{U}$ and $y_{\tilde{U}}$ the state solution in $\mathcal{T}_{\tilde{U}}(\Omega_{k})$.

We summarize our approach in algorithm 2 for the unregularized procedures. The regularized and smoothed procedures work analogously by modifying the state, adjoint and shape derivative equations. The calculations of $p_{\gamma,c},p_{c}$ are straightforward and need not the additional steps outlined in before and in algorithm 2 for the unregularized $p$. In the design of algorithm 2 a safeguarding technique is employed. This stems from the fact, that the limit of shape derivatives $D\mathcal{J}$ from eq. 32 is in general not the true shape derivative of the initial problem, see remark 14. Hence an additional testing the convergence criterion for the fully regularized shape derivative $\mathcal{D}\mathcal{J}_{\gamma,c}$ is performed after the convergence by $D\mathcal{J}$. If no convergence is detected by $D\mathcal{J}_{\gamma,c}$, $D\mathcal{J}_{\gamma,c}$ will be used to calculate a further descent direction, as the latter is a true shape derivative by theorem 2. Further, the safeguarding acts as a safety when the adjoint limit object $p_{k}$ is flawed due to erroneous allocation of the active set $A_{k}$ as discussed with eq. 47 in the beginning of this section. The smoothed model is not prone to this effect, hence acting as a substitute model for further gradient calculations.

In our calculations the safeguard was never activated by not finding a descent direction during the linesearch procedure, indicating that the shape derivative limiting object $D\mathcal{J}$ is acting appropriately for a shape derivative substitute, making the safeguard for this purpose obsolete. Still, the safeguarding is activated at convergence for coarse grids or imprecise calculations of the state $y_{k}$, indicating a non-neglectable difference in $\|D\mathcal{J}_{\gamma,c}\|$ and $\|D\mathcal{J}\|$. This is only due to false allocation of the active set $A_{k}$, resulting in inaccurate $p_{k}$ and $\|D\mathcal{J}(\Omega_{k})\|$. For grids with maximum cell diameter $0.01$ or less and error tolerance $\varepsilon_{\text{state}}<1.e-7$ for the state calculation, the errors in active set allocation are sufficiently small for $\|D\mathcal{J}(\Omega_{k})\|\approx\|D\mathcal{J}_{\gamma,c}(\Omega_{k})\|$ and the safeguarding to not be activated at convergence in all our examples.

Our findings concerning convergence of the various shape optimization approaches, using the unregularized approach for various $\varphi_{\text{adj}}$, as well as regularized approaches with different parameters $\gamma,c>0$, are displayed in fig. 8 for $\varphi_{1}=0.5$ and in fig. 9 for $\varphi_{2}=5e^{-x_{1}-1}$. Morphed shapes arising during the optimization procedure are plotted in fig. 7 for the unregularized approach using $\varepsilon_{\text{adj}}=10^{-9}$. It can be seen in the plots that there are vanishing difference between approaches using fully regularized calculation with sufficiently large $\gamma$ and $c$, regularized ones with large $c$ and the unregularized one. For smaller regularization parameters $\gamma$ and $c$, the solved state and adjoint equations begin to differ from the original problem and, thus, slow down convergence, and for very small $\gamma$ and $c$ no convergence at all.

The convergence behavior of the unregularized method strongly depends on the selection of the active set. When the state solution $y$ is not calculated with sufficient precision the numerical errors lead to misclassification of vertex indices. Hence wrong Dirichlet conditions are incorporated in the adjoint system, creating errors in the adjoint. This makes the gradient sensitive to error for smaller $\varepsilon_{\text{adj}}$, as can be seen by the slight roughness of the target graphs in fig. 8 and fig. 9 for $\varepsilon_{\text{adj}}=10^{-9}$ and $\varepsilon_{\text{adj}}=0.01$. In order to compensate this, the condition for checking active set indices (47) can be relaxed by increasing $\varepsilon_{\text{adj}}$. This increases likelihood of correctly classifying the true active indices, while also increasing likelihood of misclassification of inactive indicies. Such a relaxation can lead to errors in the adjoint increasing with $\varepsilon_{\text{adj}}$ and, thus, trading convergence speed for robustness, also visible in fig. 8 and fig. 9. Of course, this gets less feasible for highly oscillatory obstacle $\varphi$ and state $y$, as well as state solves with high tolerance $\varepsilon_{\text{\text{state}}}$.

In order to circumvent this, it is obviously sufficient to decrease error tolerance $\varepsilon_{\text{state}}$ of the state calculation. An exemplary result of this can be seen in fig. 9 under unregularized*, where we decreased the error tolerance to $\varepsilon_{\text{state}}=4.e-5$. Nevertheless, additional decrease of $\varepsilon_{\text{state}}$ comes with more computational cost, whereas with increase of $\varepsilon_{\text{adj}}$ the robustness is paid by loss of convergence speed.

It is worth to mention that implementing the unregularized state and adjoint becomes especially numerically exploitable with higher resolution meshes and more strongly binding obstacles $\varphi$, i.e., larger active sets $A$. This is possible by sparse solvers due to the incorporation of Dirichlet conditions on the active set, as we have proposed, or by a fat boundary method as in [34]. This is especially advantageous for large systems resulting from fine resolution grids, as exploitability of sparsity and accuracy of our method increase at the same time.

So in contrast to the method proposed in [15], where performance slows down for more active obstacle $\varphi$, we do not notice unusual slowdown in performance with the methods proposed in this article, and even offer possibility to actually benefit numerically from more binding obstacle $\varphi$.

5 Conclusion

Shape optimization for variational inequalities is more challenging than both, elliptic shape optimization and optimal control for variational inequalities. In this paper, we derive optimality conditions for shape optimization in the context of variational inequalities in the flavor of optimal control problems. Regularized variants are studied and limiting conditions derived. This gives rise to highly efficient optimization algorithms. In the future general investigations of necessary optimality criteria for VI constrained shape optimization like C-stationarity are conceivable. Also large-scale multidimensional computational comparisons of the presented method in comparison to other state-of-the-art methods is of particular interest.

Acknowledgement

The authors are indebted to Leonhard Frerick (Trier University) for many helpful comments and discussions about functional analytical aspects of convergence. Furthermore, the paper profited form remarks of Gerd Wachsmuth on an earlier version. This work has been partly supported by the German Research Foundation (DFG) within the priority program “Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and Hierarchical Optimization” SPP 1962/1 and SPP 1962/2 under contract numbers Schu804/15-1 and WE 6629/1-1, by the DFG research training group 2126 Algorithmic Optimization, and by the BMBF (Bundesministerium für Bildung und Forschung) within the collaborative project GIVEN (FKZ: 05M18UTA).

Appendix A Proof of existence of shape derivatives in theorem 2

What follows is a proof for the existence of shape derivatives $D\mathcal{J}_{\gamma,c}$ for all $\gamma,c>0$ under the assumptions of theorem 2.

Proof.

For the proof of existence $D\mathcal{J}_{\gamma,c}$ for all $\gamma,c>0$, we will employ the so called averaged adjoint approach, as found in [26, Ch. 7, Thm. 7.1]. We roughly follow a proof found in [26, Ch. 7, Thm. 7.2], but have to give a proof ourselves, since the situation in this paper differs from the one mentioned previously, e.g. as we are not having a bounded Nemetskii operator as the non-linearity in the semi-linear state equation.

Let $\gamma,c>0$. By definition, the Lagrangian function corresponding to our problem is given by

[TABLE]

We have to proof the assumptions (H0) - (H3) needed to apply the averaged adjoint theorem as found in [26, Chapter 7, Thm. 7.1] in order to guarantee shape differentiability. For convenience, we do not state all the lengthy assumptions here, and refer the interested reader to [26, Chapter 7, Thm. 7.1]. Let us fix a deformation vector field $V$ and denote the domains deformed by $\Omega_{t}$ for deformation parameters $t\in[0,\tau]$ and some $\tau>0$ small enough, such that the corresponding deformations $\mathcal{T}_{t}$ is bijective. For the rest of the existence proof, we will drop $\gamma,c$ as the subscripts of $y_{\gamma,c}$ and $p_{\gamma,c}$ for readability purposes, still knowing we are in the fully regularized situation. Define

[TABLE]

for the deformed domain $\Omega_{t}$ resulting from application of the deformation $\mathcal{T}_{t}$ in direction $V$ parametrized by $t\in[0,\tau]$.

The first assumption (H0) concerns well behavedness of eq. 51. First notice that the function in eq. 51 is both differentiable in the state $y$ and adjoint $p$, which can be seen after applying the transformation theorem to eq. 51. Thus for the set

[TABLE]

we have

[TABLE]

with $y^{t}=y_{t}\circ\mathcal{T}_{t}$ being the retraction of the unique solution $y_{t}\subset H^{1}_{0}(\Omega_{t})$ of the fully regularized state equation eq. 11 on the deformed domain $\Omega_{t}$ as by the use of Minty-Browder’s Theorem as portrayed in the proof of proposition 1.

Further, the Lagrangian in direction $V$ as a function of averaged states inserted

[TABLE]

is absolutely continuous in $s\in[0,1]$. For this, we make use of the fact that

[TABLE]

for bounded $\Omega\subset{\mathbb{R}}^{n}$ with a constant depending on $p$ and $q$. So absolute continuity is satisfied, as we have existing derivatives of the integrand due to 1 (i) and integrability of the integrand due to

[TABLE]

for all $z\in L^{2}(\Omega)$ and all $z\in H^{1}_{0}(\Omega)$ as by 1 (ii), $\Omega$ being bounded and $H^{1}_{0}(\Omega)\hookrightarrow L^{4}(\Omega)\hookrightarrow L^{2}(\Omega)$ with constant $C$ by eq. 22 and eq. 55. Further, the directional derivative mapping

[TABLE]

is integrable for all $\tilde{p}\in H^{1}_{0}(\Omega)$, since

[TABLE]

by Hölder’s Inequality, 1 (iii) and $a_{t}(.,.)$ being the bilinear form defined by retraction of $a(.,.)$ from $\Omega_{t}$ to $\Omega$ bound with constants $B_{t}>0$. This, and the easy to verify fact that $G$ is a affine linear function in $p$, gives us (H0). We remind the careful reader, that the Jacobians created by retraction of $\Omega_{t}$ to $\Omega$ are to be implicitly included in scalarproducts and norms above for the calculations to be valid. We don’t explicitly state these for readability.

Next we introduce the set of so called averaged adjoints

[TABLE]

We manipulate the averaged adjoint equation found in eq. 58 by interchanging integrals

[TABLE]

which is an elliptic PDE with an additional positive $L^{\infty}(\Omega)$ coefficient function term for the zero’th order terms, where we again omitted explicit statement of Jacobians. By 1 $(iii)$, the additional coefficient term for the zero’th order terms in the averaged adjoint equation eq. 59 satisfies

[TABLE]

which results in coercivity and boundedness of the corresponding bilinear form of the averaged adjoint equation. This lets us apply the Lemma of Lax-Milgram, resulting in existence of a unique solution for the averaged adjoint equation we will denote by $q^{t}\in H^{1}_{0}(\Omega)$ for all $t\in[0,\tau]$. Thus we have the identity $Y(t,y^{t},y^{0})=\{q^{t}\}\subset H^{1}_{0}(\Omega)$, which together with eq. 53 ensures condition (H2).

We also notice, that the derivatives of $\frac{\partial}{\partial t}G$ exist and can be explicitly calculated after application of the transformation theorem, giving us (H1).

To apply the averaged adjoint theorem from [25], it remains to address (H3), which is satisfied in our case by application of [45, Lemma 4.1], if for the unique solutions of the state- and adjoint equation $y^{0}\in X(0)$ and $q^{0}\in Y(0,y^{0},y^{0})$ and a given sequence $(t_{n})_{n\in\mathbb{N}}\subset[0,\tau]$ converging to zero, we can find a subsequence $(t_{n_{k}})_{k\in\mathbb{N}}\subseteq(t_{n})_{n\in\mathbb{N}}$ with $q^{t_{n_{k}}}\in Y(t_{n_{k}},y^{t}_{n_{k}},y^{0})$ such that

[TABLE]

We will mimic parts of the argumentation found in the proof of [45, Theorem 5.1] accustomed to our situation, which is slightly different than the one found in [45, Theorem 5.1] or [26, Theorem 7.2].

Consider the solutions $y^{0}\in X(0)$ and $q^{0}\in Y(0,y^{0},y^{0})$ and a sequence

$(t_{n})_{n\in\mathbb{N}}\subset[0,\tau]$ converging to zero.

First notice that by monotony of the Nemetskii operator eq. 15 of the concerning semilinear state equation we have

[TABLE]

This in turn, together with the coercivity of the retracted bilinearform $a_{t}(.,.)$ with constant $K_{t}>0$, eq. 56 and choosing $y^{t}\in X(t)\subset H^{1}_{0}(\Omega)$ as a testfunction gives us

[TABLE]

again omitting Jacobians. Dividing by $\|y^{t}\|_{H^{1}_{0}(\Omega)}$, using the convergence $\mathcal{T}_{t}$ to the identity for $t\downarrow 0$ and by taking a supremum we achieve

[TABLE]

bounding the norms by a constant $0<M<\infty$ independent of $t\in[0,\tau]$. Recognize that the norms still implicitly depend on $t$, since Jacobians are to be included.

In the same line of argumentation we can confirm the boundedness of $\|q^{t}\|_{H^{1}_{0}(\Omega)}$. For this, we apply the first inequality of eq. 60 to get

[TABLE]

By finally using eq. 64 we arrive at

[TABLE]

As we have established bound eq. 64, we can choose a subsequence $(t_{n_{k}})_{k\in\mathbb{N}}\subseteq(t_{n})_{n\in\mathbb{N}}$, such that $y^{t_{n_{k}}}\rightharpoonup z$ weakly in $H^{1}_{0}(\Omega)$ for $k\rightarrow\infty$ and some $z\in H^{1}_{0}(\Omega)$.

Further, using the convergence of the retracted functions $\bar{\lambda}_{t}$ and $\varphi_{t}$ in $L^{2}(\Omega)$ for $t\downarrow 0$, we can uniformly and independently of $t_{n_{k}}$ bound

[TABLE]

by using eq. 56 and eq. 64. By having boundedness of the coercive bilinearforms $a_{t}(.,.)$, eq. 64 and smoothness in $\text{max}_{\gamma}$ by 1 (i) we are able to apply Lebesgue’s dominated convergence theorem to the retracted state equations, giving us $y_{t_{n_{k}}}\rightharpoonup y^{0}$ weakly in $H^{1}_{0}(\Omega)$ due to the unique solution guaranteed by Minty-Browder’s theorem.

Applying the same routine due to eq. 65, we can choose a subsequence of $\{t_{n_{k}}\}$, which we will again call $\{t_{n_{k}}\}$ by abuse of notation, such that $q^{t_{n_{k}}}\rightharpoonup u$ weakly in $H^{1}_{0}(\Omega)$ for some $u\in H^{1}_{0}(\Omega)$. Then uniform boundedness eq. 60, the previously established weak convergence $y_{t_{n_{k}}}\rightharpoonup y^{0}$ and eq. 65 yield applicability of Lebegue’s theorem for inserted $t_{n_{k}}$ in eq. 59. For $k\rightarrow\infty$, the limit equation of eq. 59 is the fully regularized adjoint equation eq. 12, which has a unique solution by Lax-Milgram’s lemma. Whence $q^{t_{n_{k}}}\rightharpoonup q^{0}=p_{\gamma,c}\in Y(0,y^{0},y^{0})$ weakly in $H^{1}_{0}(\Omega)$ with the previously established weak convergence of $q^{t_{n_{k}}}$ and continuity of $\text{sign}_{\gamma}$ by 1 (i).

Now we have found a subsequence $\{t_{n_{k}}\}\subseteq\{t_{n}\}$, such that $q^{t_{n_{k}}}\rightharpoonup q^{0}$ weakly in $H^{1}_{0}(\Omega)$. Using the transformation theorem, $G(t,y^{0},q^{t_{n_{k}}})$ from eq. 51 can be stated as an integral in $\Omega$ with integrands being differentiable in $t\in[0,\tau]$. The derivative $\frac{\partial}{\partial t}G(t,y^{0},q^{t_{n_{k}}})$ is weakly continuous in it’s first and last argument, hence the weak convergence $q^{t_{n_{k}}}\rightharpoonup q^{0}$ implies eq. 61, which is condition (H3) from [26, Theorem 7.1]. All assumptions (H0)-(H3) for the averaged adjoint theorem [26, Theorem 7.1] are satisfied, finally guaranteeing existence of shape derivatives $D\mathcal{J}_{\gamma,c}$ for all $\gamma,c>0$. ∎