- Email: [email protected]
An improved Monte Carlo method in structural reliability
Reliability Engineering and System Safety 24 (1989) 275-292
An Improved Monte Carlo Method in Structural Reliability
Attila Csenki Centre for Software Reliability, The City University, Northampton Square, London ECIV 0HB, UK (Received 15 July 1988; accepted 20 July 1988)
A BSTRA C T A Monte Carlo technique is presented for obtaining accurate estimates of component failure probabilities based on small sample sizes. It is an improvement on a recent method by Harbitz and is based on the exclusion from sampling of a maximum hypersphere contained in the safe region as opposed to that of the fl-sphere. After finding the maximum hypersphere in the standard normal space of basic variables a two-stage sampling procedure involving conditional distributions is employed.
NOTATION 8
g, gt, h tl W*, Z* B
C
F,G U. M
Centre of hypersphere Limit state functions Number of basic variables Estimate of failure probability Design points Sampling region Acceptance ratio Distribution functions Distribution function of the central chi-square distribution with n degrees of freedom Indicator function of the set {...} Orthogonal n x n matrix
275 Reliability Engineering and System Safety 0951-8320/89/$03.50 ~, 1989 Elsevicr Science Publishers Ltd, England. Printed in Great Britain
27~
Q
R" S Z I , Z2 . . . .
/I O , , O 2 .... /-k v P
(1)
I: " :[ iI'i I •
.4 ttila ('senki
Probability of failure Probability measure n-Dimensional Euclidean space Hypersphere Basic variables in the standard normal space Reliability index R a n d o m angles Measures R a n d o m radius Standard normal distribution function Density of the central chi-square distribution with n degrees of freedom Euclidean norm Maximum norm
Superscripts C
T
Complement Matrix transposition ! INTRODUCTION
Given a limit state function g(Z 1.... , Z.) in terms of the standard normal basic variables Z = {Z, ..... Z.) r the expression for the failure probability Pr associated with g can be written as (see for example, Ref. 1) Pf = P(~IZ) _< O)
{I)
An approximate value of the failure probability in eqn (1) is based on the reliability index fl which is defined as the distance of the failure surface
Iz R"fg(z) = 0} from the origin, fl can be calculated by an iterative method I and Pf is then approximated as Pf "- (1)~- lt) (2) where (l) denotes the standard normal distribution function. The approximation given by eqn (2) can be improved by subsequently employing a Monte Carlo method whereby the fl-sphere with centre at the origin, which of course is contained in the safe region, is excluded from samplingf The class of procedures considered here can in outline be described as follows. Take a (measurable) region B which contains the failure region
Iz R"l z) _< O}
B
An improved Monte Carlo method in structural reliability
277
Then eqn (1) can be rewritten as Pf = P(g(Z) _< 0, Z e B ) = P(g(Z) ___0 [ Z e B ) P ( Z E B )
(3)
The first term on the r.h.s, of eqn(3) denotes the conditional failure probability given that ZeB. It will be estimated by employing a suitable Monte Carlo procedure while the second term in eqn (3) is to be calculated numerically. The method resulting if B is assumed to be the whole ndimensional space is the crude Monte Carlo method. 3 If B is assumed to be the complement of the fl-sphere the method is as given in Ref. 2. In general, N
~6N = N - l ~ ' l~z,)~o I P(Z~ B)
(4)
i=1
is an unbiased estimator of Pf if z~ ..... z N are N independent observations from the same n-dimensional parent distribution Q defined by Q(A) = P(Z~A IZeB),
A ~ R"
(5)
(I~ I denotes the indicator function of the set {...}). The variance of the estimator ,fin is given by [ p ( Z e B ) ] 2 N - 2N Varo(l:~zl_
N-
ef[e(z
a) - er]
Hence, the coefficient of variation of the estimator fin is given by COVe(/~N) = N-1/2 [P(ZeB) - Pf] ~/2Pr- 1/2
(6)
Equation (6) shows that in order to achieve a small coefficient of variation of/~.~ the region B has to be selected such that it approximates the failure region from above as closely as possible. Furthermore, in order to calculate the estimator PN from eqn (4) we have to restrict ourselves to domains B the probability content of which can be evaluated numerically. Finally, a method has to be available for obtaining independent samples according to the distribution Q. This paper presents the theory for regions B which are complements of hyperspheres with the centre not necessarily at the co-ordinate origin. This entails seeking a hypersphere S = / V in the safe region with maximum probability content (or near so). Subsequent sampling from Q as given by eqn (5) yields a value for/~N. It appears that such an optimization effort can result in a further order of magnitude reduction of the number N of limit state function evaluations when compared with Ref. 2 while maintaining the same coefficient of variation.
,4tti& ('senki
278
It should be noted that the method presented in this paper is an importance sampling technique 3 whereby the crude Monte Carlo estimator, i.e. the indicator function of the failure region is replaced by the ratio of that function and I a which then after normalization is a new estimator mimicking the r a n d o m behaviour of the old one. (Since B contains the failure region the above ratio prior to normalization coincides with the old estimator.) The optimization procedure is intended to minimize the random variation of the new estimator, while the sampling procedure allows independent samples from Q to be generated. The paper is concluded with an example in which results by the present method are compared with those obtained when Harbitz's method 2 is employed.
2 THE OPTIMIZATION PROCEDURE An algorithm for finding a hypersphere contained in the safe region with (near) m a x i m u m probability content can be formulated by means of one of the well-known constrained optimization methods, e.g. the Method of Steepest Descent. 4 The problem then reduces to formulating subroutines which evaluate the objective function and the constraint(s). The situation is illustrated in Fig. 1 for n = 2.
2.1 The objective function f For a given point a = (a I ..... an) r in the safe region R(a) denotes the distance of a from the failure surface and S(a; r) denotes the hypersphere with centre a and radius r R(a) = min {l[z - a[! Ig(z) = 0}
Sta;r)= {z~R" I IIz-all _
(7)
(8)
The objective function f evaluated at a is then simply given by / ( a ) = P(Z~S(a; R(a)))
(9)
R(a) is easily obtained by applying a usual level 2 reliability algorithm t to the 'new' limit state function h defined by h(w) = g(w + a) and by writing the point z* at which the minimum is attained in eqn (7) as
z* -- w* + a where w* is the design point associated with the limit state function h.
An improved Monte Carlo method in structural reliability
279
Fig. 1. The maximum hypersphere. The standard normal probability content of S(a;r) is equal to the distribution function evaluated at r 2 of a non-central chi-square distribution with n degrees of freedom and non-centrality parameter IIa II2. Even though methods are available to compute the distribution function of a quadratic form of a standard normal random vector 5-7 there seems to be no algorithm which is both elementary and does not use truncation of an infinite series of functions at some stage of the calculation. A simple, straightforward method of evaluation of the non-central chi-square distribution function is provided in Appendix 1. It is based on the introduction of polar co-ordinates and subsequent reduction of the dimension by means of a recurrence relation. The resulting one-dimensional integral is in terms of the standard normal distribution function • for which any of the usual high-precision numerical approximations a'9 may be taken to perform a numerical integration. 2.2 The constraints The objective function f(a) defined by eqn (9) should ideally be maximized over the entire safe region, i.e. the inequality constraint g(a) > 0
(10)
280
Attila Csenki
should be applied. In most cases, however, no maximum will exist over the whole safe region and therefore eqn (10) has to be replaced by a constraint confining a to a compact subset of the safe region over which the existence of a maximum follows from the fairly weak assumption o f a c o n t i n u o u s f The overall optimization procedure may be carried out as a sequence of optimization steps with the current region being defined as the optimum sphere from the previous step. The region for the first step may be defined as the//-sphere. The optimization procedure can be terminated if in any one step a local maximum is found, i.e. the point obtained is inside the region being considered, or after a predetermined number of steps. Using the above approach involves inequality constraints the methods for which are computationally time-consuming and not always certain to converge. Therefore, it is preferable to use range constraints. Suitable range constraints are defined by any hypercube contained in the previous optimum sphere S(a; r). Such a hypercube is given by, for example,
[z~R"l!lz-ail ~ <0"7rn-1"2} If available, library optimization routines such as the NAG 1o may be used.
3 THE SAMPLING PROCEDURE Having obtained a (near) optimum hypersphere S(a; R) in the safe region it remains to estimate the first term on the r.h.s, of eqn (3) by a Monte Carlo method. Since S will in general not be centred at the origin the conditional random radius and direction are not independent and thus the sampling procedure employed by Harbitz 2 is not applicable. In this section a composition method sampling procedure 3 is described which first delivers a sample t 1/2 of the random radius p and subsequently yields a sample (0~ ..... 0,_ t) of the angles (O~,..., O,_ ~) representing the random direction. The latter sampling is conditioned on p2 = t and both sampling steps are conditioned on the event that the optimum sphere is excluded.
3.1 Change of co-ordinates prior to sampling New co-ordinates will be introduced for which the centre of S lies on the positive section of the last co-ordinate axis (see Fig. 1). The objective of the sampling procedure is to estimate the first term on the r.h.s, of eqn (3), that is e(g(Z) < OI l i t - a l l > R)
(11)
An improved Monte Carlo method in structural reliability
281
In Appendix 2 an orthogonal matrix M is constructed such that M r a = (0,...,0, Ilall) x
(12)
and hence the probability in eqn (11) can be reformulated to give
a)ll > R)
P(g(MMrZ) < 01 IIMr(Z -
= P ( g l ( V ) < 0 l IIV - ( 0 , . . . , 0 , Ilall)rll > R)
(13)
where gx = g o M is the g function written in terms of the v co-ordinates with V = MTZ being a standard normal n-dimensional random vector. Sampling will take place in the v-space conditioned on the event that the (transformed) maximum hypersphere is excluded from sampling. 3.2 Sampling in the new space of variables
The random variable V can be represented in polar co-ordinates as n-I 11"-"1r
VI=p]
Ic°sOi & & i=l
n-j T--T
V~=psinO._~+l]
IcosO/,
j = 2 ..... n
(14)
i=1
with independent random variables O, ®1 ..... ® . - 1 such that pZ has a chisquare distribution with n degrees of freedom and the angles ®~ take values in [ - 7t/2, n/2] and [0, 2n] for i < n - I and i = n - 1 respectively. 11 Samples of V will be obtained from eqn (14) by first sampling p conditioned on the event as given in eqn(13) and then sampling the random angles ~)1 .... , O n - t . 3.2.1 Sampling the random radius p
Denoting the conditioning event in eqn (13) by E the (conditional) density of T = p2 is obtained from P(t < T < t + f t l E ) = P({t < T < l + f t } n E ) / P ( E ) = P(t < T < t + fit, T - - 2llall I/". + Ila[I2 > R2)/p(E)
= P(t < T < t + fit, T - 2Y~/Zllall sin®x + Ilall 2 > R2)/p(E) = Z2(t)L.(t;y, R ) f t / P ( E ) + o(f/), (fit ~,0)
(15)
282
A t t i & (5"enki
where
L , ( t ; y , R ) = P ( t - 2 t ~ " Z y s i n O ~ + y 2 > R2),
.1'= Ilal! > 0
and ,~2 stands for the density of the central chi-square distribution with n degrees o f freedom. Notice that the independence of T a n d ®~ has been used to arrive at eqn. (15). Observing that the density of O~ is proportional to COS"- 2 0 (Ref. 11) we get
L.(t; y, R) =
0
if0_
andR>y
1
if0
1
if t > t 2
0.5 + Oo/rt
if n = 2,
t t < t < 12
if n>_3,
t x < t <_ t 2
I and R < v (16)
0o
cn
cos"-20d0 d - hi2
with tt(2) = 0 ' ( + ) R ) 2
0 o = arcsin (max { - 1,rain { I,(}})
=(t + y2-
(17)
R2)y - ~t-~12/2
(18)
and a normalizing constant c, whose inverse is given by rt.2
c.- i =
S
cos" - 2 0 dO 1 3 5...(n-3) 24...(n-21
~
ifn>3iseven
=
(19) 24...(n-3) 2 ,1 3 5 . . . 1 n - 2 )
if n _> 3 is odd
The integral in eqn (16) together with the above formula (19) for c, may be obtained from the following recurrence relation (n = 2, 3 .... )
r' d
cos" 0 dO =
" n..2
[
sin 0o cos"- 1 0o + (n - 1)
f,o
cos"- 2 0 dO
- n/2
]/. s
Three cases will be distinguished when sampling from the distribution o f T given E whose distribution function will be denoted by G.
Case 1: y = 0. The maximum hypersphere is the fl-sphere and G is given by G(t) = {0 [H.(t) - H.(RZ)]/[1 - H.(R2)]
if 0 < t < R 2 if t > R 2
An improved Monte Carlo method in structural reliability
283
where /4. stands for the central chi-square distribution function with n degrees of freedom. Samples o f T c a n be obtained by the inverse method by observing that G - l(u) = H.-- l(u + (1 - u)H,(R2)),
0 < u< 1
C a s e 2: O < y < R. The co-ordinate origin lies in the m a x i m u m hypersphere but does not coincide with its centre. Using eqns (15) and (16) it can be seen that
P ( T > t2 I E) = [1 - H . ( t 2 ) ] / P ( E ) = 7
(20)
say, and P(tt < T < t 2 1 E ) =
1 -~
The distribution function G can be represented as a convex combination o f distribution functions G 1 and G2 as follows G(t) = ~Gt(t) + (1 - ~)G2(t)
(21)
where G ~(t) =
G2(t) =
I
I
- ~ 1,2,~,)(s)gZ.(s)/['? P( E)] ds
(22)
_ I [,,.,21(s)7~.2(s)L.(s; y, R)/[(1 - "~,)P(E)] ds ~t2
Since 7 is known from eqn (20) samples according to G are obtained from independent samples ~ and ~2 drawn from G~ and G 2 respectively when they are selected independently with probability 7 and (1 -~,) respectively, i.e. if q is the realization of an independent B(1,~)-variable then the distribution function of I'1~ 1 --~ (1 -- Y])~2
is G. ~1 can be generated by the inverse method by observing that Gl(t ) = [ ( H . ( t ) - H.(t2))/(1 - H.(t2))]Ii,2.~:)(t)
and hence G~-l(u) = H.-- l(u + (1 - u)H.(t2) ) The sample ¢2 is obtained by the rejection technique 3 as follows. Independent samples u I and u2 are drawn from uniform distributions on [tl,t2] and [0,m] respectively, where m is defined by m = max {X2(t)L.(t;y, R) it I < t < / 2 }
Attila C~'enki
284
If 2 u2 < ~.(u~)L.(u,; y, R)
then ~2 = ul is accepted as a sample from G 2 but sampling of ul and u 2 is repeated otherwise. The acceptance ratio C of this sampling procedure is given by C=
Z~(t)L.(t;y, R I d t / [ m ( t 2 - tl)] f"~I
= (1 - "/)PIE)/(4yRm) = (P(E) + H.(t2) - 1)/(4yRm)
The expected number K o f independent rectangular samples necessary in order to generate one sample o f T is K = I + 7 + 2(1 - ),)/C = 1 + (1 + 8 y R m - H.(t2))/P(E) Case 3: y > R. The m a x i m u m hypersphere does not contain the coordinate origin. This unusual case will occur very seldom and it signifies that a level 2 reliability method is not an appropriate one to be used. Nevertheless, the present Monte Carlo method might be useful to improve upon the level 2 result considerably. The sampling procedure is similar to the previous case and the formulae have to be modified as follows. Equation (21 ) is valid with 7 = P(T<_ I l or T > t2[E)
= (1 + H.(t~) - H . ( t z ) ) / P ( E ) Gl(t) =
It0.,, J~q,2..~,(s)z2(s)/[7 P(E)] ds ...j.
and G2 is again defined by eqn (22). ~l will be generated by the inverse method by observing that if t < ta
H.(t)
]+ H.(tl)Gl(t) =,
H.(t2)
H.(tl) 1 + H.(tt)-
H.(t2)
H.(tz) + H . ( t l ) 1 + H . ( t , ) - H.It2)
H.(t)-
if t 1 < t_< t2 if t > t 2
and hence G~- l(u) = {H.- t(u[1 + H.(ta) - H.(tz)]) H.-~(u + (1 + u ) [ H . ( t 2) - H.(tl)])
if 0 < u < u o if Uo < u < 1
An improved Monte Carlo method in structural reliability
285
where
Uo = H.(tl)/[l + H . ( t l ) - H.(t2)] Sampling the distribution Gz is identical to the previous case except that the acceptance ratio C now becomes
C = (P(E) + H.(tz) - H.(tl) - 1)/(4yRm) The expected number of independent rectangular samples needed to obtain one sample of T is now K = 1 + (1 + H.(t 1) + 8yRm - H.(t2))/P(E) C and K are measures of the effectiveness of the above sampling procedure. Because of the rather complicated structure of the density of G 2 it cannot be expected to find a transformation in a closed form of the random radius which would be analogous to the one that significantly increases the number of successful samples when the /J-sphere is excluded from samplingfl Nevertheless, the limit state function need not be calculated at the sampling stage, and therefore in cases where it is time-consuming to evaluate g the total computing time will still be dominated by the time spent in the limit state calculations even though a certain proportion of the random numbers drawn to generate samples of T will have to be discarded. The present procedure has the advantage of a reduced average number o f limit state function evaluations as implied by eqn (6).
3.2.2 Sampling the angles (91 ..... ®, - t Let P r =,(19 ~.IE)
(23)
denote the conditional distribution of 19 given T = t conditioned on the event E. For a given measurable subset A of the (n - l)-dimensional space the above conditional probability is thus obtained as the R a d o n - N i k o d y m derivative of a measure v with respect to a probability measure p (Ref. 12}, i.e.
v(D) = f o
PT='(O~A ] E)p(dt)
(241
where v and ,u are defined for the Borel subsets of the real line by
v(D) = P(TED, O E A JE)
(25)
p(D) = P(T6DI E) Assuming A to be an (n - l)-fold Cartesian product A = AI × "'" x A._, and
.Itttla ('~'enLi
286
using the conditional density o f Tgiven E ( e q n (15)) and the {unconditional) independence o f T, ®1 ..... ®.._ ,, eqn (25) m a y be evaluated as follows v(O)=
r e o } ~ ~ .~®ieAij.t~{ T - 2Tl"2vsin®~ + y2 > R z}
..:p(E)
i-: I n.- 1
I-] P(®i6Ai) f, _ ~=2 P(E) Jt~ P(®leAl"t n
2tl'2v® 1sin "
+ "v2 >
RE)p(T6dt)
1
=-([ I P(~)'~AiJ) ft) P(OIEAI'I-2tI'2")'Sin~)I )'. R) ............. +)'2> i-2
× (L.(t; y,
R)Z~.(t)/P(E))dt n--I
fD P(~)l ~A~,_t - 2tl'2)'sin ~)l_+ y2 > R2) ~ = --L.(t;y, R) P(®i6Ai)P(Te d t l E ) i.'2
(26) C o m p a r i n g eqns (24) and (26) shows that the conditional probability in eqn (23) is a p r o d u c t measure whose marginal distributions coincide with the respective (unconditional) marginal distributions of O except for the first one which has the distribution function
F.(s; t, 3'. R) =
P(® 1 -< s I t - 2t 1,2 v sin ® i + v2 > R2)
(27)
Several cases have to be distinguished when sampling from F.(.; t,y, R). The case n = 2 will be examined first, l f y --=-0 then, o f course, R > y and hence for L. to be non-zero t > t~ ( = t 2 = R 2) m a y be assumed (see eqn (16)). Thus, the conditioning event in eqn (27) is the full event and therefore F 2 is the unconditional distribution function of ®1 which is uniformly distributed on [0,2n]. If v > 0 then the conditioning event in e q n ( 2 7 ) i s equal to {sinO~ < if} with ~ defined by eqn (18). ~ < - 1 cannot occur since it would imply t < t~ and R > v and hence from eqn (16) L 2 = 0. The conditional distribution o f ®1 is therefore uniform on a set Ao with ['[0. 2~z] A0=~[0, arcsin~]w[~z-arcsin~,27t ] [ I x - arcsin ~. 2~ + arcsin C]
ifC > 1 if0<~
/ f n > 3 then the unconditional distribution o f (-), is concentrated on [ - ~/2, x/2] with density c. cos" - 2 0 (see eqn (19)). y = 0 again implies that F,, ~s equal to the unconditional distribution function o f O~. Equally, .v > 0 and
An improved Monte Carlo method in structural reliability
287
> 1 implies the same, while ~ < - 1 can be excluded if y > 0 for reasons already explained above. Hence, for y > 0 and s between - r e / 2 and 0o (defined by eqn (17)) we have
F.(s; t,y, R) = P(®I <- s)/P(O1 < 0o) =
cos"- 2 0 dO -
/foo
n/2
-
cos"- 2 0 dO
n/2
Samples from/7, may easily be obtained by the inverse method by applying the N e w t o n - R a p h s o n scheme which will converge for the starting point sl = min {0,0o} since F. is convex for s < s~ and concave for s >_ s~. Finally, samples of ®i for n > 3 (i = 2 .... ,n - 1) are obtained as follows. O . _ ~ is uniformly distributed on [0, 2hi. O I (2 _
4 AN EXAMPLE A general F O R T R A N 77 program has been written and applied to an example in Ref. 1. A propped cantilever with a total (deterministic) length of 7-5 m is supported at a distance of 1= 5.0m from its built-in end. It is subjected to a G u m b e l distributed end point load p with mean 4.0 kN and standard deviation 1.0kN. Both Young's modulus E and the moment of inertia I are assumed to be normally distributed with means and standard
TABLE
1
Results of the O p t i m i z a t i o n Procedure
Maximum hypersphere
fl-sphere Centre Radius R Probability y
(0, 0, 0)1 3.3220 1 -1.1534 x 10 -2 0
],~
.--
C K m
• -
(-3"7366× 4"748 2
101
1 1"6101 × 10
1'4351 1'566 5 1"974 2 1-1131 2"9922
× × × ×
10-5 10- 1 10 ~ 10 -4
3
1"3426,3"4269x 10-1)T
288
Attila ( :~enki TABLE 2
Results of the Sampling Procedure Number o[ g-eraluations N
101 5 x 101 102 5 × 102 103 5 × 10"~ 104 5 × 104 10-~
Number t)f samph's in .[ailure region . . . . . . . . . . . . . . . . . . . . . . .
fl-sphere
Max. sphere
0 2 3 23 55 277 561 2907 5888
2 20 41 211 414 2082 4 215 20955 41 851
I(# j6N fl-sphere
Max. sphere
Samph, value of K (Max. sphere)
0"0000 4'6134 3"460 1 5'3054 6'343 4 6"3900 6"4703 6"7056 6"7909
3"2202 6"4406 6-601 6 6"7948 6.6660 6.704 7 6-786 8 6"748 1 6'7386
11.20 11.44 10"58 10"76 11.31 11.03 11-11 11.10 11.12
deviations 2 x 107kN/m 2, 5 × 106kN/m 2 and 10-'~m '~, 2 x 1 0 - S m 4 respectively. Failure is assumed to occur if the deflection criterion 5 p l a / ( 4 8 E l ) < I,/30
is violated. The three basic variables p , E and I are assumed to be independent. Results o f the analysis are shown in Tables 1 and 2. The level 2 failure probability (4-47 × 10 -4) is some 34% away from the actual failure probability which is 6.74 x 10 -4. Table 2 shows that in order to achieve the same accuracy in the fiN-estimate by excluding the fl-sphere from sampling as is obtained from,fiN based on the m a x i m u m hypersphere for, say, N = I00, the number o f limit state function evaluations has to be increased to well above 104 . The program for the above example uses the N A G program library 1° and is capable o f analysing problems with a wide range of assumptions. It is available from the a u t h o r upon request.
5 CONCLUSIONS A Monte Carlo method has been presented which is based on the exclusion from sampling o f a m a x i m u m hypersphere in the standard normal space of basic variables. The number o f limit state function evaluations needed to achieve a given accuracy of the estimate o f the failure probability by this method is orders of magnitude smaller than that required when the fl-sphere
An improved Monte Carlo method in structural reliability
289
is excluded from sampling. An example has been given demonstrating the merits of the technique.
REFERENCES 1. Thoft-Christensen, P. & Baker, M. J., Structural Reliability Theory and its Applications. Springer-Verlag, Berlin, Heidelberg, New York, 1982. 2. Harbitz, A., An efficient sampling method for probability of failure calculation. Structural Safety, 3 (1986) 109-15. 3. Hammersley, J. M. & Handscomb, D. C., Monte Carlo Methods. Chapman and Hall, London, 1979. 4. Siddal, J. N., Optimal Engineering Design, Marcel Dekker, New York, Basel, 1982. 5. Rice, S. O., Distribution of quadratic forms in normal random variables-evaluation by numerical integration. S l A M J. Sci. Stat. Comput., 1 (1981) 438-48. 6. Fiessler, B., Neumann, H.-J. & Rackwitz, R., Quadratic limit states in structural reliability. J. Engng Mech. Div., ASCE, Proc. Paper 14739, 105(EM4) (1979) 661-76. 7. Narula, S. C. & Desu, M. M., Computation of probability and non-centrality parameter of a non-central chi-squared distribution, algorithm AS170. Appl. Statist., 30 (1981) 349-52. 8. Rosenblueth, E., On computing normal reliabilities. Structural Safety, 2 (1985) 165-7; Errata in Structural Safety, 3 (1985) 67. 9. Hart, J. F., Cheney, E. W., Lawson, C. L., Maehly, H. J., Mesztenyi, C. K., Rice, J. R., Thacher, H. G. & Witzgall, C., Computer Approximations. Robert E. Krieger Publ. Co., Huntington and New York, 1978. 10. Numerical Algorithms Group (NAG), FORTRAN library manual, Mark 11. Oxford, 1984. 11. Kendall, M. G., A Course in the Geometry of n Dimensions. Gri~n's Statistical Monographs & Courses, Charles Griffin, London, 1961. 12. R6nyi, A., Probability Theory. North-Holland, Amsterdam, London, 1970. 13. Abramowitz, M. & Stegun, I. A. (eds), Handbook of Mathematical Functions. Dover, New York, 1965. 14. Leon, S. J., Linear Algebra with Applications. Macmillan, New York, 1980.
A P P E N D I X 1: T H E S T A N D A R D N O R M A L P R O B A B I L I T Y C O N T E N T O F A N O F F S E T H Y P E R S P H E R E IN n > 2 DIMENSIONS The standard normal probability content o f a hypersphere S(a; R) defined by eqn (8) can be written as P(ZES(a; R)) = P(Z~S(c; R))
(28)
290
.4ttila (.~enki
where c denotes the n-dimensional vector whose only non-zero entry is the last c o m p o n e n t which is equal to the Euclidean norm ofa. The explicit form o f a suitable orthogonal transformation justifying eqn (28)is immaterial for our present purposes but it is nevertheless given in Appendix 2 because it will be used in the sampling procedure. Equation (28) simply can be viewed as a consequence o f the invariance o f the standard normal density under orthogonal transformations. The above probability thus can be written as
P(Z~S(e; R))= f "" f (2n)-"'2 exp,- f'u + ellZ/2)d(u, ..... u.) lul <_ R
d
d
'"' < ~
d(ut ..... u.)
=
... )
-n,'2
In
(2r0-";2 exp [-(r
2 + 2r taii s i n Ot
d-n,'2dO
2) times
Jr-Ilal[2)/2]
n 2 ×
r"-t [ [cos.-l-i0id0._l...d01dr
(29)
i~l
where the last step is by transition to polar co-ordinates. ~t Notice that the domains o f integration for 01 in eqn (29) are from 0 to 2rt and from - rr,.'2 to 7r/2 for n = 2 and n _> 3 respectively. Equation (29) thus can be rewritten as
~:.:~.
cos"- 20J.(O) dO
(30)
where {0 xo =
if n = 2 rt"2
if n > 3
xt =
2n rr/2
if n = 2 if n_>3
e.=
1 2rt
ifn = 2 ifn>_3
n--2
:c. -- (2rr) '"'2
[1F
cos"- i -iOd 0
-~.2
i=2
J.(O) = I "r" i exp [ 3o
(r 2 +
2r [ia I1sin 0 + tl a !12)/'2] dr
(31)
An improved Monte Carlo method in structural reliability
291
(For n = 2, 3 the product in eqn (31) is understood to be unity.) It can be shown that (2n)-t ~.0t. =
if n = 2
2 ~n/2~- I n - 1
n - 2)!
-- 1
2'-"/2rt
-
i f n > 4 is e v e n
if n is o d d
3
I \-2--/
Finally, integration by parts shows that the following recurrence relation holds true for J.(O) J2(O) = ao + bo + CoJl(O) J.(O) = (n - 2)J. _ 2(0) + CoJ. - 1(0) + bo Rn - 2
n> 3
(32)
with ao = exp ( - IIa II2/2) bo = - e x p [ - ( R 2 +
2RIlall
sin0+
Ila 112)/2-1
Co = - II a 11 sin 0
In order to evaluate Jn(O) from eqn (32) only J,(O) needs to be known which can be written as Jl(O) = (2n) 1'2 exp [ -(tl a i12 cos 2 0)/2] x ((1)(R + !',all s i n 0 ) - ~(liall sin 0))
(33)
Taking in eqn (33) a high-precision approximation for the standard normal distribution function tl)s'° enables us to obtain J~ and hence J, by means of the recurrence relation eqn(32) without numerical integration. The calculation of the n-dimensional integral in eqn (28) is therefore reduced to a simple one-dimensional integration. It should be noted that a simpler approximation o f • which is usually sufficiently accurate for standard reliability calculations 13 proved inadequate for the purposes of this probability evaluation.
A P P E N D I X 2: C O N S T R U C T I O N OF T H E M A T R I X M An orthogonal matrix M satisfying eqn (12) is required. If a is non-zero (which may be assumed since otherwise the identity matrix is suitable) then the G r a m - S c h m i d t orthogonalization process t4 can be applied as follows. The basis to be orthogonalfzed is obtained from the standard basis of the ndimensional space by replacing its ith element by a with i t { 1..... n} chosen
202
.,! ttila ( ".~cnk i
such that a i ¢ O. D e n o t i n g the e l e m e n t s o f this basis by Pl . . . . . p._ 1,a the c o l u m n v e c t o r s o f M arc defined by the following r e c u r r e n c e relation m n --
mj :
a ii
! Ui
1 a,
i I Uj,
j = n --
where
uj=pi- ~ (p~m,)mk k-j~]
1, n - - 2 . . . . .
1
An Improved Monte Carlo Method in Structural Reliability
Attila Csenki Centre for Software Reliability, The City University, Northampton Square, London ECIV 0HB, UK (Received 15 July 1988; accepted 20 July 1988)
A BSTRA C T A Monte Carlo technique is presented for obtaining accurate estimates of component failure probabilities based on small sample sizes. It is an improvement on a recent method by Harbitz and is based on the exclusion from sampling of a maximum hypersphere contained in the safe region as opposed to that of the fl-sphere. After finding the maximum hypersphere in the standard normal space of basic variables a two-stage sampling procedure involving conditional distributions is employed.
NOTATION 8
g, gt, h tl W*, Z* B
C
F,G U. M
Centre of hypersphere Limit state functions Number of basic variables Estimate of failure probability Design points Sampling region Acceptance ratio Distribution functions Distribution function of the central chi-square distribution with n degrees of freedom Indicator function of the set {...} Orthogonal n x n matrix
275 Reliability Engineering and System Safety 0951-8320/89/$03.50 ~, 1989 Elsevicr Science Publishers Ltd, England. Printed in Great Britain
27~
Q
R" S Z I , Z2 . . . .
/I O , , O 2 .... /-k v P
(1)
I: " :[ iI'i I •
.4 ttila ('senki
Probability of failure Probability measure n-Dimensional Euclidean space Hypersphere Basic variables in the standard normal space Reliability index R a n d o m angles Measures R a n d o m radius Standard normal distribution function Density of the central chi-square distribution with n degrees of freedom Euclidean norm Maximum norm
Superscripts C
T
Complement Matrix transposition ! INTRODUCTION
Given a limit state function g(Z 1.... , Z.) in terms of the standard normal basic variables Z = {Z, ..... Z.) r the expression for the failure probability Pr associated with g can be written as (see for example, Ref. 1) Pf = P(~IZ) _< O)
{I)
An approximate value of the failure probability in eqn (1) is based on the reliability index fl which is defined as the distance of the failure surface
Iz R"fg(z) = 0} from the origin, fl can be calculated by an iterative method I and Pf is then approximated as Pf "- (1)~- lt) (2) where (l) denotes the standard normal distribution function. The approximation given by eqn (2) can be improved by subsequently employing a Monte Carlo method whereby the fl-sphere with centre at the origin, which of course is contained in the safe region, is excluded from samplingf The class of procedures considered here can in outline be described as follows. Take a (measurable) region B which contains the failure region
Iz R"l z) _< O}
B
An improved Monte Carlo method in structural reliability
277
Then eqn (1) can be rewritten as Pf = P(g(Z) _< 0, Z e B ) = P(g(Z) ___0 [ Z e B ) P ( Z E B )
(3)
The first term on the r.h.s, of eqn(3) denotes the conditional failure probability given that ZeB. It will be estimated by employing a suitable Monte Carlo procedure while the second term in eqn (3) is to be calculated numerically. The method resulting if B is assumed to be the whole ndimensional space is the crude Monte Carlo method. 3 If B is assumed to be the complement of the fl-sphere the method is as given in Ref. 2. In general, N
~6N = N - l ~ ' l~z,)~o I P(Z~ B)
(4)
i=1
is an unbiased estimator of Pf if z~ ..... z N are N independent observations from the same n-dimensional parent distribution Q defined by Q(A) = P(Z~A IZeB),
A ~ R"
(5)
(I~ I denotes the indicator function of the set {...}). The variance of the estimator ,fin is given by [ p ( Z e B ) ] 2 N - 2N Varo(l:~zl_
N-
ef[e(z
a) - er]
Hence, the coefficient of variation of the estimator fin is given by COVe(/~N) = N-1/2 [P(ZeB) - Pf] ~/2Pr- 1/2
(6)
Equation (6) shows that in order to achieve a small coefficient of variation of/~.~ the region B has to be selected such that it approximates the failure region from above as closely as possible. Furthermore, in order to calculate the estimator PN from eqn (4) we have to restrict ourselves to domains B the probability content of which can be evaluated numerically. Finally, a method has to be available for obtaining independent samples according to the distribution Q. This paper presents the theory for regions B which are complements of hyperspheres with the centre not necessarily at the co-ordinate origin. This entails seeking a hypersphere S = / V in the safe region with maximum probability content (or near so). Subsequent sampling from Q as given by eqn (5) yields a value for/~N. It appears that such an optimization effort can result in a further order of magnitude reduction of the number N of limit state function evaluations when compared with Ref. 2 while maintaining the same coefficient of variation.
,4tti& ('senki
278
It should be noted that the method presented in this paper is an importance sampling technique 3 whereby the crude Monte Carlo estimator, i.e. the indicator function of the failure region is replaced by the ratio of that function and I a which then after normalization is a new estimator mimicking the r a n d o m behaviour of the old one. (Since B contains the failure region the above ratio prior to normalization coincides with the old estimator.) The optimization procedure is intended to minimize the random variation of the new estimator, while the sampling procedure allows independent samples from Q to be generated. The paper is concluded with an example in which results by the present method are compared with those obtained when Harbitz's method 2 is employed.
2 THE OPTIMIZATION PROCEDURE An algorithm for finding a hypersphere contained in the safe region with (near) m a x i m u m probability content can be formulated by means of one of the well-known constrained optimization methods, e.g. the Method of Steepest Descent. 4 The problem then reduces to formulating subroutines which evaluate the objective function and the constraint(s). The situation is illustrated in Fig. 1 for n = 2.
2.1 The objective function f For a given point a = (a I ..... an) r in the safe region R(a) denotes the distance of a from the failure surface and S(a; r) denotes the hypersphere with centre a and radius r R(a) = min {l[z - a[! Ig(z) = 0}
Sta;r)= {z~R" I IIz-all _
(7)
(8)
The objective function f evaluated at a is then simply given by / ( a ) = P(Z~S(a; R(a)))
(9)
R(a) is easily obtained by applying a usual level 2 reliability algorithm t to the 'new' limit state function h defined by h(w) = g(w + a) and by writing the point z* at which the minimum is attained in eqn (7) as
z* -- w* + a where w* is the design point associated with the limit state function h.
An improved Monte Carlo method in structural reliability
279
Fig. 1. The maximum hypersphere. The standard normal probability content of S(a;r) is equal to the distribution function evaluated at r 2 of a non-central chi-square distribution with n degrees of freedom and non-centrality parameter IIa II2. Even though methods are available to compute the distribution function of a quadratic form of a standard normal random vector 5-7 there seems to be no algorithm which is both elementary and does not use truncation of an infinite series of functions at some stage of the calculation. A simple, straightforward method of evaluation of the non-central chi-square distribution function is provided in Appendix 1. It is based on the introduction of polar co-ordinates and subsequent reduction of the dimension by means of a recurrence relation. The resulting one-dimensional integral is in terms of the standard normal distribution function • for which any of the usual high-precision numerical approximations a'9 may be taken to perform a numerical integration. 2.2 The constraints The objective function f(a) defined by eqn (9) should ideally be maximized over the entire safe region, i.e. the inequality constraint g(a) > 0
(10)
280
Attila Csenki
should be applied. In most cases, however, no maximum will exist over the whole safe region and therefore eqn (10) has to be replaced by a constraint confining a to a compact subset of the safe region over which the existence of a maximum follows from the fairly weak assumption o f a c o n t i n u o u s f The overall optimization procedure may be carried out as a sequence of optimization steps with the current region being defined as the optimum sphere from the previous step. The region for the first step may be defined as the//-sphere. The optimization procedure can be terminated if in any one step a local maximum is found, i.e. the point obtained is inside the region being considered, or after a predetermined number of steps. Using the above approach involves inequality constraints the methods for which are computationally time-consuming and not always certain to converge. Therefore, it is preferable to use range constraints. Suitable range constraints are defined by any hypercube contained in the previous optimum sphere S(a; r). Such a hypercube is given by, for example,
[z~R"l!lz-ail ~ <0"7rn-1"2} If available, library optimization routines such as the NAG 1o may be used.
3 THE SAMPLING PROCEDURE Having obtained a (near) optimum hypersphere S(a; R) in the safe region it remains to estimate the first term on the r.h.s, of eqn (3) by a Monte Carlo method. Since S will in general not be centred at the origin the conditional random radius and direction are not independent and thus the sampling procedure employed by Harbitz 2 is not applicable. In this section a composition method sampling procedure 3 is described which first delivers a sample t 1/2 of the random radius p and subsequently yields a sample (0~ ..... 0,_ t) of the angles (O~,..., O,_ ~) representing the random direction. The latter sampling is conditioned on p2 = t and both sampling steps are conditioned on the event that the optimum sphere is excluded.
3.1 Change of co-ordinates prior to sampling New co-ordinates will be introduced for which the centre of S lies on the positive section of the last co-ordinate axis (see Fig. 1). The objective of the sampling procedure is to estimate the first term on the r.h.s, of eqn (3), that is e(g(Z) < OI l i t - a l l > R)
(11)
An improved Monte Carlo method in structural reliability
281
In Appendix 2 an orthogonal matrix M is constructed such that M r a = (0,...,0, Ilall) x
(12)
and hence the probability in eqn (11) can be reformulated to give
a)ll > R)
P(g(MMrZ) < 01 IIMr(Z -
= P ( g l ( V ) < 0 l IIV - ( 0 , . . . , 0 , Ilall)rll > R)
(13)
where gx = g o M is the g function written in terms of the v co-ordinates with V = MTZ being a standard normal n-dimensional random vector. Sampling will take place in the v-space conditioned on the event that the (transformed) maximum hypersphere is excluded from sampling. 3.2 Sampling in the new space of variables
The random variable V can be represented in polar co-ordinates as n-I 11"-"1r
VI=p]
Ic°sOi & & i=l
n-j T--T
V~=psinO._~+l]
IcosO/,
j = 2 ..... n
(14)
i=1
with independent random variables O, ®1 ..... ® . - 1 such that pZ has a chisquare distribution with n degrees of freedom and the angles ®~ take values in [ - 7t/2, n/2] and [0, 2n] for i < n - I and i = n - 1 respectively. 11 Samples of V will be obtained from eqn (14) by first sampling p conditioned on the event as given in eqn(13) and then sampling the random angles ~)1 .... , O n - t . 3.2.1 Sampling the random radius p
Denoting the conditioning event in eqn (13) by E the (conditional) density of T = p2 is obtained from P(t < T < t + f t l E ) = P({t < T < l + f t } n E ) / P ( E ) = P(t < T < t + fit, T - - 2llall I/". + Ila[I2 > R2)/p(E)
= P(t < T < t + fit, T - 2Y~/Zllall sin®x + Ilall 2 > R2)/p(E) = Z2(t)L.(t;y, R ) f t / P ( E ) + o(f/), (fit ~,0)
(15)
282
A t t i & (5"enki
where
L , ( t ; y , R ) = P ( t - 2 t ~ " Z y s i n O ~ + y 2 > R2),
.1'= Ilal! > 0
and ,~2 stands for the density of the central chi-square distribution with n degrees o f freedom. Notice that the independence of T a n d ®~ has been used to arrive at eqn. (15). Observing that the density of O~ is proportional to COS"- 2 0 (Ref. 11) we get
L.(t; y, R) =
0
if0_
andR>y
1
if0
1
if t > t 2
0.5 + Oo/rt
if n = 2,
t t < t < 12
if n>_3,
t x < t <_ t 2
I and R < v (16)
0o
cn
cos"-20d0 d - hi2
with tt(2) = 0 ' ( + ) R ) 2
0 o = arcsin (max { - 1,rain { I,(}})
=(t + y2-
(17)
R2)y - ~t-~12/2
(18)
and a normalizing constant c, whose inverse is given by rt.2
c.- i =
S
cos" - 2 0 dO 1 3 5...(n-3) 24...(n-21
~
ifn>3iseven
=
(19) 24...(n-3) 2 ,1 3 5 . . . 1 n - 2 )
if n _> 3 is odd
The integral in eqn (16) together with the above formula (19) for c, may be obtained from the following recurrence relation (n = 2, 3 .... )
r' d
cos" 0 dO =
" n..2
[
sin 0o cos"- 1 0o + (n - 1)
f,o
cos"- 2 0 dO
- n/2
]/. s
Three cases will be distinguished when sampling from the distribution o f T given E whose distribution function will be denoted by G.
Case 1: y = 0. The maximum hypersphere is the fl-sphere and G is given by G(t) = {0 [H.(t) - H.(RZ)]/[1 - H.(R2)]
if 0 < t < R 2 if t > R 2
An improved Monte Carlo method in structural reliability
283
where /4. stands for the central chi-square distribution function with n degrees of freedom. Samples o f T c a n be obtained by the inverse method by observing that G - l(u) = H.-- l(u + (1 - u)H,(R2)),
0 < u< 1
C a s e 2: O < y < R. The co-ordinate origin lies in the m a x i m u m hypersphere but does not coincide with its centre. Using eqns (15) and (16) it can be seen that
P ( T > t2 I E) = [1 - H . ( t 2 ) ] / P ( E ) = 7
(20)
say, and P(tt < T < t 2 1 E ) =
1 -~
The distribution function G can be represented as a convex combination o f distribution functions G 1 and G2 as follows G(t) = ~Gt(t) + (1 - ~)G2(t)
(21)
where G ~(t) =
G2(t) =
I
I
- ~ 1,2,~,)(s)gZ.(s)/['? P( E)] ds
(22)
_ I [,,.,21(s)7~.2(s)L.(s; y, R)/[(1 - "~,)P(E)] ds ~t2
Since 7 is known from eqn (20) samples according to G are obtained from independent samples ~ and ~2 drawn from G~ and G 2 respectively when they are selected independently with probability 7 and (1 -~,) respectively, i.e. if q is the realization of an independent B(1,~)-variable then the distribution function of I'1~ 1 --~ (1 -- Y])~2
is G. ~1 can be generated by the inverse method by observing that Gl(t ) = [ ( H . ( t ) - H.(t2))/(1 - H.(t2))]Ii,2.~:)(t)
and hence G~-l(u) = H.-- l(u + (1 - u)H.(t2) ) The sample ¢2 is obtained by the rejection technique 3 as follows. Independent samples u I and u2 are drawn from uniform distributions on [tl,t2] and [0,m] respectively, where m is defined by m = max {X2(t)L.(t;y, R) it I < t < / 2 }
Attila C~'enki
284
If 2 u2 < ~.(u~)L.(u,; y, R)
then ~2 = ul is accepted as a sample from G 2 but sampling of ul and u 2 is repeated otherwise. The acceptance ratio C of this sampling procedure is given by C=
Z~(t)L.(t;y, R I d t / [ m ( t 2 - tl)] f"~I
= (1 - "/)PIE)/(4yRm) = (P(E) + H.(t2) - 1)/(4yRm)
The expected number K o f independent rectangular samples necessary in order to generate one sample o f T is K = I + 7 + 2(1 - ),)/C = 1 + (1 + 8 y R m - H.(t2))/P(E) Case 3: y > R. The m a x i m u m hypersphere does not contain the coordinate origin. This unusual case will occur very seldom and it signifies that a level 2 reliability method is not an appropriate one to be used. Nevertheless, the present Monte Carlo method might be useful to improve upon the level 2 result considerably. The sampling procedure is similar to the previous case and the formulae have to be modified as follows. Equation (21 ) is valid with 7 = P(T<_ I l or T > t2[E)
= (1 + H.(t~) - H . ( t z ) ) / P ( E ) Gl(t) =
It0.,, J~q,2..~,(s)z2(s)/[7 P(E)] ds ...j.
and G2 is again defined by eqn (22). ~l will be generated by the inverse method by observing that if t < ta
H.(t)
]+ H.(tl)Gl(t) =,
H.(t2)
H.(tl) 1 + H.(tt)-
H.(t2)
H.(tz) + H . ( t l ) 1 + H . ( t , ) - H.It2)
H.(t)-
if t 1 < t_< t2 if t > t 2
and hence G~- l(u) = {H.- t(u[1 + H.(ta) - H.(tz)]) H.-~(u + (1 + u ) [ H . ( t 2) - H.(tl)])
if 0 < u < u o if Uo < u < 1
An improved Monte Carlo method in structural reliability
285
where
Uo = H.(tl)/[l + H . ( t l ) - H.(t2)] Sampling the distribution Gz is identical to the previous case except that the acceptance ratio C now becomes
C = (P(E) + H.(tz) - H.(tl) - 1)/(4yRm) The expected number of independent rectangular samples needed to obtain one sample of T is now K = 1 + (1 + H.(t 1) + 8yRm - H.(t2))/P(E) C and K are measures of the effectiveness of the above sampling procedure. Because of the rather complicated structure of the density of G 2 it cannot be expected to find a transformation in a closed form of the random radius which would be analogous to the one that significantly increases the number of successful samples when the /J-sphere is excluded from samplingfl Nevertheless, the limit state function need not be calculated at the sampling stage, and therefore in cases where it is time-consuming to evaluate g the total computing time will still be dominated by the time spent in the limit state calculations even though a certain proportion of the random numbers drawn to generate samples of T will have to be discarded. The present procedure has the advantage of a reduced average number o f limit state function evaluations as implied by eqn (6).
3.2.2 Sampling the angles (91 ..... ®, - t Let P r =,(19 ~.IE)
(23)
denote the conditional distribution of 19 given T = t conditioned on the event E. For a given measurable subset A of the (n - l)-dimensional space the above conditional probability is thus obtained as the R a d o n - N i k o d y m derivative of a measure v with respect to a probability measure p (Ref. 12}, i.e.
v(D) = f o
PT='(O~A ] E)p(dt)
(241
where v and ,u are defined for the Borel subsets of the real line by
v(D) = P(TED, O E A JE)
(25)
p(D) = P(T6DI E) Assuming A to be an (n - l)-fold Cartesian product A = AI × "'" x A._, and
.Itttla ('~'enLi
286
using the conditional density o f Tgiven E ( e q n (15)) and the {unconditional) independence o f T, ®1 ..... ®.._ ,, eqn (25) m a y be evaluated as follows v(O)=
r e o } ~ ~ .~®ieAij.t~{ T - 2Tl"2vsin®~ + y2 > R z}
..:p(E)
i-: I n.- 1
I-] P(®i6Ai) f, _ ~=2 P(E) Jt~ P(®leAl"t n
2tl'2v® 1sin "
+ "v2 >
RE)p(T6dt)
1
=-([ I P(~)'~AiJ) ft) P(OIEAI'I-2tI'2")'Sin~)I )'. R) ............. +)'2> i-2
× (L.(t; y,
R)Z~.(t)/P(E))dt n--I
fD P(~)l ~A~,_t - 2tl'2)'sin ~)l_+ y2 > R2) ~ = --L.(t;y, R) P(®i6Ai)P(Te d t l E ) i.'2
(26) C o m p a r i n g eqns (24) and (26) shows that the conditional probability in eqn (23) is a p r o d u c t measure whose marginal distributions coincide with the respective (unconditional) marginal distributions of O except for the first one which has the distribution function
F.(s; t, 3'. R) =
P(® 1 -< s I t - 2t 1,2 v sin ® i + v2 > R2)
(27)
Several cases have to be distinguished when sampling from F.(.; t,y, R). The case n = 2 will be examined first, l f y --=-0 then, o f course, R > y and hence for L. to be non-zero t > t~ ( = t 2 = R 2) m a y be assumed (see eqn (16)). Thus, the conditioning event in eqn (27) is the full event and therefore F 2 is the unconditional distribution function of ®1 which is uniformly distributed on [0,2n]. If v > 0 then the conditioning event in e q n ( 2 7 ) i s equal to {sinO~ < if} with ~ defined by eqn (18). ~ < - 1 cannot occur since it would imply t < t~ and R > v and hence from eqn (16) L 2 = 0. The conditional distribution o f ®1 is therefore uniform on a set Ao with ['[0. 2~z] A0=~[0, arcsin~]w[~z-arcsin~,27t ] [ I x - arcsin ~. 2~ + arcsin C]
ifC > 1 if0<~
/ f n > 3 then the unconditional distribution o f (-), is concentrated on [ - ~/2, x/2] with density c. cos" - 2 0 (see eqn (19)). y = 0 again implies that F,, ~s equal to the unconditional distribution function o f O~. Equally, .v > 0 and
An improved Monte Carlo method in structural reliability
287
> 1 implies the same, while ~ < - 1 can be excluded if y > 0 for reasons already explained above. Hence, for y > 0 and s between - r e / 2 and 0o (defined by eqn (17)) we have
F.(s; t,y, R) = P(®I <- s)/P(O1 < 0o) =
cos"- 2 0 dO -
/foo
n/2
-
cos"- 2 0 dO
n/2
Samples from/7, may easily be obtained by the inverse method by applying the N e w t o n - R a p h s o n scheme which will converge for the starting point sl = min {0,0o} since F. is convex for s < s~ and concave for s >_ s~. Finally, samples of ®i for n > 3 (i = 2 .... ,n - 1) are obtained as follows. O . _ ~ is uniformly distributed on [0, 2hi. O I (2 _
4 AN EXAMPLE A general F O R T R A N 77 program has been written and applied to an example in Ref. 1. A propped cantilever with a total (deterministic) length of 7-5 m is supported at a distance of 1= 5.0m from its built-in end. It is subjected to a G u m b e l distributed end point load p with mean 4.0 kN and standard deviation 1.0kN. Both Young's modulus E and the moment of inertia I are assumed to be normally distributed with means and standard
TABLE
1
Results of the O p t i m i z a t i o n Procedure
Maximum hypersphere
fl-sphere Centre Radius R Probability y
(0, 0, 0)1 3.3220 1 -1.1534 x 10 -2 0
],~
.--
C K m
• -
(-3"7366× 4"748 2
101
1 1"6101 × 10
1'4351 1'566 5 1"974 2 1-1131 2"9922
× × × ×
10-5 10- 1 10 ~ 10 -4
3
1"3426,3"4269x 10-1)T
288
Attila ( :~enki TABLE 2
Results of the Sampling Procedure Number o[ g-eraluations N
101 5 x 101 102 5 × 102 103 5 × 10"~ 104 5 × 104 10-~
Number t)f samph's in .[ailure region . . . . . . . . . . . . . . . . . . . . . . .
fl-sphere
Max. sphere
0 2 3 23 55 277 561 2907 5888
2 20 41 211 414 2082 4 215 20955 41 851
I(# j6N fl-sphere
Max. sphere
Samph, value of K (Max. sphere)
0"0000 4'6134 3"460 1 5'3054 6'343 4 6"3900 6"4703 6"7056 6"7909
3"2202 6"4406 6-601 6 6"7948 6.6660 6.704 7 6-786 8 6"748 1 6'7386
11.20 11.44 10"58 10"76 11.31 11.03 11-11 11.10 11.12
deviations 2 x 107kN/m 2, 5 × 106kN/m 2 and 10-'~m '~, 2 x 1 0 - S m 4 respectively. Failure is assumed to occur if the deflection criterion 5 p l a / ( 4 8 E l ) < I,/30
is violated. The three basic variables p , E and I are assumed to be independent. Results o f the analysis are shown in Tables 1 and 2. The level 2 failure probability (4-47 × 10 -4) is some 34% away from the actual failure probability which is 6.74 x 10 -4. Table 2 shows that in order to achieve the same accuracy in the fiN-estimate by excluding the fl-sphere from sampling as is obtained from,fiN based on the m a x i m u m hypersphere for, say, N = I00, the number o f limit state function evaluations has to be increased to well above 104 . The program for the above example uses the N A G program library 1° and is capable o f analysing problems with a wide range of assumptions. It is available from the a u t h o r upon request.
5 CONCLUSIONS A Monte Carlo method has been presented which is based on the exclusion from sampling o f a m a x i m u m hypersphere in the standard normal space of basic variables. The number o f limit state function evaluations needed to achieve a given accuracy of the estimate o f the failure probability by this method is orders of magnitude smaller than that required when the fl-sphere
An improved Monte Carlo method in structural reliability
289
is excluded from sampling. An example has been given demonstrating the merits of the technique.
REFERENCES 1. Thoft-Christensen, P. & Baker, M. J., Structural Reliability Theory and its Applications. Springer-Verlag, Berlin, Heidelberg, New York, 1982. 2. Harbitz, A., An efficient sampling method for probability of failure calculation. Structural Safety, 3 (1986) 109-15. 3. Hammersley, J. M. & Handscomb, D. C., Monte Carlo Methods. Chapman and Hall, London, 1979. 4. Siddal, J. N., Optimal Engineering Design, Marcel Dekker, New York, Basel, 1982. 5. Rice, S. O., Distribution of quadratic forms in normal random variables-evaluation by numerical integration. S l A M J. Sci. Stat. Comput., 1 (1981) 438-48. 6. Fiessler, B., Neumann, H.-J. & Rackwitz, R., Quadratic limit states in structural reliability. J. Engng Mech. Div., ASCE, Proc. Paper 14739, 105(EM4) (1979) 661-76. 7. Narula, S. C. & Desu, M. M., Computation of probability and non-centrality parameter of a non-central chi-squared distribution, algorithm AS170. Appl. Statist., 30 (1981) 349-52. 8. Rosenblueth, E., On computing normal reliabilities. Structural Safety, 2 (1985) 165-7; Errata in Structural Safety, 3 (1985) 67. 9. Hart, J. F., Cheney, E. W., Lawson, C. L., Maehly, H. J., Mesztenyi, C. K., Rice, J. R., Thacher, H. G. & Witzgall, C., Computer Approximations. Robert E. Krieger Publ. Co., Huntington and New York, 1978. 10. Numerical Algorithms Group (NAG), FORTRAN library manual, Mark 11. Oxford, 1984. 11. Kendall, M. G., A Course in the Geometry of n Dimensions. Gri~n's Statistical Monographs & Courses, Charles Griffin, London, 1961. 12. R6nyi, A., Probability Theory. North-Holland, Amsterdam, London, 1970. 13. Abramowitz, M. & Stegun, I. A. (eds), Handbook of Mathematical Functions. Dover, New York, 1965. 14. Leon, S. J., Linear Algebra with Applications. Macmillan, New York, 1980.
A P P E N D I X 1: T H E S T A N D A R D N O R M A L P R O B A B I L I T Y C O N T E N T O F A N O F F S E T H Y P E R S P H E R E IN n > 2 DIMENSIONS The standard normal probability content o f a hypersphere S(a; R) defined by eqn (8) can be written as P(ZES(a; R)) = P(Z~S(c; R))
(28)
290
.4ttila (.~enki
where c denotes the n-dimensional vector whose only non-zero entry is the last c o m p o n e n t which is equal to the Euclidean norm ofa. The explicit form o f a suitable orthogonal transformation justifying eqn (28)is immaterial for our present purposes but it is nevertheless given in Appendix 2 because it will be used in the sampling procedure. Equation (28) simply can be viewed as a consequence o f the invariance o f the standard normal density under orthogonal transformations. The above probability thus can be written as
P(Z~S(e; R))= f "" f (2n)-"'2 exp,- f'u + ellZ/2)d(u, ..... u.) lul <_ R
d
d
'"' < ~
d(ut ..... u.)
=
... )
-n,'2
In
(2r0-";2 exp [-(r
2 + 2r taii s i n Ot
d-n,'2dO
2) times
Jr-Ilal[2)/2]
n 2 ×
r"-t [ [cos.-l-i0id0._l...d01dr
(29)
i~l
where the last step is by transition to polar co-ordinates. ~t Notice that the domains o f integration for 01 in eqn (29) are from 0 to 2rt and from - rr,.'2 to 7r/2 for n = 2 and n _> 3 respectively. Equation (29) thus can be rewritten as
~:.:~.
cos"- 20J.(O) dO
(30)
where {0 xo =
if n = 2 rt"2
if n > 3
xt =
2n rr/2
if n = 2 if n_>3
e.=
1 2rt
ifn = 2 ifn>_3
n--2
:c. -- (2rr) '"'2
[1F
cos"- i -iOd 0
-~.2
i=2
J.(O) = I "r" i exp [ 3o
(r 2 +
2r [ia I1sin 0 + tl a !12)/'2] dr
(31)
An improved Monte Carlo method in structural reliability
291
(For n = 2, 3 the product in eqn (31) is understood to be unity.) It can be shown that (2n)-t ~.0t. =
if n = 2
2 ~n/2~- I n - 1
n - 2)!
-- 1
2'-"/2rt
-
i f n > 4 is e v e n
if n is o d d
3
I \-2--/
Finally, integration by parts shows that the following recurrence relation holds true for J.(O) J2(O) = ao + bo + CoJl(O) J.(O) = (n - 2)J. _ 2(0) + CoJ. - 1(0) + bo Rn - 2
n> 3
(32)
with ao = exp ( - IIa II2/2) bo = - e x p [ - ( R 2 +
2RIlall
sin0+
Ila 112)/2-1
Co = - II a 11 sin 0
In order to evaluate Jn(O) from eqn (32) only J,(O) needs to be known which can be written as Jl(O) = (2n) 1'2 exp [ -(tl a i12 cos 2 0)/2] x ((1)(R + !',all s i n 0 ) - ~(liall sin 0))
(33)
Taking in eqn (33) a high-precision approximation for the standard normal distribution function tl)s'° enables us to obtain J~ and hence J, by means of the recurrence relation eqn(32) without numerical integration. The calculation of the n-dimensional integral in eqn (28) is therefore reduced to a simple one-dimensional integration. It should be noted that a simpler approximation o f • which is usually sufficiently accurate for standard reliability calculations 13 proved inadequate for the purposes of this probability evaluation.
A P P E N D I X 2: C O N S T R U C T I O N OF T H E M A T R I X M An orthogonal matrix M satisfying eqn (12) is required. If a is non-zero (which may be assumed since otherwise the identity matrix is suitable) then the G r a m - S c h m i d t orthogonalization process t4 can be applied as follows. The basis to be orthogonalfzed is obtained from the standard basis of the ndimensional space by replacing its ith element by a with i t { 1..... n} chosen
202
.,! ttila ( ".~cnk i
such that a i ¢ O. D e n o t i n g the e l e m e n t s o f this basis by Pl . . . . . p._ 1,a the c o l u m n v e c t o r s o f M arc defined by the following r e c u r r e n c e relation m n --
mj :
a ii
! Ui
1 a,
i I Uj,
j = n --
where
uj=pi- ~ (p~m,)mk k-j~]
1, n - - 2 . . . . .
1