- Email: [email protected]
Reliability analysis using object-oriented constrained optimization
Structural Safety 26 (2004) 69–89 www.elsevier.com/locate/strusafe
Reliability analysis using object-oriented constrained optimization Bak K. Lowa,*, Wilson H. Tangb a
School of Civil and Environmental Engineering, Nanyang Technological University, Block N1, #01b-40, 50 Nanyang Avenue, Singapore 639798, Singapore b Department of Civil and Structural Engineering, Hong Kong University of Science & Technology, Clear Water Bay, Kowloon, Hong Kong
Abstract A practical procedure is presented for reliability analysis involving correlated nonnormals. Equivalent normal transformation is performed as in the first order reliability method; the correlation matrix is used as it is without orthogonal transformation. Performance functions appear as simple spreadsheet cell objects but may embody substantial user-created program codes. The multidimensional equivalent dispersion ellipsoid is likewise a single cell object that contains various spreadsheet routines for equivalent normal transformation and matrix operations. The procedure can be viewed as the constrained optimization of the equivalent hyperellipsoid (yielding the reliability index) in the original space, and is automatic and efficient in the ubiquitous spreadsheet platform. Different probability distributions can be selected at ease and conveniently transformed to equivalent normals using a short program code given in the paper. Illustrations are presented using a three-variate simple performance function and various combinations of ten common distributions. This is followed by two relatively complicated examples, namely an asymmetrically loaded beam on Winkler medium, and a complex strut with implicit performance function. The probabilities of failure inferred from reliability indices are compared with Monte Carlo simulations. The simplicity, transparency and flexibilities of the object-oriented constrained optimization approach are demonstrated. # 2003 Elsevier Ltd. All rights reserved.
1. Introduction The Hasofer–Lind [1] index and the first order reliability method (FORM) have been welldocumented in Ditlevsen [2], Shinozuka [3], Ang and Tang [4], Madsen et al. [5], and Tichy [6], * Corresponding author. Tel.: +65-6790-5270; fax: +65-6791-0676. E-mail address: [email protected] (B.K. Low). 0167-4730/03/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0167-4730(03)00023-7
70
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
for example. The continuing interest in FORM is also evident in Melchers [7], Haldar and Mahadevan [8], Nowak and Collins [9]. The potential inadequacies of FORM in some special cases have been recognized, and more refined alternatives proposed, in Der Kiureghian et al. [10], Chen and Lind [11], Wu and Wirsching [12], and Zhao and Ono [13,14]. On the other hand, Rackwitz [15] noted the usefulness and accuracy of FORM in most applications; some counterexamples were summarized, their rather extreme nature and artificiality noted, but their healthy role also admitted. With respect to the Hasofer–Lind index and FORM, a practical procedure was presented in Low and Tang [16], based on the perspective of an expanding equivalent ellipsoid tangent to the limit state surface in the original space of the random variables. The Rackwitz–Fiessler [17] equivalent normal transformation was used, but the concepts of coordinate transformation and frame-of-reference rotation were not required. Correlation was accounted for by setting up the quadratic form directly. Iterative searching and partial derivatives were automatic using constrained optimization in the ubiquitous spreadsheet platform. The present paper extends Low and Tang [16], by testing robustness and accuracy for various nonnormal distributions and more complicated performance functions, and by providing enhanced operational advantage and versatility. A simple 3-variate case in the next section elucidates the principles of object-oriented constrained optimization, and prepares the background for the more efficient and flexible strategy that follows. Different combinations of ten common distributions are used. Two relatively complicated cases are also analyzed, namely an asymmetrically loaded beam on Winkler elastic medium, and a complex strut the performance function of which requires iterative program calculations. The results obtained by the spreadsheet object-oriented constrained optimization approach are compared with Monte Carlo simulations.
2. Constrained optimization of dispersion ellipsoid in original space The matrix formulation of the Hasofer–Lind index is: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ min ðx mÞT C1 ðx mÞ x2F
ð1aÞ
or, equivalently: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi xi mi T 1 xi mi ¼ min ½R x2F i i
ð1bÞ
where x is a vector representing the set of random variables, m the mean values, C the covariance matrix, R the correlation matrix, and F the failure domain. The index is, in classical explanation, the minimum distance from the mean value point to the limit state surface, in units of directional standard deviations. An alternative interpretation was presented in Low and Tang [16], where the perspective of an expanding ellipsoid (Fig. 1) led to a simple method of computing the Hasofer– Lind index in the original space of the variables.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
71
Fig. 1. Design point and equivalent normal dispersion ellipsoids illustrated in the plane.
One may note that the quadratic form in Eq. (1a), ðx mÞT C1 ðx mÞ; appears also as a negative exponent of the multivariate normal probability density function. Hence, to minimize (or 2 in the multivariate normal density function) is to maximize the value of the multivariate normal probability density function, and to find the smallest ellipsoid tangent to the limit state surface is equivalent to finding the most probable failure point. This perspective is consistent with Shinozuka [3] that ‘‘the design point x* is the point of maximum likelihood if x is Gaussian, whether or not its components are uncorrelated.’’ When the random variables are nonnormals and correlated, the approach to be presented below uses the Rackwitz–Fiessler [17] equivalent normal transformation, but without the need to diagonalize the correlation matrix (which is equivalent to rotating the frame of reference). Also, the iterative computations of the equivalent normal mean (mN) and equivalent normal standard deviation ( N) for each trial design point are automatic during the constrained optimization search.
3. Three-variate correlated nonnormals: a simple example The case in Fig. 4 of Low and Tang [16] for correlated nonnormals is shown in Fig. 2, with a more succinct set-up. The random variables Y and Z are correlated lognormals, and M is Type 1
72
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Fig. 2. Reliability analysis of correlated nonnormals using constrained optimization. Random variables Y and Z are correlated lognormals; M is Type 1 Extreme (Gumbel).
Extreme (Gumbel). The performance function is g(x)=Y*ZM. The original mean () and standard deviations () of Y, Z and M are (40, 5), (50, 2.5), and (1000, 200), respectively. The square root of the quadratic form in Eq. (1b) is computed directly using array formula: ‘‘ ¼ sqrtðmmultðtransposeðnxÞ; mmultðminverseðcrmatÞ; nxÞÞÞ”
ð2Þ
(Array formulas are created by pressing ‘‘Enter’’ while holding down ‘‘Ctrl’’+‘‘Shift’’). In Eq. (2), ‘‘mmult’’, ‘‘transpose’’ and ‘‘minverse’’ are Excel’s built-in functions, each being a container of VBA codes for matrix operations. The nx vector in Eq. (2) contains functions (via mN and N) of Excel routines ‘‘normsinv’’ (cell D20), which returns the inverse of the standard normal cumulative distribution, and ‘‘normdist’’ (cell E7). These two functions are used to transform the Type1 Extreme variate M into equivalent normal distribution, based on Rackwitz– Fiessler [17] two-parameter equivalent normal transformation: F1 ½F ðxÞ N ð3Þ Equivalent normal standard deviation: ¼ f ðxÞ Equivalent normal mean :
mN ¼ x N F1 ½F ðxÞ
ð4Þ
where x is the original nonnormal variate, F1 ½: is the inverse of the cumulative probability (CDF) of a standard normal distribution, F(x) is the original nonnormal CDF evaluated at x,
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
73
{.} is the probability density function (pdf) of the standard normal distribution, and f(x) is the original nonnormal probability density ordinates at x. Therefore, the single spreadsheet object for (cell N7) is a container of several Excel function objects and substantial program codes. For lognormals Y and Z, closed form equivalent normal transformation is available and has been used, as shown by the equations for mN and N of Y and Z in Fig. 2. For Type 1 Extreme variable M, Excel’s built-in routines normsinv (for 1) and normdist (for {.}) have been used to automatically update mN and N (cells D7 and E7) by Eqs. (3) and (4) during the iterative constrained optimizational search for the design point on the limit state surface. The design point (x* values) was obtained by calling Excel’s built-in optimization program Solver, to minimize the cell, by changing the x* values, subject to the constraint that the performance function g(x)=0. Prior to Solver search, the x* values had values equal to the mean values () of the original random variables. Iterative numerical derivatives and directional search for the design point x* were automatic. The value of the reliability index in Fig. 2 implies a probability of failure approximately equal to 0.00385. This compares favorably with a failure rate of 0.00382 from Monte Carlo simulation using the commercial software @RISK 2001 (http://www.palisade.com). The sample size was 200,000. The @RISK distribution functions used were RiskLogNorm(40,5) for variable Y, RiskLogNorm(50, 2.5) for variable Z, and RiskExtvalue(u, 1/) for M, i.e., RiskExtvalue(910, 155.9). The correlation matrix was modeled. Fig. 2 is concise, instructive and relatively transparent for computing the Hasofer–Lind index involving correlated nonnormals. It is ideal for pedagogic purposes. Nevertheless, a clearer and more convenient user interface is desirable so that various nonnormal distributions can be selected at ease without the need to modify the template for different probability distributions. This is developed in the next section. For reliability index computations in the following sections, Microsoft Excel XP (2002) and its built-in optimization program have been used throughout. Alternative optimization programs (e.g., http://www.Solver.com and http://www.palisade.com/ html/evolver.html, including genetic and evolutionary algorithms, are available (at a cost) as Excel add-ins. If desired, these spreadsheet add-ins may be used instead of the standard Excel Solver. The procedures and codes in the next section are unaffected by the type of optimization add-ins used. Should an alternative spreadsheet platform evolve in future, it is simple to modify the short code (Fig. 4 below) to conform to the language syntax of the new platform. For the examples below, Excel XP with its built-in standard Solver is adequate, and has the advantage of being ubiquitous.
4. Reliability analysis using object-oriented constrained optimisation The more succinct and versatile alternative to Fig. 2 is shown in Fig. 3. The key feature is that all distribution-specific equations are relegated to a short user-created (or user-defined, in VBA parlance) program code shown in Fig. 4, which is called by cells beneath the headings ‘‘mN’’ and ‘‘ N’’ in Fig. 3. The VBA (Visual Basic for Applications) program code in Fig. 4 is created via Tools/Macro/VisualBasicEditor/Insert/Module on the Excel worksheet menu. It is easily transferable to other workbooks.
74
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The main program is EqvN, which returns the equivalent normal mean (if code=1) or the equivalent normal standard deviation (if code=2). Its arguments include ‘‘DistributionName’’ and ‘‘Paralist’’, the latter being a vector of up to 4 parameters, depending on the type of probability distribution. In the fourth line of the program in Fig. 4, the ‘‘Select Case’’ control structure alters the flow of execution to one of several code segments, depending on the input ‘‘DistributionName’’. The functions Trim and UCase are VBA’s built-in functions: ‘‘Trim’’ returns a text with both leading and trailing spaces removed, ‘‘UCase’’ then renders the text uppercase throughout. Function EqvTransform in Fig. 4, called by six distribution types (from EXTVALUE1 to BETADIST) in the main function EqvN, implements the Rackwitz-Fiessler equivalent normal transformation according to Eqs. (3) and (4), based on the cumulative probability distribution (CDF) and probability density function (pdf) of the original nonnormals at x. Function betapdf in Fig. 4, called by BETADIST in the main program EqvN, returns the value of the density function of the beta distribution. The three functions in Fig. 4 make liberal calls to Excel’s built-in functions, using the syntax ‘‘Application.GammaDist(. . .)’’, Application. Weibull(. . .)’’, ‘‘Application.NormDist(.)’’, ‘‘Application.NormSInv(.)’’, where ‘‘Application’’ stands for Excel. The use of these Excel ‘‘objects’’ (each a container of VBA program codes) results in much simplicity, clarity, and brevity of the program codes and structure. The Uniform (or rectangular) distribution can be obtained [18] as a special case of BetaDist (1,1,min,max). If desired, it can also be created (only a few statement lines are required), and be one of the ‘‘Cases’’ in Fig. 4.
Fig. 3. Constrained optimization of equivalent normal ellipsoid in the original space of correlated nonnormals. This approach is more user-friendly and versatile than Fig. 2.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
75
Fig. 4. User-defined VBA codes. These functions were used in cells J5:K7 of Fig. 3 for equivalent normal transformation.
76
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The PERT distribution (an option in @RISK; also described in Vose [19]) is also a special case of Beta distribution. It is defined by the same three parameters (min, mode, max) as the Triangular distribution, but has some advantages over the latter. The PERT distribution shown below can be inserted to the code in Fig. 4: Case ‘‘PERTDIST’’: a=para1: Mode=para2: c=para3 If x<=a Then x=a+del If x>=c Then x=c-del mean=(a+4* Mode+c)/6 If Mode=mean Then f=6 Else f=(2* Modeac)/(Mode-mean) a1=(meana)*f/(ca): a2=a1*(cmean)/(meana) CDF=Application.BetaDist(x, a1, a2, a, c): pdf=betapdf(x, a1, a2, a, c) EqvN=EqvTransform (x, CDF, pdf, code) The straightforward procedure for setting up the template in a worksheet (Fig. 3) that invokes the user-created function EqvN (of Fig. 4) is as follows: 1. The probability distribution types are entered in a vertical column. In Fig. 3, the distribution names are entered in cells B5 to B7. (The handy data validation feature, available since Excel 97, can also be used to select from a drop-down list of user pre-specified distribution types.) 2. For each distribution type, the corresponding values of the distribution’s parameters are entered in a row. In Fig. 3, the parametric values are entered in cells D5 to G7. The meanings of Para1 to Para4 are explained in cells D11–G20. 3. The mean values of the random variables are entered as initial values in column I under the heading ‘‘x ”. The mean values for each distribution type can be determined from the distribution parameters, as explained in cells I11 to I20. In Fig. 3, the initial values of x (prior to Solver search) were 40, 50 and 1000. 4. The formula ‘‘=EqvN(DistributionName, paralist, x, 1)’’ is entered in cell J5, where ‘‘DistributionName’’ is input by selecting cell B5, paralist is input by selecting the entire range D5:G5 (including the blank F5 and G5 cells), x is input by selecting cell I5, and the argument ‘‘1’’ informs function EqvN (of Fig. 4) that the computed Rackwitz-Fiessler equivalent normal mean is to be returned to cell J5. The formula thus entered in cell J5 is autofilled down the column to the last random variable, in this case cell J7. The formulas of EqvN for up to 50 random variables can thus be promptly set up within seconds. 5. Similarly, the formula ‘‘=EqvN(DistributionName, paralist, x, 2)’’ is entered in cell K5. The argument ‘‘2’’ signifies to function EqvN of Fig. 4 that the equivalent normal standard deviation N is to be computed and its value returned to cell K5. This formula is likewise autofilled downwards to the last random variable. 6. A correlation matrix is created. In Fig. 3, this is shown in the range M5:O7. 7. Enter and autofill formulas for the nx column:=(x*mN)/ N. 8. The array formula for the index, Eq. (1b), is entered in cell S5: ‘‘=sqrt(mmult(transpose(nx),mmult(minverse(crmat),nx)))’’ (followed by pressing ‘‘Enter’’ while holding down the ‘‘Ctrl’’ and ‘‘Shift’’ keys.)
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
77
9. The formula (or user-defined VBA functions) for the limit state surface g(x) is created. In this case the formula is simply ‘‘=Y*ZM’’, where the values of Y, Z, and M are input by selecting from the x values in column I. 10. The constrained optimization program Solver that resides in Microsoft Excel is invoked, to set (cell S5) equal to minimum, by changing (automatically) the x* values, subject to the constraint that g(x)=0. Numerous additional constraints (including formula constraints) can be specified when necessary. It is desirable to activate the Solver option ‘‘Use Automatic Scaling’’ when the changing cells and target cell have large differences in magnitude. The other Solver options may be left at their default settings. For complicated problems (e.g. Figs. 6–9 in later sections), the main developmental effort (restricting attention to computations only) is Step 9. The other steps are straightforward. In particular, the two equations EqvN(. . ., 1) and EqvN(. . ., 2) are set-up promptly in steps 4 and 5 and autofilled down the columns, without the need to remember distribution-specific equations, since these are now decided at the program level (Fig. 4). (The two-parameter Weibull distribution is given as Weibull (,) in Excel 2002, @RISK 2001, and Vose [19]; as Weibull ( ,) in Hahn and Shapiro [18]; and as Weibull ( ,) in Evans et al. [20]. Comparisons of the equations show that the symbol in Evans et al. [20] is the shape parameter and hence equal to the (not ) in Excel 2002, @RISK 2001, and Vose [19]. The in Evans et al. [20] is the scale parameter, the same as the (not ) in Hahn and Shapiro [18]. To make explicit the meanings of the parameters, the equations for Gamma and Weibull distributions are displayed in Fig. 3. The symbol l is used in place of , since the latter denotes reliability index. For beta distribution, the notations can be ascertained by comparing the formulas for the mean).
5. Testing robustness for various nonnormal distributions and comparisons with monte-carlo simulations Fig. 3 used only the lognormal and Type 1 Extreme distribution options. In this section the other distribution types are combined arbitrarily to test robustness and accuracy, by comparing with Monte Carlo simulation using @RISK. For ease of illustrations, three correlated nonnormals Z1, Z2, and Z3 are considered, with performance function g(z)=Z1*Z2Z3, the same as in Figs. 2 and 3. However, no particular physical meanings need be associated with these variables, so that one is free to assume various distributions, and to assign positive or negative correlations. The set-up of Fig. 3 is convenient for experimenting with different combinations of distributions: one need only input (or select from a drop-down list) new distribution names in cells B5 to B7, modify the parameters in cells D5 to G7, initialize x* values to mean values (according to the relationships shown in cells I11 to I20), input changes (if any) to the correlation matrix, and invoke Solver. If more random variables are to be treated (as in Figs. 7 and 8 later), it is only necessary to add the corresponding number of rows beneath row number 7, and expand the correlation matrix accordingly. If Weibull distribution is used, with parameters and l, its mean value (for initial input in x* column) can be computed using the formula ‘‘=l/*EXP(GAM*EXP(GAMMALN(1/))’’.
78
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Of the seven cases shown in Fig. 5, the upper three cases share the same correlation matrix, but differ in their probability distributions. The lower four cases share the same correlation matrix which differs from the upper three cases. Because the means, standard deviations and distributions are not entirely identical, the lower four cases are not comparable, hence arbitrary notations (PQR, STU, GHJ, and VWZ) have been used to denote the variables Z1, Z2, and Z3 in the lower four cases. It is assumed that Z1 and Z2 are positively correlated with coefficient 0.5, Z1 and Z3 negatively correlated with coefficient 0.5, and Z2 and Z3 positively correlated with coefficient 0.4, as shown. The Monte Carlo failure probabilities shown in the last column of Fig. 5 were obtained using @RISK version 4.0. The correlation of input distributions in @RISK is based on the Spearman rank order correlations instead of the Pearson product moment correlations. The advantages and disadvantages of rank order correlation and its use in generating random numbers for correlated nonnormal distributions were discussed in Altman [21], Vose [19], and @RISK manual, for example. For the cases in hand, the failure probabilities based on reliability indices agree remarkably well with those from Monte Carlo simulations, as shown in Fig. 5.
Fig. 5. Comparisons of failure probabilities inferred from reliability indices and from Monte Carlo simulations, for various combinations of distributions.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
79
Case 3 in Fig. 5, in which all the three random variables have Type 1 Extreme distributions, requires additional constraints be specified during Solver optimization, to define the permissible bounds of the random variables. The reason for this is the nature of negative the double expo nential of the Type 1 Extreme distribution function: Fðx Þ ¼ exp exp½ðx aÞ=b ; which leads to rapid variation of F(x*). The design point and index were obtained successfully by restricting Solver’s search to the following range (obtained by re-arrangement of the distribution function F(x*)): a bln½lnðCDFmin Þ 4 x 4 a bln½lnðCDFmax Þ
ð5Þ
with CDFmin=1010 and CDFmax=(11010), and where a and b are Type 1 Extreme distribution parameters, given by a ¼ mean 0:450 StDev and b ¼ 0:7797 StDev: For this case various random initial trial values of x* (within the bounds specified by Eq. (5)) were also explored, and all found to yield the same design point (35.31, 48.19, 1702) and the same index of 2.718. The difference between the probabilities of failure based on reliability index (0.329%) and that based on Monte Carlo simulations (0.303% ) may therefore be attributed to the approximate nature of equivalent normal transformation and the possible deviation of the limit state surface from the assumed hyperplane. The good agreements shown in Fig. 5 may not always be expected, for there are cases (Der Kiureghian et al. [10], Wu and Wirsching [12], Zhao and Ono [14], for example) where the agreements were less remarkable than those between more refined methods (e.g., SORM) and Monte Carlo simulations. It is therefore of interest to investigate the extent of agreement when more complicated performance functions and greater number of correlated nonnormal distributions are involved. This is done after the next two sections.
6. Comparing expanding ellipsoid, Hasofer–Lind method, and FORM The cell-object-oriented constrained optimization approach in the previous two sections achieves the same reliability index values as the widely documented Hasofer–Lind method for uncorrelated and correlated normal variates and the first order reliability method (FORM) for correlated nonnormals. The differences between the present approach and the Hasofer–Lind method and FORM are in the alternative interpretation of the Hasofer–Lind index and in the computational approach. When the random variables are normally distributed, the present approach seeks the smallest multinormal dispersion ellipsoid that is tangent to the limit state surface (Fig. 1) in the original space of the random variables. The vector [(ximi)/ i] in Eq. (1b) is kept intact, thus retaining the original coordinate system. This is in contrast to the Hasofer–Lind method (e.g., [8], pp. 198–202, or [4], pp. 359–361) which (for uncorrelated normals) defines the reduced variable xi0 =(ximi)/ i, and uses it to transform the original limit state g(x)=0 to the reduced limit state g(x0 )=0 in the transformed or reduced coordinate system. The index is then interpreted as the minimum distance from the origin of the axes in the reduced coordinate system to the reduced limit state surface. If the normal variables are correlated, the Hasofer–Lind method also requires a rotation of frame of reference to convert the correlated normals into uncorrelated variables via orthogonal transfor-
80
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
mation ([8], pp. 231–234, [4], pp. 354–368). In the process the correlation matrix is reduced to a diagonal matrix using eigenvalue and eigenvector techniques. This is in contrast with the present approach which leaves the correlation matrix intact. The present approach seeks the expanding dispersion ellipsoid which touches the limit state surface at the most probable point of failure (MPP) in the original coordinate system. As discussed earlier in connection with Eq. (1a) and (1b), that the tangent point is the most probable failure point is evident from the established form of the probability density function of the multinormal distribution: 1 1 T 1 exp ðx mÞ C ðx mÞ ð6aÞ f ðxÞ ¼ n 2 ð2 Þ2 jCj0:5 1 1 2 ¼ exp n 2 ð2 Þ2 jCj0:5
ð6bÞ
where is defined by Eq. (1a) or (1b), without the ‘‘min’’. Since minimizing (or 2) means maximizing f(x), minimizing subject to x being on the limit state surface means the same as finding the most probable failure point on the limit state surface. This intuitive and visual understanding of the design point is an alternative to the more mathematical approach in [3]. For correlated nonnormals, the present approach seeks the equivalent normal dispersion hyperellipsoid tangent to the limit state surface in the original coordinate system: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T N x i mN 1 xi mi i ð7Þ ½R ¼ min x2F iN iN N where mN i and i are the mean and standard deviation of the equivalent normal distribution based on the Rackwitz–Fiessler two-parameter equivalent normal transformation [Eqs. (3) and (4)], a transformation used also in FORM. However, unlike FORM, the equivalent ellipsoid approach does not need to diagonalize the correlation matrix or the covariance matrix, nor does it need to repeatedly reformulate the reduced limit state function in the transformed space. The integrated, expedient and transparent nature of the present approach can perhaps be appreciated by noting that the operational procedure (from the user-point of view) for correlated nonnormals differs from that of the correlated normals only in the distribution names entered in the first column of Fig. 5, and the corresponding parameters (Para1–Para4) characterizing the nonnormal distributions. It may be mentioned that the expanding ellipsoid approach [based on Eq. (1b)] and the expanding equivalent ellipsoid approach [based on Eq. (7)] in the original space of the variables are not vulnerable to numerical ill-conditioning even when the random variables show several orders of magnitude N difference in their magnitude, by virtue of the form (ximi)/ i in Eq. (1b) and (ximN i )/ i in Eq. (7). The significant operational advancement and improved clarity of the present ellipsoidal approach over the Low and Tang [16] ellipsoidal approach can be appreciated by comparing the reliability analyses of Fig. 5 (which invokes the codes of Fig. 4) with the reliability analysis of Fig. 2 (which implements the Rackwitz–Fiessler transformation directly on the worksheet). Had the seven cases of Fig. 5 been analyzed by the approach of Fig. 2, a new template would have to be devised for each different nonnormal distribution to be treated. Even increasing the number of
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
81
variates would require tedious multi-cell equation changes. It may also be pointed out that had one followed the widely documented classical FORM procedure of seeking the shortest distance in transformed space, with eigenvalues and eigenvectors for orthogonal transformation, and iteratively reformulating the reduced limit state function for each trial equivalent normal transformation, the procedure and codes (in the same EXCEL platform) would have been far more complicated and less transparent than the ellipsoid approach presented herein, even though the final computed index will be the same by both approaches.
7. On the need for positive definite correlation matrix Although the correlation coefficient between two random variables has a range 14 ij41, one is not totally free in assigning any values within this range for the correlation matrix. This was explained in Ditlevsen p85 [2] where a 33 correlation matrix with 12= 13= 23=1 implies a contradiction and is hence inconsistent. A full-rank (i.e., no redundant variables) correlation matrix has to be positive definite. One may note that Hasofer–Lind index is defined as the squareroot of a quadratic form [Eq. (1a) or (1b)], and unless the quadratic form and its covariance matrix and correlation matrix are positive-definite, one is likely to encounter the square root of a negative entity, which is meaningless whether in the dispersion ellipsoid perspective or in the perspective of minimum distance in the reduced variate space. In Fig. 5, the correlation matrix of the last four cases becomes singular when 23= 32=0.5. The quadratic form evaluates to a negative value when 23= 32 > 0.5. This means that given
12=0.5 and 13=0.5, the correlation coefficient 23 has to be smaller than 0.5. On the other hand, if 12= 13=0.8, the correlation coefficient 23 has to be greater than 0.3 in order for R to be positive definite. In the latter case it can be appreciated intuitively that when both the random variables Z2 and Z3 are highly and positively correlated with Z1 (via 12= 13=0.8), then Z2 and Z3 should also be positively correlated with each other to a certain extent ( >0.3). Mathematically, a symmetric matrix R is positive definite if and only if all the eigenvalues of R are positive. For present purposes, to detect inconsistent correlation matrix, one can create the p equation QForm=uTR1u side by side with = (QForm). During Solver’s searches, if the cell Qform shows negative value, it means that the quadratic form is not positive definite and hence the correlation matrix is inconsistent. The Monte Carlo simulation program @RISK will also display a warning of ‘‘Invalid correlation matrix’’ and offer to correct the matrix if simulation is run involving an inconsistent correlation matrix. That a full-rank correlation or covariance matrix should be positive definite is an established theoretical requirement [29–32], although not widely known. This positive-definiteness requirement is therefore not a limitation of the paper’s proposed procedure. Further reasons for positive definiteness follow. The joint p.d.f. of an n-dimensional multivariate normal distribution is given by Eq. (6a). The p.d.f. has a peak at mean-value point and decreasing ellipsoidal probability contours surrounding m only if the quadratic form and hence the covariance matrix and its inverse are positive definite. (By lemmas of matrix theory: C positive definite implies C1 is positive definite, and also means the correlation matrix R and its inverse are positive definite.)
82
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Without assuming any particular probability distributions, it is documented in books on multivariate statistical analysis (e.g. [29–32]) that a covariance matrix (and hence a correlation matrix) cannot be allowed to be indefinite, because this would imply that a certain weighted sum of its variables has a negative variance, which is inadmissible, since all real-valued variables must have nonnegative variances. (The classical proof can be viewed as reductio ad absurdum, i.e, a method of proving the falsity of a premise by showing that the logical consequence is absurd.) An internet search using the Google search engine with ‘‘positive definite correlation matrix’’ as search terms also revealed thousands of relevant web pages, showing that many statistical computer programs (in addition to the @RISK simulation software) will check for positive definiteness of correlation matrix and offer to generate the closest valid matrix to replace the entered invalid one. For example, one such web site titled ‘‘Not Positive Definite Matrices–Causes and Cures’’ at http://www.gsu.edu/mkteer/npdmatri.html refers to the work of Wothke [32].
8. Reliability analysis of finite beam on Winkler medium Beams on elastic foundations can be analyzed using various analytical and numerical methods (Hetenyi [22], Selvadurai [23], Boresi et al. [24]). In this section reliability analysis will be performed for a beam on Winkler medium, to test robustness and accuracy of the present approach when relatively complicated nonlinear performance functions are involved. The beam has a length (l) of 15 m, a width (b) of 1 m. Six random variables are considered: the moment of inertia Ib (in m4) and Young’s modulus Eb (MN/m2) of the beam, the subgrade modulus ks (MN/m3), and the magnitudes (in MN) of three concentrated loads W1, W2, and W3, acting at 1.5, 7.5, and 13.5 m respectively from the left end of the beam. Failure is considered to have occurred if the bending moment reaches 1 MNm anywhere along the beam. The analytical solutions summarized in Young and Budynas [25] and shown in Fig. 6 are incorporated in a user-created function BMoment (. . .) whose input parameters include the load vector (W1, W2, W3), position vector (1.5, 7.5, 13.5 m), x coordinate and lamda (equal to (bks/4EbIb)0.25). The loads W1, W2, and W3 are assumed to have a coefficient of variation equal to 0.15. The assumed distribution types, correlation and other input are as shown in Fig. 7. The computed b index is 2.383. The reliability-based probability of failure (0.86%) compared very well with simulation (Pf=0.85%) using @RISK with Latin Hypercube sampling. Reliability analysis using Solver in Excel XP obtained the b index in 10 seconds (on a Pentium4 PC) after 9 trial solutions. In using Solver, the ‘‘automatic scaling’’ option has been activated.
9. Reliability analysis of complex strut involving implicit performance function The deterministic analysis of a strut with complex supports was discussed in Coates et al. [26]. The member is initially straight (Fig. 8) with a pin support at 2, an elastic support at 1 (rotational stiffness l1) which provides a restoring moment M1=l11, and a support at 3 (stiffness k3) which provides a reaction force (k3v3) proportional to the vertical displacement v3. The problem is to determine the smallest value of axial force P which will cause the strut to become elastically unstable. Coates et al. [26] showed that the problem reduces to that of finding a value of P, referred
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
83
to as Pcrt, which would make the determinant of a 7 by 7 matrix (functions of P, L, a, E, I, l1 and k3) vanish. It was also noted that a numerical procedure would be necessary to determine the determinant of the matrix by sequential variation of the axial load P until a zero value is reached. For the reliability analysis here, the critical load Pcrit can be obtained in a spreadsheet cell (Fig. 8) that calls the user-defined VBA function Pcritical (L, a, E, I, l1, k3) shown in Fig. 9. The performance function (PerFunc in Fig. 8) is ‘‘=Pcrt.P’’, where P refers to the first value of the x* column. Hence the single cell object ‘‘PerFunc’’ contains the implicit and iterative user-created routine of Fig. 9. For the distributions and input in Fig. 8, the initial (mean) values of the x* column were (700, 1000, 500, 200000, 200, 500, 10). Solver in Excel XP (2002) was then invoked to minimize the
Fig. 6. Bending moment equations for finite-length beams on Winkler elastic medium. (Ref. Young and Budynas, 2002).
84
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
cell, by changing the x* values, subject to the constraint that PerFunc=0. The Solver option ‘‘Use Automatic Scaling’’ was activated. Solver obtains the solution in 35 s (in a Pentium4 PC) after 14 trial solutions. The probability of failure (0.40%) based on the index of 2.652 compares reasonably well with the 0.36% failure probability from simulation with 250,000 trials. Note that in this case where there is five orders of magnitude difference between random variables E and k3, use of ‘‘Automatic Scaling’’ is essential for efficiency. The design point (x*) and the nx values in Fig. 8 indicate that the parameters l1 and k3 hardly deviate from their mean values of 500 and 10 respectively. The implication is that the performance function and hence Pcritical is not sensitive to l1 and k3 at their current mean values and
Fig. 7. Reliability analysis of asymmetrically loaded beam on Winkler medium, and comparison with Monte Carlo simulation.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
85
levels of uncertainty. These rather unexpected findings were further investigated in four analyses in which the initial values of l1 and k3 were set at two standard deviations above and below their respective mean values, while the other five parameters have initial values at their respective mean values. The indices (vary between 2.6513 and 2.6524) and the design point x* values obtained after constrained optimization were virtually identical to those of Fig. 8, thereby providing further support to the finding that Pcritical is insensitive to l1 and k3. Coates et al. [26] also investigated sensitivities, deterministically, by plotting curves of Pcritical versus k3, for values of l1 ranging from 0 to 1020 Nmm/ rad, and found that very large changes in l1 had only very small influence on Pcritical. In comparison, for the mean values and standard deviations of Fig. 8, reliability analysis reveals the insensitivities of Pcritical to l1 and k3 with much less effort than deterministic contour plotting, besides providing additional information, including a good estimate of the probability of strut failure.
Fig. 8. Reliability analysis of a strut with complex supports. Performance function=(Pcrt.P), where Pcrt. is programcalculated by the code of Fig. 9.
86
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
10. Limitations and flexibilities of object-oriented constrained optimization in spreadsheet There is no closed form equation for the inverse of the standard normal cumulative distribution. Numerical procedures such as that incorporated in Excel’s NormSInv is necessary. In Excel 2000 and earlier versions, NormSInv(CDF) returns correct values only for CDF between about 0.00000035 and 0.99999965, i.e., within about 5 standard deviations from the mean. There are undesirable and idiosyncratic sharp kinks at both ends of the validity range. In EXCEL XP (2002), the NormSInv (CDF) function has been refined, and returns correct values for a much broader band of CDF: 1016
11. Summary and conclusions An expanding ellipsoid perspective (Low and Tang [16]) provides an alternative and intuitive interpretation of the Hasofer–Lind index. The perspective led to a practical spreadsheet optimization procedure for performing reliability analysis. This paper extends the 1997 approach by providing a flexible user interface (Figs. 3 and 4) in which numerous different nonnormal distributions can be selected with ease for reliability analysis, and by testing the robustness and
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
87
accuracy of the approach for two relatively complicated problems: a beam on elastic foundations, and a strut with complex supports. The failure probabilities based on reliability indices are in good agreements with failure probabilities from Monte Carlo simulations (Figs. 5, 7, and 8). The main features of the procedure presented herein for reliability analysis involving correlated nonnormals are: The design point and reliability index are obtained using object-oriented constrained optimization of equivalent hyperellipsoids in the original space of the random variables.
Fig. 9. An example of program-calculated performance function that is incorporated in a single cell object in Fig. 8.
88
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The approach obtains the same result as FORM, but is operationally more direct and transparent. There is no need to diagonalize the correlation matrix. The platform is ubiquitous. The Rackwitz-Fiessler equivalent normal transformations are executed in a short program code (Fig. 4). This renders possible a consistent and uncluttered user-interface (Figs. 3, 7, and 8). To specify different nonnormal distributions, it is only necessary to name the distributions and provide its parameters. Performance functions can be implicit (e.g. containing VBA codes) and involving various correlated nonnormal variates, as in Figs. 8 and 9, for example. In general, for reliability analysis involving correlated nonnormal distributions, the accuracy of the present approach may depend on the nature of the performance functions, the number of random variables, the validity range of the NormSInv function (which is used for equivalent normal transformation), and the robustness of the nonlinear constrained optimization program. For the cases of correlated nonnormals analyzed in this study, Excel XP (with its refined NormSInv function) and its built-in standard Solver are found to be efficient and accurate (as far as FORM is concerned). The main advantage of the object-oriented constrained optimization approach for reliability analysis is in the modular nature of its constituents: The deterministic model that underlies the performance function g(x) is set-up (or programmed in VBA) separately and appears as one or more spreadsheet cell objects. The nonlinear constrained optimization program Solver is then invoked to minimize the cell object (itself a conglomerate of program routines), subject to constraints on the performance function cell object and other cell objects. The platform can be other than Microsoft Excel. There are also other choices of potentially more robust spreadsheetbased optimization programs, though less ubiquitous and at a cost.
References [1] Hasofer AM, Lind NC. Exact and invariant second-moment code format. J of Engrg Mech ASCE, New York 1974;100(1):111–21. [2] Ditlevsen O. Uncertainty modeling: with applications to multidimensional civil engineering systems. New York: McGraw-Hill; 1981. [3] Shinozuka M. Basic analysis of structural safety. Journal of Structural Engineering, ASCE 1983;109(3):721–40. [4] Ang HS, Tang WH. Probability concepts in engineering planning and design, v2: decision, risk, and reliability. New York: John Wiley; 1984. [5] Madsen HO, Krenk S, Lind NC. Methods of structural safety. NJ: Prentice-Hall: Englewood Cliffs; 1986. [6] Tichy M. Applied methods of structural reliability. Dordrecht, Boston: Kluwer Academic; 1993. [7] Melchers RE. Structural reliability analysis and prediction. 2nd edition. John Wiley; 1999. [8] Haldar A, Mahadevan S. Probability, reliability and statistical methods in engineering design. New York: John Wiley; 1999. [9] Nowak AS, Collins KR. Reliability of structures. McGraw-Hill; 2000. [10] Der Kiureghian A, Lin HZ, Hwang SJ. Second-order reliability approximations. Journal of Engineering Mechanics, ASCE 1987;113(8):1208–25. [11] Chen X, Lind NC. Fast probability integration by three-parameter normal tail approximations. Structural Safety 1983;1:269–76.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
89
[12] Wu YT, Wirsching PH. New algorithm for structural reliability estimation. Journal of Engineering Mechanics, ASCE 1987;113(9):1319–35. [13] Zhao, YG, Tetsuro Ono. A general procedure for first/second-order reliability method (FORM/SORM). Structural Safety 1999;21:95–112. [14] Zhao YG, Tetsuro Ono. Moment methods for structural reliability. Structural Safety 2001;23:47–75. [15] Rackwitz R. Reliability analysis-a review and some perspectives. Structural Safety 2001;23:365–95. [16] Low BK, Tang WH. Efficient reliability evaluation using spreadsheet. Journal of Engineering Mechanics, ASCE 1997;123(7):749–52. [17] Rackwitz R, Fiessler B. Structural reliability under combined random load sequences. Computers & Structures 1978;9:484–94. [18] Hahn GJ, Shapiro SS. Statistical models in engineering. John Wiley & Sons; 1967. [19] Vose D. Risk analysis: a quantitative guide. 2nd edition. John Wiley & Sons; 2000. [20] Evans M, Hastings N, Peacock B. Statistical distributions. 3rd edition. New York: Wiley; 2000. [21] Altman DG. Practical statistics for medical research. London: Chapman and Hall; 1991. [22] Hetenyi M. Beams on elastic foundations. Ann Arbor, Michigan: Univ of Mich Press; 1974. [23] Selvadurai APS. Combined foundations. In Henry FDC, editor. The Design and Construction of Engineering Foundations (2nd ed.). London: Chapman & Hall; 1986, chapter 5. [24] Boresi AP, Schmidt RJ, Sidebottom OM. Advanced mechanics of materials. 5th edition. New York: Wiley; 1993. [25] Young WC, Budynas RG. Roark’s formulas of stress and strain. 7th edition. McGraw-Hill; 2002. [26] Coates RC, Coutie MG, Kong FK. Structural analysis. 3rd edition. London: Chapman & Hall; 1994. [27] Low BK. Probabilistic slope analysis involving generalized slip surface in stochastic soil medium. In Proceedings of the 14th South East Asian Geotechnical Conference, 10-14 December 2001, Hong Kong, A.A.Balkema Publishers; 2001. p.825–830. [28] Low BK, Teh CI, Tang WH. Stochastic nonlinear p-y analysis of laterally loaded piles. Proceedings of the Eight International Conference on Structural Safety and Reliability, ICOSSAR ‘01, Newport Beach, California, 17–22 June 2001. 8 pages. [29] Chatfield C, Collins AJ. Introduction to multivariate analysis. London: Chapman and Hall; 1989 (reprint). [30] Graybill FA. Matrices with applications in statistics. 2nd ed. Belmont, CA: Wadsworth; 1983. [31] Searle SR. Matrix algebra useful for statistics. New York: John Wiley; 1982. [32] Wothke W. Nonpositive definite matrices in structural modelling. In: Bollen KA, Long JS, editors. Testing structural equation models. Newbury Park, CA: Sage Publications, Inc; 1993. p. 257–93.
Reliability analysis using object-oriented constrained optimization Bak K. Lowa,*, Wilson H. Tangb a
School of Civil and Environmental Engineering, Nanyang Technological University, Block N1, #01b-40, 50 Nanyang Avenue, Singapore 639798, Singapore b Department of Civil and Structural Engineering, Hong Kong University of Science & Technology, Clear Water Bay, Kowloon, Hong Kong
Abstract A practical procedure is presented for reliability analysis involving correlated nonnormals. Equivalent normal transformation is performed as in the first order reliability method; the correlation matrix is used as it is without orthogonal transformation. Performance functions appear as simple spreadsheet cell objects but may embody substantial user-created program codes. The multidimensional equivalent dispersion ellipsoid is likewise a single cell object that contains various spreadsheet routines for equivalent normal transformation and matrix operations. The procedure can be viewed as the constrained optimization of the equivalent hyperellipsoid (yielding the reliability index) in the original space, and is automatic and efficient in the ubiquitous spreadsheet platform. Different probability distributions can be selected at ease and conveniently transformed to equivalent normals using a short program code given in the paper. Illustrations are presented using a three-variate simple performance function and various combinations of ten common distributions. This is followed by two relatively complicated examples, namely an asymmetrically loaded beam on Winkler medium, and a complex strut with implicit performance function. The probabilities of failure inferred from reliability indices are compared with Monte Carlo simulations. The simplicity, transparency and flexibilities of the object-oriented constrained optimization approach are demonstrated. # 2003 Elsevier Ltd. All rights reserved.
1. Introduction The Hasofer–Lind [1] index and the first order reliability method (FORM) have been welldocumented in Ditlevsen [2], Shinozuka [3], Ang and Tang [4], Madsen et al. [5], and Tichy [6], * Corresponding author. Tel.: +65-6790-5270; fax: +65-6791-0676. E-mail address: [email protected] (B.K. Low). 0167-4730/03/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0167-4730(03)00023-7
70
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
for example. The continuing interest in FORM is also evident in Melchers [7], Haldar and Mahadevan [8], Nowak and Collins [9]. The potential inadequacies of FORM in some special cases have been recognized, and more refined alternatives proposed, in Der Kiureghian et al. [10], Chen and Lind [11], Wu and Wirsching [12], and Zhao and Ono [13,14]. On the other hand, Rackwitz [15] noted the usefulness and accuracy of FORM in most applications; some counterexamples were summarized, their rather extreme nature and artificiality noted, but their healthy role also admitted. With respect to the Hasofer–Lind index and FORM, a practical procedure was presented in Low and Tang [16], based on the perspective of an expanding equivalent ellipsoid tangent to the limit state surface in the original space of the random variables. The Rackwitz–Fiessler [17] equivalent normal transformation was used, but the concepts of coordinate transformation and frame-of-reference rotation were not required. Correlation was accounted for by setting up the quadratic form directly. Iterative searching and partial derivatives were automatic using constrained optimization in the ubiquitous spreadsheet platform. The present paper extends Low and Tang [16], by testing robustness and accuracy for various nonnormal distributions and more complicated performance functions, and by providing enhanced operational advantage and versatility. A simple 3-variate case in the next section elucidates the principles of object-oriented constrained optimization, and prepares the background for the more efficient and flexible strategy that follows. Different combinations of ten common distributions are used. Two relatively complicated cases are also analyzed, namely an asymmetrically loaded beam on Winkler elastic medium, and a complex strut the performance function of which requires iterative program calculations. The results obtained by the spreadsheet object-oriented constrained optimization approach are compared with Monte Carlo simulations.
2. Constrained optimization of dispersion ellipsoid in original space The matrix formulation of the Hasofer–Lind index is: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ min ðx mÞT C1 ðx mÞ x2F
ð1aÞ
or, equivalently: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi xi mi T 1 xi mi ¼ min ½R x2F i i
ð1bÞ
where x is a vector representing the set of random variables, m the mean values, C the covariance matrix, R the correlation matrix, and F the failure domain. The index is, in classical explanation, the minimum distance from the mean value point to the limit state surface, in units of directional standard deviations. An alternative interpretation was presented in Low and Tang [16], where the perspective of an expanding ellipsoid (Fig. 1) led to a simple method of computing the Hasofer– Lind index in the original space of the variables.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
71
Fig. 1. Design point and equivalent normal dispersion ellipsoids illustrated in the plane.
One may note that the quadratic form in Eq. (1a), ðx mÞT C1 ðx mÞ; appears also as a negative exponent of the multivariate normal probability density function. Hence, to minimize (or 2 in the multivariate normal density function) is to maximize the value of the multivariate normal probability density function, and to find the smallest ellipsoid tangent to the limit state surface is equivalent to finding the most probable failure point. This perspective is consistent with Shinozuka [3] that ‘‘the design point x* is the point of maximum likelihood if x is Gaussian, whether or not its components are uncorrelated.’’ When the random variables are nonnormals and correlated, the approach to be presented below uses the Rackwitz–Fiessler [17] equivalent normal transformation, but without the need to diagonalize the correlation matrix (which is equivalent to rotating the frame of reference). Also, the iterative computations of the equivalent normal mean (mN) and equivalent normal standard deviation ( N) for each trial design point are automatic during the constrained optimization search.
3. Three-variate correlated nonnormals: a simple example The case in Fig. 4 of Low and Tang [16] for correlated nonnormals is shown in Fig. 2, with a more succinct set-up. The random variables Y and Z are correlated lognormals, and M is Type 1
72
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Fig. 2. Reliability analysis of correlated nonnormals using constrained optimization. Random variables Y and Z are correlated lognormals; M is Type 1 Extreme (Gumbel).
Extreme (Gumbel). The performance function is g(x)=Y*ZM. The original mean () and standard deviations () of Y, Z and M are (40, 5), (50, 2.5), and (1000, 200), respectively. The square root of the quadratic form in Eq. (1b) is computed directly using array formula: ‘‘ ¼ sqrtðmmultðtransposeðnxÞ; mmultðminverseðcrmatÞ; nxÞÞÞ”
ð2Þ
(Array formulas are created by pressing ‘‘Enter’’ while holding down ‘‘Ctrl’’+‘‘Shift’’). In Eq. (2), ‘‘mmult’’, ‘‘transpose’’ and ‘‘minverse’’ are Excel’s built-in functions, each being a container of VBA codes for matrix operations. The nx vector in Eq. (2) contains functions (via mN and N) of Excel routines ‘‘normsinv’’ (cell D20), which returns the inverse of the standard normal cumulative distribution, and ‘‘normdist’’ (cell E7). These two functions are used to transform the Type1 Extreme variate M into equivalent normal distribution, based on Rackwitz– Fiessler [17] two-parameter equivalent normal transformation: F1 ½F ðxÞ N ð3Þ Equivalent normal standard deviation: ¼ f ðxÞ Equivalent normal mean :
mN ¼ x N F1 ½F ðxÞ
ð4Þ
where x is the original nonnormal variate, F1 ½: is the inverse of the cumulative probability (CDF) of a standard normal distribution, F(x) is the original nonnormal CDF evaluated at x,
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
73
{.} is the probability density function (pdf) of the standard normal distribution, and f(x) is the original nonnormal probability density ordinates at x. Therefore, the single spreadsheet object for (cell N7) is a container of several Excel function objects and substantial program codes. For lognormals Y and Z, closed form equivalent normal transformation is available and has been used, as shown by the equations for mN and N of Y and Z in Fig. 2. For Type 1 Extreme variable M, Excel’s built-in routines normsinv (for 1) and normdist (for {.}) have been used to automatically update mN and N (cells D7 and E7) by Eqs. (3) and (4) during the iterative constrained optimizational search for the design point on the limit state surface. The design point (x* values) was obtained by calling Excel’s built-in optimization program Solver, to minimize the cell, by changing the x* values, subject to the constraint that the performance function g(x)=0. Prior to Solver search, the x* values had values equal to the mean values () of the original random variables. Iterative numerical derivatives and directional search for the design point x* were automatic. The value of the reliability index in Fig. 2 implies a probability of failure approximately equal to 0.00385. This compares favorably with a failure rate of 0.00382 from Monte Carlo simulation using the commercial software @RISK 2001 (http://www.palisade.com). The sample size was 200,000. The @RISK distribution functions used were RiskLogNorm(40,5) for variable Y, RiskLogNorm(50, 2.5) for variable Z, and RiskExtvalue(u, 1/) for M, i.e., RiskExtvalue(910, 155.9). The correlation matrix was modeled. Fig. 2 is concise, instructive and relatively transparent for computing the Hasofer–Lind index involving correlated nonnormals. It is ideal for pedagogic purposes. Nevertheless, a clearer and more convenient user interface is desirable so that various nonnormal distributions can be selected at ease without the need to modify the template for different probability distributions. This is developed in the next section. For reliability index computations in the following sections, Microsoft Excel XP (2002) and its built-in optimization program have been used throughout. Alternative optimization programs (e.g., http://www.Solver.com and http://www.palisade.com/ html/evolver.html, including genetic and evolutionary algorithms, are available (at a cost) as Excel add-ins. If desired, these spreadsheet add-ins may be used instead of the standard Excel Solver. The procedures and codes in the next section are unaffected by the type of optimization add-ins used. Should an alternative spreadsheet platform evolve in future, it is simple to modify the short code (Fig. 4 below) to conform to the language syntax of the new platform. For the examples below, Excel XP with its built-in standard Solver is adequate, and has the advantage of being ubiquitous.
4. Reliability analysis using object-oriented constrained optimisation The more succinct and versatile alternative to Fig. 2 is shown in Fig. 3. The key feature is that all distribution-specific equations are relegated to a short user-created (or user-defined, in VBA parlance) program code shown in Fig. 4, which is called by cells beneath the headings ‘‘mN’’ and ‘‘ N’’ in Fig. 3. The VBA (Visual Basic for Applications) program code in Fig. 4 is created via Tools/Macro/VisualBasicEditor/Insert/Module on the Excel worksheet menu. It is easily transferable to other workbooks.
74
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The main program is EqvN, which returns the equivalent normal mean (if code=1) or the equivalent normal standard deviation (if code=2). Its arguments include ‘‘DistributionName’’ and ‘‘Paralist’’, the latter being a vector of up to 4 parameters, depending on the type of probability distribution. In the fourth line of the program in Fig. 4, the ‘‘Select Case’’ control structure alters the flow of execution to one of several code segments, depending on the input ‘‘DistributionName’’. The functions Trim and UCase are VBA’s built-in functions: ‘‘Trim’’ returns a text with both leading and trailing spaces removed, ‘‘UCase’’ then renders the text uppercase throughout. Function EqvTransform in Fig. 4, called by six distribution types (from EXTVALUE1 to BETADIST) in the main function EqvN, implements the Rackwitz-Fiessler equivalent normal transformation according to Eqs. (3) and (4), based on the cumulative probability distribution (CDF) and probability density function (pdf) of the original nonnormals at x. Function betapdf in Fig. 4, called by BETADIST in the main program EqvN, returns the value of the density function of the beta distribution. The three functions in Fig. 4 make liberal calls to Excel’s built-in functions, using the syntax ‘‘Application.GammaDist(. . .)’’, Application. Weibull(. . .)’’, ‘‘Application.NormDist(.)’’, ‘‘Application.NormSInv(.)’’, where ‘‘Application’’ stands for Excel. The use of these Excel ‘‘objects’’ (each a container of VBA program codes) results in much simplicity, clarity, and brevity of the program codes and structure. The Uniform (or rectangular) distribution can be obtained [18] as a special case of BetaDist (1,1,min,max). If desired, it can also be created (only a few statement lines are required), and be one of the ‘‘Cases’’ in Fig. 4.
Fig. 3. Constrained optimization of equivalent normal ellipsoid in the original space of correlated nonnormals. This approach is more user-friendly and versatile than Fig. 2.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
75
Fig. 4. User-defined VBA codes. These functions were used in cells J5:K7 of Fig. 3 for equivalent normal transformation.
76
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The PERT distribution (an option in @RISK; also described in Vose [19]) is also a special case of Beta distribution. It is defined by the same three parameters (min, mode, max) as the Triangular distribution, but has some advantages over the latter. The PERT distribution shown below can be inserted to the code in Fig. 4: Case ‘‘PERTDIST’’: a=para1: Mode=para2: c=para3 If x<=a Then x=a+del If x>=c Then x=c-del mean=(a+4* Mode+c)/6 If Mode=mean Then f=6 Else f=(2* Modeac)/(Mode-mean) a1=(meana)*f/(ca): a2=a1*(cmean)/(meana) CDF=Application.BetaDist(x, a1, a2, a, c): pdf=betapdf(x, a1, a2, a, c) EqvN=EqvTransform (x, CDF, pdf, code) The straightforward procedure for setting up the template in a worksheet (Fig. 3) that invokes the user-created function EqvN (of Fig. 4) is as follows: 1. The probability distribution types are entered in a vertical column. In Fig. 3, the distribution names are entered in cells B5 to B7. (The handy data validation feature, available since Excel 97, can also be used to select from a drop-down list of user pre-specified distribution types.) 2. For each distribution type, the corresponding values of the distribution’s parameters are entered in a row. In Fig. 3, the parametric values are entered in cells D5 to G7. The meanings of Para1 to Para4 are explained in cells D11–G20. 3. The mean values of the random variables are entered as initial values in column I under the heading ‘‘x ”. The mean values for each distribution type can be determined from the distribution parameters, as explained in cells I11 to I20. In Fig. 3, the initial values of x (prior to Solver search) were 40, 50 and 1000. 4. The formula ‘‘=EqvN(DistributionName, paralist, x, 1)’’ is entered in cell J5, where ‘‘DistributionName’’ is input by selecting cell B5, paralist is input by selecting the entire range D5:G5 (including the blank F5 and G5 cells), x is input by selecting cell I5, and the argument ‘‘1’’ informs function EqvN (of Fig. 4) that the computed Rackwitz-Fiessler equivalent normal mean is to be returned to cell J5. The formula thus entered in cell J5 is autofilled down the column to the last random variable, in this case cell J7. The formulas of EqvN for up to 50 random variables can thus be promptly set up within seconds. 5. Similarly, the formula ‘‘=EqvN(DistributionName, paralist, x, 2)’’ is entered in cell K5. The argument ‘‘2’’ signifies to function EqvN of Fig. 4 that the equivalent normal standard deviation N is to be computed and its value returned to cell K5. This formula is likewise autofilled downwards to the last random variable. 6. A correlation matrix is created. In Fig. 3, this is shown in the range M5:O7. 7. Enter and autofill formulas for the nx column:=(x*mN)/ N. 8. The array formula for the index, Eq. (1b), is entered in cell S5: ‘‘=sqrt(mmult(transpose(nx),mmult(minverse(crmat),nx)))’’ (followed by pressing ‘‘Enter’’ while holding down the ‘‘Ctrl’’ and ‘‘Shift’’ keys.)
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
77
9. The formula (or user-defined VBA functions) for the limit state surface g(x) is created. In this case the formula is simply ‘‘=Y*ZM’’, where the values of Y, Z, and M are input by selecting from the x values in column I. 10. The constrained optimization program Solver that resides in Microsoft Excel is invoked, to set (cell S5) equal to minimum, by changing (automatically) the x* values, subject to the constraint that g(x)=0. Numerous additional constraints (including formula constraints) can be specified when necessary. It is desirable to activate the Solver option ‘‘Use Automatic Scaling’’ when the changing cells and target cell have large differences in magnitude. The other Solver options may be left at their default settings. For complicated problems (e.g. Figs. 6–9 in later sections), the main developmental effort (restricting attention to computations only) is Step 9. The other steps are straightforward. In particular, the two equations EqvN(. . ., 1) and EqvN(. . ., 2) are set-up promptly in steps 4 and 5 and autofilled down the columns, without the need to remember distribution-specific equations, since these are now decided at the program level (Fig. 4). (The two-parameter Weibull distribution is given as Weibull (,) in Excel 2002, @RISK 2001, and Vose [19]; as Weibull ( ,) in Hahn and Shapiro [18]; and as Weibull ( ,) in Evans et al. [20]. Comparisons of the equations show that the symbol in Evans et al. [20] is the shape parameter and hence equal to the (not ) in Excel 2002, @RISK 2001, and Vose [19]. The in Evans et al. [20] is the scale parameter, the same as the (not ) in Hahn and Shapiro [18]. To make explicit the meanings of the parameters, the equations for Gamma and Weibull distributions are displayed in Fig. 3. The symbol l is used in place of , since the latter denotes reliability index. For beta distribution, the notations can be ascertained by comparing the formulas for the mean).
5. Testing robustness for various nonnormal distributions and comparisons with monte-carlo simulations Fig. 3 used only the lognormal and Type 1 Extreme distribution options. In this section the other distribution types are combined arbitrarily to test robustness and accuracy, by comparing with Monte Carlo simulation using @RISK. For ease of illustrations, three correlated nonnormals Z1, Z2, and Z3 are considered, with performance function g(z)=Z1*Z2Z3, the same as in Figs. 2 and 3. However, no particular physical meanings need be associated with these variables, so that one is free to assume various distributions, and to assign positive or negative correlations. The set-up of Fig. 3 is convenient for experimenting with different combinations of distributions: one need only input (or select from a drop-down list) new distribution names in cells B5 to B7, modify the parameters in cells D5 to G7, initialize x* values to mean values (according to the relationships shown in cells I11 to I20), input changes (if any) to the correlation matrix, and invoke Solver. If more random variables are to be treated (as in Figs. 7 and 8 later), it is only necessary to add the corresponding number of rows beneath row number 7, and expand the correlation matrix accordingly. If Weibull distribution is used, with parameters and l, its mean value (for initial input in x* column) can be computed using the formula ‘‘=l/*EXP(GAM*EXP(GAMMALN(1/))’’.
78
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Of the seven cases shown in Fig. 5, the upper three cases share the same correlation matrix, but differ in their probability distributions. The lower four cases share the same correlation matrix which differs from the upper three cases. Because the means, standard deviations and distributions are not entirely identical, the lower four cases are not comparable, hence arbitrary notations (PQR, STU, GHJ, and VWZ) have been used to denote the variables Z1, Z2, and Z3 in the lower four cases. It is assumed that Z1 and Z2 are positively correlated with coefficient 0.5, Z1 and Z3 negatively correlated with coefficient 0.5, and Z2 and Z3 positively correlated with coefficient 0.4, as shown. The Monte Carlo failure probabilities shown in the last column of Fig. 5 were obtained using @RISK version 4.0. The correlation of input distributions in @RISK is based on the Spearman rank order correlations instead of the Pearson product moment correlations. The advantages and disadvantages of rank order correlation and its use in generating random numbers for correlated nonnormal distributions were discussed in Altman [21], Vose [19], and @RISK manual, for example. For the cases in hand, the failure probabilities based on reliability indices agree remarkably well with those from Monte Carlo simulations, as shown in Fig. 5.
Fig. 5. Comparisons of failure probabilities inferred from reliability indices and from Monte Carlo simulations, for various combinations of distributions.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
79
Case 3 in Fig. 5, in which all the three random variables have Type 1 Extreme distributions, requires additional constraints be specified during Solver optimization, to define the permissible bounds of the random variables. The reason for this is the nature of negative the double expo nential of the Type 1 Extreme distribution function: Fðx Þ ¼ exp exp½ðx aÞ=b ; which leads to rapid variation of F(x*). The design point and index were obtained successfully by restricting Solver’s search to the following range (obtained by re-arrangement of the distribution function F(x*)): a bln½lnðCDFmin Þ 4 x 4 a bln½lnðCDFmax Þ
ð5Þ
with CDFmin=1010 and CDFmax=(11010), and where a and b are Type 1 Extreme distribution parameters, given by a ¼ mean 0:450 StDev and b ¼ 0:7797 StDev: For this case various random initial trial values of x* (within the bounds specified by Eq. (5)) were also explored, and all found to yield the same design point (35.31, 48.19, 1702) and the same index of 2.718. The difference between the probabilities of failure based on reliability index (0.329%) and that based on Monte Carlo simulations (0.303% ) may therefore be attributed to the approximate nature of equivalent normal transformation and the possible deviation of the limit state surface from the assumed hyperplane. The good agreements shown in Fig. 5 may not always be expected, for there are cases (Der Kiureghian et al. [10], Wu and Wirsching [12], Zhao and Ono [14], for example) where the agreements were less remarkable than those between more refined methods (e.g., SORM) and Monte Carlo simulations. It is therefore of interest to investigate the extent of agreement when more complicated performance functions and greater number of correlated nonnormal distributions are involved. This is done after the next two sections.
6. Comparing expanding ellipsoid, Hasofer–Lind method, and FORM The cell-object-oriented constrained optimization approach in the previous two sections achieves the same reliability index values as the widely documented Hasofer–Lind method for uncorrelated and correlated normal variates and the first order reliability method (FORM) for correlated nonnormals. The differences between the present approach and the Hasofer–Lind method and FORM are in the alternative interpretation of the Hasofer–Lind index and in the computational approach. When the random variables are normally distributed, the present approach seeks the smallest multinormal dispersion ellipsoid that is tangent to the limit state surface (Fig. 1) in the original space of the random variables. The vector [(ximi)/ i] in Eq. (1b) is kept intact, thus retaining the original coordinate system. This is in contrast to the Hasofer–Lind method (e.g., [8], pp. 198–202, or [4], pp. 359–361) which (for uncorrelated normals) defines the reduced variable xi0 =(ximi)/ i, and uses it to transform the original limit state g(x)=0 to the reduced limit state g(x0 )=0 in the transformed or reduced coordinate system. The index is then interpreted as the minimum distance from the origin of the axes in the reduced coordinate system to the reduced limit state surface. If the normal variables are correlated, the Hasofer–Lind method also requires a rotation of frame of reference to convert the correlated normals into uncorrelated variables via orthogonal transfor-
80
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
mation ([8], pp. 231–234, [4], pp. 354–368). In the process the correlation matrix is reduced to a diagonal matrix using eigenvalue and eigenvector techniques. This is in contrast with the present approach which leaves the correlation matrix intact. The present approach seeks the expanding dispersion ellipsoid which touches the limit state surface at the most probable point of failure (MPP) in the original coordinate system. As discussed earlier in connection with Eq. (1a) and (1b), that the tangent point is the most probable failure point is evident from the established form of the probability density function of the multinormal distribution: 1 1 T 1 exp ðx mÞ C ðx mÞ ð6aÞ f ðxÞ ¼ n 2 ð2 Þ2 jCj0:5 1 1 2 ¼ exp n 2 ð2 Þ2 jCj0:5
ð6bÞ
where is defined by Eq. (1a) or (1b), without the ‘‘min’’. Since minimizing (or 2) means maximizing f(x), minimizing subject to x being on the limit state surface means the same as finding the most probable failure point on the limit state surface. This intuitive and visual understanding of the design point is an alternative to the more mathematical approach in [3]. For correlated nonnormals, the present approach seeks the equivalent normal dispersion hyperellipsoid tangent to the limit state surface in the original coordinate system: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T N x i mN 1 xi mi i ð7Þ ½R ¼ min x2F iN iN N where mN i and i are the mean and standard deviation of the equivalent normal distribution based on the Rackwitz–Fiessler two-parameter equivalent normal transformation [Eqs. (3) and (4)], a transformation used also in FORM. However, unlike FORM, the equivalent ellipsoid approach does not need to diagonalize the correlation matrix or the covariance matrix, nor does it need to repeatedly reformulate the reduced limit state function in the transformed space. The integrated, expedient and transparent nature of the present approach can perhaps be appreciated by noting that the operational procedure (from the user-point of view) for correlated nonnormals differs from that of the correlated normals only in the distribution names entered in the first column of Fig. 5, and the corresponding parameters (Para1–Para4) characterizing the nonnormal distributions. It may be mentioned that the expanding ellipsoid approach [based on Eq. (1b)] and the expanding equivalent ellipsoid approach [based on Eq. (7)] in the original space of the variables are not vulnerable to numerical ill-conditioning even when the random variables show several orders of magnitude N difference in their magnitude, by virtue of the form (ximi)/ i in Eq. (1b) and (ximN i )/ i in Eq. (7). The significant operational advancement and improved clarity of the present ellipsoidal approach over the Low and Tang [16] ellipsoidal approach can be appreciated by comparing the reliability analyses of Fig. 5 (which invokes the codes of Fig. 4) with the reliability analysis of Fig. 2 (which implements the Rackwitz–Fiessler transformation directly on the worksheet). Had the seven cases of Fig. 5 been analyzed by the approach of Fig. 2, a new template would have to be devised for each different nonnormal distribution to be treated. Even increasing the number of
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
81
variates would require tedious multi-cell equation changes. It may also be pointed out that had one followed the widely documented classical FORM procedure of seeking the shortest distance in transformed space, with eigenvalues and eigenvectors for orthogonal transformation, and iteratively reformulating the reduced limit state function for each trial equivalent normal transformation, the procedure and codes (in the same EXCEL platform) would have been far more complicated and less transparent than the ellipsoid approach presented herein, even though the final computed index will be the same by both approaches.
7. On the need for positive definite correlation matrix Although the correlation coefficient between two random variables has a range 14 ij41, one is not totally free in assigning any values within this range for the correlation matrix. This was explained in Ditlevsen p85 [2] where a 33 correlation matrix with 12= 13= 23=1 implies a contradiction and is hence inconsistent. A full-rank (i.e., no redundant variables) correlation matrix has to be positive definite. One may note that Hasofer–Lind index is defined as the squareroot of a quadratic form [Eq. (1a) or (1b)], and unless the quadratic form and its covariance matrix and correlation matrix are positive-definite, one is likely to encounter the square root of a negative entity, which is meaningless whether in the dispersion ellipsoid perspective or in the perspective of minimum distance in the reduced variate space. In Fig. 5, the correlation matrix of the last four cases becomes singular when 23= 32=0.5. The quadratic form evaluates to a negative value when 23= 32 > 0.5. This means that given
12=0.5 and 13=0.5, the correlation coefficient 23 has to be smaller than 0.5. On the other hand, if 12= 13=0.8, the correlation coefficient 23 has to be greater than 0.3 in order for R to be positive definite. In the latter case it can be appreciated intuitively that when both the random variables Z2 and Z3 are highly and positively correlated with Z1 (via 12= 13=0.8), then Z2 and Z3 should also be positively correlated with each other to a certain extent ( >0.3). Mathematically, a symmetric matrix R is positive definite if and only if all the eigenvalues of R are positive. For present purposes, to detect inconsistent correlation matrix, one can create the p equation QForm=uTR1u side by side with = (QForm). During Solver’s searches, if the cell Qform shows negative value, it means that the quadratic form is not positive definite and hence the correlation matrix is inconsistent. The Monte Carlo simulation program @RISK will also display a warning of ‘‘Invalid correlation matrix’’ and offer to correct the matrix if simulation is run involving an inconsistent correlation matrix. That a full-rank correlation or covariance matrix should be positive definite is an established theoretical requirement [29–32], although not widely known. This positive-definiteness requirement is therefore not a limitation of the paper’s proposed procedure. Further reasons for positive definiteness follow. The joint p.d.f. of an n-dimensional multivariate normal distribution is given by Eq. (6a). The p.d.f. has a peak at mean-value point and decreasing ellipsoidal probability contours surrounding m only if the quadratic form and hence the covariance matrix and its inverse are positive definite. (By lemmas of matrix theory: C positive definite implies C1 is positive definite, and also means the correlation matrix R and its inverse are positive definite.)
82
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
Without assuming any particular probability distributions, it is documented in books on multivariate statistical analysis (e.g. [29–32]) that a covariance matrix (and hence a correlation matrix) cannot be allowed to be indefinite, because this would imply that a certain weighted sum of its variables has a negative variance, which is inadmissible, since all real-valued variables must have nonnegative variances. (The classical proof can be viewed as reductio ad absurdum, i.e, a method of proving the falsity of a premise by showing that the logical consequence is absurd.) An internet search using the Google search engine with ‘‘positive definite correlation matrix’’ as search terms also revealed thousands of relevant web pages, showing that many statistical computer programs (in addition to the @RISK simulation software) will check for positive definiteness of correlation matrix and offer to generate the closest valid matrix to replace the entered invalid one. For example, one such web site titled ‘‘Not Positive Definite Matrices–Causes and Cures’’ at http://www.gsu.edu/mkteer/npdmatri.html refers to the work of Wothke [32].
8. Reliability analysis of finite beam on Winkler medium Beams on elastic foundations can be analyzed using various analytical and numerical methods (Hetenyi [22], Selvadurai [23], Boresi et al. [24]). In this section reliability analysis will be performed for a beam on Winkler medium, to test robustness and accuracy of the present approach when relatively complicated nonlinear performance functions are involved. The beam has a length (l) of 15 m, a width (b) of 1 m. Six random variables are considered: the moment of inertia Ib (in m4) and Young’s modulus Eb (MN/m2) of the beam, the subgrade modulus ks (MN/m3), and the magnitudes (in MN) of three concentrated loads W1, W2, and W3, acting at 1.5, 7.5, and 13.5 m respectively from the left end of the beam. Failure is considered to have occurred if the bending moment reaches 1 MNm anywhere along the beam. The analytical solutions summarized in Young and Budynas [25] and shown in Fig. 6 are incorporated in a user-created function BMoment (. . .) whose input parameters include the load vector (W1, W2, W3), position vector (1.5, 7.5, 13.5 m), x coordinate and lamda (equal to (bks/4EbIb)0.25). The loads W1, W2, and W3 are assumed to have a coefficient of variation equal to 0.15. The assumed distribution types, correlation and other input are as shown in Fig. 7. The computed b index is 2.383. The reliability-based probability of failure (0.86%) compared very well with simulation (Pf=0.85%) using @RISK with Latin Hypercube sampling. Reliability analysis using Solver in Excel XP obtained the b index in 10 seconds (on a Pentium4 PC) after 9 trial solutions. In using Solver, the ‘‘automatic scaling’’ option has been activated.
9. Reliability analysis of complex strut involving implicit performance function The deterministic analysis of a strut with complex supports was discussed in Coates et al. [26]. The member is initially straight (Fig. 8) with a pin support at 2, an elastic support at 1 (rotational stiffness l1) which provides a restoring moment M1=l11, and a support at 3 (stiffness k3) which provides a reaction force (k3v3) proportional to the vertical displacement v3. The problem is to determine the smallest value of axial force P which will cause the strut to become elastically unstable. Coates et al. [26] showed that the problem reduces to that of finding a value of P, referred
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
83
to as Pcrt, which would make the determinant of a 7 by 7 matrix (functions of P, L, a, E, I, l1 and k3) vanish. It was also noted that a numerical procedure would be necessary to determine the determinant of the matrix by sequential variation of the axial load P until a zero value is reached. For the reliability analysis here, the critical load Pcrit can be obtained in a spreadsheet cell (Fig. 8) that calls the user-defined VBA function Pcritical (L, a, E, I, l1, k3) shown in Fig. 9. The performance function (PerFunc in Fig. 8) is ‘‘=Pcrt.P’’, where P refers to the first value of the x* column. Hence the single cell object ‘‘PerFunc’’ contains the implicit and iterative user-created routine of Fig. 9. For the distributions and input in Fig. 8, the initial (mean) values of the x* column were (700, 1000, 500, 200000, 200, 500, 10). Solver in Excel XP (2002) was then invoked to minimize the
Fig. 6. Bending moment equations for finite-length beams on Winkler elastic medium. (Ref. Young and Budynas, 2002).
84
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
cell, by changing the x* values, subject to the constraint that PerFunc=0. The Solver option ‘‘Use Automatic Scaling’’ was activated. Solver obtains the solution in 35 s (in a Pentium4 PC) after 14 trial solutions. The probability of failure (0.40%) based on the index of 2.652 compares reasonably well with the 0.36% failure probability from simulation with 250,000 trials. Note that in this case where there is five orders of magnitude difference between random variables E and k3, use of ‘‘Automatic Scaling’’ is essential for efficiency. The design point (x*) and the nx values in Fig. 8 indicate that the parameters l1 and k3 hardly deviate from their mean values of 500 and 10 respectively. The implication is that the performance function and hence Pcritical is not sensitive to l1 and k3 at their current mean values and
Fig. 7. Reliability analysis of asymmetrically loaded beam on Winkler medium, and comparison with Monte Carlo simulation.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
85
levels of uncertainty. These rather unexpected findings were further investigated in four analyses in which the initial values of l1 and k3 were set at two standard deviations above and below their respective mean values, while the other five parameters have initial values at their respective mean values. The indices (vary between 2.6513 and 2.6524) and the design point x* values obtained after constrained optimization were virtually identical to those of Fig. 8, thereby providing further support to the finding that Pcritical is insensitive to l1 and k3. Coates et al. [26] also investigated sensitivities, deterministically, by plotting curves of Pcritical versus k3, for values of l1 ranging from 0 to 1020 Nmm/ rad, and found that very large changes in l1 had only very small influence on Pcritical. In comparison, for the mean values and standard deviations of Fig. 8, reliability analysis reveals the insensitivities of Pcritical to l1 and k3 with much less effort than deterministic contour plotting, besides providing additional information, including a good estimate of the probability of strut failure.
Fig. 8. Reliability analysis of a strut with complex supports. Performance function=(Pcrt.P), where Pcrt. is programcalculated by the code of Fig. 9.
86
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
10. Limitations and flexibilities of object-oriented constrained optimization in spreadsheet There is no closed form equation for the inverse of the standard normal cumulative distribution. Numerical procedures such as that incorporated in Excel’s NormSInv is necessary. In Excel 2000 and earlier versions, NormSInv(CDF) returns correct values only for CDF between about 0.00000035 and 0.99999965, i.e., within about 5 standard deviations from the mean. There are undesirable and idiosyncratic sharp kinks at both ends of the validity range. In EXCEL XP (2002), the NormSInv (CDF) function has been refined, and returns correct values for a much broader band of CDF: 1016
11. Summary and conclusions An expanding ellipsoid perspective (Low and Tang [16]) provides an alternative and intuitive interpretation of the Hasofer–Lind index. The perspective led to a practical spreadsheet optimization procedure for performing reliability analysis. This paper extends the 1997 approach by providing a flexible user interface (Figs. 3 and 4) in which numerous different nonnormal distributions can be selected with ease for reliability analysis, and by testing the robustness and
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
87
accuracy of the approach for two relatively complicated problems: a beam on elastic foundations, and a strut with complex supports. The failure probabilities based on reliability indices are in good agreements with failure probabilities from Monte Carlo simulations (Figs. 5, 7, and 8). The main features of the procedure presented herein for reliability analysis involving correlated nonnormals are: The design point and reliability index are obtained using object-oriented constrained optimization of equivalent hyperellipsoids in the original space of the random variables.
Fig. 9. An example of program-calculated performance function that is incorporated in a single cell object in Fig. 8.
88
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
The approach obtains the same result as FORM, but is operationally more direct and transparent. There is no need to diagonalize the correlation matrix. The platform is ubiquitous. The Rackwitz-Fiessler equivalent normal transformations are executed in a short program code (Fig. 4). This renders possible a consistent and uncluttered user-interface (Figs. 3, 7, and 8). To specify different nonnormal distributions, it is only necessary to name the distributions and provide its parameters. Performance functions can be implicit (e.g. containing VBA codes) and involving various correlated nonnormal variates, as in Figs. 8 and 9, for example. In general, for reliability analysis involving correlated nonnormal distributions, the accuracy of the present approach may depend on the nature of the performance functions, the number of random variables, the validity range of the NormSInv function (which is used for equivalent normal transformation), and the robustness of the nonlinear constrained optimization program. For the cases of correlated nonnormals analyzed in this study, Excel XP (with its refined NormSInv function) and its built-in standard Solver are found to be efficient and accurate (as far as FORM is concerned). The main advantage of the object-oriented constrained optimization approach for reliability analysis is in the modular nature of its constituents: The deterministic model that underlies the performance function g(x) is set-up (or programmed in VBA) separately and appears as one or more spreadsheet cell objects. The nonlinear constrained optimization program Solver is then invoked to minimize the cell object (itself a conglomerate of program routines), subject to constraints on the performance function cell object and other cell objects. The platform can be other than Microsoft Excel. There are also other choices of potentially more robust spreadsheetbased optimization programs, though less ubiquitous and at a cost.
References [1] Hasofer AM, Lind NC. Exact and invariant second-moment code format. J of Engrg Mech ASCE, New York 1974;100(1):111–21. [2] Ditlevsen O. Uncertainty modeling: with applications to multidimensional civil engineering systems. New York: McGraw-Hill; 1981. [3] Shinozuka M. Basic analysis of structural safety. Journal of Structural Engineering, ASCE 1983;109(3):721–40. [4] Ang HS, Tang WH. Probability concepts in engineering planning and design, v2: decision, risk, and reliability. New York: John Wiley; 1984. [5] Madsen HO, Krenk S, Lind NC. Methods of structural safety. NJ: Prentice-Hall: Englewood Cliffs; 1986. [6] Tichy M. Applied methods of structural reliability. Dordrecht, Boston: Kluwer Academic; 1993. [7] Melchers RE. Structural reliability analysis and prediction. 2nd edition. John Wiley; 1999. [8] Haldar A, Mahadevan S. Probability, reliability and statistical methods in engineering design. New York: John Wiley; 1999. [9] Nowak AS, Collins KR. Reliability of structures. McGraw-Hill; 2000. [10] Der Kiureghian A, Lin HZ, Hwang SJ. Second-order reliability approximations. Journal of Engineering Mechanics, ASCE 1987;113(8):1208–25. [11] Chen X, Lind NC. Fast probability integration by three-parameter normal tail approximations. Structural Safety 1983;1:269–76.
B.K. Low, W.H. Tang / Structural Safety 26 (2004) 69–89
89
[12] Wu YT, Wirsching PH. New algorithm for structural reliability estimation. Journal of Engineering Mechanics, ASCE 1987;113(9):1319–35. [13] Zhao, YG, Tetsuro Ono. A general procedure for first/second-order reliability method (FORM/SORM). Structural Safety 1999;21:95–112. [14] Zhao YG, Tetsuro Ono. Moment methods for structural reliability. Structural Safety 2001;23:47–75. [15] Rackwitz R. Reliability analysis-a review and some perspectives. Structural Safety 2001;23:365–95. [16] Low BK, Tang WH. Efficient reliability evaluation using spreadsheet. Journal of Engineering Mechanics, ASCE 1997;123(7):749–52. [17] Rackwitz R, Fiessler B. Structural reliability under combined random load sequences. Computers & Structures 1978;9:484–94. [18] Hahn GJ, Shapiro SS. Statistical models in engineering. John Wiley & Sons; 1967. [19] Vose D. Risk analysis: a quantitative guide. 2nd edition. John Wiley & Sons; 2000. [20] Evans M, Hastings N, Peacock B. Statistical distributions. 3rd edition. New York: Wiley; 2000. [21] Altman DG. Practical statistics for medical research. London: Chapman and Hall; 1991. [22] Hetenyi M. Beams on elastic foundations. Ann Arbor, Michigan: Univ of Mich Press; 1974. [23] Selvadurai APS. Combined foundations. In Henry FDC, editor. The Design and Construction of Engineering Foundations (2nd ed.). London: Chapman & Hall; 1986, chapter 5. [24] Boresi AP, Schmidt RJ, Sidebottom OM. Advanced mechanics of materials. 5th edition. New York: Wiley; 1993. [25] Young WC, Budynas RG. Roark’s formulas of stress and strain. 7th edition. McGraw-Hill; 2002. [26] Coates RC, Coutie MG, Kong FK. Structural analysis. 3rd edition. London: Chapman & Hall; 1994. [27] Low BK. Probabilistic slope analysis involving generalized slip surface in stochastic soil medium. In Proceedings of the 14th South East Asian Geotechnical Conference, 10-14 December 2001, Hong Kong, A.A.Balkema Publishers; 2001. p.825–830. [28] Low BK, Teh CI, Tang WH. Stochastic nonlinear p-y analysis of laterally loaded piles. Proceedings of the Eight International Conference on Structural Safety and Reliability, ICOSSAR ‘01, Newport Beach, California, 17–22 June 2001. 8 pages. [29] Chatfield C, Collins AJ. Introduction to multivariate analysis. London: Chapman and Hall; 1989 (reprint). [30] Graybill FA. Matrices with applications in statistics. 2nd ed. Belmont, CA: Wadsworth; 1983. [31] Searle SR. Matrix algebra useful for statistics. New York: John Wiley; 1982. [32] Wothke W. Nonpositive definite matrices in structural modelling. In: Bollen KA, Long JS, editors. Testing structural equation models. Newbury Park, CA: Sage Publications, Inc; 1993. p. 257–93.