R Schöbi

B Sudret

S Marelli

Rare event estimation using Polynomial-Chaos-Kriging

Keywords: Adaptive experimental design, design enrichment, Kriging, meta-modeling, PC-Kriging, rare events

Structural reliability analysis aims at computing the probability of failure of systems whose performance may be assessed by using complex computational models (e.g. expensive-to-run finite element models). A direct use of Monte Carlo simulation is not feasible in practice, unless a surrogate model (such as Kriging, a.k.a Gaussian process modeling) is used. Such meta-models are often used in conjunction with adaptive experimental designs (i.e. design enrichment strategies), which allows one to iteratively increase the accuracy of the surrogate for the estimation of the failure probability while keeping low the overall number of runs of the costly original model.

In this paper we develop a new structural reliability method based on the recently developed Polynomial-Chaos Kriging (PC-Kriging) approach coupled with an active learning algorithm known as AK-MCS. We formulate the problem in such a way that the computation of both small probabilities of failure and extreme quantiles is unified. We discuss different convergence criteria for both types of analyses, and show in particular that the original AK-MCS stopping criterion may be over-conservative. We finally elaborate a multi-point enrichment algorithm which allows us to add several points in each iteration, thus fully exploiting high-performance computing architectures.

The proposed method is illustrated on three examples, namely a two-dimensional case which allows us to underline the advantages of our approach compared to standard AK-MCS. Then the quantiles of the 8-dimensional borehole function are estimated. Finally the reliability of a truss structure (10 random variables) is addressed. In all case, accurate results are obtained with about 100 runs of the original model.

Introduction

Advanced computer simulations of complex physical systems (e.g. finite element models) are common in modern engineering to assess and optimize their performance. At the same time, awareness is growing on the concepts of structural reliability and robust design, hence making the quantification and propagation of uncertainties a key issue [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical modelscontributions to structural reliability and stochastic spectral methods[END_REF]De Rocquigny et al., 2008).

Due to the high cost of repeatedly evaluating complex computational models, analyses with classical sampling techniques such as Monte Carlo (MC) are often intractable. In this context, meta-modeling techniques allow one to develop fast-to-evaluate surrogate models from a limited collection of evaluations of the original computational model, which is called the experimental design. Commonly used techniques include Kriging [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], Polynomial Chaos Expansions [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF], support vector machines [START_REF] Gunn | Support vector machines for classification and regression[END_REF] and neural networks [START_REF] Hurtado | Neural-network-based reliability analysis: a comparative study[END_REF].

Two cases are generally distinguished in the context of rare event estimation, namely the estimation of failure probabilities (so-called structural reliability analyses) and that of extreme quantiles. Both quantities describe important statistics of the response of the computational model. In the context of engineering and high reliability targets, failure probabilities are generally low (in the order of 10 -3 to 10 -6 ) and the quantiles extreme (e.g. 99% quantile), therefore structural reliability falls into the broad class of rare events estimation.

Rare events estimation has been in turn often related to the concept of experimental design enrichment (a.k.a. adaptive experimental design). The core idea is to start calibrating a meta-model with a limited set of computational runs and then to add new model evaluations iteratively to increase the prediction accuracy of the rare event statistics of interest. Kriging (a.k.a. Gaussian process modeling) has been introduced by [START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF] and Schueremans and Van Gemert (2005b) to the field of structural reliability. The use of adaptive experimental design algorithms has been investigated by [START_REF] Bichon | Efficient global reliability analysis for nonlinear implicit performance functions[END_REF], [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF]Ran-jan et al. (2012) who combine an adaptive Kriging meta-model with Monte Carlo simulation, and [START_REF] Echard | A combined importance sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models[END_REF] and [START_REF] Cadini | An improved adaptive Kriging-based importance technique for sampling multiple failure regions of low probability[END_REF], who modify it for the case of importance sampling. Under the name "stepwise uncertainty reduction" [START_REF] Arnaud | Evaluation d'un risque d'inondation fluviale par planification séquentielle d'expériences[END_REF], [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF] and [START_REF] Jala | Plans d'expériences adaptatifs pour le calcul de quantiles et application à la dosimétrie numérique[END_REF] discuss an adaptive Kriging approach in a Bayesian context. [START_REF] Bourinet | Assessing small failure probabilities by combined subset simulation and support vector machines[END_REF] use instead support vector machines in the context of subset simulations. [START_REF] Dani | Stochastic linear optimization under bandit feedback[END_REF] and [START_REF] Srinivas | Information-theoretic regret bounds for Gaussian process optimization in the bandit setting[END_REF] use upper confidence bounds. [START_REF] Dai | Wavelet density-based adaptive importance sampling method[END_REF] discuss the use of wavelet-density-based adaptive importance sampling. Finally, the use of adaptive Kriging to build up a quasi-optimal importance sampling density (thus leading to unbiased estimates of the failure probability) has been recently proposed in [START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF] and [START_REF] Dubourg | Metamodel-based importance sampling for reliability sensitivity analysis[END_REF].

In this paper we combine Polynomial-Chaos-Kriging (PC-Kriging) [START_REF] Schöbi | PC-Kriging: a new metamodelling method combining polynomial chaos expansions and Kriging[END_REF][START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF][START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions -Application to uncertainty analysis in computational dosimetry[END_REF] with design enrichment to solve reliability analysis and extreme quantile problems. We focus on the best choice of additional samples in the experimental design as well as on the option of adding multiple samples in each iteration of the algorithm. Adding multiple samples in each iteration would allow one to take advantage of distributed computing facilities. This paper is organized as follows: after the introduction of a proper formalism for structural reliability analysis, the principles of Polynomial-Chaos-Kriging [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF] are given.

An existing adaptive experimental design algorithm is then presented and extended to PC-Kriging and the concept of distributed computing. It is shown how the same algorithm can be used to estimate both failure probabilities and extreme quantiles with only slight modification. The application of the resulting framework is finally demonstrated on a lowdimensional analytical function and on two realistic engineering applications.

Structural reliability analysis

Computational model

In this paper, a computational model M is defined as a mapping of the M -dimensional input vector of parameters x to the one-dimensional output scalar y, i.e. M :

x ∈ D X ⊂ R M → y = M(x) ∈ R.
The uncertainties in the input parameters are represented by a random vector X with joint probability density function (PDF) f X . The components of X = {X 1 , . . . , X M } are assumed independent throughout this paper for the sake of simplicity. Moreover, when the components are dependent, it is always possible to recast the problem in terms of independent variables using the Nataf or Rosenblatt transform [START_REF] Lemaire | Structural reliability[END_REF]Lebrun and Dutfoy, 2009a,b). The model response then becomes a random variable Y obtained by propagating the input uncertainty in X through the computational model:

Y = M(X).
(1)

It is further assumed that the computational model M is a black-box, i.e. for each realization

x of the input random vector, only the corresponding response y = M(x) is accessible. In other words we assume that the inner structure and mechanisms of the computational model are not observable. This is typically the case for e.g. a finite element model in which, generally speaking, the governing equations cannot be solved analytically.

Estimation of a failure probability

In structural reliability analysis, the failure probability P f is defined as the probability that a model response Y is smaller than a given threshold value y 0 :

P f = P (Y ≤ y 0 ) = P (M(X) ≤ y 0 ) . (2)
The failure probability may be recast as the following integral:

P f = D f f X (x)dx s.t. D f = {x ∈ D X : M(x) ≤ y 0 }, (3)
where f X is the joint probability density function of the input vector X and D f is the so-called failure domain. Note that failure is defined as in Eq. ( 2) and (3) to ensure consistency with the following sections. In other literature, failure is often defined as a probability of exceedence, e.g. P (M(X) ≥ y adm ) which is also cast as P f = P (g(X) ≤ 0) where g(X) = y adm -M(X).

The integral in Eq. ( 3) has an implicit integration domain D f and cannot be solved analytically in the general case. A Monte Carlo estimation of the failure probability can be formulated to circumvent this limitation. Assume a reasonably large sample of X denoted by S = {x i , ..., x n }. The failure probability can be estimated by:

P f = n f n = 1 n n i=1 I M(x i )≤y 0 , (4)
where n f is the number of failure samples x i ∈ D f , n = |S| is the total number of samples and I M(x)≤y 0 is an indicator function for failure such that I = 1 for M(x) ≤ y 0 and I = 0 otherwise.

The accuracy of this estimate can be measured directly with the theoretical coefficient of variation for Monte Carlo samplings:

CoV P f = 1 -P f n • P f . ( 5
)
The main limitation of pure Monte Carlo approaches, namely the loss of efficiency for small failure probabilities, is apparent in Eq. ( 5). The required number of samples to obtain a target coefficient of variation in Eq. ( 5) rapidly increases with decreasing failure probability. For instance, to achieve a target coefficient of variation CoV P f = 0.1 given a failure probability

P f = 10 -k , n ≈ P f • CoV P f 2 -1
= 10 k+2 samples are needed.

Note that the indicator function I in Eq. ( 4) transforms the integration problem in Eq. (3) into a classification problem: identifying the failure (M(x) ≤ y 0 ) and safe (M(x) > y 0 ) domains.

Generally, it is sufficient to approximate the boundary between those domains, the limit state surface (M(x) = y 0 ), to accurately estimate the associated failure probability.

Estimation of a quantile

Quantile estimation consists in determining a quantile q α so that the probability of M(X)

being smaller than q α equals a preset value α, i.e.:

P (M(X) ≤ q α ) = α, (6)
where α = (0, 1). This constitutes an inverse problem compared to classical failure probability estimation shown in Eq. (2).

The quantile estimation problem can also be solved by Monte Carlo simulation. Consider again a reasonably large sample set S = {x 1 , ..., x n } and the corresponding response values y = {y 1 , ..., y n } = {M(x 1 ), . . . , M(x n )}. Assume the response values are ranked in ascending order and denote them by y (1) ≤ y (2) ≤ . . . ≤ y (n) . Then the estimator of the quantile q α reads:

q α = y nα , ( 7
)
where nα is the largest integer smaller than nα.

Analogously to failure probability estimation, the quantile estimation problem can also be cast as a classification problem:

find q α s.t. α ≈ 1 n n i=1 I M(x i )≤ qα . (8)

Limit state surface

Both the estimation of failure probabilities (Eq. ( 4)) and of quantiles (Eq. ( 8)) may be treated as classification problems dividing the input domain into failure and safe regions. We therefore introduce a limit state parameter a which determines whether a response y belongs to the failure (y ≤ a) or to the safe domain (y > a). For the estimation of failure probabilities a = y 0 , whereas for the estimation of quantiles a = q α . Due to those similarities we will denote failure probability and quantiles generically as statistics of interest in the remaining of this paper.

3 Meta-modeling

Kriging

When the computational model M is an expensive-to-evaluate function, Eq. ( 4) and Eq. ( 8) may become intractable. The behavior of the computational model M can then be approximated by a meta-model, such as Kriging (a.k.a. Gaussian process modeling). Kriging interprets the computational model as a realization of an underlying Gaussian process [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]:

M(x) ≈ M (K) (x) = β T • f (x) + σ 2 Z(x, ω), (9)
where β T • f (x) is the mean value of the Gaussian process (a.k.a. trend ), σ 2 is the process variance, Z(x, ω) is a zero-mean, unit-variance stationary Gaussian process and ω ∈ Ω denotes an elementary event in the probability space (Ω, F, P). The Gaussian process is characterized by an autocorrelation function R = R(|xx |; θ) and its hyper-parameters θ. In this paper, we consider two types of Kriging models: ordinary Kriging and universal Kriging. Ordinary

Kriging assumes a constant trend (f (x) = 1) whereas universal Kriging assumes a trend dependent on the input variables, e.g. f (x) = {1, x 1 , . . . , x M } for a linear trend.

Consider a set of N samples of the input vector X = {χ (1) , . . . , χ (N ) }, called experimental design, and the corresponding response of the exact computational model

Y = Y (1) = M(χ (1) ) . . . , Y (N ) = M(χ (N )
) . The estimation of the Kriging parameters {β, σ 2 } is computed by generalized least-squares solution [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]:

β(θ) = F T R -1 F -1 F R -1 Y, (10)
σ 2 (θ) = 1 N (Y -F β) T R -1 (Y -F β) , (11)
where

R ij = R(|χ (i) -χ (j) |; θ)
is the correlation matrix of the experimental design points and

F ij = f j (χ (i)
) is the information matrix. In case the hyper-parameters θ are unknown, their optimal values can be estimated through maximum likelihood estimation.

The prediction of the response value y of an arbitrary input sample point x is a Gaussian random variable characterized by [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]:

µ Y (x) = f (x) T β + r(x) T R -1 (Y -Fβ) , (12)
σ 2 Y (x) = σ 2 y   1 -f (x) T r(x) T   0 F T F R   -1   f (x) r(x)      , (13)
where µ Y (x) and σ 2 Y (x) are the prediction mean value and variance, and

r i (x) = R(|x - χ (i) |; θ)
is the correlation between the new sample x and the experimental design point

χ (i) ∈ X .

PC-Kriging

Framework

Polynomial-Chaos-Kriging (PC-Kriging) consists of a universal Kriging model whose trend is modeled by a set of orthogonal polynomials with respect to the input variables [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF]:

M(x) ≈ M (PCK) (x) = τ ∈T t τ ψ τ (x) + σ 2 Z(x, ω), (14)
where Due to the assumed independence of the input variables, the multivariate polynomials can be composed by tensor product of univariate orthogonal polynomials:

τ ∈T t τ ψ τ (x)
ψ τ (x) = M i=1 ψ (i) τ i (x i ), (15)
where

ψ (i) τ i (x i ) is a univariate polynomial of degree τ i in the i-th variable.
For each X i the orthonormality in coherency with the input distributions is given by:

ψ (i) j , ψ (i) k = D i ψ (i) j (x) ψ (i) k (x) f X i (x)dx = δ jk , (16)
where

ψ (i) j , ψ (i)
k are two univariate polynomials in the i-th variable, D i is the support of the distribution of X i , f X i is the marginal probability density function, and δ jk = 1 for j = k and δ jk = 0 otherwise. A summary of common univariate orthonormal polynomials can be found in [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF].

Optimal PC-Kriging

There are several ways to determine the most efficient set of polynomials T and combine it with universal Kriging. Amongst them is the so-called Optimal-PC-Kriging (OPC-Kriging) approach [START_REF] Schöbi | PC-Kriging: a new metamodelling method combining polynomial chaos expansions and Kriging[END_REF][START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF], which is presented briefly in this paper.

The information required for the algorithm comprises the input marginals, the autocorrelation function R, an experimental design X and the corresponding response values Y. The best sparse set of multivariate polynomials (summarized in the index set T ) is obtained by leastangle regression (LARS) [START_REF] Efron | Least angle regression[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]. LARS results in a ranked set of P polynomials. From this ordered set, P Kriging models are calibrated with an increasing number of polynomials each. The first trend is the polynomial chosen first in the LARS algorithm, while the trend of the P -th meta-model includes all P polynomials. The various meta-models obtained (with a polynomial trend having 1, 2, . . . , P terms respectively) are ranked according to their respective leave-one-out (LOO) cross-validation error and the meta-model with the minimal LOO error is chosen as the OPC-Kriging meta-model.

Characteristics

PC-Kriging combines a global with a local approximation of the behavior of the computational model. The global behavior is modeled by the PCE trend whereas the local behavior is modeled by the Gaussian process. This allows one to meta-model a computational model more accurately than PCE or other Kriging meta-models taken separately, see details in [START_REF] Schöbi | PC-Kriging: a new metamodelling method combining polynomial chaos expansions and Kriging[END_REF] and [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF].

Model accuracy

Universal Kriging-based meta-models generally approximate the behavior of the original model M most accurately close to the points in the experimental design. These points, however, are not always optimal for the estimation of the statistics of interest of Y . As an example, in order to compute a small failure probability P f = P (Y ≤ a) w.r.t. a given threshold a, it is of interest to concentrate the experimental design in the region close to the limit state surface (y ≈ a), which is usually unknown a-priori.

Consider the case of a given initial experimental design consisting of N 0 samples. By adding samples in a guided way, the accuracy of the estimation of the statistics of interest can be enhanced more efficiently than adding samples randomly. This is the main idea behind adaptive experimental design algorithms that are now presented.

Adaptive experimental design algorithm 4.1 Main algorithm

We adopt a slightly modified version of the Adaptive-Kriging-Monte-Carlo-Simulation (AK-MCS) algorithm introduced in [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF] for enriching the experimental design. An overview of the algorithm is shown in the flowchart in Figure 1. The white boxes mark the information required to start the algorithm, whereas the blue boxes represent the computational tasks. The main steps can be summarized as follows:

1. An initial experimental design X is generated by Latin-hypercube sampling and the corresponding exact response Y is computed.

2. A meta-model M is calibrated based on {X , Y}. In this paper, PC-Kriging meta-models are used, whereas ordinary Kriging meta-models are used in [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF].

3. A large set of candidate samples S = {x 1 , ..., x n } is generated from X and the response values are predicted by the meta-model M.

4. The limit state surface is estimated based on the current meta-model. 5. All candidate samples x ∈ S are ranked according to an enrichment criterion.

6. Then the selection step determines the sample(s) χ * to be added to the experimental design of the meta-model, so that X ← {X , χ * }.

Initial experimental design

In order to have a well-behaved initial experimental design, a space-filling sampling method is required. In this paper the Latin-hypercube sampling algorithm is used. In particular, the unit hyper-cube [0, 1] M is sampled and the samples are mapped to the variable space by the inverse cumulative distribution function of the marginals. Other sampling strategies include quasi-random sequences, such as Sobol' sequences.

Candidate selection

Single sample selection

A learning function (LF) estimates the expected value of information gained with respect to the statistics of interest when point x is added to the experimental design X of the metamodel. A variety of learning functions has been proposed in the literature including the Expected Feasibility Function (EFF) [START_REF] Bichon | Efficient global reliability analysis for nonlinear implicit performance functions[END_REF], the Expected Improvement (EI) functions [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Ginsbourger | Distance-based Kriging relying on proxy simulation for inverse conditioning[END_REF], the Gaussian process upper confidence bounds (GP-UCB) [START_REF] Srinivas | Information-theoretic regret bounds for Gaussian process optimization in the bandit setting[END_REF] and the U -function [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF]. In the following, we discuss the U -function more in depth.

The U -function is based on the concept of misclassification and the very nature of the Gaussian process meta-model. To each sample point x ∈ S corresponds a non-zero probability that the prediction mean µ Y (x) > a (safe domain) whereas the true value satisfies M(x) ≤ a (failure domain) or vice versa. Due to the fact that the prediction in PC-Kriging is a Gaussian random variable described by the first two moments {µ Y (x), σ Y (x)}, the probability of misclassification P m can be written as [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]:

P m (x) = min Φ µ Y (x) -a σ Y (x) , Φ a -µ Y (x) σ Y (x) ≡ Φ - |µ Y (x) -a| σ Y (x) . ( 17
)
The probability of misclassification is maximized when the fraction tends to zero, i.e. -|µ Y (x)-a|/σ Y (x) ∼ 0 → P m (x) = 0.5, and is small when either the prediction mean is far away from the limit state parameter a and/or the prediction variance is small. Note that the probability measure used to define misclassification is the one associated to the Gaussian nature of the predictor N (µ Y (x), σ Y (x)) in each point x. It shall not be confused with the probability P (•) in Eq. ( 2) which corresponds to the input random vector X.

In this formalism, the U -function is defined as the "reliability index" attached to the probability of misclassification [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF]:

U (x) = |µ Y (x) -a| σ Y (x) . ( 18
)
The optimal sample to enrich the experimental design is the one minimizing the U -function among x ∈ S, thus maximizing the probability of misclassification:

x * = arg min x∈S U (x) ≡ arg max x∈S P m (x). ( 19
)
Due to this similarity, both the probability of misclassification and the U -function can be used as learning functions.

Multiple sample selection

In cases where parallel computing is available, it can be beneficial to add multiple points at the same time even to speed up the overall computation even if the procedure may be slightly suboptimal. Assuming that K samples can be determined at the beginning of each iteration, the corresponding model responses could be computed simultaneously on K independent CPUs.

There are several ways to sample multiple points from the set of candidates S. Generally speaking though, attractive candidates are found close to the limit state surface. Thus we define a lower and upper boundary of the limit state surface defined by µ Y (x) = a which takes into account the prediction uncertainty in the Kriging model. The lower boundary reads:

µ Y (x) -k • σ Y (x) = a, (20)
and the upper boundary:

µ Y (x) + k • σ Y (x) = a, (21)
where k sets the confidence level typically equal to 1.96 = Φ -1 (97.5%). In other words, e.g., the lower boundary estimates the limit state surface assuming that the real value of every

sample x ∈ X is µ Y (x) -k • σ Y (x) instead of the mean value µ Y (x).
Analogously, we define the "mean" failure domain in terms of the prediction mean value:

D 0 f def = x ∈ D X : µ Y (x) ≤ a , (22)
and the corresponding lower and upper bounds to the failure domain:

D - f def = x ∈ D X : µ Y (x) + k • σ Y (x) ≤ a , (23)
D + f def = x ∈ D X : µ Y (x) -k • σ Y (x) ≤ a , (24)
so that D - f ⊂ D 0 f ⊂ D + f . The lower and upper bounds can be interpreted as the least and most conservative estimate of the failure domain, respectively.

The limit state margin M f is defined as the intersection between the lower and the upper boundaries of the failure domain [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]:

M f def = D + f ∩ D - f . (25)
The limit state margin is a natural region where to focus for candidate samples for the design enrichment. Considering a large set of samples of the input vector S = {x 1 , . . . , x n }, we define the following sets corresponding to the domains in Eq. ( 22)-( 25):

S 0 f def = x ∈ S : µ Y (x) ≤ a , (26)
S + f def = x ∈ S : µ Y (x) -k • σ Y (x) ≤ a , (27)
S - f def = x ∈ S : µ Y (x) + k • σ Y (x) ≤ a , (28)
S M def = S + f ∩ S - f , (29)
where S - f ⊂ S 0 f ⊂ S + f . These sets are the discrete corresponding counterparts to the domains defined in Eq. ( 22)-( 25). The margin set S M contains points of interest for enriching the experimental design since they lie close to the true limit state surface. Virtually any point in S M could be added to X and a simple approach would be to sample the K different points randomly. However a better coverage of M f can be obtained using clustering techniques.

To account for the relative importance of the samples in S M , a weighted K-means clustering algorithm is used [START_REF] Zaki | Data Mining and Analysis: fundamental Concepts and Algorithms[END_REF]. In this algorithm samples with high information value according to the learning function have larger weights, whereas samples with low information value have low weights. The weights are set equal to the probability of misclassification (see Eq. ( 17)) of each sample x ∈ S M . By definition of the probability of misclassification, the weights are bounded on [0, 0.5]. The additional K samples χ (i) are then determined as the samples x ∈ S M closest to the nuclei of each ith cluster.

To showcase the process of selecting additional samples, we consider the following function with two-dimensional input vector:

y = M(x) = 20 -(x 1 -x 2 ) 2 -8 • (x 1 + x 2 -4) 3 , ( 30
)
where the input variables have uniform distributions X i ∼ U(-5, 5), i = 1, 2. As an example we compute the failure probability of the function associated with the criterion {M(x) ≤ 0},

i.e. a = 0.

The function is meta-modeled with an ordinary Kriging model. Figure 2 It can be observed that the K = 5 samples in Figure 2(f) cover the domain of the limit state margin well and are close to the limit state surface determined by the prediction mean value.

This leads to an efficient selection of additional samples for enriching the experimental design.

Interestingly, however, the optimal sample in terms of the probability of misclassification is not part of the K = 5 sample set. This implies that each of the K = 5 samples is suboptimal in terms of the probability of misclassification.

The optimal sample (filled diamond marker in Figure 2(f)) depends highly on the sampling of the limit state margin. In principle any sample with µ Y (x) = 0 could be optimal in terms of the probability of misclassification. This implies that the choice of a single sample with respect to the maximal probability of misclassification is non-unique in theory. However it is unique and well-defined in the discrete optimization problem Eq. ( 19).

Stopping criterion

In [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF], the convergence measure (a.k.a. stopping criterion) is based on the accuracy of the meta-model around the limit state surface rather than on the estimation of the statistics of interest. Thus, we propose two stopping criteria for the estimation of failure probabilities and quantiles designed to maximize the accuracy of the statistics of interest while minimizing computational costs.

The stability of the estimate of the statistics of interest can be measured by the size of the limit state margin M f and consequently by the associated values of the upper and lower boundaries of the limit state surface. When the boundaries are close to each other (i.e. a small limit state margin), then the estimate of the statistics of interest is accurate. Therefore we define the stopping criterion for estimating failure probabilities as:

P + f -P - f P 0 f ≤ P f , (31)
for two consecutive iteration steps where P f = 5% is used in applications. The upper and lower bound failure probabilities are defined as:

P ± f def = P µ Y (X) ∓ k σ Y (X) ≤ a (32)
and:

P 0 f def = P µ Y (X) ≤ a . (33)
The stopping criterion defined in Eq. ( 31) can be adapted to the context of quantile estimation as follows:

q + α -q - α Σ Y ≤ qα , (34)
for two consecutive iteration steps where qα = 5% practical applications. The quantiles q ± α are computed from a large Monte Carlo sampling of the bounds µ Y (x)±k σ Y (x), where Σ Y is the empirical standard deviation of the sample Y val = µ Y (x i ), i = 1, . . . , n . The normalization in Eq. ( 34) makes use of the standard deviation Σ Y rather than by q 0 α in analogy with Eq. ( 31). Indeed, unlike P 0 f in Eq. ( 31), q 0 α can be any real number in R, which makes it inappropriate for normalization.

Applications

Four-branch function

Problem statement

The four-branch function is a common benchmark in structural reliability analysis that describes the failure of a series system with four distinct component limit states. Its mathematical formulation reads [START_REF] Waarts | Structural reliability using finite element methods: an appraisal of DARS: Directional Adaptive Response Surface Sampling[END_REF]Schueremans and Van Gemert, 2005a,b):

f 1 (x) = min                3 + 0.1 (x 1 -x 2 ) 2 -x 1 +x 2 √ 2 3 + 0.1 (x 1 -x 2 ) 2 + x 1 +x 2 √ 2 (x 1 -x 2 ) + 6 √ 2 (x 2 -x 1 ) + 6 √ 2                (35)
where the input variables are modeled by two independent Gaussian random variables X i = N (0, 1).

The failure event is defined as f 1 (x) ≤ 0, i.e. the failure probability is P f = P (f 1 (X) ≤ 0).

Failure probability estimation

Setting The adaptive experimental design algorithm is initiated with N 0 = 12 Latinhypercube samples (LHS) and a candidate Monte Carlo population S of n = 10 6 samples.

The performances of PC-Kriging (PCK) are compared to those of ordinary Kriging (OK).

The trend in PC-Kriging consists of Hermite polynomials of maximal degree of eight and a Gaussian autocorrelation function is chosen for the Gaussian process. Both meta-modeling techniques are calibrated using the Matlab-based toolbox UQLab [START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF].

The adaptive experimental design algorithm is tested in both single and multiple (K = 6) samples mode.

Visualization Figure 3 and 4 visualize several iterations of the algorithm for OK and PCK, respectively, for single sample enrichment. The gray dots represent the candidate MC population S. The empty squares mark the initial experimental design, whereas the blue filled circles mark the additional samples. The solid black line represents the exact limit state surface.

Both variants discover the four distinct failure modes in the computational model within the first 50 iterations. The main difference lies in the limit state exploration pattern. When using ordinary Kriging in Figure 3, the distinct failure mechanisms are explored one-by-one. As soon as a failure mechanism is characterized, the adaptive sampling algorithm moves to the next. At iteration 50 all four failure mechanisms are discovered. In the case of PC-Kriging, all four failure mechanisms are already discovered after 10 iterations (Figure 4(a)). Thus, PC-Kriging converges on average faster to an accurate estimate of the failure probability as compared to ordinary Kriging. This is to be expected due to the more accurate global trend in PC-Kriging. Figure 5 illustrates the convergence of the estimate of the failure probability by showing P ± f (k = 2) and P 0 f . Further, the final iteration is marked where the iterative algorithm is stopped due to the stopping criterion in Eq. (31). Figure 5 confirms the faster convergence of the PC-Kriging meta-models. Note that for the Kriging meta-model, Figure 5 shows four plateaus in P 0 f which correspond to the exploration of the four branches of the limit state function.

In order to test the statistical significance of the convergence, the analysis is replicated 50 times with different initial experimental designs. Figure 6 illustrates the convergence of the failure probability estimates P 0 f based on ordinary Kriging and PC-Kriging up to iteration 100, respectively. The solid line represents the mean value of 50 independent replications of the same analysis with different initial Latin-hypercube experimental designs. The dashed lines represent the 5% and 95% quantiles of the 50 independent replications, i.e. the 90% confidence interval (CI). The results show that PC-Kriging converges faster than ordinary Kriging. Note that the variation in P f after convergence (iterations ≥ 80) originates from the finite size of S and the corresponding variance in Monte Carlo simulations (see also Eq. ( 5)).

Parametric study In this paragraph we discuss the candidate selection algorithm and the choice of stopping criterion. The candidate selection algorithm compares the selection of the The full failure probability analysis is replicated 50 times with different initial LHS experimental design in order to assess the statistical significance (e.g. standard deviation of the resulting P f ). The analysis is focused on the total number of computational model evaluations and on the accuracy of the estimate of the failure probability. The total number of computational model runs is N tot = N 0 + N 1 , where N 1 is the number of samples added by the adaptive algorithm.

Table 1 summarizes the results for the four-branch function. The results of the adaptive experimental design algorithm are presented for both stopping criteria (Eq. (31) and Eq. ( 36)) and meta-models (ordinary Kriging and PC-Kriging). The reference value for the failure probability is computed by Monte Carlo simulation with 10 8 samples.

Considering the cases of adding a single candidate, all configurations estimate accurately The results in Table 1 show that the stopping criterion in Eq. ( 31) leads to accurate results despite the smaller number of computational model runs N tot when compared to the stopping criterion in Eq. ( 36). This indicates that Eq. ( 31) is more suitable to the estimation of failure probabilities than Eq. (36). It is interesting though that PC-Kriging converges faster than ordinary Kriging with the stopping criterion in Eq. ( 31), as opposed to the one in Eq. ( 36). The Table 2 summarizes the results and compares the estimate of the quantiles q α with a reference solution from a Monte Carlo simulation (n M C = 10 8 ). The three α values are accurately estimated with the adaptive algorithm despite the small initial experimental design of N 0 = 12 samples. The coefficient of variation of the results increases with decreasing α due to the decreasing size of the failure domain.

It is interesting to notice that the smaller the α value is, the fewer samples are required to estimate the quantile. This phenomenon can be explained by considering the failure domain defined by Eq. ( 35). For small values of α, indeed the failure domain is also small. Hence to meta-model the limit state surface, a smaller number of experimental design points is required.

α MC E [ q α ] CoV [ q α ] E [N tot ] CoV [N 1 ] 0.

Borehole model

Definition

The first realistic engineering problem considered is the so-called borehole-function, which describes the water flow through a borehole. This benchmark function has been discussed in papers such as [START_REF] Harper | Sensitivity/uncertainty analysis of a borehole scenario comparing Latin hypercube sampling and deterministic sensitivity approaches[END_REF], [START_REF] Morris | Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[END_REF], An and Owen (2001), [START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions -Application to uncertainty analysis in computational dosimetry[END_REF]. It is a fast-to-evaluate function [START_REF] Harper | Sensitivity/uncertainty analysis of a borehole scenario comparing Latin hypercube sampling and deterministic sensitivity approaches[END_REF] depending on an eight-dimensional input vector x = [r w , r, T u , H u , T l , H l , L, K w ] T :

v(x) = 2πT u (H u -H l )
ln(r/r w ) 1 +

Conclusions

In this paper, a unified framework for the estimation of both failure probabilities and quantiles is introduced in the context of rare event estimation. The combination of PC-Kriging and design enrichment (adaptive experimental design) based on the probability of misclassification increases the accuracy in the estimation of the statistics of interest.

An efficient strategy to add multiple samples at each iteration of the adaptive experimental design algorithm is introduced to enable the use of available high-performance computing resources.

The adaptive experimental design algorithm is also equipped with a new stopping criterion which monitors the convergence of the statistics of interest better than existing stopping criteria. This results in further reduction of the computational resources needed.

The performance of the proposed algorithm is assessed and illustrated through a benchmark analytical function, which shows the algorithm's capabilities in a variety of settings. The results show that the proposed strategy provides an efficient way to estimate failure probabilities as well as quantiles. Two engineering applications (the water flow through a borehole and the deflection of a truss structure modeled by finite elements) confirm the strengths of the adaptive experimental design algorithm in the context of realistic high-dimensional engineering problem settings, and thus its significance for the engineering practice.

  is a trend defined by a set of |T | multivariate orthonormal polynomials ψ τ (x) indexed by the multi-index τ = {τ 1 , . . . , τ M } and t τ are the corresponding coefficients.

Figure 1 :

 1 Figure 1: Flowchart of the adaptive algorithm for rare event estimation

  (a) displays the initial experimental design ("•" for failure samples and "+" for safe samples), the exact failure domain (thin line) and the failure region estimated from the prediction mean values µ Y (x) (thick black line). Figures 2(b)-2(d) show the lower and upper bounds of the failure domain and the domain of the limit state margin, respectively. Finally, Figure 2(e) shows the probability of misclassification in the limit state margin and Figure 2(f) shows the selected candidates. The gray area in Figure 2(f) represents the samples in the limit state margin, the large black diamond represents the optimal sample (with highest probability of misclassification) and the hollow black diamonds represent the K = 5 samples obtained by weighted K-means clustering.

Figure 2 :

 2 Figure 2: Design enrichment when adding K = 5 samples -definition of failure domains and illustration of the candidate selection algorithm.

Figure 3 :

 3 Figure 3: Four-branch function -iterations of the adaptive experimental design algorithm with ordinary Kriging meta-models.

Figure 4 :Figure 5 :

 45 Figure 4: Four-branch function -iterations of the adaptive experimental design algorithm with PC-Kriging meta-models.

Figure 6 :

 6 Figure 6: Four-branch function -convergence curves of the adaptive experimental design algorithms (statistics of 50 replications of the analysis with different initial Latin-hypercube experimental design).

  x) is the fluid water flow measured in m 3 /year, r w is the radius of the borehole, r the radius of influence, T u the transmissivity of the upper aquifer, H u the potentiometric head of

Table 1 :

 1 Four-branch function -results using different meta-models, candidate selection algorithms and stopping criteria. Coefficients of variation CoV are computed based on the 50 replications of the analysis with different initial experimental designs.

	Method Enrich. Stop. crit. E P f

Table 2 :

 2 Four-branch function -quantile estimation results. MC is based n M C = 10 8 samples.

	Coefficients of variation CoV are computed based on 50 replications of the analysis with
	different initial experimental designs.

Table 6 :

 6 Truss structure -results of the failure probability estimate based on a single run of the adaptive experimental design algorithm for different meta-modelling techniques and

	enrichment modes.
	u adm	Method Enrichment P f (CoV [P f ])	β	N tot
	10 cm MC FORM --OK single PCK single OK K = 6 PCK K = 6	10 6 251 12 + 135 = 147 12 + 158 = 170 1.72 12 + 26 • 6 = 168 4.29 • 10 -2 (0.5 %) 1.72 1.91 2.81 • 10 -2 1.71 4.32 • 10 -2 1.71 4.32 • 10 -2 4.31 • 10 -2 4.32 • 10 -2 1.71 12 + 31 • 6 = 198
	12 cm MC FORM --OK single PCK single OK K = 6 PCK K = 6	10 6 236 12 + 164 = 176 12 + 157 = 169 2.96 12 + 27 • 6 = 174 1.55 • 10 -3 (2.5 %) 2.96 3.17 7.57 • 10 -4 2.96 1.53 • 10 -3 2.96 1.52 • 10 -3 1.53 • 10 -3 1.53 • 10 -3 2.96 12 + 25 • 6 = 162
	14 cm MC FORM --OK single PCK single OK K = 6 PCK K = 6	10 6 231 12 + 110 = 122 12 + 63 = 75 3.99 12 + 27 • 6 = 174 3.6 • 10 -5 (16.7 %) 3.97 4.21 1.29 • 10 -5 3.96 3.7 • 10 -5 3.96 3.7 • 10 -5 3.4 • 10 -5 3.2 • 10 -5 4.00 12 + 11 • 6 = 78
	adaptive experimental design search rather than using FORM algorithm.

 31) 4.458 • 10 -3 1.5 % 12 + 14.4 • 6 = 98.4 8.9 % ( * ) single run results from [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF] explanation lies in the fact that the stopping criterion in Eq. ( 31) measures the convergence of the statistics of interest whereas Eq. ( 36) assesses the overall accuracy of the approximation of the Kriging surrogate.

For the case of PCK and the stopping criterion in Eq. ( 31), Table 1 includes also the results obtained by adding K = 6 candidates at each iteration. The comparison of the total number of computational model runs N tot to the case of adding a single point in each iteration shows that the addition of multiple points at each iteration requires a larger N tot : indeed in the case of single point enrichment, the total number of model runs is N tot = 73, whereas adding K = 6 points in each step leads to N tot = 98. Considering however that the runs of the computational model can be performed in parallel (K = 6 on six independent CPUs), the total time for the K-point enrichment algorithm is approximately four times smaller.

Quantile estimation

Consider now the quantile estimation problem related to the four-branch function. The quantiles are estimated for α = {0.01, 0.001, 0.0001}. PC-Kriging is used, a single sample is added at each iteration out of a candidate Monte Carlo population of n = 10 6 samples and the algorithm is stopped with the stopping criterion in Eq. ( 34).

the upper aquifer, T l the transmissivity of the lower aquifer, H l the potentiometric head of the lower aquifer, L the length of the borehole and K w the hydraulic conductivity of the soil.

The uncertainties in the input vector are modeled as independent random variables whose properties are summarized in Table 3. For the lognormal distribution, the parameters are the mean and standard deviation of the natural logarithm of the variable. For the other variables, they describe the range of uniform distributions.

Table 3: Borehole model -definition of the probabilistic model of the input variables. For the uniform distributions the parameters denote the range, whereas for the lognormal distribution, the parameters denote the mean and standard deviation of the natural logarithm of the variable. Variable Distribution Parameters

Quantile estimation

The statistics of interest are the quantiles of the water flow v(x) corresponding to α = {0.99, 0.999, 0.9999}.

An initial experimental design of N 0 = 12 Latin-hypercube samples in used. The number of samples in the candidate Monte Carlo population is set to n = 10 6 and the remaining settings are the same as in the four-branch function.

The results of the adaptive algorithm and the reference values are reported in Table 4 for the various values of α, the meta-modeling method (PCK and OK) and enrichment strategies (single or multiple) used.

The prediction of the quantiles is accurate for both meta-modeling techniques (OK and PCK) requires more model evaluations than PC-Kriging (PCK) for all values of α. The addition of K = 6 samples at a time induces slightly more model evaluations N tot but at the same time reduces the total computational time due to the parallel evaluation of K = 6 samples.

Two-dimensional truss structure

Definition

For the purpose of illustrating a realistic structural engineering application, a two-dimensional truss structure is analyzed. This specific truss example structure has been presented and discussed previously in [START_REF] Lee | Response surface augmented moment method for efficient reliability analysis[END_REF], [START_REF] Sudret | Quasi random numbers in stochastic finite element analysis -Application to global sensitivity analysis[END_REF], [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF] and [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF].

Consider the simply supported two-dimensional truss structure sketched in Figure 7 consisting of 23 bars and 13 nodes. The geometry is known deterministically whereas the material properties and the loadings are modeled stochastically. Ten stochastic input variables form the input vector X:

where E 1 , E 2 are the Young's moduli of the linear elastic material, A 1 , A 2 are the cross-sections of the bars and P 1 , . . . , P 6 are the vertical loads acting on the nodes of the upper part of the structure. The horizontal bars have the properties {E 1 , A 1 } whereas the diagonal bars have the properties {E 2 , A 2 }. The input variables are modeled by the probability distributions summarized in Table 5. It is assumed that the input variables are statistically independent.

A finite element model of the structure is used to calculate the mid-span deflection, denoted by u(x), as a function of the ten variables in the input vector x, which is defined as positive in the direction indicated in Figure 7.

P 1 , . . . , P 6 [N] Gumbel 5.0 • 10 4 7.5 • 10 3

Failure probability estimation

The failure probability is defined as the probability that the mid-span deflection u is larger than an admissible value u adm , i.e. P f = P (u admu(X) ≤ 0).

The adaptive experimental design algorithm is started with an initial experimental design of N 0 = 12 Latin-hypercube samples and the candidate Monte Carlo population S has a size of n = 10 6 samples. A single sample or K = 6 samples are added to the experimental design at each iteration.

Table 6 summarizes the results of the failure probability estimation. The results are compared to a reference plain Monte Carlo simulation with n M C = 10 6 samples and to the first order reliability method (FORM) in terms of failure probability and Hasofer-Lind reliability index [START_REF] Zhang | Two improved algorithms for reliability analysis[END_REF]. For comparison purposes the generalized reliability index β = Φ -1 ( P f ) is also given for the various sampling-based approaches.

The values of P f shown in Table 6 are consistent with the reference for all the adaptive experimental design configurations. The number of runs of the computational model, however, differ considerably. For the cases of u adm = 10 cm and 12 cm the number of exact model runs are similar for both ordinary Kriging (OK) and PC-Kriging (PCK). In the case of u adm = 0.14 cm, PCK requires significantly fewer runs than OK, for both single and multiple sample selection modes. Similarly to the borehole example, multiple sample selection results in more computational model runs, but overall lower computational times (up to 5-times smaller) due to high performance computing.

The well-established FORM gives an estimate worse than the adaptive experimental design algorithm despite the higher cost. It is thus preferable in this application to conduct an