G Allaire

C Dapogny

A LINEARIZED APPROACH TO WORST-CASE DESIGN IN PARAMETRIC AND GEOMETRIC SHAPE OPTIMIZATION

The purpose of this article is to propose a deterministic method for optimizing a structure with respect to its worst possible behavior when a 'small' uncertainty exists over some of its features. The main idea of the method is to linearize the considered cost function with respect to the uncertain parameters, then to consider the supremum function of the obtained linear approximation, which can be rewritten as a more 'classical' function of the design, owing to standard adjoint techniques from optimal control theory. The resulting 'linearized worst-case' objective function turns out to be the sum of the initial cost function and of a norm of an adjoint state function, which is dual with respect to the considered norm over perturbations. This formal approach is very general, and can be justified in some special cases. In particular, it allows to address several problems of considerable importance in both parametric and shape optimization of elastic structures, in a unified framework.

Introduction

As idealized visions of reality, most optimization frameworks assume a complete knowledge of the parameters and data of the underlying model. Unfortunately, for a lot of reasons, physical parameters involved in realistic applications are hardly ever known with such exactness, and the feasibility and optimality of the solution to an optimization problem can be tremendously jeopardized by slight variations in the data. In this spirit, an example of the less compliant microstructure for an elastic composite material submitted to a particular set of traction loads is given in [START_REF] Cherkaev | Optimal design for uncertain loading conditions[END_REF], which is infinitely compliant when submitted to any load with a different orientation (see also [4] for a discussion about the sensitivity of linear programming problems with uncertain parameters).

In the more specific context of shape optimization of elastic structures, which is at stake in this paper, optimization problems may be plagued by (at least) three very different types of uncertainties:

• Uncertainties about the location, magnitude and orientation of the body forces or surface loads exerted on shapes: not to mention the fact that they are generally known through error-prone measurements, these external stresses are affected by the outer medium, which may itself undergo unknown perturbations. • Uncertainties about the elastic material's properties: changes in the conditions in the considered medium (temperature, humidity, etc...) may alter the material's stiffness. On a different note, the material's properties could also be perturbed during its manufacturing process, in which small inclusions of 'parasitic' phases may accidentally occur. • Uncertainties about the geometry of the shape itself: due to wear, or to the manufacturing process, the geometry of the shape may cease to be (or may not be from the beginning) the one which was initially forecast.

Robust design has been paid much attention in shape and topology optimization (and in optimization in general). Depending on the available information regarding the uncertainties, the question has been studied from two distinct perspectives.

On the one hand, many authors assume the knowledge of a probabilistic distribution (which is often obtained via statistical studies) as for the behavior of perturbations around an unperturbed state; see for instance [START_REF] Frangopol | Life-cycle reliability-based optimization of civil and aerospace structures[END_REF][START_REF] Kharmanda | Reliability-based topology optimization[END_REF] and references therein. Then, the mean value of the considered objective criterion [START_REF] Conti | Shape optimization under uncertainty -a stochastic programming approach[END_REF][START_REF] Kalsi | A Comprehensive Robust Design Approach for Decision Trade-Offs in Complex Systems Design[END_REF], or a weighted sum of its mean value and standard deviation [START_REF] Chen | A new level-set based approach to shape and topology optimization under geometric uncertainty[END_REF][START_REF] Dunning | Robust Topology Optimization: Minimization of Expected and Variance of Compliance[END_REF] are minimized to guarantee a fine performance, which is relatively independent of the perturbations. Other approaches, referred to as reliability-based approaches (see the overview in [START_REF] Choi | Reliability-based Structural Design[END_REF]) put the emphasis on ensuring that constraints are satisfied in all the possible perturbed configurations. Such a requirement is often modeled under the form of upper bound constraints on so-called failure probabilities. Regardless of the particular considered model, such probabilistic methods generally prove very costly, since they imply repeated evaluations of the mean value or standard deviation of the objective function, or probabilities of violation of constraint functions. This is generally achieved by (expensive) sampling methods (e.g. Monte-Carlo methods), which may be sped up by relying on first or second order approximations of the constraint functions; see e.g. [START_REF] Choi | Reliability-based Structural Design[END_REF] about the so-called First-Order Reliability-based Methods (FORM), and [START_REF] Jung | Reliability-Based Topology Optimization of Geometrically Nonlinear Structures with Loading and Material Uncertainties[END_REF] for an example of use in the context of structural optimization.

On the other hand, when no information is available on the expected perturbations but for bounds on their magnitude, so-called worst-case designs approaches are preferred. The problem can then be rephrased either as that of minimizing the worst value of the objective function among all the perturbed designs [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF][START_REF] De Gournay | Shape and topology optimization of the robust compliance via the level set method[END_REF][START_REF] Lassila | Boundary control and shape optimization for the robust design of bypass anastomoses under uncertainty[END_REF], or guaranteeing that constraints are fulfilled by every such design [START_REF] Guo | Robust structural topology optimization considering boundary uncertainties[END_REF] (the latter approach being called confidence optimization). Such problems are generally formulated as bilevel 'min-max' problems [START_REF] Guo | Confidence structural robust design and optimization under stiffness and load uncertainties[END_REF]: a first problem consists in finding the worst-case perturbation for a given design, then a second one is about finding the optimal design with respect to this supremum criterion. Of course, worst-case design problems are very difficult and computationally expensive in full generality (see for instance [START_REF] Lassila | Boundary control and shape optimization for the robust design of bypass anastomoses under uncertainty[END_REF] for a reduced-basis method adapted to a worst-case design problem). They cannot be solved without resorting to approximations, except in a few very specific cases; see [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF][START_REF] De Gournay | Shape and topology optimization of the robust compliance via the level set method[END_REF] about the important problem of finding the less compliant shape under uncertainties on the body forces. On the theoretical side, worst-case functionals have recently been studied, for instance in [START_REF] Herzog | Weak lower semi-continuity of the optimal value function and applications to worst-case robust optimal control problems[END_REF], where conditions are given for such min-max problems to admit at least one solution. In [START_REF] Hlaváček | On topological derivatives for elastic solids with uncertain input data[END_REF], it is shown that the concept of topological derivative is robust when the linear elasticity system at play undergoes 'small' perturbations.

The aim of this paper is to propose a unified framework for the worst-case design of elastic structures, with respect to 'small' perturbations on the applied body forces and surface loads, on the material's properties and on the geometries of the structures. More precisely, a formal and rather inexpensive approach is presented for the minimization of the maximum value of a given criterion under the assumption of 'small amplitude' perturbations. Note that we do not attempt to tackle the so-called 'confidence worst-case design' approach, which would ask constraints to be satisfied for all the perturbed designs, but claim that the same philosophy would allow to deal with them. The starting ingredient of our approach is to take advantage of the smallness of the perturbations and thus to linearize the cost function with respect to perturbations around the unperturbed configuration. The maximum value of the resulting linear functional among all possible perturbations can be explicitly computed, and is then optimized. This idea is quite natural and it has already been used in some specific cases (for example, [START_REF] Guo | Robust structural topology optimization considering boundary uncertainties[END_REF] considers the case of compliance minimization under geometric uncertainties). Our approach is however systematic and very general. It is formal, for there is a priori no guarantee that the supremum of the linearized cost function should be close to the real worst value of the original criterion. Yet, we shall see that, in some cases, it can be rigorously justified and admits physical interpretation. Besides, in general, it should provide valuable help in discerning 'trends' towards robustness with respect to perturbations of various kinds.

This article is organized as follows: Section 2 opens the discussion with a presentation of the general philosophy of the proposed method in a formal, abstract framework. This method is then carried out in Section 3, in the simpler situation of (parametric) optimization of the thickness of an elastic plate, which already features almost all the salient features of the approach; then, it is used in shape optimization in Section 4. Then, some numerical examples and discussions are proposed in Section 5 to appraise the efficiency of the proposed method, and some technical details are supplied concerning the proposed implementation.

General setting and main notations

This very informal section presents the generic worst-case optimization problem addressed in this paper, and exposes the main ideas of the proposed approach to deal with it. In the meantime, some notations are introduced.

Let H be a set of admissible designs among which the 'best' element is sought, with respect to a prescribed criterion or cost function C(u(h)), depending on h ∈ H via the solution (or state) u(h) to a system where h acts as a parameter, say:

(1)

A(h)u(h) = b.

In this system, the right-hand side b represents the data (typically the applied loads). In the applications ahead, H will stand for either a space of design parameters -like the properties of an elastic material, or the thickness of a plate (see Section 3) -or for a space of geometric shapes (see Section 4).

The system (1) may undergo perturbations that affect the state u(h), thus spoiling the performance of the corresponding design h ∈ H. The set P of such perturbations is assumed to be a Banach space, with norm ||.|| P , and we only assume that the expected perturbations have 'small' maximum norm m > 0.

To keep things simple, let us assume that these perturbations only affect the right-hand side of (1), i.e. b = b(δ), δ ∈ P. As we shall see, this is the case when the optimal thickness (resp. shape) of an elastic plate (resp. structure) is sought under uncertainties over the applied external body forces or traction loads. The state u = u(h, δ) is now solution to:

(2)

A(h)u(h, δ) = b(δ).
By convention, the case δ = 0 accounts for the unperturbed situation, and we shall indifferently denote by u(h, 0) or u(h) the corresponding state when the context is clear. The associated worst-case optimization problem consists in minimizing the functional J : H → R, defined as the maximum value reached by the cost function C(u(h, δ)) for all the potential perturbations δ ∈ P, i.e.

∀h ∈ H, J (h) = sup

δ∈P ||δ|| P ≤m C(u(h, δ)).
Taking advantage of the assumption that the amplitude m of the exerted perturbations is small, and since we claimed that this problem is difficult to solve as such in general, we propose to trade it for that of minimizing another functional J , obtained from J by linearizing the dependence of C(u(h, δ)) on δ before evaluating the supremum: [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF] ∀h ∈ H, J (h) = sup

δ∈P ||δ|| P ≤m C(u(h)) + dC du (u(h)) ∂u ∂δ (h, 0)(δ) ,
where dC du stands for the (total) differential of C with respect to u, and ∂u ∂δ is the (partial) differential of u with respect to δ. Note that the application δ → C(u(h)) + dC du (u(h)) ∂u ∂δ (h, 0)(δ) is affine by construction, since it only involves differentials of functions. Now, the supremum in formula (3) can be computed explicitly, in terms of the norm ||.|| Q of another Banach space Q:

(4) ∀h ∈ H, J (h) = C(u(h)) + m dC du (u(h)) ∂u ∂δ (h, 0) Q .
Here, Q refers to either the dual space of P (i.e. Q = P * ) or to the pre-dual space of P (i.e. P = Q * ). In the former case, (4) is a consequence of the definition of the definition of the dual norm, and in the latter one, it is a corollary of Hahn-Banach's theorem [5].

The resulting expression is still not explicit in terms of h, for it involves the sensitivity ∂u ∂δ (h, 0) of the solution to the state equation with respect to perturbations. However, classical techniques in optimal control theory allow to make it explicit, up to the introduction of an adjoint state p(h), which arises as the solution to an adjoint system very similar to (1):

(5)

A(h) T p(h) = dC du (u(h)),
where A(h) T stands for the adjoint operator to A(h), and we implicitly identified the differential dC du with the gradient of C with respect to some dual pairing. Indeed, differentiating with respect to δ in (2) produces a system for ∂u ∂δ (h, 0):

A(h) ∂u ∂δ (h, 0) = db dδ (0)
; using this latter equation in combination with the definition (5) of p(h) readily leads to:

dC du (u(h)) ∂u ∂δ (h, 0) = A(h) T p(h) ∂u ∂δ (h, 0) = A(h) ∂u ∂δ (h, 0) p(h) = db dδ (0)p(h).
Eventually, J rewrites as:

(6) J (h) = C(u(h)) + m db dδ (0)p(h) Q .
This last expression can be interpreted as follows: the approximate cost function J (h) is an aggregated sum of the unperturbed cost function C(u(h)) and a penalization of the perturbations || db dδ (0)p(h)|| Q , the penalization parameter m being precisely the expected magnitude of perturbations.

Of course, the above argument is very formal, since the involved expressions mix blithely the notions of differential applications and gradients, as well as the spaces associated to perturbations, state variables, etc... As we shall see, a significant part of the work consists in giving a precise meaning to this rough sketch.

Remark 1. It is no surprise that formula (6) for functional J features the adjoint state p(h); it is indeed well-known that the adjoint measures the sensitivity of the optimal cost function C with respect to perturbations on the right-hand side of the state equation [START_REF] Allaire | Conception optimale de structures[END_REF] (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF], rem. 4.14 and 5.21).

The previous analysis leaves us with a more classical state-constrained optimization problem, save that the functional to be minimized itself brings an adjoint state into play. The computation of the derivative of J can however be carried out as follows:

• Derivating the part C(u(h)) in [START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF] does not pose any further difficulty: it is merely the unperturbed cost function of the considered problem. Doing so involves the already introduced adjoint state p(h), which accounts for the sensitivity of C with respect to its argument. • Derivating the second part || db dδ (0)p(h)|| Q is a little bit more tricky, and brings into play two further adjoint states q(h) and z(h). The first one q(h) has nothing to do with p(h) and expresses the sensitivity of the new 'objective function' || db dδ (0) • || Q with respect to its argument. As we shall see, even when the unperturbed optimization problem is self-adjoint (i.e. p(h) = ±u(h)), q(h) differs from ±u(h). The second one z(h) describes the sensitivity of p(h) -that is, in some way, of dC du -and typically involves second order derivatives of C with respect to its argument. When the unperturbed problem is self-adjoint, z(h) happens to be equal to ±q(h).

Remark 2. For simplicity of the exposition, we did not consider the possibility that the problem of minimizing J (or J ) may be constrained. Actually, in the remainder of this article, we will only be considering constraints on the volumes of structures, which are especially easy to enforce in our context (see Section 5). However, there is no theoretical difficulty in extending our approach to worst-case design problems with confidence constraints, i.e. problems of the form:

min h∈H max δ∈P C(u(h, δ)) s.t. max δ∈P c i (u(h, δ)) ≤ 0, ∀i = 1, ..., N,
where c i are scalar-valued functions of the state u.

Notations. Throughout this article, we consistently denote as C the various cost functions under consideration. Note that, contrary to the basic setting presented above, such cost functions may themselves depend on the perturbations (which is not subsequently problematic). We will denote as J the associated worst-case design functional, and as J the approximate worst-case design functional, obtained by the aforementioned linearization of the cost function.

Worst-case design in parametric optimization

Description of the model problem.

Throughout Section 3, we consider a thin linear elastic plate in plane stress situation, with Lipschitz cross-section Ω ⊂ R d (d = 2 in concrete applications) and positive thickness h ∈ L ∞ (Ω) (see Figure 1).

Ω x • h(x) • Γ D Γ N g Figure 1.
Thin plate with cross section Ω and thickness h, clamped on a part of its boundary corresponding to Γ D ⊂ ∂Ω, and submitted to transversal loads, applied on a part associated to Γ N ⊂ ∂Ω.

From the plane stress assumption, the equilibrium equation of the plate can be written as a d-dimensional system posed on Ω. More specifically, assume the plate is clamped on a part of its boundary, associated to the subset Γ D of ∂Ω. Surface loads g ∈ L 2 (Γ N ) d are applied on the complementary part Γ N := ∂Ω \ Γ D of ∂Ω. Denote also as f ∈ L 2 (Ω) d the body force exerted on the plate.

The in-plane displacement function u belongs to the space H 1 Γ D (Ω) d , where we have defined:

H 1 Γ D (Ω) := u ∈ H 1 (Ω), u = 0 on Γ D .

It arises as the unique solution in H 1

Γ D (Ω) d of the linear elasticity system in Ω:

(7)    -div(hAe(u)) = f in Ω u = 0 on Γ D hAe(u)n = g on Γ N
, where e(u) = (∇u T + ∇u)/2 is the strain tensor, n : ∂Ω → S d-1 is the unit outer normal vector field to Ω, A is the material Hooke's law, defined for any e ∈ S(R d ) (the set of d × d real symmetric matrices) by [START_REF] Chen | A new level-set based approach to shape and topology optimization under geometric uncertainty[END_REF] Ae = 2µe + λtr(e)I, with the Lamé coefficients λ, µ, satisfying µ > 0 and λ + 2µ/d > 0.

Our general purpose is to optimize the thickness h of the considered plate among a set U ad ⊂ L ∞ (Ω) of admissible thickness functions, with respect to a criterion yet to be specified. Throughout Section 3, we shall simply work with:

(9) U ad = {h ∈ L ∞ (Ω), s.t. a.e. x ∈ Ω, h min ≤ h(x) ≤ h max } ,
where 0 < h min < h max are prescribed lower and upper bounds for thickness functions. Depending on the context, we may also need to impose a constraint on the volume of the plate; yet, to keep notations as generic as possible, we chose not to incorporate it into the modeling of the present section (see Section 5 for details about the numerical treatment of this constraint). Note also that it is customary to add a uniform smoothness constraint on h in order to easily obtain an existence result for the optimal design problem. Since existence of optimal designs is not the focus of the present paper, we ignore such smoothness constraints. The remainder of Section 3 is now dedicated to illustrating the general guideline of Section 2 on a series of model problems.

Worst-case design of an elastic plate under perturbations on the body forces.

Let us start our study with the optimization of the thickness of the considered plate in the worst-case scenario when the applied body forces f are perturbed as f + ξ, for small ξ ∈ L 2 (Ω) d .

Let j : R d f × R d u → R be a function of class C 2 , which complies with the following growth conditions:

(10) ∀f ∈ R d , u ∈ R d ,    |j(f, u)|≤ C(|f | 2 +|u| 2 ) |∇ f j(f, u)|≤ C(|f |+|u|), |∇ u j(f, u)|≤ C(|f |+|u|) |∇ 2 f j(f, u)|≤ C, |∇ f ∇ u j(f, u)|≤ C, |∇ 2 u j(f, u)|≤ C
, for some positive constant C > 0. This function j could also depend on the space variable x ∈ R d , without any change in the following developments and conclusions; to keep notations simple, this dependence is omitted (as will be the case in all the considered problems). Note that ∇ u j(f, u) and ∇ f j(f, u) are vectors in R d , and

∇ 2 f j(f, u), ∇ f ∇ u j(f, u) = ∇ 2 u j(f, u) are d × d matrices.
To emphasize the dependence of the elastic displacement u on the plate thickness h and the body forces f , we denote by u h,f ∈ H 1 Γ D (Ω) d the unique solution to problem (7) using these data. The cost function for a plate of thickness h with loads f is then defined as:

(11) C(h, f ) = Ω j(f, u h,f ) dx.
To set ideas, we have only assumed this cost to depend on u h,f -and not on its gradient-, as an integral expression on Ω -not on its boundary -, but it would be easy to generalize the discussion ahead to such cases (see Section 4.4 and [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF]).

Let us now fix f ∈ L 2 (Ω) d and g ∈ L 2 (Γ N ) d . As in Section 2, when no confusion is possible, we will denote indifferently u h = u h,f the solution to the unperturbed problem. The worst-case design problem in our situation reads: [START_REF] Ciarlet | Mathematical Elasticity, vol I: Three Dimensional Elasticity[END_REF] min

h∈U ad J (h), where J (h) = sup ξ∈L 2 (Ω) d ||ξ|| L 2 (Ω) d ≤m C(h, f + ξ).
Following the general guideline of Section 2, we propose to trade problem [START_REF] Ciarlet | Mathematical Elasticity, vol I: Three Dimensional Elasticity[END_REF] for a new one where C(h, f ) has been linearized with respect to f , namely: [START_REF] Conti | Shape optimization under uncertainty -a stochastic programming approach[END_REF] min

h∈U ad J (h), where J (h) = sup ξ∈L 2 (Ω) d ||ξ|| L 2 (Ω) d ≤m C(h, f ) + ∂C ∂f (h, f )(ξ) .
The main result of this section is an explicit formula for [START_REF] Conti | Shape optimization under uncertainty -a stochastic programming approach[END_REF] and the computation of its derivative with respect to the thickness h, thanks to the introductions of several adjoint problems.

Theorem 1. The functional J , defined by (13) rewrites, for any h ∈ U ad , as

(14) J (h) = Ω j(f, u h ) dx + m ||∇ f j(f, u h ) -p h || L 2 (Ω) d ,
where

p h ∈ H 1 Γ D (Ω) d
is the first adjoint state, defined as the unique solution to

(15)    -div(hAe(p)) = -∇ u j(f, u h ) in Ω, p = 0 on Γ D , hAe(p)n = 0 on Γ N . Let h ∈ U ad such that ∇ f j(f, u h ) -p h = 0 in L 2 (Ω) d .
Then J is Fréchet differentiable at h, and its directional derivative reads, for any s ∈ L ∞ (Ω), as

(16) J ′ (h)(s) = Ω D(u h , p h , q h , z h ) s dx,
where

D(u h , p h , q h , z h ) := Ae(u h ) : e(p h ) + m 2 ||∇ f j(f, u h ) -p h || L 2 (Ω) d
Ae(u h ) : e(z h ) + Ae(p h ) : e(q h ) .

In the last formula, the second, and third adjoint states q h , z h ∈ H 1 Γ D (Ω) d are defined as the unique solutions to, respectively:

(17)    -div(hAe(q)) = -2 (p h -∇ f j(f, u h )) in Ω, q = 0 on Γ D , hAe(q)n = 0 on Γ N , (18)
   -div(hAe(z)) = -2 ∇ f ∇ u j(f, u h ) T (∇ f j(f, u h ) -p h ) -∇ 2 u j(f, u h )q h in Ω, z = 0 on Γ D , hAe(z)n = 0 on Γ N .
Proof. We first prove formula [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] for J . Since the map f → u h,f is linear and continuous from L 2 (Ω) d into H 1 Γ D (Ω) d , it is obviously Fréchet differentiable. Applying the chain rule lemma to the computation of ∂C ∂f , we deduce from (13) the following formula for J :

(19) J (h) = Ω j(f, u h )dx + sup ξ∈L 2 (Ω) d ||ξ|| L 2 (Ω) d ≤m Ω ∇ f j(f, u h ) • ξ dx + Ω ∇ u j(f, u h ) • ∂u h,f ∂f (ξ) dx .
To simplify the second term, let us use the variational formula for the derivative ∂u h,f ∂f (ξ) obtained by differentiating that of u h,f with respect to f . We get, for any ξ ∈ L 2 (Ω) d :

∀v ∈ H 1 Γ D (Ω) d , Ω hAe ∂u h,f ∂f (ξ) : e(v) dx = Ω ξ • v dx.
Now introduce the first adjoint state p h ∈ H 1 Γ D (Ω) d , the unique solution of [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. It comes from the associated variational formulation that:

Ω -∇ u j(f, u h, ) • ∂u h,f ∂f (ξ) dx = Ω hAe(p h ) : e ∂u h,f ∂f (ξ) dx = Ω p h • ξ dx.
In this view, [START_REF] Frangopol | Life-cycle reliability-based optimization of civil and aerospace structures[END_REF] becomes:

(20)

J (h) = Ω j(f, u h ) dx + sup ξ∈L 2 (Ω) d ||ξ|| L 2 (Ω) d ≤m Ω ∇ f j(f, u h ) • ξdx - Ω p h • ξ dx = Ω j(f, u h ) dx + m ||∇ f j(f, u h ) -p h || L 2 (Ω) d ,
where the supremum with respect to ξ is achieved by

ξ * = m(∇ f j(f, u h ) -p h )/ ||∇ f j(f, u h ) -p h || L 2 (Ω) d .
In a second step we prove that J is differentiable with respect to the thickness h and we compute its derivative. This is a classical issue and the whole problem of worst-case design with respect to an uncertainty on the body forces has been concealed in the adjoint state p h . As expected, the first term in the right hand side of [START_REF] De Gournay | Shape and topology optimization of the robust compliance via the level set method[END_REF] is the unperturbed cost function and its differentiation is a classical result recalled in Lemma 2 below. We obtain, for any s ∈ L ∞ (Ω), [START_REF] Guo | Robust structural topology optimization considering boundary uncertainties[END_REF] ∂ ∂h

Ω j(f, u h ) dx (s) = Ω s Ae(u h ) : e(p h ) dx.
To be able to differentiate the second term

||∇ f j(f, u h ) -p h || L 2 (Ω) d , we have to assume that h is such that ∇ f j(f, u h ) -p h = 0 in L 2 (Ω) d
. This is a reasonable hypothesis, meaning that the considered plate with thickness h is not indifferent to a small perturbation on the applied body forces. Excluding such a case, and using again Lemma 2 with a function ℓ(u, p) = |∇ f j(f, u) -p| 2 , an elementary calculation yields

(22) ∂ ∂h ||∇ f j(f, u h ) -p h || L 2 (Ω) d (s) = Ω s Ae(u h ) : e(z h ) + Ae(p h ) : e(q h ) dx 2 ||∇ f j(f, u h ) -p h || L 2 (Ω) d ,
where q h and z h are solutions of ( 17) and [START_REF] Dunning | Robust Topology Optimization: Minimization of Expected and Variance of Compliance[END_REF], respectively. Eventually, combining ( 21) and ( 22) delivers the desired formula.

In the proof of Theorem 1, we used the following general lemma for differentiating functionals depending on h via u h and the adjoint state p h , solution to [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. For simplicity we shall not indicate the dependence of the integrand j with respect to f . Lemma 2. For h ∈ U ad , denote as u h ∈ H 1 Γ D (Ω) d the unique solution to problem [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]. (i) Let j : R d → R be a function of class C 1 , which fulfills the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. Consider the functional

K(h) = Ω j(u h ) dx.
Then K is Fréchet differentiable at any h ∈ U ad and its derivative reads, for any s ∈ L ∞ (Ω),

K ′ (h)(s) = Ω sAe(u h ) : e(p h ) dx, (23)
where the adjoint state p h ∈ H 1 Γ D (Ω) d is the unique solution of the system (15). (ii) Suppose moreover that j is of class C 2 and let ℓ : R d u × R d p → R be a function of class C 1 , both of them satisfying the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. Consider the functional

L(h) = Ω ℓ(u h , p h ) dx,
where p h is defined by [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. Then L is Fréchet differentiable at any h ∈ U ad and its derivative reads:

∀s ∈ L ∞ (Ω), L ′ (h)(s) = Ω s Ae(u h ) : e(z h ) + Ae(p h ) : e(q h ) dx, where q h ∈ H 1 Γ D (Ω) d is the unique solution to (24)    -div(hAe(q)) = -∇ p ℓ(u h , p h ) in Ω, q = 0 on Γ D , hAe(q)n = 0 on Γ N ,
and

z h ∈ H 1 Γ D (Ω) d is the unique solution to (25)    -div(hAe(z)) = -∇ u ℓ(u h , p h ) -∇ 2 u j(u h )q h in Ω, z = 0 on Γ D , hAe(z)n = 0 on Γ N .
The proof of Lemma 2 requires the following well-known result about the differentiability of the solution u h to [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF] with respect to the thickness h. This is a classical result in optimal control theory, the proof of which can be found e.g. in [START_REF] Allaire | Conception optimale de structures[END_REF], or [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], Chap. 5 in a more general case, which is typically obtained by applying the implicit function theorem.

Lemma 3. Let u h ∈ H 1 Γ D (Ω) d be the solution to (7). The mapping h → u h , from U ad into H 1 Γ D (Ω) d , is Fréchet differentiable.
Proof of Lemma 2. (i): This first point is again a classical result in optimal control theory. However, for the sake of completeness, we briefly recall its derivation using Céa's method [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]. By virtue of Lemma 3 the application h → u h is differentiable; as a consequence, the forthcoming argument is not only formal, as in most cases where Céa's method generally comes in handy. Introduce the Lagrangian L :

U ad × H 1 Γ D (Ω) d × H 1 Γ D (Ω) d → R, defined as: (26) L(h, u, p) = Ω j( u) dx + Ω hAe( u) : e( p) dx - Ω f • p dx - Γ N g • p ds,
and let us search for the points (u, p) where the partial derivatives of L(h, ., .) cancel, for a given h ∈ U ad . First, the cancellation of the derivative of

L at (h, u, p) ∈ U ad × H 1 Γ D (Ω) d × H 1 Γ D (Ω) d with respect to p reads: (27) ∀ p ∈ H 1 Γ D (Ω) d , ∂L ∂p (h, u, p)( p) = Ω hAe(u) : e( p) dx - Ω f • p dx - Γ N g • p ds = 0,
which is just the variational formulation for problem [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]; hence u = u h , solution to [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF].

Next, cancelling the derivative of L with respect to u leads to the variational formulation for the adjoint state:

(28) ∀ u ∈ H 1 Γ D (Ω) d , ∂L ∂u (h, u, p)( u) = Ω ∇ u j(u) • u dx + Ω hAe( u) : e(p) dx = 0,
which allows to identify p as p h defined by [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. By definition of u h , we now have:

∀h ∈ L ∞ (Ω), ∀ p ∈ H 1 Γ D (Ω) d , K(h) = L(h, u h , p
). Thus, differentiating this last expression with respect to h yields:

∀ p ∈ H 1 Γ D (Ω) d K ′ (h)(s) = ∂L ∂h (h, u h , p)(s) + ∂L ∂u (h, u h , p)( ∂u h ∂h (s)).
Now taking p = p h in the previous expression and using [START_REF] Kharmanda | Reliability-based topology optimization[END_REF] lead to:

K ′ (h)(s) = ∂L ∂h (h, u h , p h )(s).
Eventually, differentiating [START_REF] Jung | Reliability-Based Topology Optimization of Geometrically Nonlinear Structures with Loading and Material Uncertainties[END_REF] with respect to h produces the desired formula [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF].

(ii): Similarly, the application h → p h is differentiable (again by virtue of Lemma 3). Let us introduce the (different) Lagrangian L :

U ad × H 1 Γ D (Ω) d × H 1 Γ D (Ω) d → R, defined as: L(h, p, q) = Ω ℓ(u h , p) dx + Ω hAe( p) : e( q) dx + Ω ∇ u j(u h ) • q dx. Searching for a point (p, q) ∈ H 1 Γ D (Ω) d × H 1 Γ D (Ω) d at
which the derivative of L(h, ., .) with respect to q vanishes yields:

(29) ∀ q ∈ H 1 Γ D (Ω) d , ∂L ∂q (h, p, q)( q) = Ω hAe(p) : e( q) dx + Ω ∇ u j(u h ) • q dx = 0,
and we find p = p h , the solution to [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. Now cancelling the derivative of L with respect to p one finds:

(30) ∀ p ∈ H 1 Γ D (Ω) d , ∂L ∂p (h, p, q)( p) = Ω ∇ p ℓ(u h , p) • p dx + Ω hAe( p) : e(q) dx = 0,
and since p = p h , we identify q = q h the unique solution to [START_REF] Herzog | Weak lower semi-continuity of the optimal value function and applications to worst-case robust optimal control problems[END_REF]. Consequently, we have, for any h ∈ U ad , and any q ∈ H 1 Γ D (Ω) d : L(h) = L(h, p h , q). As in the previous point, differentiating this expression with respect to h, then taking q = q h , we end up with:

L ′ (h)(s) = ∂L ∂h (h, p h , q h )(s).
Now, we are left with computing this partial derivative ∂L ∂h . Since L depends on h ∈ U ad via u h too, we use point (i): for fixed functions p, q

∈ H 1 Γ D (Ω) d (independent of h), define the C 1 function m : R d → R as: m(u) = ℓ(u, p) + ∇ u j(u) • q.
Note that m actually also depends on the space variable x ∈ R d , and this dependence is omitted. Applying the first point to the map h → Ω m(u h ) dx and introducing

z h ∈ H 1 Γ D (Ω) d as the unique solution to the system    -div(hAe(z)) = -∇ u ℓ(u h , p) -∇ 2 u j(u h )q in Ω z = 0 on Γ D hAe(z)n = 0 on Γ N ,
we end up with:

∂L ∂h (h, p, q)(s) = ∂ ∂h Ω ℓ(u h , p) dx + Ω ∇ u j(u h ).q dx (s) + Ω sAe(p) : e(q) dx = Ω sAe(u h ) : e( z h ) dx + Ω sAe(p) : e(q) dx.
Note that, in general, z h = z h since p = p h or q = q h . Now we take p = p h and q = q h , so that z h = z h , and we get the desired formula for L ′ (h)(s).

Remark 3. Using Céa's method is not the only possibility (neither the simpler) for proving Theorem 1.

Alternatively, one can differentiate the variational problems defining u h and p h , apply the chain rule lemma for the derivative of J with respect to the thickness h and simplify the resulting expression by using the two new adjoint states q h and z h . This is rather simple in the present context of parametric optimization (see for instance the first step of the proof of Theorem 1). Yet, it is much more cumbersome when it comes to doing the same in geometric shape optimization. To emphasize the similarities between both settings, we thought it better to prove Lemma 2 using Céa's method, which will be the right tool in Section 4 too.

Example 1. For the sake of simplicity, suppose that no surface loads are applied -g = 0 (however, the argument would adapt mutatis mutandis to the general case; see [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF]) -, and that we are interested in the compliance as a cost function, i.e. j(f, u) = f • u in [START_REF] Choi | Reliability-based Structural Design[END_REF]. The various derivatives of j are: 15) is p h = -u h (the unperturbed problem is self-adjoint, as is well-known in this case). Then, J , defined by [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF], has the following expression:

∇ u j(f, u) = f, ∇ f j(f, u) = u, ∇ 2 u j(u) = 0, ∇ f ∇ u j(f, u) = I. For any h ∈ U ad , the unique solution p h ∈ H 1 Γ D (Ω) d to problem (
(31) ∀h ∈ U ad , J (h) = Ω f • u h dx + 2m ||u h || L 2 (Ω) d . Furthermore, q h ∈ H 1 Γ D (Ω) d is the unique solution to the problem:    -div(hAe(q)) = 4u h in Ω q = 0 on Γ D hAe(q)n = 0 on Γ N ,
and from [START_REF] Dunning | Robust Topology Optimization: Minimization of Expected and Variance of Compliance[END_REF], one acknowledges that z h = -q h , which gives the straightforward expression for the gradient of J (h), for any h ∈ U ad such that u h = 0:

∀s ∈ L ∞ (Ω), J ′ (h)(s) = - Ω s Ae(u h ) : e(u h ) + m 2 ||u h || L 2 (Ω) d
Ae(u h ) : e(q h ) dx.

Remarks 1.

• Interestingly enough, formula [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF] expresses the fact that -at first order -optimizing the worst-case compliance when uncertainties around body forces are expected translates into a penalization of the unperturbed compliance by the norm of the displacement u h of the structure. • As we already mentioned in the introduction, in this particular case where the cost function C is the compliance, the study of the exact worst-case functional J defined in ( 12) can be addressed without linearization of C, resorting to more involved techniques [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF][START_REF] De Gournay | Shape and topology optimization of the robust compliance via the level set method[END_REF].

There are many variants of the perturbation setting considered in Theorem 1. Let us describe briefly another example for which Theorem 1 and Lemma 2 can easily be extended (the reader is refered to [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] for other examples).

Example 2 (Localization and restriction of the direction of perturbations). Still assuming that g = 0, let η ∈ L ∞ (Ω) d be a fixed vector field on Ω; the most important case we have in mind is when η is 0 everywhere except on a small portion of Ω, where it is a constant unit direction. The underlying idea is to investigate perturbed body forces of the kind

(f + ξη), ξ ∈ L 2 (Ω), ||ξ|| L 2 (Ω) ≤ m.
Still considering the compliance as an objective function, for any h ∈ U ad , the worst-case functional J of interest is in this case:

(32) J (h) = sup ξ∈L 2 (Ω) ||ξ|| L 2 (Ω) ≤m C(h, f + ξη) = Ω (f + ξη) • u f +ξη dx .
Similarly to [START_REF] Conti | Shape optimization under uncertainty -a stochastic programming approach[END_REF], the linearized worst-case functional is, for any h ∈ U ad ,

J (h) = sup ξ∈L 2 (Ω) ||ξ|| L 2 (Ω) ≤m C(h, f ) + ∂C ∂f (h, f )(ξη)
Theorem 4. The functional J rewrites:

(33) J (h) = Ω f • u h dx + 2m ||η • u h || L 2 (Ω) .
Let h ∈ U ad be such that η • u h = 0 in L 2 (Ω). Then J is Fréchet-differentiable at h, and its differential reads:

∀s ∈ L ∞ (Ω), J ′ (h)(s) = - Ω s Ae(u h ) : e(u h ) + m 2 ||η • u h || L 2 (Ω) Ae(u h ) : e(q h ) dx,
where the adjoint state

q h ∈ H 1 Γ D (Ω) d is the unique solution to:    -div(hAe(q)) = 4(η • u h )η in Ω, q = 0 on Γ D , hAe(q)n = 0 on Γ N .
3.3. Parametric optimization of a worst-case scenario problem under geometric uncertainty.

We now investigate perturbations of a different nature, searching for the optimal thickness h ∈ L ∞ (Ω) of the considered plate when robustness is expected with respect to uncertainties over the thickness of the plate itself. As we have already mentioned in the introduction, such a problem typically occurs in the case of mechanical parts which are likely to undergo high stress during their use, thus to wear out, or of mechanical parts whose manufacturing process is especially error-prone.

More precisely, let j, k : R d → R be two functions of class C 2 , satisfying the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. For any h ∈ U ad , denote as u h the solution to problem [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF] and introduce the cost of the design with thickness h:

C(h) := Ω j(u h ) dx + Γ N k(u h ) ds.
Modeling uncertainties over the geometry (i.e. thickness) of the plate itself demands first to address an important issue concerning the perturbed designs: if h ∈ U ad and s ∈ L ∞ (Ω) is a 'small' perturbation over h, the thickness (h + s) of the corresponding perturbed design may not belong to U ad (although it is still uniformly bounded away from 0 and ∞). However, we believe this is part of the modeling, for designs generally end up perturbed in an accidental way, and there is no particular reason that a perturbed design should still fulfill any imposed constraint. Furthermore, not enforcing that (h + s) should belong to U ad allows for an easier mathematical study.

Let m < h min the maximum expected amplitude of the uncertainty over the thickness h. We consider the optimization problem: [START_REF] Pironneau | FreeFem++ version 2[END_REF] min

h∈U ad J (h), where J (h) = sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m C(h + s).
As alluded to in Section 3.1, this problem is very difficult to tackle in such form. However, in the particular situation where the cost C(h) is the compliance of the structure, it turns out almost trivial, meaning that the worst case in (34) can be found analytically.

Proposition 5. Suppose that the cost function C(h) is the compliance, that is:

∀h ∈ U ad , C(h) = Ω hAe(u h ) : e(u h ) dx = Ω f • u h dx + Γ N g • u h ds.
Then, the exact worst-case functional J is equivalently given by J (h) = C(h -m). Simply put, the worst case with respect to the compliance, when there is an uncertainty of maximum amplitude m over h, is the corresponding structure with thickness (h -m), which is the less rigid (thinner) perturbed structure.

Proof. This is a simple consequence of the elastic energy minimization principle. One has:

C(h) = -2 inf u∈H 1 Γ D (Ω) d 1 2 Ω hAe(u) : e(u) dx - Ω f • u dx - Γ N g • u ds .
Hence, by inverting the two maximizations, we get

J (h) = sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m sup u∈H 1 Γ D (Ω) d 2 Ω f • u dx + 2 Γ N g • u ds - Ω (h + s)Ae(u) : e(u) dx = sup u∈H 1 Γ D (Ω) d 2 Ω f • u dx + 2 Γ N g • u ds - Ω (h -m)Ae(u) : e(u) dx
, which allows us to conclude.

In the general, non trivial setting (i.e. when C(h) is not the compliance), we propose to reformulate our optimization problem according to the general principle of Section 2: [START_REF] Schevenels | Robust topology optimization accounting for spatially varying manufacturing errors[END_REF] min

h∈U ad J (h), where J (h) = sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m C(h) + ∂C ∂h (h)(s) ,
and the following result makes it possible to build a gradient-based algorithm for this simplified minimization problem.

Theorem 6. The functional J , defined as [START_REF] Schevenels | Robust topology optimization accounting for spatially varying manufacturing errors[END_REF] rewrites:

∀h ∈ U ad , J (h) = Ω j(u h ) dx + Γ N k(u h ) ds + m ||Ae(u h ) : e(p h )|| L 1 (Ω) ,
where p h ∈ H 1 Γ D (Ω) d is the first adjoint state, defined as the unique solution to

(36)    -div(hAe(p)) = -∇ u j(u h ) in Ω p = 0 on Γ D hAe(p)n = -∇ u k(u h ) on Γ N .
Moreover, J is differentiable at any h ∈ U ad such that the set

(37) E h := {x ∈ Ω, Ae(u h ) : e(p h ) = 0}
has zero Lebesgue measure, and its differential at such a point reads, for any s ∈ L ∞ (Ω),

J ′ (h)(s) = Ω s Ae(u h ) : e(p h ) + m (Ae(p h ) : e(q h ) + Ae(u h ) : e(z h )) dx,
the second and third adjoint states q h , z h ∈ H 1 Γ D (Ω) d being respectively defined as the unique solutions to:

(38)    -div(hAe(q)) = div(εAe(u h )) in Ω q = 0 on Γ D hAe(q)n = -εAe(u h )n on Γ N , (39)
   -div(hAe(z)) = -∇ 2 u j(u h )q h + div(εAe(p h )) in Ω z = 0 on Γ D hAe(z)n = -εAe(p h )n -∇ 2 u k(u h )q h on Γ N , and ε ∈ L ∞ (Ω) is defined as ε = sgn (Ae(u h ) : e(p h )).
Proof. The proof is close in essence to that of Theorem 1 and we illustrate another way to get the desired expressions, without relying on Céa's method (see Remark 3). Knowing that h → u h and h → p h are differentiable as functions from U ad into H 1 Γ D (Ω) d (see Lemma 3), we achieve a variational formulation for their derivatives by differentiating the ones of ( 7) and [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF], then introducing the adjoint states q h and z h with simple algebraic manipulations (in the parametric setting).

• Expression of J as a function of h using an adjoint state. Computing ∂C ∂h , we deduce from [START_REF] Schevenels | Robust topology optimization accounting for spatially varying manufacturing errors[END_REF] the following formula for J (h), for any h ∈ U ad ,

J (h) = Ω j(u h ) dx + Γ N k(u h ) ds + sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m Ω ∇ u j(u h ) • ∂u h ∂h (s) dx + Γ N ∇ u k(u h ) • ∂u h ∂h (s) ds .
Then, using the variational formulation of the adjoint problem [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF] with the test function ∂u h ∂h (s), we deduce

J (h) = Ω j(u h ) dx + Γ N k(u h ) ds + sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m - Ω hAe(p h ) : e ∂u h ∂h (s) ds .
On the other hand, differentiating the variational formulation of ( 7) yields for any

s ∈ L ∞ (Ω) (40) ∀v ∈ H 1 Γ D (Ω) d , Ω sAe(u h ) : e(v) dx = - Ω hAe ∂u h ∂h (s) : e(v) dx.
Therefore, we obtain

J (h) = Ω j(u h ) dx + Γ N k(u h ) ds + sup s∈L ∞ (Ω) ||s|| L ∞ (Ω) ≤m Ω sAe(u h ) : e(p h ) ds .
Taking s(x) to be of the same sign as Ae(u h ) : e(p h )(x) and of absolute value equal to m yields the maximizer in the above expression and the desired expression for J .

• Computation of the Fréchet derivative of J . Assume now that the set E h defined by [START_REF] Soko | Introduction to shape optimization: shape sensitivity analysis[END_REF] has zero Lebesgue measure, and introduce ε := sgn (Ae(u h ) : e(p h )).

We then have, for all s ∈ L ∞ (Ω),

J ′ (h)(s) = ∂ ∂h Ω j(u h ) dx + Γ N k(u h ) ds h (s) +m Ω ε Ae ∂u h ∂h (s) : e(p h ) dx + Ω ε Ae(u h ) : e ∂p h ∂h (s) dx ,
where we used Lemma 7. As usual, the first term actually rewrites, from the definition of p h by [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF],

∂ ∂h Ω j(u h ) dx + Γ N k(u h ) ds h (s) = Ω sAe(u h ) : e(p h ) dx.
As for the second term, differentiating with respect to h directly in the variational formulation satisfied by p h , we get, for all s ∈ L ∞ (Ω), and

any v ∈ H 1 Γ D (Ω) d , Ω hAe ∂p h ∂h (s) : e(v) dx = - Ω sAe(p h ) : e(v) dx - Ω ∇ 2 u j(u h ) ∂u h ∂h (s) • v dx - Γ N ∇ 2 u k(u h ) ∂u h ∂h (s) • v ds .
Comparing the above variational formulation for ∂p h ∂h (s) with the test function q h and the variational formulation for q h (solution to [START_REF] Wang | On projection methods, convergence and robust formulations in topology optimization[END_REF]) with the test function ∂p h ∂h (s), we get:

Ω ε Ae(u h ) : e ∂p h ∂h (s) dx = - Ω h Ae(q h ) : e ∂p h ∂h (s) dx = Ω sAe(p h ) : e(q h ) dx + Ω ∇ 2 u j(u h ) ∂u h ∂h (s) • q h dx + Γ N ∇ 2 u k(u h ) ∂u h ∂h (s) • q h ds.
Now, similarly, we compare the variational formulation for z h (solution to [START_REF] Wang | A level set method for structural topology optimization[END_REF]) with the test function ∂u h ∂h (s) and the variational formulation (40) for ∂u h ∂h (s) with the test function z h . This yields, for all s ∈ L ∞ (Ω):

Ω ε Ae ∂u h ∂h (s) : e(p h ) dx + Ω ∇ 2 u j(u h ) ∂u h ∂h (s) • q h dx + Γ N ∇ 2 u k(u h ) ∂u h ∂h (s) • q h ds = - Ω hAe(z h ) : e ∂u h ∂h (s) dx = Ω sAe(z h ) : e(u h ) dx,
thus ending the proof.

In the course of the proof of Theorem 6, we used the following lemma for the differentiation of the L 1 norm (see [START_REF] Pinkus | On L 1 -approximation[END_REF], Theorem 2.1 for a proof):

Lemma 7. Let Ω ⊂ R d an
open, bounded domain, and denote as f :

L 1 (Ω) → R the L 1 norm function: ∀u ∈ L 1 (Ω), F (u) = Ω |u(x)| dx.
f is then convex, and its subgradient ∂f (u) ⊂ L ∞ (Ω) at any point u ∈ L 1 (Ω) reads:

λ ∈ ∂f (u) ⇔    λ(x) = 1 if u(x) > 0 λ(x) = -1 if u(x) < 0 λ(x) ∈ [-1, 1] if u(x) = 0 . As a consequence, if u ∈ L 1 (Ω) is such that {x ∈ Ω, u(x) = 0} is of null Lebesgue measure, F is Fréchet- differentiable at u,
and its differential dF (u) reads:

∀v ∈ L 1 (Ω), dF (u)(v) = Ω sgn(u(x))v(x) dx.
Example 3. Consider once again the model case of the minimization of the compliance, assuming g = 0 (for simplicity). Then, j(u) = f • u, k(u) = 0, and one has ∇ u j(u) = f and ∇ 2 u j(u) = 0. Since p h ∈ H 1 Γ D (Ω) d is the unique solution to [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF], one easily finds that p h = -u h . Then, the approximate worst-case functional J rewrites:

J (h) = Ω f • u h dx + m ||Ae(u h ) : e(u h )|| L 1 (Ω) .
Assume the thickness h ∈ U ad is such that the set E h := {x ∈ Ω, Ae(u h ) : e(u h ) = 0} has zero Lebesgue measure. Then we have ε = sgn (-Ae(u h ) : e(u h )) = -1 in Ω. In this view, it is then easy to show that the second and third adjoint states reduce to q h = u h and z h = -u h . Eventually, the differential of J reads:

J ′ (h)(s) = -(1 + 2m) Ω sAe(u h ) : e(u h ) dx.
Remark 4. The results of Example 3 have to be compared with those of Proposition 5 which is concerned with the worst-case compliance (without linearization). In particular, it allows us to give a mathematically rigorous approximation result. Indeed, Proposition 5 gives the exact formula J (h) = C(h -m) from which we deduce by a simple computation that, for a given h ∈ U ad ,

∀s ∈ L ∞ (Ω), J ′ (h)(s) = J ′ (h)(s) + O(m 2 ),
where the O(m 2 ) is uniform with respect to s.

3.4.

Worst-case design with uncertainties over the elastic material's properties.

The proposed method can also be used to deal with perturbations on elastic material's properties. Consider a plate made of a material with Lamé moduli that are functions λ, µ ∈ L ∞ (Ω) such that:

(41) µ(x) ≥ γ > 0, λ + 2µ d ≥ γ + 2γ d > 0, a.e x ∈ Ω.
Denote by A λ,µ e = 2µe + λtr(e) the corresponding (possibly inhomogeneous) Hooke's tensor and by u h,λ,µ ∈ H 1 Γ D (Ω) d the unique solution to the elasticity system (7) when h is the thickness of the plate and the Lamé moduli of the constituent material are λ, µ. We shall also denote as u h := u h,λ,µ , and A := A λ,µ , when no ambiguity is possible.

We investigate perturbations (λ + α, µ + β) over the Lamé coefficients of magnitude (i.e. of L ∞ -norm) m < γ, around a given reference state (λ, µ). These coefficients are confessedly not the physically relevant properties of the material over which perturbations should be considered: for instance, the impact of an increase in temperature on the material's properties is certainly better transcribed in terms of the Young's modulus and Poisson ratio. Nevertheless, we will focus on perturbations on the Lamé coefficients, so to keep expressions as light as possible; Young's modulus and Poisson ratio being analytical functions of the Lamé coefficients of the material, this last case would be no more difficult.

Let j : R d → R and k : R d → R be two functions of class C 2 , satisfying the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. The cost of the plate is defined by:

C(h, λ, µ) = Ω j(u h,λ,µ ) dx + Γ N k(u h,λ,µ ) ds.
The corresponding worst-case objective function is, in this context:

J (h) = sup α,β∈L ∞ (Ω) ||α|| L ∞ (Ω) ≤m ||β|| L ∞ (Ω) ≤m C(h, λ + α, µ + β).
By linearization, we approximate J by (42)

J (h) = C(h, λ, µ) + sup α,β∈L ∞ (Ω) ||α|| L ∞ (Ω) ≤m ||β|| L ∞ (Ω) ≤m ∂C ∂λ (h, λ, µ)(α) + ∂C ∂µ (h, λ, µ)(β) .
The next result gives an explicit formula for J , and a computation of its derive.

Theorem 8. The functional J defined by (42) reduces to:

(43) J (h) = Ω j(u h ) dx + Γ N k(u h ) ds + 2m||he(u h ) : e(u h )|| L 1 (Ω) +m||hdiv(u h )div(p h )|| L 1 (Ω) ,
where p h ∈ H 1 Γ D (Ω) d is the first adjoint state, defined as the unique solution to [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF]. Moreover, J is differentiable at any h ∈ U ad such that the set

E h := {x ∈ Ω, Ae(u h ) : e(p h ) = 0 or div(u h )div(p h ) = 0}
has zero Lebesgue measure and its derivative at such a point reads:

∀s ∈ L ∞ (Ω), J ′ (h)(s) = Ω s D(u h , p h , q h , z h ) dx,
where

D(u h , p h , q h , z h ) = Ae(u h ) : e(p h )+m 2ε e e(u h ) : e(p h )+ε d div(u h )div(p h )+Ae(p h ) : e(q h )+Ae(u h ) : e(z h ) ,
with ε e = sgn (e(u h ) : e(p h )), ε d = sgn (div(u h )div(p h )), the second and third adjoint states q h , z h ∈ H 1 Γ D (Ω) d being respectively defined as the unique solutions to some boundary value problems which are best formulated in terms of the associated variational equations:

∀v ∈ H 1 Γ D (Ω) d , Ω hAe(q h ) : e(v) dx = - Ω h 2ε e e(u h ) : e(v) + ε d div(u h )div(v) dx, ∀v ∈ H 1 Γ D (Ω) d , Ω hAe(z h ) : e(v) dx = - Ω h 2ε e e(p h ) : e(v) + ε d div(p h )div(v) dx - Ω ∇ 2 u j(u h )q h • v dx - Γ N ∇ 2 u k(u h )q h • v ds.
Proof. The proof is in essence identical to that of Theorem 4 and we refer to [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] for details, if necessary.

Remark 5. When the cost function C is the compliance of the plate, a result analogous to Proposition 5 holds true. Namely, one can compute exactly the worst-case functional J which is

J (h) = C(h, λ -m, µ -m).
The conclusions of Theorem 8 about the derivative of the approximate worst-case functional J can then be drawn closer to this exact expression for J . After some computations, we find, as in Remark 4, that for any

h ∈ U ad : ∀s ∈ L ∞ (Ω), J ′ (h)(s) = J ′ (h)(s) + O(m 2 ).
4. Worst-case design in shape optimization 4.1. Description of the model problem.

We are now interested in the design of shapes, that is bounded domains Ω ⊂ R d , with at least Lipschitz regularity. A shape Ω is submitted to body forces f ∈ H 1 (R d ) d and traction loads g ∈ H 2 (R) d on part of its boundary Γ N ⊂ ∂Ω, while it is clamped on another part Γ D ⊂ ∂Ω. Neither of these parts of the boundary is subject to optimization, and only the free boundary Γ := ∂Ω \ (Γ D ∪ Γ N ) is optimizable. The displacement of this shape is the unique solution in H 1 Γ D (Ω) d to the linear elasticity system posed in Ω:

(44)

       -div(Ae(u)) = f in Ω, u = 0 on Γ D , Ae(u)n = g on Γ N , Ae(u)n = 0 on Γ,
where A is the material's Hooke's law, defined by [START_REF] Chen | A new level-set based approach to shape and topology optimization under geometric uncertainty[END_REF]. In accordance with this setting, the set U ad of admissible domains is:

U ad = Ω ⊂ R d is open, Lipschitz and bounded, Γ D ∪ Γ N ⊂ ∂Ω .
As for representing variations of shapes, we rely on Hadamard's boundary variation method [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Soko | Introduction to shape optimization: shape sensitivity analysis[END_REF]: for a shape Ω ⊂ R d , we consider variations of the form:

Ω θ := (I + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1.
This brings about the following notion of differentiation with respect to the domain.

Definition 1. A functional J(Ω) of the domain is shape differentiable at Ω if the underlying function θ → J((I + θ)(Ω)), from W 1,∞ (R d , R d ) into R is Fréchet differentiable at θ = 0.
The shape derivative J ′ (Ω) of J at Ω is then the corresponding Fréchet differential, so that the following asymptotic expansion holds in the vicinity of 0

∈ W 1,∞ (R d , R d ): (45) J(Ω θ ) = J(Ω) + J ′ (Ω)(θ) + o(θ), where lim θ→0 |o(θ)| ||θ|| W 1,∞ (R d ,R d ) = 0.
To guarantee that all the considered variations of shapes belong to U ad , the set of deformations of shapes is restricted from

W 1,∞ (R d , R d ) to Θ ad ⊂ W 1,∞ (R d , R d )
, where:

Θ ad = θ ∈ W 1,∞ (R d , R d ) s.t. θ(x) = 0 for x ∈ Γ D ∪ Γ N .
Remark 6. Throughout Section 4, we shall consider shape derivatives of functionals which depend on the domain via the solution u Ω ∈ H 1 Γ D (Ω) d to the linear elasticity system (44), or similar adjoint systems. A rigorous justification of the shape differentiability of such functionals typically requires the shape differentiability of u Ω (or the featured adjoint states). This is a classical result (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF], [START_REF] Delfour | Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization[END_REF], [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF], [START_REF] Soko | Introduction to shape optimization: shape sensitivity analysis[END_REF]) which requires some higher smoothness of u Ω and of the adjoints. Owing to the regularity theory for the linear elasticity system [START_REF] Ciarlet | Mathematical Elasticity, vol I: Three Dimensional Elasticity[END_REF], it is indeed the case when the data at hand, namely Ω, f, g are smooth enough. We shall implicitly and systematically assume that these data are smooth although we do not give precise regularity assumptions.

4.2.

Worst-case design in shape optimization under uncertainties over the applied body forces.

In this section, we optimize the shape of a structure Ω with respect to a worst-case scenario when the body forces are subject to small perturbations. For the sake of simplicity, we limit ourselves to a single example of cost function, which is easily generalized to more complex situations (as well as to perturbations over the surface loads).

Let j : R d f × R d u → R be a function of class C 2 , satisfying the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. For any shape Ω ∈ U ad and any body force term f , let u Ω,f be the corresponding elastic displacement, solution to problem (44). The cost associated to this configuration is:

C(Ω, f ) = Ω j(f, u Ω,f ) dx.
We now fix a particular body force term f ∈ H 1 (R d ) d , and introduce the worst-case optimization problem at stake in the section: (46) min

Ω∈U ad J (Ω), where J (Ω) = sup ξ∈L 2 (R d ) d ||ξ|| L 2 (R d ) d ≤m C(Ω, f + ξ).
Note that the allowed perturbations ξ over the body force term are 'naturally' assumed to belong to L 2 (R d ) d , whereas the unperturbed body force term f is chosen in H 1 (R d ) d (which is a technical assumption when it comes to performing shape differentiation). As proposed in Section 2, the following linearized version of problem (46) lends itself to an easier analysis:

(47) min

Ω∈U ad J (Ω), where J (Ω) = sup ξ∈L 2 (R d ) d ||ξ|| L 2 (R d ) d ≤m C(Ω, f ) + ∂C ∂f (Ω, f )(ξ) .
The main result of this section is the following theorem which gives an explicit formula for J as well as its shape gradient.

Theorem 9. The functional J , defined in (47), rewrites:

(48) ∀Ω ∈ U ad , J (Ω) = Ω j(f, u Ω ) dx + m||∇ f j(f, u Ω ) -p Ω || L 2 (Ω) d ,
where p Ω ∈ H 1 Γ D (Ω) d is the first adjoint state, defined as the unique solution to

(49)    -div(Ae(p)) = -∇ u j(f, u Ω ) in Ω, p = 0 on Γ D , Ae(p)n = 0 on Γ ∪ Γ N .
Moreover, J is differentiable at any Ω ∈ U ad such that ∇ f j(f, u Ω ) -p Ω = 0 in L 2 (Ω) d and its derivative at such a point reads, for any θ ∈ Θ ad ,

J ′ (Ω)(θ) = Γ (j(f, u Ω ) + Ae(u Ω ) : e(p Ω ) -p Ω • f ) θ • n ds + m 2||∇ f j(f,uΩ)-pΩ|| L 2 (Ω) d Γ |∇ f j(f, u Ω ) -p Ω | 2 +∇ u j(f, u Ω ) • q Ω -z Ω • f θ • n ds + m 2||∇ f j(f,uΩ)-pΩ|| L 2 (Ω) d Γ Ae(u Ω ) : e(z Ω ) + Ae(p Ω ) : e(q Ω ) θ • n ds,
where the second and third adjoint states q Ω , z Ω ∈ H 1 Γ D (Ω) d are respectively defined as the unique solutions to the following systems:

(50)    -div(Ae(q)) = -2 (p Ω -∇ f j(f, u Ω )) in Ω, q = 0 on Γ D , Ae(q)n = 0 on Γ N , (51)
   -div(Ae(z)) = -∇ 2 u j(f, u Ω )q Ω -2∇ f ∇ u j(f, u Ω ) T (∇ f j(f, u Ω ) -p Ω ) in Ω, z = 0 on Γ D , Ae(z)n = 0 on Γ N .
Following the steps of Section 3, the proof of Theorem 9 requires the following lemma (which is an equivalent of Lemma 2 in the shape optimization context). As in Section 3, the dependence of the integrand function j on the f variable is omitted for the sake of simplicity.

Lemma 10. For any Ω ∈ U ad , denote by u Ω ∈ H 1 Γ D (Ω) d the unique solution to (44). (i) Let j : R d → R be a function of class C 1 and define a functional K of the domain Ω as:

∀Ω ∈ U ad , K(Ω) = Ω j(u Ω ) dx.
Then, K is shape differentiable at any Ω ∈ U ad , and its derivative reads, for any θ ∈ Θ ad ,

K ′ (Ω)(θ) = Γ j(u Ω ) + Ae(u Ω ) : e(p Ω ) -p Ω • f θ • n ds, where p Ω ∈ H 1 Γ D (Ω) d is
∀Ω ∈ U ad , L(Ω) = Ω ℓ(u Ω , p Ω ) dx, where p Ω ∈ H 1 Γ D (Ω) d is the unique solution to (52)    -div(Ae(p)) = -b(u Ω ) in Ω, p = 0 on Γ D , Ae(p)n = 0 on Γ ∪ Γ N .
Then, L is shape differentiable at any Ω ∈ U ad , and its shape derivative reads, for any θ ∈ Θ ad ,

L ′ (Ω)(θ) = Γ ℓ(u Ω , p Ω ) + b(u Ω ) • q Ω + Ae(u Ω ) : e(z Ω ) + Ae(p Ω ) : e(q Ω ) -z Ω • f θ • n ds,
where q Ω , z Ω ∈ H 1 Γ D (Ω) d are respectively defined as the unique solutions to the systems:

(53)    -div(Ae(q)) = -∇ p ℓ(u Ω , p Ω ) in Ω, q = 0 on Γ D , Ae(q)n = 0 on Γ ∪ Γ N , (54)
   -div(Ae(z)) = -∇ u b(u Ω ) • q Ω -∇ u ℓ(u Ω , p Ω ) in Ω, z = 0 on Γ D , Ae(z)n = 0 on Γ ∪ Γ N .
Proof. (i): This is a classical result in shape optimization (see e.g. Theorem 3.6 in [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]).

(ii): Let us once again rely on Céa's method. Introduce the Lagrangian L :

U ad ×H 1 Γ D (R d ) d ×H 1 Γ D (R d ) d → R, defined as: L(Ω, p, q) = Ω ℓ(u Ω , p) dx + Ω Ae( p) : e( q) dx + Ω b(u Ω ) • q dx,
and let us search for the points (p, q) ∈ H

1 Γ D (R d ) d × H 1 Γ D (R d ) d
where the partial derivatives of L(Ω, ., .) vanish, for a particular shape Ω ∈ U ad .

• The partial derivative of L with respect to q at (Ω, p, q) reads:

∀ q ∈ H 1 Γ D (R d ) d , ∂L ∂q (Ω, p, q)( q) = Ω Ae(p) : e( q) dx + Ω b(u Ω ) • q dx.
Canceling this expression yields precisely the variational formulation of system (52) which defines the 'adjoint' state p Ω . Thus, we deduce p = p Ω .

• The derivative of L with respect to p evaluated at (Ω, p, q) reads:

∀ p ∈ H 1 Γ D (R d ) d , ∂L ∂p (Ω, p, q)( p) = Ω ∇ p ℓ(u Ω , p) • p dx + Ω Ae( p) : e(q) dx.
Canceling this expression yields precisely the variational formulation of system (53) which defines the 'second adjoint' state q Ω . Thus, we deduce q = q Ω . • Eventually, one has L(Ω) = L(Ω, p Ω , q), for any Ω ∈ U ad and any q ∈ H 1 Γ D (R d ). As usual, differentiating this equality with respect to Ω (knowing that p Ω is shape differentiable according to Remark 6), then taking q = q Ω in the resulting expression yields:

∀θ ∈ Θ ad , L ′ (Ω)(θ) = L ′ (Ω, p Ω , q Ω )(θ).
This last (partial) shape derivative L ′ (Ω, p, q)(θ) can now be computed using point (i), since it depends on Ω only via u Ω , with the following choice of integrand j, for fixed functions p, q ∈ H 1 Γ D (R d ) d : j(u Ω ) = ℓ(u Ω , p) + Ae(p) : e(q) + b(u Ω ) • q, followed by the specification of p and q as p Ω and q Ω respectively. Introducing the third adjoint state z Ω ∈ H 1 Γ D (Ω) d as the unique solution to (54) (which is precisely the adjoint of point (i) for our above choice of j), we deduce the desired formula for L ′ (Ω)(θ).

We are now in a position to carry out the proof of Theorem 9.

Proof of Theorem 9. The proof follows the same steps as that of Theorem 1 and is divided into two steps.

• Expression of J as a function of Ω using an adjoint state.

Starting from formula (47), we first compute the derivative ∂C ∂f :

∀ξ ∈ L 2 (Ω) d , ∂C ∂f (Ω, f )(ξ) = Ω ∇ f j(f, u Ω ) • ξ + ∇ u j(f, u Ω ) • ∂u Ω,f ∂f (ξ) dx.
The first part ∇ f j(f, u Ω ) • ξ of the integrand does not pose any problem since it is explicitly linear in ξ. To simplify the second one, we use the adjoint system (49) for p Ω , the variational formulation of which is

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(p Ω ) : e(v) dx = - Ω ∇ u j(f, u Ω ) • v dx.
On the other hand, since the map f → u Ω,f is linear and continuous from L 2 (Ω) d into H 1 Γ D (Ω) d , it is obviously Fréchet differentiable. Differentiating with respect to f the variational formulation of (44) (giving the solution u Ω,f ) leads to

∀v ∈ H 1 Γ D (Ω) d , Ω Ae ∂u Ω,f ∂f (ξ) : e(v) dx = Ω ξ • v dx.
Combining those two relations yields

Ω ∇ u j(f, u Ω ) • ∂u Ω,f ∂f (ξ) dx = - Ω Ae(p Ω ) : e ∂u Ω,f ∂f (ξ) dx = - Ω ξ • p Ω dx.
A simple maximization with respect to ξ ∈ L 2 (Ω) d delivers the sought expression for J (Ω).

• Shape sensitivity analysis of J .

Differentiating the first term in (48) with respect to the domain is a straightforward application of point (i) in Lemma 10. Its shape derivative reads:

∀θ ∈ Θ ad , Ω j(f, u Ω ) dx ′ (θ) = Γ j(f, u Ω ) + Ae(u Ω ) : e(p Ω ) -p Ω • f θ • n ds.
As for the second term, we use point (ii) in Lemma 10 with ℓ(u, p) = |∇ f j(f, u) -p| 2 , so that

∇ u ℓ(u, p) = 2 ∇ f ∇ u j(u) T • (∇ f j(f, u) -p) , ∇ p ℓ(u, p) = 2(p -∇ f j(f, u)).
Doing so entails, for all θ ∈ Θ ad :

Ω |∇ f j(f, u Ω ) -p Ω | 2 dx ′ (θ) = Γ |∇ f j(f, u Ω ) -p Ω | 2 +∇ u j(f, u Ω ) • q Ω -z Ω • f θ • n ds + Γ Ae(u Ω ) : e(z Ω ) + Ae(p Ω ) : e(q Ω ) θ • n ds
where q Ω , z Ω ∈ H 1 Γ D (Ω) d are defined by ( 50) and (51). These two identities lead to the desired formula for J ′ (Ω)(θ).

Example 4. As an illustration of Theorem 9, assume that g = 0 and consider the case of the compliance as a cost function. We have: j(f, u) = f • u, therefore:

∇ f j(f, u) = u , ∇ u j(f, u) = f , ∇ f ∇ u j(f, u) = I , ∇ 2 u j(f, u) = 0.
In this context, it is easy to see that p Ω = -u Ω . Thus, J admits the following expression:

J (Ω) = Ω f • u Ω dx + 2m||u Ω || L 2 (Ω) d .
Furthermore, if u Ω = 0 in L 2 (Ω) d , J is shape differentiable at Ω. Easy computations allow us to deduce that the adjoint state q Ω ∈ H 1 Γ D (Ω) d is the unique solution to the system:

   -div(Ae(q)) = 4u Ω in Ω, q = 0 on Γ D , Ae(q)n = 0 on Γ ∪ Γ N ,
and that z Ω = -q Ω . The shape derivative of J at Ω then takes the form, for any θ ∈ Θ ad ,

J ′ (Ω)(θ) = Γ 2u Ω • f -Ae(u Ω ) : e(u Ω ) θ • n ds + m 4||uΩ|| L 2 (Ω) d Γ 4|u Ω | 2 +2q Ω • f -2Ae(u Ω ) : e(q Ω ) θ • n ds.

Worst-case design in shape optimization under uncertainties on the Lamé moduli of the material.

This section is intended as the mirror image of Section 3.4 in the context of shape optimization: we investigate the worst-case design of an objective functional with respect to perturbations on the Lamé coefficients of the elastic material filling Ω. The Lamé coefficients are functions λ, µ ∈ L ∞ (R d ) satisfying the uniform bounds in (41) (see Section 3.4). The solution to problem (44), posed on a shape Ω filled with a material having Lamé coefficients λ, µ, is denoted as u Ω,λ,µ .

Let j : R d → R be a function of class C 2 , satisfying the growth conditions [START_REF] Cherkaev | Principal Compliance and Robust Optimal Design[END_REF]. For any functions λ, µ ∈ L ∞ (R d ) satisfying (41), the cost of the structure Ω filled with a material with such Lamé coefficients is defined as:

C(Ω, λ, µ) = Ω j(u Ω,λ,µ ) dx.
Let us now fix reference values of the Lamé coefficients λ, µ. For the sake of simplicity, when the context is clear, we still denote u Ω = u Ω,λ,µ . Considering perturbations of magnitude m < γ over λ, µ, the corresponding worst-case objective function J is then:

∀Ω ∈ U ad , J (Ω) = sup α,β∈L ∞ (R d ) ||α|| L ∞ (R d ) ≤m ||β|| L ∞ (R d ) ≤m C(Ω, λ + α, µ + β),
and the approximated objective function J reads:

∀Ω ∈ U ad , J (Ω) = C(Ω, λ, µ) + sup α,β∈L ∞ (R d ) ||α|| L ∞ (R d ) ≤m ||β|| L ∞ (R d ) ≤m ∂C ∂λ (Ω, λ, µ)(α) + ∂C ∂µ (Ω, λ, µ)(β) .
We have the following result, whose proof (similar to that of Theorem 9) is omitted.

Theorem 11. The considered functional J rewrites, for any Ω ∈ U ad :

(55)

J (Ω) = Ω j(u Ω ) dx + 2m||e(u Ω ) : e(p Ω )|| L 1 (Ω) +m||div(u Ω )div(p Ω )|| L 1 (Ω) ,
where p Ω ∈ H 1 Γ D (Ω) d is the adjoint state, solution of (49). Moreover, J is differentiable at any shape Ω ∈ U ad such that the set

E Ω := {x ∈ Ω, e(u Ω ) : e(p Ω )(x) = 0 or div(u Ω )div(p Ω )(x) = 0}
has zero Lebesgue measure, and its shape derivative at such a point is:

∀θ ∈ Θ ad , J ′ (Ω)(θ) = Γ D(u h , p h , q h , z h ) θ • n ds,
where

D(u h , p h , q h , z h ) = j(u Ω ) + Ae(u Ω ) : e(p Ω ) -p Ω • f + m 2|e(u Ω ) : e(p Ω )|+|div(u Ω )div(p Ω )|+∇ u j(u Ω ) • q Ω + Ae(p Ω ) : e(q Ω ) + Ae(u Ω ) : e(z Ω ) ,
where q Ω , z Ω ∈ H 1 Γ D (Ω) d are respectively defined as the unique solutions to the variational problems:

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(q) : e(v) dx = - Ω 2ε e e(u Ω ) : e(v) + ε d div(u Ω )div(v) dx, ∀v ∈ H 1 Γ D (Ω) d , Ω Ae(z) : e(v) dx = - Ω 2ε e e(p Ω ) : e(v) + ε d div(p Ω )div(v) + ∇ 2 u j(u Ω )q Ω • v dx,
and ε e = sgn (e(u Ω ) : e(p Ω )), and ε d = sgn (div(u Ω )div(p Ω )).

Remark 7. The observations of Remark 5 on the compliance case can be straightforwardly extended from the parametric to the shape optimization setting (up to some extra computations).

4.4.

Worst-case design in shape optimization under geometric uncertainties.

The purpose of this section is to compute a shape derivative for a given functional of the domain J(Ω) which is robust with respect to uncertainties on the boundaries of shapes themselves. Let us first specify what we intend by shape optimization under geometric uncertainty. We assume that perturbations only affect the free boundary Γ of Ω. To that effect, we introduce a cutoff function χ :

R d → R and two open neighborhoods O 1 O 2 of Γ D ∪ Γ N in R d , enjoying the following properties: (56)
χ is smooth and nonnegative over R d ,

   χ(x) = 0 for x ∈ O 1 , χ(x) ∈ (0, 1] for x ∈ O 2 \ O 1 , χ(x) = 1 for x ∈ O 2 .
If m > 0 is the expected magnitude of perturbations over the geometry, we are interested in perturbations of Ω ∈ U ad of the form (see Figure 2):

(57) (I + χV )(Ω), V ∈ W 1,∞ (R d , R d ), ||V || L ∞ (R d ) d ≤ m.
Note that an other way to describe this problem consists in assuming perturbations of Ω of the form:

(58) (I + χvn)(Ω), v ∈ W 1,∞ (R d ), ||v|| L ∞ (R d ) ≤ m,
where n = n Ω denotes (an extension to R d of) the normal vector field to ∂Ω (the Ω -index is meant to emphasize its dependence on Ω and will be omitted when the situation is clear). As we shall observe, both descriptions are equivalent as far as we are concerned. However, the former one (57), which is the one retained in the following, proves more convenient from a mathematical viewpoint, since it features independent sets for admissible shapes and admissible perturbations of them; on a different note, perturbed shapes in the sense of (58) are 'less regular' than the unperturbed one. The chosen description for perturbations over the geometry slightly differs from that adopted in other contributions on the topic, for instance:

Ω V (I + χV )(Ω) Γ N Γ D
• in [START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF], [START_REF] Wang | On projection methods, convergence and robust formulations in topology optimization[END_REF], the only retained possibility is that Ω may suffer from a (small) uniform 'shrinking' or 'thickening', i.e. perturbations of shapes are of the form (58), with constant v; a filtering approach is used to incorporate this uncertainty into the objective function (this approach is generalized in [START_REF] Schevenels | Robust topology optimization accounting for spatially varying manufacturing errors[END_REF] to account for random perturbations). • Closer to the present work, in [START_REF] Chen | A new level-set based approach to shape and topology optimization under geometric uncertainty[END_REF], the authors also use perturbations of the form (58) with a scalar field v varying in a random fashion over the boundary of the shape, following a Gaussian probability distribution with 0 mean value.

• Eventually, in [START_REF] Guo | Robust structural topology optimization considering boundary uncertainties[END_REF], perturbations of a shape Ω are of the form (58), with v being bounded in L 2 (R d )norm, and 'small' in the sense that the discrepancy between the volumes of the perturbed and unperturbed shapes is 'small'.

The approaches considered in [START_REF] Guo | Robust structural topology optimization considering boundary uncertainties[END_REF][START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF] fit into our general framework, but we believe that the above setting is more suited to realistic uncertainties over the manufacturing process.

Let C(Ω) be a cost functional; the associated worst-case scenario functional is:

J (Ω) = sup V ∈W 1,∞ (R d ,R d ) ||V || L ∞ (R d ) d ≤m C((I + χV )(Ω)),
where χ is the cutoff function defined by ( 56).

In what follows, we will focus on several particular cases as regards the form of C(Ω); the presented techniques could easily be generalized to different problems. The first investigated example will be that of the compliance C(Ω) of a shape Ω:

(59) C(Ω) = Ω Ae(u Ω ) : e(u Ω ) dx = Ω f • u Ω dx + Γ N g • u Ω ds;
then we will turn to a functional J(Ω), which depends only on u Ω (not on e(u Ω )) by means of a smooth enough function j : R d → R (we have in mind the least-square discrepancy criterion, with respect to a target displacement):

(60)

J(Ω) = Ω j(u Ω ) dx.
Finally, we will consider the case of a functional S(Ω), which depends only on the stress tensor σ(u Ω ) := Ae(u Ω ), through a smooth function j : S(R d ) → R (we have in mind a L p -norm of the stress, as studied in [START_REF] Allaire | Minimum stress optimal design with the level set method[END_REF]):

(61) S(Ω) = Ω j(σ(u Ω )) dx,
where σ(u Ω ) := Ae(u Ω ) is the stress tensor associated to the displacement u Ω .

As in the corresponding context of parametric optimization (see Section 3.3), the case of the compliance as a cost function is especially simple, as confirmed by the following proposition, whose proof unrolls along the lines of Proposition 5. Proposition 12. Assume that no body forces are applied to the structures under optimization: f = 0, and that the cost function C(Ω) is the compliance, that is:

C(Ω) = C(Ω) = Ω Ae(u Ω ) : e(u Ω ) dx = Γ N g.u Ω ds.
Then, for any shape Ω ∈ U ad , the exact worst-case functional J reads:

J (Ω) = sup V ∈W 1,∞ (R d ,R d ) ||V || L ∞ (R d ) d ≤m C((I + χV )(Ω)) = C((I -mχn Ω )(Ω)).
Simply put, the most compliant shape among all the perturbed designs of a shape Ω according to (57) is the thinnest of all.

However, for a general objective function C(Ω), the worst-case functional J (Ω) is not explicit and we approximate it by J (Ω), defined as:

(62) ∀Ω ∈ U ad , J (Ω) = sup V ∈W 1,∞ (R d ,R d ) ||V || L ∞ (R d ) d ≤m (C(Ω) + C ′ (Ω)(χV )).
Before stating the results of interest, let us set some more notations. If Ω ∈ U ad is any shape, we denote as τ a local basis of tangent vectors to Γ, so that (τ, n) is a local orthonormal frame of R d . Any matrix M ∈ S(R d ) can be decomposed into this basis as:

M = M τ τ M τ n M nτ M nn
where M τ τ stands for the (d -1) × (d -1) tangential minor of M, M τ n is the vector of the (n -1) first tangential components of the normal column Mn, M nτ is the row vector of the (n -1) first tangential components of the normal row n T M, and M nn = Mn • n. We eventually denote as div Γ the tangential divergence operator defined on Γ.

Theorem 13. Consider the worst-case functional J (Ω), defined by (62), when the cost function is C(Ω) = J(Ω), defined by (60). Then J rewrites:

(63) J (Ω) = Ω j(u Ω ) dx + m Γ χ |j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω | ds,
where p Ω ∈ H 1 Γ D (Ω) d is the adjoint state, solution of (49). Furthermore, J is shape differentiable at any Ω ∈ U ad such that the set

E Γ := {x ∈ Γ, (j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω ) (x) = 0}
has zero (surface) Lebesgue measure. At such a point, its shape derivative reads, for any θ ∈ Θ ad , (64)

J ′ (Ω)(θ) = Γ χ j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω θ • n ds +m Γ ∂ ∂n + κ χ |j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω | θ • n ds +m Γ Ae(p Ω ) : e(q Ω ) + Ae(u Ω ) : e(z Ω ) -f • z Ω θ • n ds,
where κ is the mean curvature of ∂Ω and the second and third adjoint states q Ω , z Ω are respectively defined as the unique solutions in H 1 Γ D (Ω) d to the systems:

(65)        -div(Ae(q)) = 0 in Ω, q = 0 on Γ D , Ae(q)n = 0 on Γ N , Ae(q)n = εχ f + div Γ (εχ (σ(u Ω )) τ τ ) on Γ, and (66)
       -div(Ae(z)) = -∇ 2 u j(u Ω )q Ω in Ω, z = 0 on Γ D , Ae(z)n = 0 on Γ N , Ae(z)n = -εχ ∇ u j(u Ω ) + div Γ (εχ (σ(p Ω )) τ τ ) on Γ, where ε := sgn (j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω ).
Remark 9. As already said in Remark 6, computing shape derivatives require smoothness of the data Ω, f, g. It is clearly necessary to simply give a meaning to formula (63). More regularity is even required for the shape derivative (64). We shall not dwell on the precise assumptions and simply ask the data Ω, f, g to be smooth enough.

To prove Theorem 13 we need the following technical lemma. Lemma 14. (i) For any Ω ∈ U ad , the shape derivative of the functional J(Ω), defined by (60), reads:

∀θ ∈ Θ ad , J ′ (Ω)(θ) = Γ j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω θ • n ds, where p Ω ∈ H 1 Γ D (Ω) d
is the adjoint state, solution of (49).

(ii) Let ℓ : R d u × R d p × R e → R be a smooth enough function which vanishes in a neighborhood of Γ D ∪ Γ N , and define the functional L(Ω) as:

(67) L(Ω) = Γ ℓ(u Ω , p Ω , Ae(u Ω ) : e(p Ω )) ds,
where p Ω is defined by system (49). Then L is shape differentiable and its shape derivative reads:

(68) ∀θ ∈ Θ ad , L ′ (Ω)(θ) = Γ D(u Ω , p Ω , q Ω , z Ω ) θ • n ds, with (69)
D(u Ω , p Ω , q Ω , z Ω ) = ∂ ∂n + κ ℓ(u Ω , σ(u Ω ) τ τ : e(p Ω ) τ τ ) + Ae(p Ω ) : e(q Ω ) + Ae(u Ω ) : e(z Ω ) -f • z Ω ,
and the second an third adjoint states q Ω , z Ω ∈ H 1 Γ D (Ω) d are respectively defined as the unique solutions to:

(70)            -div(Ae(q)) = 0 in Ω q = 0 on Γ D Ae(q)n = 0 on Γ N Ae(q)n = -∇ p ℓ(u Ω , p Ω , Ae(u Ω ) : e(p Ω )) +div Γ ∂ℓ ∂e (u Ω , p Ω , Ae(u Ω ) : e(p Ω ))(σ(u Ω )) τ τ on Γ , and: (71)
           -div(Ae(z)) = -∇ 2 u j(u Ω )q Ω in Ω z = 0 on Γ D Ae(z)n = 0 on Γ N Ae(z)n = -∇ u ℓ(u Ω , p Ω , Ae(u Ω ) : e(p Ω )) +div Γ ∂ℓ ∂e (u Ω , p Ω , Ae(u Ω ) : e(p Ω ))(σ(p Ω )) τ τ on Γ .
Proof. (i): This is a very classical result in shape optimization (see e.g. [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]).

(ii): We need to assume that u Ω and p Ω enjoy more regularity than the sole 'natural'

H 1 (Ω) d regularity, e.g u Ω ∈ H 2 (Ω) d ∩ H 1 Γ D (Ω) d and p Ω ∈ H 2 (Ω) d ∩ H 1 Γ D (Ω) d
, so that the very definition of L(Ω) makes sense. To differentiate L(Ω), it proves convenient to introduce the Lagrangian functional L defined on

U ad × H 2 (R d ) d ∩ H 1 Γ D (R d ) d 4 by: L(Ω, u, z, p, q) = Γ ℓ( u, p, Ae( u) : e( p)) ds + Ω Ae( u) : e( z) dx - Ω f • z dx - Γ N g • z ds + Ω Ae( p) : e( q) dx - Ω -∇ u j( u) • q dx,
which incorporates as constraints both variational formulations for u Ω and p Ω , the respective solutions to (44) and (49). As usual, we look for points (u, z, p, q) ∈ H

2 (R d ) d ∩ H 1 Γ D (R d ) d 4
where the partial derivatives of L(Ω, ., .) vanish, for a particular shape Ω ∈ U ad .

• As before, canceling the partial derivative of L with respect to z at (Ω, u, z, p, q) imposes that u should satisfy:

∀ z ∈ H 2 (R d ) d ∩ H 1 Γ D (R d ) d , Ω Ae(u) : e( z) dx = Ω f • z dx + Γ N g • z ds. As H 2 (R d ) is dense in H 1 (R d ),
this is equivalent to the fact that u = u Ω .

• Similarly, canceling the derivative of L with respect to q at (Ω, u, z, p, q) imposes that p should satisfy:

∀ q ∈ H 2 (R d ) d ∩ H 1 Γ D (R d ) d , Ω Ae(p) : e( q) dx = - Ω ∇ u j(u Ω ) • q dx.
For the same reason, this implies that p = p Ω .

• Let us now study the partial derivative of L with respect to p at (Ω, u, z, p, q). It reads, for all p ∈ H , because ℓ ≡ 0 on Γ D ∪ Γ N . Under our standing assumption that u = u Ω is smooth enough, canceling this last expression yields a well-defined variational problem for q, which admits as unique solution in

2 (R d ) d ∩ H 1 Γ D (R d ) d : (72) ∂L ∂p (Ω, u, z, p, q)( p) = Γ ∇ p ℓ(u,
∀ p ∈ H 2 (R d ) d ∩ H 1 Γ D (R d ) d ,
H 2 (R d ) d ∩ H 1 Γ D (R d ) d
(owing to the regularity theory for linear elasticity, see [START_REF] Ciarlet | Mathematical Elasticity, vol I: Three Dimensional Elasticity[END_REF]) q = q Ω , defined by (70). • The study of the partial derivative of L with respect to u at (Ω, u, z, p, q) unrolls in the same way. It reads, for all u ∈ H

2 (R d ) d ∩ H 1 Γ D (R d ) d : ∂L ∂u
(Ω, u, z, p, q)( u) =

+ Ω Ae( u) : e(z) dx + Ω ∇ 2 u j(u)q • u dx ,
where the first equality holds because Ae(p)n = 0 on Γ, and the second one follows again from integration by parts on Γ, with ℓ ≡ 0 on ∂Γ. Under our assumption that u = u Ω is smooth enough, canceling this last expression yields a well-defined variational problem for z, which admits z = z Ω as unique solution in

H 2 (R d ) d ∩ H 1 Γ D (R d ) d
, defined by (71). Eventually, for any domain Ω ∈ U ad , and any fixed functions q, z ∈ H

2 (R d ) d ∩ H 1 Γ D (R d ) d , one has: L(Ω) = L(Ω, u Ω , z, p Ω , q),
whence, differentiating this expression with respect to Ω, and evaluating at q = q Ω and z = z Ω ,

∀θ ∈ Θ ad , L ′ (Ω)(θ) = L ′ (Ω, u Ω , z Ω , p Ω , q Ω )(θ),
and the desired formula (68) follows.

Proof of Theorem 13. First, using point (i) in Lemma 14, we know that J is shape differentiable. A maximizer 62) is then easily obtain: only its trace on Γ matters and it is simply given by V = m ε n on Γ, where ε := sgn (j(u Ω ) + Ae(u Ω ) : e(p Ω ) -f • p Ω ). This immediatly gives formula (63).

V ∈ W 1,∞ (R d , R d ) in formula (
Then, using point (ii) in Lemma 14, with ℓ(u, p, e) = χ |j(u) + e -f • p| produces formula (64).

Remark 10. Actually, the proposed approach in this paper almost requires to differentiate the considered cost function twice with respect to the domain. Hence, the formulae derived in this section could also be used in the context of a second-order algorithm for shape optimization (see [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] for more details, or [START_REF] Delfour | Anatomy of the Shape Hessian[END_REF] for further discussions around the shape Hessian).

We conclude this tour with the study of the stress-based cost function S(Ω) defined by (61). The following result is proved in the exact same way as Theorem 13, and the proof is therefore omitted.

Theorem 15. Consider the worst-case functional J (Ω), defined by (62), when the cost function is C(Ω) = S(Ω), defined by (61). Then J rewrites:

(73) J (Ω) = Ω j(σ(u Ω )) ds + m Γ χ |j(σ(u Ω )) + Ae(u Ω ) : e(p Ω ) -f • p Ω | ds,
where the adjoint state p Ω ∈ H 1 Γ D (Ω) d is defined as the unique solution to:

(74)    -div(Ae(p)) = div(A ∂j ∂σ (σ(u Ω ))) in Ω, p = 0 on Γ D , Ae(p)n = -A ∂j ∂σ (σ(u Ω ))n on Γ ∪ Γ N .
Furthermore, J is shape differentiable at any Ω ∈ U ad such that the set

E Γ := {x ∈ Γ, (j(σ(u Ω )) + Ae(u Ω ) : e(p Ω ) -f • p Ω ) (x) = 0}
has zero (surface) Lebesgue measure. At such a point, its shape derivative reads, for all θ ∈ Θ ad ,

J ′ (Ω)(θ) = Γ χ j(σ(u Ω )) + Ae(u Ω ) : e(p Ω ) -f • p Ω θ • n ds +m Γ ∂ ∂n + κ χ |j(σ(u Ω )) + Ae(u Ω ) : e(p Ω ) -f • p Ω | θ • n ds +m Γ Ae(p Ω ) : e(q Ω ) + Ae(u Ω ) : e(z Ω ) -f • z Ω θ • n ds , (75)
where the second and third adjoint states q Ω , z Ω are respectively defined as the unique solutions in H 1 Γ D (Ω) d to the following variational problems:

(76) ∀v ∈ H 1 Γ D (Ω) d , Ω Ae(q) : e(v) dx = - Γ εχ ((σ(u Ω )) τ τ : e(v) τ τ -f • v) ds, and (77)
∀v ∈ H 1 Γ D (Ω) d , Ω Ae(z) : e(v) dx = - Γ εχ A ∂j ∂σ (σ(u Ω )) τ τ : e(v) τ τ + (σ(p Ω )) τ τ : e(v) τ τ ds, - Ω ∂ 2 j ∂σ 2 (σ(u Ω ))Ae(q Ω ) : Ae(v) dx
where ε := sgn (j(σ(u Ω )) + Ae(u Ω ) : e(p Ω ) -f • p Ω ).

Numerical results

All computations are performed within the FreeFem++ environment [START_REF] Pironneau | FreeFem++ version 2[END_REF], except otherwise mentioned. We give approximate CPU times for the test cases which are ran on a MacBook Pro with a 2.66 GHz Intel Core 2 Duo processor. 5.1. Worst-case optimization problems in parametric structural optimization.

Let us start with the parametric structural optimization setting, and test the derivatives computed in Section 3. We reuse the notations introduced then: in every case, a cost functional C of the thickness (and of perturbation parameters) is considered, and the corresponding approximate worst-case functional J is minimized using either a steepest-descent algorithm, or an augmented Lagrangian algorithm (see [START_REF] Nocedal | Numerical optimization[END_REF], §17.4).

No efforts have been made to minimize the CPU time by using a better optimization algorithm. In both examples below, we impose bounds h min = 0.1 and h max = 1 over admissible thickness functions, and the initial design of the plate is described by a uniform thickness h = 0.5. The elastic material filling the plate is characterized by its (normalized) Young's modulus and Poisson ration ν given by: (78) E = 1, ν = 0.3.

5.1.1.

Uncertainties on the applied body forces in parametric optimization.

This first example illustrates the results of Section 3.2, and more precisely those of Example 2. The situation is as depicted in Figure 3: the plate of size 2 × 1 is clamped on its bottom-left and bottom-right sides, and its cost, when its thickness is h and when submitted to body forces f is its compliance:

C(h, f ) = Ω f • u h,f dx,
(no surface loads are applied). The plate is equipped with a triangular computational mesh having 10128 vertices (thus, twice as many triangles). The unperturbed state is associated to the following distribution of forces: f = (0, -1) near the centre of the bottom side of the plate (red spot on Figure 3 are located on the bottom side on the plate, between the regions where it is clamped and where the given force f is applied (grey areas on Figure 3, top). In such a case, only one adjoint is necessary to compute a gradient (see Theorem 4). The approximate worst-case functional J defined by ( 33) is considered for minimization, and so that the problem is not trivial, a volume constraint is added, using a fixed Lagrange multiplier ℓ = 5.10 -4 . The considered minimization problem thus becomes:

min h∈U ad J (h) + ℓ Ω h dx .
For increasing values of m, 100 iterations of a gradient-based steepest descent algorithm are performed and the resulting shapes and convergence histories are reported in Figures 3 and4.

Predictably, this simple setting does not really allow us to compare the resulting shapes for various values of the perturbation parameter m. Indeed, since the volume Lagrange multiplier is kept fixed for all value of m, the resulting optimal volume is increasing with m,.

To better capture the behavior of the optimal design in terms of the perturbation parameter m, we turn to a more 'realistic' context, where a volume constraint

Vol(h) :=

Ω h dx = V T ,
is enforced owing to an augmented Lagrangian method in the course of minimizing J . The same test case is run with a target volume V T = 0.7; 150 iterations prove necessary to achieve convergence of the algorithm, and each computation takes about 6 -7 minutes (except for the one associated to the value m = 0, which does not involve any computation of an adjoint state). Results and convergence histories are reported on Figures 5 and6. We now illustrate the proposed model for dealing with geometric uncertainty in parametric optimization of Section 3.3. Recall that, in this context, the case of the compliance as a cost function is almost trivial, as assessed by Proposition 5. Hence, we consider another example, depicted in Figure 7, which is the optimization of a gripping mechanism. The considered plate is fixed on a part Γ D ⊂ ∂Ω, and submitted to surface loads g ∈ L 2 (Γ N ) 2 on another part Γ N ⊂ ∂Ω. g equals (0, -1) on the upper part of Γ N , and (0, 1) on its lower part. The plate is endowed with a triangular mesh of 12382 vertices. In this case, we chose as cost function C(h) a least-square discrepancy criterion between the solution u h to [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF], and a target displacement u 0 , that is:

∀h ∈ U ad , C(h) = Γ T |u h -u 0 | 2 ds,
where Γ T is another non optimizable subset of ∂Ω, disjoint from Γ D and Γ N . We chose u 0 = (0, -1) on the upper part of Γ T , and u 0 = (0, 1) on its lower part.

To help the optimization algorithm of the associated approximate worst-case functional J to reach a 'light' optimal desing, a very small volume constraint is imposed by using a fixed Lagrange multiplier ℓ = 0.0003, which is a mere numerical token (the cost criterion does not vary monotonically with the volume of the structure in this case).

For increasing values of m, we perform 100 iterations of a steepest-descent algorithm which requires three adjoints, accoridng to the conclusions of Theorem 6. The optimal designs are depicted in Figures 7 and8. Each computation takes about 10 minutes (except for the one associated to m = 0, which only involves the computation of one adjoint state at each iteration of the process). As in the previous test-case, one observes that the performances of the obtained shapes in their unperturbed states worsen as m grows (which is coherent, since the larger the value of m, the lower the importance of this unperturbed problem in the balance expressed by J between perturbed and unperturbed states). Note also that the overall looks of the obtained shapes change significantly around the value m ≈ 0.0365, at which the optimization algorithm follows a different path towards a local minimum of J .

Details around the numerical implementation.

As far as numerical simulations are concerned, shape optimization of elastic structures differs from its parametric counterpart mainly regarding the difficulty to account for the evolutions of shapes during the process. To deal with this issue, we rely on the level set method, as was originally suggested in [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF], which roughly speaking consists in describing every shape Ω ⊂ R d by means of a scalar function φ : R d → R enjoying the properties:

∀x ∈ R d ,    φ(x) < 0 if x ∈ Ω φ(x) = 0 if x ∈ ∂Ω φ(x) > 0 if x ∈ c Ω .
The main asset of this change in perspectives lies in that the motion of a domain Ω(t), t ∈ [0, T ] evolving in time, driven by a normal velocity V (t, x) is translated in terms of a corresponding level set function φ(t, .) into the following Hamilton-Jacobi equation:

(79) ∂φ ∂t + V |∇φ|= 0 on [0, T ] × R d .
In the situation of Section 4 (whose notations are reused here), the minimization of a functional J (Ω) of the domain is considered, whose shape derivative is of the form:

∀θ ∈ Θ ad , J ′ (Ω)(θ) = Γ D(u Ω , p Ω , q Ω , z Ω ) θ • n ds,
for some algebraic combination D(u Ω , p Ω , q Ω , z Ω ) of solutions to (state and adjoints) elasticity systems posed on Ω. Then, V is set to 0 on Γ D ∪ Γ N and to -D(u Ω , p Ω , q Ω , z Ω ) on the free boundary Γ. In numerical practice, the whole space R d is reduced to a large working domain D ⊂ R d , which encloses all the considered shapes, and comes equipped with a fixed simplicial mesh T . The (state or adjoints) linear elasticity systems posed on a given shape Ω, involved in the expressions of V , cannot be computed exactly since Ω is only known by means of an associated level set function (i.e. no mesh of Ω is available). The Ersatz material approach [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF] is then used to transfer a linear elasticity system posed on Ω to one posed on D by filling the void D \ Ω with a very soft material with Hooke's law εA, ε ≪ 1 (ε = 1.e -3 in our examples).

While the elasticity analysis are performed in FreeFem++, the Hamilton-Jacobi equation ( 79) and the redistanciation of φ are solved using a characteristic method as described in [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF][START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF].

For the sake of completeness, the computational times of two representative computations are provided, in Sections 5.2.2 (whose model involves one adjoint state) and 5.2.4 (whose model involves three adjoint states). 5.2.2. Shape optimization under uncertainties about the applied loads.

∈ L 2 (Γ N ) d is its compliance: C(Ω, f ) = Ω f • u Ω,f dx + Γ N g • u Ω,f ds.
In a first example we consider the situation depicted in Figure 10: a mast is clamped on a part Γ D of its boundary and traction loads g = (0, -1) are applied on Γ N , near the bottom-left and bottom-right parts of its arms. In the unperturbed state, no body forces are applied (f = 0). Perturbations are expected as vertical body forces (0, ξ), of amplitude ||ξ|| L 2 (R d ) ≤ m, which are located on near the bottom of the arms of the mast (blue areas in Figure 10).

We minimize the corresponding worst-case scenario functional J with respect to the shape for different values of parameter m, using a fixed Lagrange multiplier ℓ = 1 to impose a volume constraint. 200 iterations of a gradient algorithm are used, and results are displayed on Figure 10.

One observes that, once again, the shapes tend to thicken as m grows, but also notices interesting changes in trends in the layout of the structure. Once again, to better appraise this feature, we run the very same example, using an augmented Lagrangian method to enforce a volume constraint Vol(Ω) = V T , where V T is a target volume (in this example, V T = 2000). Each computation (except for the one associated to m = 0) takes about 25 minutes, for a computational mesh composed of 11257 vertices. The results are reported on Figure 11, and confirm our initial guess (see also Figure 12 for convergence histories). A second example is the benchmark optimal bridge test case, as described in Figure 14 (top). The same computations (fixed Lagrange multiplier followed by fixed volume constraint) are applied on this example. The bridge is clamped on two sides of its boundary, and vertical body forces f = (0, -10) are applied at the middle of the bottom of the structure (yellow box). Vertical perturbations of amplitude lower than m are expected to occur on the blue areas. First, a minimization procedure is carried out, using a fixed Lagrange multiplier ℓ = 0.2 for the volume constraint, for several values of m, and results are to be seen on Figure 14. The results of the subsequent step, to get optimal shapes with the same target volume V T = 0.75 are displayed on Figure 15. In both cases, the computational mesh is worth 9205 vertices.

Remark 11. Some of the 'optimal' shapes displayed turn out to be non symmetric, whereas the setting of the corresponding test case is. This is mainly because no particular attention has been paid about this feature; in particular, the meshes of the (symmetric) bounding boxes are triangular, and not symmetric. The proposed approach is evaluated in the context of Section 4.3; we seek the optimal design of a force inverter: the considered shapes are clamped on the upper and lower parts of their left-hand side, a surface load g = (-0.1, 0) is applied at the centre of this left-hand side, and should exhibit a prescribed displacement u 0 = (1, 0) in a (non optimizable) area located at the centre of their right-hand side (see the details on Figure 16).

In this context, the cost of a shape Ω ∈ U ad , when filled with a material with Lamé coefficients λ, µ is:

C(Ω, λ, µ) = Ω k(x)|u Ω,λ,µ -u 0 | 2 dx,
where k is the characteristic function of the area where the target displacement should be reached. We are in search of a shape Ω ∈ U ad which minimizes this cost, when perturbations |α|≤ m, |β|≤ m are expected over the 'reference' Lamé coefficients λ, µ associated to (78), and this leads us to consider the functional J of formula (48).

As was the case in the example of Section 5.1.2, the performance of a shape has nothing to do with its weight. Nevertheless, for purely numerical purposes, we add a very small penalization, with respective parameters ℓ v = 5.e -3 and ℓ c = 0.02 on the volume and compliance of the shapes to the expression of J . The first additional term helps in removing the small 'islands' (i.e. disconnected parts obtained after topological (From left to right, top to bottom): optimal shape for the worst-case optimal bridge example, for m = 0, 0.2, 0.5, 1, 1.5, 2. The same volume constraint V = V T = 0.75 is imposed in all six cases. changes occurred), while the second one makes it easier to obtained a connected structure (which is difficult, since shapes tend to develop very small parts in the course of the process in order to gain flexibility).

Figure 16 shows the shapes obtained after 400 iterations of a gradient algorithm performed on a mesh of 8401 vertices, for several values of m, and Figure 17

  the adjoint state, solution of (49). (ii) Let b : R d → R d and ℓ : R d u × R d p → R be two functions of class C 1 . Note that ∇ u b is a d × d matrix, and that ∇ u ℓ and ∇ p ℓ are vectors in R d . Introduce the functional L, defined as:

Figure 2 .

 2 Figure 2. Perturbation (I + χV )(Ω) of a domain Ω

  p, Ae(u) : e(p)) • p + ∂ℓ ∂e (u, p, Ae(u) : e(p))Ae(u) : e( p) ds + Ω Ae( p) : e(q) dx .This last expression is not so convenient when it comes to achieving an unambiguous definition of q by means of a variational formulation over the spaceH 1 Γ D (R d ) d ,because the term p → Γ ∂ℓ ∂e (u, p, Ae(u) : e(p))Ae(u) : e( p) ds is not a continuous linear form over H 1 Γ D (R d ) d . However, we already identified u as u Ω . In particular, u = u Ω satisfies Neumann homogeneous boundary conditions Ae(u)n = 0 over Γ. This allows for a simplification of the nasty term in (72):

Γ∇

  u ℓ(u, p, Ae(u) : e(p)) • u + ∂ℓ ∂e (u, p, Ae(u) : e(p))(Ae(p)) τ τ : e( u) τ τ ds + Ω Ae( u) : e(z) dx + Ω ∇ 2 u j(u)q • u dx = Γ ∇ u ℓ(u, p, Ae(u) : e(p)) • u -div Γ ∂ℓ ∂e (u, p, Ae(u) : e(p))(Ae(p)) τ τ • u ds

Figure 3 .

 3 Figure 3. Minimization of the compliance with a fixed Lagrange multiplier ℓ = 5.10 -4 and uncertainties over the body forces; from left to right, top to bottom, m = 0 (with boundary conditions), 0.05, 0.1, 0.2, 0.35, 0.5, 1, 1.5, 2, 5.

Figure 4 .

 4 Figure4. Convergence histories for the problem of compliance minimization under uncertainties over the body forces, using the same Lagrange multiplier for all examples (cf. Figure3).

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Minimization of the compliance with uncertainties over the body forces, with an imposed volume V T = 0.7; from left to right, top to bottom, m = 0, 0.05, 0.1, 0.2, 0.35, 0.5, 1, 1.5, 2, 5, with a target volume V T = 0.7.

Figure 8 .

 8 Figure 8. Convergence histories for the parametric optimization example under uncertainties over the thickness (cf. Figure7).

Figure 9 .

 9 Figure 9. Optimization of the shape of a mast with respect to the worst-case compliance with perturbations under the body forces; setting of the test case.

Figure 10 .

 10 Figure 10. (From left to right): Optimal shapes of a mast for compliance minimization, with uncertainties on the body forces of amplitude m = 0, 0.6, 2, 3. The same Lagrange multiplier for the volume constraint is used in all cases.

Figure 11 .

 11 Figure 11. (From left to right): optimal shapes of a mast for compliance minimization, under uncertainties on the body forces of amplitude m = 0, 0.6, 2, 3; all the shapes have the some volume V T = 2000.

Figure 12 .Figure 13 .

 1213 Figure 12. Convergence history for approximate worst-case compliance (left) and for the volume (right) in the (worst-case) optimal mast test case.

Figure 14 .

 14 Figure 14. (From left to right, top to bottom): Optimal shape of a bridge under perturbations over the body forces of amplitude m = 0, 0.2, 0.5, 1, 1.5, 2. The same Lagrange multiplier for the volume constraint is used in all cases.

Figure 15 .

 15 Figure 15. (From left to right, top to bottom): optimal shape for the worst-case optimal bridge example, for m = 0, 0.2, 0.5, 1, 1.5, 2. The same volume constraint V = V T = 0.75 is imposed in all six cases.

  displays the corresponding deformed shapes. Note

Figure 16 .

 16 Figure 16. (From left to right, top to bottom): details of the test-case, optimal shape for the worst case force inverter test case, with perturbations over the Lamé coefficients of the material of magnitude m = 0, 0.001, 0.002, 0.003, 0.0045, 0.0075, 0.01, 0.02, 0.1.

Figure 17 .

 17 Figure 17. (From left to right, top to bottom): deformed configurations of the optimal shapes of Figure16, with m = 0, 0.001, 0.002, 0.003, 0.0045, 0.0075, 0.01, 0.02, 0.1. The bounding box of the optimal shapes is displayed in red.

1 Figure 18 .

 118 Figure 18. Convergence histories for the (worst-case) force inverter test case.

Figure 20 .

 20 Figure 20. (From left to right, top to bottom): deformed configurations of the optimal shapes for the worst-case gripping mechanism test case, with m = 0, 0.001, 0.002, 0.004, 0.005, 0.007, 0.009, 0.01, 0.02 (the bounding box of the optimal shapes is displayed in red).

Figure 23 .Figure 24 .

 2324 Figure 23. (From left to right, top to bottom): stress distribution (||σ|| 2 ) for m = 0, 0.005, 0.01, 0.015, 0.02, 0.05 in the optimal L-Beams displayed on Figure 22, using p = 2.

Figure 25 .

 25 Figure 25. (From left to right): optimal shape for m = 0, 0.01, 0.02, for the L-Beam example, using p = 5.

Figure 26 .Figure 27 .

 2627 Figure 26. (From left to right): stress distribution (||σ|| 2 ) for m = 0, 0.01, 0.02 for the L-Beam example, using p = 5.

  Ae(u) : e( p) = (Ae(u)) τ τ : e( p) τ τ a.e. on Γ. Using this information in (72) together with an integration by parts on Γ yields, for all p ∈ H 2 (R d ) d ∩

	H 1 Γ D (R d ) d :
	∂L ∂p	(Ω, u, z, p, q)( p) =	Γ	∇ p ℓ(u, p, Ae(u) : e(p)) • p -div Γ	∂ℓ ∂e	(u, p, Ae(u) : e(p))(Ae(u)) τ τ • p ds
			+

Ω

Ae( p) : e(q) dx

Acknowledgements. The authors wish to thank G. Michailidis (Centre de mathématiques appliquées, Ecole Polytechnique) for numerous fruitful discussions. Part of this work has been supported by the RODIN project (FUI AAP 13). G. A. is a member of the DEFI project at INRIA Saclay Ile-de-France.

that the two last ones are not 'inverting' the displacement of the shape. The convergence histories for these computations are reported on Figure 18.

Shape optimization under geometric uncertainties.

We eventually look into the setting of Section 4.4, where geometric uncertainties are considered, first in the context of the device of a gripping mechanism, as illustrated in Figure 19: the shapes of interest are clamped on the top and bottom parts of their left-hand side, and a small horizontal force g = (0.1, 0) is applied at the centre of this side, with the goal that the jaws (corresponding to the blue area in Figure 19) will close, namely comply with a target displacement u 0 , equal to (0, -0.2) on the upper part, and (0, 0.2) on the lower part. The cost C(Ω) of a shape Ω ∈ U ad reads:

where k is the characteristic function of the area near the jaws.

As perturbations of magnitude m on the geometry of shapes are expected, we aim at optimizing the approximate worst-case functional J associated to this problem, defined by formula (63). Small constraints over the volume and compliance of shapes are incorporated using fixed Lagrange multipliers ℓ v = 0.003, and ℓ c = 1, serving the same purposes as in the force inverter test case, and 200 iterations of the usual gradient-based algorithm are performed.

Several results are displayed on Figure 19, corresponding to different values of m. Each computation is performed on a computational mesh of 10288 vertices, and takes about 35 minutes, except for the one associated to m = 0. The corresponding displacements are shown on Figure 20, and it is easily seen that, as expected, the performances of the unperturbed shapes are less and less efficient as m increases, and the 'optimal shape' obtained for the value m = 0.1 is not functional, since the jaws do not tighten when the shape is deformed.

On a different note, enforcing robustness with respect to geometric uncertainties seems to urge the joints of the shape (i.e. the very thin parts which allow to invert the displacement caused by the applied surface loads g) to thicken. Intuitively, this goes in the same direction as imposing manufacturing constraints on the resulting shapes (see [START_REF] Michailidis | Manufacturing constraints and multi-phase shape and topology optimization via a level set method[END_REF] for a more extensive discussion on this issue).

The proposed approach for addressing geometric uncertainties is eventually applied to a case where the stress of structures is at stakes. Our shapes are now L-shaped beams, clamped on their upper part, and submitted to traction loads g = (0, -1) on a portion of their right-hand side (see the details on Figure 21). The cost C(Ω) of a shape Ω ∈ U ad is now related to the stress σ(u Ω ) induced by its displacement as:

where p ≥ 2, k is a characteristic function which equals 1 everywhere on the working domain except near the area Γ N where loads are applied, and ||.|| is the Frobenius norm for matrices. The worst-case design associated to this cost function is investigated, when uncertainties over the geometry of the shape of maximum amplitude m are expected. The approximate worst-case function J defined by (73) is minimized, using parameter p = 2.

An augmented Lagrangian method is used to impose a target volume V T = 0.8 on shapes, and all the computations take place on a mesh of 7296 vertices; see the results on Figure 22, the stress distribution in the resulting structures in Figure 23, and the convergence histories in Figure 24. Eventually, the same procedure is applied for the value p = 5 (and still increasing values for m). As expected, the resulting optimal shapes are more 'rounded' in the vicinity of the reentrant corner, where a stress singularity develops. See Figure 25 for results, Figure 26 for the stress distribution in the shapes, and Figure 27 for convergence histories.