Global Reliability Sensitivity Analysis Based on Maximum Entropy and 2-Layer Polynomial Chaos Expansion

To optimize contributions of uncertain input variables on the statistical parameter of given model, e.g., reliability, global reliability sensitivity analysis (GRSA) provides an appropriate tool to quantify the effects. However, it may be difficult to calculate global reliability sensitivity indices compared with the traditional global sensitivity indices of model output, because statistical parameters are more difficult to obtain, Monte Carlo simulation (MCS)-related methods seem to be the only ways for GRSA but they are usually computationally demanding. This paper presents a new non-MCS calculation to evaluate global reliability sensitivity indices. This method proposes: (i) a 2-layer polynomial chaos expansion (PCE) framework to solve the global reliability sensitivity indices; and (ii) an efficient method to build a surrogate model of the statistical parameter using the maximum entropy (ME) method with the moments provided by PCE. This method has a dramatically reduced computational cost compared with traditional approaches. Two examples are introduced to demonstrate the efficiency and accuracy of the proposed method. It also suggests that the important ranking of model output and associated failure probability may be different, which could help improve the understanding of the given model in further optimization design.


Introduction
Sensitivity analysis (SA) is one of the most important methodologies dealing with optimization in engineering practice. It aims to quantify importance in the output of a model with regards to some input parameters or design variables [1][2][3][4]. With the results of SA, important variables of given model are chosen for further study when less important variables may be treated as constant. This could help reduce the dimension of inputs (number of input variables) and improve the model understanding.
Generally, SA methodologies can be classified into local sensitivity analysis (LSA) and global sensitivity analysis (GSA) [1,5]. The former only concentrate on the influence of single input parameter on model output by gradient estimating in a small neighborhood around a normal point. It is often used to provide deterministic direction to be close to the optimal object [6]. Meanwhile, the latter one tries to quantify the output uncertainty when all input parameters are varying over entire value domain. And it does not depend on the choice of a nominal point [5,7]. Thus, GSA is suitable for uncertainty quantification (UQ) and reliability based design optimization (RBDO) [8]. 2 of 22 Numerous methodologies are developed for GSA, among which variance-based methods are widely used [1,5]. Variance-based methods take the well-known variance decomposition formula V(Y) = V(E(Y|X)) + E(V(Y|X) as their foundation. This formula has variance separated by using the variance of conditional expectation and associated residual variance. Thus, the contribution or importance of X on Y could be further studied. Many researchers take Sobol's indices because they could provide a clear, accurate and robust indication of the selected uncertain inputs both quantitatively and qualitatively [9]. Although Sobol's indices are traditionally applied for non-statistical model output, some researchers have extended Sobol's indices to identify uncertainty effects of inputs on reliability. Morio only studied the effects of some variable distribution parameters, such as variable means, but it did not consider the variable uncertainties themselves [10]. Borgonovo proposed the moment-independent global sensitivity indices which are derived from conditional probability density function (PDF) [11,12]. Lu accepted this idea and further developed the variance-based importance measurement for global reliability sensitivity analysis (GRSA) with associated global reliability sensitivity indices [13,14]. Cui and Li proved the equivalence relation between these global reliability sensitivity indices and Sobol's indices when taking the indicator function as the assumed model output.
Unlike traditional global sensitivity indices calculation, the global reliability sensitivity indices aim at calculating importance of input variables on reliability that is a statistical parameter. These indices are specifically concerned about the variance of conditional expectation of indicator function according to variance decomposition formula. However, the indicator function maps every single model output to only two real numbers 0 or 1 discretely, so the analytical method and surrogate method, which are popular with uncertainty analysis, have not been developed for GRSA because it is difficult to obtain an analytical formula of the indicator function. In this way, the Crude Monte Carlo Simulation (MCS) is always employed to calculate such measurement [15,16], which leads that the computation cost is unacceptable when the given model is complex. Although Li et al. used the importance sampling technologies to improve sampling efficiency [13], quite a few samples are inevitable.
Fortunately, the conditional expectation of indicator function is equal to the conditional failure probability (CFP) that could be regarded as a continuous random variable. Thus, the variance of conditional expectation of indicator function could be replaced by the variance of CFP. In this way, if CFP could be assumed as the model output, there might be an unknown expression between CFP and corresponding conditional variable. So, the surrogate method could be used to approximate the CFP. Then, the statistical properties of CFP, such as mean and variance, may be calculated by surrogate model efficiently.
As for the issue of the statistical properties of a surrogate model, polynomial chaos expansions (PCE) technology provides a powerful tool for this application [17,18]. The PCE is a surrogate modelling technology with orthogonal polynomial basis functions, and it has statistical moments calculation simplified. In other words, variance information of surrogate model could be readily obtained by post-processing PCE's coefficients analytically [19,20]. Therefore, it is suitable to make a surrogate model for CFP and calculate its variance. However, a challenge for PCE method in such application is that quite a few specific probability values are required. In addition, since the CFP is dependent on a set of conditional variables, each global reliability sensitivity index corresponding to different conditional variables requires a specific surrogate model. This may be time consuming when computing several indices. To circumvent this problem, the maximal entropy (ME) method is introduced since it could compute probabilities precisely and analytically only with prior moment information of model output [21][22][23][24]. Considering the orthogonal characteristics, the PCE could not only provide the first two order moments directly, but also provide higher order moments by PCE multiplication algorithm without additional sampling. In this way, the combination of PCE and ME could make GRSA efficiently. To the author's knowledge, this idea has not been developed for global reliability sensitivity indices calculation. Therefore, a new 2-layer PCE algorithm coupling PCE and ME is developed. That is, the first layer PCE is to generate a surrogate model by traditional algorithm. The second layer PCE is to construct several surrogate models for CFPs with the help of ME, which are referred as reliability polynomial chaos expansions (R-PCEs). The moment information that ME requires is provided by the first layer PCE when the conditional variables are given and fixed. Then, the variance of each CFP comes out by analytically post processing R-PCE coefficients, followed by global reliability sensitivity indices.
The presentation of the work is structured as follows: Section 2 provides a brief introduction of GSA and GRSA respectively. Section 3 presents the proposed 2-layer PCE framework and detailed calculation procedure. A numerical example and an engineering example are presented to test the performance of the proposed method in Section 4, while Section 5 summarizes the main contributions of this paper. The symbols used in this paper are listed in Table A1 in Appendix A.

Global Sensitivity Analysis and Sobol's Indices
In this section, the traditional variance based GSA is briefly reviewed [16,19]. Let us consider a model y = g(x) with n-dimensional independent identical distribution random input x = (x 1 , . . . , x n ) ∈Ω n and a scalar output y. Denote f x (x) as the joint probability density function (PDF), then we have The Sobol's decomposition of g(x) can be expressed as: where: and so on x~[ i] in Equation (2) means the vector without the element of i, that is, x~[ i] = (x 1 , . . . , x i−1 , x i+1 , . . . , x n ). Furthermore, the expansion items above are orthogonal to each other, that is: And the variance of output y = g(x) is defined as: Sobol's proposed a variance-based representation to measure the importance of input variable. With Equations (1), (3) and (4), the Sobol's indices are given as: where V g i 1 ,··· ,i s x i 1 , · · · , x i s is called as the partial variance, which is distinguished from V(g(x)) above. Generally, S i is the first-order sensitivity index, S ij is the second-order sensitivity index, etc. There are 2 n − 1 Sobol's sensitivity indices, and they are subjected to the following formula: Moreover, since S i represents the single contribution to the output uncertainty only due to an input x i , it is also called as the main effect. And the total effect index S T i has been defined to measure the total contribution of an input x i . It could be calculated as follows: where S~[ i] is the sum of all S i 1 ,··· ,i s without subscribe index i. Traditionally, the Sobol's indices could be solved by MCS [16], which may be computationally demanding. Sudret proposed a PCE-based method to mitigate this problem with much less computational costs [18].
Moreover, φ k is the kth one-dimensional orthogonal polynomials from the Askey scheme.
ψ α ξ i 1 , · · · , ξ i s is the multi-dimensional PCE basis, which is constructed by tensor products of the corresponding one-dimensional polynomials, namely, the highest degree of orthogonal polynomials in Equation (8) is a predefined PCE degree p. Therefore, we have |α| ≤ p. It should be noted that p could be determined by considering the balance of accuracy and efficiency. Specifically, PCEs with small degrees p and p + 1 could be obtained respectively, which is followed by comparing the difference (usually the first two order moments). If the difference is smaller than a pre-defined threshold, the PCE with p + 1 degree is adopted. Otherwise, the PCE with p + 2 degree is used to repeat this procedure. Generally, Equation (8) can be rewritten as: where coefficient c j and expansion base ψ j (ξ) are corresponding to Equation (8) sequentially, N = n + p p = (n+p)! n!p! is the number of PCE items. Thus, the main effect index S i and the total effect index S T i could be solved directly with PCE coefficients as: For more details see Appendix B.

Global Reliability Sensitivity Analysis
If the given model g(x) above is the limit-state function, the input x falls in the failure region if g(x) returns a negative value. Denote P f as the failure probability, then we have: where is the indicator function of the failure domain, and E(I F ) is the mean of indicator function. Similarly, the CFP is defined as follows: where x I is the vector of conditional input variable with subscript I = {i 1 , . . . , is the conditional indicator function of the failure domain of the conditional limit-state function g(x|x I ) and it can be obtained by assuming x I at a fixed realization value while other inputs are free. Cui and Lü put forward the moment-based importance measure δ P I of the basic variable on the distribution as [14]: where δ P I eflects the effect of conditional variable x I on failure probability when x I takes its value according to its PDF. Furthermore, Cui proved that: In this way, δ P I is equal to the Sobol's index when I F is taken as the model output. Therefore, Cui proposed that the main effect index and the total effect index of x i can be defined as: where V(I F ), V(E(I F |x i )) and V(E(I F |x~[ i] )) are expressed as follows: The MCS method and associated improved Importance Sampling (IS) method for solving S (R) i and S (R) T i are presented in [27], where the IS method uses the advanced first-order second-moment method [28] to help construct importance sampling density functions within MCS framework. However, they are both time-consuming due to the sampling principle, and the computational burden becomes worse when given model are complex. In this way, a new algorithm is proposed to mitigate this problem.

PCE-Based GRSA Method
Take the limit-state function in Section 2.2 into consideration. As mentioned above, the model output for global reliability sensitivity indices calculation is I F , and it is difficult to compute I F related value except for MCS. This is due to the fact that I F only maps every single model output to two discrete real numbers 0 or 1. Either analytical methods or surrogate methods are efficient when they are applied to continuous model output rather than the discrete case. However, E(I F |x i ) and E(I F |x~[ i] ) could be replaced by two CFP values P f (x i ) and P f (x~[ i] ) respectively, which could be treated as continuous random variables. In this way, the key point of global reliability sensitivity index calculation is to compute variance values V(P f (x i )) and V(P f (x~[ i] )). Equations (17) and (18) are replaced by Equations (22) and (23): Since the PCE technology provides an efficient surrogate modeling tool to determinate the mean as well as the variance of the model output, a PCE-based method is developed for this implement. However, several probability values required during surrogate modeling procedure may cause computational demanding. So, the analytical ME method is integrated with PCE to solve this problem. Thus, the surrogate model of CFP could be developed, which is referred as R-PCE in this work.
The PCE-based GRSA method involves a 2-layer PCE surrogate model construction process. The first layer PCE is to generate traditional surrogate model approximating the limit-state function. It is used to obtain mean and variance when given conditional variables are fixed in specific values. And the second layer is to construct R-PCEs for different CFPs. To obtain probability response during the surrogate modeling of R-PCE, the ME method is taken to analytically compute these probability values with the first layer PCE. Besides, since each CFP is dependent on a different conditional variable, each global reliability sensitivity index has its own R-PCE formula. The algorithm framework is shown in Figure 1, and the whole procedure consists of four phases.
The first phase is to prepare inputs or data for GRSA.

C-PCEn
Index n

Index 1
Analytic solution of global reliability indices Index 2 Index n . . .

C-PCE2
Index 2 Analytic solution of CFP by Maximum entropy method The 2 nd layer PCE construction (R-PCE) The mulitiplication of C-PCE In the second phase, the PCE technology is an efficient and accurate tool to approximate the given model with finite samples. Thus, the original model is replaced by the 1st-layer PCE in the following procedures. Besides, the first and second moment of the model output are obtained directly without traditional Mote Carlo sampling.
In the third phase, the PCE technology is developed to make the surrogate model for CFP. Unlike the surrogate model in the second phase, quite a few failure probabilities of given model require to be calculated. Unlike the model output which may be obtained by simulating directly, it is usually difficult to obtain these statistical values. Fortunately, this problem could be overcome by the ME method. According to ME, the probability distribution that best represents the current state of knowledge, which is usually presented in terms of the first n-order moments, is the one having the largest entropy [24]. Under this principle, the probability distribution of the model output could be analytically obtained with moment information based on the 1st-layer PCE. The details of this procedure are presented as follows.
First of all, for each global reliability sensitivity index, a set of conditional variable samples is generated by the probabilistic collocation method (PCM). Second, due to the regularity of orthogonal polynomials, the conditional variable is removed from the 1st-layer PCE. Also, the specific formula of the 1st-layer PCE with other variables could remain the form of linear combination of orthogonal polynomials, so does the statistical characteristic. For simplicity, we call it conditional PCE (C-PCE). In the second phase, the PCE technology is an efficient and accurate tool to approximate the given model with finite samples. Thus, the original model is replaced by the 1st-layer PCE in the following procedures. Besides, the first and second moment of the model output are obtained directly without traditional Mote Carlo sampling.
In the third phase, the PCE technology is developed to make the surrogate model for CFP. Unlike the surrogate model in the second phase, quite a few failure probabilities of given model require to be calculated. Unlike the model output which may be obtained by simulating directly, it is usually difficult to obtain these statistical values. Fortunately, this problem could be overcome by the ME method. According to ME, the probability distribution that best represents the current state of knowledge, which is usually presented in terms of the first n-order moments, is the one having the largest entropy [24]. Under this principle, the probability distribution of the model output could be analytically obtained with moment information based on the 1st-layer PCE. The details of this procedure are presented as follows.
First of all, for each global reliability sensitivity index, a set of conditional variable samples is generated by the probabilistic collocation method (PCM). Second, due to the regularity of orthogonal polynomials, the conditional variable is removed from the 1st-layer PCE. Also, the specific formula of the 1st-layer PCE with other variables could remain the form of linear combination of orthogonal polynomials, so does the statistical characteristic. For simplicity, we call it conditional PCE (C-PCE). Thus, the associated moment evaluations are also available. Third, the ME method is taken to solve failure probabilities with regards to specific conditional variable. Fourth, a set of conditional variable samples and associated conditional failure probabilities are prepared for following 2nd-layer PCE (R-PCE) construction, and this procedure is similar to PCE in the second phase. In this way, R-PCE corresponding to specific conditional variable is generated.
Finally, with the coefficients of R-PCE, we could obtain variance of each CFP and calculate global reliability sensitivity index. During the procedure, the number of sampling of both the 1st-layer PCE and the 2nd-layer PCE of proposed method is limited according to PCE sample principle. This makes proposed method efficient. The following section will provide details of how to construct the proposed 2-layer PCE.

Conditional Polynomial Chaos Expansion and Moments Calculation
Without loss of generality, PCE formula of given limit-state function is the same as Equation (8).
If the value of ξ I is given as the conditional variable, Equation (8) could be simplified as: where |I| is the number of elements in subscript I, {i 1 , . . . , i s }

of 24
nditional variable is removed from the 1st-layer PCE. Also, the specific formula with other variables could remain the form of linear combination of orthogonal s the statistical characteristic. For simplicity, we call it conditional PCE (C-PCE). moment evaluations are also available. Third, the ME method is taken to solve with regards to specific conditional variable. Fourth, a set of conditional variable ted conditional failure probabilities are prepared for following 2nd-layer PCE (Rand this procedure is similar to PCE in the second phase. In this way, R-PCE ecific conditional variable is generated. e coefficients of R-PCE, we could obtain variance of each CFP and calculate global index. During the procedure, the number of sampling of both the 1st-layer PCE E of proposed method is limited according to PCE sample principle. This makes fficient. The following section will provide details of how to construct the E.

omial Chaos Expansion and Moments Calculation
generality, PCE formula of given limit-state function is the same as Equation (8). If en as the conditional variable, Equation (8) could be simplified as:    (24) is called as Cand variance could be obtained by the following equations: expression, C-PCE could be rewritten as: I represents that ∀i k ∈ {i 1 , · · · , i s }, i k / ∈ I. The coefficients are subject to: and: The total number of coefficients is N = n + p − |I| p = (n+p−|I|)! (n−|I|)!p! . Equation (24) is called as C-PCE. Thus, the mean and variance could be obtained by the following equations: To simplify the expression, C-PCE could be rewritten as: where c j and ψ α (ξ|ξ I ) are corresponding to Equation (28) sequentially, notation "ξ|ξ I " in ψ α (ξ|ξ I ) means that expansion base does not include conditional variable ξ I . The mean and variance of C-PCE are given as Equations (30) and (31) respectively: Specifically, for main effect indices, since we have I = {i}, Equation (29) can be further rewritten as: and for total effect indices, since we have I = {1, . . . , n}/{i}, Equation (29) can be further rewritten as: The C-PCE could provide the first two order moments accurately. However, there are some difficulties for high order moment calculation. Consider the fact that the C-PCE consists of orthogonal bases, the products of two C-PCEs could also be expanded as a linear combination of orthogonal bases [29]. Therefore, the orthogonal characteristic of PCE multiplication remains. The generalization form of orthogonal bases are presented in [29].
Take the C-PCE with Hermite polynomial bases for example, suppose two C-PCEs u = ∑ |α|≤p α c α ψ α (ξ|ξ I ) and v = ∑ |β|≤p β c β ψ β (ξ|ξ I ). If E |uv| 2 < ∞, then the product of u and v has the PCE formula as: where: The subscripts α, β, r, and θ in Equations (34) and (35) are tuples associated with C-PCE terms, e.g., α = (α 1 , α 2 , . . . , α n ). We say β ≤ θ if β i ≤ θ i for all i = 1, 2, . . . , n. The operation of these subscripts, such as + or −, is also defined as component-wise. Especially, the factorial of tuples is defined like Then, the mean of uv could be obtained as: Based upon the preparation above, we could calculate the high-order moments of PCE by replacing u and v with two PCEs respectively. That is, u and v in Equation (34) could be replaced by g PCE (ξ|ξ I ), g 2 PCE (ξ|ξ I ), g 3 PCE (ξ|ξ I ), . . . , respectively as Equation (37): where Entropy 2018, 20, x 10 of 24 I  I  I  I   I  I  I  I   I  I  I where ⌊•⌋ is rounded down and ⌈•⌉ is rounded up. Thus, the high order moments could be computed analytically.

The Maximum Entropy Method
Let fy(y ) be the unknown PDF of the output of limit-state function and assume H( fy(y)) to be the entropy of fy(y) . Then we have: where Ωy is the domain of y. According to the principle of ME, the failure probability distribution which best represents the current information, such as the moments, is the one with the largest entropy [24]. Thus, the ME method for estimating fy(y) can be formulated as: where NME is the number of moments, and μk is the kth moment, i.e., μ1 = E(y), μ2 = E(y 2 ), etc. It can be proved that the optimal solution for the PDF formula is given as: Specifically, if the first two moments of the output are available, we have: V y E y E y and the failure probability can be calculated as: In the proposed method, if taking Equations Error! Reference source not found. and Error! Reference source not found. into Equations Error! Reference source not found. and where ⌊•⌋ is rounded down and ⌈•⌉ is rounded up. Thus, the high order moments could be computed analytically.

The Maximum Entropy Method
Let fy(y ) be the unknown PDF of the output of limit-state function and assume H( fy(y)) to be the entropy of fy(y) . Then we have: where Ωy is the domain of y. According to the principle of ME, the failure probability distribution which best represents the current information, such as the moments, is the one with the largest entropy [24]. Thus, the ME method for estimating fy(y) can be formulated as: where NME is the number of moments, and μk is the kth moment, i.e., μ1 = E(y), μ2 = E(y 2 ), etc. It can be proved that the optimal solution for the PDF formula is given as: Specifically, if the first two moments of the output are available, we have: V y E y E y and the failure probability can be calculated as: In the proposed method, if taking Equations Error! Reference source not found. and Error! Reference source not found. into Equations Error! Reference source not found. and is rounded up. Thus, the high order moments could be computed analytically.

The Maximum Entropy Method
Let f y (y) be the unknown PDF of the output of limit-state function and assume H(f y (y)) to be the entropy of f y (y). Then we have: where Ω y is the domain of y. According to the principle of ME, the failure probability distribution which best represents the current information, such as the moments, is the one with the largest entropy [24]. Thus, the ME method for estimating f y (y) can be formulated as: H f y (y) s.t. y∈Ω y f y (y)dy = 1 y∈Ω y y k f y (y)dy = µ k , k = 1, 2, · · · , N ME , where N ME is the number of moments, and µ k is the kth moment, i.e., µ 1 = E(y), µ 2 = E(y 2 ), etc. It can be proved that the optimal solution for the PDF formula is given as: Thus, λ i in Equation (40) should satisfy the following equation: Specifically, if the first two moments of the output are available, we have: where V(y) = E y 2 − E(y) 2 and the failure probability can be calculated as: In the proposed method, if taking Equations (27) and (28) into Equations (42) and (43), the failure probability of C-PCE could be obtained directly without additional sampling.
For more general cases, the first four moments are usually used for PDF calculation because they are corresponding to mean, variance, skewness, kurtosis respectively, which are clear in physics meaning [30]. The 3rd moment and the 4th moment could be obtained efficiently by Equation (37) instead of traditional sampling method. Besides, since the ME method may not have the closed form solution, some optimization techniques, such as Newton method or its improved versions [31,32], could be utilized to solve the failure probability.

Global Reliability Sensitivity Indices Calculation
Let ξ denote a set of conditional variable samples, where M R is the number of samples.
As described above, when the PCE for limit-state function is given, if taking each conditional variable sample into Equation (24), it returns the mean and variance of C-PCE. Then, the CFP is calculated by ME method. In this way, we could get a set of surrogate model relationship between . Similar to the traditional PCE procedure as Equation (A7) in Appendix B, we could solve the coefficients of R-PCE. Then, the result of R-PCE is: where c R,j is jth coefficient of R-PCE, and the subscript "R" means the coefficient is used for GRSA.
ξ I is independent standardized random variables corresponding to original variables x I , N R is the total number of coefficients, which is determined by: where p R is the predefined degree of R-PCE. Then, the CFP is: And the variance of CFP is: Substituting Equations (46) and (47) into Equations (22) and (23), we could calculate the main effect and the total effect of global reliability sensitivity indices.

Implementation
For a given limit-state function g = g(x), we assume here that the input random variables are independent. With preparation mentioned above, the complete algorithm description is introduced as follows: Step 1: Represent distributions of the input into the associated standardized random variable. If the distribution types are not unified, we could map different distributions into the normal distribution. The following computation will be based upon standardized random variables.
Step 2: Compute PCE for limit-state function by Equation (A7) in Appendix B, and it refers as g PCE (ξ).
Step 3: Based upon PCE above, compute R-PCE. First of all, generate a set of conditional variable samples. Then, integrate g PCE (ξ) and ME method to obtain the CFP value P f (ξ|ξ I ), where I = {i} PCE (ξ|ξ I ). With the moments given by C-PCE, we could calculate the P f ξ (j) I with ME method analytically. Third, the R-PCE comes out as Equation (44).
The step-by-step procedure of algorithm for the main effect indices and the total effect indices by the proposed 2-layer PCE method are presented in Algorithms A1 and A2 in Appendix C respectively. The total evaluations of the limit-state function is 2(n+p)! n!p! + 1, where n and p are input variable dimension and the PCE degree respectively. If taking times of PCE evaluation in R-PCE construction process into consideration, the computation costs are where p R is the degree of R-PCE. This computation burden is usually 1-2 times less than IS procedure and 2-3 times less than MCS procedure, which could be seen in the following test examples. And the efficiency of the proposed method could be further improved by adopting sparse algorithms [6,17,33]. However, such work is not the scope of this paper.

Test Examples
Implementation of the proposed method for global reliability sensitivity indices calculation is illustrated by two test examples. The first example is a numerical example, while the second one is a structural analysis of a bridge crane. The results obtained by the proposed method are compared with solutions obtained by the full-scale MCS and the IS method from the perspective of accuracy and efficiency respectively.

A Numerical Example
Consider a limit-state function: where the input variables follow the normal distribution or the lognormal distribution, i.e., lnx 1~N (0.5,0.1), lnx 2~N  , when the PCE degree for the limit-state function is 2, and the R-PCE degree for CFP is 7. The results are compared with the results obtained by the full-scale MCS method and the IS method in Table 1. Moreover, the total times of limit-state function evaluation by different methods are also shown in the last column. It could be seen that the relative error of both methods are less than 3%, and about half of the results calculated by the proposed method is better than those obtained by the IS method. Moreover, the computational cost of proposed method is reduced sufficiently. Also, the main effect indices and the total effect indices of GRSA are schematically shown in Figure 2. We could see that both the main effect indices S (R) i and the total effect indices S (R) T i produce the same importance ranking: x 1 > x 3 > x 2 , and there are strong interaction effects among these input variables since the total effects are much larger than the main effects. This suggests that the risk may be reduced by decreasing the uncertainty of x 1 .
Entropy 2018, 20, x 13 of 23 these input variables since the total effects are much larger than the main effects. This suggests that the risk may be reduced by decreasing the uncertainty of x1.

An Engineering Example
The second example is a bridge crane with respect to its rigidity design in Figure 3. The crane is made by universal beam with movable sheaves and hook. When the weight is in the middle of the beam, the maximal deflection  max of such a crane beam could be computed as: where q is the magnitude of a normally distributed load that represents the dead-weight of crane beam, E and I denote Young's modulus and inertia moment of universal beam respectively, L is the span length, and Fp is maximal weight including listed load, sheaves and hook. Two distribution types are taken for this example, namely, the normal distribution and the lognormal distribution. Besides, two different variances are also considered to test this example, which are labeled as case1 and case 2 respectively. The detail distribution parameters are given in Table 2.

An Engineering Example
The second example is a bridge crane with respect to its rigidity design in Figure 3. The crane is made by universal beam with movable sheaves and hook. When the weight is in the middle of the beam, the maximal deflection ω max of such a crane beam could be computed as: where q is the magnitude of a normally distributed load that represents the dead-weight of crane beam, E and I denote Young's modulus and inertia moment of universal beam respectively, L is the span length, and F p is maximal weight including listed load, sheaves and hook. Two distribution types are taken for this example, namely, the normal distribution and the lognormal distribution. Besides, two different variances are also considered to test this example, which are labeled as case1 and case 2 respectively. The detail distribution parameters are given in Table 2.
Entropy 2018, 20, x 13 of 23 these input variables since the total effects are much larger than the main effects. This suggests that the risk may be reduced by decreasing the uncertainty of x1.

An Engineering Example
The second example is a bridge crane with respect to its rigidity design in Figure 3. The crane is made by universal beam with movable sheaves and hook. When the weight is in the middle of the beam, the maximal deflection  max of such a crane beam could be computed as: where q is the magnitude of a normally distributed load that represents the dead-weight of crane beam, E and I denote Young's modulus and inertia moment of universal beam respectively, L is the span length, and Fp is maximal weight including listed load, sheaves and hook. Two distribution types are taken for this example, namely, the normal distribution and the lognormal distribution. Besides, two different variances are also considered to test this example, which are labeled as case1 and case 2 respectively. The detail distribution parameters are given in Table 2.    According to design requirement, the admissible ratio of deflection to span length is lower than 0.0018. Thus, the limit-state function is given as: The results obtained by the full-scale MCS method, the IS method and the proposed method are given in Tables 3 and 4. The full-scale MSC results with large sample size are seen as the accurate results. In the calculation procedure of the proposed method, the PCE degree for the limit-state function and the PCE degree for CFP are 5 and 7 respectively.
It is shown that both the IS method and the proposed method could approximate S (R) i and S (R) T i with much less computational burden, and the latter one presents better efficiency than the former one. In case 1, the average relative error of global reliability sensitivity indices obtained by proposed method is about 0.0196, while that obtained by the IS method is about 0.0283. In case 2, the average relative error of global reliability sensitivity indices obtained by proposed method is about 0.0206, while that obtained by the IS method is about 0.0295. The computational costs of our method are about two orders of magnitude less than the IS method, and four orders of magnitude less than the MCS method. Thus, the proposed method may be better than the other two methods.   is given as E > L > I > F p > q. Obviously, the influence of Young's modulus on failure probability is higher than other random variables, which is different from case 1. Moreover, strong interaction effects also exist in this case. However, the interaction effects make more contributions than those in case 1 even when input variances are smaller compared with those of case 1. Therefore, the importance ranking should be paid more attention during further optimization design.
Comparing the results of GRSA with those of GSA, it suggests that the importance of input variables on the failure probability is not equivalent to that on the model output. So, both global reliability sensitivity indices and global sensitivity indices may be considered in engineering practice when choosing importance variables. And the proposed method could provide a useful tool for such purpose.

Conclusions
Global reliability sensitivity analysis aims at measuring the importance statistical parameter with regards to the input variable. The variance of expectation of conditional indicator function in global reliability sensitivity indices is difficult to calculate. This paper develops a PCE-based GRSA method to solve the problem. The proposed method takes conditional failure probability to replace the expectation of conditional indicator function, and then utilizes a 2-layer PCE construction technology to obtain the variance of conditional failure probability. The efficiency of proposed method is illustrated by a numerical example and an engineering model example. Both main effect indices and total effect indices present good performance in accuracy and significantly less computational burden in limit-state function evaluations compared with the traditional full-scale MCS method and IS method. Through global reliability sensitivity analysis, we could determine which uncertainties of inputs variables should be controlled and pave the way for further design improvement or optimization.
Acknowledgments: This work was supported by the National Natural Science Foundation of China (NSFC: 61304218) and Beijing Natural Science Foundation (BNSF: 3153027). The authors would like to extend their sincere thanks to anonymous reviewers for their valuable comments.
Author Contributions: Jianbin Guo conceived the framework and structured the whole paper; Jianyu Zhao performed the experiments, wrote and revised the paper; Shengkui Zeng checked the results and revised the paper; Shaohua Du helped write programs.

Conflicts of Interest:
The authors declare no conflict of interest.  The global sensitivity indices of maximal deflection of cases 1 and 2 are also computed in Table 5 by traditional GSA method mentioned in Section 2.1, and they are plotted in Figures 4b and 5b respectively. It is shown that the results of global sensitivity indices and global reliability sensitivity indices may be different. First of all, although the importance ranking based on GSA of case 1 is the same as that based on GRSA, the importance ranking obtained by GSA of case 2 is L > E > I > Fp > q, which is different from that obtained by GRSA. Secondly, the total effect indices of both cases are similar to the main effect indices, which means that there are few interaction effects among input variables. This is also different from the results of GRSA.

Appendix A
Comparing the results of GRSA with those of GSA, it suggests that the importance of input variables on the failure probability is not equivalent to that on the model output. So, both global reliability sensitivity indices and global sensitivity indices may be considered in engineering practice when choosing importance variables. And the proposed method could provide a useful tool for such purpose.

Conclusions
Global reliability sensitivity analysis aims at measuring the importance statistical parameter with regards to the input variable. The variance of expectation of conditional indicator function in global reliability sensitivity indices is difficult to calculate. This paper develops a PCE-based GRSA method to solve the problem. The proposed method takes conditional failure probability to replace the expectation of conditional indicator function, and then utilizes a 2-layer PCE construction technology to obtain the variance of conditional failure probability. The efficiency of proposed method is illustrated by a numerical example and an engineering model example. Both main effect indices and total effect indices present good performance in accuracy and significantly less computational burden in limit-state function evaluations compared with the traditional full-scale MCS method and IS method. Through global reliability sensitivity analysis, we could determine which uncertainties of inputs variables should be controlled and pave the way for further design improvement or optimization.
Acknowledgments: This work was supported by the National Natural Science Foundation of China (NSFC: 61304218) and Beijing Natural Science Foundation (BNSF: 3153027). The authors would like to extend their sincere thanks to anonymous reviewers for their valuable comments.
Author Contributions: Jianbin Guo conceived the framework and structured the whole paper; Jianyu Zhao performed the experiments, wrote and revised the paper; Shengkui Zeng checked the results and revised the paper; Shaohua Du helped write programs.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Notation.

Appendix A
n the dimension of a variable x x = (x 1 , x 2 , · · · , x n ), the input of given model ξ ξ = (ξ 1 , ξ 2 , · · · , ξ n ), the independent standardized orthogonal random variable y the output of given model , the joint PDF of x = (x 1 , x 2 , · · · , x n ) f y( y) the PDF of y E(·) expectation V(·) variance S Ti the main effect of GSA S i the total effect of GSA α k ≤ p and α = (α 1 , α 2 , · · · , α n ), the base of PCE I F the indicator function P f the failure probability λ i the unknown parameter of ME formula N ME the number of unknown parameter of ME formula µ k the kth moment of y y ≈ g PCE (ξ) = N−1 ∑ j=0 c j ψ j (ξ), f or ξ = (ξ 1 , · · · , ξ n ), where coefficient c j and expansion base ψ j (ξ) are corresponding to Equation (A4) sequentially.
Let ξ i M i=1 denote a set of the random variable samples, and it could be determined by the probabilistic collocation method (PCM). Again Let g ξ i M i=1 denote the corresponding set of model output, where M is the number of samples. Denote c = (c 0 , · · · , c N−1 ) T , thus an approximationĉ is given by: where ψ = ψ 0 ξ (i) · · · ψ N ξ (i) denotes the coefficient matrix and y ξ 1:M = g ξ 1 , · · · , g ξ M T . Since ψ T ψ may be ill-conditioned, M is suggested to be selected as M = 2(N + 1) [18]. Due to the orthogonality of the basis, the mean value and the variance of y in Equation (A6) can be calculated as: It could be proved that the PCE-based main effect index S i and total effect index S T i are derived as [6]: Therefore, the PCE could obtain the expectation and variance of the model output efficiently. Moreover, compared with the traditional MCS method which requires a large number of samples [16], the PCE method could calculate Sobol's indices with quit a few computation costs. In this way, PCE is widely employed for GSA application.