Metamodel-based importance sampling for the simulation of rare events
MMetamodel-based importance sampling forthe simulation of rare events
V. Dubourg , , F. Deheeger , B. Sudret , Clermont Universit´e, IFMA, EA 3867, Laboratoire de M´ecanique et Ing´enieries, BP 10448F-63000 Clermont-Ferrand Phimeca Engineering, Centre d’Affaires du Z´enith, 34 rue de Sarli`eve, F-63800 Cournon d’Auvergne
ABSTRACT: In the field of structural reliability, the Monte-Carlo estimator is considered as the referenceprobability estimator. However, it is still untractable for real engineering cases since it requires a high numberof runs of the model. In order to reduce the number of computer experiments, many other approaches knownas reliability methods have been proposed. A certain approach consists in replacing the original experiment bya surrogate which is much faster to evaluate. Nevertheless, it is often difficult (or even impossible) to quantifythe error made by this substitution. In this paper an alternative approach is developed. It takes advantage ofthe kriging meta-modeling and importance sampling techniques. The proposed alternative estimator is finallyapplied to a finite element based structural reliability analysis.1 INTRODUCTIONReliability analysis consists in the assessment of thelevel of safety of a system. Given a probabilisticmodel (a random vector X with probability densityfunction (PDF) f ) and a performance model (a func-tion g ), it makes use of mathematical techniques inorder to estimate the system’s level of safety in theform of a failure probability. A basic approach, whichmakes reference, is the Monte-Carlo simulation tech-nique that resorts to numerical simulation of theperformance model through the probabilistic model.Failure is usually defined as the event F = { g ( X ) ≤ } ,so that the failure probability is defined as follows: p f = P ( F ) = P ( { g ( X ) ≤ } ) = (cid:90) g ( x ) ≤ f ( x ) d x (1)Introducing the failure indicator function g ≤ beingequal to one if g ( x ) ≤ and zero otherwise, the failureprobability turns out to be the mathematical expecta-tion of this indicator function with respect to the jointprobability density function f of the random vector X . The Monte-Carlo estimator is then derived fromthis convenient definition. It reads: (cid:98) p f MC = (cid:98) E f (cid:2) g ≤ ( X ) (cid:3) = 1 N N (cid:88) k =1 g ≤ (cid:16) x ( k ) (cid:17) (2)where (cid:110) x (1) , . . . , x ( N ) (cid:111) , N ∈ N ∗ is a set of samplesof the random vector X . According to the central limit theorem, this estimator is asymptotically unbi-ased and normally distributed with variance: σ (cid:98) p f MC = V (cid:2)(cid:98) p f MC (cid:3) = p f (cid:0) − p f (cid:1) N (3)In practice, this variance is compared to the unbiasedestimate of the failure probability in order to decidewhether it is accurate enough or not. The coefficientof variation is defined as δ (cid:98) p f MC = σ (cid:98) p f MC / (cid:98) p f MC . Given N , this coefficient dramatically increases as soon asthe failure event is too rare ( p f → ) and proves thatthe Monte-Carlo estimation technique intractable forreal world engineering problems for which the perfor-mance function involves the output of an expensive-to-evaluate black box function – e.g. a finite elementcode. Note that this remark is also true for too fre-quent events ( p f → ) as the coefficient of variation of − (cid:98) p f MC exhibits the same property.In order to reduce the number of simulation runs,a large set of other approaches known as reliabilitymethods have been proposed. They might be classi-fied as follows.A first approach consists in replacing the orig-inal experiment by a surrogate which is muchfaster to evaluate. Among such approaches, thereare the well-known first and second order re-liability methods ( e.g. Ditlevsen and Madsen,1996; Lemaire, 2009), quadratic response surfaces(Bucher and Bourgund, 1990) and the more re-cent meta-models such as support vector machines a r X i v : . [ s t a t . M E ] A p r Hurtado, 2004; Deheeger and Lemaire, 2007), neu-ral networks (Papadrakakis and Lagaros, 2002) and kriging (Kaymaz, 2005; Bichon et al., 2008). Never-theless, it is often difficult or even impossible to quan-tify the error made by such a substitution.The other approaches are the so-called variance re-duction techniques . In essence, these techniques aimsat favoring the Monte-Carlo simulation of the fail-ure event F in order to reduce the estimation vari-ance. These approaches are more robust because theydo not rely on any assumption regarding the func-tional relationship g , though they are still too compu-tationally demanding to be implemented for industrialcases. For an extended review of these techniques,the reader is referred to the book by Rubinstein andKroese (2008).In this paper an hybrid approach is developed. Itis based on both margin meta-models (defined here-after) and the importance sampling technique. It isthen applied to an academic structural reliability prob-lem involving a linear finite-element model and a two-dimensional random field.2 ADAPTIVE PROBABILISTICCLASSIFICATION USING MARGINMETA-MODELS
A meta-model means to a model what the model it-self means to the real-world. Loosely speaking, it is the model of the model . As opposed to the model, itsconstruction does not rely on any physical assump-tion about the phenomenon of interest but rather onstatistical considerations about the coherence of somescattered observations that result from a set of experi-ments. This set is usually referred to as a design of ex-periments (DOE): X = { x , . . . , x m } . It should be care-fully selected in order to retrieve the largest amount ofstatistical information about the underlying functionalrelationship over the input space D x . Here, we attemptto build a meta-model for the failure indicator func-tion g ≤ . In the statistical learning theory (Vapnik,1995) this is referred to as a classification problem.Hereafter, we define a margin meta-model as ameta-model that is able to give a probabilistic predic-tion of the response quantity of interest whose spread( i.e. variance) depends on the lack of informationbrought by the DOE. It is thus reducible by bring-ing more observations into the DOE. In other words,this is an epistemic (reducible) source of uncertainty.To the authors’ knowledge, there exist only two fam-ilies of such margin meta-models: the probabilisticsupport vector machines P -SVM by Platt (1999) andGaussian-Process- (or kriging-) based classification(Santner et al., 2003). The present paper makes use ofthe kriging meta-model, but the overall concept couldeasily be extended to P -SVM. The theoretical aspectsof the kriging prediction are briefly introduced in thefollowing subsection before it is applied to the classi-fication problem of interest. 2.1 Gaussian-process based prediction
The Gaussian-Process based prediction (also knownas kriging ) theory is detailed in the book by Santneret al. (2003). In essence, kriging assumes that the per-formance function g is a sample path of an underlyingGaussian stochastic process G that would read as fol-lows: G ( x ) = f ( x ) T β + Z ( x ) (4)where: • f ( x ) T β denotes the mean of the GP which cor-responds to a classical linear regression modelon a given functional basis { f i , i = 1 , . . . , p } ∈L ( D x , R ) ; • Z ( x ) denotes the stochastic part of the GP whichis modelled as a zero mean, constant variance σ G ,stationary Gaussian process with a given sym-metric positive definite autocorrelation model. Itis fully defined by its autocovariance functionwhich reads ( (cid:0) x , x (cid:48) (cid:1) ∈ D x × D x ): C GG (cid:0) x , x (cid:48) (cid:1) = σ G R (cid:0)(cid:12)(cid:12) x − x (cid:48) (cid:12)(cid:12) , (cid:96) (cid:1) (5)where (cid:96) is a vector of parameters defining R .The most widely used class of autocorrelation func-tions is the anisotropic squared exponential model: R (cid:0)(cid:12)(cid:12) x − x (cid:48) (cid:12)(cid:12) , (cid:96) (cid:1) = exp (cid:32) n (cid:88) k =1 − (cid:12)(cid:12)(cid:12)(cid:12) x k − x (cid:48) k (cid:96) k (cid:12)(cid:12)(cid:12)(cid:12) (cid:33) (6)The best linear unbiased estimation of G at point x is shown (Santner et al., 2003; Severini, 2005, Chap.8) to be the following Gaussian random variate: (cid:98) G ( x ) = G ( x ) |{ g ( x ) , . . . , g ( x m ) }∼ N (cid:16) µ (cid:98) G ( x ) , σ (cid:98) G ( x ) (cid:17) (7)with moments: µ (cid:98) G ( x ) = f ( x ) T (cid:98) β + r ( x ) T R − (cid:16) Y − F (cid:98) β (cid:17) (8) σ (cid:98) G ( x ) = σ G (cid:32) − (cid:20) f ( x ) r ( x ) (cid:21) T (cid:20) T F R (cid:21) − (cid:20) f ( x ) r ( x ) (cid:21)(cid:33) (9)where we have introduced r , R and F such that: r i ( x ) = R ( | x − x i | , (cid:96) ) , i = 1 , . . . , m (10) R ij = R (cid:0)(cid:12)(cid:12) x i − x j (cid:12)(cid:12) , (cid:96) (cid:1) , i = 1 , . . . , m, j = 1 , . . . , m (11) F ij = f i (cid:0) x j (cid:1) , i = 1 , . . . , p, j = 1 , . . . , m (12)At this stage it can easily be proven that µ (cid:98) G ( x i ) = g ( x i ) and σ (cid:98) G ( x i ) = 0 for i = 1 , . . . , m , thus meaning thekriging surrogate is an exact interpolator .iven a choice for the regression and correlationmodels, the optimal set of parameters β ∗ , (cid:96) ∗ and σ ∗ G can then be inferred using the maximum likeli-hood principle applied to the unique sparse obser-vation of the GP sample path grouped in the vector y = (cid:104) g ( x i ) , i = 1 , . . . , m (cid:105) . This inference problem turnsinto an optimization problem that can be solved an-alytically for both β ∗ and σ ∗ G assuming (cid:96) ∗ is known.Thus, the problem is solved in two steps: the max-imum likelihood estimation of (cid:96) ∗ is first solved by aglobal optimization algorithm which in turns gives theoptimal values of β ∗ and σ ∗ G .2.2 Probabilistic classification function
A surrogate-based reliability analysis simply consistsin replacing the performance function g by its meta-model µ (cid:98) G . Again, this meta-model may be a first-or second-order Taylor expansion of the limit-statesurface g = 0 at a so-called design point (Ditlevsenand Madsen, 1996, FORM/SORM), a polynomial re-sponse surface (Bucher and Bourgund, 1990), a neu-ral networks based prediction (Papadrakakis and La-garos, 2002), or a kriging based prediction (Bichonet al., 2008). This surrogate may not be fully accurateand it is difficult or even impossible to quantify theerror made by substitution on the final failure proba-bility of interest.In this paper as in the work by Picheny (2009),it is proposed to use the complete probabilistic pre-diction provided by the kriging meta-model insteadof the sole mean prediction ( e.g. as in Bichon et al.,2008). Indeed, since the probabilistic distribution ofthe prediction is fully characterized, the probabilitythat the prediction is negative may be expressed inclosed-form and reads as follows: P (cid:104) (cid:98) G ( x ) ≤ (cid:105) = Φ (cid:32) − µ (cid:98) G ( x ) σ (cid:98) G ( x ) (cid:33) (13)Note that this latter quantity is not the sought failureprobability p f , this is simply the probability that theprediction (cid:98) G at some deterministic vector x is nega-tive.Picheny (2009) proposes then to use this proba-bilistic classification function as a surrogate for thereal failure indicator function, and uses crude Monte-Carlo simulation to estimate the failure probability. Adifferent use is proposed in the next section.Figure 1 illustrates the concepts introduced in thissection on a basic structural reliability example fromDer Kiureghian and Dakessian (1998). This exam-ple involves two independent standard Gaussian ran-dom variates X and X , and the performance functionreads: g ( x , x ) = b − x − κ ( x − e ) (14)where b = 5 , κ = 0 . and e = 0 . . In subfigure 1(a),the limit-state function g ( x ) = 0 is represented by the black dash-dot line. The red minusses ( g ≤ ) andblue plusses ( g > ) represent the initial DOE fromwhich the kriging meta-model is built. The meanprediction’s limit-state µ (cid:98) G ( x ) = 0 is represented bythe dashed black line. It can be seen that the meta-model is not fully accurate since the green triangle x (among others) is misclassified. Indeed, x is safe ac-cording to the real performance function g , but it failsaccording to the mean prediction of the meta-model µ (cid:98) G . The probabilistic classification function makes asmoother decision possible: x fails with a 60% prob-ability w.r.t. the epistemic uncertainty in the randomprediction (cid:98) G ( x ) ∼ N ( µ (cid:98) G ( x ) , σ (cid:98) G ( x )) . Note also thatthe red and blue points in the DOE fails with prob-abilities 100% and 0% (safe) respectively due to theinterpolating property of the kriging metamodel. Sub-figure 1(a) is the one-dimensional illustration of thethree classification strategies for the vector x . Thedeterministic decision function is an heaviside func-tion centered in zero, and the probabilistic classifica-tion is a smoother Gaussian cumulative density func-tion.2.3 Refinement of the probabilistic classificationfunction
In this subsection, a strategy is proposed in order torefine the probabilistic classification function so thatit tends towards the real indicator function g ≤ .First, let the margin of uncertainty M be defined asfollows: M = (cid:110) x : − k σ (cid:98) G ( x ) ≤ (cid:98) G ( x ) ≤ + k σ (cid:98) G ( x ) (cid:111) (15)where k might be chosen as k = Φ − (97 . . meaning a 95% confidence interval onto the predic-tion of the limit-state surface is chosen. Such a 95%confidence margin is illustrated in Subfigure 1(a) asthe area bounded below by the blue line (2.5% confi-dence level) and above by the red line (97.5% confi-dence level). The points that are located in this marginhave an uncertain sign, the others being either failedor safe with a confidence level greater than 97.5%.We also define the probability that a point x ∈ D x belongs to this margin of uncertainty. Due to theGaussian nature of the prediction, this probabilitymay also be expressed in closed-form and reads asfollows: P [ x ∈ M ] = Φ (cid:32) k σ (cid:98) G ( x ) − µ (cid:98) G ( x ) σ (cid:98) G ( x ) (cid:33) − Φ (cid:32) − k σ (cid:98) G ( x ) − µ (cid:98) G ( x ) σ (cid:98) G ( x ) (cid:33) (16)Then, finding the point that maximizes this quan-tity on the support of the PDF of X will finally bringthe best improvement point in the DOE. Starting withthis statement, many authors in the kriging literature x x x g ( x ) = [ c G ( x ) ] = . % µ c G ( x ) = [ c G ( x ) ] = . % [ c G ( x ) ] = . % [ c G ( x ) ] (a) Limit-states and probability contours. − . σ c G ( x ) g ( x ) µ c G ( x ) + . σ c G ( x ) g ( x ) =0 . c G ( x ) 0] =0 . µ c G ( x ) =1 .
00 0 . .
025 0
DeterministicProbabilistic (b) Classification given x = x .Figure 1: Comparison of the three classification strategies on a two-dimensional example from Der Kiureghian and Dakessian (1998). decide to use global optimization algorithms in or-der to find the best improvement point. For instance,Bichon et al. (2008) use a different criterion namedthe expected feasibility function , and Lee and Jung(2008) use the constraint boundary sampling crite-rion. Note also that an equivalent concept is used byHurtado (2004); Deheeger and Lemaire (2007); De-heeger (2008); Bourinet et al. (2010) for SVM.In this paper, as in Dubourg et al. (2011), a slightlydifferent strategy is proposed in order to add severalpoints in the DOE. The proposed criterion P [ x ∈ M ] ismultiplied by a weighting density function w so that C ( x ) = P ( x ∈ M ) w ( x ) (17)can itself be regarded as a PDF up to an unknown butfinite normalizing constant. The weighting density w can either be chosen as the original PDF of X , or, asit is proposed here, the uniform PDF on a sufficientlylarge confidence region of the original PDF. Such aconfidence region might be difficult to define for anygiven PDF, but as it is usually done in structural relia-bility (Ditlevsen and Madsen, 1996), the original ran-dom vector X can be transformed into a probabilisti-cally equivalent standard Gaussian random vector U for which the confidence region is simply an hyper-sphere with radius β . The reader is referred to Lebrunand Dutfoy (2009) for a recent discussion on suchmappings U = T ( X ) . In that given space, β can beeasily selected as e.g. β = 8 which corresponds to themaximal generalized reliability index (Ditlevsen andMadsen, 1996) that can be justified numerically, andthe sought uniform PDF is simply defined in terms ofthe following indicator function: w ( u ) ∝ √ u T u ≤ β ( u ) (18)Markov-chain Monte-Carlo simulation techniques( e.g. the slice sampling technique proposed by Neal, 2003) might be used in order to generate N (say N = 10 ) samples from the pseudo-PDF C . These sam-ples are expected to be highly concentrated aroundthe maxima of the criterion P ( x ∈ M ) , and thus in thevicinity of the predicted limit-state µ (cid:98) G ( x ) = 0 wherethe sign of (cid:98) G is the most uncertain. This large candi-date population can then be reduced to a smaller onethat condensate its statistical properties by means ofa K -means clustering algorithm (MacQueen, 1967).The K ( K being given) cluster centers uniformly spanthe margin M and may be added to the DOE in orderto enrich the prediction of the performance functionin the vicinity of the limit-state and thus reduce themargin of uncertainty.3 META-MODEL-BASED IMPORTANCESAMPLINGPicheny (2009) proposes to use the probabilistic clas-sification function as a surrogate for the real indica-tor function, so that the failure probability is rewrittenfrom its definition in Eq. (1) as follows: p f ε = (cid:90) P (cid:104) (cid:98) G ( x ) ≤ (cid:105) f ( x ) d x ≡ E f (cid:104) P (cid:104) (cid:98) G ( X ) ≤ (cid:105)(cid:105) (19)It is argued here that this latter quantity does notequal the failure probability of interest because itsums the aleatory uncertainty in the random vector X and the epistemic uncertainty in the prediction (cid:98) G . Thisis the reason why p f ε will hereafter be referred to asthe augmented failure probability . As a matter of fact,even if the epistemic uncertainty in the prediction canbe reduced ( e.g. by enriching the DOE as proposed insection 2.3), it is impossible to quantify the contribu-tion of each source of uncertainty a posteriori .his remark motivates the approach introducedin this section where the probabilistic classificationfunction is used in conjunction with the importancesampling technique in order to build a new estimatorof the failure probability.3.1 Importance sampling
According to Rubinstein and Kroese (2008), impor-tance sampling (IS) is the most efficient variance re-duction technique. This technique consists in comput-ing the mathematical expectation of the failure indica-tor function according to a biased PDF which favorsthe failure event of interest. This PDF is called the instrumental density .Given an instrumental density h , such that h dom-inates g ≤ f , the definition of the failure probabilityof Eq. (1) may be rewritten as follows: p f = (cid:90) g ( x ) ≤ f ( x ) h ( x ) h ( x ) d x ≡ E h (cid:20) g ≤ ( X ) f ( X ) h ( X ) (cid:21) (20)where the subscript h on the expectation operator isadded to recall that X is therefore distributed accord-ing to h . Note that the domination requirement of h over g ≤ f simply means that: ∀ x ∈ D x , h ( x ) = 0 ⇒ g ≤ ( x ) f ( x ) = 0 (21)so that the so-called likelihood ratio (cid:96) ( x ) = f ( x ) /h ( x ) is finite for any given x ∈ D x .The latter definition of the failure probability easilyleads to the establishment of the importance samplingestimator which reads as follows: (cid:98) p f IS = (cid:98) E f (cid:2) g ≤ ( X ) (cid:3) = 1 N N (cid:88) k =1 g ≤ (cid:16) x ( k ) (cid:17) f ( x ( k ) ) h ( x ( k ) ) (22)where (cid:110) x (1) , . . . , x ( N ) (cid:111) , N ∈ N ∗ is a set of samples from h . According to the central limit theorem, this estima-tion is unbiased and its quality may be measured bymeans of its variance of estimation which reads: Var (cid:2)(cid:98) p f IS (cid:3) =1 N − N N (cid:88) k =1 g ≤ ( x ( k ) ) (cid:32) f ( x ( k ) ) h ( x ( k ) ) (cid:33) − (cid:98) p f IS (23)Rubinstein and Kroese (2008) show that this vari-ance is zero (optimality of the IS estimator) when theinstrumental PDF is chosen as: h ∗ ( x ) = g ≤ ( x ) f ( x ) (cid:82) g ≤ ( x ) f ( x ) d x = g ≤ ( x ) f ( x ) p f (24)However this instrumental PDF is not implementablein practice because it involves the sought failure prob-ability in its denominator. There exists infinitely manyPDF h that allows to significantly reduce the varianceof estimation though. 3.2 A meta-model-based approximation of theoptimal instrumental PDF
Different strategies have been proposed in order tobuild quasi-optimal instrumental PDF suited for spe-cific estimation problems. For instance, Melchers(1989) uses a standard normal PDF centered onto the most probable failure point (MPFP) in the space ofthe independent standard Gaussian random variables U = T ( X ) in order to estimate a failure probability.Although this approach may lose accuracy as soon asthe MPFP is not unique. Cannamela et al. (2008) usea kriging prediction of the performance function g inorder to build an instrumental PDF suited for the es-timation of extreme quantiles of the random variate G = g ( X ) .Here, it is proposed to use the probabilistic classi-fication function in Eq. (13) as a surrogate for the realindicator function in the optimal instrumental PDF inEq. (24). The proposed quasi-optimal PDF thus readsas follows: (cid:99) h ∗ ( x ) = P (cid:104) (cid:98) G ( x ) ≤ (cid:105) f ( x ) (cid:82) P (cid:104) (cid:98) G ( x ) ≤ (cid:105) f ( x ) d x ≡ (cid:98) G ≤ ( x ) f ( x ) p f ε (25)where p f ε is the augmented failure probability whichhas been already defined in Eq. (19). This quasi-optimal instrumental PDF is compared to its impracti-cal optimal counterpart in Figure 2 using the exampleof Section 2.2.3.3 The meta-model-based importance samplingestimator
Choosing the proposed quasi-optimal instrumentalPDF in Eq. (25) in the importance sampling defini-tion of the failure probability in Eq. (20) leads to thefollowing new definition: p f = (cid:90) g ≤ ( x ) f ( x ) (cid:99) h ∗ ( x ) (cid:99) h ∗ ( x ) d x (26) = p f ε (cid:90) g ≤ ( x ) P (cid:104) (cid:98) G ( x ) ≤ (cid:105) (cid:99) h ∗ ( x ) d x (27) ≡ p f ε α corr (28)where we have introduced: α corr ≡ E (cid:99) h ∗ g ≤ ( X ) P (cid:104) (cid:98) G ( X ) ≤ (cid:105) (29)This means that the failure probability is now definedas the product between the augmented failure proba-bility p f ε and a correction factor α corr . This correctionfactor is defined as the expected ratio between the real x h ∗ ( x ) g ( x ) = f ( x ) (a) The optimal instrumental PDF. x x d h ∗ ( x ) g ( x ) = f ( x ) (b) A quasi-optimal PDF.Figure 2: Comparison of the instrumental PDF on the two-dimensional example from Der Kiureghian and Dakessian (1998). indicator function g ≤ and the probabilistic classifi-cation function P [ (cid:98) G ( • ) ≤ . Thus, if the kriging pre-diction is fully accurate, the correction factor is equalto unity and the failure probability is identical to theaugmented failure probability (optimality of the pro-posed estimator). On the other hand, in the more gen-eral case where the kriging prediction is not fully ac-curate, the correction factor modifies the augmentedfailure probability accounting for the epistemic uncer-tainty in the prediction.The two terms of the latter definition of the failureprobability may now be estimated using Monte-Carlosimulation: (cid:98) p f ε = 1 N ε N ε (cid:88) k =1 P (cid:104) (cid:98) G ( x ( k ) ) ≤ (cid:105) (30) (cid:98) α corr = 1 N corr N corr (cid:88) k =1 g ≤ ( x ( k ) ) P (cid:104) (cid:98) G ( x ( k ) ) ≤ (cid:105) (31)where the first N ε -sample is generated from the orig-inal PDF f , and the second N corr -sample is generatedfrom the quasi-optimal instrumental PDF (cid:99) h ∗ . Accord-ing to the central limit theorem, these two estimatesare unbiased and normally distributed. Their respec-tive variance of estimation denoted by σ ε and σ corr arenot given here but they might be easily derived.To generate samples from (cid:99) h ∗ , it is proposed touse a Markov chain Monte-Carlo simulation tech-nique which is applicable to a broad class of improperPDF for which the normalizing constant is not knownand thus for the instrumental PDF of interest (cid:99) h ∗ ( x ) ∝ P [ (cid:98) G ( x ) ≤ f ( x ) . The work presented here makes useof the slice sampling technique (Neal, 2003).Finally, the final estimator of the failure probabilitysimply reads as follows: (cid:98) p f metaIS = (cid:98) p f ε (cid:98) α corr (32) The calculation of the coefficient of variation of the fi-nal estimator (cid:98) p f metaIS will be detailed in a forthcomingpaper.4 RELIABILITY ANALYSIS OF AN 8-HOLEPLATEThis structural reliability example is inspired fromDeheeger and Lemaire (2007). It concerns the relia-bility analysis of a × mm 8-hole plate illus-trated in Figure 3. The diameter of the holes ∅ is setequal to 10 mm. Its left end is clamped both horizon-tally and vertically while its right end is subjected toa distributed line load with magnitude q = 100 MPa.Plain stress is assumed and the material is supposedto have a linear elastic behavior. The Poisson coeffi-cient ν is set equal to 0.3. Due to the boundary con-ditions the Poisson effect is not the same on all theplate though. The Young’s modulus is modeled by anhomogeneous lognormal random field with a mean µ E = 200 000 MPa, a coefficient of variation δ E = 25% and assuming an isotropic squared exponential auto-correlation function with a 20 mm correlation length (cid:96) . The two-dimensional random field is represented bya translated Karhunen-Loeve expansion discretizedby means of a wavelet-Galerkin strategy proposed byPhoon et al. (2002). The stochastic model involves20 independent standard Gaussian random variatesgrouped in the vector X to simulate the random field.The mechanical model is solved with Code Aster(eDF, R&D Division, 2006) in order to retrieve themaximal Von Mises stress in the plate P . The perfor-mance function is then defined as follows: g ( x ) = σ − max p ∈P { σ Von Mises ( p ) } (33)with respect to an arbitrary threshold σ = 450 MPa. l a m p e d q = M P a (a) Mesh, boundary conditions and loads. C l a m p e d q = M P a (b) One sample path of the Young modulus random field.Figure 3: Illustration of the 8-hole plate example from Deheeger and Lemaire (2007). DOE MPFP Simulations P f estimate C. o. V.Subset (ref.) - - 25 000 1.70 × − multi-FORM - 1 168 - 0.65 × − - meta-IS × − ¡ 10%Table 1: Reliability analyses results for the 8-hole plate example from Deheeger and Lemaire (2007). The proposed meta-model-based importance sam-pling procedure is applied to this structural reliabilityexample. First an initial kriging predictor is built forthe performance function g using a 100-point DOE.These 100 points are uniformally generated withinthe β radius hypersphere. Based on this initial pre-diction, the DOE refinement procedure introduced inSection 2.3 is used. K = 100 new points are added ateach refinement iteration. The refinement procedureis stopped after 1 000 estimations of the performancefunction. This may seem arbitrary but it is difficult toprovide another stopping criterion for the refinementprocedure – this needs further investigation. Then, theprobabilistic classification function is defined with re-spect to the latest (finest) kriging prediction and it isused to compute the proposed estimator of the failureprobability.The results are provided in Table 1. They are com-pared to a reference solution obtained by subset sim-ulation (Au and Beck, 2001), and the multi-FORMestimator from Der Kiureghian and Dakessian (1998)using FERUM v4.0 (Bourinet et al., 2009) implemen-tations of these algorithms. FERUM is a Matlab tool-box for reliability analysis published under the Gen-eral Public License. The estimate of the augmentedfailure probability is equal to (cid:98) p f ε = 2 . × − , and thecorrection factor is equal to (cid:98) α corr = 0 . . It means thatthe kriging predictor is rather accurate in that case.The probabilistic classification function is very closeto its deterministic counterpart – and so is the instru-mental importance sampling density (cid:98) h .5 CONCLUSIONStarting from the double premise that a surrogate-based reliability analyses does not permit to quantifythe substitution error, and that the existing variance reduction techniques remain time-consuming whenthe performance function involves the output of anexpensive-to-evaluate black box function, an hybridstrategy has been proposed. First, the probabilisticclassification function was introduced, this functionallows a smoother classification than its determinis-tic counterpart accounting for the epistemic uncer-tainty in the kriging prediction. Using this smootherclassification function within an importance samplingframework then allowed to derive a meta-model-based importance sampling estimator. This estimatorconverges towards the theoretically impractical opti-mal importance sampling estimator and may providea significant reduction of the estimation variance asillustrated in the example.In the present paper, the refinement procedure thatleads to the probabilistic classification function isstopped arbitrarily. Work is in progress in order to es-tablish the best trade-off between the size of the DOEand the number of simulations required to estimatethe correction factor α corr .ACKNOWLEDGEMENTSThe first author is funded by a CIFRE grant fromPhimeca Engineering S.A. subsidized by the ANRT(convention number 706/2008). The financial supportfrom the ANR through the KidPocket project is alsogratefully acknowledged.REFERENCESAu, S. & J. Beck (2001). Estimation of small failureprobabilities in high dimensions by subset simula-tion. Prob. Eng. Mech. 16 (4), 263–277.Bichon, B., M. Eldred, L. Swiler, S. Mahadevan, &. McFarland (2008). Efficient global reliabilityanalysis for nonlinear implicit performance func-tions.
AIAA Journal 46 (10), 2459–2468.Bourinet, J.-M., F. Deheeger, & M. Lemaire (2010).Assessing small failure probabilities by combinedsubset simulation and support vector machines.
Submitted to Structural Safety .Bourinet, J.-M., C. Mattrand, & V. Dubourg (2009). Areview of recent features and improvements addedto FERUM software. In
Proc. ICOSSAR’09, IntConf. on Structural Safety And Reliability, Osaka,Japan .Bucher, C. & U. Bourgund (1990). A fast and efficientresponse surface approach for structural reliabilityproblems.
Structural Safety 7 (1), 57–66.Cannamela, C., J. Garnier, & B. Iooss (2008). Con-trolled stratification for quantile estimation.
Annalsof Applied Statistics 2 (4), 1554–1580.Deheeger, F. (2008).
Couplage m´ecano-fiabiliste, SMART m´ethodologie d’apprentissage stochas-tique en fiabilit´e . Ph. D. thesis, Universit´e BlaisePascal - Clermont II.Deheeger, F. & M. Lemaire (2007). Support vectormachine for efficient subset simulations: 2SMARTmethod. In
Proc. 10th Int. Conf. on Applicationsof Stat. and Prob. in Civil Engineering (ICASP10),Tokyo, Japan .Der Kiureghian, A. & T. Dakessian (1998). Multipledesign points in first and second-order reliability.
Structural Safety 20 (1), 37–49.Ditlevsen, O. & H. Madsen (1996).
Structural reli-ability methods (Internet (v2.3.7, June-Sept 2007)ed.). John Wiley & Sons Ltd, Chichester.Dubourg, V., B. Sudret, & J.-M. Bourinet (2011).Reliability-based design optimization using krig-ing and subset simulation.
Struct. Multidisc. Op-tim. Accepted .eDF, R&D Division (2006).
Code Aster : Analyse desstructures et thermo-m´ecanique pour des ´etudes etdes recherches, V.7 . .Hurtado, J. (2004). Structural reliability – Statisticallearning perspectives , Volume 17 of
Lecture notesin applied and computational mechanics . Springer.Kaymaz, I. (2005). Application of kriging methodto structural reliability problems.
StructuralSafety 27 (2), 133–151.Lebrun, R. & A. Dutfoy (2009). An innovating anal-ysis of the Nataf transformation from the copulaviewpoint.
Prob. Eng. Mech. 24 (3), 312–320. Lee, T. & J. Jung (2008). A sampling techniqueenhancing accuracy and efficiency of metamodel-based RBDO: Constraint boundary sampling.
Computers & Structures 86 (13-14), 1463–1476.Lemaire, M. (2009).
Structural Reliability . John Wi-ley & Sons Inc.MacQueen, J. (1967). Some methods for classifi-cation and analysis of multivariate observations.In J. Le Cam, L.M. & Neyman (Ed.),
Proc. 5 th Berkeley Symp. on Math. Stat. & Prob. , Volume 1,Berkeley, CA, pp. 281–297. University of Califor-nia Press.Melchers, R. (1989). Importance sampling in struc-tural systems.
Structural Safety 6 (1), 3–10.Neal, R. (2003). Slice sampling.
Annals Stat. 31 ,705–767.Papadrakakis, M. & N. Lagaros (2002). Reliability-based structural optimization using neural net-works and Monte Carlo simulation.
Comput. Meth-ods Appl. Mech. Engrg. 191 (32), 3491–3507.Phoon, K., S. Huang, & S. Quek (2002). Simulationof second-order processes using Karhunen-Lo`eveexpansion.
Computers & Structures 80 (12), 1049–1060.Picheny, V. (2009).
Improving accuracy and compen-sating for uncertainty in surrogate modeling . Ph.D. thesis, University of Florida.Platt, J. (1999). Probabilistic outputs for support vec-tor machines and comparisons to regularized like-lihood methods. In
Advances in large margin clas-sifiers , pp. 61–74. MIT Press.Rubinstein, R. & D. Kroese (2008).
Simulation andthe Monte Carlo method . Wiley Series in Probabil-ity and Statistics. Wiley.Santner, T., B. Williams, & W. Notz (2003).
Thedesign and analysis of computer experiments .Springer series in Statistics. Springer.Severini, T. (2005).
Elements of distribution the-ory . Cambridge series in Statistical and Probabilis-tic mathematics. Cambridge University Press.Vapnik, V. (1995).