Review on Energy Resilience

Probabilistic Engineering Mechanics 34 (2013) 1–11

Contents lists available at ScienceDirect

Probabilistic Engineering Mechanics

0266-89 http://d

n Tel.: E-m

journal homepage: www.elsevier.com/locate/probengmech

Local estimation of failure probability function by weighted approach

Xiukai Yuan n

Department of Aeronautics, Xiamen University, 422 South Siming Road, Xiamen, Fujian 361005, PR China

a r t i c l e i n f o

Article history: Received 21 April 2012 Received in revised form 29 April 2013 Accepted 8 May 2013 Available online 17 May 2013

Keywords: Reliability Failure probability Importance sampling Subset Simulation Monte Carlo simulation

20/$ - see front matter & 2013 Elsevier Ltd. A x.doi.org/10.1016/j.probengmech.2013.05.001

+86 0592 2180699. ail addresses: [email protected], erudite

a b s t r a c t

In the reliability-based design of engineering systems, it is often required to evaluate the failure probability for different values of distribution parameters involved in the specification of design configuration. The failure probability as a function of the distribution parameters is referred as the ‘failure probability function (FPF)’ in this work. From first principles, this problem requires repeated reliability analyses to estimate the failure probability for different distribution parameter values, which is a computationally expensive task. A “weighted approach” is proposed in this work to locally evaluate the FPF efficiently by means of a single simulation. The basic idea is to rewrite the failure probability estimate for a given set of random samples in simulation as a function of the distribution parameters. It is shown that the FPF can be written as a weighted sum of sample values. The latter must be evaluated by system analysis (the most time-consuming task) but they do not depend on the distribution. Direct Monte Carlo simulation, importance sampling and Subset Simulation are incorporated under the proposed approach. Examples are given to illustrate their application.

& 2013 Elsevier Ltd. All rights reserved.

1. Introduction

Reliability-based design optimization (RBDO) [1–3] aims at seeking the best compromise between cost and safety while considering system uncertainties. In solving a RBDO problem, the failure probabilities of the target system under various design configurations are usually needed to be evaluated. The failure probability function, PFðθÞ, can be defined as the failure probability with regard to the set of design parameters θ ¼ ½θ1; θ2; :::; θnθ �(nθ is their number). If the whole function PFðθÞ can be obtained before- hand, the RBDO problem can be transformed into an ordinary one as a decoupled problem [4–8]. Another application of FPF is “reliability sensitivity analysis”, defined as the derivative of PFðθÞ with respect to θ. This derivative is also the main target where gradient-based methods are used for optimization.

Obtaining the FPF as a function is often a computationally challenging task in real applications. From first principles, it requires repeated reliability analyses over various design para- meter values in order to obtain the point-wise values of the FPF. It is difficult to obtain the failure probability as an explicit function with regard to the design parameters. One strategy is to construct an approximation by selecting some predefined interpolation points in the space of the design parameters. For example, Gasser [4] used a quadratic function of θ to locally approximate log½PFðθÞ�, The number of coefficients to be determined in the linear function

ll rights reserved.

[email protected]

is equal to nθ þ nθðnθ þ 1Þ=2. At least nθ þ nθðnθ þ 1Þ=2 reliability analyses are needed to obtain all the coefficients. Jensen [5] adopted a linear function when dealing with a linear system subject to stochastic excitation, for which at least nθ reliability analyses are needed.

Au [9] proposed an augmenting idea that allowed the informa- tion of the FPF to be obtained in a single simulation run. The design parameter θ is treated as uncertain with a predefined “instrumental” distribution and Bayes′ rule is used to transform the PFðθÞ function into an expression including the probability density function (PDF) of θ conditioning on the failure event F, denoted by f ðθjFÞ; the PDF is estimated by histograms. Ching [6,7] adopted a pre-defined form of the FPF and applied the maximum entropy principle to estimate the coefficients. Expositions in reliability- based design in geotechnical engineering can also be found in [10,11]. Based on the augmenting idea, Taflanidis proposed a framework for robust stochastic system design and reliability optimization [12,13].

This paper focuses on the class of RBDO where the design parameters correspond to the distribution parameters of the random variables in the problem [8]. Such problems are encoun- tered when, for example, the mean of the geometrical size of a structure such as thickness, length, height, etc are the design parameters in RBDO. Under this limited scope, advantage could be taken by analyzing the mathematical structure of the problem. A ‘weighted approach’ is proposed that allows the FPF to be obtained through a single traditional reliability analysis, thereby avoiding repeated reliability analyses. It is shown that the FPF can be expressed as a weighted sum of sample values. The influence of

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–112

the design parameters is encapsulated in the weights. The sample values must be evaluated through system analyses but they do not depend on the design parameters.

This paper is outlined as follows. The problem context is presented first. The weighted approach is then developed based on three different reliability analysis methods. Examples are given to illustrate the proposed method.

2. Problem context

Let x ¼ ½x1; x2; :::; xn� be the set of the uncertain variables of a system (n is their number) and θ ¼ ½θ1; θ2; :::; θnθ � be the set of design parameters. In this work, it is assumed that θ corresponds to the parameters that specify the probability distribution of x, e.g. the mean or standard deviation. Failure of the system is defined as F ¼ x : gðxÞ≤0

� � where gðxÞ is the limit state function. For a given

design θ, the failure probability PFðθÞ is formulated as

PFðθÞ ¼ PðFjθÞ ¼ Z

IFðxÞf ðxjθÞdx ð1Þ

where the integral is over Rn; IFðxÞ is the indicator function of the failure region F: if x∈F, IFðxÞ ¼ 1, and IFðxÞ ¼ 0 otherwise; f ðxjθÞ is the joint probability density function (PDF) of x conditional on θ. For design sensitivity purposes, it is often required to evaluate PFðθÞ for a number of values of θ within some admissible regions. The quantity PFðθÞ as a function of θ is referred to as the ‘failure probability function’ (FPF), which is the main target in this work.

3. Weighted approach based on simulation

The basic idea behind the proposed approach is that, when the design parameters correspond to the distribution parameters of the random variables in the problem, they only affect the prob- ability term in the failure probability integral and are therefore ‘uncoupled’ from the indicator function that involves system analysis. As a result, it is possible to make use of the samples generated from a ‘nominal value’ of the design parameters to infer the failure probability corresponding to other values of the design parameters. The actual mechanism by which this is done depends on the simulation method adopted, however. In what follows, the development shall be focused on three conventional simulation methods, namely, direct Monte Carlo simulation, importance sampling simulation and Subset Simulation.

3.1. Weighted approach based on direct Monte Carlo simulation

Consider evaluating PFðθ0Þ for a given ‘nominal’ design θ0 in direct Monte Carlo simulation. The failure probability is viewed as an expectation where x is distributed as the original PDF f ðxjθ0Þ, i.e.,

PFðθ0Þ ¼ Z

IFðxÞf ðxjθ0Þdx ¼ Eθ0 ½IFðxÞ� ð2Þ

Here, Eθ0 ½⋅� denotes the expectation under the PDF f ðxjθ0Þ. Using f ðxjθ0Þ to generate a set of samples, xðjÞ (j ¼ 1; 2; :::; N, N is the number of samples), an estimate of PFðθ0Þ is given by:

P̂Fðθ0Þ ¼ 1 N

∑ N

j ¼ 1 IFðxðjÞÞ ð3Þ

For general θ, the corresponding PDF is f ðxjθÞ. Instead of generat- ing another set of samples from f ðxjθÞ, the FPF value PFðθÞ can be estimated using samples generated from f ðxjθ0Þ via a change of

probability distribution

PFðθÞ ¼ Z

IFðxÞf ðxjθÞdx ¼ Z

IFðxÞ f ðxjθÞ f ðxjθ0Þ

f ðx θ0Þdx ��

¼ Eθ0 ½IFðxÞrMCSðx; θ; θ0Þ� ð4Þ where

rMCSðx; θ; θ0Þ ¼ f ðxjθÞ f ðxjθ0Þ

ð5Þ

is a weighting factor, being the ratio of two PDFs. It can be seen that the FPF is expressed as the expectation of the weighting factor under the PDF f ðxjθ0Þ. Thus, through a single Monte Carlo simula- tion analysis with f ðxjθ0Þ, the FPF PFðθÞ can be obtained. Here f ðxjθ0Þ is referred as an “instrumental” PDF because it is merely chosen to generate random samples for inferring FPF values for different design parameter values.

Based on the same set of samples generated by f ðxjθ0Þ, the estimate of PFðθÞ is given by

P̂FðθÞ ¼ 1 N

∑ N

j ¼ 1 IFðxðjÞÞ

f ðxðjÞjθÞ f ðxðjÞjθ0Þ

ð6Þ

It should be noted that Eqs. (6) and (3) share the same set of samples xðjÞ (j ¼ 1; 2; :::; N) and indicator function values IFðxðjÞÞ. This means that for different values of θ, the values of the indicator function IFðxðjÞÞ are the same and do not need to be evaluated repeatedly in order to estimate PFðθÞ. Only the weighting factor needs to be re-evaluated, but the computational cost is minimal as it does not involve system analyses.

For proper choice of the instrumental PDF (see later), standard analysis indicates that the function estimate P̂FðθÞ is unbiased. Its variance and c.o.v. (coefficient of variation¼standard deviation/ mean) can be estimated by

Var P̂FðθÞ h i

≈ 1

N−1 1 N

∑ N

j ¼ 1 IFðx jð ÞÞ

f ðx jð ÞjθÞ f ðx jð Þjθ0Þ

� �2 −½P̂FðθÞ�2

( ) ð7Þ

Cov½P̂FðθÞ� ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Var½P̂FðθÞ�

q =E½P̂FðθÞ� ð8Þ

3.2. Weighted approach based on importance sampling

In importance sampling method [14–17], the failure probability under f ðxjθ0Þ in Eq. (2) is rewritten as

PFðθ0Þ ¼ Z

IFðxÞ f ðxjθ0Þ HðxÞ HðxÞdx ¼ EH IFðxÞ

f ðxjθ0Þ HðxÞ

� � ð9Þ

where HðxÞ is the importance sampling function; EH½⋅� denotes the expectation under the PDF HðxÞ. Using HðxÞ to generate a set of samples, xðjÞ (j ¼ 1; 2; :::; N), an unbiased estimate of PFðθ0Þ is given by

P̂Fðθ0Þ ¼ 1 N

∑ N

j ¼ 1 IFðxðjÞÞ

f ðxðjÞjθ0Þ HðxðjÞÞ ð10Þ

For general θ, the FPF value PFðθÞ can be expressed as an expecta- tion under HðxÞ

PFðθÞ ¼ Z

IFðxÞ f ðxjθÞ HðxÞ HðxÞdx ¼ EH IFðxÞ

f ðxjθÞ HðxÞ

� � ¼ EH IFðxÞrISðx; θÞ½ �

ð11Þ where

rISðx; θÞ ¼ f ðxjθÞ HðxÞ ð12Þ

is the weighting factor. Here HðxÞ is naturally the instrumental PDF. Through a single importance sampling analysis with HðxÞ, the FPF PFðθÞ can be obtained.

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–11 3

Based on the same set of samples generated from HðxÞ, xðjÞ (j ¼ 1; 2; :::; N), PFðθÞ can be estimated by virtue of Eq. (11)

P̂FðθÞ ¼ 1 N

∑ N

j ¼ 1 IFðxðjÞÞ

f ðxðjÞjθÞ HðxðjÞÞ ð13Þ

Again, Eqs. (10) and (13) share the same set of samples xðjÞ

(j ¼ 1; 2; :::; N) and indicator function IFðxðjÞÞ values. These quanti- ties need not be evaluated repeatedly. The FPF is explicitly expressed as a function of the design parameters θ based on the samples generated from one traditional importance sampling procedure.

For proper choice of the instrumental PDF, the estimate P̂FðθÞ is unbiased. Its variance and c.o.v. can be estimated by

Var P̂FðθÞ h i

≈ 1

N−1 1 N

∑ N

j ¼ 1 IFðx jð ÞÞ

f ðx jð ÞjθÞ Hðx jð ÞÞ

� �2 −½P̂FðθÞ�2

( ) ð14Þ

Cov½P̂FðθÞ� ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Var½P̂FðθÞ�

q =E½P̂FðθÞ� ð15Þ

The weighted Monte Carlo simulation can be considered as a special case of the weighted importance sampling when HðxÞ ¼ f ðxjθ0Þ.

3.3. Weighted approach based on Subset Simulation

Subset Simulation [18,19] is an adaptive procedure based on Markov chain Monte Carlo for investigating rare events. For the original PDF f ðxjθ0Þ the failure probability is given by

PFðθ0Þ ¼ PðFmÞ ¼ PðF1Þ ∏ m−1

i ¼ 1 PðFiþ1jFiÞ ð16Þ

where F1⊃F2⊃:::⊃Fm ¼ F is a sequence of failure events, Fi ¼ gðxÞ4bi

� � (i ¼ 1; 2; :::; m); b1; b2; :::; bm−1 are some intermediate

threshold values which are adaptively chosen so that the corre- sponding sample estimates for the probabilities PðF1Þ, PðF2jF1Þ,…, PðFm−1jFm−2Þ are all equal to p0, e.g., p0 ¼ 0:1 for convenience which is used in this contribution. The conditional distribution of the samples in the k-th (k ¼ 1; 2; :::; m) level is denoted by f ½xjðFk; θ0Þ� ¼ IFk ðxÞf ðxjθ0Þ=PðFkÞ. The failure probability PFðθ0Þ in Eq. (16) can be rewritten as

PFðθ0Þ≈pm−10 PðFmjFm−1Þ ¼ pm−10 ⋅EFm−1;θ0 ½IFm ðxÞ� ð17Þ

where EFm−1;θ0 ½⋅� denotes the expectation under the PDF f ½xjðFm−1; θ0Þ� ¼ IFm−1 ðxÞf ðxjθ0Þ=PðFm−1Þ.

Suppose Ni samples are generated in the i-th (i ¼ 1; 2; :::; m−2) level and N samples are generated in the (m−1)th level, i.e., xðjÞ

(j ¼ 1; 2; :::; N). The reason for a non-uniform number of samples at different levels will be explained later in Section 4. The failure probability PFðθ0Þ can be estimated as

P̂Fðθ0Þ ¼ pm−10 PðFmjFm−1Þ ¼ pm−10 ⋅ 1 N

∑ N

j ¼ 1 IFm ðxðjÞÞ ð18Þ

For general θ, the corresponding PDF is f ðxjθÞ. The FPF PFðθÞ can be obtained using the samples generated from f ðxjθ0Þ. This can be reasoned as follows. For the first level F1, the failure probability of F1 under f ðxjθÞ can be expressed by

PF1 ðθÞ ¼ Z

IF1 ðxÞf ðx; θÞdx ¼ Z

IF1 ðxÞ f ðxjθÞ f ðxjθ0Þ

f ðx θ0Þdx ��

¼ Eθ0 IF1 ðxÞ f ðxjθÞ f ðxjθ0Þ

� � ð19Þ

Similarly, for the k-th (k ¼ 1; 2; :::; m) level Fk, the conditional probability PðFkjFk−1Þ under f ðxjθÞ, denoted by PFkjFk−1 ðθÞ, can be

expressed as

PFkjFk−1 ðθÞ ¼ Z

IFk ðxÞ f ½xjðFk−1; θÞ� f ½xjðFk−1; θ0Þ�

f x ðFk−1; θ0Þ �� �dx�

¼ EFk−1;θ0 IFk ðxÞ IFk−1 ðxÞf ðxjθÞ IFk−1 ðxÞf ðxjθ0Þ

⋅ PFk−1;θ0 PFk−1 ðθÞ

� �

¼ PFk−1;θ0 PFk−1 ðθÞ

⋅EFk−1;θ0 IFk ðxÞ f ðxjθÞ f ðxjθ0Þ

� � ð20Þ

where PFk−1 ðθ0Þ ¼ PðF1ÞPðF2jF1Þ⋯PðFk−1jFk−2Þ≈pk−10 is the failure probability of Fk−1 under f ðxjθ0Þ. Setting k ¼ m in Eq. (20) yields

PFmjFm−1 ðθÞPFm−1 ðθÞ ¼ PFm−1;θ0 ⋅EFm−1;θ0 IFm ðxÞ f ðxjθÞ f ðxjθ0Þ

� � ð21Þ

That is

PFðθÞ ¼ PFm ðθÞ ¼ PFmjFm−1 ðθÞPFm−1 ðθÞ ¼ PFm−1;θ0 ⋅EFm−1;θ0 ½IFm ðxÞrSSðx; θ; θ0Þ� ð22Þ

where

rSSðx; θ; θ0Þ ¼ f ðxjθÞ f ðxjθ0Þ

ð23Þ

is the weighting factor. It can be seen that the FPF is expressed as the expectation of the weighting factor under the PDF f ðxjθ0Þ. Thus the FPF PFðθÞ can be obtained through a single Subset Simulation analysis under f ðxjθ0Þ. Here f ðxjθ0Þ is the instrumental PDF.

Using the same set of samples generated in the (m−1)th level, xðjÞ (j ¼ 1; 2; :::; N), the FPF with respect to θ can be estimated by

P̂FðθÞ ¼ pm−10 ⋅ 1 N

∑ N

j ¼ 1 IFm ðxðjÞÞ

f ðxðjÞjθÞ f ðxðjÞjθ0Þ

ð24Þ

As before, the estimation of the FPF based on Eq. (24) does not require repeated evaluations of the limit state function.

4. Choice of instrumental PDF

In the context of RBDO, θ is the set of design parameters, varying over the design intervals, i.e., ½θL; θU�. The selection of the instrumental PDF is crucial for the estimation of the FPF. It should take into consideration all the situations where f ðxjθÞ is under different values of θ in the sense that the statistical error of the FPF is reasonably controlled over the whole design intervals. The proper choice of the instrumental PDF is discussed in the following.

4.1. Weighted Monte Carlo simulation and weighted Subset Simulation

Generally, the nominal value θ0 can be set in the center of the design intervals, that is θ0 ¼ ðθL þ θUÞ=2. For the case that the mean values of the random variables are the design parameters, the c.o.v. of x, Cx, can be taken as deterministic. The mean and standard deviation of the instrumental PDF, θ0 and r0, can be simply chosen as θ0 ¼ ðθL þ θUÞ=2 and r0 ¼ ρCxθ0, respectively, where ρ is a factor, e.g., ρ ¼ 1−10. In order to cover all the important regions under different parameter values, a relatively large standard deviation value should be used. The bigger ρ is, the wider region the instrumental PDF covers. However, an excessively large ρ may lead to poor performance.

In weighted Subset Simulation, the FPF is calculated based on the weighting factors for the samples that lie in the (m−1)th failure level as shown in Eq. (24). The conditional probabilities affect the approximation of FPF, but only through the coefficients of the FPF expression. In this sense, at a given computational cost, more samples should be spent on the (m−1)th level which affects the key terms that contain the variable θ (design parameter) so that

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–114

a more accurate approximation of FPF over different θ values can be obtained. While a proper number of samples are spent in the levels prior to the (m−1)th level to avoid ergodic problems, more samples are spent in the (m−1)th level to increase accuracy.

In order to illustrate how the instrumental PDF works, consider a one-dimensional Gaussian variable x with mean value taken as the design parameter θ∈½θL; θU� ¼ ½3; 7� and a c.o.v. of Cx ¼ 0:1. The parameters of instrumental PDF f ðxjθ0Þ are set as θ0 ¼ 5 and ρ ¼ 3. The PDFs of f ðxjθÞ under different θ values and the PDF of f ðxjθ0Þ are shown in Fig. 1. It can be seen that the instrumental PDF f ðxjθ0Þ, whose sampling center is set in the center of the design interval, approximately covers all the regions of the cases when θ changes from 3 to 7.

4.2. Weighted importance sampling

There are many schemes for constructing the importance sampling function, such as those based on the Most Probable Point (MPP) [14,15] or adaptive pre-samples [16,17]. The impor- tance sampling function based on MPP is investigated here, as it is commonly used in practice.

Gaussian importance sampling function is commonly used and easy to apply. The random variables can be classified into two groups, one related to the design parameters, denoted by xd ¼ ½x1; x2; :::; xnd �, and the other containing the remaining ones, denoted by xr ¼ ½xndþ1; xndþ2; :::; xndþnr �. The instrumental PDF HðxÞ can be represented as

HðxÞ ¼ ϕðxdjxnd; r ðHÞ xd Þϕðxrjx

n

r ; r ðHÞ xr Þ ð25Þ

1 2 3 4 5 6 7 8 9 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

x

f( x|

)

=3

=5

=7

0 =5, =3

Fig. 1. Instrumental PDF.

L L

L L

L

L

P P P P P

Fig. 2. Composite

where ϕð⋅Þ denotes the Gaussian PDF, xnd(xnr ) and r ðHÞ xd (r

ðHÞ xr ) are the

sampling center and the standard deviation for random variable xd(xr), respectively. Usually the sampling center is set at the MPP. Here it can be set as follows: first, set θ0 ¼ ðθL þ θUÞ=2, and then determine the corresponding MPP, i.e., xn ¼ ½xnd; xnr �. The standard deviations for HðxÞ can be set as rðHÞxd ¼ ρCxd xnd, r

ðHÞ xr ¼ rxr where ρ is

a factor, e.g., ρ ¼ 1−10. In order to cover all the important regions under different parameter values, a relatively large standard deviation value can be selected.

From Eqs. (6), (13) and (24), it seems that the estimator of the FPF under the weighted approach can be applied globally when the samples are sufficiently spread off. However, in reality this is not easy to achieve. The proposed approach based on the instru- mental PDF for the construction of FPF appears to be a ‘local’ approximation. At a given computational effort, the approximation does not guarantee that it is valid over arbitrary size of the design parameter domain. When the design value θ is too far from the sampling center, the samples generated by the instrumental PDF may not cover the corresponding important region of failure domain, and the proposed approach may fail to obtain a satisfac- tory estimation of FPF in the vicinity of θ. Of course, it can be improved by increasing the computational effort, but it is usually not feasible in problems of practical interest. So the limitation of the size of the design intervals still exists. This will be investigated by different examples in the following section.

5. Illustrative examples

In this section, several examples are considered to investigate the performance of the weighted Monte Carlo simulation (weighted MCS or WMCS), the weighted importance sampling (weighted IS or WIS) and the weighted Subset Simulation (weighted SS or WSS). First, a static reliability problem with explicit limit state function is considered. In the second example, a structure is modeled with finite elements and the response function is only implicitly given. A multi-story building with uncertainties in the structural parameters and dynamic load is also considered, which is a challenging problem with an extremely large number of uncertain parameters.

5.1. Example 1

Consider a composite beam shown in Fig. 2. This example is modified from [20]. The beam is A (mm) wide by B (mm) high by L (mm) long with the Young's modulus Ew. It has an aluminum plate with Yong's modulus Ea and a net section fastened to its bottom face which is C (mm) wide by D (mm) high. Along the beam, six external vertical forces, P1, P2, P3, P4, P5 and P6 (kN), are applied at six different locations, L1, L2, L3, L4, L5 and L6. A total of

D C

A

B

L

P

beam.

95 96 97 98 99 100

101 102 103 104 105

190

195

200

205

210 -8

-6

-4

-2

0

1

2

lo g1

0[ P F

( 1,

2)

]

190 195

200 205

210

95

100

105 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

1

2

C ov

[P F(

1,

2) ]

Fig. 3. Two-dimension FPF PFðθ1; θ2Þ and its c.o.v. with 106 samples by weighted MCS (solid lines) for Example 1. Exact values (circles) are obtained by direct MCS with 106 samples for each value.

Table 1 Distribution information of random variables for Example 1.

No. Random variable Mean C.o.v. Distribution

1 A (mm) θ1 ¼ μA 0.01 Normal 2 B (mm) θ2 ¼ μB 0.01 Normal 3 C (mm) 80 0.01 Normal 4 D (mm) 20 0.01 Normal 5 L1(mm) 200 0.01 Normal 6 L2 (mm) 400 0.01 Normal 7 L3 (mm) 600 0.01 Normal 8 L4 (mm) 800 0.01 Normal 9 L5 (mm) 1000 0.01 Normal

10 L6 (mm) 1200 0.01 Normal 11 L (mm) 1400 0.01 Normal 12 P1 (kN) 15 0.2 Gumbel 13 P2 (kN) 15 0.2 Gumbel 14 P3 (kN) 15 0.2 Gumbel 15 P4 (kN) 15 0.2 Gumbel 16 P5 (kN) 15 0.2 Gumbel 17 P6 (kN) 15 0.2 Gumbel 18 Ea (GPa) 70 0.01 Normal 19 Ew (GPa) 8.75 0.01 Normal

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–11 5

19 random variables are considered in this example

x ¼ ½A; B; C; D; L1; L2; L3; L4; L5; L6; L; P1; P2; P3; P4; P5; P6; Ea; Ew� ð26Þ The distribution information of the random variables is given in

Table 1. Failure is defined as maximum binding normal stress of the beam exceeding the allowable tensile stress S (strength). The limit state function is then given by

gðxÞ ¼ S−smaxðxÞ ð27Þ where S ¼ 0.0195 GPa is the allowable tensile stress; smaxðxÞ is the maximum stress given by

smaxðxÞ ¼ max skðxÞ : k ¼ 1; :::; 6 � �

ð28Þ and skðxÞ is the stress of the cross-section where Pk is applied. It can be explicitly obtained as

s1ðxÞ ¼ ½ðL1=LÞ∑6i ¼ 1PiðL−LiÞ�YmaxðxÞ

WðxÞ ;

skðxÞ ¼ ½ðLk=LÞ∑6i ¼ 1PiðL−LiÞ−∑k−1i ¼ 1PiðLk−LiÞ�YmaxðxÞ

WðxÞ ðk ¼ 2; :::; 6Þ

ð29Þ where

YmaxðxÞ ¼ 0:5AB2 þ DCðB þ DÞEa=Ew

AB þ DC Ea=Ew ð30Þ

WðxÞ ¼ AB 3

12 þ AB YmaxðxÞ−

B 2

� �2 þ CD

3Ea 12Ew

þ CDEa Ew

D 2 þ B−YmaxðxÞ

� �2 ð31Þ

The mean values of A and B are taken as design parameters, i.e., θ ¼ ½θ1; θ2� and θ1 ¼ μA∈[95, 105] (mm), θ2 ¼ μB∈[190, 210] (mm).

5.1.1. Weighted MCS In the weighted MCS, N¼106 samples are used to obtain the

two-dimension FPF with an instrumental PDF f ðxjθ0Þ, where θ0 ¼ ½100; 200� (mm) and ρ ¼ 3. The results are shown in Figs. 3–5. In each figure, the exact results calculated by direct MCS are also shown with circles. Here 10 values are uniformly selected along each design interval, thus the direct MCS method needs to analyze 102 times to obtain the discrete values of failure probabilities with 106 samples for each analysis.

Figs. 4 and 5 show the variation of the FPF along each individual dimension. It can be seen that as the parameter θ1 varies from 95 to 105, the FPF varies from about 10−2 to 10−3. As

the parameter θ2 varies from 190 to 210, the FPF varies from about 10−1 to 10−5, spanning several orders of magnitude. The results of the direct MCS and the weighted MCS are consistent with each other. Note that the same number of samples has been used for each single run of the direct MCS (for each value of the design parameters) and the weighted MCS. In the latter repeated relia- bility analyses are avoided and reduced to one single simulation run. It can also be seen that the c.o.v.s of the weighted MCS are mostly slightly higher than those of the direct MCS. This may be due to the fact that the information of samples in a single run is spread over a range of values as ρ ¼ 3 is selected. This seems inevitable and is a trade-off between ‘efficiency’ and ‘robustness’.

5.1.2. Weighted IS In the implementation of the weighted IS, the instrumental PDF

HðxÞ is chosen as the Gaussian density function where the sampling center is taken as the design point when ½θ1; θ2�¼[100, 200] (mm). The design point is approximately calculated as xn ¼ [99.78, 199.09, 79.96, 19.99, 200.03, 400.14, 600.34, 799.22, 998.82, 1198.40, 1406.90, 15.932, 17.796, 16.864 18.723, 17.497, 16.271, 69.970, 8.7538]. The standard deviations of A and B in the instrumental PDF are set as 3 times of the original ones, i.e., ρ ¼ 3 and ½sðHÞA ; s

ðHÞ B � ¼ 3½sA; sB� ¼ [3, 6] (mm). For the other vari-

ables the standard deviations are taken the same values as the original ones (see details in Eq. (25)).

In the simulation, N¼8000 samples have been used to obtain the two-dimension FPF. The results are shown in Fig. 6. The results

95 96 97 98 99 100 101 102 103 104 105 10

-3

10-2

10-1

θ 1

P F (θ

1)

WMCS 95% Confidence interval 95% Confidence interval Exact Sampling center

95 96 97 98 99 100 101 102 103 104 105 0

0.02

0.04

0.06

0.08

0.1

0.12

θ 1

C ov

(P F (θ

1) )

WMCS Exact Sampling center

Fig. 4. FPF PFðθ1; θ2 ¼ 200Þ (solid line) and its 95% confidence band (dashed line) by weighted MCS with 106 samples for Example 1. Exact values (circles) are obtained by direct MCS with 106 samples for each value.

190 192 194 196 198 200 202 204 206 208 210 10-5

10-4

10-3

10-2

10-1

100

P F (θ

2)

WMCS 95% Confidence interval 95% Confidence interval Exact Sampling center

190 192 194 196 198 200 202 204 206 208 210 0

0.05

0.1

0.15

0.2

0.25

C ov

(P F (θ

2) )

WMCS Exact Sampling center

θ2

θ2

Fig. 5. FPF PFðθ1 ¼ 100; θ2Þ (solid line) and its 95% confidence band (dashed line) by weighted MCS with 106 samples for Example 1. Exact values (circles) are obtained by direct MCS with 106 samples for each value.

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–116

of the direct MCS with 106 samples for each FPF value are also shown using circles. It can be seen that the results of the weighted IS are consistent with those of the direct MCS over the whole region. The c.o.v.s are below 40% in the most part of the design region.

Figs. 7 and 8 show the variation of the FPF along each individual dimension. The results of the direct MCS and weighted IS are also consistent with each other. The c.o.v.s of the weighted IS are mostly below 30% except in the tail of the design region which is a little far from the sampling center.

5.1.3. Weighted SS In the weighted SS method, Ni ¼ 1000 and N ¼ 3000 (a total of

4000 samples as m¼1) are used to obtain the two-dimension FPF with an instrumental PDF f ðxjθ0Þ, where θ0 ¼ ½100; 200� (mm) and ρ ¼ 3. The results are shown in Fig. 9. Note that the confidence band of FPF is not given here, as the c.o.v. of the failure probability by the weighted SS cannot be obtained explicitly. In the figure, the exact results by the direct MCS with 106 samples are also shown with circles. It can be seen that the results of the weighted SS are mostly consistent with those of the direct MCS over the whole region. Fig. 10 shows the variation of the FPF along each individual

dimension. From the figure, it can be seen that the results by the weighted SS method are mostly consistent with the exact values except in the tail of the design region where the failure probability values are small.

As mentioned above, the weighted approach is a local approx- imation. From this example, it seems that the weighted approach can obtain the satisfactory estimation of the FPF (ranges from 10−1

to 10−5 ) when θ is in a range of about θ0 75Cxθ0 where θ0 is the center of the design intervals, and moderate computational effort is used, i.e., 106 samples for WMCS, 8000 for WIS and 4000 for WSS, respectively.

5.2. Example 2

A three-box flap structure is shown in Fig. 11. It consists of 16 plates and 28 bar elements. The 28 bars can be classified by different directions, i.e., X, Y and Z, into three kinds, the length of bar of each kind is the same, denoted by L1, L2 and L3 respectively. The section area for all the bars is the same, denoted by A. The modulus of elasticity of bars and plates is E and the Poisson ratio is taken as constant 0.3; P is the external load applied on the nodes; T is the thickness of the plates. The finite element model of the three-box structure is shown in Fig. 12.

95 96 97 98 99 100 101 102

103 104 105190

195

200

205

210 -7

-6

-5

-4

-3

-2

-1

0

θ1

θ2

lo g 1

0[ P

F (θ

1, θ 2

)]

190 195

200 205

210

95

100

105 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

θ1 θ2

C ov

[P F

(θ 1,

θ 2 )]

Fig. 6. Two-dimension FPF PFðθ1; θ2Þ with 8000 samples by weighted IS (solid line) for Example 1. Exact values by direct MCS with 106 samples are also shown with circles.

95 96 97 98 99 100 101 102 103 104 105 10-4

10-3

10-2

10-1

θ1

P F (θ

1)

WIS 95% Confidence interval 95% Confidence interval Exact Sampling center

95 96 97 98 99 100 101 102 103 104 105 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

θ1

C ov

(P F (θ

1) )

WIS Exact Sampling center

Fig. 7. FPF PFðθ1; θ2 ¼ 200Þ (solid line) and its 95% confidence band (dashed line) by weighted IS with 8000 samples for Example 1.

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–11 7

There are seven Gaussian random variables, i.e., L1, L2, L3, A, E, P and T. Their distribution information is shown in Table 2. Failure is defined as the displacement of Node 16 exceeding the allowable displacement. The limit state function is given by

gðxÞ ¼ Dt−dðL1; L2; L3; A; E; P; TÞ ð32Þ

where Dt ¼3 mm is the allowable displacement; dðL1; L2; L3; A; E; P; TÞ is the displacement of Node 16 which can be obtained by finite element analysis.

Here the mean value of the thickness T is viewed as a design parameter, i.e., θ ¼ μT . The design interval for θ is [1, 6] (mm).

In this Example, the same instrumental PDF f ðxjθ0Þ, where θ0 ¼ 3:5 (mm) and ρ ¼ 3, is used for weighted MCS and weighted SS. In the weighted MCS, N¼10,000 are used to analyze and the results are shown in Fig. 13. In the weighted SS, Ni ¼1000 and N¼3000 are used to analyze and the results are shown in Fig. 14. Direct MCS is also applied to calculate the failure probability which is taken as the exact value. For each value of design parameter, the direct MCS uses 104 samples to obtain the corresponding value of the failure probability.

It can be seen from the figures that as θ varies from 1 to 6 (mm), the failure probability varies from about 10−1 to 10−4, spanning several orders of magnitude. The FPF estimates obtained by the weighted MCS and the weighted SS are consistent with the exact values (dots). Note that the computation of each exact value requires one reliability analysis.

From this example, it is also seen that the weighted approach can obtain satisfactory estimation of FPF when θ is in about θ0 75Cxθ0.

5.3. Example 3

Consider a six-story steel frame as shown in Fig. 15. The floors are assumed to be rigid. It is assumed to be linear-elastic and classically damped. The damping ratio of the i-th modal is denoted as ζi(i ¼ 1; 2; :::; 6). Passive dampers, which are modeled as linear viscous damping elements, are also installed as diagonal bracings. The damping coefficient of the damper in the i-th story (i ¼ 1; 2; :::; 6) is denoted as Cdi whose value can be adjusted.

A total of 24 structural random variables are considered in this example, which include the original modal damping ratios of the building, ζi(i ¼ 1; 2; :::; 6), the damping coefficient of each damper Cdi(i ¼ 1; 2; :::; 6), the mass of each story, m1, m2, m3, m4, m5 and m6, and the stiffness of each story, k1, k2, k3, k4, k5 and k6. The distribution information of the random variables is given in Table 3.

The structure is subject to stochastic ground acceleration €agðtÞ which is modeled as Gaussian white noise WðtÞ. The failure is defined as the peak interstory drift ratio over all stories of the structure exceeding a threshold level b. Three damage levels,

190 192 194 196 198 200 202 204 206 208 210 10-5

10-4

10-3

10-2

10-1

100

θ2

P F (θ

2)

WIS 95% Confidence interval 95% Confidence interval Exact Sampling center

190 192 194 196 198 200 202 204 206 208 210 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

θ2

C ov

(P F (θ

2) )

WIS Exact Sampling center

Fig. 8. FPF PFðθ1 ¼ 100; θ2Þ (solid line) and its 95% confidence band (dashed line) by weighted IS with 8000 samples for Example 1.

95 96 97 98 99 100 101 102 103 104 105190

195

200

205

210 -10

-8

-6

-4

-2

0

θ1

θ2

lo g 1

0[ P

F (θ

1, θ

2) ]

Fig. 9. Two-dimension FPF PFðθ1; θ2Þ (solid line) by weighted SS with N ¼ 2000 and exact results are shown with circles for Example 1.

95 96 97 98 99 100 101 102 103 104 105 10-3

10-2

10-1

θ1

190 192 194 196 198 200 202 204 206 208 210 10-6

10-5

10-4

10-3

10-2

10-1

100

θ2

P F (θ

1) P

F (

θ 2 )

Fig. 10. FPF (solid line) by weighted SS with N ¼ 2000 for Example 1 (exact results are shown with dots). a) PF (y1, y2 ¼ 200) and b) PF(y1¼100,y2).

x

y

z

PP

PP

PP

Fig. 11. Three-box structure.

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–118

b¼0.2% (Full Operational), b¼0.5% (Operational), b¼1.5% (Life- Safety) are considered for reliability assessment. Here the case for b¼1.5% is considered first. A duration of T¼10 s and time interval of Δt¼0.02 s are assumed, this gives nt ¼ T=Δt¼500 input random variables in the discrete approximation for WðtÞ at tk ¼ kΔt (k ¼ 1; 2; :::; nt), i.e., WðtkÞ ¼ ZðtkÞ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2πS=Δt

p where S ¼ 0:01 m2=s3 is

the spectral intensity and Zk ¼ ZðtkÞ are independent identical

Fig. 12. Finite element model of three-box structure.

Table 2 Distribution information of three-box structure.

Variable A L1 L2 L3 E P T

Mean 100 mm2 600 mm 200 mm 400 mm 71 GPa 1.5 kN θ (mm) C.o.v. 0.1 0.1 0.1 0.1 0.1 0.2 0.1

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 10-5

10-4

10-3

10-2

10-1

100

P F (θ

)

WMCS 95% Confidence interval 95% Confidence interval Sampling center Exact

θ

Fig. 13. FPF PFðθÞ and its 95% confidence band by weighted MCS with 10,000 samples and ρ¼3 (exact results are showed by dots).

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

10-4

10-3

10-2

10-1

100

P F (θ

)

WSS Exact Sampling center

Fig. 14. FPF PFðθÞ (solid line) by weighted SS with Ni ¼1000, N¼3000 and ρ¼3 (exact results are showed by dots).

Fig. 15. Six-story steel frame structure with passive dampers.

Table 3 Distribution information of structural random variables.

Variable Mean C.o.v. Distribution

Cdi(i ¼ 1; 2; :::; 6) (103 Ns/mm) θ 0.1 Normal ζi(i ¼ 1; 2; :::; 6) (%) 2 0.2 Normal m1 (10

3 kg) 283 0.05 Normal m2 (10

3 kg) 263 0.05 Normal m3 (10

3 kg) 256 0.05 Normal m4 (10

3 kg) 255 0.05 Normal m5 (10

3 kg) 247 0.05 Normal m6 (10

3 kg) 215 0.05 Normal k1 (kN/mm) 123 0.1 Normal k2 (kN/mm) 367 0.1 Normal k3 (kN/mm) 246 0.1 Normal k4 (kN/mm) 246 0.1 Normal k5 (kN/mm) 175 0.1 Normal k6 (kN/mm) 175 0.1 Normal

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–11 9

distributed Gaussian random variables. When all the structural variables are fixed at the mean values, the natural frequencies of the first two modes of vibration are computed to be 1.08 Hz and 2.97 Hz, respectively.

Here the mean values of damping coefficients of dampers are viewed as the design parameters whose effects on system relia- bility are studied, i.e., θi ¼ μCdi (i ¼ 1; 2; :::; 6). The design interval for θi is [0.6, 1.5] (10

3 Ns/mm). In the context of seismic design, the mean values of different dampers, θi(i ¼ 1; 2; :::; 6), are assigned the same value, θi ¼ θ. In the following, the proposed weighted IS and weighted SS are applied to obtain the six-dimension FPF, but only the variation of the FPF with respect to θ is shown.

In this example, the instrumental PDF for the weighted IS, HðxÞ, is selected as follows: for variables Cdi, the Gaussian distribution is used with the sampling center CðHÞdi ¼ 1 (10

3 Ns/mm) and the standard deviation sðHÞCdi ¼ 0:1ρC

ðHÞ di (10

3 Ns/mm), respectively. For variables ζi, mi and ki(i ¼ 1; 2; :::; 6), the original PDFs of ζi, mi and ki are used respectively. For variables Zk, the importance sampling

0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2

10-4

10-3

10-2

10-1

( 103 Ns/mm)

( 103 Ns/mm)

WIS 95% Confidence interval 95% Confidence interval Exact Sampling center

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

10-4

10-3

10-2

WIS 95% Confidence interval 95% Confidence interval Exact Sampling center

P F (θ

) P

F (θ

)

Fig. 16. FPF PFðθÞ and its 95% confidence band by weighted IS with 8000 samples (exact results are showed by circles). a) r =1 and b) r =3.

0.7 0.8 0.9 1 1.1 1.2 1.3

10-4

10-3

10-2

10-1

P F (θ

)

WSS Exact Sampling center

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 10-6

10-5

10-4

10-3

10-2

10-1

( 103 Ns/mm)

( 103 Ns/mm)

WSS Exact Sampling center

P F (θ

)

Fig. 17. FPF PF ðθÞ by weighted SS with N¼8000 (solid line) (exact results are showed by circles). a) r =1 and b) r =3.

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–1110

function proposed in [15] for linear structure subjected to random excitation is adopted here. The instrumental PDF for the weighted SS, f ðxjθ0Þ, is selected as follows: for variables Cdi, the Gaussian distribution is used with the sampling center CðHÞdi ¼ 1 (10

3 Ns/mm) and the standard deviation sðHÞCdi ¼ 0:1ρ C

ðHÞ di (10

3 Ns/mm), respec- tively. For other variables, the original PDFs are used. For both WIS and WSS, ρ ¼ 1 and 3 are used, respectively.

Fig. 16 shows the estimate FPF versus design parameter θ and the 95% confidence band obtained by the weighted IS with N¼8000 samples. Fig. 17 shows the estimate FPF versus design parameter θ obtained by the weighted SS with Ni ¼1000 and N¼8000. It is seen that both the results of the weighted IS and weighted SS are approximately consistent with the exact values (denoted by circles) when θ is near the intended sampling region, but they give a poor estimation when θ is far from the center for both ρ ¼ 1 and 3.

6. Conclusions

A weighted approach has been proposed to efficiently estimate the FPF of a system which is the function of distribution para- meters. The FPF with respect to the design parameters can provide important information for structural design under uncertainty. The key idea of the weighted approach is to rewrite the failure probability estimate for a given set of random samples in simula- tion as a function of distribution parameters. The proposed

method relies on the fact that the important regions of the failure domains for different parameter values overlap each other, and so they could be explored by only one properly selected PDF, the instrumental PDF. One attractive feature of this approach is that it only requires a single reliability analysis run under a pre-selected instrumental PDF. Three reliability analysis methods have been incorporated into the proposed approach, i.e., weighted MCS, weighted IS and weighted SS.

Numerical examples have been considered to investigate the feasibility and efficiency of the proposed approach. The static exam- ples (with implicit or explicit limit state function) show that for general problem, the proposed approach can give a satisfactory estimation of the FPF when θ is in about θ0 75Cxθ0 where θ0 is the center of design intervals. The dynamic example shows that for high dimensional problem, the corresponding interval is about θ0 7Cxθ0.

The proposed approach may be used for RBDO and reliability sensitivity analysis. Current work focuses on combining the proposed approach and sequential approximation optimization framework for RBDO.

Acknowledgments

Part of work in this paper was performed during the author′s postdoctoral research under the supervision of Prof. Siu-Kui Au

X. Yuan / Probabilistic Engineering Mechanics 34 (2013) 1–11 11

formerly at the City University of Hong Kong. The author would like to acknowledge financial support from the Research Grants Council of the Hong Kong Special Administrative Region (General Research Fund 9041550, CityU 110210) and the National Natural Science Foundation of China (Grant no. 51105309).

References

[1] Schuëller G, Jensen HA. Computational methods in optimization considering uncertainties—an overview. Computer Methods in Applied Mechanics and Engineering 2008;198(1):2–13.

[2] Enevoldsen I, Sorensen JD. Reliability-based optimization in structural engi- neering. Structural Safety 1994;15(3):169–96.

[3] Valdebenito MA, Schuëller GI. A survey on approaches for reliability-based optimization. Structural and Multidisciplinary Optimization 2010;42 (5):645–63.

[4] Gasser M, Schueller GI. Reliability-based optimization of structural systems. Mathematical Methods of Operations Research 1997;46(3):287–307.

[5] Jensen HA. Structural optimization of linear dynamical systems under sto- chastic excitation: a moving reliability database approach. Computer Methods in Applied Mechanics and Engineering 2005;194(12–16):1757–78.

[6] Ching J, Hsieh YH. Approximate reliability-based optimization using a three- step approach based on subset simulation. Journal of Engineering Mechanics —ASCE 2007;133(4):481–93.

[7] Ching J, Hsieh YH. Local estimation of failure probability function and its confidence interval with maximum entropy principle. Probabilistic Engineer- ing Mechanics 2007;22(1):39–49.

[8] Zou T, Mahadevan S. A direct decoupling approach for efficient reliability- based design optimization. Structural and Multidisciplinary Optimization 2006;31(3):190–200.

[9] Au SK. Reliability-based design sensitivity by efficient simulation. Computers and Structures 2005;83(14):1048–61.

[10] Wang Y. Reliability-based economic design optimization of spread founda- tions. Journal of Geotechnical and Geoenvironmental Engineering 2009;135 (7):954–9.

[11] Wang Y, Au SK, Kulhawy FH. Expanded reliability-based design approach for drilled shafts. Journal of Geotechnical and Geoenvironmental Engineering— ASCE 2011;137(2):140–9.

[12] Taflanidis AA, Beck JL. An efficient framework for optimal robust stochastic system design using stochastic simulation. Computer Methods in Applied Mechanics and Engineering 2008;198(1):88–101.

[13] Taflanidis AA, Beck JL. Stochastic Subset Optimization for reliability optimiza- tion and sensitivity analysis in system design. Computers & Structures 2009;87(5–6):318–31.

[14] Melchers RE. Importance sampling in structural systems. Structural Safety 1989;6(1):3–10.

[15] Au SK, Beck JL. First excursion probabilities for linear systems by very efficient importance sampling. Probabilistic Engineering Mechanics 2001;16 (3):193–207.

[16] Bucher CG. Adaptive sampling—an iterative fast Monte Carlo procedure. Structural Safety 1988;5:119–26.

[17] Ang GL, Ang AH-S, Tang WH. Optional importance sampling density estimator. Journal of Engineering Mechanics—ASCE 1992;118(6):1146–63.

[18] Au SK, Beck JL. Estimation of small failure probability in high dimensions by subset simulation. Probabilistic Engineering Mechanics 2001;16(4):263–77.

[19] Au SK, Ching J, Beck JL. Application of subset simulation method to reliability benchmark problems. Structural Safety 2007;29(3):183–93.

[20] Huang B, Du X. Probabilistic uncertainty analysis by mean-value first-order Saddlepoint approximation. Reliability Engineering and System Safety 2008;93(2):325–36.