Achille Murangira

Timothée Miguel Munoz Zuniga

M Munoz Zuniga

T Perdrizet

Structural reliability assessment through metamodel based

Keywords: Reliability analysis, Surrogate model, Importance sampling, Sufficient dimension reduction, Kernel dimension reduction

Introduction

Structural reliability analysis aims at assessing the probability of occurrence of an extreme event related to a given structure. In the usual setting, we are given a limit state function (l.s.f.) g which describes the safety level of the structure for a given input vector X ∈ R d . The failure domain F, corresponds to the set of inputs for which the performance function g is negative, i.e. F = {X ∈ R d | g(X) ≤ 0}, and the failure probability is thus:

P f = P(g(X) ≤ 0) = E(1 g(X)≤0 ) (1)
This probability may readily be estimated through standard Monte Carlo simulation, however since P f is often less than 10 -5 it takes roughly 10 7 -10 8 evaluations of G to obtain an estimate with a coefficient of variation (c.o.v.) less than 5 %. This is clearly problematic for many engineering applications since most of the time G depends on the output of an expensive computer model. Analytical approximations based on the most probable failure point have long been the practical alternative to Monte Carlo sampling. This has led to the popular FORM/SORM approximations [START_REF] Ditlevsen | Structural Reliability Methods[END_REF]. However these approximations are often poor whenever there exists multiple failure modes or when the non-linearity are extreme. Moreover, the FORM/SORM methods do not provide any confidence bounds on the failure probability estimate. Importance sampling (IS) [START_REF] Hammersley | Monte Carlo Methods[END_REF]) is a standard variance reduction Monte Carlo method that can * Corresponding author be used to estimate P f by sampling from an importance density whose main contribution takes place near the limit state surface. A preliminary design point or multiple failure mode computation can be used to design an importance density such as a mixture of standard distributions centred on the failure modes [START_REF] Au | Reliability of uncertain dynamical systems with multiple design points[END_REF] or a single design point [START_REF] Melchers | Importance sampling in structural systems[END_REF]. Still, the number of evaluations of the limit state function required to reach a given accuracy may still be important, limiting the applicability of such a scheme. Subset simulation [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF] eliminates the need to design an importance density by estimating the failure probability as a product of intermediate probabilities which are evaluated by Monte Carlo Markov Chain simulations. However, the incurred computing cost is usually prohibitive in many industrial reliability cases. As a result, methods based on surrogate modelling have been devised to limit the number of computer model evaluations. AK-MCS [START_REF] Echard | AK-MCS: An active learning reliability method combining kriging and monte carlo simulation[END_REF], and AK-IS [START_REF] Echard | A combined importance sampling and kriging reliability method for small failure probabilities with time-demanding numerical models[END_REF] are based on the active kriging principle, that is the performance function is replaced by a kriging metamodel which is iteratively refined so as to provide accurate kriging predictions in the vicinity of the limit state surface.

While often superior to a simple FORM analysis, simple substitution of G by the metamodel gives generally non-consistent failure probability estimators. In the subset simulation framework, a similar approach has been proposed by [START_REF] Bourinet | Assessing small failure probabilities by combined subset simulation and Support Vector Machines[END_REF], where a SVM classifier is built to emulate the intermediate limit state surfaces.

To obtain consistent importance sampling estimators, the Meta-IS [START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF] method was introduced, where the sought probability is evaluated via an approximation of the optimal importance density (OID) based on a kriging emulator. The corresponding IS estimator of P f is cast as a product of the failure probability relative to the metamodel -the augmented failure probability -and a correction factor based on the limit state function G. Simulations on mechanical reliability problems have demonstrated the efficiency of this methodology both in terms of accuracy and simulation budget. However, it is known that fitting a kriging metamodel can break down in high dimensional settings (d ≥ 20 -30), or demand an unreasonable amount of code evaluations to obtain sufficient accuracy in the failure region . Furthermore, in the methodology of Meta-IS, sampling from the approximate OID is accomplished by resorting to an MCMC algorithm, most of which are known to suffer from the curse of dimensionality. On the other hand dimension reduction techniques have been the object of much research in computational statistics in the last couple of decades. In this paper, we lay out a methodology for importance sampling based on the sufficient dimension reduction (SDR) framework [START_REF] Chiaromonte | Sufficient dimension reduction and graphics in regression[END_REF] for regression. SDR is based on the assumption that the output statistical dependency on the input X can be described entirely by projecting X on a lower dimensional subspace. Based on this assumption, we propose to build a kriging metamodel in the dimension reduction subspace, making the Meta-IS algorithm tractable. Assuming that the input vector is Gaussian, we show that the augmented failure probability can be estimated through Monte Carlo simulation in the dimension reduction subspace. As for the correction factor which is based on a sample from an approximation to the optimal importance density, we demonstrate how sample generation can be achieved efficiently since MCMC sampling is only performed in the dimension reduction subspace. We also discuss the case of application of SDR for Gaussian process emulation of complex models with non-Gaussian inputs. The paper is organized as follows. Section 2 recalls the basics of metamodelbased importance sampling. Section 3 gives an overview of dimension reduction tools, more specifically Kernel Dimension Reduction (KDR) [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF] and gradient-based KDR [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]. Section 4 is devoted to metamodel based importance sampling with sufficient dimension reduction. Finally, in section 5 we numerically illustrate these ideas on an academic example and on a relevant industrial case study which consists in the assessment of the structural reliability of a wind turbine, where the wind and wave processes are are modelled as Gaussian processes leading to a high dimensional (d ≥ 100) case study.

Importance sampling with a kriging metamodel

Throughout this article, the input is a real d dimensional random vector X = (X 1 , . . . , X d ) T with density q . As previously stated g denotes the limit state function.

Kriging models, also known as Gaussian process regression models, are flexible and efficient surrogates to complex computer codes. As most metamodels, they rely on an initial design of experiments (DoE) y * = x 1 , . . . , x N D . The popularity of Kriging stems among other things from the availability of prediction uncertainty estimates, given by the Kriging variance. This makes it possible to devise refinement strategies aimed at increasing some measure of accuracy in the region of interest. For reliability purposes this region is in the vicinity of the limit state surface defined by {x ∈ R d | g(x) = 0}.

Gaussian process models

Kriging is based on the assumption that the performance function G is a sample from a Gaussian process (GP) G so that

G(x) = f T (x)β + Z(x) (2)
where f T β is the GP mean and Z a zero-mean stationary Gaussian process.

f = f 1 , . . . , f p T is a vector of basis functions ∈ L 2 (R d , R) and β a constant vector in R d . Z is parametrized by its autocovariance function C(x, x ′ ) = σ 2 G R θ (x -x ′ ) (3)
where σ 2 G is the GP variance and θ is a vector of parameters of the autocorrelation function R. Equation 2 is known as universal kriging (UK) as it allows for a non constant mean of the process G. Prediction at a previously unobserved input x is based on the best linear unbiased prediction (BLUP) of G(x) given the observations y = g(x 1 , . . . , x N D ) at the DoE y * . The BLUP at x, denoted Ĝ(x), is a normal random variable N(m Ĝ, σ 2 Ĝ) where

m Ĝ(x) = f T (x) β + r T (x)R -1 (y -F β) (4) σ 2 Ĝ (x) = σ 2 G (1 -r T (x)R -1 r(x) + v T (x)(F T R -1 F) -1 v(x)) (5)
are respectively the prediction mean and variance, R is the correlation matrix of the DoE defined by

R i j = R θ (x i , x j ), i, j = 1, . . . , N D . r = R θ (x -x i ) N i=1
is the cross-correlation vector between the prediction and the observations while F is the matrix defined by

F i j = f j (x i ), 1 ≤ i ≤ N D , 1 ≤ j ≤ p. The vector β is the solution to a generalized least-squares problem β = (F T R -1 F) -1 F T R -1 y (6) and v(x) = F T R -1 r(x) -f(x).

Basics of metamodel based importance sampling (Meta-IS)

In this section, we give an overview of the Meta-IS algorithm originally due to [START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF]. The method is founded on the use of a quasi-optimal importance density which mimics the intractable optimal importance density.

Standard importance sampling

Importance sampling is a well known variance reduction method where the quantity of interest is an expectation of a integrable function. In the context of reliability methods, the expectation of interest is

P f = E q (1 g(X)≤0 ) = R d 1 g(x)≤0 q(x)dx.
Let q be a density such that:

1 g(x)≤0 q(x) 0 =⇒ q(x) 0 (7) E q 1 g(X)≤0 q(X) q(X) < ∞ (8)
IS stems from the equality

P f = E q (1 g(X)≤0 ) = E q 1 g(X)≤0 q(X)
q(X) . Given an i.i.d. sample x (i) , i = 1, . . . , n from q, the importance sampling estimator of P f reads

PIS f = 1 n n i=1 1 g(x (i) )≤0 q(x (i) ) q(x (i) ) (9)
The above estimator is unbiased and its variance is given by

Var q( PIS f ) = 1 n E q 1 g(X)≤0 q 2 (X) q2 (X) - P 2 f n (10)
which can readily be estimated by [START_REF] Rubinstein | Simulation of stochastic processes by spectral representation[END_REF], for the following optimal proposal density qopt (x) = 1 g(x)≤0 q(x) P f (12) which is generally difficult to sample from because it depends on the complex limit state function g. In Meta-IS samples are drawn from a quasi-optimal density obtained by replacing the indicator function with a a continuous probabilistic classification function linked to the kriging predictor.

σ2 IS = 1 n n i=1 1 g(x (i) )≤0 q 2 (x (i) ) q2 (x (i) ) - ( PIS f ) 2 n (11) An important result is that Var q( PIS f ) = 0

Probabilistic classification function and quasi-optimal importance density

Since the Kriging prediction Ĝ(x) at an unknown location x is a normal random variable, it's straightforward to compute the following probabilistic classification function

π(x) = P 0 ( Ĝ(x) ≤ 0) ( 13
)
where P 0 is the probability measure associated with the Gaussian process G. It turns out that

π(x) =        Φ(- m Ĝ (x) σ Ĝ (x) ) if x y * 1 G(x)≤0 if x ∈ y * ( 14
)
where Φ is the standard normal cumulative distribution function (cdf). Now by swapping π(x) and 1 G(x)≤0 in (12), one obtains the following quasi-optimal importance density [START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF]:

q * (x) = π(x)q(x) P f,ǫ (15)
where P f,ǫ = π(x)q(x)dx is the augmented failure probability.

One would be tempted to use the simpler density qmeta (x) = 1 m Ĝ (x)≤0 q(x) P f but this choice does not guarantee that the proposal density dominates the unnormalized pdf 1 g(x)≤0 q(x) as required by condition (7). Note that this is condition is always satisfied for q * no matter the Kriging model.

The Meta-IS failure probability estimator

Having defined the quasi-optimal instrumental density, the failure probability estimate may now be broken down as follows Dubourg et al. ( 2011)

P IS f = 1 g(x)≤0 q(x) q * (x) q * (x)dx = P f,ǫ 1 g(x)≤0 π(x) q * (x)dx (16) = P f,ǫ α corr ( 17
)
where

α corr = E q * ( 1 g(X)≤0 π(X)
) is a factor that corrects for the bias in the augmented probability estimator which relies solely on the metamodel. The value α corr actually depends on how accurate the Kriging predictor is in the failure domain F. It's clear that whenever the metamodel is highly accurate in the failure domain, the ratio

1 g(X)≤0
π(X) is close to 1 and so is α corr . In case of lower accuracy, there are regions of F with high uncertainty on the sign of the Kriging prediction, i.e. π(x) is far from 1. This yields a correction factor also different from 1. Let x (i) , i = 1, . . . , N meta be an i.i.d. sample from the prior density q, then

P f,ǫ = 1 N meta N meta i=1 π(x (i) ) ( 18
)
is an unbiased and consistent estimator of P f,ǫ . The estimation error is quantified by the Monte Carlo variance

σ2 ǫ = 1 N meta (N meta -1) N meta i=1 π(x (i) ) -πN meta 2 (19)
where πN meta , γN corr is the respective sample average. From a practical standpoint, sampling from q is usually straightforward whereas obtaining draws from q * requires adequate algorithms such as Monte Carlo Markov Chain simulation. Let t ( j) , j = 1, . . . be a chain with stationary distribution q * , then from the ergodic theorem 1 N N j=1

1 g(t ( j) )≤0 π(t ( j) ) converges a.s. to α corr as N → ∞. A biased finite sample estimate of α corr is then αcorr = 1 N corr N corr +b j=b+1 1 g(t ( j) )≤0 π(t ( j) ) ( 20
)
where b is a burnin parameter. The standard squared error can be assessed as

σ2 corr = 1 N corr         1 N corr N corr +b j=b+1 h(t j ) -αcorr 2         (1 + γcorr ) (21)
where h(t j ) =

1 g(t ( j) )≤0 π 2 (t ( j) )
and γ corr is an estimator of 2 ∞ k=0 Corr(h(t (0) ), h(t (k) )), which can be estimated as detailed in [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]. In the original paper [START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF], samples from the quasi-optimal density are obtained through a modified Metropolis-Hastings scheme [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF].

To mitigate the inflation of the estimation variance due to the dependence in the chain, thinning can be performed prior to the correction factor estimation, that is taking one in k draws in the simulated chain. This of course increases the length of the simulated chain which must be kN corr long in order to retain N corr samples. From a computational perspective, the MCMC sampler only evaluates the probabilistic classification function, and is therefore expected to be relatively efficient.

2.2.4. Efficient splitting estimator for the augmented failure probability While the estimation of the augmented failure probability P f,ǫ only resorts to the Kriging predictor, its computation in a high dimensional through standard Monte Carlo simulation can still be improved. One idea due to [START_REF] Sudret | Enhancing meta-model-based importance sampling by subset simulation[END_REF] is to use a splitting estimator of P f,ǫ , much like in subset simulation which is base on the following equality:

P f,ǫ = P ǫ,1 s k=2 P ǫ,k|k-1 (22)
where

P ǫ,1 = E q (π 1 (X)) (23) P ǫ,k|k-1 = E qk * π (k) (X) π (k-1) (X) (24)
where π (k) (x) = Φ(-

l k -m Ĝ (x)
σ Ĝ (x) ) and qk * (x) ∝ π (k-1) (x)q(x) is the quasi-optimal density for a failure level l k-1 . The thresholds l 1 > . . . > l s = 0 should be such that the intermediate probabilities P ǫ,k|k-1 are easy to estimate (around 0.1-0.3 for instance). They can be fixed a priori or by running a subset simulation algorithms on the Kriging mean of the limit state function and using the corresponding intermediate levels as the l k thresholds. Notice that P ǫ,1 can evaluated by standard Monte Carlo efficiently since it is a relatively large probability and q. Estimation of P ǫ,k|k-1 and its standard error σ ǫ,k|k-1 can be achieved by running a Monte Carlo Markov Chain targeting qk * . Assuming δǫ,k = σǫ,k|k-1 Pǫ,k|k-1 is the coefficient of variation for the estimation of P ǫ,k|k-1 , a conservative estimate of the overall c

.o.v. of P f,ǫ is δǫ = s i=1 δ 2 ǫ,k (25)

Kriging metamodel refinement

Before computing both quantities P f,ǫ and α corr , it is paramount to have a sufficiently accurate Kriging metamodel. Otherwise, the approximation q * of the optimal importance density would be poor and yield high estimation variance. To obtain an accurate emulator, one has to iteratively enrich the initial design DoE with new limit state function evaluations. This initial DoE can be any relevant space filling design. In reliability analysis, iterative refinement methods are usually geared towards accurate approximation of a target region as evidenced by [START_REF] Vasquez | A sequential bayesian algorithm to estimate a probability of failure[END_REF], [START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF] or [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]. The sign of the performance function G which delineates the failure domain is the most uncertain in this target region. It is usually determined by optimizing some in-fill criterion [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]. However, the in-fill criterion might feature multiple optima or cause the optimizer to be stuck on a local plateau. Therefore, unless one uses a multi-start scheme (which wouldn't guarantee finding all local optima anyway), most iterative in-fill criteria do not allow for multiple points to be added to the DoE. This has motivated sampling based design enhancements as in [START_REF] Bourinet | Assessing small failure probabilities by combined subset simulation and Support Vector Machines[END_REF] which rely on a so-called margin shrinking concept: initially these were based on adding multiple points in the margin of a support vector margin (SVM) classifier, which acts as a surrogate to the function G. The idea is to cluster points generated from the density of the inputs given that they belong to the margin. In practice after generating a sufficiently large number of samples (through accept-reject for instance), k-means clustering could be used to obtain any given number of points. The Kriging metamodel enhancement proceeds in a similar fashion by drawing a large number of samples from the weighted margin probability density defined as

h(x) ∝ w(x)q(x) (26)
where

w(x) = P 0 ( Ĝ(x) ∈ [-t 1-α σ Ĝ (x), t 1-α σ Ĝ(x)]) (27) = Φ t 1-α σ Ĝ (x) -m Ĝ(x) σ Ĝ(x) -Φ -t 1-α σ Ĝ(x) -m Ĝ(x) σ Ĝ(x) ( 28
)
is the margin probability [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF] and t 1-α is the 1α quantile of the standard normal pdf. Drawing samples from h yields points with a high uncertainty on the sign of their Kriging predictor: these are good candidates to be added to the DoE. In practice, an MCMC algorithm (e.g. slice sampling) can provide a sufficiently large sample with stationary distribution h. The k-means algorithm then provides N add clusters from these samples, which are evaluated on the performance function and added to the DoE. The Kriging model is then updated on the basis of this enhanced design.

The model refinement is iterated as long as a stopping criterion is not met. In the Meta-IS algorithm, the usual criterion is a leave-one-out estimate of the correction factor which is defined as

α LOO = 1 N D N D i=1 1 G(x i )≤0 π -i (x i ) ( 29
)
where

π -i (x) = Φ - mG D\x i (x i ) σG D\x i (x i )
is the classification function obtained by removing observation x i from the DoE D. As stated in section 2.2.3, a factor close to 1 signifies accuracy of the classifier based on the Kriging classification function π. The DoE enhancement can therefore be terminated whenever 0.1 ≤ α LOO ≤ 10 and the number of experiments N D is higher than a pre-specified threshold (typically a few tenths/hundreds). It's also necessary to impose a maximum of computer model evaluations during the initial space-filling and refinement phase, as there is no guarantee that α LOO will get closer to 1 with a reasonable amounts of l.s.f. evaluations.

Sufficient dimension reduction for regression and classification

3.1. Motivations and principle of sufficient dimension reduction While flexible enough to cover a wide range of reliability analysis problems, the metamodel based importance sampling procedure can sometimes be intractable when dealing with high dimensional inputs. Indeed, the Kriging metamodel usually needs a DoE with sufficient size to cover the input space. Furthermore, the Gaussian process training which requires learning the kernel hyperparameters might be faced with Hessian matrix singularity during marginal likelihood optimization, especially when the kernel is parametrized by a scalar parameter for each input variable (e.g. the length-scale parameters for the squared-exponential (SE) kernel). Attempts have been made recently to propose adequate representations for high dimensional Kriging. For instance, additive Kriging [START_REF] Durrande | Additive covariance kernels for high-dimensional gaussian process modeling[END_REF] can be used by assuming an additive model for the emulator, i.e. that model is a sum of univariate metamodels. The resulting kernel function is simply the sum of the kernel submodels. This was shown to be more effective than using standard separable kernels on a few test cases involving between 10 and 50 variables. Another relevant approach is to incorporate variable selection into the Gaussian process regression framework as suggested by Yi et al. (67).

In this work, we take another approach which is motivated by the literature in supervised dimension reduction, specifically the sufficient dimension reduction (SDR) framework [START_REF] Chiaromonte | Sufficient dimension reduction and graphics in regression[END_REF]. The idea of SDR is to find a low rank matrix B such that regression of the output Y on X can be performed by replacing X by B T X. Without requirements on the rank r of matrix B, B always exists since B = I d does the trick. However, for many regression models, it is frequent for the output to be only a function of a linear combination of inputs as well as an optional noise independent of the inputs. For instance, if Y depends only on a subset of covariates, B is a sparse matrix with zero rows on the irrelevant predictors. Indeed, assuming without loss of generality that

X T = (X T A , X T -A ) where X A (resp. X -A ) is the subset of rel- evant (resp. irrelevant) components, then B = I |A| O d-|A|×|A| sat- isfies the SDR requirements. More formally, S (B) = span(B) is dimension reduction subspace if Y ⊥ ⊥ X | B T X (30)
where B ∈ R d×r and ⊥ ⊥ stands for independence. Equivalently, (30) could be rewritten as

Y = g r (B T X, ε) (31)
where ε is a noise term independent of X and g r an unknown function. An object of interest is the central subspace (CS) which is defined as the intersection of all dimension reduction subspace and therefore has minimal dimension: it exists under mild conditions (see [START_REF] Cook | Regression Graphics[END_REF]. The interest of SDR for metamodeling is obvious : if we could find a matrix B with sufficiently low rank r, a surrogate model of the link function g r can be built. Typically, in structural reliability of offshore structures we might be dealing with over 100 input variables. If r ≤ 10 usual GP metamodels can then be applied to model the response Y = g(X) but in the subspace S B . Because, we have no prior knowledge on the noise term, we may obtain a coarse metamodel by using the simplified model

Y = g r (B T X) (32)
In this case, the surrogate model is fitted to the link function g r .

To deal with the fact that model (32) does not take into account the noise in the original dimension reduction model ( 7), we suggest to use a nugget term in the corresponding metamodel. Let us stress that this is equivalent to Kriging with noisy observations when considering data points outside the DoE. However, for points in the experimental design, the Kriging metamodel with nugget effect interpolates the observations which apparently contradicts the original model ( 31). This shouldn't however be a major problem since the whole objective is to design an emulator that is close enough to the original limit state function. The fact that we are fitting a metamodel not to the true limit state function but to the link function g r in the simplified regression model (32) shouldn't necessarily hamper the metamodel's performance if the regression noise ε is not predominant. We call the importance sampling procedure with SDR subspace metamodelling Meta-ISDR, where DR stands for dimension reduction. It will be outlined in 4.

SDR estimation methods

Let us get back at dimension reduction subspace estimation. The first arguably popular SDR method is sliced inverse regression [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] which finds a basis of the central subspace under restrictive conditions on the marginal distributions of X. Although very efficient computationally, it is difficult to use in the context of metamodelling for reliability purposes. Indeed, we are interested in reducing the dimension for the regression of the output Y in a region that is large enough to part of the failure domain. Therefore, it is necessary to obtain a training set {x i tr } N tr i=1 consisting in a random samples that fill the input space but also populate the failure domain. Because SIR constrains the training set to satisfy restrictive conditions on the marginal distributions of the training sample, it seems difficult to find such a sampling distribution that also explores efficiently both the safety and failure domains. Also, SIR tends to miss some SDR directions in some special cases [START_REF] Cook | SAVE: A method for dimension reduction and graphics in regression[END_REF]. Other methods such as Sliced Average Variance Estimation (SAVE) [START_REF] Cook | SAVE: A method for dimension reduction and graphics in regression[END_REF] and contour regression [START_REF] Li | Contour regression: a general approach to dimension reduction[END_REF] have similar limitations. More recently, kernel methods have lead to the kernel dimension reduction(KDR) algorithm [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF] and its gradient based counterpart (gKDR, [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]): these algorithms impose no strict conditions on the distribution of the input vector X. As such, they lends themselves well to the task at hand since the user has much more latitude in the choice of the sampling distribution. These are therefore preferable methods to learn a dimension reduction subspace in very general settings.

Standard and gradient based kernel dimension reduction 3.3.1. Kernel dimension reduction

The basis of KDR is to express the conditional independence criterion (30) as an optimization problem that reaches its minimum at the dimension reduction matrix. The optimization problem is defined through conditional covariance operators on reproducing kernel Hilbert spaces (RKHSs) that capture conditional independence. From a practical standpoint these operators are estimated from Gram matrices of the data and a specific functional of these Gram matrices is minimized to yield an estimate of a SDR matrix. Let (X, Y) be a random variable on X×Y with distribution probability P XY where (X, B X ) and (Y, B Y ) are measure spaces. Let k X and k Y be measurable positive definite kernels on X and Y respectively, with associated reproducing kernel Hilbert spaces H X and

H Y . Furthermore, assume E(k X (X, X)) < ∞ and E(k Y (Y, Y)) < ∞.
The conditional independence criterion of (30) in KDR relies on the cross-covariance operator Σ YX : H X → H Y defined as [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF]

g, Σ YX f H Y = E f (X) -f )(g(Y) -ḡ) (33)
for all f ∈ H X and g ∈ H Y . The covariance operator Σ XX is defined similarly. The conditional covariance operator Σ YY|X is defined as

Σ YY|X = Σ YY -Σ 1/2 YY V YX V XY Σ 1/2 YY ( 34
)
where V YX , V XY are bounded operators such that

Σ YX = Σ 1/2 XX V XY Σ 1/2 YY and Σ XY = Σ 1/2 YY V YX Σ 1/2
XX and the conditional crosscovariance is sometimes abusively expressed as

Σ YY|X = Σ YY - Σ 1/2 YX Σ -1 XX Σ 1/2 XY .
Let B ∈ R d×r a matrix such that span(B) is a dimension reduction subspace and B T B = I r . Letting k d a positive kernel and k B X a positive definite kernel on R r defined by k

B X (x, x ′ ) = k d (B T x, B T x ′ ), define the (cross-)covariance opera- tors w.r.t. kernel k B X as Σ B YX , Σ B XY , Σ B YY and Σ B XX .
The conditional cross covariance operator Σ B YY|X associated to k B X is then defined as in (34). Under mild conditions that are satisfied in particular if all involved kernels are Gaussian RBFs, the following fundamental relationship links the dimension reduction subspace and conditional cross-covariance operators [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF].

Σ YY|X = Σ B YY|X ⇔ Y ⊥ ⊥ X | B T X (35)
In addition to (35), we have Σ B YY|X ≥ Σ YY|X (with the order of self-adjoint operators). Thus, in order to find the SDR matrix B, the KDR algorithm minimizes Tr(Σ B YY|X ) subject to B T B = I r . The practical algorithm takes as inputs n i.i.d. samples (x i , y i ) from P XY and solves min

B T B=I r Tr G Y (G B X + nε n I n ) -1 (36)
where G Y and G X are the centred Gram matrices defined by

G X = HK X H, G Y = HK Y H, (K X ) i j = k X (x i , x j ), (K Y ) i j = k X (y i , y j ), 1 ≤ i, j ≤ n and H = I n -11 T n .
The parameter ε n is a regularization term that facilitates matrix inversion. In KDR this non-convex optimization problem is solved by a steepestdescent with algorithm line search (see [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF] for more details).

Gradient-based kernel dimension reduction(gKDR)

The previous kernel dimension reduction procedure is efficient and consistent under non-restrictive assumptions on the joint and marginal distributions of X and Y. However, the required non-convex optimization step makes it somewhat CPU intensive depending on the size of the training samples (x i , y i ). Moreover, kernel parameter selection and the choice of the regularization term is needed and cross-validation for parameter selection is not possible because of the computing cost of the optimization procedure. The gKDR method [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]) solves these issues with a much faster algorithm for dimension reduction subspace identification. It relies on the fact that the derivative ∂E(Y | X = x)/∂x is contained in the SDR subspace. A non-parametric estimator of this derivative is obtained through covariance operators and the dimension reduction matrix B is identified as the solution of an eigenproblem with few matrices inversions. This fast procedure enables regularization and kernel parameter selection through cross-validation and can be applied to large training sets in high dimensions. For a detailed description of gKDR, we refer the reader to the original paper by [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]. Note that both KDR and gKDR require a priori knowledge of the SDR subspace dimension r which is seldom the case in practical applications. In the case of gKDR, we suggest a crossvalidation procedure in 1 based on the mean-squared error of the regression estimate of Y w.r.t. B T r X where B r is the dimension reduction matrix of rank r.

Algorithm 1 CV procedure for dimension selection Require: samples (x i , y i ), i = 1, . . . , n, K cv : number of CV folds, r max : maximum SDR subspace dimension, kernels k X ,k Y Split dataset into K cv folds of similar size C 1 , . . . , C K cv for r = 1 to r max do for j = 1 to K cv do -Estimate Br using (x i , y i ), i C j -Estimate the mean-squared error e 2 r, j of k-NN regression of Y on BT r X using validation fold (x i , y i ), i ∈ C j end for Estimate the regression error of Y on B T r X by e 2 r =

1 K cv K cv j=1 e 2 r, j
end for r ← arg min 1≤r≤r max e 2 r is the SDR dimension estimate Compute B = Br using the whole dataset As a matter of fact, algorithm 1 can be made more efficient since gKDR defines the columns of the dimension reduction subspace matrix estimate B as the eigenvectors corresponding to the r largest eigenvalues of some symmetric n × n matrix (see [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]). This implies that if r < s, then the columns of Br are the first r columns of Bs , possibly up to a sign. It is therefore sufficient to only compute B for the largest candidate dimension for all cross-validation folds leading to the faster algorithm 2.

Algorithm 2 Fast CV procedure for dimension selection Require: samples (x i , y i ), i = 1, . . . , n, K cv : number of CV folds, r max : maximum SDR subspace dimension, kernels

k X ,k Y Split dataset into K cv folds of similar size C 1 , . . . , C K cv for j = 1 to K cv do Compute B( j) = B( j)
r max ∈ R d×r max using (x i , y i ), i C j end for for r = 1 to r max do for j = 1 to K cv do B( j) r ← B( j) :,1:r matrix of first r columns Estimate the mean-squared error e 2 r, j of k-NN regression of Y on B( j) r T X using validation fold (x i , y i ), i ∈ C j end for Estimate the regression error of Y on B T r X by e 2 r =

1 K cv K cv j=1 e 2 r, j
end for r ← arg min 1≤r≤r max e 2 r is the SDR dimension estimate Compute B = Br using the whole dataset

Importance sampling with dimension reduction subspace metamodel

As stated in 3.1, the aim of sufficient dimension reduction is to enable Kriging metamodelling in a space of much lower dimension, given by the dimension reduction subspace. We now consider a Kriging metamodel of the limit state function output Y by using a Gaussian process model GP(f T β, C) for the link function g r of the dimension reduction model in 32. The ensuing Kriging predictor is denoted Ĝr . For reasons that will be clear in the sequel, we restrict ourselves to limit state functions that are defined for random inputs X ∈ R d .

Quasi-optimal importance density with reduced metamodel

Consider the natural extension of the original quasi-optimal importance density defined as qr * (x) = π r (B T x)q(x)

π r (B T x)q(x)dx (37)
where π r is the probabilistic classification function of the Kriging metamodel in the SDR subspace

π r (z) = Φ - m Ĝr (z) σ Ĝr (z) (38)
The ensuing augmented failure probability estimates and correction factor using the IS density qr * then read

P f,ǫ = E q (π r (B T X)) (39) α corr = E qr * 1 g(X)≤0 π r (B T X) (40)

Reduced Kriging metamodel refinements

Regarding the sampling based Kriging refinement strategy outlined in 2.2.5, a similar approach may readily be applied, that is we define the weighted margin probability density

h r (x) ∝ w r (B T x)q(x) (41)
where w r is defined as in ( 28). The refinement strategy using a reduced metamodel may now proceed as follows: given a sample x i , i = 1, . . . , n drawn from a Markov chain targeting h r , N add cluster centers are determined via a k-means algorithm and then projected onto the dimension reduction subspace using matrix B. These projected points can then be added to the DoE. However this runs the risk that two cluster centers have very close projections which is inefficient since the metamodel actually emulates the link function in the simplified regression model ( 32). The same remark can be made in the initial space filling design.

To avoid this, we suggest another strategy. The DoE for the reduced Kriging model is a set of input-evaluation pairs in R r × R in which the evaluated function is the link function g r of the SDR model. Because g r is unknown, we should use the true limit state function g. To do so, a mapping from R r to R d is necessary. The simplified model states that Y = g r (B T X). On the other hand Y = g(X). Hence if (z * , y * = g r (z * )) is a point to be added to the DoE, letting x * = (B T ) # z * where B T (B T ) # = I d ensures that y * = g(x * ) = g r (B T (B T ) # z * ) = g r (z * ). Therefore, (z * , y * ) can be added to the design of experiments by evaluating y * = g((B T ) # z * ). In the case where the input domain X is not R d , there is no guarantee that x * = (B T ) # z * ∈ X, therefore the following tweak can be adopted:

1. let z i = B T x i , i = 1, . . . , n be the projected samples 2. find N add cluster centers {c j } N j i=1 of the z i dataset 3. let i j be the index of the closest projections z i to c j : add (z i j , y j * = g r (z i j )) by evaluating y j * = g(x i j )

The case of Gaussian inputs

The augmented failure probability and correction factor (39) and (40) normally require ordinary Monte Carlo or MCMC sampling in a d-dimensional space. Admittedly, (39) doesn't resort to expensive limit state function evaluations but a small value of the target failure probabibility may require a substantial amounts of sampling in R d and calls to the Kriging predictor. More importantly, the estimation of the correction factor requires sampling from qr * which is achieved through an MCMC algorithm such as the Metropolis-Hastings (MH) algorithm. Tuning the MH algorithm in high dimensions can prove cumbersome though approaches such as the modified Metropolis-Hastings algorithm [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF]) may prove successful. Let us assume the input X to be a standard multivariate Gaussian variable. The estimation of ( 39) and ( 40) is made more efficient as shown below.

Lemma 1. (i) The augmented failure probability P f,ǫ can be expressed as P f,ǫ = E(π r (Z)) where Z is an r dimensional standard normal variable.

(ii) Let B a = [B, B ⊥ ]
where the columns of B ⊥ form an orthonormal basis of span(B)

⊥ . Let W 2 ∼ N(0 d-r×1 , I d-r ), W 1 ∼ p W 1 (w 1 ) = π r (w 1 )ϕ r (w 1 ) P f ,ǫ
where ϕ r is the standard rdimensional multinormal pdf and W

= [W T 1 , W T 2 ] T . Then X = B -T a W is distributed according to qr * (x).
The first result regarding the augmented failure probability computation is a straightforward consequence of the fact that Z = B T X ∼ N(0, I r ) since B T B = I r . The proof of (ii) is deferred to the appendix. The practical implications of these results are that:

• the estimation of P f,ǫ requires sampling of a standard distribution in R r instead of R d .

• the estimation of the correction factor boils down to MCMC sampling in a space of expected much lower dimension r which voids some of the limitations of high dimensional MCMC. The second part of the sampling procedure is a straightforward standard Gaussian generation which makes the whole procedure faster than performing MCMC sampling in R d directly where d ≫ r potentially.

• note that lemma 1 can also be used to sample the weighted margin probability density: the method is strictly the same if one replaces the probabilistic classification function with the margin probability of the reduced Kriging model. Also, given that we only need the projection of cluster centers of the samples onto the dimension reduction subspace, this implies that we only need to simulate a chain with stationary unnormalized distribution w r (w 1 )ϕ r (w 1 ).

The general case

Assuming Gaussian inputs yields a simple importance sampling scheme after linear dimension reduction. In the general case, one might simply consider the limit state function in the standard Gaussian space U, which is denoted g 0 . The U space formulation is classically achieved through an isoprobabilistic transform T : X ∈ X → T (X) = U ∈ U so that g 0 (U) = g T -1 (U) . In this paper, T we consider the Nataf transform defined as T = T 2 • T 1 , where

T 1 : X → Z = T 1 (X) =             Φ -1 (F 1 (X 1 ))
. . .

Φ -1 (F d (X d ))             (42) T 2 : Z → T 2 (Z) = L -1 0 Z
where L 0 L T 0 = R 0 is the linear correlation matrix of Z, Φ is the standard normal cdf and F 1 , . . . , F d are the marginal CDFs of X. One could simply try dimension reduction on the response model Y = g 0 (U), but because of the non-linear nature of the Nataf transform, there is no guarantee that an SDR type hypothesis of the form Y ⊥ ⊥ U | BT U holds for some B even if Y ⊥ ⊥ X | B T X. Nevertheless, the U space formulation being convenient for practical reliability analysis, to achieve dimension reduction, we suggest to simply backtransform to the X space and perform linear dimension reduction. A quasi-optimal importance density in the U space is then obtained as

qr * (u) = π r (B T T -1 (u))ϕ d (u) P f,ǫ (43)
with

P f,ǫ = π r (B T T -1 (u))ϕ d (u)du.
The correction factor reads

α corr = E qr * 1 g 0 (U)≤0 π r (B T T -1 (U)) (44)
Now, because of the non-linear transformation, drawing samples from ( 43) is not as simple as in the case of Gaussian inputs X. A chain targeting qr * has to be simulated in R d either via relevant MCMC method. For this purpose, we suggest using a modified Metropolis-Hastings algorithm [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF].

As far as the sampling based refinement criterion is concerned, we redefine the weighted margin probability in a similar fashion, that is

h r (u) = w r (B T T -1 (u))ϕ d (u) (45)
Since this density may have modes that are far apart, slice sampling is recommended as an alternative to standard Metropolis Hastings which typically struggles in this setting, as noted by [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]. Once samples are obtained, we backtransform them to the physical input domain X and proceed as outlined in 2.2.5.

Alternative probabilistic classification function

A common challenge in estimating the correction factor is that the quasi-optimal density qr * may be multimodal in which case a standard Metropolis-Hastings scheme may be trapped in a mode and fail to provide an accurate sample from the target density. One way to deal with this is to use an ad-hoc MCMC scheme such as population MCMC [START_REF] Liu | Monte Carlo Strategies for Scientific Computing[END_REF]. Instead of changing the sampling scheme, we suggest modifying the quasi-optimal density by altering the probabilistic classification function π r . The modification consists in replacing the standard normal cdf with a logistic function L in the definition of π r :

π r (z) = L - m Ĝr (z) σ Ĝr (z) ( 46
)
where L(x) = 1 1+exp (-x) . With this definition, π r no longer corresponds to the probability that the Kriging predictor falls into the failure domain, however it is a reasonable approximation. The main advantage is that L decreases less slowly to zero as x → ∞, hence the resulting importance density qr * has modes that are separated by valleys where the density achieves larger values. For this choice of IS density, the standard MH scheme has better mixing properties and thus provides more reliable estimates of the correction factor in multimodal settings.

The Meta-ISDR algorithm

We hereafter sketch the Meta-ISDR algorithm, assuming X to be a standard Gaussian vector. The general case may be obtained by considering the modifications stated in 4.4. Note that we assume a preliminary single/multiple failure points search but this is optional. The main parameters of the algorithm, along with sensible default parameters in brackets are:

• x * 1 , . . . , x * m : m ≥ 1 most probable failure points (optional)
• N S DR : number of samples for the SDR basis estimation

• r max : maximum SDR subspace dimension

• K cv (5): number of CV folds for the reduced dimension estimation

• k X ,k Y : Gaussian kernels used in SDR.

• N max : maximum number of limit state function evaluations

• R max : radius of the hypersphere enclosing the design of experiments (10 5 ) number of samples generated for the spacefilling design

• N init D (2d): initial size of DoE • N f ill
• N min D , N max D
minimum/maximum number of points in final design of experiments

• N add : number of points added to DoE during refinement step

• α min LOO (0.1), α max LOO (10): min/max of leave-one-out criteria for the metamodel quality assessment

• N aug : number of samples used for the estimation of the augmented failure probability P f,ǫ

• δ ǫ : target coefficient of variation (c.o.v.) for P f,ǫ

• δ α : target coefficient of variation (c.o.v.) for α corr Although the random input vector is not bounded, the radius parameter R max , which defines the domain of the Kriging metamodel, can be set by choosing R max such that P( X 2 > R 2 max ) is much smaller than P f . For instance if P f >> 10 -b for some positive integer b then solving P(

X 2 > R 2 max ) = 10 -b yields R max = q χ 2 d ,1-10 -b
where q χ 2 d ,1-10 -b is the quantile of order 1 -10 -b of the chi-squared distribution with d degrees of freedom. The complete algorithm is given in table 3.

Numerical illustrations

We now investigate the performance of the metamodel-based importance sampling procedure with dimension reduction on a simple analytical setting and a realistic high dimensional industrial application consisting in estimating the failure probability of a floating wind turbine under stationary wind and wave loads.

Academic example

We consider an analytical limit state function whose dimensionality can be varied to illustrate the impact of increasing number of inputs. The performance function reads

g(X) = d + aσ √ d - d i=1 X i (47)
where the X i are d i.i.d. lognormal variables with mean 1 and standard deviation σ = 0.2 and a = 1. We consider 3 cases with increasing dimension, d = 2, 50, 100. It is obvious that a sufficient dimension reduction subspace is spanned by b

= 1 • • • 1 T since g(X) = d + aσ √ d -b
T X so we can expect sizeable performance improvements for high dimensions for the Meta-ISDR method. The strategy presented here is confronted with standard reliability methods, namely FORM, standard Monte Carlo simulation, importance sampling with a proposal centered on the design point (IS-FORM) and the standard Meta-IS. Note that, with the exception of FORM, we only consider methods that yield consistent probability estimation. The design point in FORM was computed using a sequential quadratic programming algorithm as described in [START_REF] Liu | Optimization algorithms for structural reliability[END_REF] As can be expected for d = 2, the dimension reduction is not worth it since it comes at an additional cost of limit state function evaluations: from table 5.1 the total calls to g for the MetaIS-SDR is 5 times that of standard MetaIS. For higher dimensions however, it appears that the metamodelling in a reduced space limits the number of limit state function evaluations: this reduction occurs at the Kriging model construction but is even more noticeable for the correction factor estimation. This is probably because a more accurate Kriging model is obtained in the one-dimensional SDR subspace leading to an importance sampling density closer to the optimal one hence achieving lower estimation variance.

Reliability assessment of a wind turbine in stationary conditions

We now turn to an industrial case study where the goal is the reliability analysis of an onshore wind turbine under wind loading. We use a 5MW wind turbine mode developed by the National Renewable Energy Laboratory (NREL). The turbine has a 90m hub height and implements a blade and generator control strategy. The mechanical response analysis to the wind sample n 1 ∈ [10 4 , 5 × 10 4 ] points z i , i = 1, . . . , n 1 from h r (defined in ( 41)) using procedure from lemma 1.

7:

add N add points to DoE according to the procedure outlined in 4.2.

8:

if N D > N max D then break 9:
end if 10: end while Failure probability estimation 1: Augmented probability estimation: compute P f,ǫ , δǫ using standard Monte Carlo or the splitting algorithm in 2.2.4 2: Correction factor estimation: run an MCMC of length N corr = b + N max -N S DR -N D targeting qr * as described in lemma 1.

Compute αcorr and δα = σcorr αcorr according to ( 20) and ( 21) 3: Failure probability estimate:

P f = αcorr P f,ǫ , δ = δ2 α + δ2 ǫ + δ2 α δ2 ǫ
inflow is obtained using the FAST software which provides extreme and fatigue loads for a a wide array of turbines.

The random wind field is described in terms of its (u, v, w) coordinate where u is a vector pointing towards the mean wind flow, (v, w) completes the orthogonal basis and thus mean wind speeds in the v and w directions are zero. For the purpose of our analysis, we neglect the turbulent wind flow components in the v and w directions and consider the wind speed along u as the random input load.

Wind speed model

The wind speed process {X(t), t ≥ 0} was modelled as a stationary Gaussian process according to a spectral expansion model (Shinozuka and Deodatis, 1991)

X(t) = U 10 + n i=1 (u i σ i cos(ω i t) -ūi σ i sin(ω i t)) ( 48
)
u i , ūi are standard independent normal variables, ω i are the frequencies with increment dω i = w i+1w i and σ 2 i = S (ω)dω i where S (ω) is the power spectrum density (p.s.d.) of {X(t), t ≥ 0}. U 10 is the 10-minute mean wind speed. The spectrum was estimated from real measurements on the Danish coastal site Hornsrev. To achieve near-stationarity, we extracted 10-minute wind speed time series corresponding to a turbine rated wind speed of U 10 = 11.5 m/s and turbulence intensity I = 6%. These measurements where broken down into 5 minutes segments. A p.s.d. estimate was then obtained by averaging the periodograms of each segment. The resulting spectrum is shown in figure 1. For simulation purposes, the spectrum is discretized into n = 50 harmonics so that the random vector U = [u 1 , . . . , u n , ū1 , . . . , ūn ] describing X is 100-dimensional.

Time variant reliability analysis

The objective of this reliability analysis problem is to evaluate the 10-minute failure probability of the turbine. Failure is hereafter characterized as the exceedance of a distance threshold s by the tower-top displacement due to fore-aft forces. Letting {Y(t), t ≥ 0} the mechanical response, we seek the following quantity

P f (T ) = P(∃t ∈ [0, T ] | Y(t) ≥ s) ( 49
)
where T = 10 minutes. Equivalently, we may introduce a timevarying limit state function g(t, U) = s -Y(t) where U are the coefficients of the spectral expansion model, which yields

P f (T ) = P(∃t ∈ [0, T ] | g(t, U) ≤ 0).
It is well known that a related quantity is the expected number of downcrossings E(N + (0, T )) (of level 0 by the limit state g) which is such that

P f (T ) ≤ P(g(0, U) ≤ 0) + E(N + (0, T )) (50)
Furthermore, E(N + (0, T )) = ν + (t)dt, where

ν + (t) = lim ∆t→0 P (g(t, U) > 0, g(t + ∆t, U) ≤ 0) ∆t (51)
is the outcrossing rate. This is the actual quantity of interest in the usual time-varying reliability framework. In most random vibration models, the response Y is modelled thanks to a linearized stochastic differential equation. The non linearity and non Gaussian nature of Y is then due to the loading force term (see [START_REF] Jensen | Extreme wave and wind predictions[END_REF]. Assuming wind speed stationarity, one does not automatically obtain output stationarity because of the coupling that goes on between the loading force and the response. However, simulations show that for stationary wind conditions, the response can generally be considered stationary after a transitory period which corresponds to the memory of the system. This means ν + is independent of the time t for sufficiently large t. Therefore, its computation can be obtained by considering the limit state function at a fixed time t 0 . In practice, the simulation length must be sufficiently high so as to discard the initial transient part of the response signal. We however acknowledge that the stationarity assumption of the response assuming stationary inputs must be assessed more carefully but this is left out for a future study.

An approximation of the outcrossing rate may be obtained by considering an integration time-step ∆t << 1 so that

ν+ = P (g(t 0 , U) > 0, g(t 0 + ∆t, U) ≤ 0) ∆t (52)
In our settings, we consider a threshold on the tower-top displacement s = 0.4 m and ∆t = 0.01. The simulation length is one minute.

Instantaneous failure probability

To illustrate the performance of the MetaIS-DR algorithm with respect to other simulation-based reliability algorithms, we first look at the instantaneous failure probability at time t 0 = 60 s. As for the academic examples, we also implement for comparison purposes a design point based method, an importance sampler, and subset simulation [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF]):

• We considered a multiple failure point (MFP) search based on the approach of Der Kiureghian and Dakessian (1998): this method has the ability to find several failure points which contribute significantly to the failure probability. This algorithm essentially modifies the limit state function in the vicinity of a previously found design point thus redefining a new constrained optimization problem geared towards finding a failure point far enough from previous solutions. The constrained optimization method used in our simulations was the sequential quadratic approximation (SQA) algorithm, a gradient free method which proved efficient in terms of number of limit state function evaluations. The MFP search was done with a limit of 5 distinct design points.

• The importance sampler used a mixture proposal where each component was a unit covariance Gaussian distribution centered on a failure points identified by the MFP approach

• The subset simulation algorithm was ran using the original algorithm by [START_REF] Au | Estimation of small probabilities of failure in high dimensions by subset simulation[END_REF]. The conditional probabilities where set to p 0 = 0.1 and a modified Metropolis-Hastings chain of length N subset = 1000 was generated at each subset level, using a thinning interval of length 3.

The Meta-ISDR algorithm was implemented with the following setup:

• a maximum dimension of the SDR subspace r max = 6

• the dimension reduction subspace was estimated with the KDR algorithm using N KDR = 1000 training samples. The training set was constructed by drawing 50% samples from a uniform distribution in the hypersphere in R d of radius R max = √ q χ 2 ,1-10 -8 and 50% samples from a Gaussian mixture distribution centered on the design points. Note that this R max parameter corresponds to a lower bound on the instantaneous failure probability of 10 -8 (see section 4.6).

• a Kriging model with constant mean β and squared exponential covariance function

• initial DoE of size N init D = 10r + M: Latin Hypercube Sampling of size 10 × r + m "design points", where r is the estimated reduced dimension For all methods, the target coefficient of variation was 7%. For the Meta-ISDR algorithm, a maximum number of N max = 11000 limit state function evaluations was set, including the number of calls due to the preliminary multiple failure points.

The multiple failure points algorithm found two significant design points. A FORM approximation based on the lowest reliability index was then performed, while both design points were used in the importance sampler and MetaIS-DR algorithm. The results in table 5.2.3 indicate the good performance of the MetaIS-DR algorithm as it achieves both the target coefficient of variation while also requiring the least amount of

+ f,ǫ (Q) α corr (Q) Multi-FORM 2.50 × 10 -5 N/A 4657 N/A N/A MetaIS-DR
1.99 × 10 -5 0.16 10657 2.07 × 10 -3 9.63 × 10 -3

Table 3: Short-term outcrossing rate estimation limit state function evaluations. The standard MetaIS requires 37% more limit state evaluations. The performance of the subset simulation algorithm in terms of estimation error is due to the limited length of Markov chains simulated at each level, which isn't enough to mitigate the chain's autocorrelation. A better error can be achieved by increasing the thinning interval for instance but this results in running a longer chain and more evaluations of the limit state.

Outcrossing rate estimation

As previously stated, the object of interest for time-variant reliability analysis in stationary conditions is the outcrossing rate ν + . To evaluate this quantity, let us rewrite (52) in a more convenient fashion ν+ = P (g(t 0 , U) > 0, g(t 0 + ∆t, U) ≤ 0) ∆t = P min (-g(t 0 , U), g(t 0 + ∆t, U)) ≤ 0 ∆t -P g(t 0 , U) = 0, g(t 0 + ∆t, U) ≤ 0 = P min (-g(t 0 , U), g(t 0 + ∆t, U)) ≤ 0 ∆t = P H t 0 ,t 0 +∆t (U) ≤ 0 ∆t (53)

where H t 0 ,t 0 +∆t (U) = min (-g(t 0 , U), g(t 0 + ∆t, U)) In (53), we implicitly assume that the response process {Y(t), t ≥ 0} and therefore the limit state g(t 0 , U), admit a density leading to P (g(t 0 , U) = 0, g(t 0 + ∆t, U) ≤ 0) = 0. The formulation (53) lends itself to standard reliability analysis algorithms and is used to compute numerical outcrossing rate estimates. The FORM outcrossing rate was computed using an approximation due to [START_REF] Koo | Design-point excitation for non-linear random vibrations[END_REF] for i.i.d. standard Gaussian inputs u i :

ν + (t 0 ) = 1 2πβ e -1 2 β 2 n i=1 u * 2 i (t 0 ) + ū * 2 i (t 0 ) ω 2 i ( 54
)
where β = u * (t 0 ) is the reliability index and u * (t 0 ) is the design point of the limit state function g ( t 0 , •). Surprisingly, the FORM approximation for this case study is more accurate than the corresponding instantaneous failure probability computed in the previous subsection and is comparable to the estimation of MetaIS-DR. We however expect the FORM outcrossing rate estimation to be inaccurate for more pronounced non-linearities.

Conclusion

A concern of Monte Carlo based methods for structural reliability assessment such as standard importance sampling or subset simultation is their high computational cost as the dimension exceeds a few tenths and/or the failure probability gets below 10 -4 . Recently, the MetaIS algorithm, which combines a Kriging metamodel of the failure surface and an importance sampling procedure, enables more efficient estimation by sampling a so-called quasi-optimal density which acts as a surrogate to the optimal IS density. In high dimensional settings however, the Kriging metamodel construction demands a non-negligible number of limit state function evaluations so as to achieve sufficient accuracy in the vicinity of the failure region. This number influences directly the achievable variance reduction by the IS scheme. Considering that in some reliability problems, the performance function depends on a projection of the input variables on a lower-dimension subspace, we have suggested to build a surrogate to the limit state function in this reduced subspace. The proposed approach leverages recent sufficient dimension reduction techniques to find this subspace. A cross-validation type procedure is suggested in order to infer the dimension of the reduced subspace. The MetaIS algorithm is then cast into this framework and yields particularly efficient MCMC sampling for Gaussian distributed inputs. Its applicability in dimensions up to 100 is demonstrated on a well known academic example which illustrates the impact on constructing a metamodel in the reduced subspace on the efficacy of the probability estimator for a given confidence level. Finally, an industrial case study focused on the extreme response prediction of a wind turbine shows a notable reduction in the computational cost, compared to existing approaches which yield consistent estimators (subset simulation, standard importance sampling and standard MetaIS). where w T = (w T 1 , w T 2 ). By the same argument, if W has density p W (w) = 1 P f ,ǫ π r (w 1 )ϕ r (w 1 )ϕ d-r (w 2 ) then X = B -T a W ∼ qr * . It is clear from (A.3) that to sample W = (W T 1 , W T 2 ) T from p W , suffice it to sample W 2 ∼ N(0 d-r×1 , I d-r ), and W 1 ∼ p W 1 (w 1 ) = π r (w 1 )ϕ r (w 1 ) P f ,ǫ .

Figure 1 :

 1 Figure 1: One-sided power spectral density of wind speed process

• 300 •

 300 minimum/maximum size of the DoE N min D = 60, N max D = number of points added to DoE during refinement N add = 2r

=E

  Appendix A. Quasi-optimal density sampling in Meta-ISDRWe provide the proof of lemma 1. Let B ∈ R d×r a matrix such that B T B = I r and span(B) is a dimension reduction subspace for the regression ofY on X . Let B a = [B, B ⊥ ] ∈ R d×d where the columns of B ⊥ form an orthonormal basis of span(B) ⊥ . Now assume X ∼ N(0, I d ). X has density q(x) = ϕ d (x) = (2π) -d 2 exp(-x 2 2) where ϕ d is the d dimensional multinormal pdf. If x B and x B ⊥ are the orthogonal projections of x on span(B) and span(B) respectively, thenx = x B + x B ⊥ = B(B T x) + B ⊥ (B T ⊥ x) and x 2 = x B 2 + x B ⊥ 2 = B T x 2 + B T ⊥ x 2 . Hence, q(x) = (2π) -d 2 exp(ϕ r (B T x)ϕ d-r (B T ⊥ x)and the quasi-optimal density is equal toqr * (x) = π r (B T x)ϕ r (B T x)ϕ d-r (B T ⊥ x) P f,ǫNow let X ∼ qr * and consider the mapping X → W = B T a X. By a change of variable, for any continuous bounded functionψ from R d to R d π r (B T B -T a w)ϕ r (B T B -T a w)ϕ d-r (B T ⊥ B -T a w) (A.2)= 1 P f,ǫ π r (w 1 )ϕ r (w 1 )ϕ d-r (w 2 ) (A.3)

Table 1 :

 1 . Failure probability estimation: analytical limit state function(47)

	Method	FORM	Monte Carlo IS-FORM	MetaIS	MetaIS-DR
	d = 2
	N P f c.o.v.	18 3.81 × 10 -3 -	5 × 10 5 4.93 × 10 -3 2%	8818 4.90 × 10 -3 4.88 × 10 -3 4.94 × 10 -3 103 518 2% 1.3% 2.1%
	d = 50
	N P f c.o.v.	255 1.54 × 10 -3 -	1.4 × 10 6 1.89 × 10 -3 1.9%	14000 1.95 × 10 -3 1.87 × 10 -3 1.85 × 10 -3 3255 1026 1.9% 2% 1.9%
	d = 100
	N P f c.o.v.	505 3.74 × 10 -5 -	1.5 × 10 6 1.73 × 10 -3 1.8%	18000 1.76 × 10 -3 1.76 × 10 -3 1.71 × 10 -3 5505 1783 1.9% 2% 1.9%

  Algorithm 3 Meta-ISDR algorithm for standard Gaussian inputs SDR Basis estimation:1: Sample x i S DR ∼ N(0, I d ), i = 1, . . . ,N S DR and compute y i S DR = g(x i S DR ) 2: Using gKDR algorithm find a dimension reduction subspace estimate B, r ← gKDR(x i S DR , y i S DR , k X , k Y , K cv , r max )

	Initial DoE:
	1: Add design points x * 1 , . . . , x * m to DoE (optional) 2: Space-filling design: sample N f ill points in the hypersphere of radius R max and add N init D -m corresponding cluster centers x init j to initial DoE by evaluating y init j = g(x init j ), j = 1, . . . , N init D -m 3: Initial metamodel: fit a Kriging model to initial DoE and compute αLOO
	4: Metamodel refinement:
	5: while N D < N min D or αLOO [α min LOO , α max LOO ] do
	6:

Preprint submitted to Structural SafetyDecember 3, 2015