M Moustapha

B Sudret

J.-M Bourinet

B Guillaume

Quantile-based optimization under uncertainties using adaptive Kriging surrogate models

Keywords:

Uncertainties are inherent to real-world systems. Taking them into account is crucial in industrial design problems and this might be achieved through reliability-based design optimization (RBDO) techniques. In this paper, we propose a quantile-based approach to solve RBDO problems. We first transform the safety constraints usually formulated as admissible probabilities of failure into constraints on quantiles of the performance criteria. In this formulation, the quantile level controls the degree of conservatism of the design. Starting with the premise that industrial applications often involve high-fidelity and time-consuming computational models, the proposed approach makes use of Kriging surrogate models (a.k.a. Gaussian process modeling).

Thanks to the Kriging variance (a measure of the local accuracy of the surrogate), we derive a procedure with two stages of enrichment of the design of computer experiments (DoE) used to construct the surrogate model. The first stage globally reduces the Kriging epistemic uncertainty and adds points in the vicinity of the limit-state surfaces describing the system performance to be attained. The second stage locally checks, and if necessary, improves the accuracy of the quantiles estimated along the optimization iterations. Applications to three analytical examples and to the optimal design of a car body subsystem (minimal mass under mechanical safety constraints) show the accuracy and the remarkable efficiency brought by the proposed procedure.

Introduction

In engineering design, one often seeks to lower the product cost while ensuring its integrity.

These are by construction two conflicting objectives. Optimization has therefore been used as an automatic procedure to find a good trade-off. The optimal solution usually lies at the boundary of the feasible space. However uncertainties are ubiquitous to engineering systems whether arising from modeling approximations or input parameters inherent variability. They make any optimal design solution likely to depart from its real-world counterpart. Such discrepancy may turn a feasible solution into an unfeasible one. It is therefore of prime importance to account for uncertainties during optimization. This is generally achieved through robust and reliabilitybased design optimization (respectively RDO and RBDO). In the former, emphasis is put on the cost function. The designer actually searches for a design that is immune to the inputs uncertainties. The cost function is in this case replaced by robustness measures which include worst-case scenarios or moment-based criteria [START_REF] Trosset | Taguchi and robust optimization[END_REF]. [START_REF] Beyer | Robust optimization: A comprehensive survey[END_REF] and [START_REF] Baudoui | Optimisation robuste multiobjectifs par modèles de substitution[END_REF] give a comprehensive review of such techniques. On the other hand, reliabilitybased design optimization rather seeks to balance the cost and the safety requirements by moving the solution away from the boundary of the admissible space. The work presented in this paper is concerned with the latter approach.

Following the notations in [START_REF] Dubourg | Reliability-based design optimization using Kriging and subset simulation[END_REF]

where a cost function c is minimized with respect to design variables d. This minimization task is carried out under a set of constraints divided into two groups respectively denoted by soft and hard constraints. The n s soft constraints f j are simple analytical functions, often bounding the design space while the n h hard constraints g k are actually the system performance functions.

They rely on the mechanical model M k used to predict the structural behavior. In our case, they result from a finite element model and may be written as g k = ḡk -M k , where ḡk is a threshold not to be exceeded by the structural response which is computed from a simulation model x → M k (x) (usually a time-consuming finite element model). When safety requirements are of interest, performance may be measured in terms of a failure probability. To this end, random variables accounting for the uncertainties in the inputs are introduced and denoted respectively by X ∼ f X (x) for the design variables and Z ∼ f Z for the environmental variables.

The former notation means that the distribution of X is conditioned on the design parameters.

Typically, design parameters d are nominal dimension and f X (x) models the uncertainties due to manufacturing tolerances. Environmental variables Z may for instance be parameters of the crash protocol such as the impact speed in crashworthiness design. By propagating these uncertainties to the output, the failure probability for a given design d reads:

P f k (d) = P (g k (W ) ≤ 0) = g k (w)≤0 f W (w) dw, (2)
where W = {X|d, Z} T ∼ f W is a vector gathering all the input parameters of the mechanical model that governs the structure's behavior. The solution of this integral is generally not tractable because the failure domain defined by {w : g k (w) ≤ 0} has an implicit definition. One rather resorts to approximation or simulation methods [START_REF] Madsen | Methods of structural safety[END_REF]. Within the former group, the first-order reliability method (FORM) is the most widely used [START_REF] Ditlevsen | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF][START_REF] Lemaire | Structural reliability[END_REF][START_REF] Hasofer | Exact and invariant second-moment code format[END_REF]. It consists in mapping the random variables into the standard normal space where the limit-state surface is linearly approximated. Therefore, the failure probability can be equivalently expressed by the so-called reliability index [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF]. Curvatures of the limit-state surface may be handled by the second-order reliability method (SORM). As for the simulation methods, the most straightforward one is crude Monte Carlo simulation (MCS) where the failure probability is estimated by the relative occurrence of failed samples. The accuracy of the estimate depends on the number of samples. For extremely small probabilities of failure, the required number of samples for an accurate estimate becomes relatively high, typically 10 6-8 , which makes the approach not affordable. Variance-reduction techniques have been introduced in order to by-pass this limitation [START_REF] Asmussen | Stochastic simulation: algorithms and analysis[END_REF].

Applied to rare events simulations, such techniques include importance sampling [START_REF] Au | A new adaptive importance sampling scheme for reliability calculations[END_REF][START_REF] Melchers | Importance sampling in structural systems[END_REF] and subset simulation [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. The former proceeds by sampling from an instrumental distribution which puts a higher weight to the failure domain and afterwards correct the introduced bias appropriately. The latter splits the failure domain into nested auxiliary domains such that the failure probability can be estimated by the product of larger ones, the latter being easier to evaluate by simulation. However, the computational cost of all these techniques is in the order of 10 3-4 (for each design d), and can thus not be used within an optimization loop.

Indeed, the solution of the RBDO problem relies on the estimate of the failure probability for different values of the design parameters. Many techniques exist and may be classified into two-level, mono-level and decoupled approaches [START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF][START_REF] Aoues | Benchmark study of numerical methods for reliabilitybased design optimization[END_REF]. Two-level approaches, which basically consist of two nested loops, are among the most straightforward to implement. The outer loop explores the design space and the inner one solves the reliability analysis for any given design. Usually, the inner loop resorts to FORM approximations as in the so-called reliability index [START_REF] Enevoldsen | Reliability-based optimization in structural engineering[END_REF] and performance measure [START_REF] Tu | A performance measure approach in reliability-based structural optimization[END_REF][START_REF] Tu | A new study on reliability-based design optimization[END_REF] approaches (respectively RIA and PMA). Simulation techniques may also be used in the inner loop, as we show in the sequel. The mono-level approach transforms the double-loop problem into a single-loop one by introducing optimality criteria for the FORM problem. [START_REF] Kuschel | Two basic problems in reliability-based structural optimization[END_REF] propose an equivalent formulation based on RIA while [START_REF] Agarwal | An inverse-measure-based unilevel architecture for reliability-based design optimization[END_REF] rather rely on PMA. Finally, the decoupled approaches transform the double-loop into a sequence of deterministic problems. A well-known example is sequential optimization and reliability analysis (SORA) proposed in Du and Chen (2004).

As introduced above, all theses methods rely on repeated evaluations of the mechanical models, i.e. during the outer optimization loop and more intensively in the reliability analysis steps. This limits their range of applications to engineering problems of practical interest. This issue is even more dramatic in the design of complex industrial systems which relies on highfidelity models and henceforth time-consuming simulations. Surrogate modeling, a technique in which the mechanical model is replaced by a well calibrated easy-to-evaluate analytical function, has been extensively used in the past decade to alleviate the computational burden. For instance, support vector machines have been used for structural reliability assessment in [START_REF] Hurtado | Neural-network-based reliability analysis: a comparative study[END_REF]; [START_REF] Bourinet | Assessing small failure probabilities by combined subset simulation and support vector machines[END_REF]; [START_REF] Deheeger | Support vector machine for efficient subset simulations: 2 SMART method[END_REF]. Polynomial chaos expansion were considered in [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]; [START_REF] Blatman | Reliability analysis of a pressurized water reactor vessel using sparse polynomial chaos expansions[END_REF]; [START_REF] Hu | Adaptive-sparse polynomial chaos expansion for reliability analysis and design of complex engineering systems[END_REF]. Kriging proposed some approaches which rely on locally or globally refined Kriging approximations.

Likewise, the present work considers Kriging because it provides not only an approximation of the original mechanical model but also gives a built-in error estimate. This enables adaptive techniques that further reduce the computational cost.

In contrast to most of the literature in RBDO, we are not interested in this paper in highly reliable designs, for which the probability of failure has to be computed by one of the methods mentioned above. Our goal is rather to develop a conservative optimal design methodology. To this aim, we first introduce a quantile-based design optimization procedure, in which the hard constraints are formulated on quantiles of the performance criteria, instead of target probabilities of failure. At each iteration of the design optimization, we therefore evaluate quantiles through Monte Carlo simulation instead of solving a reliability problem. This approach is justified by the "degree of conservatism" targeted in our applications to car body mass optimization: 95%-quantiles are indeed considered as sufficiently conservative, which makes their evaluation relatively easy. As the performance criteria are obtained from time-consuming simulation (e.g.

frontal impact of a full car body or a subsystem), quantile evaluation must rely on surrogate models. In this respect, the goal of the paper is to propose a quantile-based design optimization methodology which is relying on adaptive Kriging surrogate models.

The paper is organized as follows: in Section 2 we introduce the quantile-based optimization and prove its formal equivalence with the classical RBDO setting. The basics of Kriging are then summarized in Section 3. In Section 4, an original two-stage strategy of enrichment of the experimental designs used in Kriging is proposed, as a means to reduce the overall computational burden (e.g. , at most a few hundreds runs of the time-consuming computational model) to regions of the design space that are relevant for optimization. Finally, Section 5 presents four examples: the three first involve analytical constraints and allows us to validate our approach against benchmark results. The final example is a real case study which deals with the mass optimization of a car body subsystem under crashworthiness constraints.

2 Formulation of the quantile-based optimization procedure

Equivalence between RBDO and quantile-based formulation

Prior to formulating the quantile-based procedure, let us consider the reliability-based design optimization problem in Eq. (1). By explicitly introducing the computational model of interest M which describes the system performance, the following equivalence holds:

P (g (X (d) , Z) ≤ 0) ≤ Pf ⇔ P (M (X (d) , Z) ≥ ḡ) ≤ Pf , ⇔ P (M (X (d) , Z) ≤ ḡ) ≥ 1 -Pf , (3)
where ḡ is an upper threshold on the system mechanical response.

From the last expression, we can introduce the following quantile as an alternative way of measuring the failure probability:

Q α (d; M (X (d) , Z)) = inf {q ∈ R : P (M (X (d) , Z) ≤ q) ≥ α} , (4)
where α = 1 -Pf .

The computed quantile may henceforth be used as a measure of reliability given a target failure probability. Considering Eq. (3) and Eq. ( 4), the following equivalence holds:

P (g (X (d) , Z) ≤ 0) ≤ Pf ⇔ Q α (d; M (X (d) , Z)) ≤ ḡ, (5)
where the value of α is directly related to the target failure probability.

This equivalence between the failure probability and the quantile estimation is illustrated in Figure 1 where the distributions of a mechanical response in two configurations are shown. In the upper panel, the quantile corresponding to the target failure probability is lower than the constraint threshold. This corresponds to a safe design since the probability that M (X (d) , Z)

is greater that ḡ is smaller than Pf = 1α. In contrast, the lower panel shows an unsafe design since the associated quantile is higher than the threshold ḡ, meaning that the probability of failure P (M (X (d) , Z) ≥ ḡ) is greater than Pf .

Following the previous developments, the RBDO problem of Eq. ( 1) may eventually be recast as:

d * = arg min d∈D c (d) subject to:    f j (d) ≤ 0, {j = 1, . . . , n s } , Q α k (d; M k (X (d) , Z)) ≤ ḡk , {k = 1, . . . , n h } , (6)
where

α k = 1 -Pf k .

Monte Carlo estimate of the quantile

To solve the optimization problem in Eq. ( 6), the quantile must be estimated in each iteration for the current design d (i) . In this paper, we consider crude Monte Carlo sampling. The following steps describe the numerical procedure:

1. Sample the Monte Carlo set needed to evaluate the quantile: j) , z (j) , j = 1, . . . , N ,

C q d (i) = x (
where X ∼ f X|d (i) , Z ∼ f Z and N is the size of the Monte Carlo sample set.

2. Compute the set of associated responses for each mechanical model:

Y k = y (j) k = M k x (j) , z (j) , j = 1, . . . , N (8)
3. Sort them in ascending order such that

y k (1) ≤ y k (2) ≤ . . . ≤ y k (N )
4. Retrieve the quantile corresponding to the k-th constraint by:

Q α k d (i) ; M k X d (i) , Z ≡ q α k d (i) = y k ( N α k ) , (9)
where t denotes the floor function yielding the largest integer smaller than t.

To apply this approach, the Monte Carlo sample set in Eq. ( 7) needs to be large enough so that the computed quantile is accurate. For our application, where the target failure probability is 1% -10%, we choose N = 10, 000. As this simulation is embedded in the iterative process of optimization, the number of calls to the mechanical model may reach hundreds of thousands.

Pf = 1 -α PDF of M(X(d), Z) M(x , z) q 0.5 ḡ q α (a) Safe design: P (M (X (d) , Z) ≥ ḡ) ≤ Pf ⇔ Qα (d; M (X (d) , Z)) ≤ ḡ Pf = 1 -α PDF of M(X(d), Z) M(x , z) q 0.5 ḡ q α (b) Unsafe design: P (M (X (d) , Z) ≥ ḡ) > Pf ⇔ Qα (d; M (X (d) , Z)) > ḡ
When a high-fidelity model is involved, such a large number of calls is not affordable. We therefore couple the proposed approach to a well-known surrogate modeling technique, namely Kriging.

3 Kriging (a.k.a Gaussian process modeling)

Surrogate models have been increasingly used as proxies of time-consuming functions in the past decade. In the computer experiments setting, such a function is considered to be a black-box i.e. only pairs of inputs/outputs are known with respect to a limited set of observations. This set constitutes the design of experiments and reads, for a given model M:

D = {(x i , y i ) , x i ∈ R s , y i = M (x i ) , i = 1, . . . , n} , (10)
where x i is an s-dimensional input, y i is the corresponding scalar output and n is the number of available observations in the design of experiments.

Kriging a.k.a. Gaussian process modeling [START_REF] Santner | The design and analysis of computer experiments[END_REF] is one particular emulator which considers the function M to approximate as a realization of a stochastic process, which may be cast as:

M (x) = p j=1 β j f j (x) + Z (x) , (11)
where the first summand is the deterministic part referred to as the trend. It reads as a linear combination of a vector of p weight coefficients β = {β j , j = 1, . . . , p} and a set of function basis f = {f j , j = 1, . . . , p}. The second summand is a zero-mean stationary Gaussian process. It is completely defined by its auto-covariance function Cov [Z (x) , Z (x )] = σ 2 R (x, x ; θ), where σ 2 is the constant variance of the Gaussian process, R is the auto-correlation function whose hyperparameters are gathered in the vector θ.

The calibration of the Kriging model involves making a few choices that can be motivated by some prior knowledge on the function to approximate. The first one is the choice of the mean trend. In this work, we consider an unknown constant trend. This results in the socalled ordinary Kriging. The second one is the choice of the auto-correlation function which encodes assumptions such as the degree of regularity of the underlying process. A wide family of auto-correlation functions have been used in the literature. Here, we consider the Matérn 5/2 auto-correlation family, defined in the one-dimensional case by:

R (x, x ; l) = 1 + √ 5 |x -x | l + 5 3 (x -x ) 2 l 2 exp - √ 5 |x -x | l , ( 12
)
where l is the so-called characteristic length scale. The multi-dimensional case is obtained by tensor product of the above equation:

R (x, x ; θ) = s i=1 R x i , x i ; θ i . (13)
Once these choices are made, the Kriging predictor at the point x is assumed to follow a

normal distribution M (x) ∼ N µ M (x) , σ 2 M (x) : µ M (x) = f T (x) β + r T (x) R -1 y -F T β , σ 2 M (x) = σ 2 1 -r T (x) R -1 r (x) + u T (x) F T R -1 F -1 u (x) , (14)
where

β = F T R -1 F -1
F T R -1 y is the generalized least-square estimate of the weight coefficients β, r (x) is a vector of cross-correlations between the point x and each point of the design of experiments, F is the information matrix whose components are f j (x i ) , i = {1, . . . , n} , j = {1, . . . , p} and u = F T R -1 r (x)f (x) has been introduced for the sake of clarity. Beside the prediction given by µ M (x), Kriging features a measure its own accuracy through the prediction variance σ 2 M (x). Confidence intervals on the prediction can then be derived since the distribution of the prediction M (x) is Gaussian by assumption. More importantly, this has supported the development of infill-sampling criteria used for the adaptive refinement of Kriging models.

Eventually, one has to estimate the hyperparameters of the auto-correlation functions to completely define the Kriging predictor. This is achieved through automatic calibration following techniques such as cross-validation or maximum likelihood estimation. The latter is used in this work and boils down to the following optimization problem:

θ = arg min θ∈R d ψ (θ) = σ 2 (θ) det R (θ) 1 n , (15)
where ψ (θ) is the so-called reduced likelihood function and d is the number of parameters in θ [START_REF] Koehler | Computer experiments[END_REF][START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF].

The accuracy of the solution of the optimization problem in Eq. ( 15) is crucial as it conditions the quality of the Kriging predictor. General-purpose algorithms such as genetic algorithm or BFGS are often used. Available softwares such as DiceKriging [START_REF] Roustant | DiceKriging, DiceOptim: two R packages for the analysis of computer experiments by Kriging-based metamodeling and optimization[END_REF] in R, UQLab [START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF][START_REF] Lataniotis | UQLab user manual Kriging[END_REF] or ooDace [START_REF] Couckyut | ooDace toolbox A Matlab Kriging toolbox: Getting started[END_REF] in MATLAB make use of such algorithms and more generally provide a framework to train a Kriging model.

Kriging-based optimization 4.1 Construction of a Kriging model in the augmented reliability space

In this section, the optimization problem in Eq. ( 6) is solved while each performance function M k is replaced by a Kriging model M k as introduced above. This simply means that the performance functions for now on read:

g k (x, z) = ḡk -M k (x, z) , k = {1, . . . , n h } . ( 16
)
Computing the quantile with respect to this surrogate model becomes an extremely cheap operation. The expensive part is the initial building of the surrogate model which requires to set and evaluate a design of experiments. Given the possible number of iterations before convergence is achieved, building one surrogate model for each reliability analysis (i.e. the quantile computation for a given design in our case) would be quite cumbersome. Instead we advocate for the use of a single Kriging model as already proposed in other contributions. The idea is to build the surrogate model in a unique space that embeds both the design and random variables. In [START_REF] Kharmanda | Efficient reliability-based design optimization using a hybrid space with application to finite element analysis[END_REF], this space is called the hybrid design space and is defined as the tensor product between the design and random variables. This needlessly increases the dimension of the space where the surrogate model is built and may be problematic when it comes to space-filling design of experiments. From another perspective, [START_REF] Au | Reliability-based design sensitivity by efficient simulation[END_REF] efficiently computes the failure probability in the so-called augmented reliability problem for a given design considering a space where the design variables are artificially considered as random. [START_REF] Taflanidis | Stochastic subset optimization for optimal reliability problems[END_REF] use this augmented reliability problem to construct a stochastic optimization problem. Eventually, [START_REF] Dubourg | Reliability-based design optimization using Kriging and subset simulation[END_REF]; [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF] propose an augmented reliability space following the ideas in the two above contributions. In these works, the size of the augmented reliability space remains equal to that of the original reliability problem. This is because they consider that the uncertainty in the random variables is simply augmented by the choice of the design points. In practice, the design and environmental variables are treated separately in the so-called confidence regions which span a sufficiently large space such that any point sampled during the analysis is extremely likely to fall within the space of definition of the surrogate model. For the design variables, this region is hyper-rectangular and consists of the design space with extended bounds. The confidence region for the environmental variables is a hypersphere in the standard normal space (i.e. , after transforming these variables into standard Gaussian ones) with a sufficiently large radius to account for extreme realizations of the random variables.

The augmented space is henceforth considered as the tensor product between these two spaces.

In this paper, we propose an augmented space which is quite close to that defined in [START_REF] Dubourg | Reliability-based design optimization using Kriging and subset simulation[END_REF]. We indeed treat separately the design and environmental variables. However, in our case, the environmental variables are not defined in a hypersphere. The reasons for this are twofold. First, due to its very formulation, the reliability analysis we perform does not need any mapping to the standard Gaussian space. Second, the non-linear mapping from the unit hypersphere (where the space-filling design of experiments is sampled) to the physical space may add complexity and non linearity to the function that is eventually surrogated. To avoid this, we rather consider a hypercube. Since the surrogate models are built in the unit hypercube, the mapping to the physical space is simply linear.

The augmented space is therefore the tensor product between two hyperrectangular confidence regions X×Z, where X refers to the design variables and Z to the environmental parameters.

The former is defined by:

X = s d i=1 q - di , q + di , (17)
where s d is the number of design variables and q - di and q + di are respectively quantiles associated to the lower and upper bounds of the design variables. They are defined in such a way that the confidence region spans a space sufficiently large to contain with high probability (e.g. 99-99.9% in the application) all realizations of X|d sampled during the optimization procedure. They read as follows:

q - di = F -1 Xi|di (α di /2) q + di = F -1 Xi|di (1 -α di /2) , (18)
where X i follows the marginal distribution f Xi|di , F -1 Xi|di is the associated inverse CDF, d - i and d + i are respectively the lower and upper bounds of the design variable d i , and α di is the probability of sampling outside the augmented space. In applications we select α di = 2.7 10 -3 for each variable, which corresponds to µ ± 3σ for a Gaussian variable.

In the same fashion, the confidence region for the environmental variables is defined by:

Z = sz j=1 q - zj , q + zj , (19)
where s z is the number of environmental variables and the bounding quantiles are defined by:

q - zj = F -1 Zj α zj /2 , q + zj = F -1 Zj 1 -α zj /2 , ( 20
)
where Z j follows the marginal distribution f Zj whose inverse CDF is F -1 Zj and α zj is the probability of sampling outside the augmented space in the direction of Z j , again in the order of 10 -3 in applications.

To illustrate the augmented space defined in this paper, we consider a problem where the design space is one-dimensional: D = [d -, d + ]. For the RBDO problem, the design variable is supposed random with distribution d ∼ N d, σ 2 d . We also assume that the RBDO problem features a unique environmental random variable defined by Z ∼ N µ z , σ 2

z . An augmented space associated to this problem is shown in Figure 2. The design space is the blue line and the augmented space is the gray area. The distributions of the design and environmental variables are also plotted.

z x d - d + q - z q + z µ z q - d q + d f X|d - f X|d + f Z α z /2 α z /2 X × Z D Figure 2:
Illustration of the augmented space for a two-dimensional problem with both random design and environmental variables.

Note that the proposed framework naturally encompasses the two following cases:

• When the analysts disregard uncertainties in the design parameters (e.g. , ignore manufacturing tolerances), q d -and q d + are simply set to d -and d + respectively.

• When no environmental variables are considered, the augmented space reduces to X.

Once the augmented space is defined for a specific problem, one may build a single global surrogate model. This may be achieved by a space-filling design of experiments, i.e. by using uniformly distributed samples so as to cover the entire space. The built Kriging model may henceforth be used for any reliability analysis during the RBDO. This one-shot approach is theoretically possible but would require the Kriging model to be accurate in the entire space.

However, during the optimization only a subset of the space is actually of interest, i.e. regions in the vicinity of the limit-state surface and those where the objective function decreases. These two issues can be dealt with using so-called adaptive design of experiments. We propose in this paper a two-stage enrichment scheme where each stage is geared toward achieving one of the two above goals.

Adaptive design of experiments: A short literature review

Adaptive design of experiments have been developed from the premise that only a limited region of the space is of interest to the designer during an optimization analysis. Thus, instead of densely filling the space so as to have an evenly accurate model in the entire space, the optimization starts with a not so accurate model built upon a scarcely sampled design of experiments.

Enrichment is then made so as to improve the surrogate accuracy in regions that matter. [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] proposed the efficient global optimization (EGO) scheme relying on an expected improvement function which focuses on sequentially updating a Kriging model so as to converge to a global minimum. From the same idea, numerous authors have proposed infill sampling criteria to achieve the same goal. As for RBDO, emphasis is rather put on the vicinity of the limit-state surface in order to accurately estimate failure probabilities. A first family of infill sampling criteria comes from EGO techniques as they are mere adaptation, e.g. adjusted expected improvement [START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF], expected violation [START_REF] Audet | A surrogate-modelbased method for constrained optimization[END_REF] or expected improvement for contour approximation [START_REF] Ranjan | Sequential experiment design for contour estimation from complex computer codes[END_REF]. On the other hand, Bichon et al.

(2008) introduced a so-called efficient global reliability analysis (EGRA) where an expected feasibility function is used to improve the surrogate model in the vicinity of the limit-state surface.

Similarly, [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF] used in his PhD thesis work the margin probability function. In the present work, we will focus on the deviation number developed by [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF] for their Active Kriging Monte Carlo simulation technique (AK-MCS). In AK-MCS, some candidates to enrichment are considered among a Monte Carlo set. The point that is most likely to improve the Kriging model is defined as the one that minimizes the following U -function:

U (x) = ḡ -µ M (x) σ M (x) . (21)
In practice, points that tend to minimize this function are those which are close to the constraint threshold i.e. µ M (x) → ḡ (otherwise put, g (x) → 0), or those for which the Kriging variance is high (σ M (x) → ∞), thus implying that the Kriging model may lack of accuracy there because of the DoE scarcity.

In the sequel, we adapt this function for contour estimation with respect to a quantile which is referred to as the global stage of enrichment.

Proposed global stage of enrichment

This first stage of enrichment is aimed at revealing regions of the space where the constraints, as defined in terms of quantiles, are likely to be violated. We call it global as this enrichment spans the entire augmented space just as the AK-MCS defined above. There is however one difference in our setting. In contrast to AK-MCS, the constraint is defined with respect to d in the design space but the Kriging model is built in the augmented space. The idea with the proposed approach is to find the pair of points in the augmented space that most likely leads to an improvement of the quantile estimation in the design space. The following steps are completed to achieve this task:

1. Sample candidates for enrichment in the design space:

C = d (1) , d (2) , . . . , d (m) (22)
2. For each design d (i) , i = {1, . . . , m}:

(a) Sample the Monte Carlo set required to compute the quantile:

C (i) q = {(x j , z j ) , j = 1, . . . , N } (23)
(b) Compute the associated quantile q α d (i)

(c) Identify the point in the augmented space that is associated to the quantile, i.e.

x (i) α , z (i) α = (x, z) ∈ C (i) q : q α d (i) = µ M (x, z) (24)
(d) Compute the modified deviation number:

U d (i) ≡ U x (i) α , z (i) α = ḡ -µ M x (i) α , z (i) α σ M x (i) α , z (i) α (25)
3. The next best point to add to the design of experiments is therefore defined as:

(x next , z next ) = arg min (xα,zα)∈Cα U (d) , (26)
where

C α = x (i) α , z (i) α
, i = 1, . . . , m .

To illustrate this enrichment scheme, let us consider the mathematical function from Janusevskis and Le Riche (2013) which reads:

M (d, z) = 1 3 z 4 -2.1z 2 + 4 z 2 + dz + 4d 2 d 2 -1 , (27)
where d ∈ [-1, 1]. It is considered as a performance function for an RBDO problem where the constraint threshold is set to ḡ = 0.5. The probabilistic model consists of the random design variable X ∼ N d, 0.05 2 and the random environmental variable Z ∼ N 0.5, 0.05 2 . Figure 3 shows the various iterations of the enrichment procedure. The left panel shows the augmented space with contour of ϕ(-U ), where ϕ is the standard Gaussian PDF (U is conveniently mapped for proper illustration). The contour g (x, z) = 0 is plotted as black dotted line and the small black crosses form the set C α . In the right panel, the true and estimated quantiles (q α and q α ) are respectively plotted in blue and black lines. The threshold ḡ = 0.5 is represented by the red dotted line. The blue triangles are the initial DoE. As enrichment is carried out, the red squares are added to the DoE. At each iteration, the best next point corresponds to the black diamond. From this example, we can see that the points added in the augmented space actually corresponds to those where q α (d) → ḡ. Figure 3: Enrichment with the mathematical function. In the left panel, the augmented space with contours of the enrichment functions and the set C α shown as small crosses. In the right panel, the true and estimated quantiles are shown (resp. blue and black curve). Triangles and squares respectively stand for initial and enrichment points.

In this example, the enrichment was stopped when min d∈C U ≥ 2. This criterion actually means that there is only 5% of chance of mistaking a safe design for a failed one (and vice-versa) w.r.t. all the points in C. This is quite a conservative stopping criterion. We do not need such an accuracy in the entire design space. Since the next step is optimization, we may go further in reducing the computational budget by saving model evaluations to regions that actually improve the objective function. For this reason we propose a second stage of enrichment as explained in the sequel.

Local stage of enrichment

In order to further reduce the number of calls to the original computational model, we stop the first and global stage of enrichment earlier and proceed to a local enrichment which is coupled with optimization. In fact, the idea is to have a roughly accurate surrogate model that reveals the different regions of the space where the approximated limit-state function is close to zero.

To this end, we may relax the criterion min U ≥ 2 to a certain proportion of the enrichment candidates rather than all of them. The criterion may therefore be written as:

η = Card (C 2 ) /Card (C) ≤ η. ( 28
)
where C 2 = {d ∈ C : U (d) ≤ 2}. Note that the original criterion corresponds to η = 0. We may consider a relaxed criterion by setting η = 0.30 for instance.

Assuming that the first stage of enrichment has been stopped with respect to the above criterion, there is residual epistemic uncertainty to the Kriging model. This uncertainty can be monitored during optimization and dealt with by updating the Kriging model only when necessary. To achieve this goal, we may consider a local accuracy measure associated to the quantile estimates, as they ultimately define the constraints of interest.

Following the idea in [START_REF] Dubourg | Reliability-based design optimization using Kriging and subset simulation[END_REF] where bounds on failure probabilities were developed, we define the following lower and upper bounds, q - α and q + α , which are quantiles computed with respect to µ M -2σ M and µ M + 2σ M . Since the standard deviation is positive the following relationship holds:

q - α (d) ≤ q α (d) ≤ q + α (d) for any d ∈ D. (29)
The spread of this interval is a good measure of the local Kriging accuracy for the quantile estimation. The following local accuracy criterion may henceforth be derived:

η q (d) = q + α (d) -q - α (d) ḡ ≤ ηq , ( 30
)
where ηq is a pre-defined threshold. In the case where ḡ = 0, we may replace the denominator by Var Y M CS , where Var Y M CS is the variance of Kriging prediction over a large Monte Carlo set sampled in the augmented space.

The surrogate model is considered to be accurate enough for the quantile estimation at the design d (i) if this relationship holds. If in contrast η q > ηq , then a local enrichment is made. To this end, candidates for enrichment are selected among the Monte Carlo set C

(i)

q . The following deviation number is computed on this set [START_REF] Schöbi | PC-Kriging: A new meta-modelling method and its applications to quantile estimation[END_REF][START_REF] Schöbi | Rare event estimation using Polynomial-Chaos-Kriging[END_REF]:

U (x, z) = q α d (i) -µ M (x, z) σ M (x, z) . ( 31
)
The best next point is the one that minimizes this function. This point corresponds to a certain

x (i) α , z (i) α
from Eq. ( 24). By iteratively adding points in this fashion, it is expected that the quantile will be more and more accurately estimated.

Implementation of the proposed procedure

We now consider the implementation of the whole procedure. Prior to that, let us specify two additional points that most often characterize the actual problems we intend to address, i.e. the possibility of adding many points per iterations and the presence of multiple constraints. The first point may be interesting when one has computational resources that allow for distributed computations. It may also be argued that there is not one single point that is likely to improve the surrogate model but many points located in disjoint regions. In such a case, multiple enrichment points allow us to reach them simultaneously. In order to add K points in the DoE, we consider a weighted K-means clustering of the candidates for enrichment, where each point is weighted by ϕ(-U ). This way, regions with small values of U are favored. Finally, K clusters centers are chosen as the next points to add in the DoE.

As for the case of multiple constraints, many techniques exist. We may, for instance, rank the constraints and enrich sequentially starting with the most important one. This is not an optimal scheme. Fauriat and Gayton (2014) proposed a composite criterion which focuses on the most violated constraints. However, the notion of "most violated" is not adequate when the constraints are defined on completely different scales. In this work, we thus consider a composite criterion where, for each enrichment candidate, the constraint with minimum value of U is taken, that is:

U comp (x, z) = min l∈{1,...,n h } U l (x, z) = ḡl -µ M l (x, z) σ M l (x, z) , U comp (x, z) = min l∈{1,...,n h } U l (x, z) = q α l d (i) -µ M l (x, z) σ M l (x, z) . ( 32
)
Considering all these developments, the pseudo-code in Algorithm 1 summarizes the proposed procedure. Here we consider that the first stage of enrichment has been already performed. The selected optimization algorithm is the (1+1)-CMA-ES (Covariance matrix adaptation -evolution strategy) for constrained problems [START_REF] Arnold | A (1+1)-CMA-ES for constrained optimisation[END_REF]. This is a stochastic global search algorithm which relies on multivariate normal distributions to search candidates with increased fitness as iterations grows. It also accounts for constraints by decreasing the likelihood to sample in the direction of previously unfeasible sampled points. Such a global search algorithm is quite convenient for the proposed procedure since only one parent generates one offspring, thus allowing us to check the quantile accuracy for the offspring before moving on. The entire procedure is illustrated in Figure 4.

Initialize the DoE and the Kriging models

Proceed to global enrichment (Section 4.3)

i = 0 Build C (i) q
Compute the quantile q α (d (i) ), the bounds q ± α (d (i) ) and the criterion η

(i) q Is η (i) q ≤ ηq ?
Update DOE following Eq. ( 32)

Update Kriging models Run one iteration of CMA-ES:

d (i+1) = d (i) + ν (i)
Convergence achieved ?

End Draw samples

i = i +
C (i) q = {(x 1 , z 1 ) , . . . (x N , z N )} in the augmented space where X ∼ f X|d (i) and Z ∼ f Z 4:
for l = 1 to n h do 5:

for j = 1 to N do 6:

y j = µ M l (x j , z j ) 7: y - j = µ M l (x j , z j ) -2σ M l (x j , z j ) 8: y + j = µ M l (x j , z j ) + 2σ M l (x j , z j ) 9:
end for 10:

q α l d (i) = quantile { y} N j=1 , α l Estimated quantile 11: q - α l d (i) = quantile { y -} N j=1 , α l
Lower bound of the quantile 12:

q + α l d (i) = quantile { y + } N j=1 , α l
Upper bound of the quantile 13:

end for 14:

if (q + αq - α ) /ḡ > ηq then 15:

for k = 1 to N do 16:

for l = 1 to n h do 17:

U l (x k , z k ) = µ M l (x k , z k ) -q α l /σ M l (x k , z k ) 18:
end for 19:

U comp (x k , z k ) = min l∈{1,...,n h } U l (x k , z k ) 20:
end for 21:

if K == 1 then

The point that minimizes Ucomp is chosen 22: end if

(x next , z next ) = arg min {U comp (x k , z k )} N

27:

Update the DoE D and the Kriging models M l 28:

d (i) ← d (i) + ν (i)
Explore the next design point using (1 + 1)-CMA-ES

29:

i ← i + 1 30:

Check convergence of the optimization algorithm 31: end while

Application examples

The proposed methodology is now validated with four application examples. The first three are analytical problems whose solutions are available in the literature. The last one is related to the lightweight design of an automotive body structure under crashworthiness constraints.

The following settings are common to all the problems. L 2 -discrepancy-based optimized Latin hypercube is used to generate the initial designs of experiments. Anisotropic Kriging with Matérn 5/2 autocorrelation function and a constant trend is considered as the default surrogate model.

Column under compression

This first example, introduced in Dubourg ( 2011), is concerned with a column of rectangular cross-section b × h submitted to a compressive load F ser . The aim is to minimize the crosssectional area while avoiding buckling. Buckling may occur here if the service load is higher than critical Euler force which reads:

F cr = π 2 EI L 2 , ( 33
)
where L is the length of the column, E is the Young's modulus of its constitutive material and

I = bh 3 /12 (b > h) is the column area moment of inertia.
The deterministic optimization problem then reads: 150,350] 2 bh subject to:

d * = arg min d∈[
   f (d) = h -b ≤ 0, g (d, z) = F ser -k π 2 Ebh 3 12L 2 , ( 34
)
where k is a parameter which accounts for model uncertainty in the Euler force (i.e. it represents the effect of imperfections in the beam geometry and may be viewed as a model correction factor with respect to the ideal Euler force) and z = {k, E, L} T is the vector of environmental variables.

Uncertainties are considered by introducing the probabilistic model as described in Table 1. With all parameters being lognormally distributed, an analytical solution can be derived [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]:

b * = h * = 12F ser π 2 exp λ k + λ E -2λ L + Φ -1 (P f ) ζ 2 k + ζ 2 E + 4ζ 2 L , (35)
where

ζ • = ln (1 + δ 2 • ) and λ • = ln (µ • )-1 2 ζ 2
• are respectively the scale and location parameters of the lognormal distribution. By setting the target probability of failure to 5%, i.e. α = 0.95, the analytical solution, b * = h * = 238.45 mm.

To apply the methodology on this five-dimensional problem, we start with a scarce initial design of 10 points and set the global accuracy threshold in Eq. ( 28) to η = 0.15. Only two enrich- We then start the optimization by setting a simulated-annealing-like threshold ηq with three levels which are respectively 1, 0.5 and 0.1. The idea is to start with a relaxed threshold in the early iterations where CMA-ES is exploring and gradually reduce it as iterations grow and CMA-ES starts exploring identified local minima. With this optimally tuned scheme, six points are added to the DoE. The found solution is b * = h * = 239.12 mm, has 0.28% discrepancy with the analytical solution. Note that the exact solution can be reached, should we increase the number of iterations of CMA-ES or refine the solution by a gradient-based algorithm. Figure 5 illustrates the convergence of CMA-ES. In total, only 18 points were necessary to achieve convergence. By comparison, a oneshot approach with a DoE of size 18 does not systematically converge to the reference solution.

Additional points are needed most of the time. The blue and red dots fall respectively in the feasible and unfeasible sets. The green ones are the successive best sample points. In the right panel, the evolution of the local accuracy criterion is shown with respect to the number of iterations. In the two figures, the points corresponding to enrichment have been circled in cyan.

Two-dimensional problem

This analytical example has been widely used for benchmark purposes in the related literature (Du and Chen, 2004;[START_REF] Shan | Reliable design space and complete single-loop reliability-based design optimization[END_REF][START_REF] Liang | A single-loop method for reliability-based design optimization[END_REF][START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]. The optimization problem consists in minimizing the sum of the design parameters under three non-linear limit state functions whose deterministic formulation reads:

d * = arg min d∈[0,10] 2 d 1 + d 2 s.t.:                    g 1 (d) = d 2 1 d 2 20 -1 ≤ 0 g 2 (d) = (d 1 + d 2 -5) 2 30 + (d 1 -d 2 -12) 2 120 -1 ≤ 0 g 3 (d) = 80 (d 2 1 + 8d 2 + 5) -1 ≤ 0 . (36)
In order to solve the RBDO problem, we consider the following setting. The two design variables are considered as random:

X i ∼ N d i , 0.6 2 , i = {1, 2}.
The target failure probability is Pfi = 1.35 • 10 -3 and thus corresponds to

β i = 3 for i = {1, 2, 3}.
We start the procedure with a 10-point experimental design. Considering η = 0.3, five points are added during the first stage of enrichment. Figure 6a shows the convergence of this enrichment stage. In Figure 6b, the contours of the limit-state with respect to the current Kriging models in the augmented space X = [-1.8, 11.8] 2 are shown. In this figure, the black rectangle corresponds to the bounds of the design space D = [0, 10] 2 and the initial and added points are shown respectively as blue triangles and red squares. 2 compares the results presented here with those reported in the literature for benchmark. All the selected methods provide a solution with a good accuracy. The difference rather lies in their costs. In these results, the single loop and the reliable design space (RDS) approaches require a relatively small number of functions evaluation despite they do not rely upon surrogate models. The two cases considering surrogate models (Meta-RBDO and Quantile-RBDO) are among the best in terms of model evaluations. In this example, the quantile-based approach we propose is on average the less expensive one. Since the initial design is random, the analysis is replicated 50 times. The number of calls to the true models varies between 11 and 23, all of them leading to good results. On average the number of calls is 14.6 and among them only three are above 20.

Bracket structure

This mechanical example consists of the two-member bracket structure illustrated in Figure 8 ( [START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF]. The two-members are pin-joined at the point B and a vertical load P is applied on the right end of the member CD at a distance L of its hinge. The aim is to minimize its weight while considering two failure modes: • The bending stress in the member CD whose maximum value σ b is required to be smaller than the yield stress σ y :

g 1 (d, z) = σ y -σ b , (37)
where σ b = 6M B /w CD t 2 , with M B = P L/3 + ρgw CD tL 2 /18. Here w CD and t are the cross-sectional dimensions of CD, ρ is the unit mass of its constitutive material and g is the gravity intensity.

• The compression force F AB that must be lower than the critical Euler force F b :

g 2 (d, z) = F b -F AB with: F b = π 2 EI L 2 AB = π 2 Etw 3 AB 12 (2L/3sinθ) 2 , F AB = 1 cosθ 3P 2 + 3ρgw CD tL 4 , (38)
where w AB and L AB are respectively the width and length of AB and θ is its inclination angle.

The RBDO application problem as in [START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF] and [START_REF] Dubourg | Reliability-based design optimization using Kriging and subset simulation[END_REF] consists in minimizing the weight of the structure, given the following objective function:

c (d) = ρtL 4 √ 3 9 w AB + w CD , (39)
where d = {w AB , w CD , t} T ∈ D = [5, 30] 3 represents the set of design parameters.

The target reliability index for this problem is set to β 1 = β 2 = 2, thus corresponding to a component failure probability of 0.0227. The RBDO problem therefore reads:

d * = arg min d∈D c d subject to:    P g 1 X (d) , Z ≤ 0 ≤ Pf1 P g 2 X (d) , Z ≤ 0 ≤ Pf2 , (40)
where Pfi = Φ (-β 1 ) ≈ 0.0227 and the functions g 1 , g 2 and c are respectively given by Eqs. ( 37)

-( 39).

The probabilistic model associated to this problem is shown in Table 3. Table 4 shows the bounds of the augmented space in which the training points are sampled. The surrogate model is built in the unit hypercube following a linear mapping from this augmented space. For the solution of this problem, we start with an initial design of 50 points. The threshold for the first stage of enrichment is set to η = 0.30. A total of 60 enrichments points have been added to reach the required accuracy through 6 iterations of 10 points each. The optimization is then initiated starting from d (0) = {6.1, 20.2, 26.9} T which is also the initial design in the benchmark references [START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF]; [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF] and corresponds to the optimal deterministic solution. In the second stage, K = 3 points are added per enrichment, thus leading to 15 additional points in the DoE. Convergence of the CMA-ES algorithm is illustrated in Figure 9 where the evolution of the cost function with respect to the iteration number is

shown. The green circles highlight the points that were feasible and improved the current best design. The CMA-ES algorithm is stopped after 150 iterations and the solution locally refined through a gradient-based approach using the final Kriging model. The overall number of calls to the original model is 125 for this illustrated case.

As in the previous example, we replicate the optimization 20 times because of the random nature of the initial experimental design. The number of calls varies between 80 and 170, with the maximum clearly being an outlier. On average, this number of calls is 107. In the light of this result, the proposed procedure is more efficient than the approaches from the two benchmark references as shown in Table 5. In this table, the brute force approach refers to a solution that is found by a quantile-based procedure directly relying on the true mechanical models rather than surrogates. Beside, the resulting weight saving is higher in our approach. This may be explained by the fact that we use a global optimization algorithm rather than a gradient-based one as was done in the two references.

Sidemember subsystem

This final application is related to the lightweight design of an automotive body structure under crashworthiness constraints. This involves finding the best distribution of the metal sheet thicknesses which allows one to satisfy frontal impact-related constraints. These constraints are evaluated by finite element crash simulations which are extremely time-consuming, i.e. 24 hours for a single model run on distributed CPUs. The use of surrogate models is therefore the only alternative in order to perform such an optimization. In this application, we consider the so-called sidemember subsystem which is a subset of the front end of a vehicle. This subsystem actually has the same behavior in frontal impact as a full vehicle, yet requires reduced computational time (10 to 15 minutes on a cluster of 48 CPUs). The sidemember subsystem is illustrated in Figure 10. Five parts are considered for optimization as shown in the figure . To account for noise which is inherent to frontal impact, some parameters of the crash protocol are considered as random. These are respectively the initial speed and the position of the barrier:

V ∼ U (34, 35) km/h and P ∼ N (0, 2) mm. The two constraints that are considered for this problem are the maximum wall contact force that should not be larger than ḡ1 = 170 kN and the maximum sidemember compression which should be kept below ḡ2 = 525 mm. In order to obtain a conservative design with respect to uncertainties, the quantile-based optimization procedure is applied with a quantile level α 1 = α 2 = 0.95. The associated RBDO problem reads as follows:

d * = arg min d∈D c d subject to:    P (ḡ 1 -M 1 (d, Z) ≤ 0) ≤ 0.05 P (ḡ 2 -M 2 (d, Z) ≤ 0) ≤ 0.05, (41)
where Z = {V, P } T and M 1 and M 2 are the outputs of the finite element model giving respectively the maximum wall force and maximum sidemember compression.

The initial design and the bounds of the augmented space associated to this problem are given in Table 6. Since the design variables are deterministic X reduces to the design space. to a weight saving of 1.08 kg, that is 11.5% of the initial weight, which is considered significant in car manufacturing. The thicknesses associated to the initial and optimal solutions are shown in Figure 11. The corresponding quantile constraints for this solution are q α1 (d * ) = 155.62 kN and q α2 (d * ) = 523.12 mm, which are below the thresholds. The validation of the reliability of this solution with respect to the finite element model is not possible due to the large cost of a single run. Instead, we focus here on the accuracy of the Kriging surrogates in the vicinity of the optimal design. To this end, we estimate quantiles with the original and surrogate models from a set of Monte Carlo samples of size 100. To account for this reduced size, we consider the mean value of the quantiles estimates obtained from 500 bootstrap replicates. Each bootstrap replicate consists in sampling with replacement 100 points from the available data. The resulting relative error is still biased because of the small size of the Monte Carlo set. However, this allows us to give a flavor on the ability of the surrogate model to approximate the quantile in the vicinity of the found solution. Table 7 compares the results by considering the finite element model ( q FE α ) on the one hand and the Kriging model ( q KRG α ) on the other. The two responses are quite close for each output, showing that the Kriging models were accurate enough (at least locally) for the purpose of quantile estimation.

Conclusion

The aim of this paper is to propose a quantile-based, conservative optimization procedure for structures in an uncertain environment. Furthermore, structures whose behavior is simulated by high-fidelity and expensive-to-evaluate models are considered. Such simulations are timeconsuming. Surrogate modeling approaches are therefore introduced as computationally costless approximations of these models.

The optimization problem is first posed in the framework of reliability-based design optimization (RBDO). After a brief review of the most-widely used techniques to solve a RBDO problem, we formulate a new quantile-based approach of optimal design. This approach is motivated by the relatively high target failure probabilities that can be accepted in the applications under consideration in the field of car body design. These probabilities of failure will be estimated by crude Monte Carlo sampling.

Kriging, with its basic equations, is introduced for the purpose of surrogate modeling. To further reduce the computational burden associated to building the Kriging surrogate model, a two-stage enrichment of the design of computer experiments is proposed in an augmented space that combines both design variables and uncertain environmental variables. The first stage, which is global, aims at reducing the overall Kriging epistemic uncertainty by adding points in the vicinity of the limit-state surfaces. The second stage, which is local, is embedded in the optimization procedure. At each iteration, the accuracy of the estimated quantiles is checked.

Enrichment of the design of experiments is made locally only when the accuracy is not sufficient.

This allows us to direct the experimental design points to regions of the space that decrease significantly the cost function while ensuring that the performance criteria are fulfilled.

Three applications are considered to validate the proposed procedure. The first one is a five-dimensional example related to a beam buckling problem, whose analytical solution can be computed. This allows us to validate the proposed method against exact results. The second and third problems respectively involve three non-linear limit state functions and a bracket structure.

The optimal solutions obtained from different approaches are already available in the literature.

The application of the proposed procedure shows increased efficiency as the number of calls to the original computational model is reduced. For the bracket structure, a better solution in terms of the cost function is even found compared to the best results available in the literature.

Finally, we apply the methodology to an industrial problem related to the lightweight design of an automotive sidemember subsystem under frontal impact. A reliable solution is found within a reasonable number of calls to the expensive finite element model. All these applications feature relatively low-dimensional problems. Applications to high-dimensional cases, say s > 20, is still a challenging task and require further work.

  , a reliability-based design optimization may be formulated as follows: d * = arg min d∈D c (d) subject to:    f j (d) ≤ 0, {j = 1, . . . , n s } , P (g k (X (d) , Z) ≤ 0) ≤ Pf k , {k = 1, . . . , n h } ,

  (a.k.a. Gaussian process modeling) has been successfully used for reliability analysis in Echard et al. (2011); Picheny et al. (2010); Bichon et al. (2008); Balesdent et al. (2013). For the specific task of RBDO, conservative surrogate models which rely on Kriging or polynomial response surfaces were considered in Viana et al. (2010); Picheny et al. (2008). From another perspective, Dubourg et al. (2011); Chen et al. (2015); Lee et al. (2011); Li et al. (2016), for instance, have

Figure 1 :

 1 Figure 1: Comparison of a safe and an unsafe design with respect to a quantile Q α corresponding to a target failure probability Pf = 1α.

  Augmented space (left) and design space (right): iteration #10

Figure 4 :

 4 Figure 4: Flowchart of the optimization procedure with the two stages of enrichment.

  , z next ) obtained from weighted K-means clustering with weight ϕ (-U comp )

Figure 5 :

 5 Figure 5: Convergence of the column under compression problem. The left panel shows the evolution of CMA-ES in the design space. The blue and red dots fall respectively in the feasible and unfeasible sets. The green ones are the successive best sample points. In the right panel, the evolution of the local accuracy criterion is shown with respect to the number of iterations. In the two figures, the points corresponding to enrichment have been circled in cyan.

  Limit-state surface after the first stage of enrichment

Figure 6 :Figure 7 :

 67 Figure 6: Illustration of the first stage of enrichment for the two-dimensional problem

Figure 8 :

 8 Figure 8: A sketch of the bracket structure (as illustrated in Dubourg (2011))

Figure 9 :

 9 Figure 9: Convergence of the CMA-ES algorithm for the bracket structure.

  Figure 10: Sidemember subsystem

Figure 11 :

 11 Figure11: Comparison of the initial and optimal designs with respect to each parameter for the 5-part sidemember subsystem. The blue color stands for the nominal design and the green color for the optimal one.

Table 1 :

 1 Probabilistic model for the column under compression example.

	Parameter Distribution	Mean (µ)	COV (δ%)
	k	Lognormal	0.6	10
	E (MPA)	Lognormal	10, 000	5
	L (mm)	Lognormal	3, 000	1
	F ser (N)	-	1.4622 × 10 6	-
	ment points are necessary to reach the required global accuracy.

Table 2 :

 2 Results comparison for the Choi problem.

	Method Brute force PMA 1	d * 1 3.45 3.30 6.75 d * c (d * ) 2 3.43 3.29 6.72	g-calls ≈ 10 6 1, 551
	SORA 2	3.44 3.29 6.73	151
	Single loop 3	3.43 3.29 6.72	19
	RDS 1	3.44 3.28 6.72	27
	Meta-RBDO 4	3.46 3.27 6.74 20(20/10/10)
	Quantile-RBDO 3.44 3.29 6.73	17
	1 As calculated in Shan and Wang (2008)
	2 As calculated in Du and Chen (2004)
	3 As calculated in Liang et al. (2004)
	4 As calculated in Dubourg (2011)

Table 3 :

 3 Parameters of the variables defining the probabilistic model for the bracket structure problem: d = {w AB , w CD , t} T are the design variables and z = {P, E, σ y , ρ, L} T are the environmental variables.

	Parameter	Distribution Mean COV (δ%)
	Width of AB (w AB in m)	Normal	w AB	0.05
	Width of CD (w CD in m)	Normal	w CD	0.05
	Thickness (t in m)	Normal	t	0.05
	Applied load (P in kN)	Gumbel	100	0.15
	Young's modulus (E in GPa)	Gumbel	200	0.08
	Yield stress (σ y in MPa) Unit mass (ρ in kg/m 3 )	Lognormal Weibull	225 7860	0.08 0.10
	Length (L in m)	Normal	5	0.05

Table 4 :

 4 Bounds of the augmented space for the bracket structure problem.

	Parameter	Lower bound Upper bound
	Width of AB (w AB in m)	4.25	34.5
	Width of CD (w CD in m)	4.25	34.5
	Thickness (t in m)	4.25	34.5
	Applied load (P in kN)	15.98	109.58
	Young's modulus (E in GPa)	110.38	224.31
	Yield stress (σ y in MPa) Unit mass (ρ in kg/m 3 )	176.49 4760.20	285.01 9576.3
	Length (L in m)	4.25	5.75

Table 5 :

 5 Comparative results for the bracket structure. The PMA result comes from[START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF] and Meta-RBDO from[START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF].

	Design method	Weight (kg) w AB (cm) w CD (cm) t (cm)	g-calls
	Brute force PMA 1	1357 1673	5.35 6.08	7.40 15.68	30.00 20.91	≈ 10 6 2340
	Meta-RBDO 2	1584	5.80	12.80	23.30 160(160/90)
	Quantile-based RBDO	1364	5.57	7.28	30.00	125

1 As computed by

[START_REF] Chateauneuf | Structural design optimization considering uncertainties[END_REF]

2 As calculated by

[START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]

Table 6 :

 6 Bounds of the augmented space and initial design for the sidemember subsystem.

	Lower Upper	1.5 2.5	1.5 2.5	2 3	1.5 2.5	0.6 1.2	34 35	-6 6
	Initial	2	2	2.5	2	0.9	-	-
			Lower bulkhead
		Rear side-member
		Forward side-member
		Wheel arch
	Forward side-member base

Param. d 1 (mm) d 2 (mm) d 3 (mm) d 4 (mm) d 5 (mm) V (km/h) P (mm)

Table 7 :

 7 Quantiles of the performance criteria g 1 and g 2 computed from the Kriging model ( q KRG α ) (resp. the original model ( q FE α )) obtained from 100 Monte Carlo samples, averaged over 500 bootstrap replicates.

	Criterion	g 1 (kN) g 2 (mm)
	Original model q FE α Kriging model q KRG α	150.66 148.02	527.81 528.04
	Error (%)	1.75	0.04