A novel adaptive importance sampling algorithm based on Markov chain and low-discrepancy sequence

Aerospace Science and Technology 29 (2013) 253–261

Contents lists available at SciVerse ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

A novel adaptive importance sampling algorithm based on Markov chain and low-discrepancy sequence Xiukai Yuan a,∗ , Zhenzhou Lu b , Changcong Zhou b , Zhufeng Yue c a b c

Department of Aeronautics, Xiamen University, Xiamen 361005, PR China School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, PR China School of Mechanics, Civil Engineering & Architecture, Northwestern Polytechnical University, Xi’an 710072, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 11 January 2012 Received in revised form 18 March 2013 Accepted 22 March 2013 Available online 29 March 2013 Keywords: Reliability Importance sampling Markov chain Low-discrepancy sequence Monte Carlo simulation

A novel adaptive importance sampling method is proposed to estimate the structural failure probability. It properly utilizes Markov chain algorithm to form an adaptive importance sampling procedure. The main concept is suggesting the proposal distributions of Markov chain as the importance sampling density. Markov chain states can adaptively populate the important failure regions thus the importance sampling based on them will yield an efficient and accurate estimate of the failure probability. Compared with existent methods, it does not need the solution of the design point(s) or the pre-sampling in the failure region. Various examples are given to demonstrate the advantages of the proposed method. © 2013 Elsevier Masson SAS. All rights reserved.

1. Introduction In structural reliability analysis, the failure probability can be defined as the following multifold probability integral:

 PF =

 ···

 f (x) dx =

F

 ···

I F (x) f (x) dx

(1)

Rn

where x = [x1 , x2 , . . . , xn ]T is the n-dimensional vector of random variables; f (x) is the joint probability density function (PDF) of x and I F (x) is the indicator function of the failure region F = {x: g (x)  0} (where g (x) is the limit state function): I F (x) = 1 if g (x)  0 and I F (x) = 0 otherwise. In order to evaluate the failure probability P F , it is necessary to compute an integral over n-dimension space. However, in practice, it is difficult to evaluate P F when the number of the dimensions n is large or the failure region F has a complicated shape. In addition, the efficiency of the computing method is the main concern when dealing with a practical problem, as reliability analysis with finite element model is usually quite time-consuming. Several methods have been developed in order to solve the integral. Among which Monte Carlo simulation (MCS) [25] method is one of the most widely used methods for handling complex

*

Corresponding author. E-mail addresses: [email protected] (X. Yuan), [email protected] (Z. Lu). 1270-9638/$ – see front matter © 2013 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ast.2013.03.008

problems. However, it shows a poor computational efficiency in dealing with problems of small failure probability or problems with costly finite element analyses. Hence, many variance reduction techniques have been developed to address this disadvantage, such as importance sampling [1,4,5,8,10–13,15,18,20, 21], subset simulation [6,7], and line sampling technique [23,24], etc. Importance sampling techniques have been developed over the past few decades, and have been one of the most prevalent approaches in the structural reliability analysis. Its basic idea is to shift underlying distribution towards the failure region to gain information from rare events more efficiently. It draws samples of the random variables from a pre-selected distribution which is concentrated in the important failure region. The success of the method relies on the proper choice of the importance sampling density (ISD). The optimal choice for the ISD is practically infeasible since it requires knowledge of failure probability a priori. Several methods have been developed to construct the ISD, which is the approximate optimal ISD or a different one. At present, the most prevalent choices are those based on design point(s) [5,8,18] or kernel density estimators [1,4,10–13,15,20,21]. A design point is usually obtained by solving a constrained optimization problem. When multiple design points exist, the search for them may be difficult as it requires more sophisticated algorithms for the optimization problem. Thus in an importance sampling process, the search for design point(s) may occupy a great portion of the total computational effort.

254

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

The kernel sampling density (KSD) is constructed based on the samples in the failure region which are simulated according to f (x). It tends to the optimal ISD as the number of samples increases, thus it can be served as a good choice for the ISD. However, it needs pre-sampling in the failure region F . The points in F can be generated using standard rejection sampling from the original PDF [1], however this is extremely inefficient in cases where the failure probability is small. Another method which adopts Markov chain Metropolis algorithm [19] can overcome this drawback [4]. In this method, a number of points in F are simulated as the states of a Markov chain and their PDF asymptotically tends to the optimal ISD as the number of Markov chain steps increases. Subsequently, a KSD estimator can be constructed based on these points. The total computational cost of this approach mainly consists of the simulation of Markov chain samples in order to construct the KSD and the subsequent simulation of importance sampling in the estimation of the failure probability. Some other adaptive importance sampling strategies have also been developed, i.e., the non-parametric adaptive ISD [20] and the cross-entropy (CE) method [10,11,13]. In [20], a non-parametric adaptive importance sampling procedure was proposed which adapts the auxiliary PDF (weight and kernels) by using an algorithm at each iteration to approach the optimal ISD. The cross-entropy method [10,11,13] can be seen as an adaptive importance sampling procedure. The cross-entropy is a measure of closeness between two sampling distributions. The basic idea of the cross-entropy method is to choose the ISD in a specified class of densities such that the cross-entropy between the optimal ISD and implementing ISD is minimal. These methods provide alternative and adaptive ways of approximating the optimal ISD. Other strategies can also be found in literatures to further improve the performance of importance sampling, i.e., adopting the splitting strategy as given in [21] or the surrogate model (Kriging model) as in [15]. Recently, quasi-Monte Carlo (QMC) has gained more and more popularity in studying multidimensional integration [9,14,22,26]. Compared with the traditional MCS method, quasi-Monte Carlo method replaces the pseudo-random numbers used in MCS by low-discrepancy sequences. Sobol [26] drew the conclusion that applying quasi-Monte Carlo is more promising than Monte Carlo algorithms and Bosserta et al. [9] also pointed out that the numerical integration converges faster with a higher accuracy by using quasi-random numbers compared with using pseudo-random numbers. The quasi-Monte Carlo method which is also called the low-discrepancy sampling method has been applied in many fields, i.e. computational physics, statistics and mathematical finance. Recently, Nie and Ellingwood [22] and Dai and Wang [14] applied it to structural reliability analysis. Various low-discrepancy sequences were investigated in both Refs. [22] and [14]. Dai and Wang successfully combined the low-discrepancy sequences with traditional importance sampling technique to further improve the computational efficiency. In this paper, a novel importance sampling approach, called Markov chain importance sampling method (MCIS), is proposed to efficiently handle structural reliability analysis. The innovative part of this contribution is that the proposed method modify the Markov chain Metropolis algorithm into an adaptive importance sampling procedure, thus the proposed method combines the advantages of both the Markov chain Metropolis algorithm and the importance sampling method. The basic idea of the proposed method is to take the proposal distribution of Markov chain as the ISD, and the Markov chain time steps are selected as the importance sampling centers. According to the property of the Metropolis algorithm, Markov chain states can adaptively populate the important region, thus the proposed importance sampling based on these states will yield an efficient and accurate estimate of the failure probability. The most desirable feature is that it can complete

the importance sampling with only an initial point, no more information is needed, e.g. design point(s) or pre-samples. Furthermore, the low-discrepancy sequence is adopted in the proposed method to further improve the accuracy. 2. Original Markov chain Metropolis algorithm Markov chain, also called Metropolis–Hastings algorithm, is a very useful and powerful tool for simulating samples according to an arbitrary PDF and it has become more and more widely used in structural reliability analysis [4,6]. In the Metropolis–Hastings algorithm, the samples are simulated as the states of a special Markov chain whose limiting stationary distribution is equal to the target PDF. This means that the PDF of the samples, which are Markov chain states, converges to a unique limit, the target PDF, as the number of Markov steps increases. In structural reliability analysis, this algorithm can be used to efficiently simulate the samples in the failure region and the samples can adaptively populate the important failure region. Taking the optimal sampling density f (x| F ) = I F (x) f (x)/ P F as the stationary distribution, the Markov chain states {x1 , . . . , x M } (M is the total number of Markov chain time steps) are asymptotically distributed as f (x| F ). See [4, 6,19] and [3] for more details of applying the Metropolis–Hastings algorithm to reliability analysis. The original Markov chain procedure for simulating samples {x1 , . . . , xM } with limit stationary distribution equal to target PDF f (x| F ) is briefly described below. (1) Select an initial state of Markov chain. A ‘seed’ point x1 should be selected in the failure region F either randomly or deterministically to start the sequence. It may be obtained as a point simulated according to the PDF f (x), or be assigned based on engineering judgment as a point that leads to system failure. In most structural systems of practical interest the determination of the point that lead to failure does not present serious difficulties. Even if the current sample is not distributed as f (x| F ), the limit distribution property of Markov chain guarantees that the distribution of simulated samples will tend to f (x| F ) as the number of Markov chain steps increases [3,4]. (2) Determine the ( j + 1)-th time step state x j +1 of Markov chain. Suppose the Markov chain is in j-th time step x j , then generate a new point ε from a chosen PDF, f ∗ (ε |x j ), called the ‘proposal distribution’, which is a PDF for ε depending on x j . Calculate the ratio r = f (ε | F )/ f (x j | F ). If r > 1, then ε is accepted as the ( j + 1)-th state, that is x j +1 = ε , otherwise ε is accepted as the ( j + 1)-th state with probability r and x j is accepted as the ( j + 1)-th state with the remaining probability 1 − r, i.e., x j +1 = x j , where a repeated state is obtained. (3) Repeat step (2) until M Markov chain states, {x1 , . . . , x M }, are obtained. 3. Markov chain importance sampling method The proposed method is to modify the original Markov chain Metropolis algorithm to establish an adaptive importance sampling procedure. It takes the proposal distribution of Markov chain as the ISD, and the Markov chain time step states are selected as the sampling centers, thus it can properly combine the advantages of both the Markov chain Metropolis algorithm and the importance sampling method. The Markov chain states can adaptively populate the importance failure region, thus the importance sampling centered on these states owns high efficiency in exploring the important failure region. The procedure of the proposed method is illustrated firstly. It stems form the original Markov chain algorithm, but several mod-

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

255

ifications are made here to establish an adaptive importance sampling procedure. (1) Select an initial point as the first sampling center. An initial point x1 in the failure region is selected to start the Markov chain, and it is also the first importance sampling center. (2) Generate samples based on j-th sampling center. Based on the j-th ( j = 1, . . . , M) importance sampling center x j which is also the j-th time step of Markov chain, the proposal distribution f ∗ (ε |x j ) is taken as the j-th ISD component to generate (1)

(N j )

(2)

a set of samples, i.e. {x j , x j , . . . , x j }, where N j is the number of samples for the j-th sampling. Note that this is different form the original Metropolis–Hastings algorithm, where only one sample is generated at each time step. These samples are not only the important samples but also can be seen as the candidate states of ( j + 1)-th Markov chain time step state. Thus it seems that based on x j , not one but N j Markov chains are generated. (3) Determine the next sampling center. The next step of the proposed method is to choose the next sampling center which is also the next time step of Markov chain. After the proposal distribution generates N j sample at a time in j-th sampling, the ( j + 1)-th sampling center will be selected from these samples. The strategy is illustrated as follows. First, use the Metropolis algorithm to determine the ‘next state’ of Markov chain for each sample of the j-th set of sam(1)

In the procedure described above, based on each sampling center (Markov chain time step state) x j , a number of samples, N j , are generated by using the proposal distribution (ISD) f ∗ (ε |x j ). Thus the total ISD H (x) used in the proposed method is actually the weighted sum of the proposal distributions which can be given by:

H (x) =

(k)

puted. If rk > 1, x j

is accepted as the ‘pseudo’ ( j + 1)-th state,

(k)

= x j , otherwise x(jk) is accepted as the ‘pseudo’ ( j + 1)-th

ε (jk+)1 = x(jk) , and x( j) is accepted with (k) the remaining probability 1 − rk , ε j +1 = x( j ) . As a result, there is state with probability rk ,

one ( j + 1)-th Markov chain state corresponding to each sample (k) x j (k = 1, 2, . . . , N j ), thus a total of N j ‘pseudo’ ( j + 1)-th Markov (1)

(2)

(N j )

chain states, {ε j +1 , ε j +1 , . . . , ε j +1 }, are obtained. Second, choose one of these ‘pseudo’ ( j + 1)-th states as the next sampling center. The next (( j + 1)-th) sampling center is also the ( j + 1)-th time step of the modified Markov chain. Note that repeated samples may be obtained in Markov chain, slowing down the sample coverage of the importance regions of the failure domains. Therefore, among these N j ( j + 1)-th ‘pseudo’ Markov (1)

(2)

(N j )

chain states, {ε j +1 , ε j +1 , . . . , ε j +1 }, the one that is different from current state x j , should be preferred to be selected as the next sampling center. This leads to acceleration of the convergence of (k) (k) Markov chain, i.e., x j +1 = ε j +1 if ε j +1 = x j (1  k  N j ). If there (1)

(2)

(N )

are several such points, i.e., {ε j +1 , ε j +1 , . . . , ε j +d1 } (N d  N j is the number of points which are different from current state), then one of them can be randomly selected as the ‘real’ ( j + 1)-th state which is the next sampling center x j +1 . (4) Calculate the failure probability. Repeat steps (2)–(3) until M Markov chain time steps (sampling centers) are reached {x1 , x2 , . . . , x M }. And based on each time step, N j samples are generated. Then the failure probability of the problem can be calculated as follows. The schematic diagram of the proposed procedure for twodimensional case is shown in Fig. 1. When the low-discrepancy sequence is used in step (2) to generate samples, it is called ‘Markov chain importance sampling method with quasi-Monte Carlo simulation’, and denoted as ‘MCIS-QMC’. When the pseudo-random numbers are used, it is called ‘Markov chain importance sampling method’ and denoted as ‘MCIS’.

M 

α j h j ( x) =

j =1

(N j )

(2)

ples {x j , x j , . . . , x j }. Note that these states are called ‘pseudo’ ( j + 1)-th states here, as just one of them is selected as the ‘real’ next time step of Markov chain. That is, for each sample, i.e., (k) (k) x j (k = 1, 2, . . . , N j ), the ratio rk = f (x j | F )/ f (x j | F ) is com-

ε(jk+)1

Fig. 1. Schematic representation of the procedure for the proposed method.

M 

α j f ∗ ( x| x j )

(2)

j =1

where h j (x) = f ∗ (x|x j ) is the j-th ISD component, α j is the corresponding weighted coefficient of the j-th ISD component. In order M to guarantee that H (x) is a PDF, the condition j =1 α j = 1 is imposed. When each ISD component generate the same number of samples, then α j = 1/ M. Using the weighted importance sampling function H (x) in Eq. (2), the failure probability can be rewritten as:

 PF =

 ···

I F ( x) Rn

=

M 



=

H ( x)



α j Eh j

α j h j (x) dx

j =1

 ···

αj

j =1 M 

M f (x) 

I F ( x) Rn

I F (x) f (x)

f (x) H ( x)

 h j (x) dx

 (3)

H (x)

j =1

where E h j [·] denotes expectation with respect to the ISD component h j (x). Suppose that the M sampling centers are used, and based on (1)

(2)

(N j )

each center x j ( j = 1, . . . , M ), N j samples, {x j , x j , . . . , x j }, are generated by the j-th ISD, i.e. h j (ε ) = f ∗ (ε |x j ). Then according to Eq. (3) the unbiased estimate of failure probability Pˆ F can be given by

Pˆ F =

M 



αj

j =1

where Pˆ F h j =

N

(k) (k) j 1  I F (x j ) f (x j )

Nj 1 Nj

k =1

N j k=1

(k)

H (x j ) (k)

 =

M 

α j Pˆ F h j

(4)

j =1

(k)

I F (x j ) f (x j ) (k)

H (x j )

is the estimate of the failure

probability evaluated by the samples of the j-th ISD. 4. Low-discrepancy sampling The low-discrepancy sampling method can significantly increase the accuracy and efficiency of MCS for structural reliability analysis by replacing the pseudo-random numbers with lowdiscrepancy sequences. In this paper, the low-discrepancy sampling

256

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

method is incorporated into Markov chain importance sampling simulation in order to further improve the accuracy and efficiency of the proposed method. In reliability analysis, the success of MCS method lies in the samples used, i.e. the better the samples approximate the uniform distribution the more precise the integral estimate is. However, the distribution samples are usually obtained by the pseudo-random numbers which commonly achieve a poor representation in a finite number of samples. Due to the Koksma–Hlawka inequality [16], a sequence of points with the smallest discrepancy is the best set of integration points. These low-discrepancy sequences have the best uniformity on the integration domain thus a faster rate of convergence and a precise integral estimate can be obtained. The low-discrepancy sampling method overcomes the MCS approximation by using deterministic sample points that fill the space as fully as possible, as it has deterministic and smaller error bound, thus, it implies that low-discrepancy sampling method has faster convergence rate than MCS method. See [14,16,27] and [2] for more details of the low-discrepancy sequences. There are several classes of sequences of points that are well suited to multivariate integration [14], e.g. lattice points and digital sequences. The Good lattice point set and the Hua–Wang point set belong to lattice points. The Halton and Hammersley sequences are the representative of the digital sequences which are popularly used [2,27]. In this paper, the Halton sequence is used.

5.3. Choice of parameter α j The parameter α j ( j = 1, . . . , M ) is the weight coefficient of the j-th ISD component and its value affects the estimate of the failure probability to a certain extent. As the same number of samples are generated for each ISD component in the proposed method, α j should be assigned the same value, that is, α j = 1/ M. 5.4. Statistical property of the estimate Here, the mean and standard deviation of the estimate are investigated. Note that, the proposed method is different from the traditional importance sampling method, as the samples obtained by the proposed method are correlated samples. That is, each sampling center is correlated with the samples generated in the following time steps. First, the estimate obtained by the proposed method is unbiased as it seen form Eq. (6)

E ( Pˆ F ) =

M 

5.1. Choice of proposal PDF The proposal distribution f ∗ (ε |x j )

=

M 

( j = 1, . . . , M ) is a PDF for

distribution and the uniform distribution over an n-dimensional hypercube are possible candidates. As the proposal PDF is taken as the importance sampling PDF and the Gaussian distribution is a popular choice for ISD, herein the Gaussian distribution is chosen as the proposal PDF. In general, a product of n-dimensional Gaussian PDFs with mean x ji (i = 1, . . . , n) and standard deviation σi (i = 1, . . . , n) is given by:

i =1

√

1

2πσi

Nj

(k)

(k)



α j Eh j



(k)

H (x j )

k =1

I F (x) f (x)



H (x)

j =1

Var( Pˆ F ) ≈

 exp −

2

(εi − x ji ) 2σ i

(5)

The parameter σi affects the distance the next sample point can depart from the actual one, thus it governs the efficiency in generating samples by Markov chain. Generally speaking, too large standard deviation in the proposal PDF may lead to too many rejections in the Metropolis algorithm and yield a highly-correlated chain. In this paper, the parameter of σi is selected to be equal to the standard deviation of xi .

M 



α

2 j Var

j =1

ε centered at x j , and the proposal PDF is chosen to be symmetric with respect to ε and x j , for example, the n-dimensional Gaussian

f ∗ (ε |x j ) =

αjE

j =1

N

j 1  I F (x j ) f (x j )

= PF

(6)

The variance Var( Pˆ F ) of the estimate Pˆ F in Eq. (4) can’t be derived analytically. However, it can be approximated by Eq. (7) when the correlation between the samples is regardless.

5. Implementation issues

n



=

M 

N

(k) (k) j 1  I F (x j ) f (x j )

Nj

k =1



(k)

H (x j )

α 2j Var( Pˆ F h j )

(7)

j =1

where

Var[ Pˆ F h j ] ≈



1 Nj −1



N

(k)

(k)

j 1  I F (x j ) f (x j )

Nj

k =1

(k)

H j (x j )

2

 − ( Pˆ F h j )

j = 1, . . . , M

2

(8)

Then the coefficient of variation (c.o.v.) of the estimate Pˆ F , Cov( Pˆ F ), can also be given as follow:

Cov( Pˆ F ) ≈

Var( Pˆ F ) Pˆ F

 M  α 2j Var( Pˆ F h j )  = Pˆ 2 j =1

F

 M   2  = α 2j Cov( Pˆ F h j )

(9)

j =1



5.2. Choices of M and N j

where Cov( Pˆ F h j ) =

Var( Pˆ F h j )/ Pˆ F is the c.o.v. of Pˆ F h j .

For convenience, the same value can be assigned to N j , that is, N 1 = N 2 = · · · = N M = N, then the total number of samples (evaluations of limit state function) used in an analysis is N T = M × N. For a fixed amount of computational effort N T , it is suggested to select a small value for N, i.e. 5 or 10, and a relatively large value for M (i.e., 100 or more) to guarantee the convergence of Markov chain as M represents the number of time steps of Markov chain. This will lead to the reduction of the variance of the estimate Pˆ F .

5.5. Comparison with other existent methods It should be noted that there are some obvious differences between the proposed method and that by Au and Beck [4]. The first one is the way of generation of samples. In the method by Au [4], the original Markov chain is used to pre-sample in the failure region F and then a KSD for importance sampling is constructed based on these samples, where the estimation of kernel density parameters is involved. In the proposed method here, the original

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

Markov chain algorithm is modified and the proposal distribution is taken as the ISD to generate not one but N j samples at a time. It seems that both the modified Markov chain algorithm and importance sampling are carried out at the same time, that is, there is no need to pre-sample and construct a KSD for importance sampling. Second, in the proposed method, the sample point which is distinct from the current Markov chain state is prior to be selected as the next time step of Markov chain, thus it seems that Markov chain samples converges to the important failure region more quickly. Since the convergence of Markov chain is crucial for the performance of ISD, the proposed method is more likely to perform better than the method of Au and Beck while the same number of Markov chain time steps is used. This can be seen from the examples given in the following section. The third difference is that the low-discrepancy sequence is adopted in the proposed method to further improve the accuracy of the method. The proposed method contains multiple ISDs, and for each ISD component, only a few samples are generated, thus the low-discrepancy sequence is used to achieve a better estimate with a few samples due to the merits of low-discrepancy sequence in expediting convergence and increasing solution stability for structural reliability. Thus, the proposed method can be seen as an improvement of the method of Au and Beck. Meanwhile, the importance sampling method based on the kernel sampling density (IS-KSD) with fixed width is used in following section for comparison. Note that the proposed method is also different with the method proposed in [3], although both of them are related with importance sampling and Markov chain algorithm. Actually, the differences between them are obvious. First, the aim of the method proposed in [3] is probabilistic analysis for failure analysis but not the failure probability estimation which is the aim of the proposed method in this contribution. In [3], it studies the statistic of the system behavior after the reliability analysis has been made by importance sampling. Second, the procedure in [3] is used to efficiently generate the conditional samples for failure analysis. When the reliability analysis is carried out by importance sampling, it make use of the samples obtained form importance sampling, and transform the importance samples into the conditional samples by Markov chain algorithm and then use them to yield estimate of conditional expectations of interest for failure analysis, thus it is obviously different from the proposed method in this work. 6. Illustrative examples The proposed importance sampling methods, MCIS and MCISQMC, are applied to different problems to investigate their accuracy, efficiency and robustness. Note that only one simulation run cannot demonstrate the performance of the proposed methods completely. Therefore, in order to show the robustness and stability of the proposed methods, multiple independent runs are needed. In some examples, 100 independent runs are carried out to investigate the statistical performance of the proposed methods. The direct Monte Carlo simulation (MCS), subset simulation (Subset), and the importance sampling method based on the kernel sampling density (IS-KSD) with fixed width [4] are also applied to these examples for comparison. For each method with 100 runs, a number of 100 failure probability estimates, Pˆ F ,i (i = 1, . . . , 100), can be obtained. Based on these estimates, the corresponding statistical values are computed, including the mean values and c.o.v.sof the 100 1 estimates, and the mean values of relative error ε¯ = 100 i =1 εi (where

εi = | Pˆ F ,i − P F |/ P F ).

6.1. Example 1: Series system of nonlinear limit state functions This example considers a series system with nonlinear limit state functions given by [4]:

g (x) = min( g 1 , g 2 )

257



g 1 (x) = c − 1 − x2 + exp − g 2 ( x) =

c2 2

− x1 x2

x21 10



 +

x1

4

5 (10)

where x = [x1 , x2 ], x1 and x2 are independent, identically and standard normally distributed. This system has two design points ∗(1) = which have the same√distance √ from the origin, given by x [0, c ] and x∗(2) = [c / 2, c / 2]. Both of them are important in contributing to the failure probability. Herein, three cases, (a), (b) and (c) are considered corresponding to c = 3, 4 and 5, respectively. The proposed methods are applied to this example which has multiple important regions to show their robustness and efficiency. The ‘exact’ failure probabilities taken from [4] for the cases c = 3, 4, 5 are 2.53 × 10−3 , 6.81 × 10−5 and 6.97 × 10−7 , respectively. In order to show the performance of the proposed method visually, a total of N T = M × N = 100 × 5 samples are used to complete the analysis for each case and the samples for cases (a), (b) and (c) simulated in MCIS method are shown in Fig. 2(a)–(c) respectively. In Fig. 2, the circles denote the Markov chain time steps which are also the sampling centers of importance sampling, and the dots denote the importance samples. It can be seen that the Markov chain time steps and the simulated samples well populate the two important regions which are relatively independent of the value of c and hence the probability level of the problem. The estimates obtained by MCIS for cases (a), (b) and (c) with 500 samples are 2.6044 × 10−3 , 6.6131 × 10−5 and 6.5357 × 10−7 , respectively. Similar results can also be obtained by the MCIS-QMC method, thus it can be seen that the proposed methods can efficiently handle the problem of multiple failure modes and also the problem of small failure probability. Compared with direct MCS method, they obviously own higher efficiency while dealing with small failure probability problems, i.e. case (c), as the MCS method need at least 1/ P F samples in order to evaluate the failure probability P F . Compared with traditional importance sampling method based on design point(s), the proposed methods do not need the solution of design point(s) which is sometimes extremely complicated. Compared with importance sampling method based on kernel density, the proposed methods do not need the extra presamples in the failure region and the computational effort in the estimation of the kernel density. Table 1 shows the statistical results obtained by different methods in 100 independent runs. In each run of the proposed methods (MCIS and MCIS-QMC) M = 100 and N = 10 are adopted, thus a total of N T = M × N = 100 × 10 = 1000 samples are used. In the ISKSD method, M = 100 samples are used to construct the KSD, and then 1000 importance samples are generated from the constructed KSD to complete the analysis. Meanwhile, in order to make thorough comparison, a different setting is also used for the IS-KSD method, say, M = 500 pre-samples and 500 importance samples are adopted. The direct MCS method and subset simulation are also applied where MCS uses 104 , 105 and 107 samples for these cases, respectively and the subset simulation uses 1000 samples for each intermediate level and conditional probability is chosen to be 0.1. The statistical results obtained from 100 runs, i.e., the mean values of the failure probability estimates, the c.o.v.s and mean values of relative errors, are listed in Table 1. Note that smaller value of c.o.v. implies faster convergence rate of the method. And smaller value of relative error implies higher accuracy of the method. It can be seen that the mean values of c.o.v.s and relative errors of the proposed methods (MCIS and MCIS-QMC) for all cases are all smaller than those of the other methods. It can be drawn that the proposed methods performance better than the IS-KSD and owns

258

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

(a)

(b)

(c) Fig. 2. The samples of MCIS method for cases (a), (b) and (c) in example 1 (the circles denote the Markov chain states and the dots denote the importance samples).

high efficiency and accuracy. Furthermore, it can also be seen that the MCIS-QMC has the smallest mean values of c.o.v.s and relative errors, thus the MCIS-QMC method performances better than MCIS. This indicates that the low-discrepancy sequence has further improved the MCIS method. For the three cases of different probability levels (10−3 , 10−5 and 10−7 ) in this example, the number of samples used by the proposed methods and the values of c.o.v.s of estimates do not vary much, it seems that the computational cost needed by the proposed methods is nearly independent of the probability level. Thus the proposed methods own extremely high efficiency in handling problems with small failure probabilities. 6.2. Example 2: Parallel system of linear limit state functions This five-dimensional example is studied by Au and Beck [4], and it represents a parallel system with limit state functions given by:

g (x) = max( g 1 , g 2 , g 3 , g 4 ) g i (x) = 2.5 + 0.25 cos(π i /5) − xi − xi +1 ,

i = 1, . . . , 4

(11)

where xi , i = 1, . . . , 5 are independent, identically and standard normally distributed. The exact probability of failure was computed to be P F = 0.0001779 by direct MCS with 107 samples. Table 2 shows the statistical results based on 100 independent runs. The direct MCS method with 105 samples and subset simulation with 1000 × 4 samples are applied. The proposed methods, MCIS and MCIS-QMC, use a total of N T = M × N = 200 × 5 = 1000 samples to compute the failure probability for each run. The ISKSD method firstly uses M = 200 samples to construct the ISD, and then uses 1000 importance samples to calculate the failure probability. In order to compare them thoroughly, a different setting for the IS-KSD method is also selected, i.e., M = 500 samples to construct the ISD and 500 importance samples to carry out the analysis, which use the same computational cost with the proposed methods. The direct MCS method is applied with 105 samples. The subset simulation with 1000 samples for each intermediate level is also applied. Good convergence is achieved for the proposed methods with the c.o.v.s below 0.2 which are smaller than those of IS-KSD and subset. Except the direct MCS method, the MCIS-QMC obtains the estimates of failure probabilities which have the smallest mean values of the c.o.v.s and relative errors. The superiority of the proposed methods in accuracy and efficiency is demonstrated.

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

259

Table 1 Statistical performance based on 100 runs for example 1. Number of samples

Mean of P F

C.o.v. of P F

Mean of relative error

Case 1: C = 3 Exacta [4] MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC

104 1000 × 3 100 + 1000 500 + 500 100 × 10 100 × 10

2.53 × 10−3 3.4100 × 10−3 3.4962 × 10−3 2.4388 × 10−3 2.5578 × 10−3 2.5011 × 10−3 2.6192 × 10−3

1.5698 × 10−1 2.8840 × 10−1 2.5633 × 10−1 1.3192 × 10−1 6.9482 × 10−2 5.0665 × 10−2

3.5771 × 10−1 4.3036 × 10−1 1.2921 × 10−1 6.7024 × 10−1 5.6478 × 10−2 5.0940 × 10−2

Case 2: C = 4 Exact MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC

105 1000 × 5 100 + 1000 500 + 500 100 × 10 100 × 10

6.81 × 10−5 8.9100 × 10−5 9.1150 × 10−5 6.2665 × 10−5 6.7058 × 10−5 6.6656 × 10−5 6.8555 × 10−5

3.1018 × 10−1 8.4562 × 10−1 4.5838 × 10−1 1.1970 × 10−1 8.7755 × 10−2 5.9962 × 10−2

3.9944 × 10−1 7.6545 × 10−1 2.3768 × 10−1 7.9225 × 10−2 6.7771 × 10−2 2.3355 × 10−2

Case 3: C = 5 Exact MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC

107 1000 × 7 100 + 1000 500 + 500 100 × 10 100 × 10

6.97 × 10−7 9.2300 × 10−7 9.5468 × 10−7 2.7508 × 10−7 6.4073 × 10−7 6.7596 × 10−7 7.0204 × 10−7

3.2259 × 10−1 2.1692 3.9501 × 10−1 2.1537 × 10−1 1.5531 × 10−1 8.5524 × 10−2

4.4017 × 10−1 1.0188 6.1403 × 10−1 1.4359 × 10−1 1.0406 × 10−1 6.5849 × 10−2

Number of samples

Mean of P F

C.o.v. of P F

Mean of relative error

107 105 1000 × 4 200 + 1000 500 + 500 200 × 5 200 × 5

1.779 × 10−4 1.5720 × 10−4 1.7537 × 10−4 1.0807 × 10−4 1.6883 × 10−4 1.6022 × 10−4 1.8197 × 10−4

1.2661 × 10−1 3.1562 × 10−1 6.3743 × 10−1 2.4167 × 10−1 1.7284 × 10−1 1.5823 × 10−1

1.2651 × 10−1 2.5299 × 10−1 4.8019 × 10−1 1.6670 × 10−1 1.5041 × 10−1 1.2490 × 10−1

a

The ‘exact’ failure probability is taken from [4].

Table 2 Statistical performance based on 100 runs for example 2.

Exacta MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC a

The ‘exact’ probability of failure was computed by MCS with 107 samples.

Table 3 Statistical performance based on 100 runs for example 3.

Exact [17] MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC

Number of samples

Mean of P F

C.o.v. of P F

Mean of relative error

104 1000 × 2 200 + 1000 500 + 500 200 × 5 200 × 5

1.23 × 10−2 1.2180 × 10−2 1.2483 × 10−2 1.1433 × 10−2 1.1816 × 10−2 1.1069 × 10−2 1.1392 × 10−2

9.5204 × 10−2 1.8042 × 10−1 3.0938 × 10−1 2.6720 × 10−1 1.7532 × 10−1 8.0830 × 10−2

7.5935 × 10−2 1.4855 × 10−1 2.2200 × 10−1 1.9168 × 10−1 1.3224 × 10−1 8.9857 × 10−2

6.3. Example 3: A noisy boundary case with non-normal variables This example is to test the efficiency of the proposed method with respect to noisy failure boundaries [17]. The limit state function is

g (x) = x1 + 2x2 + 2x3 + 2x4 − 5x5 − 5x6

+ 0.001

6 

cos(100xi )

(12)

i =1

where all the variables xi (i = 1, . . . , 6) are lognormally distributed, the mean values and standard deviations of x1 − x4 are 120 and 12, respectively, the mean and standard deviation of x5 are 50 and 15 respectively, and the mean and standard deviation of x6 are 40

and 12 respectively. The exact failure probability of this example is 0.0123 [17]. Table 3 shows the statistical performance based on 100 independent runs. The proposed methods complete the analysis with N T = M × N = 200 × 5 = 1000 samples. Two different settings are used for the IS-KSD, i.e., 200 pre-samples plus 1000 importance samples and 500 pre-samples plus 500 importance samples, respectively. The direct MCS method is applied with 105 samples. The subset simulation with 1000 samples for each intermediate level is also applied. It can be seen that the proposed methods achieve the smaller c.o.v.s and relative errors of the failure probability estimates than those of subset and IS-KSD methods. The applicability of the proposed methods for noisy failure boundaries problem are demonstrated.

260

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

Table 4 Statistical performance based on 100 runs for example 4.

Exact [28] MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC

Number of samples

Mean of P F

C.o.v. of P F

Mean of relative error

105 1000 × 3 200 + 2000 1000 + 1000 200 × 10 200 × 10

4.85 × 10−3 4.8913 × 10−3 5.0331 × 10−3 3.1221 × 10−3 4.6076 × 10−3 4.2260 × 10−3 4.5314 × 10−3

4.4582 × 10−2 2.3153 × 10−1 3.6840 × 10−1 1.8662 × 10−1 8.9911 × 10−2 8.7978 × 10−2

3.7423 × 10−2 1.8769 × 10−1 3.9173 × 10−1 1.5144 × 10−1 1.3487 × 10−1 8.5637 × 10−2

Fig. 3. The finite element model of the inside flap of a aircraft.

6.4. Example 4: Series system of non-normal variables

Table 5 The distribution information of the random variables for example 5.

A plane frame structure which is a series system from [28] is considered in this example, based on the principle of virtual work three limit state functions are given as

Beam:

g beam = x2 + 2x3 + x4 − hx7

Sway:

g sway = x1 + x2 + x4 + x5 − hx6

Frame:

g frame = x1 + 2x3 + 2x4 + x5 − hx6 − hx7

(13)

where plastic moment capacities x1 –x5 and loads x6 , x7 are lognormally distributed independent variables, and the mean values and standard deviations of x1 –x5 are 134.9 and 13.49, respectively, x6 has mean 50 and standard deviation 15, and x7 has mean 40 and standard deviation 12. The constant parameter h = 5. The exact failure probability of this example is 0.00485 [28]. For this series system with non-normal variables, the proposed methods use N T = M × N = 200 × 10 = 2000 samples to calculate the failure probability, and the statistical results based on 100 independent runs are listed in Table 4. Two settings are used for the IS-KSD method: 200 pre-samples plus 2000 importance samples and 1000 pre-samples plus 1000 importance samples. The direct MCS method is applied with 105 samples. The subset simulation with 1000 samples for each intermediate level is also applied. Both of the proposed methods achieve good convergences with the mean values of the c.o.v.s of the failure probability estimates below 0.1. From the smallest mean values (except those of the direct MCS method) of c.o.v.s and relative errors of the estimates listed in Table 4, it also can be seen that the proposed MCIS-QMC obtains the most stable and accurate estimate. 6.5. Example 5: Reliability of a aircraft flap An aircraft inside flap is considered here and its finite element model is shown in Fig. 3. The structure suffers from aerodynamic load which is transformed into concentrated load F applied to the nodes of the finite element model. Taking the most dangerous working condition into consideration, failure is defined as the maximum of the strain of all the nodes exceeding an admissible maximal strain. The limit state function is then given by:

g (x) = D − f (t 1 , t 2 , t 3 , E 1 , G 1 , E 2 , G 2 , F )

(14)

Variable

Mean

C.o.v.

Standard deviation

Distribution

t 1 (cm) t 2 (cm) t 3 (cm) E 1 (MPa) G 1 (MPa) E 2 (MPa) G 2 (MPa) F

2 2 4 70 380 26458.6 72 450 27236.8 0

0.02 0.02 0.02 0.02 0.02 0.02 0.02 –

0.04 0.04 0.08 1407.6 529.172 1449 544.736 0.2

Normal Normal Normal Normal Normal Normal Normal Normal

Table 6 The results by different methods for example 5.

Direct MCS Subset IS-KSD IS-KSD MCIS MCIS-QMC a

Number of samples

PF

Cov( P F )

Relative error

104 1000 × 3 200 + 1000 500 + 500 200 × 5 200 × 5

3.5 × 10−3 3.3603 × 10−3 2.3185 × 10−3 4.9581 × 10−3 3.3110 × 10−3 3.2965 × 10−3

0.1687 0.3665 0.8780 0.4743 0.0498a 0.0502a

– 3.99% 33.76% 41.66% 5.40% 5.81%

The c.o.v.s are approximated by Eq. (10).

where t 1 , t 2 , t 3 are the thicknesses of three kinds of beams of the flap, E 1 ( E 2 ) and G 1 (G 2 ) are the elastic modulus and shear modulus, respectively. As an instrumental random variable F represents the randomness of the load applied to the nodes. For the load applied in the node i, the value is F i = (1 + F ) F i0 where F i0 is a constant nominal value. All variables are mutual independent normal variables whose distribution information is given in Table 5. The failure probability results calculated by different methods are given in Table 6. The direct MCS method is used with 104 samples to obtain an estimate which is seen as the ‘exact’ value. In order to compare with the proposed methods, the subset simulation and IS-KSD are also applied. As seen from Table 6, the IS-KSD uses 200 pre-samples and 1000 importance samples to analyze and obtains an estimate with c.o.v. bigger than 0.8. Another trial of IS-KSD with 500 pre-samples and 500 importance samples obtains an estimate with c.o.v. still bigger than 0.4. Subset simulation uses 1000 samples for each level and a total of 3000 samples are used to obtain an estimate with c.o.v. bigger than 0.3. The proposed methods (MCIS and MCIS-QMC) obtain the results with only N T = M × N = 200 × 5 = 1000 samples which use the

X. Yuan et al. / Aerospace Science and Technology 29 (2013) 253–261

lowest computational cost among these methods. The corresponding c.o.v.s of estimates are also approximately computed by using Eq. (9) and the values are both below 0.1. The relative errors by different methods are also listed in the table. It is seen that the IS-KSD has bigger relative error (exceed 30%) for both settings and the proposed methods achieve accurate results with relative errors below 10%. The advantages of the proposed method in efficiency and precision are shown. 7. Conclusions A novel adaptive importance sampling algorithm has been proposed to compute the failure probability in structural reliability analysis. It utilizes the Markov chain Metropolis algorithm to establish a novel adaptive multi-point importance sampling procedure. The proposal distributions are taken as not only the distribution for candidate state generation but also the ISD for importance sampling, that is, not one but a set of samples is generated at a time by the proposal distribution, and then the next Markov chain time step, which is also the next sampling center of importance sampling, is properly chosen from these samples. As the number of Markov chain steps increases, Markov chain can adaptively explore the important region which significantly contributes to the failure probability, thus the importance sampling based on the Markov chain time steps leads to an efficient and accurate estimate of the failure probability. Compared with traditional importance sampling based on design point(s), it does not need the solution of design point(s). Compared with existent adaptive sampling schemes, it does not need the pre-sampling process and the construction of the ISD. In addition, low-discrepancy sequence is adopted in importance sampling in order to further improve the accuracy and robustness of the method. Numerical examples have been given to demonstrate the efficiency, robustness and accuracy of the proposed methods. The results show that the proposed method is relatively independent of the probability level and the geometry of failure region of the problem. Both MCIS-QMC and MCIS method can obtain desirable the failure probability estimate efficiently, i.e., they can achieve a c.o.v. of about 0.1–0.2 of the failure probability estimate over several orders of magnitude with about 1000–2000 samples. Furthermore, the MCIS-QMC method with low-discrepancy sampling exhibits higher accuracy and robustness than MCIS. Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant Nos. 51105309 and 51175425) and the Aviation Science Foundation (Grant No. 2011ZA53015). These supports are gratefully acknowledged. References [1] G.L. Ang, A.H.-S. Ang, W.H. Tang, Optimal importance sampling density estimator, Journal of Engineering Mechanics (ASCE) 118 (1992) 1146–1163.

261

[2] W.O. Aszkenazy, Hammersley/Halton sequence generation, http://www.math. rwth-achen.de/mapleAnswers/html/456.html, 2002. [3] S.K. Au, Probabilistic failure analysis by importance sampling Markov chain simulation, Journal of Engineering Mechanics (ASCE) 130 (2004) 303–311. [4] S.K. Au, J.L. Beck, A new adaptive important sampling scheme, Structural Safety 21 (1999) 135–158. [5] S.K. Au, J.L. Beck, First excursion probabilities for linear systems by very efficient importance sampling, Probabilistic Engineering Mechanics 16 (2001) 193–207. [6] S.K. Au, J.L. Beck, Estimation of small failure probability in high dimensions by subset simulation, Probabilistic Engineering Mechanics 16 (2001) 263–277. [7] S.K. Au, J. Ching, J.L. Beck, Application of subset simulation method to reliability benchmark problems, Structural Safety 29 (3) (2007) 183–193. [8] S.K. Au, C. Papadimitriou, J.L. Beck, Reliability of uncertain dynamical systems with multiple design points, Structural Safety 21 (2) (1999) 113–133. [9] J. Bosserta, M. Feindta, U. Kerzelb, Fast integration using quasi-random numbers, Nuclear Instruments and Methods A 559 (2006) 232–236. [10] Z. Botev, D.P. Kroese, An efficient algorithm for rare-event probability estimation, combinatorial optimization, and counting, Methodology and Computing in Applied Probability 10 (2008) 471–505. [11] Z. Botev, D.P. Kroese, T. Taimre, Generalized cross-entropy methods with applications to rare-event simulation and optimization, Simulation 83 (2007) 785– 806. [12] C.G. Bucher, Adaptive sampling – an iterative fast Monte Carlo procedure, Structural Safety 5 (1988) 119–126. [13] M. Caserta, M. Cabo Nodar, A cross entropy based algorithm for reliability problems, Journal of Heuristics 15 (2009) 479–501. [14] H. Dai, W. Wang, Application of low-discrepancy sampling method in structural reliability analysis, Structural Safety 31 (2009) 55–64. [15] B. Echard, N. Gayton, M. Lemaire, N. Relun, A combined importance sampling and Kriging reliability method for small failure probabilities with timedemanding numerical models, Reliability Engineering and System Safety 111 (2013) 232–240. [16] K.T. Fang, Y. Wang, Number-Theoretic Methods in Statistics, Chapman & Hall, London, 1994. [17] P.L. Liu, A. Der Kiureghian, Optimization algorithms for structural reliability, Structural Safety 9 (1991) 161–177. [18] R.E. Melchers, Importance sampling in structural systems, Structural Safety 6 (1989) 3–10. [19] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, Equations of state calculations by fast computing machines, Journal of Chemical Physics 21 (1953) 1087–1092. [20] J. Morio, Non-parametric adaptive importance sampling for the probability estimation of a launcher impact position, Reliability Engineering and System Safety 96 (1) (2011) 178–183. [21] J. Morio, R. Pastel, F. Le Gland, Missile target accuracy estimation with importance splitting, Aerospace Science and Technology 25 (2013) 40–44. [22] J. Nie, B.R. Ellingwood, A new directional simulation method for system reliability. Part I: Application of deterministic point sets, Probabilistic Engineering Mechanics 19 (2004) 425–436. [23] H.J. Pradlwarter, M.F. Pellissetti, C.A. Schenk, G.I. Schuëller, A. Kreis, S. Fransen, A. Calvi, M. Klein, Realistic and efficient reliability estimation for aerospace structures, Computer Methods in Applied Mechanics and Engineering 194 (2005) 1597–1617. [24] H.J. Pradlwarter, G.I. Schuëller, P.S. Koutsourelakis, D.C. Charmpis, Application of line sampling simulation method to reliability benchmark problems, Structural Safety 29 (2007) 208–221. [25] R.Y. Rubinstein, Simulation and the Monte Carlo Method, Wiley, 2007. [26] I.M. Sobol, On quasi-Monte Carlo integrations, Mathematics and Computers in Simulation 47 (1998) 103–112. [27] X. Wang, F.J. Hichernell, Randomized Halton sequences, Mathematical and Computer Modelling 32 (2000) 887–899. [28] K. Ziha, Descriptive sampling in structural safety, Structural Safety 17 (1995) 33–41.