An a posteriori error estimator for shape optimization: application to EIT

In this paper we account for the numerical error introduced by the Finite Element approximation of the shape gradient to construct a guaranteed shape optimization method. We present a goal-oriented strategy inspired by the complementary energy principle to construct a constant-free, fully-computable a posteriori error estimator and to derive a certified upper bound of the error in the shape gradient. The resulting Adaptive Boundary Variation Algorithm (ABVA) is able to identify a genuine descent direction at each iteration and features a reliable stopping criterion for the optimization loop. Some preliminary numerical results for the inverse identification problem of Electrical Impedance Tomography are presented.


Introduction
A classical approach to shape optimization and shape identification problems relies on the application of the Steepest Descent Method using the so-called shape gradient. In most applications, the shape gradient depends on the solution of a PDE which can only be solved approximately. Thus stopping criteria based on its norm are never fulfilled if the a priori given tolerance is too small. In this paper, we evaluate the discretization error in the shape gradient relying on the well-established results on goal-oriented error estimators [1] and we propose a guaranteed shape optimization strategy. Several works in the literature have focused on Adaptive Finite Element Method for shape optimization problems ( [2], [3], [4], [5]) but to the best of our knowledge, all the proposed approaches used error estimators only to drive mesh adaptivity without accounting for the error in the shape gradient itself.
After briefly recalling the framework of shape optimization and shape identification problems in section 2 and the notion of shape gradient in section 3, we introduce the Adaptive Boundary Variation Algorithm for shape optimization problems (Section 4). In section 5 we present some preliminary numerical simulations for the inverse problem of Electrical Impedance Tomography and a summary of our results ends the paper.

Shape optimization and shape identification problems
Let us consider an open domain Ω ⊂ R d (d ≥ 2) with Lipschitz boundary ∂Ω. Let V Ω be a separable Hilbert space depending on Ω, we define u Ω ∈ V Ω to be the solution of a state equation where a Ω (·, ·) is a continuous and coercive bilinear form on V Ω and F Ω (·) is a continuous linear form on V Ω , both of them depending on Ω. Under these assumptions, problem (1) has a unique solution.
We introduce a cost functional J(Ω) = j(Ω, u Ω ) which depends on the domain Ω itself and on the solution u Ω of the state equation. We consider the following shape optimization problem where U ad is the set of admissible domains in R d . Within this framework, problem (2) may be viewed as a PDE-constrained optimization problem where we aim at minimizing j(Ω, u) under the constraint u = u Ω .

The case of Electrical Impedance Tomography
Let D be a domain featuring an inclusion Ω. We define the conductivity k Ω as a piecewise constant function equal to k I inside the inclusion Ω and k E otherwise (k I , k E > 0). We are interested in identifying the shape and the location of the inclusion, fitting given boundary measurements g and U D of the flux and the potential. Thus we introduce u Ω,N , u Ω,D ∈ H 1 (D) such that u Ω,D = U D on ∂D, solutions of respectively a Neumann and a Dirichlet Boundary Value Problem: The inverse identification problem of Electrical Impedance Tomography may be written as in (2) and we seek the open subset that minimizes the Kohn-Vogelius functional [6]: 3. Differentiation with respect to the shape Let X be a Banach space and θ ∈ X be an admissible smooth deformation of Ω. The objective functional J is said to be X-differentiable at Ω ∈ U ad if there exists a continuous linear form dJ(Ω) on X such that ∀θ ∈ X we have J((Id +θ)Ω) = J(Ω) + ⟨dJ(Ω), θ⟩ + o(θ). Let us introduce the Lagrangian functional, defined for every admissible open set Ω and every and the adjoint state Let us now introduce a diffeomorphism φ : For every admissible set Ω, and for all We admit that u → u φ is a one-to-one map between V Ω and V φ(Ω) . We define the material derivative ∂L ∂φ of L as the derivative of the map φ → L(Ω φ , u φ , p φ ) evaluated in φ = Id. Hence, from the fast derivation method of Céa [7] the shape gradient reads as The case of Electrical Impedance Tomography In order to compute the shape gradient of (4), we need to determine the adjoint solutions p Ω,N and p Ω,D associated with the states u Ω,N and u Ω,D : the Kohn-Vogelius problem is self-adjoint and we get that p Ω,N = u Ω,N − u Ω,D and p Ω,D = 0. The most common approach in the literature is based on an Eulerian viewpoint and leads to a surface expression of the shape gradient Let us now introduce the following operator: (7), we get a volumetric expression of the shape gradient We refer to [8] for more details on the differentiation of the Kohn-Vogelius functional and its application to the identification of discontinuities of the conductivity parameter.

The Adaptive Boundary Variation Algorithm for shape optimization
In algorithm 1, we present the Boundary Variation Algorithm, that is the application of the classical Steepest Descent Method to a shape optimization problem. After computing the solution of the state equation, we derive the expression of the shape gradient and we compute the corresponding descent direction to deform the computational domain.
This method relies on the computation at each iteration of a direction θ such that the shape gradient of the objective functional in this direction is strictly negative, that is we seek θ such that ⟨dJ(Ω), θ⟩ < 0. Let us now introduce the discretized solutions u h Ω and p h Ω of the state and adjoint equations using the Finite Element Method. We consider for both u h Ω and p h Ω the approximation in the space of linear Lagrangian Finite Element functions. In the same fashion, we define the Finite Element discretization of the shape gradient as The discretized descent direction θ h ∈ X is computed by solving problem (12) using linear Lagrangian Finite Elements.
(θ h , δθ h ) X + ⟨d h J(Ω), δθ h ⟩ = 0 ∀δθ h ∈ X (12) We remark that due to the numerical approximation, it is not sufficient that ⟨d h J(Ω), θ h ⟩ < 0 for J(Ω) to decrease along the direction θ h . Moreover, if the approximation error is too large the stopping criterion (Algorithm 1 -step 4) will usually not be fulfilled. In this paper we present a strategy to bypass this issue by providing a computable upper bound of the error E h in the approximation of the shape gradient and by seeking a certified direction θ h such that ⟨d h J(Ω), θ h ⟩ + |E h | < 0.

Numerical error in the shape gradient
Let us write the numerical error E h in the shape gradient as follows In order to compute an upper bound of the error (13), we use the well-established framework of a posteriori error estimators for a Quantity of Interest [1]. For this purpose we introduce two adjoint problems, each of which is associated with one term on the right-hand side of (13) and we seek r Ω , s Ω ∈ V Ω such that ∀δr, δs ∈ V Ω a * Ω (r Ω , δr) = By plugging (14) into (13) and exploiting Galerkin orthogonality and Cauchy-Schwarz inequality, we may derive the upper bound E for the numerical error in the shape gradient: where |||·||| Ω is the energy-norm induced by the bilinear form a Ω (·, ·). We remark that the criterion ⟨d h J(Ω), θ h ⟩ + E < 0 ensures that θ h is a certified genuine descent direction, whether it is the solution of equation (12) or not. In particular, this also stands when the latter problem is only solved approximately.

Certified Steepest Descent Method
We may now present the resulting Adaptive Boundary Variation Algorithm arising from the coupling of the classical Boundary Variation Algorithm for shape optimization with the goaloriented estimator of the error in the shape gradient.

The case of Electrical Impedance Tomography
Since the Kohn-Vogelius problem is self-adjoint, (15) reduces to where r Ω,N ∈ H 1 (D) and r Ω,D ∈ H 1 0 (D) are the solutions of the following adjoint problems introduced to construct the estimator in the Quantity of Interest: k Ω ∇r Ω,i · ∇δr i + r Ω,i δr i To fully compute the estimate (16), we use the complementary energy principle. For the state equations, we introduce the dual flux variables σ Ω,N , σ Ω,D ∈ H(div) = {σ ∈ L 2 (D; R d ) : div σ ∈ L 2 (D)} such that σ Ω,N · n = g on ∂D which satisfy ∀δσ N , δσ D ∈ H(div) such that δσ N · n = 0 on ∂D Let σ h Ω,N and σ h Ω,D be the dual fluxes discretized using the lowest order Raviart-Thomas Finite Elements. We obtain the following bound for the energy-norm of the error in the state equations: In a similar fashion, we introduce the dual fluxes ξ Ω,N , ξ Ω,D ∈ H(div) such that ξ Ω,N · n = 0 on and via their RT 0 approximations ξ h Ω,i 's, we derive the upper bound of the adjoint errors (21) Eventually, by plugging (19) and (21) into (16), we are able to explicitly compute the bound E.

Numerical results
We apply the ABVA to the identification of an unknown inclusion Ω inside a circular domain D using one boundary measurement. The deformation of the inclusion (Algorithm 2 -step 6) is performed by moving the whole computational mesh. In figure 1c, we present the reconstructed interface after 25 iterations of the procedure; figure 1d confirms that the ABVA identifies a genuine descent direction at each iteration, generating a sequence of minimizing shapes such that the objective functional is monotonically decreasing and provides a guaranteed stopping criterion. Moreover, in figure 1e we observe that the number of Degrees of Freedom remains small until the reconstructed interface approaches the real inclusion thus coarse meshes prove to be reliable during the initial iterations of the algorithm. In figure 2, we report the results of a more complex case in which two inclusions are sought. We use ten boundary measurements and the ABVA is still able to identify a genuine descent direction at each iteration ( Fig. 2b) but the reconstructed interface appears to be less precise (Fig. 2a). This phenomenon depends on the severe ill-posedness of the problem and the diffusive nature of the state equation is responsible for the loss of the information in the middle of the domain, even when using several measurements. This results in the ability of the algorithm to reconstruct the inclusions near the boundary of the domain whereas the inner interfaces appear more difficult to identify and with a less precise outcome. As previously observed, figure 2c confirms that the ABVA allows to use the same meshes for several iterations increasing the number of DOFs only when the descent direction is no more validated.

Conclusion
Introducing an a posteriori error estimator for the shape gradient provides quantitative information to define a reliable stopping criterion for the algorithm. Nevertheless, solving the  dual flux problems to derive an upper bound of the error in the shape gradient is computationally expensive hence the ABVA may result in higher computing times than the classical BVA on fine meshes. Prompted by the above promising results, we are now investigating alternative procedures to estimate the error in order to rely only on the computation of local quantities.