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

Abstract

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.

1. 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].

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. Notation.

$n$	the dimension of a variable
$\mathit{x}$	$\mathit{x}=\left(x_1,x_2,\cdots ,x_n\right)$, the input of given model
$\mathit{\xi }$	$\mathit{\xi }=\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)$, the independent standardized orthogonal random variable
$y$	the output of given model
${\mathit{x}}_{~\left[i\right]}\left({\mathit{\xi }}_{~\left[i\right]}\right)$	the vector without the ith element
$f_x\left(x_i\right)$	the PDF of xi
$f_x\left(\mathit{x}\right)$	$f_x\left(\mathit{x}\right)={\displaystyle \prod _{i=1}^n{f}_x\left(x_i\right)}$, the joint PDF of $\mathit{x}=\left(x_1,x_2,\cdots ,x_n\right)$
$f_y\left(y\right)$	the PDF of y
$E(·)$	expectation
$V(·)$	variance
$S_{T_i}$	the main effect of GSA
$S_i$	the total effect of GSA
$S_i^{\left(R\right)}$	the main effect of GRSA
$S_{T_i}^{\left(R\right)}$	the total effect of GRSA
$p$	the degree of PCE
$p_R$	the degree of R-PCE
$c_j$	the coefficient of jth PCE item
${c^{\prime }}_j$	the coefficient of jth CPCE item
$c_{R,j}$	the coefficient of jth R-PCE item
$N$	the number of PCE items
$N^{\prime }$	the number of CPCE items
$N_R$	the number of R-PCE items
${\varphi }_k$	the kth one-dimensional orthogonal polynomial
${\psi }_{\mathit{\alpha }}\left(\mathit{\xi }\right)$	${\psi }_{\mathit{\alpha }}\left(\mathit{\xi }\right)={\displaystyle \prod _{k=1}^n{\varphi }_{{\alpha }_k}\left({\xi }_k\right)}$ where ${\displaystyle \sum _{k=1}^j{\alpha }_k\le p}$ and $\mathit{\alpha }=\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n\right)$, the base of PCE
$I_F$	the indicator function
$P_f$	the failure probability
${\lambda }_i$	the unknown parameter of ME formula
$N_{ME}$	the number of unknown parameter of ME formula
${\mu }_k$	the kth moment of $y$

in Appendix A.

2. Related Work

2.1. 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₁, …, x_n) ∈Ωⁿ and a scalar output y. Denote f_x(x) as the joint probability density function (PDF), then we have $f_x\left(\mathit{x}\right)={\displaystyle \prod _{i=1}^n{f}_x\left(x_i\right)}$. In addition, we assume that $\mathit{x}\in {\mathit{\Omega }}^n\mapsto g\left(\mathit{x}\right)\in L^2\left({\mathit{\Omega }}^n,f_x\left(\mathit{x}\right)\right)$.

The Sobol’s decomposition of g(x) can be expressed as:

$$g\left(\mathit{x}\right)=g_0+{\displaystyle \sum _{i=1}^n{g}_i\left(x_i\right)+{\displaystyle \sum _{1\le i<j\le n}{g}_{ij}\left(x_i,x_j\right)+\cdots +g_{1,2,\cdots ,n}\left(x_1,\cdots ,x_n\right)}},$$

(1)

where:

$$\begin{array}{l}{g}_0={\displaystyle {\int }_{{\mathit{\Omega }}^n}g\left(\mathit{x}\right)f_x\left(\mathit{x}\right)d\mathit{x}},\\ g_i\left(x_i\right)={\displaystyle {\int }_{{\mathit{\Omega }}^{n-1}}g\left(\mathit{x}\right)f_x\left({\mathit{x}}_{~[i]}\right)d{\mathit{x}}_{~[i]}-g_0},\\ g_{ij}\left(x_i,x_j\right)={\displaystyle {\int }_{{\mathit{\Omega }}^{n-2}}g\left(\mathit{x}\right)f_x\left({\mathit{x}}_{~[i,j]}\right)d{\mathit{x}}_{~[i,j]}-g_i\left(x_i\right)}-g_j\left(x_j\right)-g_0,\end{array}$$

(2)

and so on x_~[i] in Equation (2) means the vector without the element of i, that is, x_~[i] = (x₁, …, x_i−1, x_i+1, …, x_n).

Furthermore, the expansion items above are orthogonal to each other, that is:

$${\displaystyle {\int }_{{\mathit{\Omega }}^n}{g}_{1,2,\cdots ,i_s}\left(x_1,x_2,\cdots ,x_{i_s}\right)}{g}_{1,2,\cdots ,j_t}\left(x_1,x_2,\cdots ,x_{j_t}\right)f_x\left(\mathit{x}\right)d\mathit{x}=0\forall \left(x_1,x_2,\cdots ,x_{i_s}\right)\ne \left(x_1,x_2,\cdots ,x_{j_t}\right).$$

(3)

And the variance of output y = g(x) is defined as:

$$V\left(g\left(\mathit{x}\right)\right)={\displaystyle {\int }_{{\mathit{\Omega }}^n}{\left(g\left(\mathit{x}\right)-g_0\right)}^2}{f}_x\left(\mathit{x}\right)d\mathit{x}.$$

(4)

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:

$$S_{i_1,\cdots ,i_s}=\frac{V\left(g_{i_1,\cdots ,i_s}\left(x_{i_1},\cdots ,x_{i_s}\right)\right)}{V\left(g\left(\mathit{x}\right)\right)}=\frac{{\displaystyle {\int }_{{\mathit{\Omega }}^s}{g}_{i_1,\cdots ,i_s}^2\left(x_{i_1},\cdots ,x_{i_s}\right)}{f}_x\left(x_{i_1},\cdots ,x_{i_s}\right)d{x}_{i_1}\cdots x_{i_s}}{{\displaystyle {\int }_{{\mathit{\Omega }}^n}{\left(g\left(\mathit{x}\right)-g_0\right)}^2}{f}_x\left(\mathit{x}\right)d\mathit{x}},$$

(5)

where $V\left(g_{i_1,\cdots ,i_s}\left(x_{i_1},\cdots ,x_{i_s}\right)\right)$ 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ⁿ − 1 Sobol’s sensitivity indices, and they are subjected to the following formula:

$${\displaystyle \sum _{i=1}^n{S}_i}+{\displaystyle \sum _{1\le i<j\le n}{S}_{ij}}+\cdots S_{1,2,\cdots ,n}=1.$$

(6)

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:

$$S_{T_i}={\displaystyle \sum _{L_i}{S}_{i_1,\cdots ,i_s}}=1-S_{~[i]},\forall L_i=\left\{\left(i_1,\cdots ,i_s\right):\exists k,1\le k\le s,i_k=i\right\},$$

(7)

where S_~[i] is the sum of all $S_{i_1,\cdots ,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].

PCE, which was originally introduced by Wiener [25], employs the Hermite polynomials in the random space to approximate the Gaussian stochastic processes [19,26]. In [20], Xiu and Karniadakis developed the PCE under Wiener-Askey scheme that could be applied to non-Gaussian scenarios. It can uniformly approximate any random process with finite second-order moments. Specifically, the PCE approximation of model y = g(x) is:

$$\begin{array}{ll}y\approx g_{PCE}\left(\mathit{\xi }\right)=& c_0+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_i}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_i\right)}}+{\displaystyle \sum _{1\le i_1<i_2\le n}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{i_1{i}_2}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},{\xi }_{i_2}\right)}}+\cdots +\\ & {\displaystyle \sum _{1\le i_1<\cdots <i_s\le n}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{i_1\cdots i_s}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)}}+\cdots +{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{1,2,\cdots ,n}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_1,\cdots ,{\xi }_n\right)},\end{array}$$

(8)

where ξ = (ξ₁, ξ₂, …, ξ_n) is independent standardized orthogonal random variables associated with input x = (x₁, …, x_n), the subscript α is a tuple defined as α = (α₁, …, α_n), and ${\phi }_{i_1,\cdots ,i_s}$ is defined as a realization of α tuple in which only the indices {i₁, …, i_s} are nonzero:

$${\phi }_{i_1,\cdots ,i_s}=\left\{\mathit{\alpha }:\begin{array}{c}{\alpha }_k>0\forall k=1,2,\cdots ,n.k\in \left(i_1,\cdots ,i_s\right)\\ {\alpha }_k=0\forall k=1,2,\cdots ,n.k\notin \left(i_1,\cdots ,i_s\right)\end{array}\right\}.$$

(9)

Moreover, ${\varphi }_k$ is the kth one-dimensional orthogonal polynomials from the Askey scheme. ${\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)$ is the multi-dimensional PCE basis, which is constructed by tensor products of the corresponding one-dimensional polynomials, namely, ${\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)={\displaystyle \prod _{k=1}^{i_s}{\varphi }_{{\alpha }_k}\left({\xi }_k\right)}$. Besides, 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:

$$y\approx g_{PCE}\left(\mathit{\xi }\right)={\displaystyle \sum _{j=0}^{N-1}{c}_j}{\psi }_j\left(\mathit{\xi }\right),for\mathit{\xi }=\left({\xi }_1,\cdots ,{\xi }_n\right),$$

(10)

where coefficient c_j and expansion base ψ_j(ξ) are corresponding to Equation (8) sequentially, $N=\left(\begin{array}{c}n+p\\ p\end{array}\right)=\frac{\left(n+p\right)!}{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:

$$S_i=\frac{{\displaystyle {\sum }_{\mathit{\alpha }\in {\phi }_i}{c}_{\mathit{\alpha }}^{2}E\left[{\psi }_{\mathit{\alpha }}^2\right]}}{{\displaystyle \sum _{j=1}^{N-1}{c}_j^2}E\left({\psi }_j^2\left(\mathit{\xi }\right)\right)},$$

(11)

$$S_{T_i}={\displaystyle \sum _{\left(i_1,\cdots ,i,\cdots ,i_s\right)\in {\phi }_{i_1,\cdots ,i_s}}{S}_{i_1,\cdots ,i,\cdots i_s}}.$$

(12)

For more details see Appendix B.

2.2. 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:

$$\begin{array}{ll}{P}_f& =P\left(g\left(\mathit{x}\right)<0\right)\\ & ={\displaystyle \int \cdots {\displaystyle {\int }_{g\left(\mathit{x}\text{<}\right)}{f}_x\left(x_1,x_2,\cdots ,x_n\right)}}d{x}_{1}d{x}_2\cdots d{x}_n\\ & ={\displaystyle \int \cdots {\displaystyle {\int }_{{\mathit{\Omega }}^n}{I}_F{f}_x\left(x_1,x_2,\cdots ,x_n\right)}}d{x}_{1}d{x}_2\cdots d{x}_n\\ & =E\left[I_F\right],\end{array}$$

(13)

where $I_F=\{\begin{array}{cc}1,& g\left(\mathit{x}\right)\le 0,\\ 0,& g\left(\mathit{x}\right)>0.\end{array}$ 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:

$$P_{f|{\mathit{x}}_{\mathit{I}}}=E\left[I_F|{\mathit{x}}_{\mathit{I}}\right],$$

(14)

where x_I is the vector of conditional input variable with subscript I = {i₁, …, i_m} (1 ≤ i₁ ≤…≤ i_m ≤ n), $I_F|{\mathit{x}}_{\mathit{I}}=\{\begin{array}{cc}1,& g\left(\mathit{x}|{\mathit{x}}_{\mathit{I}}\right)\le 0,\\ 0,& g\left(\mathit{x}|{\mathit{x}}_{\mathit{I}}\right)>0.\end{array}$ 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 ${\delta }_{\mathit{I}}^P$ of the basic variable on the distribution as [14]:

$${\delta }_{\mathit{I}}^P={\displaystyle {\int }_{-\infty }^{+\infty }{\left(P_f-P_{f|{\mathit{x}}_{\mathit{I}}}\right)}^2}{f}_{\mathit{x}}\left({\mathit{x}}_{\mathit{I}}\right)d{\mathit{x}}_{\mathit{I}}=E{\left[P_f-P_{f|{\mathit{x}}_{\mathit{I}}}\right]}^2,$$

(15)

where ${\delta }_{\mathit{I}}^P$ 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:

$${\delta }_{\mathit{I}}^P=V\left[E\left(I_F|{\mathit{x}}_{\mathit{I}}\right)\right].$$

(16)

In this way, ${\delta }_{\mathit{I}}^P$ 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:

$$S_i^{\left(R\right)}=\frac{V\left(E\left(I_F|x_i\right)\right)}{V\left(I_F\right)},$$

(17)

$$S_{T_i}^{\left(R\right)}=1-\frac{V\left(E\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)\right)}{V\left(I_F\right)},$$

(18)

where V(I_F), V(E(I_F|x_i)) and V(E(I_F|x_~[i])) are expressed as follows:

$$V\left(I_F\right)=E\left(I_F^2\right)-E{\left(I_F\right)}^2=E\left(I_F\right)-E{\left(I_F\right)}^2=P_f-P_f^2,$$

(19)

$$V\left(E\left(I_F|x_i\right)\right)=E\left(E{\left(I_F|x_i\right)}^2\right)-E{\left(E\left(I_F|x_i\right)\right)}^2=E\left(E{\left(I_F|x_i\right)}^2\right)-P_f^2,$$

(20)

$$V\left(E\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)\right)=E\left(E{\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)}^2\right)-E{\left(E\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)\right)}^2=E\left(E{\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)}^2\right)-P_f^2.$$

(21)

The MCS method and associated improved Importance Sampling (IS) method for solving $S_i^{\left(R\right)}$ and $S_{T_i}^{\left(R\right)}$ 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.

3. R-PCE Framework for GRSA

3.1. 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):

$$S_i^{\left(R\right)}=\frac{V\left(E\left(I_F|{\mathit{x}}_i\right)\right)}{V\left(I_F\right)}=\frac{E\left(P_f{\left(x_i\right)}^2\right)-P_f^2}{P_f-P_f^2},$$

(22)

$$S_{T_i}^{\left(R\right)}=1-\frac{V\left(E\left(I_F|{\mathit{x}}_{~\left[i\right]}\right)\right)}{V\left(I_F\right)}=1-\frac{E\left(P_f{\left(x_{~\left[i\right]}\right)}^2\right)-P_f^2}{P_f-P_f^2}.$$

(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.

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.

3.2. 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:

$$\begin{array}{ll}{g}_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)=& {c^{\prime }}_0+{\displaystyle \sum _{i=1}^{n-\left|I\right|}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{i\notin I}}{c^{\prime }}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_i\right)}}\\ & +{\displaystyle \sum _{1\le i_1<\cdots <i_s\le n-\left|I\right|}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{\left\{i_1\cdots i_s\right\}\notin I}}{c^{\prime }}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)}}+\cdots \\ & +{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{\left\{1,2,\cdots ,n-\left|I\right|\right\}\notin I}}{c^{\prime }}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_1,\cdots ,{\xi }_{n-\left|I\right|}\right)},\end{array}$$

(24)

where |I| is the number of elements in subscript I, {i₁, …, i_s}∉I represents that $\forall i_k\in \left\{i_1,\cdots ,i_s\right\},i_k\notin \mathit{I}$. The coefficients are subject to:

$${c^{\prime }}_0=c_0+{\displaystyle \sum _{\forall I^{\prime }\in I,I^{\prime }\notin \varnothing ,\mathit{\alpha }\in {\phi }_{I^{\prime }}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\mathit{\xi }}_{I^{\prime }}\right)}=c_0+{\displaystyle \sum _{\forall I^{\prime }=\left\{i_1,\cdots ,i_k\right\}\in I,I^{\prime }\notin \varnothing ,\mathit{\alpha }\in {\phi }_{I^{\prime }}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_k}\right)}$$

(25)

and:

$${c^{\prime }}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)={\displaystyle \sum _{\forall I^{\prime }\in I,I^{\prime }\notin \varnothing ,\mathit{\alpha }\in {\phi }_{i_1\cdots i_s,I^{\prime }}}{c}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left({\xi }_{i_1},\cdots ,{\xi }_{i_s},{\mathit{\xi }}_{I^{\prime }}\right)}.$$

(26)

The total number of coefficients is $N^{\prime }=\left(\begin{array}{c}n+p-\left|I\right|\\ p\end{array}\right)=\frac{\left(n+p-\left|I\right|\right)!}{\left(n-\left|I\right|\right)!p!}$. Equation (24) is called as C-PCE. Thus, the mean and variance could be obtained by the following equations:

$$E\left(g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)={c^{\prime }}_0,$$

(27)

$$\begin{array}{ll}V\left(g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)=& {\displaystyle \sum _{i=1}^{n-\left|I\right|}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{i\notin I}}{\left({c^{\prime }}_{\mathit{\alpha }}\right)}^{2}E\left({\psi }_{\mathit{\alpha }}^2\left({\xi }_i\right)\right)}}+\cdots \\ & +{\displaystyle \sum _{1\le i_1<\cdots <i_s\le n-\left|I\right|}{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{\left\{i_1\cdots i_s\right\}\notin I}}{\left({c^{\prime }}_{\mathit{\alpha }}\right)}^{2}E\left({\psi }_{\mathit{\alpha }}^2\left({\xi }_{i_1},\cdots ,{\xi }_{i_s}\right)\right)}}+\cdots \\ & +{\displaystyle \sum _{\mathit{\alpha }\in {\phi }_{\left\{1,2,\cdots ,n-\left|I\right|\right\}\notin I}}{\left({c^{\prime }}_{\mathit{\alpha }}\right)}^{2}E\left({\psi }_{\mathit{\alpha }}^2\left({\xi }_1,\cdots ,{\xi }_{n-\left|I\right|}\right)\right)}.\end{array}$$

(28)

To simplify the expression, C-PCE could be rewritten as:

$$g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)={\displaystyle \sum _{j=0}^{N^{\prime }-1}{c^{\prime }}_j{\psi }_{\mathit{\alpha }}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)},$$

(29)

where ${c^{\prime }}_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:

$$E\left(g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)={c^{\prime }}_0,$$

(30)

$$V\left(g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)={\displaystyle \sum _{j=1}^{N^{\prime }-1}{c^{\prime }}_j^{2}E\left({\psi }_j^2\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)}.$$

(31)

Specifically, for main effect indices, since we have I = {i}, Equation (29) can be further rewritten as:

$$g_{PCE}\left({\mathit{\xi }}_{~[i]}\right)={\displaystyle \sum _{j=0}^{N^{\prime }-1}{c^{\prime }}_j{\psi }_{\mathit{\alpha }}\left({\mathit{\xi }}_{~[i]}\right)}.$$

(32)

and for total effect indices, since we have I = {1, …, n}/{i}, Equation (29) can be further rewritten as:

$$g_{PCE}\left({\xi }_i\right)={\displaystyle \sum _{j=0}^{N^{\prime }-1}{c^{\prime }}_j{\varphi }_{\mathit{\alpha }}\left({\xi }_i\right)}.$$

(33)

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={\displaystyle \sum _{\left|\mathit{\alpha }\right|\le p_{\alpha }}{c^{\prime }}_{\mathit{\alpha }}{\psi }_{\mathit{\alpha }}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)}$ and $v={\displaystyle \sum _{\left|\mathit{\beta }\right|\le p_{\beta }}{c^{\prime }}_{\mathit{\beta }}{\psi }_{\mathit{\beta }}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)}$. If $E\left({\left|uv\right|}^2\right)<\infty$, then the product of u and v has the PCE formula as:

$$uv={\displaystyle \sum _{\left|\mathit{\theta }\right|\le p_{\alpha }+p_{\beta }}{\displaystyle \sum _{0\le \mathit{\beta }\le \mathit{\theta }}{\displaystyle \sum _{\begin{array}{l}\left|\mathit{\theta }-\mathit{\beta }+\mathit{r}\right|\le p_{\alpha },\\ \left|\mathit{\beta }+\mathit{r}\right|\le p_{\beta }\end{array}}C\left(\mathit{\theta },\mathit{\beta },\mathit{r}\right)u_{\mathit{\theta }-\mathit{\beta }+\mathit{r}}{v}_{\mathit{\beta }+\mathit{r}}{\psi }_{\mathit{\theta }}\left(\mathit{\xi }\right)}}},$$

(34)

where:

$$C\left(\mathit{\theta },\mathit{\beta },\mathit{r}\right)={\left[\left(\begin{array}{c}\mathit{\theta }-\mathit{\beta }+\mathit{r}\\ \mathit{r}\end{array}\right)\left(\begin{array}{c}\mathit{\beta }+\mathit{r}\\ \mathit{r}\end{array}\right)\left(\begin{array}{c}\mathit{\theta }\\ \mathit{\theta }-\mathit{\beta }\end{array}\right)\right]}^{\frac{1}{2}}.$$

(35)

The subscripts α, β, r, and θ in Equations (34) and (35) are tuples associated with C-PCE terms, e.g., α = (α₁, α₂, …, α_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 α! = ∏_i α_i!.

Then, the mean of uv could be obtained as:

$$E\left(uv\right)={\displaystyle \sum _{\left|\mathit{r}\right|\le \min (p_{\alpha },p_{\beta })}C\left(\mathit{\theta }=\mathit{0},\mathit{\beta }=\mathit{0},\mathit{r}\right)u_{\mathit{r}}{v}_{\mathit{r}}}={\displaystyle \sum _{\left|\mathit{r}\right|\le \min (p_{\alpha },p_{\beta })}{u}_{\mathit{r}}{v}_{\mathit{r}}}.$$

(36)

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}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)$, $g_{PCE}^2\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)$, $g_{PCE}^3\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)$, …, respectively as Equation (37):

$$\begin{array}{l}E\left(y^2\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)\approx E\left(g_{PCE}^2\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)=E\left(g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\cdot g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)\\ E\left(y^3\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)\approx E\left(g_{PCE}^3\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)=E\left(y_p^2\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\cdot g_{PCE}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)\\ ⋮\\ E\left(y^k\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)\approx E\left(g_{PCE}^k\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right)=E\left(g_{PCE}^{\lfloor \frac{k}{2}\rfloor }\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\cdot g_{PCE}^{\lceil \frac{k}{2}\rceil }\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)\right),\end{array}$$

(37)

where ⌊•⌋ is rounded down and ⌈•⌉ is rounded up. Thus, the high order moments could be computed analytically.

3.3. 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:

$$H\left(f_y\left(y\right)\right)=-{\displaystyle {\int }_{y\in {\Omega }_y}{f}_y\left(y\right)\ln \left(f_y\left(y\right)\right)}dy,$$

(38)

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:

$$\begin{array}{l}\underset{f_y(·)}{\max }H\left(f_y\left(y\right)\right)\\ s.t.{\displaystyle {\int }_{y\in {\Omega }_y}{f}_y\left(y\right)dy}=1\\ {\displaystyle {\int }_{y\in {\Omega }_y}{y}^k{f}_y\left(y\right)dy}={\mu }_k,k=1,2,\cdots ,N_{ME},\end{array}$$

(39)

where N_ME is the number of moments, and μ_k is the kth moment, i.e., μ₁ = E(y), μ₂ = E(y²), etc.

It can be proved that the optimal solution for the PDF formula is given as:

$$f_y\left(y\right)=\exp \left(-{\displaystyle \sum _{i=0}^{N_{ME}}{\lambda }_i{y}^i}\right).$$

(40)

Thus, λ_i in Equation (40) should satisfy the following equation:

$$\begin{array}{l}{\displaystyle {\int }_{y\in {\Omega }_y}{e}^{-{\displaystyle \sum _{i=0}^N{\lambda }_i}}dy}=1\\ {\displaystyle {\int }_{y\in {\Omega }_y}{y}^k{e}^{-{\displaystyle \sum _{i=0}^N{\lambda }_i}}dy}={\mu }_k;k=1,2,\cdots ,N.\end{array}$$

(41)

Specifically, if the first two moments of the output are available, we have:

$${\lambda }_0^{\ast }=\frac{E{\left(y\right)}^2}{2V\left(y\right)}+\frac{1}{2}\ln \left(2\pi V\left(y\right)\right),{\lambda }_1^{\ast }=-\frac{E\left(y\right)}{V\left(y\right)},{\lambda }_2^{\ast }=\frac{1}{2V\left(y\right)},$$

(42)

where $V\left(y\right)=E\left(y^2\right)-E{\left(y\right)}^2$ and the failure probability can be calculated as:

$$P_f=\frac{1}{2}\frac{\sqrt{\pi }\exp \left(\frac{1}{4}\frac{{\left({\lambda }_1^{\ast }\right)}^2}{{\lambda }_2^{\ast }}\right)\left(erf\left(\frac{1}{2}\frac{{\lambda }_1^{\ast }}{\sqrt{{\lambda }_2^{\ast }}}\right)+1\right)\exp \left(-{\lambda }_0^{\ast }\right)}{\sqrt{{\lambda }_2^{\ast }}}.$$

(43)

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.

3.4. Global Reliability Sensitivity Indices Calculation

Let ${\left\{{\xi }_{\mathit{I}}^{\left(i\right)}\right\}}_{i=1}^{M_R}$ 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 ${\left\{{\left[\begin{array}{ccc}{\psi }_0\left({\mathit{\xi }}_{\mathit{I}}^{\left(i\right)}\right)& \cdots & {\psi }_N\left({\mathit{\xi }}_{\mathit{I}}^{\left(i\right)}\right)\end{array}\right]}^T\right\}}_{i=1}^{M_R}$ and ${\left\{P_f\left({\mathit{\xi }}_{\mathit{I}}^{\left(i\right)}\right)\right\}}_{i=1}^{M_R}$. 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:

$$P_{f|X_{\mathit{I}}}=E\left[I_F|{\mathit{x}}_{\mathit{I}}\right]=P_f\left({\mathit{x}}_{\mathit{I}}\right)=P_f\left({\mathit{\xi }}_{\mathit{I}}\right)={\displaystyle \sum _{j=0}^{N_R-1}{c}_{R,j}}{\psi }_j\left({\mathit{\xi }}_{\mathit{I}}\right),$$

(44)

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:

$$N_R=\frac{\left(\left|\mathit{I}\right|+p_R\right)!}{\left|\mathit{I}\right|!p_R!},$$

(45)

where p_R is the predefined degree of R-PCE. Then, the CFP is:

$$P_f=E\left(P_f\left({\mathit{\xi }}_{\mathit{I}}\right)\right)=c_{R,0}.$$

(46)

And the variance of CFP is:

$$V\left(P_f\left({\mathit{\xi }}_{\mathit{I}}\right)\right)={\displaystyle \sum _{j=1}^{N_R-1}{c}_{R,j}^2}E\left({\psi }_j^2\left({\mathit{\xi }}_{\mathit{I}}\right)\right).$$

(47)

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.

3.5. 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} for main effect indices and I = {1, 2, ..., n}/{i} for total effect indices. Specifically, for each sample set $\left\{{\mathit{\xi }}_{\mathit{I}}^{\left(j\right)}\right\}$, substitute it into Equations (24)–(26) to obtain C-PCE $g_{PCE}^{\left(j\right)}\left(\mathit{\xi }|{\mathit{\xi }}_{\mathit{I}}\right)$. With the moments given by C-PCE, we could calculate the $P_f\left({\mathit{\xi }}_{\mathit{I}}^{\left(j\right)}\right)$ with ME method analytically. Third, the R-PCE comes out as Equation (44).

Step 4:

Compute the global reliability sensitivity indices by Equation (22), Equation (23) and Equations (46) and (47).

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 $\frac{2\left(n+p\right)!}{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 $\frac{2\left(n+p\right)!}{n!p!}+2n{p}_R+\frac{2\left(n+p_R-1\right)!n}{\left(n-1\right)!p_R!}+3+n$, 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.

4. 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.

4.1. A Numerical Example

Consider a limit-state function:

$$g\left(\mathit{x}\right)=10-18x_1+x_2^2+x_2+x_3^2+5x_3,$$

(48)

where the input variables follow the normal distribution or the lognormal distribution, i.e., lnx₁~$N$ (0.5,0.1), lnx₂~$N$ (0.6,0.1), and lnx₃~$N$ (0.5,0.1). The proposed method is employed to compute the main effect indices $S_i^{\left(R\right)}$ and total effect indices $S_{T_i}^{\left(R\right)}$, 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. The global reliability sensitivity indices of example 1.

Method		x1	x2	x3	Times of Function Evaluations
Full-scale MCS	$S_i^{\left(R\right)}$	0.7120	0.0005274	0.005265	1.5 × 107
	$S_{T_i}^{\left(R\right)}$	0.9946	0.1130	0.2602
IS	$S_i^{\left(R\right)}$	0.7206	0.0005147	0.005327	1.2 × 105
	$S_{T_i}^{\left(R\right)}$	0.9945	0.1124	0.2583
Proposed method	$S_i^{\left(R\right)}$	0.7181	0.0005426	0.005139	182
	$S_{T_i}^{\left(R\right)}$	0.9932	0.1119	0.2625

. 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_i^{\left(R\right)}$ and the total effect indices $S_{T_i}^{\left(R\right)}$ produce the same importance ranking: x₁ > x₃ > x₂, 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₁.

4.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 ${\omega }_{\max }$ of such a crane beam could be computed as:

$${\omega }_{\max }\left(q,E,I,L,F_p\right)=\frac{L^3}{EI}\left(\frac{F_p}{48}+\frac{5qL}{384}\right),$$

(49)

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. The distribution parameters of the input variables of example 2.

Input Variables		q/(N/M)	E/(kN/m2)	I/m4	L/m	Fp
Distribution type		Normal	Normal	Normal	Lognormal	Lognormal
Mean value		662.5	2 × 1011	2.172 × 10−4	10	3 × 104
variances	Case 1	0.1	0.08	0.1	0.06	0.1
	Case 2	0.05	0.06	0.05	0.03	0.04

.

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:

$$g\left(q,E,I,L,F_p\right)=0.0018-\frac{{\omega }_{\max }}{L}=0.0018-\frac{L^2}{EI}\left(\frac{F_p}{48}+\frac{5qL}{384}\right).$$

(50)

The results obtained by the full-scale MCS method, the IS method and the proposed method are given in

Table 3. The global reliability sensitivity indices of case 1.

Method		q	E	I	L	Fp	Times of Function Evaluations
Full-scale MCS	$S_i^{\left(R\right)}$	0.001283	0.05398	0.1029	0.1446	0.05583	1.5 × 107
	$S_{T_i}^{\left(R\right)}$	0.07086	0.4701	0.5967	0.6891	0.4895
IS	$S_i^{\left(R\right)}$	0.001267	0.05103	0.09755	0.1397	0.05129	1.2 × 105
	$S_{T_i}^{\left(R\right)}$	0.06961	0.4727	0.6026	0.6985	0.4905
Proposed method	$S_i^{\left(R\right)}$	0.001118	0.05404	0.1031	0.1445	0.05591	3882
	$S_{T_i}^{\left(R\right)}$	0.07251	0.4714	0.5879	0.6865	0.4983

and

Table 4. The global reliability sensitivity indices of case 2.

Method		q	E	I	L	Fp	Times of Function Evaluations
Full-scale MCS	$S_i^{\left(R\right)}$	0.000486	0.03123	0.01150	0.02078	0.002947	1.5 × 107
	$S_{T_i}^{\left(R\right)}$	0.09433	0.8405	0.7389	0.8240	0.5270
IS	$S_i^{\left(R\right)}$	0.000502	0.03111	0.01194	0.02260	0.002806	1.2 × 105
	$S_{T_i}^{\left(R\right)}$	0.09845	0.8509	0.7356	0.8333	0.5336
Proposed method	$S_i^{\left(R\right)}$	0.000532	0.03064	0.01183	0.02047	0.003107	3882
	$S_{T_i}^{\left(R\right)}$	0.09512	0.8497	0.7298	0.8199	0.5277

. 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_i^{\left(R\right)}$ and $S_{T_i}^{\left(R\right)}$ 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.

Moreover, the global reliability sensitivity indices are also schematically shown in Figure 4a and Figure 5a. In case 1, the importance ranking based upon global reliability sensitivity indices $S_i^{\left(R\right)}$ and $S_{T_i}^{\left(R\right)}$ are given as L > I > F_p > E > q. It indicates that the reliability could be increase mostly by decreasing the uncertainty of the span length. Besides, since the total effect indices of input variables are much higher than their main effect indices, there are strong interaction effects among these variables. In case 2, the importance ranking based upon global reliability sensitivity indices $S_i^{\left(R\right)}$ and $S_{T_i}^{\left(R\right)}$ 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.

The global sensitivity indices of maximal deflection of cases 1 and 2 are also computed in

Table 5. The global sensitivity indices of example 2.

		q	E	I	L	Fp
Case 1	$S_i$	0.003538	0.15818	0.253083	0.387629	0.183083
	$S_{T_i}$	0.003718	0.16375	0.260962	0.397182	0.188962
Case 2	$S_i$	0.003182	0.317395	0.218915	0.349925	0.106503
	$S_{T_i}$	0.003227	0.319898	0.220892	0.352517	0.107553

by traditional GSA method mentioned in Section 2.1, and they are plotted in Figure 4b and Figure 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.

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.

5. 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.