An efficient reliability method combining adaptive importance sampling and Kriging metamodel

Accepted Manuscript An efficient reliability method combining Adaptive Importance Sampling and Kriging metamodel Hailong Zhao, Zhufeng Yue, Yongshou Liu, Zongzhan Gao, Yishang Zhang PII: DOI: Reference:

S0307-904X(14)00488-0 http://dx.doi.org/10.1016/j.apm.2014.10.015 APM 10168

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

21 February 2013 12 May 2014 2 October 2014

Please cite this article as: H. Zhao, Z. Yue, Y. Liu, Z. Gao, Y. Zhang, An efficient reliability method combining Adaptive Importance Sampling and Kriging metamodel, Appl. Math. Modelling (2014), doi: http://dx.doi.org/ 10.1016/j.apm.2014.10.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An efficient reliability method combining Adaptive Importance Sampling and Kriging metamodel

Hailong Zhao, Zhufeng Yue*, Yongshou Liu, Zongzhan Gao, Yishang Zhang Institute of Aircraft Reliability Engineering, School of Mechanics, Civil Engineering & Architecture, Northwestern Polytechnical University, Xi’an 710129, China

* Corresponding author: [email protected], Tel/Fax: +86-29-88431002

ABSTRACT In practice, there are many engineering problems characterized by complex implicit performance functions. Accurate reliability assessment for these problems usually requires very time-consuming computation, and sometimes it is unacceptable. In order to reduce the computational load, this paper proposes an efficient reliability method combining adaptive importance sampling and Kriging model based on the active learning mechanism. It inherits the superiorities of Kriging metamodel, adaptive importance sampling and active learning mechanism, and enables only evaluating the interested samples in actual performance function. The proposed method avoids a large number of time-consuming evaluation processes, and the important samples are mainly predicted by a well-constructed Kriging metamodel, thus the calculating efficiency is increased significantly. Several examples are given as validations, and results show that the proposed method has great advantages in the aspect of both efficiency and accuracy.

Keywords: Reliability; Importance Sampling; Kriging metamodel; Active learning; Markov chain

1. Introduction In actual engineering problems, there always exist various uncertainties (such as material properties, structure dimensions, boundary conditions), and

these

uncertainties eventually affect the structural safety. Reliability analysis is the technique which provides the estimation of the safety degree by taking these uncertainties into consideration. Usually, analytical explicit performance function cannot be obtained easily, under which the external numerical tools (such as finite element method) have to be adopted to simulate the structural behavior and evaluate the performance function. Monte Carlo Simulation (MCS) [1] is the most commonly used method when dealing with these kinds of implicit reliability performance function. However, an accurate estimation of the failure probability requires a large number of simulations. As a consequence, the computational demand is excessively large and sometimes unacceptable for dealing with small failure probability in engineering problems. To solve the computational efficiency problems, researchers have proposed a lot of improved methods. Importance Sampling (IS) [2-5] is the most commonly method for its relatively high sampling efficiency and convenient execution, of which the basic principle is sampling with importance sampling density function to get a higher efficiency. IS increases the simulation efficiency to a certain extent, but the number of calls to actual performance function is still excessive. When dealing with complex external tools, evaluating the reliability of engineering models requires extremely time-consuming processes despite the improvements in computer technology. As another way to solve the efficiency problems, surrogate models are introduced in reliability field by some researchers. For these kinds of methods, a surrogate model is firstly built to approximate the actual performance function, then by analyzing this metamodel, the reliability can be estimated. Given that the time demand for evaluating the analytical metamodel is much less than that of actual performance function, the computational efficiency can be evidently increased. Recently, response surface [6-9], neural network [10-14], support vector machine [14-17] and Kriging metamodel [18-22] are getting much attention, and they all can be classified as

surrogate models. As a well-known spatial interpolation technique, Kriging has been recently gained attentions in literature [22-25]. It is a semi-parametric interpolation technique, which gives an unbiased estimation with minimum variance. Compared with traditional response surface methods, Kriging does not need to construct a particular mathematical model, and its application is more flexible. Furthermore, Kriging presents several interesting characters [20, 26]: first, it is an accurate interpolation method; the prediction at a point belonging to the design of experiment is exactly accurate. Second, with the existence of stochastic process, Kriging provides not only predicted value but also the so-called Kriging variance, and this variance can illustrate the local uncertainty of the predicted value. The higher the Kriging variance is, the less certain the prediction will be, and vice versa. Due to the superiority listed above, Kriging metamodel has been used in structural reliability. Giunta et al. [27] compared the prediction accuracy of Kriging with classical polynomial response surface, and pointed out that Kriging shows a better fitting accuracy in the aspect of global approximation. Kaymaz [18] proposed a Kriging-based reliability method on foundation of a certain experimental design, and the failure probability is finally obtained by estimating the Kriging metamodel. Bichon et al. [19] developed a method to adaptively search additional samples, and these samples are subsequently used to update the design of experiment and make Kriging more precise. Echard et al. [20] proposed a reliability method called AK-MCS combining Kriging with MCS. With an elaborate active learning function, Kriging is updated in an iterative way and it performs better approximating ability. Several examples have demonstrated the accuracy and efficiency of AK-MCS. Inspired by AK-MCS, this paper proposes an efficient reliability method combining adaptive importance sampling and Kriging. The proposed method aims at estimating the structural safety in an efficient and accurate way. Its fundamental purpose is to carry out an importance sampling method without evaluating the actual performance function. The sign of each point is predicted by a Kriging metamodel constructed with a few evaluated points. The proposed method enables analysts to

focus only on the points locating in the interested region and having significant impact on the accuracy of Kriging metamodel. Then, the approximation of the limit state is accurate in the important region. This method gives at the same time the estimation of failure probability and its coefficient of variation without expensive evaluation of the whole important population. The proposed method alleviates the computational burden of classic numerical simulation method. When dealing with engineering problems with expensive external commercial tools, the method performs excellent in the aspect of efficiency and accuracy. This paper includes the following parts. Section 2 reviews the fundamental theory of importance sampling, Markov chain Metropolis algorithm and Kriging metamodel. Section 3 describes the detailed principles and steps of the proposed method. Section 4 proposes several academic examples to validate the proposed method. The results are also compared with some other available methods to illustrate the superiorities. Finally this paper is concluded in Section 5.

2. Fundamental theory 2.1. Importance Sampling method Generally, failure probability Pf can be expressed as follows, Pf = P ( g ( x ) ≤ 0 ) = ∫ n I F ( x ) f X ( x ) dx

(1)

R

where x is the vector, and denotes basic input variables, g ( x ) is the performance function, R n shows the n -dimension variable space, I F ( ⋅) is the denoting function,

I F ( x ) = 1 when

g (x) ≤ 0

and

IF (x) = 0

when

g ( x) > 0

as an

opposite, and f X ( x ) is the joint probability density function (PDF) of basic variables. The basic principle of importance sampling is to shift the sampling center towards the failure region to gain information more efficiently. By introducing importance the sampling density function hX ( x ) , Eq.(1) can be expressed in the following form,

Pf = ∫ n I F ( x ) R

fX ( x) hX ( x ) dx hX ( x )

(2)

Obviously, the success of importance sampling relies on the proper choice of importance sampling density (ISD). The optimal ISD can be drawn through theoretical analysis [2-4], hopt ( x ) = f ( x | F ) = I F ( x ) f X ( x ) Pf

(3)

In practice, the optimal ISD is not feasible since it requires knowledge of failure probability as a priori. One popular choice for ISD is moving sampling center to the design point (also the most probable failure point). However, in order to obtain the design point, a constrained optimization problem should be solved firstly, which may occupy a considerable portion of the total computational effort. As another good choice for forming ISD, the adaptive importance sampling is introduced by Au et al. [3]. They made use of Markov chain Metropolis algorithm to adaptively generate a scheme of important samples and avoided the solution of design point(s). Yuan et al. [5] modified Au’s methodology and proposed another adaptive importance sampling method based on Markov chain and low-discrepancy sequence. In Yuan’s method, the Markov chain states were selected as importance sampling centers. Both Markov chain algorithm and importance sampling were carried out at the same time and the important region was adaptively populated.

2.2. Markov chain Metropolis algorithm Markov chain Metropolis algorithm is considered as a powerful tool to generate samples according to arbitrary PDF. In this paper, it is used to efficiently generate samples in the failure region. In Markov chain Metropolis algorithm, samples are simulated as the states of a special Markov chain whose stationary distribution is equal to the target PDF. As the number of Markov steps increases, the PDF of Markov states converges to the target PDF. The procedure for generating samples in the failure region is briefly described below.

(1) Determine the stationary and proposal distribution

The stationary distribution, as mentioned above, is the PDF Markov chain states asymptotically distribute when the steps increase. In order to generate samples in the failure region, the stationary distribution is chosen as the optimal ISD hopt ( x ) = I F ( x ) f X ( x ) Pf . The proposal distribution is the rule determining how to

generate new point from the present state. Experience shows that its type is not crucial in the determination of efficiency, and thus the distributions which can be operated conveniently could be employed. In this paper, the proposal distribution is formed by centering the original joint probability density function to present state as classic importance sampling method does.

(2) Select an initial state In order to start the Markov chain sequence, an initial state in failure region is needed. In most structural reliability problems, the point can be assigned based on engineering judgment [3, 5], and this does not present serious difficulties.

(3) Determine the j -th step state of Markov chain The j -th step state x j is generated on the foundation of x j−1 . A candidate state ε is firstly generated from proposal distribution depending on x j −1 . Then, the ratio

r=

hopt ( ε ) hopt ( x j −1 )

=

I F ( ε ) f X ( ε ) Pf I F ( x j −1 ) f X ( x j −1 ) Pf

=

IF (ε ) f X (ε ) I F ( x j −1 ) f X ( x j −1 )

(4)

is computed. According to Metropolis-Hasting algorithm, if r ≥ 1 , ε is accepted as the j -th step state, that is, x j = ε , whereas if r < 1 , ε is accepted as x j with

probability r and x j−1 is accepted as x j with the remaining probability 1 − r , where a repeated state is obtained. Note that in this step, whether the candidate state ε is accepted or not, it has to be evaluated in actual performance function to compute the denoting function I F ( ε ) . In this paper, all the true evaluations will be used to construct the Kriging metamodel and no resource is wasted. (4) Repeat step (3) until enough Markov chain states are obtained

With the introduction of Markov chain Metropolis algorithm, the important region will be adaptively populated, which makes the proposed method more applicable for various reliability problems. 2.3. Kriging metamodel Generally speaking, Kriging metamodel consists of two parts: linear regression and stochastic process [28]. Assume that x denotes input variables and y ( x ) denotes response, Kriging is given as: y ( x ) = F ( β, x ) + z ( x )

(5)

where F ( β, x ) is the deterministic part and represents an averaged approximation of the response. It can be expressed as an ordinary polynomial regression of x . β is the vector of regression coefficients. In this paper, F ( β, x ) is simplified to a constant: F ( β, x ) = β , and all the following formulas are deduced on this basis. The second part z ( x ) is the realization of stochastic process, and it provides the approximation of local fluctuation.

z ( x)

shows

the

following statistical

characteristics, E  z ( x )  = 0

(6)

Var  z ( x )  = σ 2

(7)

Cov  z ( xi ) , z ( x j )  = σ 2  R ( xi , x j ) 

(8)

where σ 2 is the process variance. Eq.(8) defines the covariance between two arbitrary points xi and x j , R ( xi , x j ) is the spacial correlation function and can be defined as several models. The Gauss correlation function is adopted here, and it can be expressed as, ndv 2 R ( xi , x j ) = exp  − ∑ θ k xki − xkj   k =1 

(9)

where ndv is the number of design variables (also the dimension of sample points).

xki and xkj denote the k th component of x i and x j respectively. θ k is the

correlation parameter which ensures high flexibility of the metamodel, and it is usually obtained by optimization. Given a set of ns experimental samples X =  x1, x 2 ,K x ns 

T

and responses

T

Y =  y1, y2 ,K yns  , x i ∈ R ndv , yi ∈ R . The unknown parameters β and σ 2 can be estimated by: −1

βˆ = (1T R −11) 1T R −1Y T

σˆ

2

( Y − 1βˆ ) =

(

R −1 Y − 1βˆ

(10)

)

(11)

ns

where 1 is a vector filled with 1 of length ns , R is the correlation matrix with Ri , j = R ( x i , x j ) .

However, the correlation parameter θ k in Eq.(9) still has to be estimated before calculations of β and σ 2 . θ k can be obtained by optimizing a maximum ϕ .

ϕ =−

ns ln (σˆ 2 ) + ln R

(12)

2

With θk , β and σ 2 known, for an unknown point x new , the linear unbiased prediction is:

(

yˆ ( x new ) = βˆ + rT ( x new ) R −1 Y − 1βˆ

)

(13)

where rT ( x new ) is a ns -dimension row vector revealing the correlative relations between x new and experimental samples X =  x1 , x 2 ,L , x ns  . It can be expressed as,

rT ( x new ) =  R ( x new , x1 ) , R ( xnew , x 2 ) ,K , R xnew , xns   

(

)

T

(14)

Kriging variance σ y2ˆ ( x ) , which defines the variance of prediction yˆ ( x ) , can be expressed as,

(

−1

σ y2ˆ ( xnew ) =σ 2 1 + uT (1T R −11) u − rT ( xnew ) R −1r ( xnew )

)

(15)

where u = 1T R −1r ( x new ) − 1 .

σ y2ˆ ( x ) provides an important index to quantify the uncertainty of predictions and to further adjudge the fitting accuracy. Its existence supplies an approach to improve the design of experiment and to make Kriging more precise.

3. Proposed method combining Adaptive Importance Sampling and Kriging metamodel 3.1. Basic principles The method proposed in this paper can be treated as an improvement of the adaptive importance sampling. It combines Kriging metamodel with adaptive importance sampling in theory of active learning. Here, the basic fundamental is briefly summarized. The proposed method first generates a number of samples in the failure region using the above-mentioned Markov chain Metropolis algorithm. The samples introduced in Markov chain Metropolis algorithm can be sorted into two parts: N M Markov chain states and N R rejected samples. All these N M + N R samples are already carried into the actual performance function to evaluate the true responses, and they are subsequently used to construct an initial Kriging metamodel. Then, important samples (also the active learning candidate points) are generated by treating N M Markov chain states as sampling centers. Suppose there are N I samples

generated for each center, the population of these important samples is N I ⋅ N M . Their responses are efficiently predicted by Kriging rather than evaluated with the actual performance function. At the same time, in order to guarantee the predicting accuracy, an active learning mechanism is introduced as well. The active learning mechanism properly utilizes Kriging predictions and Kriging variances of important samples. It exactly determines the dangerous points which affect Kriging approximating accuracy most. Afterwards, by adding these points into previous design of experiment, the prediction accuracy of Kriging can be improved. The active

learning process runs iteratively until the convergence criteria are arrived. With the last-obtained Kriging, the N I ⋅ N M

important samples can be efficiently and

accurately predicted. Meanwhile, the failure probability can be estimated using these Kriging predictions with the theory of importance sampling reliability method. The proposed method inherits the advantages of adaptive importance sampling, active learning and Kriging metamodel, and enables evaluating the responses for the samples with more stochastic characters and high effects for Kriging approximation accuracy. Compared with other reliability methods, its computational efficiency is greatly increased while remaining similar accuracy. Detailed computational steps of this method are explained in Section 3.2.

3.2. Computational steps The proposed algorithm consists of the following steps, and its flowchart is shown in Figure 1. 1) Select an initial point in the failure region. The first step of the proposed method is to select an initial point for starting Markov chain. As mentioned above, the selection of this point which locates in the failure region does not present serious difficulties. 2) Generate more samples using Markov chain Metropolis algorithm. Based on the initial seed point, more samples which locate in the interested region can be generated using Markov chain Metropolis algorithm with one or more Markov chains. The stationary and proposal distributions are chosen as mentioned in Section 2.2. It should be noted that, with the generation of N M Markov chain states, N R samples are concomitantly rejected in the Metropolis algorithm, and the true responses of all these N M + N R samples are already evaluated in actual performance function. All these true evaluations will be used in the following steps and no excrescent calculations are introduced. The length of Markov chain is not necessary to be very long. According to our experiments, a hundred steps can be sufficient, that is

N M + N R = 100 . If the samples are generated by more than one chain, the results are

more satisfactory especially for reliability problems involving multiple design points. 3) Generate candidate important samples based on Markov chain states. In this step, the N M Markov chain states are used as importance sampling centers to generate more sample points. The sampling density function is constructed by moving the sampling center to the Markov chain states respectively. Unlike classic design point based importance sampling which has only one sampling center, the ISD of the present proposed method has N M sampling centers. Suppose N I samples are generated for each Markov chain state, there will be a total of N I ⋅ N M samples (denoted as SC ) generated in this step. For the reason that all the N M Markov chain states are located in the failure region, thus SC contains much more stochastic characters of the actual performance function. These N I ⋅ N M samples perform in two roles, one is important samples for reliability analysis, and the other one is candidate points for active learning mechanism. None of them is evaluated with the actual performance function in this step. 4) Build an initial Kriging metamodel. In order to obtain an initial Kriging metamodel, a corresponding initial design of experiment is required first. Here, we can just use previous N M + N R samples whose true responses were already evaluated to avoid unnecessary calculation wastes. Kriging metamodel can be constructed with the help of a MATLAB toolbox DACE, as this toolbox has been used in several reference articles [18, 20]. As mentioned in Section 2.3, the correlation model is chosen as Gaussian and the regression model is taken as constant. 5) Implement active learning process. With the Kriging metamodel built previously, the Kriging predictions and Kriging variances of population S C can be obtained by DACE. With these predictions and variances processed in a specified learning function referred to as U, the dangerous

point which affects the fitting accuracy can be identified. And then, its corresponding U value is compared with the stopping criterion, if the stopping criterion is satisfied, the active learning process stops and the algorithm goes to step 7, otherwise goes to step 6. It should be noted that in this step, the evaluations on U learning function do not need additional calls to the actual performance function, and it can be done only using Kriging predictions and variances. Further discussions of the learning function U including learning criterion and stopping condition can be found in following Section 3.3. 6) Update previous DOE and Kriging. If the stopping criterion is not satisfied, the previous Kriging metamodel is considered not to be precise enough to approximate the actual performance function. Consequently, previous DOE and Kriging should be updated. The dangerous point identified in step 5 is evaluated with the actual performance function firstly, and then the evaluation and the point itself are added to previous DOE. Afterwards, a new Kriging which is more precise can be rebuilt with DACE. Once the new Kriging is obtained, the algorithm goes back to step 5 to carry on the active learning process. Note that in this step, the actual performance function is evaluated only once. 7) Estimate failure probability and C.O.V. If the stopping criterion is satisfied, present Kriging is considered to be precise enough to estimate the failure probability of actual performance function. With the help of DACE, the evaluations of Kriging at population SC are predicted firstly, and then the failure probability and coefficient of variation (C.O.V.) of failure probability can be estimated. The C.O.V. is used to determine whether the importance sampling population SC is sufficient enough to ensure the convergence of the estimation of failure probability. In this paper, a C.O.V. less than 5% is considered to be precise enough. If the estimated C.O.V. is less than 5%, the algorithm goes to step 9, otherwise goes to step 8. Further information about how to estimate failure probability and C.O.V. from Kriging predictions will be discussed in following Section 3.4. 8) Update previous population.

If the estimated C.O.V. is larger than the requested criterion, population S C is considered too small to obtain an accurate result, and it needs to be enhanced. SC can be enlarged by generating another group of important samples like step 3 does and then adding them to present population. After the population is updated, the algorithm should go back to step 5 to check out whether the present Kriging is accurate enough to predict this new population, and the active learning process should carry on until stopping criterion is met again. It should be noted that in this step, no information about previous true evaluations is lost. This avoids the waste of computational resources and ensures efficiency of the method. 9) End the algorithm. If the estimated C.O.V. is small enough, the algorithm can stop and the last estimation of failure probability is considered as the final result of the method, which is accurate enough to evaluate the reliability of actual performance function.

3.3. Active learning mechanism It can be seen from above-mentioned sections that, to the proposed methodology, the active learning process is one of the important factors that affect the fitting accuracy. At the same time, lots of computations can be avoided owing to the existence of active learning process, and thus the computational efficiency is affected as well. The core of active learning mechanism is the learning function, and it identifies which points should be selected from candidate population to improve the fitting accuracy. Meanwhile, it determines when the active learning process can be stopped as the required stopping criterion is satisfied. Essentially, the proposed method can be classified as numerical simulating methods. To these kinds of methods, the signs of performance function are more important and concerned than the actual evaluations. The sign directly indicates where the sample drops: a positive value means the structure is safe in this configuration, and a negative value implies that the sample drops in the failure domain. If the signs of samples can be predicted exactly, the failure probability can be estimated

accurately as well. Compared with other samples which locate far from the limit state surface, a wrong predication for the samples near the limit state may cause the corresponding prediction to change from positive to negative (or negative to positive) and eventually to an unreliable failure probability. Whether the signs of the samples near the limit state can be correctly predicted affects the result directly. To identify these dangerous points close to the actual limit state surface and have a high potential effect on the approximating accuracy, the U learning function proposed by Echard [20] is adopted in this paper, and it can be expressed as, U ( x) =

Gˆ ( x )

σ Gˆ ( x )

(16)

where Gˆ ( x ) denotes the Kriging prediction of sample point x , and σ Gˆ ( x ) is the corresponding Kriging variance. Since Kriging makes the assumption of a Gaussian process, Φ (U ( x ) ) is the probability that x is classified in the correct group (failure of safe) according to the sign of its Kriging mean Gˆ ( x ) . The dangerous point which affects fitting accuracy most is the one which minimizes learning function U ( x ) . Meanwhile, a stopping criterion is needed to terminate the learning process. Echard [20] suggested that when the minimum evaluation of U ( x ) is bigger than 2 (in mathematic words, min U ( x ) ≥ 2, ( x ∈ SC ) ), the Kriging is considered accurate enough to approximate

the actual performance function. This stopping criterion is also adopted in this paper.

3.4. Estimate failure probability

As mentioned in Section 3.2, important samples SC is generated based on each Markov chain state respectively and each group of N I samples are generated according to a so-called sub-importance sampling density. SC is actually distributed to a weighted sum of each sub-importance sampling density, and the ISD used in the proposed method can be given by:

NM

hX ( x ) =∑ α j h j ( x )

(17)

j =1

where h j ( x ) is the j -th ISD component, α j is the weighted coefficient. In order to guarantee that hX ( x ) is a PDF, the condition

NM

∑α

j

= 1 is imposed. Since the

j =1

sample number of each ISD component is equal to N I , α j = 1/ N M . Replace the weighted importance sampling function hX ( x ) of Eq.(17) in Eq.(2), the failure probability is estimated as: Pf = ∫ n I F ( x ) R

= ∫ n IF (x) R

f X ( x) hX ( x ) dx hX ( x ) f X ( x ) NM ∑α j h j ( x )dx hX ( x ) j =1

(18)

NM  f (x)  = ∑ α j Eh j  I F ( x ) X  hX ( x )  j =1 

where Eh j [⋅] denotes expectation with respect to the j -th ISD component h j ( x ) . Then, the unbiased estimation of failure probability can be given by  1 NI NM f X ( xkj )  NM k Pˆf = ∑ α j  I x  = ∑ α j Pˆfhj ∑ F( j) k j =1 k = 1 N hX ( x j )  j =1  I

(19)

where x kj denotes the k -th sample generated according to j -th ISD component f X ( x kj ) 1 NI k ˆ is the failure probability evaluated by the h j ( x ) , and Pfh j = ∑ IF (x j ) N I k =1 hX ( x kj )

samples of h j ( x ) .

( )

The variance Var Pˆf

cannot be derived analytically. However, according to

Ref. [5], it can be approximated by Eq.(20)  1 NI NM f X ( x kj )  N M 2 2 k ˆ Var Pf ≈ ∑ α j Var   = ∑ α j Var Pˆfh j ∑ IF (x j ) j =1 hX ( x kj )  j =1  N I k =1

( )

where

( )

(20)

( )

Var Pˆfh j

2  f X ( x kj )  1  1 NI  k ≈  − Pˆfh j  ∑ IF (x j ) k N I − 1  N I k =1  h x ( ) X j    

( )

2

   

( j = 1, 2,K, N M )

(21)

Then the C.O.V. of Pˆf can be given as:

( )

Cov Pˆf ≈

( )

Var Pˆf Pˆf

(22)

4. Numerical examples In this section, several numerical examples are presented to illustrate the efficiency and accuracy of the proposed method. These examples are also analyzed with crude Monte-Carlo Simulation (MCS), classic design point-based Importance Sampling method (IS) and AK-MCS proposed by Echard [20] for comparison. 4.1. Example 1 This example is a linear reliability problem, and its limit state function is defined as g ( x ) = 2 − ( x1 + x2 ) / 2

(23)

where x1 and x2 are assumed to be independent with each other, and both of them obey the standard normal distribution. The exact solution of this problem can be obtained by Advanced First Order Reliability Method (AFORM). This solution will be used as an exact reference. All results are summarized in Table 1. In Table 1, AFORM, MCS, IS and AK-MCS correspond to the above reliability methods, and KAIS denotes the proposed Kriging-based Adaptive Importance Sampling method. The results are given in forms of failure probability Pf , the number of calls to actual performance function N call , the population number n , the

( )

C.O.V. of failure probability Cov Pˆf , the corresponding reliability index β and the error percentage ε β compared to reference reliability index. It should be noted that N call and ε β are two important indices which can be respectively used to

evaluate the efficiency and accuracy. For N call of AK-MCS, it is composed by initial design of experiment and the points introduced in active learning process, for example, “12+23” denotes 12 initial samples and 23 samples added in active learning process. Similarly, for N call of the proposed method, it is composed by Markov chain states, the samples rejected in Metropolis algorithm and the samples added in active learning process, i.e. “30+70+6” denotes 30 Markov chain states, 70 samples rejected and 6 samples added. Furthermore, for IS, it is supposed that the design point is identified prior, and the actual computational demand is larger than N call . It can be seen from Table 1, the proposed method exhibits a significant increase in computational efficiency compared with MCS and IS, and at the same time, the estimation of failure probability is very accurate. Compared with IS, although the population number n is larger, there are only 106 points carried into the actual performance function, and the rest are evaluated by Kriging constructed from these 106 points. Furthermore, for the proposed method, no more calculations are needed to find the design point or locate the importance region. As a result, the computational demand is much lower than IS. Compared with AK-MCS, the actual computational effort N call is a little larger. However, in these 106 samples, there are only 6 samples introduced in active learning process, whereas for AK-MCS, there are 23 samples added. At the same time, the population number n is much lower than that of AK-MCS, thus the proposed method still shows its great advantages in aspect of efficiency and accuracy. In order to illustrate the fitting accuracy more clearly, the actual and predicted limit state surfaces are compared in Figure 2. The Markov chain states, the rejected points, the important samples and the added points selected by active learning mechanism are also plotted in Figure 2. It can be seen that the 30 Markov chain states all locate in the failure region, and the important samples generated from these points represent the important regions very well. Meanwhile, for the active learning process, the added points all locate near the true limit state surface, which exhibits the

excellent learning ability of active learning mechanism. The Kriging metamodel constructed from these points reflects the characters of true limit state surface accurately. As Figure 2 shows, the true limit state surface (red solid line) and the predicted limit state surface (blue dotted line) are very close to each other. As a result, the Kriging metamodel can be used to accurately predict the responses of the important samples, and the number of calls to actual performance function N call can be greatly decreased as well.

4.2. Example 2 The previous example is a linear reliability problem. In this example, applicability of the proposed method to nonlinear problems will be examined. The example is taken from Ref. [18], and behaves nonlinearly around the design point. The performance function is given as g ( x ) = x13 + x12 x2 + x23 − 18

(24)

where x1 and x2 are independent and normally distributed with mean values 10 and 9.9, respectively, and a standard deviation of 5 for both. The analysis results are summarized in Table 2. The symbols in Table 2 have the same meanings with those in Table 1. Due to the nonlinearity, the failure probability and reliability index of MCS with C.O.V. 0.01 are selected as the exact solution. In the following examples, Monte-Carlo simulation result with C.O.V 0.01 will also be viewed as the exact solution. As can be seen from Table 2, although this problem behaves a certain degree of nonlinearity around the design point, the proposed method shows its excellent applicability and performs efficiently and accurately. For the proposed method, although there are only 8 samples accepted as Markov chain states, the important samples populate the important region very well and these samples reflect the statistical characters of the actual performance function perfectly. Compared with MCS and IS, the proposed method shows its high efficiency and the accuracy is acceptable as well. Whereas compared with AK-MCS, the advantage seems not so obvious, as the number of calls to actual performance function is a little

larger than AK-MCS. Considered that the added points in active learning process are less and the candidate important population is smaller, the method still shows its competitive strength. Compared with the previous example, there is almost no increase in N call . The reason lies in that the initial DOE is already sufficient for obtaining a relatively accurate Kriging metamodel, and only a small number of active learning processes can assure the fitting accuracy. Like example 1, the true and predicted limit state surfaces are also compared, and are plotted in Figure 3. It can be seen from Figure 3 that the 8 Markov chain states are mainly located near the design point, meanwhile, although this performance function behaves nonlinearly, the active learning mechanism still finds the dangerous points correctly, and the predicted limit state surface approximates the true limit state surface accurately.

4.3. Example 3 This example considers a series system with nonlinear performance functions modified from Refs. [3, 5], it is given as: g ( x ) = min ( g1 , g 2 )  x2   x  g1 ( x ) = 3 − 1 − x2 + exp  − 1  +  1   10   5  2 3 g 2 ( x ) = − x1 x2 2

4

(25)

where x1 and x2 are independent, identically and standard normally distributed. This system has two design points which have the same distance from the origin, i.e. ∗1 x ( ) = [ 0,3] and x∗( 2 ) = 3 / 2,3 / 2  . Both of them are important in the  

contribution of failure probability. The proposed method is applied to this example which has multiple important regions to show its efficiency and accuracy. The ‘exact’ failure probability is taken from Ref. [3], and is 2.53 × 10−3 . The results are

summarized in Table 3. It can be seen from Table 3 that the proposed method can efficiently handle the

problem of multiple design points with the premise of accurate result. Compared with MCS and IS, the proposed method obviously exhibits a higher efficiency, and the extra solution of design points is also avoided. Furthermore, compared with AK-MCS, the proposed method also shows its high efficiency with equivalent accuracy. Like previous examples, the samples of this example are also shown in Figure 4. It can be seen that the Markov chain states well simulate two important regions and hence a large part of candidate important samples locate in the important region and have more information. For the active learning process, the added points are mainly located near the limit state surface and thus the constructed Kriging metamodel shows its excellent fitting ability in the important region. It should be noted that for the bottom region, the Kriging model does not perform very well, but this affects the final result slightly. Considering that the probability of samples locating in that region is rather small, the imprecise fitting hardly affects the calculation result.

4.4. Example 4 The main purpose of introducing metamodel is to reduce the computational load when estimating the reliability of the structures involving implicit or time consuming performance functions. In order to test the applicability of the proposed method for problems with more random variables, the roof truss in Ref. [29] is selected in the example 4. The truss is simply illustrated as Figure 5. The top chords and compression bars of the truss are made by steel reinforced concrete, and the bottom chords and tension bars are made by steel. Assume the truss bears uniformly distributed load q , which can be transformed into nodal load P = ql / 4 . According to structural mechanics, the perpendicular deflection of truss peak node C ∆ C can be calculated as ∆ C =

ql 2 2

(

3.81 AC EC

+

1.13 AS ES

) , where AC , AS are the cross sectional areas of

reinforced concrete and steel bars, respectively, EC , ES are their corresponding elastic modulus, and l is the length of the truss. Considered the safety of the truss, the perpendicular deflection ∆ C should satisfy

the constraint ∆ C ≤ 0.03m . Hence, the performance function can be given as follows, g ( q, l , AS , AC , ES , EC ) = 0.03 − ∆ C = 0.03 −

ql 2 3.81 1.13 ( ) + 2 AC EC AS ES

(26)

All random variables are independent and normally distributed, and the statistical properties are given in Table 4. The reliability results are listed in Table 5. It can be clearly seen that the proposed method still shows high efficiency and accuracy when dealing with this truss structural problem. With the intervention of Markov chain Metropolis algorithm, active learning mechanism and Kriging metamodel, the proposed method only evaluates 215 sample points with the actual performance function. Whereas, when adopting MCS or IS, the same exact reliability estimation is obtained after 50000 or 2000 calls to the actual performance function.

4.5. Example 5 A three-span beam structure shown in Figure 6 is selected as another engineering problem. The beam subjects to a uniformly distributed load w , and the length of a single span L is determinate, L = 5m . The maximum vertical deflection of the beam can be expressed as ∆ max = 0.0069 wL4 / EI . Given the maximum allowable vertical deflection

L / 400

as constraint, the performance function can be

established: g ( w, E , I ) = L / 400 − 0.0069 wL4 / EI

(27)

where w denotes the uniformly distributed load, E denotes elastic modulus of the beam, and I denotes the moment of inertia. All three random variables are normally distributed and independent. Their corresponding statistical properties are given in Table 6. The reliability results are summarized in Table 7. The results of this problem still illustrate the computational efficiency and accuracy of the proposed method. From Table 5 and Table 7, it can be concluded that the proposed method can be used for these kinds of engineering problems with accurate results.

5. Conclusions

This paper proposes an efficient reliability method, and it can be seen as an improvement of the adaptive importance sampling method based on Kriging metamodel and active learning mechanism. The proposed method inherits the superiorities of Kriging metamodel, adaptive importance sampling and active learning mechanism. With the Markov chain Metropolis Algorithm, the sampling centers can be efficiently generated and the candidate important samples can populate the important regions very well. Due to the introduction of an active learning process, the points which may greatly affect the fitting accuracy can be precisely determined with a rigorous criteria, this makes the Kriging metamodel more accurate. With the Kriging metamodel, the evaluations of a large part of candidate samples can be efficiently predicted instead of calculating in actual time-consuming performance function, thus the calculation efficiency can be greatly increased. With application of the proposed method, there are only a few sample points evaluated with the actual performance function, while the majority of samples with little effect are predicted by Kriging metamodel. The proposed method avoids the evaluations for the samples whose exact responses are unnecessary, the boring calculations for design points and the waste of computational resources, thus the calculation efficiency greatly increases. Results of validating examples illustrate that the proposed method performs very well either in explicit reliability problems or implicit ones. Furthermore, the proposed method can also be adopted in nonlinear, high dimensional or non-normal distributed problems, and shows great prospective application.

Acknowledgements The financial supports from Chinese Aerospace Support Foundation (Grant No. NBXW0001), National Natural Science Foundation of China (Grant No. 51205312) and Foundation Research Funds of Northwestern Polytechnical University (Grant No. JCY20130123) are gratefully acknowledged, and the pertinent suggestions from anonymous reviewers are also acknowledged by authors.

References

[1] O. Ditlevsen, H.O. Madsen, Structural reliability methods, John Wiley & Sons, Chichester, 1996. [2] V. Bayer, C. Bucher, Importance sampling for first passage problems of nonlinear structures, Probab. Eng. Mech. 14 (1999) 27-32. [3] S.K. Au, J.L. Beck, A new adaptive importance sampling scheme for reliability calculations, Struct. Saf. 21 (1999) 135-138. [4] F. Zhang, Z.Z. Lu, L.J. Cui, S.S. Song, Reliability sensitivity algorithm based on stratified importance sampling method for multiple failure modes systems, Chinese J. Aeronaut. 23 (2010) 660-669. [5] X.K. Yuan, Z.Z. Lu, C.C. Zhou, et al., A novel adaptive importance sampling algorithm based on Markov chain and low-discrepancy sequence, Aerosp. Sci. Technol. 19 (2013) 253-261. [6] C.G. Bucher, U. Bourgand, A fast and efficient response surface approach for structural reliability problems, Struct. Saf. 7 (1990) 57-66. [7] M. Kleiber, J. Knabel, J. Rojek, Response surface method for probabilistic assessment of metal forming failures, Int. J. Numer. Meth. Eng. 60 (2004) 51-67. [8] I. Kaymaz, C.A. McMahon, A response surface method based on weighted regression for structural reliability analysis, Probab. Eng. Mech. 20 (2005) 11-17. [9] D.L. Allaix, V.I. Carbone, An improvement of the response surface method, Struct. Saf. 33 (2011) 165-172. [10] J. Deng, Structural reliability analysis for implicit performance function using radial basis function network, Int. J. Solids Struct. 43 (2006) 3255-3291. [11] J. Cheng, J. Zhang, C.S. Cai, et al., A new approach for solving inverse reliability problems with implicit response functions, Eng. Struct. 29 (2007) 71-79. [12] V. Papadopoulos, D.G. Giovanis, N.D. Lagaros, et al., Accelerated subset simulation with neural networks for reliability analysis, Comput. Methods Appl. Mech. Eng. 223-224 (2012) 70-80. [13] G.B. Kingston, M. Rajabalinejad, B.P. Gouldby, et al., Computational intelligence methods for the efficient reliability analysis of complex flood defence structures, Struct. Saf. 33 (2011) 64-73.

[14] J.E. Hurtado, An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory, Struct. Saf. 26 (2004) 271-293. [15] A. Basudhar, S. Missoum, Adaptive explicit decision functions for probabilistic design and optimization using support vector machines, Comput. Struct. 86 (2008) 1904-1917. [16] J.M. Bourinet, F. Deheeger, M. Lemaire, Assessing small failure probabilities by combined subset simulation and Support Vector Machines, Struct. Saf. 33 (2011) 343-353. [17] C.M. Rocco, J.A. Moreno, Fast Monte Carlo reliability evaluation using support vector machine, Reliab. Eng. Syst. Saf. 76 (2002) 237-243. [18] I. Kaymaz, Application of kriging method to structural reliability problems, Struct. Saf. 27 (2005) 133–151. [19] B.J. Bichon, M.S. Eldred, L.P. Swiler, et al., Efficient global reliability analysis for nonlinear implicit performance functions, AIAA J. 46 (10) (2008) 2459-2468. [20] B. Echard, N. Gayton, M. Lemaire, AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation, Struct. Saf. 33 (2011) 145-154. [21] B. Echard, N. Gayton, M. Lemaire, et al., A combined importance sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models, Reliab. Eng. Syst. Saf. 111 (2013) 232-240. [22] B.H. Ju, B.C. Lee. Reliability-based design optimization using a moment method and a kriging metamodel, Eng. Optim. 40 (2008) 421-438. [23] J.P.C. Kleijnen, Kriging metamodeling in simulation: A review, Eur. J. Oper. Res. 192 (3) (2009) 707-716. [24] A. Marrel, B. Iooss, F.V. Dorpe, et al., An efficient methodology for modeling complex computer codes with Gaussian processes, Comput. Stat. Data Anal. 52 (2008) 4731-4744. [25] S. Sakata, F. Ashida, M. Zako, Structural optimization using Kriging approximation, Comput. Methods Appl. Mech. Eng. 192 (2003) 923-939.

[26] K.W. Juang, W.J. Liao, T.L. Liu, et al., Additional sampling based on regulation threshold and kriging variance to reduce the probability of false delineation in a contaminated site, Sci. Total Environ. 389 (1) (2008) 20-28. [27] A.A. Giunta, L.T. Watson, A comparison of approximation modeling techniques: polynomial versus Interpolating models, AIAA-98-4758. [28] S.N. Lophaven, H.B. Nielsen, J. Sondergaard, DACE, a Matlab Kriging Toolbox, 2002. [29] S.S. Song, Z.Z. Lu, H.W. Qiao, Subset simulation for structural reliability sensitivity analysis, Reliab. Eng. Syst. Saf. 94 (2009) 658-665.

Table Captions

Table 1

Reliability results for Example 1

Table 2

Reliability results for Example 2

Table 3

Reliability results for Example 3

Table 4

Statistical properties of random variables for roof truss

Table 5

Reliability results for Example 4

Table 6

Statistical properties of random variables for three-span beam

Table 7

Reliability results for Example 5

Table 1

Reliability results for Example 1

Pf

N call

n

AFOSM

0.02275

--

MCS

0.02273

IS

( )

β

εβ (%)

--

--

2

--

440000

440000

0.01

2.0004

0.02

0.02266

24000

24000

0.01

2.0017

0.09

AK-MCS 0.02259

12+23

440000

0.01

2.0030

0.15

39000

0.01

1.9796

1.02

KAIS

0.02388 30+70+6

Cov Pˆf

MCS

Table 2

Reliability results for Example 2

Pf

N call

n

( )

Cov Pˆf

β

εβ (%)

0.005798 1720000 1720000

0.01

2.5242

--

0.005816

42000

42000

0.01

2.5231

0.04

AK-MCS 0.005783

12+33

1720000

0.01

2.5251

0.04

40000

0.01

2.5198

0.17

IS

KAIS

0.005871 8+92+10

Table 3

Reliability results for Example 3

εβ (%)

0.01

2.8033

--

160000

0.05

2.7959

0.26

2000

2000

0.05

2.7673

1.28

12+172

170000

0.05

2.8178

0.52

7000

0.05

2.7834

0.71

N call

n

MCS

0.002529

3950000

3950000

MCS

0.002588

160000

IS

0.002826

AK-MCS 0.002418 KAIS

( )

Cov Pˆf

β

Pf

0.002690 35+65+63

Table 4 Random

Statistical properties of random variables for roof truss Distribution

Standard

Coefficient of

deviation

variation

Mean variable

type

q( N / m)

Normal

20000

1400

0.07

l(m)

Normal

12

0.12

0.01

AS ( m 2 )

Normal

9.82×10-4

5.982×10-5

0.06

AC ( m 2 )

Normal

0.04

0.0048

0.12

ES ( Pa )

Normal

1×1011

6×109

0.06

EC ( Pa )

Normal

2×1010

1.2×109

0.06

Table 5

Reliability results for Example 4

εβ (%)

0.01

2.3497

--

50000

0.05

2.3503

0.03

2000

2000

0.05

2.3510

0.06

12+255

50000

0.05

2.3471

0.11

19600

0.05

2.3509

0.05

N call

n

MCS

0.009395

1060000

1060000

MCS

0.009380

50000

IS

0.009361

AK-MCS 0.009460 KAIS

( )

Cov Pˆf

β

Pf

0.009365 28+72+115

Table 6 Random

Statistical properties of random variables for three-span beam Distribution

Standard

Coefficient of

deviation

variation

Mean variable

type

w ( KN / m )

Normal

10

0.4

0.04

E ( KN / m 2 )

Normal

2×107

0.5×107

0.25

I ( m4 )

Normal

8×10-4

1.5×10-4

0.1875

Table 7

Reliability results for Example 5

εβ (%)

0.01

3.0282

--

340000

0.05

3.0334

0.17

3000

3000

0.05

3.0282

0

12+341

460000

0.05

3.0059

0.73

2100

0.05

3.0231

0.17

N call

n

MCS

0.001230

8120000

8120000

MCS

0.001209

340000

IS

0.001230

AK-MCS 0.001324 KAIS

( )

Cov Pˆf

β

Pf

0.001251 21+79+209

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6