Local Estimation of Failure Probability Function with Direct Monte Carlo Simulation

Jianye Ching¹and Yi-Hung Hsieh²

1Assistant Professor, Dept of Construction Engineering, National Taiwan University of Science and Technology, Taiwan; Tel: 886-2-27301055; Fax: 886-2-27376606;

jyching@gmail.com

2Graduate Student, Dept of Construction Engineering, National Taiwan University of Science and Technology, Taiwan;yihung.hsieh@gmail.com

Abstract

An approach is developed to locally estimate the failure probability of a system under various design values by using simple Monte Carlo simulation. Although it seems to require numerous reliability analysis runs to locally estimate the failure probability function, which is a function of the design variables, the approach only requires simple Monte Carlo simulations. It is also possible to find the confidence interval of the failure probability function as well as estimate the gradient of the logarithm of that function with respect to the design variables. The use of the new approach is demonstrated with two simulated examples. The results show that the new approach can effectively locally estimate the entire failure probability function and the confidence interval. The approach should be valuable for reliability-based design and reliability sensitivity analysis.

Introduction

Uncertainties are often encountered in civil engineering problems, and reliability analyses (e.g.: Benjamin and Cornell 1970; Thoft-Christensen and Murotsu 1986; Ang and Tang 1984; Der Kiureghian 2000; Au and Beck 2001) are the main tools of quantifying them. In fact, reliability analyses are key components of many other research areas, such as reliability-based optimization and design (Enevoldsen and

Sorensen 1994; Gasser and Schueller 1997; Royset et al. 2001; Papadrakakis and Lagaros 2002). In practice, it is often of interest to calculate reliability of the target system under various design settings. Taking reliability-based optimization as an example, let us consider the following optimization problem:

*

minc0( ) s t. . P_F( ) P_F c_i( ) 0 i 1,...,m

     (1)

where is the design parameter; P_F(is the failure probability of the target system given the design variable is ; c₀() is the objective function; {c_i():i=1,…,m}are deterministic constraints of ; P_F^* is the target probability of failure. If the entire PF(

function can be obtained beforehand, the reliability-based optimization in (1) can be transformed into an ordinary optimization problem. Otherwise, the optimization problem in (1) can be highly non-trivial because of the reliability constraint. Another example is reliability sensitivity analysis (Bjerager and Krenk 1989; Karamchandani and Cornell 1992; Enevoldsen 1994), where reliability sensitivity under design variable variation, i.e.

sensitivity of PF(with respect to ,is the main target.

However, finding the failure probability function PF(, denoted by FPF, is not trivial. It seems that numerous reliability analyses need to be conducted for various  locations. This approach has been taken by Jensen (2005), where he adopted a linear function of to locally approximate log[PF(of a deterministic linear system subject to stochastic excitation. The number of the to-be-determined coefficients in the linear function is equal to the dimension of ,denoted by n, so at least nreliability analyses are required to determine those coefficients. In Gasser and Schueller (1997), a quadratic function is assumed, so the minimum required number of reliability analyses rapidly increases to n+n(n)/2. However, for complicated problems, large number of reliability analyses can be computationally infeasible. This may explain why most examples in relevant reliability-based optimization literature are academic type, i.e.

simple systems subject to few uncertainties. There is a strong need to enhance our ability of estimating the PF( function for general nonlinear systems subject to high-dimensional uncertainties.

In Ching and Hsieh (2006), a novel approach is proposed to locally estimate the failure probability function PF(and its confidence interval by using Subset Simulation (Au and Beck 2001). In this paper, it is shown that the same problem can be easily solved by using Monte Carlo simulation (MCS) and the maximum entropy principle (Jaynes 1957). Moreover, it is possible to estimate the gradient of the logarithm of the FPF with respect to the design variables. Those features are essential for reliability-based design and reliability sensitivity computations.

Problem Definition

Let Z be the uncertain variables of the target system andbe the design variables. Given a design,the probability of failure of the target system is

( ) ( | ) ( | )

F

PF  P F  f Z dZ



 



where F is the failure event, f(Z|) is the probability density function (PDF) of Z conditioned on ,_Fis the failure domain in the Z space. The quantity PF(in (2) is exactly the FPF. It is usually of interest to estimate the FPF only in the design region D where the failure probability is small since the design is expected to be safe.

Estimation of the Failure Probability Function

The basic idea is to treatas uncertain although it is, in fact, controlled by the designer.

First, let us assign a prior probability density function (PDF) to ,denoted by f(). By the Bayes rule,

where f(|F) is the PDF of conditioned on the failure event; P(F) is the probability the failure when both Z andare considered to be uncertain:

( ) ( | ) ( )

F

P F f Z  f dZd







where_F= the failure region in the (Z,)space. Note that P(F) is different from the FPF PF(: the former is the failure probability if both Z andare uncertain, while the latter is the failure probability if only Z is uncertain. From (3), it can be seen that the FPF can be determined if P(F), f(|F) and f()are known.

It is trivial to compute f() since we have the freedom to choose it, so it is always possible to choose f() so that it is easy to evaluate. P(F) can be estimated using any efficient reliability analysis in which both Z and are treated as uncertain. Assuming that we can obtain samples of from f(|F), we can then use these samples to estimate the underlying PDF f(|F).

Monte Carlo simulation. In this section, MCS is introduced to estimate P(F) and obtain samples from f(|F). According to the Law of Large Number,

 

where Z⁽ⁱ⁾ and ⁽ⁱ⁾ are the i-th independent identically distributed (i.i.d.) sample pair drawn from the joint PDF f(Z,)=f(Z|)f(); I_F(Z⁽ⁱ⁾,⁽ⁱ⁾) is the indicator function of the failure event, i.e. IF(Z⁽ⁱ⁾,⁽ⁱ⁾)=1 if Z⁽ⁱ⁾ and ⁽ⁱ⁾ represent a failure state and IF(Z⁽ⁱ⁾,⁽ⁱ⁾)= 0 otherwise. Note that the samples satisfying failure are distributed as f(Z,|F), so the  parts of the samples are distributed as f(|F). With these samples, it is possible to estimate f(|F), as described in the next section.

Estimating f(|F): maximum entropy method. The underlying PDF f(|F) can be estimated from its samples. One option of estimating f(|F) is the kernel method, i.e. the sum of Gaussian kernels centered at the sample points. However, kernel methods are not implemented in this paper because it requires many samples to converge well, especially when the dimension of is large. According to past experience (see Gasser and Schueller (1997) and Jensen (2005)), f(|F) is usually quite smooth in the low failure probability region. In fact, in Gasser and Schueller (1997) and Jensen (2005), linear and quadratic functions ofare used to approximate log[P_F(in the low probability region.

If this approximation is appropriate, it is a pity to abandon this valuable prior information, which is the case when either kernels or histograms are chosen.

In this study, the maximum entropy method (Zellner and Highfield 1988; Ormoneit and White 1999) is adopted to estimate the PDF. The idea is to first estimate the first moment (mean) of f(|F) from the samples and then to find the PDF which maximizes the entropy subject to the first moment constraint. The rationale of this approach is that the PDF maximizing entropy is the least subjective PDF subject to the moment information. The resulting PDF estimate is usually quite smooth. The maximum entropy method of estimating f(|F) is stated as following:

max( ) log[ ( )] ( )

whereiis the i-th design variable;iis the sample mean of theisamples. Note that the variable in this optimization problem is the entire g(.) function. The optimal solution of

(6) is the maximum entropy estimate of f(|F). Differentiating the Lagrangian function with respect to g(.) obtains the following unique optimal solution:

 

wheres are the Lagrange multipliers that satisfy the following equations:

0 0

The optimal solution^*can be found by solving (8) (there are nof them) with standard Newton methods. Once the optimal solution^*is found, the maximum entropy estimator g^*() is exactly exp(-^T^*), which is usually a satisfactory approximation of f(|F) in the region where the samples are abundant regardless how linear log[f(|F)] is. This is why the result of the proposed method should be considered as a local approximation of the FPF. However, if log[f(|F)] is roughly linear in region D, which is the case for many examples, the result of the proposed approach can be satisfactory not only in the region with manysamples but also in the region with few samples. This will be seen in a later section where examples are demonstrated.

Final estimate of PF(The final estimate of PF(is

The required computational cost by the method of estimating P_F(may be much less than required by the previous methods, i.e. Gasser and Schueller (1997) and Jensen (2005), especially when dimension is high. Moreover, with the maximum entropy approach, it is possible to estimate the gradient of the logarithm of that function with respect to the design variables as well as find the confidence interval of the failure probability function (Ching and Hsieh 2006), summarized as follows:

 ²  ²

and (g^*()) is the coefficient of variation of the g^*() estimator, whose formula can be found in Ching and Hsieh (2006).

Examples

We present two numerical examples in this section: (a) an earth-fill dam subject to various uncertainties; and (b) a retaining wall subject to various uncertainties. The goal of the examples is to demonstrate the use of the proposed approach, i.e. locally estimate the FPF PF(in the design region D, where failure, design parameters and design regions D will be defined differently in each example. For all examples, is treated as uncertain and its prior PDF is taken to be uniform over the design region D.

Example 1: earth fill dam. This example consists of a cross section of an earth fill dam which is depicted in Figure 1, where the material properties are uncertain. It is proposed to model the elastic and shear moduli as lognormal random fields E(x,z) and G(x,z), respectively:

where Y(x,z) is a zero-mean, unit-variance Gaussian random field with a given autocorrelation R_Y(Δx,Δz)=exp(-Δx/10-Δz/3) (Note: the mean of Y = e^0.5and variance of Y = e²-e).

The cross section is discretized into triangular, 3-node, finite elements and the resulting system is solved under the plane-strain condition in the linear elastic regime.

The dam is subject to a fixed loading q = 30 kPa/m at the top, self weight (soil mass density = 1800 kg/m³) and the hydrostatic loading where the top level of the water is 2m below the crest. Failure is defined when the Mohr-Coulomb criterion is satisfied in at

least one element: = c+_ntan(), where is the shear stress, c = 150 KPa/m is the cohesive strength, _n is the normal stress and = 40^o is the friction angle. The design parameters (₁,₂) are the width and height of the dam, respectively, and the design region is ₁[50,120] and ₂[15,25]. Figure 2 shows the estimated FPF and the 95%

confidence intervals from the approach with MCS sample size N = 10000 for various₁ values when₂is fixed at 15, 20 and 25. The results from brute-force reliability analyses (MCS with sample size = 100,000) are plotted, and it is seen that the actual FPF should lie within the confidence interval for various₁and₂values.

Figure 1. The cross section of the earth fill dam for Example 1.

50 100

Figure 2. The PF(estimate made by the proposed approach with MCS sample size N = 10000 (solid line) and its 95% confidence interval (dashed line) for Example 1. The dots

are the results of brute-force analyses.

Example 2: retaining wall. Consider the retaining wall in Figure 3 that is subject to self weight and earthquake excitation. The uncertainties include the density and friction angle of the backfill cohesiveless soil, friction angle between the wall and soil and peak horizontal and vertical earthquake acceleration kh and kv. All uncertain variables are modeled as Gaussian random variables with means equal to [_kh _kv] = [20KN/m³ 32^o 19^o 0.2g 0.067g] and standard deviations equal to [_kh _kv] = [1.5KN/m³ 3^o 3^o 0.02g 0.0067g]. The design parameters include the height H(1), the backfill angle (₂), the rare-wall inclination angle(₃) and the weight WW(₄) of the wall. Failure is defined as the retaining wall fails in lateral sliding. The limit-state function is:

 

( , ) _aecos( ) _{W ae}sin( ) tan( )

R Z P   W P    (13)

Figure 3. The cross section of the retaining wall.

where

In this example, we examine the ability of the approach of estimating the sensitivity of log[PF(with respect to around the following design location: H = 4 m, = 15^o,

= 5^oand WW= 450 kN/m. We obtain the sensitivity estimate using the following strategy.

First, we choose the design region D to be a small rectangular centered at the target design location: H[3.5,4.5]m, [10,20]^o, [0,10]^o and W_W[400,500]kN/m. Note that -_i is the sensitivities of the logarithm of the failure probability with respect to _i. Therefore, by estimating the Lagrangian multipliers, the sensitivities of log[P_F(with respect tocan be quantified.

Table 1 shows the results of the estimated Lagrange multipliers from the approach with N = 10000. The results include the sample mean and standard deviation of the estimates obtained using 50 independent runs. It is clear that the height parameter (1) dominates the failure probability. Moreover, the following brute-force MCS approach is used to validate the results: Eight design locations are chosen: the first and second are shifted from the rectangle center in the ₁ direction by amount of +25% and -25%, respectively, of the width of the rectangular in that direction. The third to eighth locations are similarly chosen for the ₂, ₃ and ₄ directions. At each location, brute-force MCS of sample size = 1,000,000 is employed to estimate the failure probability to give eight very accurate failure probability estimates. The logarithms of the eight estimates can then be used to estimate the actual sensitivity of log[PF(with respect to . The brute-force results are also listed in Table 1. It is evident that the sensitivity estimates made by the proposed approach are similar to the actual answers.

Table 1. The sample mean and c.o.v. of the estimated Lagrange multiplier from the proposed approach for Example 2

1(H) 2() 3() 4(W_w)

Mean -2.5606 -0.2242 -0.0261 0.0099

Standard

The proposed approach is able to locally estimate the failure probability function PF(

of a system, a function depending on the deterministic system parameters . The proposed approach employs the maximum entropy principle to estimate PF(. Confidence intervals can be derived with the new approach, and the sensitivity of log[PF(can be readily estimated.

The applicability of the proposed approach is verified using two examples, showing that the approach can effectively estimate P_F(, and the derived confidence intervals mostly contain the actual PF(.

To sum up, the proposed approach of estimating PF(can be attractive because: (a)

it only requires simple Monte Carlo simulations (MCS); (b) pointwise confidence interval of the PF(estimator can be established with computational cost less than brute-force MCS simulations; (c) it is applicable to nonlinear systems subject to high dimensional uncertainties since it is based on MCS. The approach should be valuable for reliability-based design and reliability sensitivity analysis.

Reference

Ang, A.H-S. and Tang, W.H. (1984). Probability Concepts in Engineering Planning and Design. New York, NY: John Wiley & Sons.

Au, S.K. and Beck, J.L. (2001). “Estimation of small failure probability in high dimensionsby subsetsimulation.”Probabilistic Engineering Mechanics, 16(4), 263-277.

Benjamin, J.R. and Cornell, C.A. (1970). Reliability, Statistics and Decision for Civil Engineers, McGraw Hill, New York.

Bjerager,P.and Krenk S.(1989).“Parametricsensitivity in first order reliability theory.”

Journal of Engineering Mechanics, 115(7), 1577-1582.

Ching,J.and Hsieh,Y.H.(2006).“Localestimation offailureprobability function and its confidence interval with maximum entropy principle.”Probabilistic Engineering Mechanics, in press.

Der Kiureghian, A. (2000).“Thegeometry ofrandom vibrationsand solutionsby FORM and SORM.”Probabilistic Engineering Mechanics, 15(1), 81-90.

Enevoldsen, I. (1994). “Sensitivity analysis or reliability-based optimal solution.”

Journal of Engineering Mechanics, 120(1), 198-205.

Enevoldsen,I.and Sorensen,J.D.(1994).“Reliability-based optimization in structural engineering.”Structural Safety, 15(3), 169-196.

Gasser, M. and Schueller, G.I. (1997). “Reliability-based optimization of structural systems.”Mathematical Methods of Operations Research, 46(3), 287-307.

Jaynes, E.T. (1957).“Information theory and statisticalmechanics.”Physical Review, 106(4), 620-630.

Jensen, H.A. (2005). “Structural optimization of linear dynamical systems understochastic excitation: a moving reliability database approach.” Comput.

Methods Appl. Mech. Engrg., 194(16), 1757–1778.

Karamchandani,A.and Cornell,C.A.(1992).“Sensitivity estimation within firstand second order reliability methods.”Structural Safety, 11(2), 95-107.

Ormoneit, D. and White H. (1999). “An efficient algorithm to compute maximum entropy densities.”Econometric Reviews, 18(2), 127-140.

Papadrakakis, M. and Lagaros, N.D. (2002).“Reliability-based structural optimization using neural networks and Monte Carlo simulation.” Comput. Methods Appl.

Mech. Engrg., 191(32), 3491-3507.

Royset,J.O.,Kiureghian,A.D.and Polak,E.(2001).“Reliability-based optimal design of series structural systems.”Journal of Engineering Mech.., 127(6), 607-614 Thoft-Christensen, P. and Murotsu, Y. (1986). Application of Structural Systems

Reliability Theory, Springer Verlag, Berlin.

Zellner,A.and Highfield,R.A.(1988).“Calculation ofmaximum entropy distributions and approximation ofmarginalposteriordistributions.”Journal of Econometrics, 37(2), 195-209.