Computational Aspects for Interface Identification Problems with Stochastic Modelling

TL;DR

This paper addresses interface identification problems under uncertainty by formulating a shape optimization constrained by stochastic PDEs, introducing a generalized stochastic gradient method, and demonstrating its effectiveness through numerical experiments.

Contribution

It introduces a novel stochastic gradient approach for shape optimization with uncertain PDE constraints, extending classical methods to stochastic shape spaces.

Findings

The proposed algorithm effectively solves stochastic interface identification problems.

Numerical experiments confirm the method's convergence and robustness.

The approach handles uncertainty in coefficients and inputs in PDE-based shape optimization.

Abstract

In this paper, a shape optimization problem constrained by a random elliptic partial differential equation with a pure Neumann boundary is presented. The model is motivated by applications in interface identification, where we assume coefficients and inputs are subject to uncertainty. The problem is posed as a minimization of the expectation of a random objective functional depending on the uncertain parameters. A numerical method for iteratively solving the problem is presented, which is a generalization of the classical stochastic gradient method in shape spaces. Moreover, we perform numerical experiments, which demonstrate the effectiveness of the algorithm.

Equations

$$B_e=B_e(S^1,{\mathbb{R}}^2):=\mathrm{Emb}(S^1,{\mathbb{R}}^2)/\mathrm{Diff}(S^1).$$

$$f(x,\omega )=f(x,{\xi }_f(\omega )),g(x,\omega )=g(x,{\xi }_g(\omega ))\text{on }D\times \Omega .$$

$$\kappa :D\times {\Omega }_{\kappa }\to \mathbb{R},(x,\omega )\mapsto \sum _{i=0}^N{\kappa }_i(\omega )\text{mathbbm}{1}_{D_i}(x).$$

$$J(u,\xi ):=\frac{1}{2}{\int }_D(y(x,\xi )-\bar{y}(x){)}^2\text{d}x$$

$$\displaystyle\underset{\text{with }u_{i}\in B_{e}}{\min_{u=\sqcup_{i=1,\dots,N}u_{i}}}\quad\left\{j(u):=\mathbb{E}\big{[}J(u,\xi)\big{]}\right\}$$

$$\text{s.t.}y:D\times \Xi \to \mathbb{R},(x,\xi )\mapsto y(x,\xi )\text{ satisfies}$$

$$-\nabla \cdot (\kappa \nabla y)$$

$$\kappa \frac{\partial y}{\partial n}$$

$$\left\llbracket\kappa\frac{\partial y}{\partial n}\right\rrbracket=0,\quad\llbracket y\rrbracket=0,\quad\text{in }u\times\Xi,$$

$$\Vert y{\Vert }_{L^2(\Xi ,H^r(D))}^2:=\mathbb{E}[\Vert y{\Vert }_{H^r(D)}^r]<\infty .$$

$$a_{\xi }(y,p)=b_{\xi }(p),\forall p\in H_{\text{av}}^1(D)$$

$$a_{\xi }(y,p)$$

$$b_{\xi }(p)$$

$${\mathcal{L}}_{\xi }(u,y,p):=J(u,\xi )+a_{\xi }(y,p)-b_{\xi }(p),$$

$$\frac{\partial {\mathcal{L}}_{\xi }(u,y,p)}{\partial y}=\frac{\partial {\mathcal{L}}_{\xi }(u,y,p)}{\partial p}=0,$$

$$-\nabla \cdot (\kappa \nabla p)$$

$$\kappa \frac{\partial p}{\partial n}$$

$$\left\llbracket\kappa\frac{\partial p}{\partial n}\right\rrbracket=0,\quad\llbracket p\rrbracket=0,\quad\text{in }u\times\Xi.$$

$$d{\mathcal{L}}_{\xi }(u,y,p)[V]:=\underset{t\to 0^{+}}{\lim }\frac{{\mathcal{L}}_{\xi }(F_t(u),y,p)-{\mathcal{L}}_{\xi }(u,y,p)}{t},$$

$$\begin{array}{rl}dJ(u,\xi )[V]=& {\int }_D-\kappa \nabla y^T(\nabla V+\nabla V^T)\nabla p-(y-\bar{y})\nabla {\bar{y}}^{T}V-\nabla f^{T}Vp\\ & +\text{div}(V)\left(\frac{1}{2}(y-\bar{y}{)}^2+\kappa \nabla y^T\nabla p-fp\right)\text{d}x.\end{array}$$

$$dJ(u,\xi )[V]$$

$$\underset{x}{\min }\mathbb{E}[J(x,\xi )].$$

$$t^n\ge 0,\sum _{n=0}^{\infty }{t}^n=\infty ,\sum _{n=0}^{\infty }(t^n{)}^2<\infty ,$$

$$J(x^n-t^n\nabla J(x^n,{\xi }^n),{\xi }^n)\le J(x^n,{\xi }^n)-t^{n}c\Vert \nabla J(x^n,{\xi }^n){\Vert }^2$$

$$G^S(v,u)=dJ(u,\xi )[U]=a(V,U)\forall U\in H_0^1(D,{\mathbb{R}}^2),$$

$$a^{\text{elas}}(V,U):={\int }_D(\lambda \text{tr}(ϵ(V))\mathrm{id}+2\mu ϵ(V)):ϵ(U)\text{d}x,$$

$$a^{\text{elas}}(V,U)=dJ(u,\xi )[U]\forall U\in H_0^1(D,{\mathbb{R}}^2).$$

$$W_{h}\coloneqq\left\{y\in H^{1}(D)\,\big{|}\,y|_{T}\in\mathcal{P}_{1}\text{ for all }T\text{ in }D_{h}\right\}.$$

$$W_{\text{av},h}\coloneqq\left\{y\in H_{\text{av}}^{1}(D)\,\big{|}\,y|_{T}\in\mathcal{P}_{1}\text{ for all }T\text{ in }D_{h}\right\}$$

$$W_{0,h}\coloneqq\left\{V\in H_{0}^{1}({D,\mathbb{R}^{2}})\,\big{|}\,V_{i}|_{T}\in\mathcal{P}_{1}\text{ for all }i=1,2\text{ and }T\text{ in }D_{h}\right\}.$$

$$\text{Find }{y}_h\in W_{\text{av},h}:{\int }_{D_h}{\kappa }_h\nabla y_h^T\nabla p_h\text{d}x={\int }_{D_h}{f}_h{p}_h\text{d}x+{\int }_{\partial D_h}{g}_h{p}_h\text{d}s\forall p_h\in W_{\text{av},h},$$

$$\text{Find }{p}_h\in W_{\text{av},h}:{\int }_{D_h}{\kappa }_h\nabla p_h^T\nabla v_h\text{d}x={\int }_{D_h}(y_h-{\bar{y}}_h)v_h\text{d}x,\forall v_h\in W_{\text{av},h}$$

Full text

Key words. Computational methods, stochastic approximation, shape optimization, stochastic gradient method, interface identification, PDE constrained optimization under uncertainty.

AMS subject classifications. 60H35, 49Q10, 35R60, 60H15, 35R15

1 Introduction

Shape optimization is concerned with problems where an objective function is supposed to be minimized with respect to a shape, or a subset of $\mathbb{R}^{d}$. Many relevant problems contain additional constraints in the form of a partial differential equation (PDE), which describe the physical laws that the shape should obey. In applications, the material coefficients and external inputs might not be known exactly, but rather be randomly distributed according to a probability distribution obtained empirically. In this case, one might still wish to optimize over a set of these possibilities to obtain a more robust design. The random parameters are modeled with random fields, for example with a Karhunen-Loève expansion [13, 14] or polynomial chaos [30]. The PDE constraint is then replaced by a set of PDE constraints, or a random PDE constraint, accounting for all possible realizations of the stochastic space.

One challenge in shape optimization is finding the correct model to describe the set of shapes; another is finding a way to handle the lack of vector structure of the shape space. In principle, a finite dimensional optimization problem can be obtained for example by representing shapes as splines. However, this representation limits the admissible set of shapes, and the connection of shape calculus with infinite dimensional spaces [6, 26] leads to a more flexible approach. Recently, it was suggested to embed shape optimization problems in the framework of optimization on shape spaces [21, 29]. One possible approach is to cast the sets of shapes in a Riemannian viewpoint, where each shape is a point on an abstract manifold equipped with a notion of distances between shapes; see, e.g., [15, 16]. From a theoretical and computational point of view it is attractive to optimize in Riemannian shape manifolds because algorithmic ideas from [1] can be combined with approaches from differential geometry. Here, the Riemannian shape gradient can be used to solve such shape optimization problems using the gradient descent method. In the past, major effort in shape calculus has been devoted towards expressions for shape derivatives in the so-called Hadamard form, which are integrals over the surface (cf. [6, 26]). During the calculation of these expressions, volume shape derivative terms arise as an intermediate result. In general, additional regularity assumptions are necessary in order to transform the volume into surface forms. Besides saving analytical effort, this makes volume expressions preferable to Hadamard forms. In this paper, we consider the Steklov-Poincaré metric, which allows to use the volume formulations (cf. [22]).

Additional challenges arise in the stochastic setting. When the number of possible scenarios in the probability space is small, then the optimization problem can be solved over the entire set of scenarios. This approach is not relevant for most applications, as it becomes intractable if the random field has more than a few scenarios. For problems with random PDEs, either the stochastic space is discretized, or sampling methods are used. If the stochastic space is discretized, one typically relies on a finite-dimension assumption, where a truncated expansion is used as an approximation of the infinite-dimensional random field. Numerical methods include stochastic Galerkin method [3] and sparse-tensor discretization [23]. Sample-based approaches involve taking random or carefully chosen realizations of the input parameters; this includes Monte Carlo or quasi Monte Carlo methods and stochastic collocation [2]. In the stochastic approximation approach, dating back to a paper by Robbins and Monro [18], one uses a stochastic gradient in place of a gradient to iteratively minimize the expected value over a random function. Recently, stochastic approximation methods have been proposed to efficiently solve PDE-constrained optimization problems involving uncertainty [7, 11]. In this paper, we demonstrate a novel use of the stochastic gradient method on a model shape optimization problem.

The paper is structured as follows. In Section 2, a model problem will be presented; the model is an interface identification problem with uncertainty arising from PDE constraints containing random parameters and inputs. The problem will be formulated as a minimization problem over the expected deviation from a target measurement. The proposed algorithm for the numerical solution to the problem will be reviewed in Section 3. Numerical experiments demonstrating the effectiveness of the method are shown in Section 4. Finally, closing remarks are shared in Section 5.

2 Model formulation

We consider a model interface identification problem, which has been studied in the deterministic setting in a number of texts [5, 12, 26], and which we modify to allow for stochastic data. We allow for a domain containing multiple shapes, leading to a multi-interface identification problem as illustrated in Figure 1. This model is for instance used in electrical impedance tomography, where the distribution of materials with different properties is examined based on measurements on the boundary.

Let $D\subset\mathbb{R}^{2}$ be a bounded domain with boundary $\partial D$. This domain is assumed to be Lipschitz and partitioned in an open subdomain $D_{0}\subset D$ and a finite number $N\in\mathbb{N}$ of disjoint open subdomains $D_{i}\subset D$ with variable boundaries $u_{i}:=\partial D_{i}$, $\sqcup_{i=0,\dots,N}D_{i}\bigsqcup\left(\sqcup_{i=1,\dots,N}u_{i}\right)=D$, where $\sqcup$ denotes the disjoint union. For $i=1,\dots,N$, the union of all domains $D_{i}$ is denoted by $D_{\text{int}}:=\sqcup_{i=1,\dots,N}D_{i}$ and the union of all boundaries $u_{i}$ is called the interface and denoted by $u:=\sqcup_{i=1,\dots,N}u_{i}$. Figure 1 illustrates this situation. In our setting, $D$ is meant to be composed of up to $N+1$ distinct materials. One should keep in mind that $D$ depends on $u$, i.e., $D=D(u)$. If $u$ changes, then the subdomains $D_{i}\subset D$ change in a natural manner.

For the model, we concentrate on one-dimensional smooth shapes. In [16], the set of all one-dimensional smooth shapes is characterized by

[TABLE]

Here, $\mathrm{Emb}(S^{1},\mathbb{R}^{2})$ denotes the set of all embeddings from the unit circle $S^{1}$ into $\mathbb{R}^{2}$, which contains all simple closed smooth curves in $\mathbb{R}^{2}$. Note that the boundary of the subdomain already characterize a shape. Thus, we can think of one-dimensional smooth shapes as images of simple closed smooth curves in the plane of the unit circle. In (1), $\mathrm{Diff}(S^{1})$ is the set of all diffeomorphisms from $S^{1}$ into itself, which characterize all smooth reparametrizations. These equivalence classes are considered because we are only interested in the shape itself and images are not changed by reparametrizations. More precisely, shapes that have been translated represent the same shape. In contrast, shapes with different scaling are not equivalent in this shape space. In the following, we assume $u_{i}\in B_{e}$ for all $i=1,\dots,N$.

In this model, we allow different random variables accounting for source terms and material coefficients subject to uncertainty. It is assumed that probability distributions are known, for example by priorly obtained empirical samples. We allow for uncertainty in source terms and material constants by definition of a probability space $(\Omega,\mathcal{F},\mathbb{P})$, where $\mathcal{F}\subset 2^{\Omega}$ is the $\sigma$-algebra of events and $\mathbb{P}\colon\Omega\rightarrow[0,1]$ is a probability measure. In this paper, the probability space is to be understood as a product space $(\Omega,\mathcal{F},\mathbb{P})=(\Omega_{f}\times\Omega_{g}\times\Omega_{\kappa},\mathcal{F}_{f}\times\mathcal{F}_{g}\times\mathcal{F}_{\kappa},\mathbb{P}_{f}\times\mathbb{P}_{g}\times\mathbb{P}_{\kappa})$ as in [10]. We assume without loss of generality that this product space is complete. We define $f\colon D\times\Omega_{f}\rightarrow\mathbb{R}$ and $g\colon\partial D\times\Omega_{g}\rightarrow\mathbb{R}$ as a volume input and a boundary input function, respectively. We assume that these functions are finite dimensional noise in that there exist random vectors $\xi_{f}(\omega)=(\xi_{f}^{1}(\omega),\dots,\xi_{f}^{m_{1}}(\omega))$ with $\xi_{f}^{i}:\Omega\rightarrow\Xi_{f}^{i}\subset\mathbb{R}$ and $\xi_{g}(\omega)=(\xi_{g}^{1}(\omega),\dots,\xi_{g}^{m_{2}}(\omega))$ with $\xi_{g}^{i}:\Omega\rightarrow\Xi_{g}^{i}\subset\mathbb{R}$ such that

[TABLE]

On each subdomain $D_{i}$, the material is assumed to be constant, but randomly distributed. The function

[TABLE]

defines the random material coefficient over the domain $D$, where $\kappa_{i}$ are independent real-valued random variables and $\mathbbm{1}_{D_{i}}$ denotes the indicator function of the set $D_{i}$. To simplify notation, we set $\xi:=(\xi_{\kappa},\xi_{f},\xi_{g})$, $\Xi:=\Xi_{\kappa}^{\text{in}}\times\Xi_{\kappa}^{\text{out}}\times\Xi_{f}^{1}\times\cdots\times\Xi_{f}^{m_{1}}\times\Xi_{g}^{1}\times\cdots\times\Xi_{g}^{m_{2}}$ and now write $\kappa(\xi)=\kappa(\cdot,\xi)$, $f(\xi)=f(\cdot,\xi)$, and $g(\xi)=g(\cdot,\xi)$ for a given $\xi\in\Xi$.

For $J:B_{e}^{N}\times\Xi\rightarrow\mathbb{R}$, we use the notation $J(\cdot,\xi)\colon B_{e}^{N}\rightarrow\mathbb{R}$ to denote a deterministic functional for a single realization $\xi\in\Xi$. Here, $B_{e}^{N}$ stands for the cross-product of the shape manifold $B_{e}$ with itself $N$ times, one for each $u_{i}$. The expectation is then defined as the integral $\mathbb{E}[J(u,\xi)]:=\int_{\Omega}J(u,\xi(\omega))\,\textup{d}\mathbb{P}(\omega).$ For the tracking-type objective functional

[TABLE]

we consider the following problem:

[TABLE]

[TABLE]

where $D=D(u)$ depends on $u$ and $\bar{y}\colon D\to\mathbb{R}$ denotes (deterministic) observed measurements. With the tracking-type objective functional $J(\cdot,\xi)$ the model is fitted to the data measurements $\bar{y}$. The following continuity conditions are imposed explicitly for the state and flux at the interface:

[TABLE]

where the jump symbol \left\llbracket\cdot\right\rrbracket is defined on $u$ by $\llbracket v\rrbracket=v_{1}-v_{2}$ with $v_{1}=\text{tr}_{1}(v)$ and $v_{2}=\text{tr}_{2}(v)$ and $\text{tr}_{1}$, $\text{tr}_{2}$ are the trace operators defined on $u$ and $D_{0}$.

Remark 1

A regularization term such as perimeter regularization can be added to the objective functional (3). This is often necessary for analytical investigations, e.g., in order to establish that a problem is well-defined or to guarantee the existence of unique solutions. In simulations, we observed that the scaling of the perimeter penalization needed to be set so small as to be negligible, so for simplicity we leave off this term.

In order to solve model problem (4)–(7), we need to calculate its shape derivative. We use the standard notation for the Sobolev space $H^{1}(D)$, where $H_{0}^{1}(D)$ indicates the subspace of $H^{1}(D)$ containing functions with disappearing trace, and its vector valued versions $H^{1}(D,\mathbb{R}^{d}),\,H_{0}^{1}(D,\mathbb{R}^{d})$. For $H^{r}(D)$-valued random variables $y\colon H^{r}(D)\times\Xi\rightarrow\mathbb{R}$, we recall $y\in L^{2}(\Xi,H^{r}(D))$ if

[TABLE]

In addition, we define the Hilbert space $H_{\text{av}}^{1}(D):=\{v\in H^{1}(D)|\int_{D}v\,\textup{d}x=0\}$. In the following, $\xi\in\Xi$ is fixed but arbitrary. The weak formulation of the boundary value problem (5)–(7) for a fixed realization $\xi\in\Xi$ is: find $y=y(\cdot,\xi)\in H_{\text{av}}^{1}(D)$ such that

[TABLE]

with

[TABLE]

where $\kappa$ is defined as in (2). In order for the shape derivative to be well-defined, we make the following technical assumptions: $\bar{y}\in H^{1}(D)$, $f\in L^{2}(\Xi,H^{1}(D))$ and $g\in L^{2}(\Xi,L^{2}(\partial D))$.

The Lagrangian of (4)–(7) is defined as

[TABLE]

where $J(\cdot,\xi)$ is defined in (3), $a_{\xi}$ in (9) and $b_{\xi}$ in (10). For any $\xi\in\Xi$, a saddle point $(y,p)\in H_{\text{av}}^{1}(D)\times H_{\text{av}}^{1}(D)$ of the Lagrangian is given by

[TABLE]

which leads to the state equation (5)–(7) and the adjoint equation

[TABLE]

for $p\colon D\times\Xi\to\mathbb{R},\,(x,\xi)\mapsto p(x,\xi)$. As with the state equation, we have interface conditions for the adjoint equation:

[TABLE]

The design equation is given by the shape derivative. If we consider the perturbation of identity, the shape derivative of ${\mathcal{L}}$ at $u$ and fixed $\xi\in\Xi$ in the direction of a sufficiently smooth vector field $V$ is defined by

[TABLE]

where $F_{t}(u):=\{\text{id}(x)+tV(x)\colon x\in u\}$ and $F_{0}(u)=u$. Using standard techniques for calculating the shape derivative (min-max approach [6], chain rule approach [27], rearrangement method [12]), we get the shape derivative of our model problem in volume formulation:

[TABLE]

Remark 2

Note that the shape derivative arises in two equivalent notational forms:

[TABLE]

Here, $R$ is a differential operator acting linearly on the vector field $V$ and $r\in L^{1}(u)$. In this paper, we use the volume formulation.

Remark 3

In this paper, we consider the perturbation of identity. Of course, it is also possible to use, e.g., the speed method in order to transform the domain $D$. In this case, the definition of $F_{t}$ changes.

Remark 4

As mentioned in the introduction, the algorithm proposed in this work is basically a generalization of the stochastic gradient method. For this method, descent directions can be chosen as the gradient of $J(\cdot,\xi)$ for a randomly chosen $\xi$. This will be clarified in the next section. This is the reason why the shape derivative is computed just for $J(\cdot,\xi)$ and not for $j$.

Now that we have the state, adjoint and design equations, we are ready to solve our problem computationally. In the next section, we will present our proposed method for the computational solution of the problem, including the discretization scheme, optimization method, and post-processing techniques for the improvement of mesh aspect ratios.

3 Algorithmic and computational details

This section is split into two subsections. Section 3.1 is devoted to the presentation of the algorithmic aspects for the solution of the problem described in section 2. To our knowledge, this is a novel application of the stochastic gradient method to a shape optimization problem under uncertainty. In section 3.2, we focus on computational aspects and combine these with the algorithmic details of section 3.1 resulting in a stochastic shape gradient method for a discretized version of the model problem.

3.1 Algorithmic aspects

The stochastic gradient method is an algorithm for solving problems of the form

[TABLE]

The method uses iterations of the form $x^{n+1}=x^{n}-t^{n}\nabla J(x^{n},\xi^{n})$, where at iteration $n$ the stochastic gradient at $x^{n}$ for a random sample ${\xi^{n}}\in\Omega$, is given by $\nabla J(x^{n},\xi^{n})\approx\nabla\mathbb{E}[J(x^{n},\xi)]$. We make the assumption that it is possible to generate independent and identically distributed (i.i.d.) samples $\xi^{n}$ from a known probability distribution. This method dates back to [18] and has been used in many different applications, notably as well in optimization on manifolds [4] and in PDE-constrained optimization [7]. Much research has been devoted to the proper choice of the step size $t^{n}$. The least involved choice is the so-called Robbins-Monro step size,

[TABLE]

which was also used in the original work [18] to establish almost sure convergence of the sequence $\{x^{n}\}$ to a stationary point for (18). It is notable that the Armijo rule $t^{n}:=\alpha\rho^{m_{n}}$ with minimal $m_{n}\in\mathbb{N}_{0}$ satisfying

[TABLE]

with $\alpha>0$, $\rho\in(0,1)$, $c\in(0,1)$ may not converge without additional variance reduction techniques. If instead of choosing a single realization at iteration $n$, one chooses $N_{n}$ samples $\xi^{n}=(\xi^{n,1},\dots,\xi^{n,N_{n}})$ with $N^{n}\rightarrow\infty$ as $n\rightarrow\infty$, then convergence using (20) can be demonstrated under relatively weak assumptions [28, 24]. Note that this involves increasingly expensive calculations for the estimate of the objective function $J(x,\xi^{n})=1/N_{n}\sum_{l=1}^{N_{n}}J(x,\xi^{n,l})$ and likewise for the stochastic gradient $\nabla J(x,\xi^{n}).$ In numerical simulations, we will compare both of these approaches. Further step size rules can be reviewed in a nice summary by George and Powell [9].

Now, we formalize what is meant by stochastic gradient for our problem. In this work, we focus on one-dimensional smooth shapes which are elements of the shape space $B_{e}$ defined in (1). This shape space is a Riemannian manifold (cf. [16, 17]). If we want to optimize on a Riemannian shape manifold, we have to find a representation of the shape derivative with respect to the Riemannian metric under consideration, called the Riemannian shape gradient, which is required to formulate optimization methods on a shape manifold. In [22], the authors present a metric based on the Steklov-Poincaré operator, which allows for the computation of the Riemannian shape gradient as a representative of the shape derivative in volume form. Besides saving analytical effort during the calculation process of the shape derivative, this technique is computationally more efficient than using an approach which needs the surface shape derivative form. For example, the volume form allows us to optimize directly over the hold-all domain $D$ containing one or more elements $u_{i}\in B_{e}$, whereas the surface formulation would give us descent directions (in normal directions) for the boundary $u_{i}$ only, which would not help us to move mesh elements around the shape. Additionally, when we are working with a surface shape derivative, we need to solve another PDE in order to get a mesh deformation in the hold-all domain $D$. Following the ideas presented in [22], we choose the Steklov-Poincaré metric, denoted by $G^{S}$. In the setting of the shape space $B_{e}$, the mesh deformation vector $V\in H_{0}^{1}(D,\mathbb{R}^{2})$ can be viewed as an extension of a Riemannian shape gradient to the hold-all domain $D$ because of the identities

[TABLE]

where $v=(\text{tr}(V))^{T}n,u=(\text{tr}(U))^{T}n$ with $\text{tr}(\cdot)$ denotes again the trace operator on the Sobolev spaces. Here, $a(\cdot,\cdot):H_{0}^{1}(D,\mathbb{R}^{2})\times H_{0}^{1}(D,\mathbb{R}^{2})\rightarrow\mathbb{R}$ is a symmetric and coercive bilinear form. One option for the operator $a(\cdot,\cdot)$ is chosen to be the bilinear form associated with the linear elasticity problem, i.e.,

[TABLE]

where $\epsilon(U)\coloneqq\frac{1}{2}\,(\nabla U+\nabla U^{T})$, $A:B$ denotes the Frobenius inner product for two matrices $A,B$ and $\lambda,\mu$ denote the Lamé parameters. To summarize, for a fixed $\xi\in\Xi$, we need to solve the following so-called deformation equation: find $V\in H_{0}^{1}(D,\mathbb{R}^{2})$ such that

[TABLE]

The main advantage of this Steklov-Poincaré metric approach is that the identity (21) holds, meaning that the Riemannian metric $G^{S}(\cdot,\cdot)$, which is naturally defined over the interfaces, can be equivalently reformulated in terms of the bilinear form $a(\cdot,\cdot)$ over the whole domain. This last observation is the main approach we will use in the numerical solution of our model problem.

Remark 5

In general, we need the concept of the exponential map and vector transports in order to formulate optimization methods on a shape manifold. The calculations of optimization methods have to be performed in tangent spaces because manifolds are not necessarily linear spaces. This means points from a tangent space have to be mapped to the manifold in order to get a new shape-iterate, which can be performed with the help of the exponential map. However, the computation of the exponential map is prohibitively expensive in the most applications because a calculus of variations problem must be solved or the Christoffel symbols need to be known. It is much easier and much faster to use a first-order approximation of the exponential map. In [1], it is shown that a so-called retraction is such a first-order approximation and sufficient in most applications. We refer to [20], where a suitable retraction on $B_{e}$ is given.

3.2 Computational aspects

In order to numerically solve problem (4)–(7), we need to consider some kind of discretization; we use the finite element method (FEM). Let $D_{h}$ be a triangulation of the domain $D$ with a partition of subsets $D_{i,h}:=\overline{D_{i}}\cap D_{h}$ such that $D_{h}=\sqcup_{i=0,\dots,N}D_{i,h}$. Moreover, let $u_{i,h}:={D_{i,h}}\cap{D_{0,h}}$ denote the discretization of the interfaces for $i=1,\dots,N$ and $\partial D_{h}:=\partial D\cap D_{h}$. We assign $D_{\text{int},h}:=\sqcup_{i=1,\dots,N}D_{i,h}$ and $u_{h}:=\sqcup_{i=1,\dots,N}u_{i,h}.$ We consider the set,

[TABLE]

where $\mathcal{P}_{1}$ is the space of all first degree polynomials. Furthermore, let us define the spaces

[TABLE]

where $H_{\text{av}}^{1}(D):=\{v\in H^{1}(D)|\int_{D}v\,\textup{d}x=0\}$, and

[TABLE]

Once we have discretized the domain and found an approximation for the Sobolev spaces, we are able to present the discrete version of the state equation, which is given in its weak formulation as follows:

[TABLE]

where $\mathbb{I}_{h}:C(\bar{D})\rightarrow W_{h}$ is the piecewise linear interpolation of a continuous function in the space $W_{h}$ and $\kappa_{h}:=\mathbb{I}_{h}(\kappa)$, $f_{h}:=\mathbb{I}_{h}(f)$, $g_{h}:=\mathbb{I}_{h}(g)$. In the same way, the adjoint state will be computed by solving the problem

[TABLE]

where $\bar{y}_{h}:=\mathbb{I}_{h}(\bar{y})$. Furthermore, as explained in the previous section, these equations need to be solved for a random $\xi$. For readability, we have suppressed the dependence of $\kappa_{h}$, $y_{h}$, and $p_{h}$ on $\xi$.

Once we have computed the state and adjoint variables, the shape derivative will be given by the discrete counterpart of equation (17), i.e.,

[TABLE]

Finally, the discretized version of the deformation equation is given by: find $V_{h}\in W_{0,h}$ such that

[TABLE]

As mentioned in Remark 5, after the computation of the deformation field, we need to use the exponential function or an approximation of it (like a retraction mapping) in order to update the shape for the next iteration. In this context, following the ideas of [20], we have chosen a retraction that in its discretized version acts like the perturbation of identity method. That means we assign a direction given by $V_{h}$ to each of the nodes of the triangulation $D_{h}$ and move them with certain step size. In order to guarantee good behavior of the method (invertibility of the operator $F_{t}$ from the perturbation by identity), this step sizes need to be small; otherwise, the algorithm could produce low aspect ratios or overlapping elements. One of the principal reasons for this behavior is that the perturbation by identity method does not consider mesh connectivity, or how close a node is to its neighbors. In order to overcome this drawback, we use the following techniques.

As discussed in [22], an unmodified right-hand side of the discretized deformation equation leads to deformation fields causing meshes with bad aspect ratios. Thus, we set values of the discretized shape derivative to zero if the corresponding element does not intersect with the interface, i.e.,

[TABLE]

Furthermore, in order to ensure good behavior of the mesh along iterations of the algorithm, we follow ideas from [19]; for each iteration step $n$, we choose the Lamé parameters as follows: $\lambda=0$ and $\mu\in[\mu_{\min},\mu_{\max}]$ decreasing smoothly from $u$ to the outer boundary. One possible way to model this behavior is to solve the following Poisson equation

[TABLE]

In section 4.2, we will present an experiment that shows the importance of the choice for $\mu_{\min}$ and $\mu_{\max}$. Correct choices for these parameters ensure that the cells in the mesh maintain good aspect ratios throughout the optimization process.

In algorithm 1, we present the Stochastic Shape Gradient Method (SSGM) formulated specifically for the discrete version of the model problem described in section 2.

The model problem (4)-(6) is formulated for one-dimensional shapes $u_{i}\in B_{e}$, $i=1,\dots,N$. This means that we are dealing with objects in 2D. Of course, it is possible to generalize our problem to 3D shapes. For a generalization, we refer to to [25], which considers a multiple shape interface problem in 3D. One focus of this paper is the performance of 3D computations. Here, the effectiveness of the computations is guaranteed by using multigrid methods and parallelization techniques. We want to mention that our proposed method can also be applied to higher-dimensional objects. However, one needs to take into account that problems occur, which are due to the high dimension. In 3D, we have to deal with many degrees of freedom. Thus, our method needs some tuning in higher dimensions to solve the problem efficiently. For example, parallelization techniques or multigrid methods can be performed, as suggested by [25]. Since addressing this kind of question goes beyond the scope of this paper, we restrict ourselves to 2D. Although the stochastic gradient, the numerical solution of PDEs, constrained shape optimization problems had been widely studied in the past, there have not been a combination of both, and the scientific computing is challenging as the experiments in the next section show. We want to highlight that the main advantage of our method is the lower computational cost added by considering random parameters. This fact is obtained due to the low complexity of the stochastic gradient. Furthermore, up-to-the authors’ knowledge, there has not been an extensive study on the different techniques for preventing mesh destruction and bad aspect ratios like the ones presented here.

Remark 6

We want to mention that model problem (4)-(6) for $N=1$ in higher dimensions is investigated analytically in [8]. In particular, an asymptotic convergence result on Riemannian manifolds is proven.

4 Numerical experiments

The principal aim of this section is to show the behavior of algorithm 1 through numerical experiments. This section is divided into five experiments. The first one shows the applicability to a problem with multiple shapes. The second one presents an analysis of the behavior of the meshes for different values of the parameters $\mu_{\min}$ and $\mu_{\max}$ in the equation (28). The third shows convergence to different target shapes. The fourth shows the influence of the individual random variables on the outcome. Finally, we look at results obtained using different step size rules and probability distributions with high variance.

Throughout this section we choose $D=(0,1)^{2}$ and discretize as described in section 3.2 with up to 10,000 triangles. In all experiments, the volume input $f$ is set to $f\equiv 0.$ We use $\mathcal{N}(\rho,\sigma,a,b)$ to denote a truncated normal distribution with mean $\rho$, standard deviation $\sigma$, and lower and upper bounds denoted by $a$, $b$, respectively. We consider a simplified version of the model where $\kappa_{\text{int}}:=\kappa_{1}=\cdots=\kappa_{N}.$ The measurement $\bar{y}$ is generated by solving the state equation (24) on a target shape $\bar{u}$ with constant parameters $\kappa_{0}=1.5$, $\kappa_{\text{int}}=4$ and $g=10$. Unless otherwise indicated, the bounds of the Lamé parameter $\mu$ were set to $\mu_{\min}=10$ and $\mu_{\max}=25$.

In all experiments, we use a single realization $\xi^{n}$ at each iteration $n$ to generate the stochastic shape gradient. For some convergence plots, we approximate function values $\mathbb{E}[J(u^{n},\xi)]$ and the $L^{2}$-norm of the shape gradient $\mathbb{E}[\lVert V^{n}\rVert_{L^{2}(D,\mathbb{R}^{2})}]$ using additional sampling at each point $u^{n}$ by the estimates

[TABLE]

with $m$ i.i.d. samples $\{\xi^{n,1},\dots,\xi^{n,m}\}$ generated at each $n$.

4.1 Multiple shapes

The main objective of this experiment is to present the numerical solution of problem (4)- (7) for the target configuration depicted in Figure 2. Additionally, we present experiments on meshes with different sizes to observe their effect on the solution. We used the Armijo line search rule (20) with $\alpha=50$, $\rho=0.5$, and $c=10^{-4}$ and the random parameters $\kappa_{0}\sim\mathcal{N}(1.5,10^{-2},1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$ and $g\sim\mathcal{N}(10,10^{-2},9,11)$.

First, we used a mesh with approximately 3,000 elements. We let the algorithm iterate 300 times, and Figure 3 shows the initial, an intermediate and the final shapes obtained. A convergence plot for a single sample of the objective function value and stochastic shape gradient is shown in Figure 5(a). Furthermore, we use $m=100$ samples to approximate the objective function at the last iteration to find $\hat{j}_{300}\approx 0.00323.$

Second, we would like to show that the solutions can be improved just by taking a finer mesh for both the initial and target shapes, in this case with around 10,000 elements. In Figure 4, we present the obtained solutions. Moreover, the convergence plot for this experiment is showed in Figure 5(b), and the value of the objective function is approximated with $m=100$ samples with $\hat{j}_{300}\approx 0.00311$. As expected, we observe a lower objective function value with the finer mesh.

4.2 Influence of the Lamé parameters on the mesh quality

For simplicity, from now on we only will consider shapes as depicted in Figure 6(a); i.e., consisting only of one interior interface. The main goal of this experiment is to show the importance of the parameters $\mu_{\min}$ and $\mu_{\max}$ in achieving better aspect ratios of the final meshes. The initial shape was chosen to be an ellipse as depicted in Figure 6(b) on a mesh with about $10,000$ triangles. We used the Armijo line search rule (20) with $\alpha=400$, $\rho=0.5$, $c=10^{-4}$ and $\kappa_{1}\sim\mathcal{N}(1.5,10^{-2},1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$ and $g\sim\mathcal{N}(10,10^{-2},9,11)$. We let the algorithm iterate 200 times.

Figures 7 and 9 show that of the tested values, the aspect ratio for the final mesh is the best in the case where we choose $\mu_{\min}=10,\mu_{\max}=25$. These results suggest that the correct choice of these values is problem-dependent and would be an opportunity for future research. Specifically, we have improved the aspect ratios of the cells near to the moving boundary by increasing the value of $\mu_{\min}$ and the non-moving boundary by increasing the value of $\mu_{\max}$.

4.3 Different target shapes

In this experiment, we show convergence to different target shapes. We used the Armijo line search rule (20) with $\alpha=400$, $\rho=0.5$, $c=10^{-4}$ and $\kappa_{0}\sim\mathcal{N}(1.5,10^{-2},1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$ and $g\sim\mathcal{N}(10,10^{-2},9,11)$. The measurement $\bar{y}$ was generated for two different target shapes $\bar{u}$: the ellipse (Figure 11(a)) and a small circle (Figure 10(a)). The starting shape $D^{0}$, final shape $D^{200}$ and convergence plot is shown in each figure. In the convergence plots, one sees a clear decrease in the objective function value and the norm of the shape gradient as well as significant progress toward to the target shape.

4.4 Influence of the individual random variables

The main purpose of this experiment is to study the influence of the individual random variables on the performance of the algorithm. That means we will choose different values to be deterministic while other will remain random. The initial shape $D^{0}$ is depicted in Figure 12(a) and the target shape $\bar{u}$ is depicted in Figure 12(b).

The behavior of the objective function and the norm of the shape gradient are depicted in Figure 13, where we choose the following parameters.

Figure 13(a):

$\kappa_{0}\sim\mathcal{N}(1.5,10^{-2},1,2)$, while $\kappa_{\text{int}}=4$ and $g=10$.

Figure 13(b):

$\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$, while $\kappa_{0}=1.5$ and $g=10$.

Figure 13(c):

$g\sim\mathcal{N}(10,10^{-2},9,11)$, while $\kappa_{0}=1.5$ and $\kappa_{\text{int}}=4$.

Figure 13(d):

$\kappa_{0}\sim\mathcal{N}(1.5,10^{-2},1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$, and $g=10$.

Figure 13(e):

$\kappa_{0}\sim\mathcal{N}(1.5,10^{-2},1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,10^{-2},3,5)$ and $g\sim\mathcal{N}(10,10^{-2},9,11)$.

Figure 13(f):

$\kappa_{0}=1.5$, $\kappa_{\text{int}}=4$ and $g=10$.

By comparing Figures 13(b) and 13(f) we can conclude that the randomness in $\kappa_{\text{int}}$ influences the behavior of the algorithm only minimally. That agrees with the model, since this parameter affects only the moving boundary from the inside. Comparing Figures 13(a), 13(d), and 13(e), we observe that the stochastic shape gradient appears to be most sensitive to the parameter $\kappa_{0}$, which induces greater oscillations in the objective function throughout optimization. We can also see in Figure 13(c) that the stochastic shape gradient and objective function value do not appear to be very sensitive to randomness in the boundary condition $g$.

4.5 Step size rules for high variance distributions

In this experiment, we will show the influence of the step size $t^{n}$ in Algorithm 1 on convergence in the high-variance case. We will compare the Robbins-Monro step size rule (19) with the Armijo line search rule (20) as well as a modification of the latter. To demonstrate the effect of high variance on the problem, we will work with the following random variables with higher variance than in the past examples: $\kappa_{0}\sim\mathcal{N}(1.5,0.2,1,2)$, $\kappa_{\text{int}}\sim\mathcal{N}(4,0.2,3,5)$ and $g\sim\mathcal{N}(10,0.2,9,11)$. The initial shape is shown in Figure 14(a) and the target measurement $\bar{y}$ is generated on the shape in Figure 14(b). Convergence plots are shown with $m=1000$ samples for each iteration $n$.

4.5.1 Robbins-Monro step size

We work with the rule $t^{n}=\alpha n^{-0.85}$, which obviously satisfies the conditions in (19). In Figures 15 and 16, we compare results of a simulation using $\alpha=800$ and $\alpha=400$, respectively. This demonstrates that a poorly chosen constant can drastically affect convergence speed. Interestingly, at $n=200$, the better choice $\alpha=800$ has $\hat{j}_{200}\approx 0.371,$ while $\alpha=400$ has $\hat{j}_{200}\approx 0.368$, so a slightly lower objective function value. However, as one clearly sees from the plots, even with $m=1000$ samples, the estimated function values are quite noisy; therefore, we cannot conclude that the shape obtained in Figure 16 has a lower objective function value than the one in Figure 15.

4.5.2 Armijo line search

We use the Armijo line search rule (20) with $\alpha=300$, $\rho=0.5$, $c=10^{-4}$ and $N_{n}\equiv 1$ sample per iteration. In Figure 17, one sees the final shape as well as a plot of the objective function. In this simulation, the Armijo line search was successful in the first iteration for every $n$, effectively making the step size constant $t^{n}\equiv 150$. Unlike the Robbins-Monro rule, this step size rule does not reduce variance well, since large steps are taken even if one has already reached the neighborhood of the optimum. This induces large oscillations in the shape, which may compromise mesh integrity. The objective function value also oscillates heavily as a result, as one can see in the plot of the objective function. At $n=200$, $\hat{j}_{200}\approx 0.360$, where it is noted that one would need to take more samples before making a fair comparison with the previous example.

4.5.3 Armijo line search with damping

To mimic the variance reduction property of the Robbins-Monro rule, we introduce Armijo line search with damping to reduce oscillations around the optimum. This line search takes the form

[TABLE]

for a decreasing sequence $\alpha_{n}\in[\alpha_{\min},\alpha_{\max}].$ Here, we choose $\alpha_{0}=400,$ $\rho=0.5$, and $c=10^{-4}$, and damp $\alpha_{n}$ every $20^{th}$ iteration (i.e. $\alpha_{20}=0.9\alpha_{0}$, $\alpha_{21}=\cdots=\alpha_{39}=\alpha_{20}$, $\alpha_{40}=0.9\alpha_{20}$, etc.) Figure 18 shows the shape and a plot of objective function. Although the function values still oscillate, the extreme values are lower. At $n=200$, $\hat{j}_{200}\approx 0.361$. While this experiment shows some promise in the hybrid step size rule (29), it is to be noted that to the authors’ knowledge, no convergence theory for this method exists; this would be an opportunity for future research.

4.5.4 Concluding observations

Notably, all of these methods are quite sensitive to the scaling of the step size, requiring the user to tune parameters offline. This is in some sense unavoidable for iterative methods on nonconvex problems. For all step size rules, if the scaling is chosen to be too large, mesh destruction or bad aspect ratios of the mesh may appear as seen in Figures 19 and 20. For the Robbins-Monro rule, the danger of mesh destruction occurs at the beginning, when the step size is its largest. Figure 19 shows how the mesh has bad aspect ratios already at the second iteration if $t^{n}=2000\cdot{n}^{-0.85}$. Figure 20 shows a mesh with bad aspect ratios at the final iteration, where high oscillations around iteration $n=175$ caused a compromised mesh.

The Armijo rule, which is widely used for deterministic problems, should be used with care for stochastic problems. As mentioned in Section 3, current convergence theory requires that one takes more and more samples of the gradient to achieve variance reduction. The Robbins-Monro rule, which naturally reduces variance over the optimization, suffers from a lack of robustness: slight changes in the scaling of the step size lead to drastically different efficiency in the first stage of optimization, as we saw in the experiments. An idea is to try to reduce variance while still allowing for large step sizes, which we used to develop the damping procedure. This method, however, introduces an additional parameter needing to be tuned, namely the damping rate. So far, no theoretical results exist for such a procedure, but it would be an interesting topic to explore, particularly for nonconvex problems.

5 Conclusion

A stochastic shape interface problem was investigated computationally. The problem under consideration is of high importance because adding randomness makes numerical simulations for real-world problems with uncertainty in the model. In order to solve the problem, we combine the gradient descent method using the shape derivative in its volume expression with the well-known stochastic gradient, which results in a stochastic gradient descent method for shape optimization problems in shape spaces. We present various numerical experiments in this paper by applying this method to the shape interface model. Numerical solutions of the problem are presented, where the domain consists of multiple shapes. Moreover, the influence of the Lamé parameters on the mesh quality and the influence of individual random variables on the performance of the algorithm are investigated. In a further experiment, we perform computations for different target shapes, observing the convergence of a circle to an ellipse and vice versa. In a final experiment, which focuses on the various step size rules, we use an example with higher variance and compare simulations for the Robbins-Monro step size rule, the Armijo line search rule, and a hybrid method we call Armijo line search with damping. The strengths and weaknesses of the methods are compared.

Many interesting questions remain, some of which we will attempt to address in future work. For example, how to extend our method efficiently to problems with two-dimensional shapes, i.e., where the underlying domains are subsets of $\mathbb{R}^{3}$. Moreover, it is well done that with highly dimensional problems it arises the question of improving the performance of the algorithm. For example, parallelization techniques could be solution. Additionally, it would be useful to compute efficiency estimates so that a termination condition can be defined for the algorithm. As we demonstrated in the last experiment, step size plays a big role, especially when the variance is high; although no convergence theory exists for the Armijo with damping rule, it may a good choice in applications, as it appears to possess the strengths of the other two methods.

Acknowledgments

This work has been partly supported by the German Research Foundation (DFG) within the priority program SPP 1962 under contract number WE 6629/1-1 and Schu804/15-1, the Austrian Science Fund (FWF) within project no. W1260-N35, and the German Academic Exchange Service (DAAD) within the Doctoral Programm 2017/18. The authors are indebted to Laura Bittner for discussions about the model problem.