Estimation of statistical moments of performance functions using efficient importance measure

ABSTRACT This paper introduces an efficient methodology for evaluating the first four central moments of performance functions in structural reliability assessment. The procedure of the proposed method is divided into three major steps. First, a performance function is approximated as the sum of univariate and bivariate functions by the bivariate dimension-reduction method. Next, a new criterion based on importance measure analysis is proposed for dividing the bivariate functions into significant bivariate functions and non-significant bivariate functions, and then the new approximation of the performance function is obtained. Finally, formulas are derived for evaluating the first four central moments of the proposed performance function, in which the first four raw moments of the univariate and bivariate functions can be achieved by a point estimating method in a standard normal space. The effectiveness and efficiency of the proposed method are demonstrated through three numerical examples.


Introduction
The moment method is one of the effective techniques for analyzing structural reliability (Fan, Ang, and Li 2016;Wan et al. 2020), which is known for its simplicity and the avoidance of the difficulties involved in searching design point by iteration using the derivatives of the performance function (Zhao and Ono 2001). With the development of moment-based structural analysis, many proposals are able to provide accurate prediction of failure probability by using the obtained first few moments of the performance function (Rajan et al. 2016). However, it is not easy to estimate statistical moments of the performance function in practical engineering as it might be quite intricate, especially in dealing with multidimensional variables and strong nonlinearity.
Statistical moments are typically defined in terms of multi-dimensional integrals. For some simple explicit functions, such as linear sum, it is easy to obtain the moments directly from the integral-based definitions. For complicated and implicit functions, two-point and three-point estimating methods with weighted sum of functions were proposed by Rosenblueth (1975) and Gorman (1980), respectively. These methods are quite simple and do not require the computation of derivatives of performance functions; however, the accuracy of these methods in general cases is quite low, and the estimate points may be outside of the definition domains of random variables. In order to remove these weaknesses, Zhao and Ono (2000) proposed a new point estimating method in standardized normal space, in which the performance function is approximated by the sum of some one-dimensional functions. Though the evaluation accuracy of a higher moment is unoptimistic in many cases when facing strong nonlinearity, the computational efficiency of moment evaluation is greatly improved Li and Zhang 2011). A generalized dimensionreduction method was proposed by Xu and Rahman (2004), demonstrating its application to the momentbased reliability analysis by combining with the point estimating method in standardized normal space (Zhao and Lu 2011). Usually two-dimensional integrals are considered in generalized dimension-reduction methods, and the precision of estimating moments is effectively improved when solving some strong nonlinear performance functions. However, the computational efforts are not acceptable if a large number of bivariate integrals are involved (Fan, Ang, and Li 2016). Besides, excessive terms originated from bivariate integrals also increase the possibility of errors. Hence, there remains a contradiction between accuracy and efficiency in evaluating moments by dimensionreduction methods (Xu and Dang 2019). There are mainly two approaches available for improving this deficiency, one is to optimize selected points and the other is to modify a model based on cross terms. In the first approach, alternatives are sought to evaluate moments for two-dimensional integrations instead of using Gaussian spatial points, such as the sparse grid (He, Gao, and Gong 2014), number-theoretic net (Chen and Li 2008), and the cubature formulas (Cools 2003). Recently, an adaptive bivariate dimension-reduction method (ABDRM) has been proposed from another perspective (Fan, Ang, and Li 2016), in which some portions of the two-dimensional functions can be decomposed as the sum of one-dimensional functions according to the existing criterion of cross terms. After that, seeking better precision, different point-selection methods, such as the sparse grid and cubature formulae, are applied to compute two-dimensional integrals based on that adaptive bivariate dimension-reduction method (Cai et al. 2019). The cross-term-based methods provide an explicit criterion to decompose dimension-reduction integrals, thus a moment algorithm is optimized with fewer computational efforts. However, it is not easy to identify the interaction with each variable, especially for the performance functions corresponding to the assessment of the system reliability with multiple failure modes. Furthermore, even if the cross terms can be identified, the improvement in the computing efficiency cannot be guaranteed when the number of cross terms in the performance functions is large. For such cases, the criterion of decomposing twodimensional integrals needs to be investigated further.
In order to address the limitations identified above, this paper develops a new procedure for moment evaluation based on the efficient importance measure analysis (IMA). The rest of this paper is organized as follows. In section 2, the bivariate dimension-reduction method (BDRM) for approximating performance functions is briefly reviewed, and its main drawbacks are summarized. Based on IMA, the bivariate functions in the performance function are decomposed as significant ones and nonsignificant ones. The non-significant ones are approximated further using univariate functions. The proposed performance function is then expressed as a sum of these decomposed low-dimensional functions. In section 3, the formula used for estimating moment in the proposed performance function is derived. In section 4, the efficiency and accuracy of the proposed method over the point estimate method based on BDRM are investigated. Furthermore, three illustrative examples are then presented to demonstrate the efficiency and effectiveness of the proposed method in section 5. Finally, the findings of this paper are concluded in Section 6.

General approximation using bivariate dimension-reduction method
For the performance function Z = G(X), G(X) can be approximated, as follows (Cai et al. 2019), by using the inverse Rosenblatt transformation (Hohenbichler and Rackwitz 1981) and bivariate dimension-reduction method .
where i, j = 1, 2, . . ., n and i < j; n is the dimension of random variables; T À 1 � ð Þ is the inverse Rosenblatt transformation (Hohenbichler and Rackwitz 1981); μ X i represents the mean of X i ; G μ is a constant; G U i ð Þ denotes a univariate function with random variables of U i ; and G U i ; U j À � is a bivariate function with the random vector of (U i , U j ). Then, the statistical moments of a structural performance function can be computed by estimating the moment of each low-dimensional function.
Obviously, G(X) is approximated by summing up ð Þ and a constant value. With increasing n, the number of bivariate functions grows quadratically in dimensions, which generally results in excessive computational burdens and possibility of accumulated errors in the estimation of moments of bivariate functions. Therefore, it is necessary to reduce the number of bivariate functions.
A simple linear expression can be used to approximate using its corresponding univariate functions as given below : In this way, the performance function in Equation (1) can be directly degenerated, which is also known as the univariate dimension-reduction method (UDRM). However, the impact of the inherent approximation error has to be considered in Equation (3) when it is used to reduce the number of G(U i ,U j ). Note that each bivariate function in Equation (1) can be regarded as an integrated random input, and a different random input generally has a different contribution to the random output. If the bivariate function G U i ; U j À � has very little interference in the uncertainty of the output (i.e., a very small range around the output value G μ ), its corresponding approximation error in Equation (3) can be neglected in a random output. Thus, in order to further optimize the performance function in Equation (1), the global sensitivity analysis is necessary for quantifying and comparing the relative importance of the bivariate functions.

Importance measure analysis for bivariate functions
The importance measure analysis (IMA), also known as the global sensitivity analysis, can be regarded as an extension of the probabilistic local sensitivity analysis. It is defined as the contribution degree of uncertainties contained in random variables to the overall uncertainty (Sobol 2001). One of the most well-known methods for IMA is the variance-based method, referred to as the analysis of variance (Zhou, Lu, and Li 2013). In order to measure the importance of bivariate functions based on the approximation in Equation (3), we firstly investigate the corresponding univariate functions.
Based on the Sobol's decomposition theory (Sobol 1993;Shayan 2013), the significant index (λ i ) for univariate function G(U i ) can be expressed as , which can be given as follows: where ϕ(·) denotes the probability density function (PDF) of the standard normal variable. Using the point estimating method in the standard normal space (Zhao and Ono 2000), the above integral can be approximated as follows: In Equation 6(b), u a and P a are the estimating points and weights, respectively (see Appendix I for detail).
For an n-dimensional random vector, let the classification criterion � λ be defined as the mean of significant indices as follows: In practical structural performance functions, if λ i � � λ � ε (ε ¼ 10 À 4 is considered in this study), then the corresponding variable X i (or U i ), which is called as the ineffective variable, should be dropped from the variable vector X (or U) and considered as a constant parameter rather than a random variable in the performance function G(X) (or G(U)). Then, the number of variables n is reduced to n* (n* ≤ n), and Equation (7a) can then be expressed as: where n* is the number of effective random variables. Therefore, the practical criterion for classifying G(U i ) is simply defined as follows: the significant univariate function GðU þ i Þ (i.e.,λ i > � λ) and the non-significant univariate function GðU À i Þ (i.e.,). Next, consider the bivariate function GðUÞ ¼ GðU i ; U j Þ(i < j). Similar to λ i , the significant index (K i,j ) for bivariate functions GðUÞ can be expressed as According to Equation (3), the obtained variance of the one-dimensional function is utilized to approximate as Then, based on Equations (8) and (9), K i , j can be rewritten as where P i < j K i;j ¼ 1. It can be observed that there exists an explicit relationship between the proposed two significant indexes of λ i and K i , j .
Using the mean value of K i , j as a classification boundary like � λ, we get where n b is the total number of the bivariate functions G(U i , U j ). According to Equations (10) and (11), the bivariate function G(U i , U j ) can simply be separated as the significant bivariate function GðU þ i;j Þ (i.e., K i;j > K) and the non-significant bivariate function If the number of GðU þ i;j Þ is considered to be n þ b , such a decomposition rate can be defined as

Proposed approximation of performance functions using methods of dimension-reduction integration and importance measure analysis
According to the criteria in Equation (11), the performance function G(X) in Equation (1) can further be decomposed as where the sum of non-significant bivariate functions P GðU i ; U j Þ can be decomposed using Equation (3) to reduce the number of bivariate functions as follows: X 3. Statistical moments of the proposed approximation of performance functions

Moment estimation formula based on point estimating method
In order to propose a clear formula for calculating moments, the first four raw moments of H where 2,3,4) are the kth raw moments of the bivariate and univariate functions, respectively, which can be obtained by using the expressions of the point estimating method in standardized normal space. The approximate results of Usually the number of estimating points m = 7 is applied. Then, Equations (17a) and (17b) can be used to obtain the first four raw moments of H k h i G of the performance function. Finally, the first four statistical moments, i.e., the meanμ G , standard deviationσ G , skewness α 3G and kurtosisα 4G , can be derived as

Procedure of the proposed method for moment estimation
From the preceding investigation, the procedure of the proposed method for estimating the first four statistical moments of performance functions is summarized in Figure 1.

Efficiency and accuracy of the proposed method
It can be observed in Equations (16b) and (16c) that the computation of the performance functions in H after the decomposition, and hence, they share the same computational efforts. Therefore, the maximum number of computations of the performance functions in the proposed method is Upon combining with the point estimating method in Equation (1), n s ¼ n n À 1 ð Þ m 2 À 2m ð Þ=2 þ mn þ 1 points will be required for determining the first four moments of G(X) by using the m-point estimating method, which is known as PE-BDRM. Compared with PE-BDRM, the efficiency improvement index of the proposed method can be described as: Similarly, (mn+1) points are needed for the point estimating method based on UDRM (known as PE-UDRM). According to Equation (20), the proposed method approaches PE-UDRM when γ b decreases to 0%, while it is degenerated from PE-BDRM when γ b increases to 100%. Besides, Equations (17a) and (17b) are approximated equations and could have small random errors. Hence, they can be rewritten as: where ε U i; j ;k (i < j) and ε U i ;k are assumed as the random errors with the same variance, which are caused by the kth moment approximation of low-dimensional functions. Based on PE-BDRM, the total random approximated errors of the kth raw moment approximation, Θ k , can be described as In order to assess the uncertainty of Θ k , the total variance V Θ k can be given based on Equation (23) as follows: where V ε;k is the variance of ε U i ; j ;k (i < j) or ε U i ;k , which is usually a small positive number. Here, we define V Θ k � V ε;k as a variance amplification factor, which increases significantly with increasing n, thus resulting in a huge impact on the accuracy of the output in Equation (23).
Since the decomposition errors of the nonsignificant bivariate functions in Equation (14) are negligible, Θ k in Equation (22) can be modified as Note that the random error P can be eliminated by combining similar terms of ε U i ;k as they come from the same estimation of the kth raw moment of one-dimensional functions. The corresponding total variance V Θ k of the proposed method is obtained as: To make a comparison with Equation (23), suppose n* = n in Equation (25) (all the random variables are considered to be effective), and Equation (25) has the same expression as Equation (23) when γ b = 0%. The relationship between n and V Θ k =V ε;k under different decomposition rates is shown in Figure 2. It can be observed from Figure 2 that when the approximate errors of the non-significant bivariate functions are ignored in Equation (14), the proposed method provides smaller uncertainty to the kth estimation of raw moments, especially with increasing n.
In comparison to PE-UDRM, PE-BDRM can provide better accuracy when strong non-linearity or nonnormality is involved as the interaction between each variable is taken into account. However, PE-BDRM still faces the problem of both accuracy and efficiency when the number of cross terms is large in the performance functions. Hence, the proposed method, from a perspective of importance measure, finds a road to achieve a balanced computational efficiency and accuracy in the estimation of moments.

Numerical examples
Three numerical examples are investigated to demonstrate the effectiveness and advantages of the proposed method over traditional methods. The first example uses a simple nonlinear performance function for the same.
Two typical performance functions for truss structure systems are used in the second example, which shows obvious superiority of the proposed method compared with cross-term-based methods in terms of function decomposition. Dozens of random variables are considered in the third example with an implicit performance function of a six-storey three-bay ductile frame structure.

Example 1: a simple nonlinear performance function of an annular section column
An annular section column subjected to an axial compressive load is considered (Huang and Zhang 2013), which is shown in Figure 3. The performance function is given as where the random information about the involved variables is listed in Table 1. According to Equation (1), the performance function in this example can be approximated as Using Equation (4), the significant indices of each onedimensional function can be obtained as λ E = 0.732, λ P = 0.166, λ L = 0.052, λ D = 0.018, and λ T = 0.032. As it can be observed, the first two univariate functions G(U E ) and G(U P ) are significantly important than the others, and there is no ineffective variable in this case (that is, n* = n = 5). Using Equation (10), the significant indices of each bivariate function can be obtained as K E , P = 0.225, K E , L = 0.196, K E , D = 0.188, K E , T = 0.191, K P , L = 0.055, K P , D = 0.046, K P , T = 0.049, K L , D = 0.018, K L , T = 0.021, and K D , T = 0.012. According to Equation (11), the classification boundary for two-dimensional functions can be calculated as G(U P ,U L ), G(U P U D ), G(U P ,U T ), G(U L ,U D ), G(U L ,U T ), and G(U D ,U T ) are then identified as the non-significant bivariate functions. Based on Equation (13), the decomposed performance function with a decomposition rate of γ b = 60% can be rewritten as According to Equations (16a-d) and (17a-b), the first four raw moments of both one-and two-dimensional functions could be obtained, which are presented in Tables 2 and 3, respectively. According to Equations (18a-d), the first four statistical moments of Equation   (26) and the number of the computations of the performance function can be obtained as listed in Table 4 together with those determined from different methods, namely MCS, PE-BDRM, PE-UDRM and ABDRM for judging cross terms (Fan, Ang, and Li 2016). The results of MCS are taken from 10 6 runs as the benchmark values for comparison. It can be observed from Table 4 that (1) although the PE-UDRM is very efficient, its predicted skewness is far away from the benchmark value due to the existence of cross terms in Equation (26); (2) the results obtained from both the PE-BDRM and the proposed method are almost the same with the same benchmark values; and (3) the computational effort of the proposed method is reduced with an efficiency improvement index of ϴ = 54.5%; (4) because there are 6 different cross terms in the performance function, the decomposition rate of ABDRM is 40% (smaller than 60% in the proposed method), and the proposed method is more efficient in this example.

Example 2: performance functions of truss structures containing multiple failure modes
The second example considers performance functions corresponding to truss structures (see Figure 4) having multiple failure modes. Two cases of the example are considered: (a) a simple brittle truss (non-series system) and (b) an elastoplastic truss with two stories and one bay (series system). The statisticals information about the strength and load of the member is listed in Table 5.
For the truss in case (a), the corresponding performance functions for each of the failure modes are given as follows (Zhao and Ang 2003): The performance function can be defined as the minimum of the above, i.e., The proposed procedure is performed for Equation (32) and the number of effective variables is obtained as n* = 4, which means that W 1 , W 2 and W 5 can be deemed to be constant (taken their mean values). The computational results are presented in Table 6. It can be observed that: (1) the proposed method provides results almost similar to those of PE-BDRM, and the efficiency of the proposed method is greatly improved with an efficiency improvement index of ϴ = 80.25%; and (2) the results of the proposed method are closer to those of MCS than those of PE-UDRM; (3) because there are no cross terms in this performance function, ABDRM is equivalent to PE-UDRM. For the elastoplastic truss structure in Case (b), the performance functions corresponding to the eight most likely failure modes are given below as per the work of Ono et al. (1990).

Example 3: an implicit performance function of six-storey three-bay ductile frame system
The third example considers a six-storey three-bay frame under concentrated lateral loads as illustrated in Figure 5. The statistical properties of 28 random variables (n = 28) including member strengths (M 1 , . . ., M 22 ) and loads (S 1 , . . ., S 6 ) are listed in Table 9.
The performance function for this example is expressed as where λ(M, P) is the load factor, which can be obtained from the analysis of structural limit using the Compact Procedure (CP) (Aoyama 1988).
After the filtering process, 16 random variables having almost no effect on the outcome of the performance function can be eliminated (n* = 12). The first four moments obtained by different methods are presented in Table 10. It can be observed in Table 10 that: (1) PE-UDRM is found again to be very efficient. However, the predicted skewness and kurtosis are far away from those of MCS; (2) compared with PE-BDRM, the efficiency improvement of the proposed method is    Figure 5. A six-storey three-bay ductile frame structure for Example 3.

Concluding remarks
The focus of this study was to evaluate the statistical moments for structural performance functions in an accurate and efficient way. On the basis of dimensionreduction methods, a novel method combined with importance measure analysis is proposed. In this method, the approximated performance function based on BDRM is optimized through the decomposition of non-significant bivariate functions, which are identified according to the proposed classification criterion. Point estimating method is applied to calculate the first four moments of the proposed performance function, in which the computational efficiency and estimating errors are also analyzed. Three examples, including non-linear, piecewise and implicit performance functions which usually relate to complex, systematic and multivariable structural reliability problems, are investigated to illustrate the efficiency and accuracy of the proposed method. Besides, PE-UDRM and PE-BDRM are also carried out for comparison, together with MCS providing some benchmark results. The work can be concluded as: (1) PE-UDRM, without any bivariate function, usually provides poor results in estimating moments of many non-linear and non-normal performance functions as the interactions between two random variables are ignored.
(2) Although PE-BDRM considers the effectiveness of cross terms in the performance function, it still faces intricate difficulties in efficiency when the number of involved bivariate functions increases. (3) With high accuracy, the proposed method reduces the number of bivariate functions in most cases, and hence it becomes more efficient than PE-BDRM; (4) The proposed method makes the total variance of random errors, caused by the point estimating method, smaller than PE-BDRM, especially with increasing number of variables.
In summary, the proposed method maintains a good balance between PE-BDRM and PE-UDRM by using a simple importance measure analysis. The proposed method is expected to be a useful tool in structural reliability analysis. Based on the proposed performance function, further investigations can be made to improve the accuracy and efficiency when different point selection methods are used (e.g., the sparse grid stochastic collocation method and the number-theoretic method).

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
The research reported in this paper is partially supported by the National Natural Science Foundation of China (Grant Nos.51820105014, 51738001, and U1934217). The support is gratefully acknowledged.

Notes on contributors
Dong-Zhu Hu is currently a Ph. D. student at Central South University, China. His research interesting mainly focuses on reliability assessment of structural system.

Zhao-Hui
Lu is currently a professor at Beijing University of Technology and former professor at Central South University, China. His expertise includes structural reliability and timedependent reliability based risk-cost optimised maintenance strategy. Yan-Gang Zhao is currently a professor at Kanagawa University, Japan and a guest professor at Beijing University of Technology, China. His research interesting mainly focuses on structural reliability and seismic performance of engineering structures.

Appendix. The estimating points and weights in standardized normal space
According to Zhao and Ono (2000), the estimating points in a standard normal space u k (k = 1, . . ., m) and the corresponding weights P k can be expressed as . The estimating points u k and weights P k in a standard normal space (m = 7).
Estimating points u k u 1 = -u 7 u 2 = -u 6 u 3 = -u 5 u 4 -3.7504397 -2.3667594 -1.1544054 0.0 Weights P k P 1 = P 7 P 2 = P 6 P 3 = P 5 P 4 5.48269 × 10 -4 3.07571 × 10 -2 0.2401233 0.4571427 where x k and w k are the abscissas and weights for the Hermite integration with the weight function exp(-x 2 ) that can be found in Abramowitz and Stegum (1972). For convenience, the abscissas and corresponding weights of a seven-point estimate (m = 7) in a standard normal space (Zhao and Ono 2000) are listed in Table A1.