- Email: [email protected]
A new look at the response surface approach for reliability analysis
Structural Safety, 12 (1993) 205-220
205
Elsevier
A new look at the response surface approach for reliability analysis * Malur R. Rajashekhar and Bruce R. Ellingwood Department of Civil Engineering, The Johns Hopkins University, Baltimore MD 21218-2699, USA
Abstract. Closed-form mechanical models to predict the behaviour of complex structural systems often are unavailable. Although reliability analysis of such systems can be carried out by Monte Carlo simulations, the large number of structural analyses required results in prohibitively high computational costs. By using polynomial approximations of actual limit states in the reliability analysis, the number of analyses required can be minimized. Such approximations are referred to as Response Surfaces. This paper briefly describes the response surface methodology and critically evaluates existing approaches for choosing the experimental points at which the structural analyses must be performed. Methods are investigated to incorporate information on probability distributions of random variables in selecting the experimental points and to ensure that the response surface fits the actual limit state in the region of maximum likelihood. A criterion for reduction in the number of experiments after the first iteration is suggested. Two numerical examples show the application of the approach.
Key words: limit state design; probability; reliability; response surface; statistics; structural engineering
1. Introduction Structural reliability analysis requires a model of the structural system to describe its behavior under various loading conditions [1]. Structural behaviour is defined by a function g(X) of basic random variables X defining loads, strength of materials, and dimensions. Let g(X) = 0 define a limit state or failure function with g(X) <<.0 being the failure condition. The probability of failure or limit state probability is,
pf=e[g(X)
= ff (x) dx
(1)
where fx(x) is the joint probability distribution function for the vector X of basic variables. The region of integration is the domain .~ of x in which g(x) <~O. Equation (1) is difficult to evaluate for realistic structures and systems. Part of the difficulty arises in identifying the joint density function, fx(x), and in performing numerically the multi-dimensional integrations of f~(x) over the domain .~. Moreover, the limit state function, * Discussion is open until April 1994 (please submit your discussion paper to the Editor, Ross B. Corotis). 0167-4730/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
206
g(x) = 0, often is determined from finite element analysis or other approximate computational methods, making it possible to define the domain 2 only pointwise rather than in closed form (e.g., [2]). First order reliability analysis addresses the first difficulty by mapping the random variables into the space of uncorrelated standard normal variables U, linearising the limit state surface g(u) -- 0 at the point with minimum distance,/3, from the origin (denoted the minimum norm or design point), and estimating the limit state probability as, pf = (I)(--/3)
(2)
in which @(.) = standard normal probability distribution function [3-7]. However, because of the nature of first-order reliability analysis, it is most convenient if g(X) is available in closed-form. There are many practical reliability analysis situations where the limit state cannot be represented mathematically in closed form or the numerical integration cannot be performed conveniently. Monte Carlo simulation, often involving some variance reduction technique to reduce the required number of samples [7,8], can be used in such situations. If the limit state probability with respect to sidesway of a multi-storeyed multi-bay framed structure were to be calculated by simulation, an analysis of the structure would have to be performed for each random sample of loads and strengths, leading to hundreds or perhaps thousands of analyses. As an alternative, if the sidesway were to be approximated as a polynomial function of random variables, the coefficients of each term in the polynomial could be obtained by performing the structural analysis for a lesser but sufficient number of times. This polynomial may then be used as a basis for simulations; provided that the polynomial fits the actual limit state reasonably well, one could obtain a fairly accurate estimate of the probability of failure. Such a polynomial is referred to as a Response Surface.
2. Response Surfaces 2. 1. Basic ideas and assumptions In a structural system, a large number of variables influences the response of the system, e.g., the maximum crack width, deflection or stress at a particular point of interest. Suppose response variable Y depends on the input variables X1, X 2. . . . , X n. Experiments are conducted with design variables X 1, X 2 , . . . , X n a sufficient number of times to define the response surface to the level of accuracy desired. Each experiment can be represented by a point with coordinates xlj, x2j,...,x,,j in an n-dimensional space. At each point, a value of yj is observed. Although the actual response Y is a function of input variables, i.e., Y = g(X1,..., Xn), this function is generally unavailable in closed form. The basic response surface procedure is to approximate g(X) by an nth order polynomial g(X) with undetermined coefficients. Structural analysis is performed at various points x~, in order to determine the unknown coefficients in the polynomial ~(x) such that the error of approximation is minimum in the region of interest.
2. 2. Selection of the order of the polynomial and quadratic approximations The selection of the order of the approximating polynomial and points x i for experimentation requires careful consideration. The degree of if(X) should be less than or equal to the
207 degree of g(X) to get a well-conditioned system of linear equations for the unknown coefficients. If ~(X) is of much higher degree than g(X), one obtains an ill-conditioned system of equations. Moreover, higher order polynomials can exhibit erratic behaviour in the sub-domains not covered by the experiments [9]. Of course, g(X) itself is not known. U p to a certain degree, a higher order polynomial improves the accuracy of the approximation at the expense of additional computation. The rate of increase in accuracy reduces with the degree of polynomial but the computational costs increase exponentially, as a higher order polynomial involves greater number of unknown coefficients and requires correspondingly more structural analyses. For reliability estimates, one needs to have a good approximation to g(X) around the design or minimum norm point, or the region of the failure domain .q~ where the joint probability density of X is relatively large and thus contributes most to the overall failure probability. Since we neither know the actual limit state function nor the actual design point, the accuracy of the reliability estimate depends on the accuracy of the polynomial approximation in the region of the design point. With an eye toward this accuracy as well as the high cost of repeated finite element analyses, a second degree polynomial is used in the present study. A second order response surface for n input variables is described by a quadratic model,
=A +xTB +xTcx
(3)
where A, B T = [B1, B2,..., nn] , and, C ~--
Cll
"'"
syn'l
Cln] Cnn]
define the undetermined coefficients. Experiments are conducted as per the adopted design and the resulting system of equations may be put in the form,
g(x) =De +e
(4)
where d is a vector of constants A, B i, Cij, the matrix D contains constant, linear, quadratic and cross-combinations of x i and e is the error vector, the components e i of which consist of a systematic or lack of fit error resulting from approximating g by if, and a pure experimental error, assumed to be a zero mean random vector. The solution,
E(d) = (DTD)-IDTg(x)
(5)
provides the expected values of the unknown coefficients. Other polynomial interpolation schemes using Lagrangian and Hermite polynomials are possible, although no specific examples of their application for reliability analysis could be located. At a higher level of sophistication, Ditlevsen [10] has presented a random field model for stochastic interpolation between point by point measured values of a spatially distributed material property.
2.3. Experimental designs In the absence of prior information, the experiments initially can be centered around the mean values of the variables. Points are located along the axes and elsewhere. The axial points
208
define the effect of each variable and are needed to determine coefficients for the terms involving only that variable i.e., x i, x 2 etc. The points in the plane x i x j define the interaction effects of variables x~ and xj. and are needed to determine coefficients for the cross terms, i.e., x~x i etc. When including the cross terms, a particular design of experiments depending on the nature of variables is discussed subsequently in the context of a specific response surface. The commonly used orthogonal experimental designs are 2 n and 3 n factorial and fractional factorial designs [11,12]. The 2 n factorial points for a simple two variable problem are indicated in Fig. 1 by , . Fractional factorials are useful when the number of variables is large and hence the number of experiments that can be conducted is less than the number of combinations in the full factorial set. Augmented 2 n factorials are obtained by adding n 2 points at the center of the design (total points = 2 ~ + n2). The center point is indicated by • in Fig. 1. Another efficient class of designs for fitting a second order response surface is the so-called central c o m p o s i t e design. This consists of a 2 ~ factorial design, augmented by n 2 > 1 central points and 2 ~ points placed at coordinates - a and + a along the axes. These axial points are indicated by ~9 along each axis in Fig. 1. Hence the design consists of (2" + 2n + n 2) points [13]. Other experimental designs include Simplex designs for first order models and Equiradial designs for second order surfaces. Khuri and Cornell [14] have discussed various experimental designs in detail. These factorial designs, though efficient, lead to unacceptably high computational efforts with the increase in number of variables for complex systems and may become more time consuming than simulation. An iterative response surface approach for reliability analysis was presented by Bucher and Bourgund [15]. The experimental design in each iteration consists of as many locations as the total number of undetermined coefficients in the polynomial i(X) =a
bixi+
+ i=1
Y'~cix2i i=1
(6)
in which xi, i = 1,..., n are basic variables and the parameters a, bi, Ci are to be determined. In order to obtain these constants, (2n + 1) experiments have to be conducted. Such an experimental design is referred to as fully saturated, and the surface fits exactly at the experimental points. The suggested experimental points for obtaining if(x) lie along axes x i (see Fig. 1). The points chosen are mean values tz i of the random variable X i , and xi = [1.1i, Jr h p i
X2
0
(
0
XI
Fig. 1. Location of experiments for a two variable problem.
209 in which h i is an arbitrary factor and o-i, is the standard deviation of X i. The constants are determined using (2n + 1) values of g ( x ) at these points. In the next step ~(x) is used along with information on/J,i and o-i to determine the minimum norm point x o for the surface ~ ( x ) = 0, based on the assumption that x are uncorrelated Gaussian variables. Once x o is found, g ( x o ) is evaluated and a new center point x M for interpolation is chosen on a straight line from mean vector ~ to x 9 so that g ( x ) = 0 from linear interpolation, i.e.,
XM= ~ +(Xo--~)[g(~l-----g(XO) ]
(7)
This strategy is to position the new center point reasonably close to g(x) = 0. A new surface obtained using center point x m may reasonably define the actual response in the failure domain. The total number of evaluations of g ( x ) is (4n + 3). Bucher and Bourgund [15] assert that this method of updating the polynomial results in a sufficiently accurate response surface. However, the writers believe that the accuracy depends on the characteristics of the limit state being explored and hence one cycle of updating may not always be sufficient. Key issues in this method are the selection of h i (which is arbitrary) and a suitable criterion to ensure that the fitted surface if(x) does indeed represent the actual surface in the failure domain reasonably well. The following section describes some ideas and approaches that address these issues.
3. Improvements to response surface methods
3.1. Criteria for fitting O(x) to g(x) The accuracy of the response surface obtained after two cycles of iteration and 4n + 3 experiments varies with each problem analysed, as shown subsequently, and thus the accuracy of the failure probability computed on the basis of the approximate response surface is unpredictable. One obvious question is how far should the process of interpolation and experimentation continue? Figure 2 illustrates the process of successive approximations, design points and center points for a two-variable limit state. The first approximation is obtained by conducting experiments at the mean values/z i denoted by XM1 in Fig. 2 and at (ix i + h i o i ) . The minimum norm point x m is determined for the fitted surface ~l(x) = 0. By linear interpolation, the new center point XM2 for experimentation is found from Equation 6 and a new surface ~ 2 ( x ) = 0 is then obtained. For this surface, again the minimum norm point x o 2 is obtained but the distance is measured from XM2 instead of from the mean /~(=Xul). This results in a low value for this distance, which is henceforth referred to as 8. ~-Fhe process is repeated until 6 becomes very small or zero. Then the factor h i is reduced and the process is repeated until 6 becomes zero again. This can be achieved in about three to four iterations, resulting in a very good approximation for the actual surface and a value of x 9 that is close to the actual minimum norm point. However, if the process is repeated indefinitely with progressively smaller values of h i , then at some stage one obtains an ill-conditioned system of equations. Experience has shown that starting the process with h i = 2.0 and then repeated once more with h i = 1.0 gives a satisfactory
210 LEGEND : X t, X= = A two v ~ i a b l e
..... ....
X2 'x' ~'-~g.IL J
space
for
g(x) m Actual Limit State g t ( x ) =, Raeponee S u r f , a c e gel x ] = Ra=ponse SuP f a c e g~l(x] = Reeponsa S u r f , a c e -
illustration
Function Fir'st Approximation Second A p p r o x t mat t on Third Approximation
xul - Center Point t,or the t n approximation x m - D e s i g n P o i n t t,or t h e i th a p p r o x i m a t i o n +g - Ihnasur'ed C o n v e r g e n c e P a r a m e t e r t ' o r t h e i ~h a p p r o x i m a t i o n
9(x)
-
xui
X1
Fig. 2. Convergence to the minimum norm point: Approach A-0.
response surface around the minimum norm point. Once the surface is obtained, advanced simulation techniques can be used to estimate the failure probability [7,16]. This procedure is referred to as Approach A-0 in the examples to follow in Section 4. When the variables are dependent a n d / o r not Gaussian, they can be transformed into an independent standard normal space and the procedure as outlined above can be carried out to ensure the criterion based on ~. Experience has shown that this improves the conditioning of the system of equations as well.
3.2. Experiments at distribution extremes The tail regions of the probability distributions of random load and strength are most important in the estimation of failure probabilities. Hence, to obtain an accurate fit to the actual surface in the vicinity of the design point, one might think that experimental points should be chosen in the lower tails of the resistance variables and upper tails of the load variables (Figs. 3 and 4) rather than over their entire range. Three different approaches incorporating experiments at distribution extremes are considered in this section. t'.( x 2 )
f x ( Xl )
12345
Xl
12545
Fig. 3. Experiments at distribution extremes: Approach A-1.
~ x2
211
f,( xl )
f',,( ×2]
k,tktt~t --- x I
X2
Fig. 4. Experiments at distribution extremes: Approach A-2.
In the first approach (Fig. 3), the width of the sampling region for each variable is kept sufficiently large and in the first iteration, as many points are chosen in this region as the number of experiments. The experiments are conducted at points whose coordinates are rank ordered locations in the tail region of each variable. Thus the first approximation is obtained and then the minimum norm point for the first approximation is determined. For the second iteration, the width of the sampling region is reduced around this point and the distance 8 is measured as before (Section 3.1). At low values of 8, a fairly accurate response surface results, as shown subsequently in Example 1 in Section 4.1. This approach is designated henceforth as Approach A-1. In the second approach (Fig. 4), the effect of shape of the probability distributions of variables on the selection of experimental points is considered. Using a two-variable problem as an illustration, the center point for the first iteration is chosen at x = (x 1, x 2) where x 1 and x 2 are selected at k ~ and k ~ fractiles (e.g., typically the 95th percentile of load or 5th percentile of resistance) of random variables X 1 and X 2 (Fig. 4). Two other points on either side of both k n and k12 are also chosen. The experiments are conducted at ( k l l , k12) , (kol , k12) , (k21 , k12) , (k n, k02) and ( k l l , k22). These points are at the k ~ fractiles of the respective distribution functions. Once the minimum norm point for the first trial surface is obtained, sampling is concentrated around that point and the process is repeated as before. The same logic is applicable for an n-variable problem. This approach is designated henceforth as Approach A-2. As a variant of this approach, points on either side of the k ~ fractiles of X i can be chosen as a constant times its standard deviation, as done in Approach A-0, i.e., the experiments are conducted at ( k l l , k12) , ( k l l -t- hlOrl, k12) and ( k l l , k12 -k- h20"2). However, h i would typically be less than what it is when experimenting around the mean but large enough to have a well-conditioned system of equations. This approach is henceforth designated as Approach A-3. The objective in procedures A-1 through A-3 is to include distribution information in the selection of experimental points. Example 4.1 in Section 4 shows that fairly accurate response surfaces can be generated from these procedures.
3.3. Addition of terms to quadratic O(x) The accuracy of the response surface may be improved by inclusion of the cross terms in the polynomial [15,17] as follows, n ff(x)=a+
n
Ebixi.-]- E i=1
~
cijxixj
i=1 j>~i=l
(8)
212 TABLE 1 Quadratic combinations for experimental locations Resistance(Xi)* Load(Xj) Resistance(Xi)* Resistance(Xj) Load( X i ) * Resistance(Xj) Load(X/) * Load(X~)
- ( x i ) . +(x~) - (xi)* - (x j) + ( x ) . - (xj) +(x).(xp
The total n u m b e r of experiments to be conducted for each approximation would then increase to (n + 1 ) x (n + 2)/2. To gain some insight on the appropriate location of experiments, consider a two variable limit state function approximation ~(x 1, x2). There is only one cross t e r m (clzxlX 2) in this polynomial which has six constants in all. The first five locations are taken at the m e a n and four axial points. The last location for the cross term has to be selected among the four choices indicated by * in the Fig. 1. The best location depends on w h e t h e r X~ and Xj are "resistance" or "load" variables, as summarized in Table 1. The same logic should be applied for a cross term involving a greater n u m b e r of variables, by selecting points whose coordinates lie on the negative side of the "resistance" axes and positive side of the " l o a d " axes, with the origin of the axes taken at the center point. Such a combination results in the best possible fit in the failure domain among all possible combinations. This is illustrated in Example 4.2 in Section 4. The inclusion of cross terms in the approximation increases the number of experiments from (2n + 1) to (n + 1) x (n + 2)/2, an increase of n x (n - 1 ) / 2 which is significant if the reliability analysis involves a large n u m b e r of variables. However, totally neglecting the cross terms may be inappropriate as there are situations where cross terms are more important to the actual response than the squared terms of the individual variables. As a simple example to illustrate this point, the limit state of a truss m e m b e r subjected to axial force is defined by A x f y - P = 0, A, fy, and P being the area of cross section, yield stress and applied force respectively; it is obvious that a better response surface results by taking cross terms than squared terms. For efficiency, only those terms which are important in terms of accuracy should be included. As the iterations proceed, terms that are relatively unimportant can be neglected, reducing the number of experiments at each cycle. If gc is the response computed at the center point and gi is the response computed at the ith experimental point, then the square of the difference normalised with respect to gc gives the relative importance of the ith location. If,
( gi-gc l2< c. gc
(9)
]
then the corresponding term may be neglected in subsequent iterations. The neglected term may be a linear or a square or a cross term; this is more realistic and accurate than just neglecting the cross terms and should lead to better subsequent approximations. This is illustrated in Example 4.2 in the following section.
4. Numerical examples Two numerical examples are presented in this section. The cantilever beam in the first example facilitates visualization of the accuracy of successive approximations by response
213 surfaces of the actual limit state. In the second example, the non-linear dynamic response of a single degree of freedom system with random parameters is analysed and compared with a previously published response surface analysis.
4. 1. Cantilever beam A cantilever beam with a rectangular cross section subjected to uniformly distributed loading is considered. The limit state of serviceability with respect to maximum deflection at the free end being greater than 1/325, is defined by (W X b ) l 4 g =
8EI
l
+ 325
(10)
where w, b, l, E and I are respectively load per unit area, width, span, modulus of elasticity and moment of inertia of the cross section. The random variables are the load per unit area and the depth of the cross section. Assuming E and l as deterministic and equal to 2.6 × 104 MPa and 6 m respectively, the limit state equation reduces to g ( X ) = 18.46154 - 7.476923 X 101° X---~= 0 X2
(11)
where X 1 is the load in MPa and X 2 is the depth in mm. The load is assumed to be normally distributed, with a mean ~.L1 1000 N / m E and a coefficient of variation V1 -- 0.2. The beam depth is also normally distributed, with ~2 = 250 mm and V2 = 0.15. The probability of failure for this limit state is calculated according to the various procedures described in earlier sections. The "exact" solution using importance sampling with 1000 simulations is 9.6071 × 10 -3. Advanced first order reliability analysis gives a reliability index of 2.33 corresponding to a probability of failure of 9.9031 × 10 -3. The "exact" design point is located in Figures 5-9 that follow. Consistent with eqn. (11) the response surface is assumed to have the form of Equation 5, with no cross terms, and is obtained for this limit state by carrying out analyses along axial points(see Fig. 1) Curves 1 and 2 in Fig. 5 represent Bucher's algorithm with h i = 2.0. Simulation using curve 2 gives a probability of failure of 13.7538 × 10 -3. Following Approach =
8
I
I
I
7N
x
I
I
I
I
I
I
I
I
....................
I
.
I
I
~
5 - ' ~ 4
3 2
f ["
~ 0.0
Exit
C~esents
Response S u r f a c e
Limit
State
o f I th I t e r a t i o n Xo represents "exact" des1 gn p o t n t
1
0
C ~ e n t s
~ 0.2
i
i 0.4
i
i 0.6
I
l 0.8
I
I 1.0
l
i 1.2
X1/~1 F i g . 5. C o n v e r g e n c e - - A p p r o a c h
A-0 - h i = 2.0.
l
I 1.4
214 I
0.8
I
I
I
I
I
I
I
I
I
I
I
I
I
0.70,6fw :1
X
0,50,40.50.2-
LEGEND
f •
C ~ e n t s Exact L i ~ t State C ~ e s e n t s Response Surface of I ~ I t e r a t i o n represents "exact" design potnt
0.1 0.0
I 0.0
t
I
0.2
I 0.4
I
I 0.6
I
I 0.8
I
I
I
I
1.0
I
1.2
I 1.4
xl/~x Fig. 6. Convergence--ApproachA-0
-
hi
=
1.0.
A-0 to improve the response surface, successive surfaces are obtained with h i = 2.0 (curves 3 and 4 in Fig. 5). Fig. 6 shows the surfaces obtained in the second stage of iteration with h i = 1.0. Table 2 gives the distance 6, between successive center and design points. The probability of failure obtained from simulations on the final response surface (curve 6) is 9.5410 x 10-3which is very close to the exact solution. Thus Approach A-0 yields a very close result; the n u m b e r of iterations required to fit the response surface to g(x) = 0 is limit state d e p e n d e n t and cannot be generalised as two or three cycles. Progressively smaller values of h i lead to an ill-conditioned system of equations. Table 3 illustrates the effect on the solution when the response surfaces are generated by locating experiments in the tail regions of the random variables using Approach A-1. The process is initiated by selecting a relatively large range of 0.75 to 0.95 fractiles for the load variable and 0.05 to 0.25 fractiles for the resistance variable. Once the design point is obtained, the range is successively reduced until the distance 6 becomes very small. Fig. 7 illustrates the successive approximations; sampling using curve 3 of Fig. 7 yields a probability of failure of 9.6398 x 10 -3 which is very close to the exact solution of 9.6071 x 10 -3. It can be seen from Fig. 6 and Fig. 7 that this technique yields a response surface with a better fit to the exact limit state function over a wider range compared to that obtained by Approach A-0 (Table 2). The drawback of this approach lies in the selection of experimental points for the first iteration. With an increase in the n u m b e r of variables, it becomes increasingly difficult to choose the range for each variable.
TABLE 2 Successive Response Surface fits to
g(x)= 0 using Approach A-0 (cf Figs. 5 and 6)
Iteration
hi
6
XloiN/m2
X2Dimm
1 2 3 4 5 6
2.0 2.0 2.0 2.0 1.0 1.0
2.6436 0.2868 0.0169 0.0002 0.0002 0.0000
1069.29 1064.15 1073.15 1072.40 1072.38 1072.38
151.72 168.36 162.29 163.52 163.15 163.15
Experiments at XMg and
(XMg +_hio"i)
215 TABLE 3 Experiments at distribution extremes: Approach A-1 (cf Fig. 7) Iteration
Variable
kI
k2
k3
k4
1
x~ x2 x1 x2 xI x2
0.75 0.05 0.78 0.0050 0.7365 0.0120
0.80 0.10 0.83 0.0075 0.7865 0.0145
0.85 0.15 0.88 0.0100 0.8365 0.0170
0.90 0.20 0.93 0.0125 0.8865 0.0195
2 3
ks 0.95 0.25 0.98 0.0150 0.9365 0.0220
kdesign
t~
0.8810 0.0107 0.8365 0.0170 0.8413 0.0158
1.51 0.3327 0.0456
Fig. 8 shows plots of successive approximations obtained from Approach A-2 of Section 3.2, incorporating the shape of the probability distributions in selecting experimental points in the tail region. C o m p a r e d to Figs. 5-7, the response surfaces fit well over a narrower range. Table 4 gives the details of two iterations; the probability of failure obtained is 11.1508 × 10 -3 which is 16% higher than the exact value. Table 5 and Fig. 9 show details and plots of successive
.8
°.,
0
.
o.st 20.4
7
.....
~
......
....
t~ 0"-~JI- / - - I /"
C ~ ~ n t s Exact Limit State Curve / r e p r e s e n t s Response S u r f a c e of i ~ Ite;etion " X0 r e p r e s e n t s " e x a c t " d e s | g n , p o t n t
0 , 1 --'~
0
. 0.2
0.0
0 0.4
~ 0.6
0.8
1. 0
1.2
1.4
xt/~l Fig. 7. Convergence--Approach A-1.
I
0.8
I
I
I
I
I
I
I
I
I
I
I
I
I
0.70.6-
- . ....... ~ .....
r
NO.5:l ~O.4x 0.5-
=
:
=
:
=
~
~
Curve E r e p r e s e n t s Exact L i m i t S t a t e Curve I r e p r e s e n t s R e s p o n s e ~ J r f a e e of I m I t e r a t i o n represents "exect" design p o | n t
0.20.1 0.0 0.0
:
I
I
0.2
I
I
0.4
I
I
0.6
I
I
0.8 Xl/~l
I
I
I
1.0
Fig. 8. Convergence--ApproachA-2.
I
1.2
I
I
1.4
216 TABLE 4 Experiments at distribution extremes: Approach A-2 (cf Fig. 8) Iteration
Fractiles selected
Experiments at kth fractile
Design point at kth fractile
6
1
x a ~ (0.7,0.8,0.9) x 2 ---,(0.01,0.1,0.2)
(0.87,0.017)
1.0778
2
x 1 ~ (0.80,0.87,0.94) x 2 ~ (0.01,0.017,0.024)
(0.8,0.1) (0.8,0.01),(0.8,0.2) (0.7,0.1),(0.9,0.1) (0.87,0.017) (0.87,0.01),(0.87,0.024) (0.80,0.017),(0.94,0.017)
(0.87,0.017)
0.0003
TABLE 5 Experiments at distribution extremes: Approach A-3 (cf Fig. 9) Iteration
Mean point for experiment as a fraction of actual mean, x M / ~
Design point as a fraction of actual mean, x D/tz
6
1
(1.256,0.808) corresponding to (0.9,0.1) fractile (1.0673,0.6695) (1.0774,0.6559) (1.0788,0.6542) (1.0790,0.6540)
(1.0673,0.6695)
2.2291
(1.0774,0.6559) (1.0788,0.6542) (1.0790,0.6540) (1.0790,0.6540)
0.1438 0.0181 0.0026 0.0003
2 3 4 5
Experiments at (kl, k2), (k 1 + or1, k 2) and (kl, k 2 + ~2)
a p p r o x i m a t i o n s o b t a i n e d f r o m A p p r o a c h A-3 l o c a t i n g e x p e r i m e n t s at k i a n d ki _+ tri in t h e tail region. T h i s gave a p r o b a b i l i t y o f f a i l u r e o f 9.5410 × 10 -3 w h i c h c o i n c i d e d with t h a t o b t a i n e d by A p p r o a c h A - 0 ( T a b l e 2). T h e results o b t a i n e d f r o m v a r i o u s p r o c e d u r e s a r e s u m m a r i s e d in T a b l e 6. It is e v i d e n t f r o m the a b o v e discussion t h a t i n c o r p o r a t i n g i n f o r m a t i o n o n t h e p r o b a b i l i t y d i s t r i b u t i o n s o f t h e basic
:
0.8 0.7
I
:
I ..-~
0.6
I
I
I
I
I
I
I
I
......................................
I
I
xD
~0.5 0.4 x
0,5
//
0.2- / O,l I 0.0 0.0
C
~
ents Response Surface
e
0 f I~h I t e r a t i O n X0 r e p r e s e n t s
0.2
0.4
0.6
.
.
.
.
.
"exa(: " t d e s t g n p o i n t
0.8
1.0
Xl/l~t
Fig. 9. Convergence--Approach A-3.
1.2
1.4
217 TABLE 6
Comparison of various methods Methods
Probability of failure
Simulation using actual limit state model Advanced First Order Second Moment Method Bucher's Method Approach A-0 Approach A-1 Approach A-2 Approach A-3
0.0096071 0.0099031 0.0137538 0.0095410 0.0096398 0.0111508 0.0095410
variables in the selection of experimental points at which the finite element analysis for complex structures is carried out yielded fairly accurate response surfaces but did not save computational effort required to obtain a good approximation (measured, for comparable level of accuracy, by 6) to the probability of failure. Hence, it can be concluded that sampling in the tail regions of the random variables instead of the entire region does not lead to significant improvements in the response surfaces or reliability estimates. 4.2. Dynamic response of a nonlinear oscillator A non-linear undamped single degree of freedom system analysed by Bucher et al. [15,17] is illustrated in Fig. 10. The limit state is defined by
g(x) =
3r-
I Zmax I
(12)
in which Zmax is the maximum displacement response of the system and r is the displacement at which one of the springs yields. The statistical parameters of the basic variables are listed in Table 7. All variables are assumed to be independent and lognormally distributed.
I
.
C l
~
z(t)
_ ~ _ ~ . F( t ) r
C2 &l
Restoring
r
force
~
z
il
F(t)
tl
~t
Fig. 10. Nonlinear SDOF oscillator with nonlinear restoring force and pulse load for Example 2.
218 TABLE 7 Statistical data for Example 2 Number
Variable
Mean value
Standard deviation
1 2 3 4 5 6
m c1 c2 r F1 t1
1.00 1.00 0.10 0.50 1.00 1.00
0.05 0.10 0.01 0.05 0.20 0.20
Importance sampling with 1000 simulations of nonlinear dynamic response yielded a probability of failure of 3.99864 x 10 -2. Bucher et al. [15,17] obtained a probability of failure of 3.81 x 10 -2 by adaptive sampling using 1500 simulations, 3.01 x 10 -2 using a response surface without cross terms, and 3.67 x 10 -2 using a response surface with cross terms. Approach A-0, using a response surface without cross terms and with just one more cycle, gives a failure probability of 3.93365 x 10 -2. Table 8 gives details of the iterations. The effect of cross terms are considered for this example only for the first iteration as it was found that a good approximation could be obtained in the second iteration without them. W h e t h e r the cross terms may be neglected altogether depends on the complexity of the problem and the nature of the unknown limit state. Without cross terms, the first response surface approximation to the limit state gives a probability of failure of 2.61898 x 10 -2. The experimental points for cross terms are then selected in the first quadrant (see Fig. 1) and the resulting failure probability associated with the revised response surface is 3.22473 x 10 -2. In order to study the effect of location of experiments, the procedure is repeated with experimental points in the second, third and fourth quadrants. With the minimum norm distance for the approximating surfaces taken as an index of the accuracy, the result is almost same regardless of the quadrant in which the points are located, as indicated in Table 9. The suggested procedure of selecting points on the negative side of the axes for "resistance" variables and positive side for " l o a d " variables with the mean as the origin gives the best result in terms of both minimum fl and closest estimate of the failure probability. The probability of failure is 3.42005 x 10 -2 in the first iteration, which is quite comparable to the exact value of 3.99864 x 10 -2 The importance of the individual terms in the polynomial response surface is evaluated using eqn. (9). The proposed index is calculated for each term and terms for which the squared error is less than 0.01 (E) are neglected. In all, five terms (out of 21) are neglected. Treating this as a first iteration, a new surface is fitted which gives an almost identical failure probability of
TABLE 8 Results of Approach A-0 for Example 2 Iteration 1 2 3 Non linear dynamic response
6 2.08001 0.14650 0.09543
Probability of failure 2.61898 x 10 - 2 4.09679 x 10 -2 3.93365 x 10 - 2 3.99864 x 10 - 2
219 TABLE 9 First iteration accuracy Method
/3
Probability of failure by simulation
Without cross terms Cross terms (x i, x j) in I quadrant ( + , + ) II quadrant ( - , + ) III quadrant ( - , - ) IV quadrant ( + , - ) Resistance ( - ) and Load ( + ) With reduced crossterms Simulation by non-linear analysis
2.08001
2.61898 x 10 2
1.90273 1.89055 1.90991 1.91874 1.88704
3.22473 x 10- 2
1.88882
3.40402 x 10- 2
3.42005 x 10 - 2
3.99864 x 10-2
3.40402 x 10 -2. Thus the terms neglected do not have a significant effect on the response surface or probability of failure. The results are summarised in Table 9. The accuracy of the response surface obtained depends on the characteristics of the actual limit state and is affected by the selection of points at which the experiments are performed. The proposed method (Approach A-0) of continuing the response surface fitting process until the distance 6 becomes very small or zero is a reasonable way to tackle reliability analysis of complex structural systems. On the basis of the limited studies made, the suggested location of experimental points for finite element analysis seems to be logical. The index in eqn. (9) is sufficient to test the significance of each term after the first iteration so that unnecessary experiments can be reduced in the subsequent iterations.
5. Summary and conclusions An adaptive iterative procedure was presented to develop a response surface that is locally a good approximation to the actual limit state surface in the region of maximum joint probability density and can be used for structural reliability analysis. Methods were presented to enable inclusion of information on probability distributions of random variables in selecting experimental points used to generate response surfaces. However, it was found that such methods have no particular advantages with regard to computational effort and accuracy as compared to more conventional methods for selecting experimental points. The location of experiments for cross terms and a systematic procedure to reduce the number of experiments depending upon the importance of each term in the previous approximation also were suggested. It is felt that the response surface method can be used efficiently for reliability analysis of certain structural systems where behavioral models to describe their various limit states cannot be developed in closed form.
6. Acknowledgements The first author expresses special thanks to Professor A.S.R. Sai at his parent institution, the Indian Institute of Technology, Kanpur, India, for initiating him into research in the field of
220
reliability. He also expresses grateful acknowledgement to the Council for International Exchange of Scholars, Washington D.C., for sponsoring his research at the Johns Hopkins University.
References 1 R.E. Melchers, Structural Reliability : Analysis and Prediction, Ellis Horwood Limited, Chichester, 1987. 2 A. Der Kiureghian, and J.B. Ke, The Stochastic Finite Element Method in Structural Reliability, Stochastic Structural Mechanics, Lecture Notes in Engineering 31, Springer, Berlin, 1987. 3 A.M. Hasofer and N.C. Lind, Exact and invariant second moment code format, J. Engng. Mech. Div., 100 (EM1) (1974) 111-121. 4 R. Rackwitz and B. Fiessler, Structural reliability under combined random load sequences, Computers and Structures, 9 (1978) 490-494. 5 O. Ditlevsen, Principle of normal tail approximation, J. Engng. Mech. Div., 107 (EM6) (1981) 1191-1208. 6 M. Shinozuka, Basic analysis of structural safety, J. Struct. Div., 109 (3) (1983) 721-740. 7 R.E. Melchers, Importance sampling in structural systems, Structural Safety, 6 (1989) 3-10. 8 G.I. Schueiler and R. Stix, A critical appraisal of methods to determine failure probabilities, Structural Safety, 4(4) (1987) 293-309. 9 S. Engelund and R. Rackwitz, Experiences with experimental design schemes for failure surface estimation and reliability, Proceedings of the Sixth Speciality Conference on Probabilistic Mechanics and Structural and Geotechnical Reliability, ASCE, 1992, pp. 252-255. 10 O. Ditlevsen, Random field interpolation between point by point measured properties, Proceedings, 1st International Conference on Computational Stochastic Mechanics, Corfu, Greece, 17-19 September 1991. 11 R.H. Myers, Response Surface Methodology, Allyn and Bacon, Boston, 1971. 12 R.G. Petersen, Design and Analysis of Experiments, Marcel Dekker, New York, 1985. 13 L. Faravelli, Response surface approach for reliability analysis, J. of Engng. Mech., 115 (12) (December 1989) 2763-2781. 14 A.I. Khuri and J.A. Cornell, Response Surfaces, Designs and Analyses, Statistics Textbooks and Monographs series, Volume 81, Marcel Dekker, New York, 1987. 15 C.G. Bucher and U. Bourgund, A fast and efficient response surface approach for structural reliability problems, Structural Safety, 7 (1990) 57-66. 16 C.G. Bucher, Adaptive sampling--an iterative fast Monte Carlo procedure, Structural Safety, 5(3) (1988) 119-126. 17 C.G. Bucher, Y.M. Chen and G.I. Schueller, Time variant reliability analysis utilising response surface approach, in: P. Thoft-Christensen, ed., Reliability and Optimization of Structural Systems '88, (Lecture Notes in Engineering 48), Springer, Berlin, 1989, pp. 1-14.
205
Elsevier
A new look at the response surface approach for reliability analysis * Malur R. Rajashekhar and Bruce R. Ellingwood Department of Civil Engineering, The Johns Hopkins University, Baltimore MD 21218-2699, USA
Abstract. Closed-form mechanical models to predict the behaviour of complex structural systems often are unavailable. Although reliability analysis of such systems can be carried out by Monte Carlo simulations, the large number of structural analyses required results in prohibitively high computational costs. By using polynomial approximations of actual limit states in the reliability analysis, the number of analyses required can be minimized. Such approximations are referred to as Response Surfaces. This paper briefly describes the response surface methodology and critically evaluates existing approaches for choosing the experimental points at which the structural analyses must be performed. Methods are investigated to incorporate information on probability distributions of random variables in selecting the experimental points and to ensure that the response surface fits the actual limit state in the region of maximum likelihood. A criterion for reduction in the number of experiments after the first iteration is suggested. Two numerical examples show the application of the approach.
Key words: limit state design; probability; reliability; response surface; statistics; structural engineering
1. Introduction Structural reliability analysis requires a model of the structural system to describe its behavior under various loading conditions [1]. Structural behaviour is defined by a function g(X) of basic random variables X defining loads, strength of materials, and dimensions. Let g(X) = 0 define a limit state or failure function with g(X) <<.0 being the failure condition. The probability of failure or limit state probability is,
pf=e[g(X)
= ff (x) dx
(1)
where fx(x) is the joint probability distribution function for the vector X of basic variables. The region of integration is the domain .~ of x in which g(x) <~O. Equation (1) is difficult to evaluate for realistic structures and systems. Part of the difficulty arises in identifying the joint density function, fx(x), and in performing numerically the multi-dimensional integrations of f~(x) over the domain .~. Moreover, the limit state function, * Discussion is open until April 1994 (please submit your discussion paper to the Editor, Ross B. Corotis). 0167-4730/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
206
g(x) = 0, often is determined from finite element analysis or other approximate computational methods, making it possible to define the domain 2 only pointwise rather than in closed form (e.g., [2]). First order reliability analysis addresses the first difficulty by mapping the random variables into the space of uncorrelated standard normal variables U, linearising the limit state surface g(u) -- 0 at the point with minimum distance,/3, from the origin (denoted the minimum norm or design point), and estimating the limit state probability as, pf = (I)(--/3)
(2)
in which @(.) = standard normal probability distribution function [3-7]. However, because of the nature of first-order reliability analysis, it is most convenient if g(X) is available in closed-form. There are many practical reliability analysis situations where the limit state cannot be represented mathematically in closed form or the numerical integration cannot be performed conveniently. Monte Carlo simulation, often involving some variance reduction technique to reduce the required number of samples [7,8], can be used in such situations. If the limit state probability with respect to sidesway of a multi-storeyed multi-bay framed structure were to be calculated by simulation, an analysis of the structure would have to be performed for each random sample of loads and strengths, leading to hundreds or perhaps thousands of analyses. As an alternative, if the sidesway were to be approximated as a polynomial function of random variables, the coefficients of each term in the polynomial could be obtained by performing the structural analysis for a lesser but sufficient number of times. This polynomial may then be used as a basis for simulations; provided that the polynomial fits the actual limit state reasonably well, one could obtain a fairly accurate estimate of the probability of failure. Such a polynomial is referred to as a Response Surface.
2. Response Surfaces 2. 1. Basic ideas and assumptions In a structural system, a large number of variables influences the response of the system, e.g., the maximum crack width, deflection or stress at a particular point of interest. Suppose response variable Y depends on the input variables X1, X 2. . . . , X n. Experiments are conducted with design variables X 1, X 2 , . . . , X n a sufficient number of times to define the response surface to the level of accuracy desired. Each experiment can be represented by a point with coordinates xlj, x2j,...,x,,j in an n-dimensional space. At each point, a value of yj is observed. Although the actual response Y is a function of input variables, i.e., Y = g(X1,..., Xn), this function is generally unavailable in closed form. The basic response surface procedure is to approximate g(X) by an nth order polynomial g(X) with undetermined coefficients. Structural analysis is performed at various points x~, in order to determine the unknown coefficients in the polynomial ~(x) such that the error of approximation is minimum in the region of interest.
2. 2. Selection of the order of the polynomial and quadratic approximations The selection of the order of the approximating polynomial and points x i for experimentation requires careful consideration. The degree of if(X) should be less than or equal to the
207 degree of g(X) to get a well-conditioned system of linear equations for the unknown coefficients. If ~(X) is of much higher degree than g(X), one obtains an ill-conditioned system of equations. Moreover, higher order polynomials can exhibit erratic behaviour in the sub-domains not covered by the experiments [9]. Of course, g(X) itself is not known. U p to a certain degree, a higher order polynomial improves the accuracy of the approximation at the expense of additional computation. The rate of increase in accuracy reduces with the degree of polynomial but the computational costs increase exponentially, as a higher order polynomial involves greater number of unknown coefficients and requires correspondingly more structural analyses. For reliability estimates, one needs to have a good approximation to g(X) around the design or minimum norm point, or the region of the failure domain .q~ where the joint probability density of X is relatively large and thus contributes most to the overall failure probability. Since we neither know the actual limit state function nor the actual design point, the accuracy of the reliability estimate depends on the accuracy of the polynomial approximation in the region of the design point. With an eye toward this accuracy as well as the high cost of repeated finite element analyses, a second degree polynomial is used in the present study. A second order response surface for n input variables is described by a quadratic model,
=A +xTB +xTcx
(3)
where A, B T = [B1, B2,..., nn] , and, C ~--
Cll
"'"
syn'l
Cln] Cnn]
define the undetermined coefficients. Experiments are conducted as per the adopted design and the resulting system of equations may be put in the form,
g(x) =De +e
(4)
where d is a vector of constants A, B i, Cij, the matrix D contains constant, linear, quadratic and cross-combinations of x i and e is the error vector, the components e i of which consist of a systematic or lack of fit error resulting from approximating g by if, and a pure experimental error, assumed to be a zero mean random vector. The solution,
E(d) = (DTD)-IDTg(x)
(5)
provides the expected values of the unknown coefficients. Other polynomial interpolation schemes using Lagrangian and Hermite polynomials are possible, although no specific examples of their application for reliability analysis could be located. At a higher level of sophistication, Ditlevsen [10] has presented a random field model for stochastic interpolation between point by point measured values of a spatially distributed material property.
2.3. Experimental designs In the absence of prior information, the experiments initially can be centered around the mean values of the variables. Points are located along the axes and elsewhere. The axial points
208
define the effect of each variable and are needed to determine coefficients for the terms involving only that variable i.e., x i, x 2 etc. The points in the plane x i x j define the interaction effects of variables x~ and xj. and are needed to determine coefficients for the cross terms, i.e., x~x i etc. When including the cross terms, a particular design of experiments depending on the nature of variables is discussed subsequently in the context of a specific response surface. The commonly used orthogonal experimental designs are 2 n and 3 n factorial and fractional factorial designs [11,12]. The 2 n factorial points for a simple two variable problem are indicated in Fig. 1 by , . Fractional factorials are useful when the number of variables is large and hence the number of experiments that can be conducted is less than the number of combinations in the full factorial set. Augmented 2 n factorials are obtained by adding n 2 points at the center of the design (total points = 2 ~ + n2). The center point is indicated by • in Fig. 1. Another efficient class of designs for fitting a second order response surface is the so-called central c o m p o s i t e design. This consists of a 2 ~ factorial design, augmented by n 2 > 1 central points and 2 ~ points placed at coordinates - a and + a along the axes. These axial points are indicated by ~9 along each axis in Fig. 1. Hence the design consists of (2" + 2n + n 2) points [13]. Other experimental designs include Simplex designs for first order models and Equiradial designs for second order surfaces. Khuri and Cornell [14] have discussed various experimental designs in detail. These factorial designs, though efficient, lead to unacceptably high computational efforts with the increase in number of variables for complex systems and may become more time consuming than simulation. An iterative response surface approach for reliability analysis was presented by Bucher and Bourgund [15]. The experimental design in each iteration consists of as many locations as the total number of undetermined coefficients in the polynomial i(X) =a
bixi+
+ i=1
Y'~cix2i i=1
(6)
in which xi, i = 1,..., n are basic variables and the parameters a, bi, Ci are to be determined. In order to obtain these constants, (2n + 1) experiments have to be conducted. Such an experimental design is referred to as fully saturated, and the surface fits exactly at the experimental points. The suggested experimental points for obtaining if(x) lie along axes x i (see Fig. 1). The points chosen are mean values tz i of the random variable X i , and xi = [1.1i, Jr h p i
X2
0
(
0
XI
Fig. 1. Location of experiments for a two variable problem.
209 in which h i is an arbitrary factor and o-i, is the standard deviation of X i. The constants are determined using (2n + 1) values of g ( x ) at these points. In the next step ~(x) is used along with information on/J,i and o-i to determine the minimum norm point x o for the surface ~ ( x ) = 0, based on the assumption that x are uncorrelated Gaussian variables. Once x o is found, g ( x o ) is evaluated and a new center point x M for interpolation is chosen on a straight line from mean vector ~ to x 9 so that g ( x ) = 0 from linear interpolation, i.e.,
XM= ~ +(Xo--~)[g(~l-----g(XO) ]
(7)
This strategy is to position the new center point reasonably close to g(x) = 0. A new surface obtained using center point x m may reasonably define the actual response in the failure domain. The total number of evaluations of g ( x ) is (4n + 3). Bucher and Bourgund [15] assert that this method of updating the polynomial results in a sufficiently accurate response surface. However, the writers believe that the accuracy depends on the characteristics of the limit state being explored and hence one cycle of updating may not always be sufficient. Key issues in this method are the selection of h i (which is arbitrary) and a suitable criterion to ensure that the fitted surface if(x) does indeed represent the actual surface in the failure domain reasonably well. The following section describes some ideas and approaches that address these issues.
3. Improvements to response surface methods
3.1. Criteria for fitting O(x) to g(x) The accuracy of the response surface obtained after two cycles of iteration and 4n + 3 experiments varies with each problem analysed, as shown subsequently, and thus the accuracy of the failure probability computed on the basis of the approximate response surface is unpredictable. One obvious question is how far should the process of interpolation and experimentation continue? Figure 2 illustrates the process of successive approximations, design points and center points for a two-variable limit state. The first approximation is obtained by conducting experiments at the mean values/z i denoted by XM1 in Fig. 2 and at (ix i + h i o i ) . The minimum norm point x m is determined for the fitted surface ~l(x) = 0. By linear interpolation, the new center point XM2 for experimentation is found from Equation 6 and a new surface ~ 2 ( x ) = 0 is then obtained. For this surface, again the minimum norm point x o 2 is obtained but the distance is measured from XM2 instead of from the mean /~(=Xul). This results in a low value for this distance, which is henceforth referred to as 8. ~-Fhe process is repeated until 6 becomes very small or zero. Then the factor h i is reduced and the process is repeated until 6 becomes zero again. This can be achieved in about three to four iterations, resulting in a very good approximation for the actual surface and a value of x 9 that is close to the actual minimum norm point. However, if the process is repeated indefinitely with progressively smaller values of h i , then at some stage one obtains an ill-conditioned system of equations. Experience has shown that starting the process with h i = 2.0 and then repeated once more with h i = 1.0 gives a satisfactory
210 LEGEND : X t, X= = A two v ~ i a b l e
..... ....
X2 'x' ~'-~g.IL J
space
for
g(x) m Actual Limit State g t ( x ) =, Raeponee S u r f , a c e gel x ] = Ra=ponse SuP f a c e g~l(x] = Reeponsa S u r f , a c e -
illustration
Function Fir'st Approximation Second A p p r o x t mat t on Third Approximation
xul - Center Point t,or the t n approximation x m - D e s i g n P o i n t t,or t h e i th a p p r o x i m a t i o n +g - Ihnasur'ed C o n v e r g e n c e P a r a m e t e r t ' o r t h e i ~h a p p r o x i m a t i o n
9(x)
-
xui
X1
Fig. 2. Convergence to the minimum norm point: Approach A-0.
response surface around the minimum norm point. Once the surface is obtained, advanced simulation techniques can be used to estimate the failure probability [7,16]. This procedure is referred to as Approach A-0 in the examples to follow in Section 4. When the variables are dependent a n d / o r not Gaussian, they can be transformed into an independent standard normal space and the procedure as outlined above can be carried out to ensure the criterion based on ~. Experience has shown that this improves the conditioning of the system of equations as well.
3.2. Experiments at distribution extremes The tail regions of the probability distributions of random load and strength are most important in the estimation of failure probabilities. Hence, to obtain an accurate fit to the actual surface in the vicinity of the design point, one might think that experimental points should be chosen in the lower tails of the resistance variables and upper tails of the load variables (Figs. 3 and 4) rather than over their entire range. Three different approaches incorporating experiments at distribution extremes are considered in this section. t'.( x 2 )
f x ( Xl )
12345
Xl
12545
Fig. 3. Experiments at distribution extremes: Approach A-1.
~ x2
211
f,( xl )
f',,( ×2]
k,tktt~t --- x I
X2
Fig. 4. Experiments at distribution extremes: Approach A-2.
In the first approach (Fig. 3), the width of the sampling region for each variable is kept sufficiently large and in the first iteration, as many points are chosen in this region as the number of experiments. The experiments are conducted at points whose coordinates are rank ordered locations in the tail region of each variable. Thus the first approximation is obtained and then the minimum norm point for the first approximation is determined. For the second iteration, the width of the sampling region is reduced around this point and the distance 8 is measured as before (Section 3.1). At low values of 8, a fairly accurate response surface results, as shown subsequently in Example 1 in Section 4.1. This approach is designated henceforth as Approach A-1. In the second approach (Fig. 4), the effect of shape of the probability distributions of variables on the selection of experimental points is considered. Using a two-variable problem as an illustration, the center point for the first iteration is chosen at x = (x 1, x 2) where x 1 and x 2 are selected at k ~ and k ~ fractiles (e.g., typically the 95th percentile of load or 5th percentile of resistance) of random variables X 1 and X 2 (Fig. 4). Two other points on either side of both k n and k12 are also chosen. The experiments are conducted at ( k l l , k12) , (kol , k12) , (k21 , k12) , (k n, k02) and ( k l l , k22). These points are at the k ~ fractiles of the respective distribution functions. Once the minimum norm point for the first trial surface is obtained, sampling is concentrated around that point and the process is repeated as before. The same logic is applicable for an n-variable problem. This approach is designated henceforth as Approach A-2. As a variant of this approach, points on either side of the k ~ fractiles of X i can be chosen as a constant times its standard deviation, as done in Approach A-0, i.e., the experiments are conducted at ( k l l , k12) , ( k l l -t- hlOrl, k12) and ( k l l , k12 -k- h20"2). However, h i would typically be less than what it is when experimenting around the mean but large enough to have a well-conditioned system of equations. This approach is henceforth designated as Approach A-3. The objective in procedures A-1 through A-3 is to include distribution information in the selection of experimental points. Example 4.1 in Section 4 shows that fairly accurate response surfaces can be generated from these procedures.
3.3. Addition of terms to quadratic O(x) The accuracy of the response surface may be improved by inclusion of the cross terms in the polynomial [15,17] as follows, n ff(x)=a+
n
Ebixi.-]- E i=1
~
cijxixj
i=1 j>~i=l
(8)
212 TABLE 1 Quadratic combinations for experimental locations Resistance(Xi)* Load(Xj) Resistance(Xi)* Resistance(Xj) Load( X i ) * Resistance(Xj) Load(X/) * Load(X~)
- ( x i ) . +(x~) - (xi)* - (x j) + ( x ) . - (xj) +(x).(xp
The total n u m b e r of experiments to be conducted for each approximation would then increase to (n + 1 ) x (n + 2)/2. To gain some insight on the appropriate location of experiments, consider a two variable limit state function approximation ~(x 1, x2). There is only one cross t e r m (clzxlX 2) in this polynomial which has six constants in all. The first five locations are taken at the m e a n and four axial points. The last location for the cross term has to be selected among the four choices indicated by * in the Fig. 1. The best location depends on w h e t h e r X~ and Xj are "resistance" or "load" variables, as summarized in Table 1. The same logic should be applied for a cross term involving a greater n u m b e r of variables, by selecting points whose coordinates lie on the negative side of the "resistance" axes and positive side of the " l o a d " axes, with the origin of the axes taken at the center point. Such a combination results in the best possible fit in the failure domain among all possible combinations. This is illustrated in Example 4.2 in Section 4. The inclusion of cross terms in the approximation increases the number of experiments from (2n + 1) to (n + 1) x (n + 2)/2, an increase of n x (n - 1 ) / 2 which is significant if the reliability analysis involves a large n u m b e r of variables. However, totally neglecting the cross terms may be inappropriate as there are situations where cross terms are more important to the actual response than the squared terms of the individual variables. As a simple example to illustrate this point, the limit state of a truss m e m b e r subjected to axial force is defined by A x f y - P = 0, A, fy, and P being the area of cross section, yield stress and applied force respectively; it is obvious that a better response surface results by taking cross terms than squared terms. For efficiency, only those terms which are important in terms of accuracy should be included. As the iterations proceed, terms that are relatively unimportant can be neglected, reducing the number of experiments at each cycle. If gc is the response computed at the center point and gi is the response computed at the ith experimental point, then the square of the difference normalised with respect to gc gives the relative importance of the ith location. If,
( gi-gc l2< c. gc
(9)
]
then the corresponding term may be neglected in subsequent iterations. The neglected term may be a linear or a square or a cross term; this is more realistic and accurate than just neglecting the cross terms and should lead to better subsequent approximations. This is illustrated in Example 4.2 in the following section.
4. Numerical examples Two numerical examples are presented in this section. The cantilever beam in the first example facilitates visualization of the accuracy of successive approximations by response
213 surfaces of the actual limit state. In the second example, the non-linear dynamic response of a single degree of freedom system with random parameters is analysed and compared with a previously published response surface analysis.
4. 1. Cantilever beam A cantilever beam with a rectangular cross section subjected to uniformly distributed loading is considered. The limit state of serviceability with respect to maximum deflection at the free end being greater than 1/325, is defined by (W X b ) l 4 g =
8EI
l
+ 325
(10)
where w, b, l, E and I are respectively load per unit area, width, span, modulus of elasticity and moment of inertia of the cross section. The random variables are the load per unit area and the depth of the cross section. Assuming E and l as deterministic and equal to 2.6 × 104 MPa and 6 m respectively, the limit state equation reduces to g ( X ) = 18.46154 - 7.476923 X 101° X---~= 0 X2
(11)
where X 1 is the load in MPa and X 2 is the depth in mm. The load is assumed to be normally distributed, with a mean ~.L1 1000 N / m E and a coefficient of variation V1 -- 0.2. The beam depth is also normally distributed, with ~2 = 250 mm and V2 = 0.15. The probability of failure for this limit state is calculated according to the various procedures described in earlier sections. The "exact" solution using importance sampling with 1000 simulations is 9.6071 × 10 -3. Advanced first order reliability analysis gives a reliability index of 2.33 corresponding to a probability of failure of 9.9031 × 10 -3. The "exact" design point is located in Figures 5-9 that follow. Consistent with eqn. (11) the response surface is assumed to have the form of Equation 5, with no cross terms, and is obtained for this limit state by carrying out analyses along axial points(see Fig. 1) Curves 1 and 2 in Fig. 5 represent Bucher's algorithm with h i = 2.0. Simulation using curve 2 gives a probability of failure of 13.7538 × 10 -3. Following Approach =
8
I
I
I
7N
x
I
I
I
I
I
I
I
I
....................
I
.
I
I
~
5 - ' ~ 4
3 2
f ["
~ 0.0
Exit
C~esents
Response S u r f a c e
Limit
State
o f I th I t e r a t i o n Xo represents "exact" des1 gn p o t n t
1
0
C ~ e n t s
~ 0.2
i
i 0.4
i
i 0.6
I
l 0.8
I
I 1.0
l
i 1.2
X1/~1 F i g . 5. C o n v e r g e n c e - - A p p r o a c h
A-0 - h i = 2.0.
l
I 1.4
214 I
0.8
I
I
I
I
I
I
I
I
I
I
I
I
I
0.70,6fw :1
X
0,50,40.50.2-
LEGEND
f •
C ~ e n t s Exact L i ~ t State C ~ e s e n t s Response Surface of I ~ I t e r a t i o n represents "exact" design potnt
0.1 0.0
I 0.0
t
I
0.2
I 0.4
I
I 0.6
I
I 0.8
I
I
I
I
1.0
I
1.2
I 1.4
xl/~x Fig. 6. Convergence--ApproachA-0
-
hi
=
1.0.
A-0 to improve the response surface, successive surfaces are obtained with h i = 2.0 (curves 3 and 4 in Fig. 5). Fig. 6 shows the surfaces obtained in the second stage of iteration with h i = 1.0. Table 2 gives the distance 6, between successive center and design points. The probability of failure obtained from simulations on the final response surface (curve 6) is 9.5410 x 10-3which is very close to the exact solution. Thus Approach A-0 yields a very close result; the n u m b e r of iterations required to fit the response surface to g(x) = 0 is limit state d e p e n d e n t and cannot be generalised as two or three cycles. Progressively smaller values of h i lead to an ill-conditioned system of equations. Table 3 illustrates the effect on the solution when the response surfaces are generated by locating experiments in the tail regions of the random variables using Approach A-1. The process is initiated by selecting a relatively large range of 0.75 to 0.95 fractiles for the load variable and 0.05 to 0.25 fractiles for the resistance variable. Once the design point is obtained, the range is successively reduced until the distance 6 becomes very small. Fig. 7 illustrates the successive approximations; sampling using curve 3 of Fig. 7 yields a probability of failure of 9.6398 x 10 -3 which is very close to the exact solution of 9.6071 x 10 -3. It can be seen from Fig. 6 and Fig. 7 that this technique yields a response surface with a better fit to the exact limit state function over a wider range compared to that obtained by Approach A-0 (Table 2). The drawback of this approach lies in the selection of experimental points for the first iteration. With an increase in the n u m b e r of variables, it becomes increasingly difficult to choose the range for each variable.
TABLE 2 Successive Response Surface fits to
g(x)= 0 using Approach A-0 (cf Figs. 5 and 6)
Iteration
hi
6
XloiN/m2
X2Dimm
1 2 3 4 5 6
2.0 2.0 2.0 2.0 1.0 1.0
2.6436 0.2868 0.0169 0.0002 0.0002 0.0000
1069.29 1064.15 1073.15 1072.40 1072.38 1072.38
151.72 168.36 162.29 163.52 163.15 163.15
Experiments at XMg and
(XMg +_hio"i)
215 TABLE 3 Experiments at distribution extremes: Approach A-1 (cf Fig. 7) Iteration
Variable
kI
k2
k3
k4
1
x~ x2 x1 x2 xI x2
0.75 0.05 0.78 0.0050 0.7365 0.0120
0.80 0.10 0.83 0.0075 0.7865 0.0145
0.85 0.15 0.88 0.0100 0.8365 0.0170
0.90 0.20 0.93 0.0125 0.8865 0.0195
2 3
ks 0.95 0.25 0.98 0.0150 0.9365 0.0220
kdesign
t~
0.8810 0.0107 0.8365 0.0170 0.8413 0.0158
1.51 0.3327 0.0456
Fig. 8 shows plots of successive approximations obtained from Approach A-2 of Section 3.2, incorporating the shape of the probability distributions in selecting experimental points in the tail region. C o m p a r e d to Figs. 5-7, the response surfaces fit well over a narrower range. Table 4 gives the details of two iterations; the probability of failure obtained is 11.1508 × 10 -3 which is 16% higher than the exact value. Table 5 and Fig. 9 show details and plots of successive
.8
°.,
0
.
o.st 20.4
7
.....
~
......
....
t~ 0"-~JI- / - - I /"
C ~ ~ n t s Exact Limit State Curve / r e p r e s e n t s Response S u r f a c e of i ~ Ite;etion " X0 r e p r e s e n t s " e x a c t " d e s | g n , p o t n t
0 , 1 --'~
0
. 0.2
0.0
0 0.4
~ 0.6
0.8
1. 0
1.2
1.4
xt/~l Fig. 7. Convergence--Approach A-1.
I
0.8
I
I
I
I
I
I
I
I
I
I
I
I
I
0.70.6-
- . ....... ~ .....
r
NO.5:l ~O.4x 0.5-
=
:
=
:
=
~
~
Curve E r e p r e s e n t s Exact L i m i t S t a t e Curve I r e p r e s e n t s R e s p o n s e ~ J r f a e e of I m I t e r a t i o n represents "exect" design p o | n t
0.20.1 0.0 0.0
:
I
I
0.2
I
I
0.4
I
I
0.6
I
I
0.8 Xl/~l
I
I
I
1.0
Fig. 8. Convergence--ApproachA-2.
I
1.2
I
I
1.4
216 TABLE 4 Experiments at distribution extremes: Approach A-2 (cf Fig. 8) Iteration
Fractiles selected
Experiments at kth fractile
Design point at kth fractile
6
1
x a ~ (0.7,0.8,0.9) x 2 ---,(0.01,0.1,0.2)
(0.87,0.017)
1.0778
2
x 1 ~ (0.80,0.87,0.94) x 2 ~ (0.01,0.017,0.024)
(0.8,0.1) (0.8,0.01),(0.8,0.2) (0.7,0.1),(0.9,0.1) (0.87,0.017) (0.87,0.01),(0.87,0.024) (0.80,0.017),(0.94,0.017)
(0.87,0.017)
0.0003
TABLE 5 Experiments at distribution extremes: Approach A-3 (cf Fig. 9) Iteration
Mean point for experiment as a fraction of actual mean, x M / ~
Design point as a fraction of actual mean, x D/tz
6
1
(1.256,0.808) corresponding to (0.9,0.1) fractile (1.0673,0.6695) (1.0774,0.6559) (1.0788,0.6542) (1.0790,0.6540)
(1.0673,0.6695)
2.2291
(1.0774,0.6559) (1.0788,0.6542) (1.0790,0.6540) (1.0790,0.6540)
0.1438 0.0181 0.0026 0.0003
2 3 4 5
Experiments at (kl, k2), (k 1 + or1, k 2) and (kl, k 2 + ~2)
a p p r o x i m a t i o n s o b t a i n e d f r o m A p p r o a c h A-3 l o c a t i n g e x p e r i m e n t s at k i a n d ki _+ tri in t h e tail region. T h i s gave a p r o b a b i l i t y o f f a i l u r e o f 9.5410 × 10 -3 w h i c h c o i n c i d e d with t h a t o b t a i n e d by A p p r o a c h A - 0 ( T a b l e 2). T h e results o b t a i n e d f r o m v a r i o u s p r o c e d u r e s a r e s u m m a r i s e d in T a b l e 6. It is e v i d e n t f r o m the a b o v e discussion t h a t i n c o r p o r a t i n g i n f o r m a t i o n o n t h e p r o b a b i l i t y d i s t r i b u t i o n s o f t h e basic
:
0.8 0.7
I
:
I ..-~
0.6
I
I
I
I
I
I
I
I
......................................
I
I
xD
~0.5 0.4 x
0,5
//
0.2- / O,l I 0.0 0.0
C
~
ents Response Surface
e
0 f I~h I t e r a t i O n X0 r e p r e s e n t s
0.2
0.4
0.6
.
.
.
.
.
"exa(: " t d e s t g n p o i n t
0.8
1.0
Xl/l~t
Fig. 9. Convergence--Approach A-3.
1.2
1.4
217 TABLE 6
Comparison of various methods Methods
Probability of failure
Simulation using actual limit state model Advanced First Order Second Moment Method Bucher's Method Approach A-0 Approach A-1 Approach A-2 Approach A-3
0.0096071 0.0099031 0.0137538 0.0095410 0.0096398 0.0111508 0.0095410
variables in the selection of experimental points at which the finite element analysis for complex structures is carried out yielded fairly accurate response surfaces but did not save computational effort required to obtain a good approximation (measured, for comparable level of accuracy, by 6) to the probability of failure. Hence, it can be concluded that sampling in the tail regions of the random variables instead of the entire region does not lead to significant improvements in the response surfaces or reliability estimates. 4.2. Dynamic response of a nonlinear oscillator A non-linear undamped single degree of freedom system analysed by Bucher et al. [15,17] is illustrated in Fig. 10. The limit state is defined by
g(x) =
3r-
I Zmax I
(12)
in which Zmax is the maximum displacement response of the system and r is the displacement at which one of the springs yields. The statistical parameters of the basic variables are listed in Table 7. All variables are assumed to be independent and lognormally distributed.
I
.
C l
~
z(t)
_ ~ _ ~ . F( t ) r
C2 &l
Restoring
r
force
~
z
il
F(t)
tl
~t
Fig. 10. Nonlinear SDOF oscillator with nonlinear restoring force and pulse load for Example 2.
218 TABLE 7 Statistical data for Example 2 Number
Variable
Mean value
Standard deviation
1 2 3 4 5 6
m c1 c2 r F1 t1
1.00 1.00 0.10 0.50 1.00 1.00
0.05 0.10 0.01 0.05 0.20 0.20
Importance sampling with 1000 simulations of nonlinear dynamic response yielded a probability of failure of 3.99864 x 10 -2. Bucher et al. [15,17] obtained a probability of failure of 3.81 x 10 -2 by adaptive sampling using 1500 simulations, 3.01 x 10 -2 using a response surface without cross terms, and 3.67 x 10 -2 using a response surface with cross terms. Approach A-0, using a response surface without cross terms and with just one more cycle, gives a failure probability of 3.93365 x 10 -2. Table 8 gives details of the iterations. The effect of cross terms are considered for this example only for the first iteration as it was found that a good approximation could be obtained in the second iteration without them. W h e t h e r the cross terms may be neglected altogether depends on the complexity of the problem and the nature of the unknown limit state. Without cross terms, the first response surface approximation to the limit state gives a probability of failure of 2.61898 x 10 -2. The experimental points for cross terms are then selected in the first quadrant (see Fig. 1) and the resulting failure probability associated with the revised response surface is 3.22473 x 10 -2. In order to study the effect of location of experiments, the procedure is repeated with experimental points in the second, third and fourth quadrants. With the minimum norm distance for the approximating surfaces taken as an index of the accuracy, the result is almost same regardless of the quadrant in which the points are located, as indicated in Table 9. The suggested procedure of selecting points on the negative side of the axes for "resistance" variables and positive side for " l o a d " variables with the mean as the origin gives the best result in terms of both minimum fl and closest estimate of the failure probability. The probability of failure is 3.42005 x 10 -2 in the first iteration, which is quite comparable to the exact value of 3.99864 x 10 -2 The importance of the individual terms in the polynomial response surface is evaluated using eqn. (9). The proposed index is calculated for each term and terms for which the squared error is less than 0.01 (E) are neglected. In all, five terms (out of 21) are neglected. Treating this as a first iteration, a new surface is fitted which gives an almost identical failure probability of
TABLE 8 Results of Approach A-0 for Example 2 Iteration 1 2 3 Non linear dynamic response
6 2.08001 0.14650 0.09543
Probability of failure 2.61898 x 10 - 2 4.09679 x 10 -2 3.93365 x 10 - 2 3.99864 x 10 - 2
219 TABLE 9 First iteration accuracy Method
/3
Probability of failure by simulation
Without cross terms Cross terms (x i, x j) in I quadrant ( + , + ) II quadrant ( - , + ) III quadrant ( - , - ) IV quadrant ( + , - ) Resistance ( - ) and Load ( + ) With reduced crossterms Simulation by non-linear analysis
2.08001
2.61898 x 10 2
1.90273 1.89055 1.90991 1.91874 1.88704
3.22473 x 10- 2
1.88882
3.40402 x 10- 2
3.42005 x 10 - 2
3.99864 x 10-2
3.40402 x 10 -2. Thus the terms neglected do not have a significant effect on the response surface or probability of failure. The results are summarised in Table 9. The accuracy of the response surface obtained depends on the characteristics of the actual limit state and is affected by the selection of points at which the experiments are performed. The proposed method (Approach A-0) of continuing the response surface fitting process until the distance 6 becomes very small or zero is a reasonable way to tackle reliability analysis of complex structural systems. On the basis of the limited studies made, the suggested location of experimental points for finite element analysis seems to be logical. The index in eqn. (9) is sufficient to test the significance of each term after the first iteration so that unnecessary experiments can be reduced in the subsequent iterations.
5. Summary and conclusions An adaptive iterative procedure was presented to develop a response surface that is locally a good approximation to the actual limit state surface in the region of maximum joint probability density and can be used for structural reliability analysis. Methods were presented to enable inclusion of information on probability distributions of random variables in selecting experimental points used to generate response surfaces. However, it was found that such methods have no particular advantages with regard to computational effort and accuracy as compared to more conventional methods for selecting experimental points. The location of experiments for cross terms and a systematic procedure to reduce the number of experiments depending upon the importance of each term in the previous approximation also were suggested. It is felt that the response surface method can be used efficiently for reliability analysis of certain structural systems where behavioral models to describe their various limit states cannot be developed in closed form.
6. Acknowledgements The first author expresses special thanks to Professor A.S.R. Sai at his parent institution, the Indian Institute of Technology, Kanpur, India, for initiating him into research in the field of
220
reliability. He also expresses grateful acknowledgement to the Council for International Exchange of Scholars, Washington D.C., for sponsoring his research at the Johns Hopkins University.
References 1 R.E. Melchers, Structural Reliability : Analysis and Prediction, Ellis Horwood Limited, Chichester, 1987. 2 A. Der Kiureghian, and J.B. Ke, The Stochastic Finite Element Method in Structural Reliability, Stochastic Structural Mechanics, Lecture Notes in Engineering 31, Springer, Berlin, 1987. 3 A.M. Hasofer and N.C. Lind, Exact and invariant second moment code format, J. Engng. Mech. Div., 100 (EM1) (1974) 111-121. 4 R. Rackwitz and B. Fiessler, Structural reliability under combined random load sequences, Computers and Structures, 9 (1978) 490-494. 5 O. Ditlevsen, Principle of normal tail approximation, J. Engng. Mech. Div., 107 (EM6) (1981) 1191-1208. 6 M. Shinozuka, Basic analysis of structural safety, J. Struct. Div., 109 (3) (1983) 721-740. 7 R.E. Melchers, Importance sampling in structural systems, Structural Safety, 6 (1989) 3-10. 8 G.I. Schueiler and R. Stix, A critical appraisal of methods to determine failure probabilities, Structural Safety, 4(4) (1987) 293-309. 9 S. Engelund and R. Rackwitz, Experiences with experimental design schemes for failure surface estimation and reliability, Proceedings of the Sixth Speciality Conference on Probabilistic Mechanics and Structural and Geotechnical Reliability, ASCE, 1992, pp. 252-255. 10 O. Ditlevsen, Random field interpolation between point by point measured properties, Proceedings, 1st International Conference on Computational Stochastic Mechanics, Corfu, Greece, 17-19 September 1991. 11 R.H. Myers, Response Surface Methodology, Allyn and Bacon, Boston, 1971. 12 R.G. Petersen, Design and Analysis of Experiments, Marcel Dekker, New York, 1985. 13 L. Faravelli, Response surface approach for reliability analysis, J. of Engng. Mech., 115 (12) (December 1989) 2763-2781. 14 A.I. Khuri and J.A. Cornell, Response Surfaces, Designs and Analyses, Statistics Textbooks and Monographs series, Volume 81, Marcel Dekker, New York, 1987. 15 C.G. Bucher and U. Bourgund, A fast and efficient response surface approach for structural reliability problems, Structural Safety, 7 (1990) 57-66. 16 C.G. Bucher, Adaptive sampling--an iterative fast Monte Carlo procedure, Structural Safety, 5(3) (1988) 119-126. 17 C.G. Bucher, Y.M. Chen and G.I. Schueller, Time variant reliability analysis utilising response surface approach, in: P. Thoft-Christensen, ed., Reliability and Optimization of Structural Systems '88, (Lecture Notes in Engineering 48), Springer, Berlin, 1989, pp. 1-14.