Sparse grid integration based solutions for moment-independent importance measures

Probabilistic Engineering Mechanics 39 (2015) 46–55

Contents lists available at ScienceDirect

Probabilistic Engineering Mechanics journal homepage: www.elsevier.com/locate/probengmech

Sparse grid integration based solutions for moment-independent importance measures Changcong Zhou a, Zhenzhou Lu b,n, Wei Li b a b

School of Mechanics, Civil Engineering and Architecture, Northwestern Polytechnical University, Xi'an, China School of Aeronautics, Northwestern Polytechnical University, Xi'an, China

art ic l e i nf o

a b s t r a c t

Article history: Received 21 July 2013 Received in revised form 16 December 2014 Accepted 16 December 2014 Available online 18 December 2014

Different importance measures exist in the literature, aiming to quantify the contributions of model inputs to the output uncertainty. Among them, the moment-independent importance measures consider the entire model output without dependence on any of its particular moments (e.g., variance), and have attracted growing attention among both academics and practitioners. However, until now robust and efficient computational methods for moment-independent importance measures are unavailable in the literature. To properly address this issue, a new computational method based on sparse grid integration (SGI) is proposed to perform moment-independent importance analysis in an efficient framework. In the proposed method, SGI is combined with the moment method to obtain the conditional and unconditional distributions of the model output, which are further used to estimate the moment-independent importance measures. The proposed method makes full use of the advantages of the SGI technique, and is easy to implement. Numerical and engineering examples are studied in this work and the results have demonstrated that the proposed method is feasible and applicable for practical use. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Uncertainty Sensitivity analysis Importance measure Moment-independent Sparse grid

1. Introduction As the complexity of mathematical models in the engineering keeps increasing, it is necessary to take all the inputs into consideration thoroughly and find which of them are more important to the model output [1,2]. Global sensitivity analysis is such a tool to study “how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input factors” [3]. Global sensitivity indices are also known as importance measures, which generally include nonparametric techniques [4], variance based methods studied by Sobol [5], Rabitz and Alis [6], Saltelli et al. [7,8], Frey and Patil [9], and moment-independent approaches proposed by Borgonovo [10–13], Liu and Homma [14]. Saltelli [3] underlined that an importance measure should satisfy the requirement of being “global, quantitative, and model free”. Borgonovo [10] extended Saltelli's three requirements by adding “moment independent”, and proposed a new importance measure which looks at the influence of input uncertainty on the entire output distribution without reference to a specific moment of the output. In a similar manner, Liu and Homma [14] proposed another moment-independent importance measure, which describes the contribution of the n

Corresponding author. E-mail address: [email protected] (Z. Lu).

http://dx.doi.org/10.1016/j.probengmech.2014.12.002 0266-8920/& 2014 Elsevier Ltd. All rights reserved.

input uncertainty on the cumulative distribution function (CDF) of the model output, while Borgonovo's importance measure studies the contribution on the probability density function (PDF) of the model output. A key concern in performing importance analysis is to improve the computational efficiency [15], i.e. to obtain results with a modest number of model evaluations, especially in engineering cases where the models involved usually take a long processing time. Numerical simulation methods have been adopted by Borgonovo [10], Liu and Homma [14,16] to obtain the two momentindependent importance measures proposed by them respectively, of which the computational burden is considerably large as double loop samplings are involved. In this respect, more efficient algorithms should be developed. According to the work by Liu and Homma [16], the calculation of the two importance measures can be, as a matter of fact, transformed into calculations of CDFs. Simply speaking, once the unconditional and conditional CDFs of the model output are obtained, both the two moment-independent importance measures can be readily obtained. In this work, we develop a new method to obtain the two moment-independent importance measures by using the sparse grid integration (SGI) technique. The sparse grid technique [17], which is based on one-dimensional formulas and then extended to higher dimensions, has been extensively utilized for high dimensional multivariate integration [18], as well as interpolation [19]. To a certain extent, this technique avoids the “curse of dimension”

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

of conventional integration algorithms meaning that the computational cost grows exponentially with the dimension of the problem. In the SGI technique, multivariate quadrature formulas are constructed with combinations of tensor products of suitable onedimensional formulas, thus the function evaluations needed and the numerical accuracy become independent of the dimension of the problem up to logarithmic factors. Existing literature has reported the high efficiency and accuracy of sparse grid applications. The SGI technique is combined with moment methods [26–28] in the proposed framework, where the first four moments of functions will be estimated by SGI and further used to obtain the unconditional and conditional CDFs of the model output. The obtained CDFs will be finally utilized to calculate the importance measures. The remainder of the paper is organized as follows. Section 2 reviews the definitions of the two moment-independent importance measures. This is followed by a brief introduction of the mathematical formulas of the SGI technique in Section 3. In Section 4, we propose a new method based on SGI to obtain the moment-independent importance measures more efficiently. In Section 5, numerical and engineering examples are presented to demonstrate the feasibility of the proposed method. Finally, conclusions of this work are highlighted in Section 6.

2. Two moment-independent importance measures

Suppose the model is denoted by Y = g (X), where Y is the model output, and X = [X1, X2 , … , Xd ]T (d is the number of inputs) is the set of uncertain inputs. The uncertainties of the inputs are represented by probability distributions. The unconditional PDF and CDF of Y are denoted as fY (y) and FY (y), respectively. The conditional PDF and CDF of Y are denoted as fY | Xi (y) and FY | Xi (y) , which are obtained by fixing the input of interest Xi at its realization x i⁎ generated from its distribution. The conditional and unconditional distribution functions of the output are depicted in Fig. 1. Borgonovo [10] proposed his importance measure by considering the contribution of inputs to the PDF of the model output. The shift between fY (y) and fY | Xi (y) is measured by the area s (Xi ) closed by the two PDFs, which is given as follows:

s (X i ) =

The expected shift can be obtained by varying the value of Xi over its distribution range, i.e.

E Xi [s (Xi )] =

∫ |fY (y) − fY|X (y) | dy i

(1)

∫ f X (x i ) s (X i ) d x i i

(2)

where fXi (x i ) is the marginal PDF of Xi . Then Borgonovo's importance measure is defined as [10] 1

δi = 2 E Xi [s (Xi )]

(3)

In a similar manner, Liu and Homma [14] proposed a new moment-independent importance measure concerning the CDF of model output. The deviation of FY | Xi (y) from FY (y) is measured by using the area A (Xi ) closed by the two CDFs, and can be calculated as

A (X i ) =

∫ | FY (y ) − FY | X (y ) | d y i

(4)

The expected deviation of FY | Xi (y) from FY (y) can be obtained as

E Xi [A (Xi )] =

∫ f X (x i ) A (X i ) d x i i

(5)

E Xi [A (Xi )] can be used to illustrate the influence of the input Xi on the output distribution. Furthermore, for a general model of which the unconditional expectation of the output E (Y ) is not equal to 0, Liu and Homma defined that the CDF based sensitivity indicator, Si(CDF ) , can be normalized as follows:

Si(CDF) =

2.1. Definitions

47

E Xi [A (Xi )] | E (Y ) |

(6)

For models of which E (Y ) is equal to 0, we can compare the influences of the input variables on the output CDF by directly using E Xi [A (Xi )]. Both δ i and Si(CDF ) evaluate the influence of the entire distribution of the input Xi on that of the model output Y , except that the former focuses on the PDF of output while the latter on the CDF. For both measures, the effect of Xi on the output uncertainty is estimated by varying all the other inputs over their variation range. The two moment-independent importance measures represent the expected shift in the decision maker’s view on the output uncertainty provoked by Xi . The inputs with larger values of δ i or Si(CDF ) will be associated with higher importance, indicating these inputs should be paid more attention to if we want to modify the performance of the model. In fact, the two measures work regardless of whether the model is linear or nonlinear, additive or nonadditive, and can be easily extended to a group of inputs, for which more detailed information can be found in the work of Borgonovo [10,11], Liu and Homma [14]. 2.2. Computational issues about the two measures Firstly, let us talk about how to obtain the value of δ i . It can be seen from Eqs. (1) and (2) that the calculation of δ i includes two major parts, i.e. the first would be to get the unconditional and conditional PDFs of the output to obtain s (Xi ), and the second is to calculate the expectation of s (Xi ). Now we describe the computational method utilized by Borgonovo [10] to obtain δ i , which generally can be viewed as a double-loop simulation method. First, the unconditional PDF of output, fY (y), is obtained by uncertainty propagation, and varying the inputs over their variation range. Second, one generates a value for Xi , namely x i1 sampled from fXi (x i ). Given this value, the conditional PDF of output, fY | Xi (y), is obtained by propagating uncertainty in the model fixing Xi = x i1. Then s (x i1) is computed from the numerical integration of the absolute value of the difference between fY | Xi (y) and fY (y). Repeat

Fig. 1. Unconditional and conditional distribution functions of the model output.

the procedure for a certain number of times, δ i can be finally estimated. In fact, the calculation of Si(CDF ) can be implemented in a

48

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

similar way, only needing to substitute the PDFs in Borgonovo's method with CDFs. Obviously, Borgonovo's method is a numerical simulation method. Liu and Homma [16] have discussed its computational cost. Suppose the sample size is N , then N model simulations are needed to estimate the unconditional distribution fY (y) (or FY (y)). For each input and for each value of this input, N model simulations are needed to estimate its conditional distribution fY | Xi (y) (or

FY | Xi (y) ). And to compute δ i , N numbers of fY | Xi (y) (or FY | Xi (y) ) are needed. Therefore, for a model with d inputs, a total number of N + d (N × N) model runs are needed to obtain the values of δ i for all inputs. It is clear that the computational cost will be tremendous even for a modest number of N . Liu and Homma [16] proposed a CDF-based method to calculate δ i , transforming the calculation of the area between fY (y) and fY | Xi (y) to that of vertical distances at some points between FY (y) and FY | Xi (y) . The CDF-based method is rather novel, since CDF can be more easily obtained than PDF. However, although Liu and Homma [16] proposed this novel method, they did not take full advantage of it. They have admitted that the efficiency improvement of their way to calculate δ i can be negligible for computationally intensive models. Still, the CDF-based method provides the researchers a brand-new viewpoint. As a matter of fact, the work in this paper is partly enlightened by Liu and Homma's method. According to Liu and Homma's method [16], suppose there are m (m ≥ 1) intersection points between fY (y) and fY | Xi (y), i.e.

y = a1, a2, … , am . If FY (a1) − FY | Xi (a1) > 0, then s (Xi ) can be obtained as

s (Xi ) = 2{[FY (a1) − FY | Xi (a1)] − [FY (a2 ) − FY | Xi (a2 )] + ⋯ + ( − 1)(m − 1) [FY (a m ) − FY | Xi (a m )]}

(7)

dimension of the problem. However, the SGI technique, which is based on Smolyak's construction [17], can overcome this curse of dimension to a certain extent. The underlying idea of sparse grid can be traced back to the Russian mathematician Smolyak [17], who used it for numerical integration, and it has been further improved for use of interpolation and approximation [18,19]. The Smolyak's construction uses a weighted linear combination of special tensor products to reduce the integration grid size, thus the number of function evaluations and the numerical accuracy become independent of the dimension of the problem up to logarithmic factors. The SGI has been demonstrated as an efficient discretization scheme which is well suited for high dimensional problems. It is necessary to point out that only the basic concepts are reviewed in this section, as technical details of the SGI fall out of the scope of this paper. More advanced sparse grid techniques are available in the literature, such as delayed construction [20], generalized construction [21], dimension-adaptive construction [22], adaptive construction [23], etc. Readers can refer to the literature for more details. Let U1i and w1i denote the quadrature points and weights in the one-dimension space, which can be obtained by Gaussian quadrature, trapezoidal, Clenshaw–Curtis rules, etc. [21]. Then the quadrature points Udk in the d-dimension space with k-level (k ≥ 0) accuracy generated by the SGI can be established by the Smolyak algorithm in the form of tensor product as follows:

Udk =

∪

k + 1 ≤| i |≤ q

U1i1 ⊗ U1i2 ⊗ ⋯ ⊗ U1id

(9)

where q = k + d and |i| denotes the sum of the multi-indices (|i| = i1 + , … , + id ). The sum is intelligently bounded by the restrictive condition k + 1 ≤ |i| ≤ q so as to exclude those points from full grids that contribute little to the improvement of the required accuracy. By use of the Smolyak algorithm [21], the weight wl corresponding to the lth quadrature point ξl = [ξ jii1 , … , ξ jiid ]T ∈ Udk 1

If FY (a1) − FY | Xi (a1) < 0, then

d

can be obtained as follows:

s (Xi ) = 2{[FY | Xi (a1) − FY (a1)] − [FY | Xi (a2 ) − FY (a2 )] + ⋯ + ( − 1)(m − 1) [FY | Xi (a m ) − FY (a m )]}

(8)

Until now, the calculation of s (Xi ) has been successfully converted to formulations involving vertical distances between FY (y) and FY | Xi (y) at the intersection points of the two PDFs. This method takes full advantage of the relationship between PDF and CDF. Then the job is to obtain FY (y) and FY | Xi (y) , and find the intersection

Si(CDF )

also relies on the calpoints. Meanwhile, the calculation of culation of FY (y) and FY | Xi (y) , judging from its definition illustrated by Eqs. (4)–(6). That is to say, if the unconditional and conditional CDF of the model output can be obtained more efficiently and conveniently, δ i and Si(CDF ) can be obtained more easily, and estimate of PDFs is no longer needed. Readers can refer to the work of Liu and Homma [16] for detailed information about this method. In this paper, we propose a new method based on sparse grid integration (SGI) technique, combined with moment methods, to provide a new way of calculating the moment-independent importance measures. In the next section, we will give a short introduction of SGI.

3. The sparse grid integration (SGI) technique Multivariate integrals arise extensively in structural reliability analysis, most of which have to be calculated numerically. Conventional algorithms for numerical computation of such integrals are often prohibited by the so-called “curse of dimension” [21], as the computational burden grows exponentially with the

⎛ d − 1⎞ i ωl = ( − 1)q−| i| ⎜ ⎟ (w 1 … w ijdid ) ⎝q − |i|⎠ ji1

(10)

Then, for the input–output function g (X) with d-dimension inputs X = [X1, X2 , … , Xd ]T discussed in Section 2, the integration of its expectation can be computed with SGI by the following equation with up to (2k + 1)-order polynomial accuracy

∫

Ndk

g ( x) f X ( x) d x ≈

Ndk

∑ wl g (T −1(ξl )) = ∑ wl g (x l ) l= 1

l= 1

(11)

where fX (x) is the joint PDF of the inputs and x = [x1, x2, … , x d ]T is an observation of X , Ndk is the number of all the quadrature points in Udk . T −1 (ξl ) is the inverse transformation of T (x l ) which transforms x space to ξ space, and the corresponding input for ξl is x l . In this work, we choose to use the Gaussian quadrature, which is widely used in the existent literature. Since the quadrature points and weights based on Gaussian quadrature are selected in the standard normal space, i.e. the ξ space, it is necessary to perform the isoprobabilistic transformation that maps the inputs from the original space into the standard normal space. In the original space, the input variables can be dependent or non-normally distributed. Generally, there are two ways to perform the isoprobabilistic transformation, i.e. Rosenblatt transformation or Nataf transformation [29,30]. Researchers have argued that the isoprobabilistic transformation will possibly introduce errors, and still work on this issue [31]. The isoprobabilistic transformation can be avoided if the quadrature points and weights are obtained according to the corresponding probability density function [24,25]. Besides, transformation from dependent variables to

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

independent ones also may potentially lead to significant degradation of the convergence properties of the solutions. This can be seen as a major limitation of the sparse grid methods. In this work, we use the Rosenblatt transformation when necessary. And for normally distributed inputs x of which the mean and standard deviation are μ X and σ X respectively, there exists the analytical transformation x l = T −1 (ξl ) = μ X + ξl Cσx , where Cσx is a d
d-dimension diagonal matrix with the elements of σx on the main diagonal. From the above discussion, it can be seen that the SGI technique is based on one-dimensional quadrature but extends it to higher dimensions in a more careful way than the full tensor product rule. SGI will exhibit its overwhelming advantage when applied to high dimensional problems (e.g. dZ 5), which has been proved by abundant research work. In this work, we will employ the SGI technique to estimate the first four moments of the model output, which will be further utilized to establish the conditional and unconditional CDFs.

4. A new method to solve moment-independent importance measures From the discussion in Section 2.2, s (Xi ), i.e. the shift between fY (y) and fY | Xi (y), can be estimated once FY (y) and FY | Xi (y) are obtained. So it is for A (Xi ). After s (Xi ) and A (Xi ) are obtained, their expectations have to be calculated to finally obtain δ i and Si(CDF ) . So, we conclude that there are two main issues to solve: the first is to obtain FY (y) and FY | Xi (y) with Xi fixed; and the second is to find a way to perform the estimate of the expectation. In this section, we will detail our method to properly address these issues.

4.1. Estimate of the unconditional CDF of model output

49

2 3(α4Z − 1) βSM + α3Z (βSM − 1)

βFM =

(5α32Z

− 9α4Z + 9) (1 − α4Z )

(16)

where β SM is the second moment reliability index and is obtained by

α1Z α2Z

βSM =

(17)

In Eqs. (16) and (17), αlZ (l = 1, 2, 3, 4) are the first four moments of Z (X , y), respectively. In fact, from Eq. (13) a clear relationship between the first four moments of Z (X , y) and those of g (X) can be established as

⎧ α1Z ⎪α ⎪ 2Z ⎨ ⎪ α3Z ⎪ ⎩ α4Z

= α1g − y = α2g = α3g = α4g

(18)

where αlg (l = 1, 2, 3, 4) denote the first four moments of g (X), respectively. It should be pointed out that for the probabilistic model defined by g (X), its central moments will not change once the distribution information of the inputs are known. Then from Eq. (18) it can be seen that the central moments of Z (X , y) keep changeless except that the α1Z varies with the changing value of y . After αlg (l = 1, 2, 3, 4) are obtained, an explicit mapping relationship between the values of y and FY (y) can be established from the propagation through Eqs. (15)–(18), which is denoted as

FY (y ) = φ (y )

(19)

where φ (⋅) denotes the mapping relationship. This means FY (y) can be seen as a univariate function of y if the values of αlg (l = 1, 2, 3, 4) are known. Now, we employ the SGI technique to estimate the values of αlg (l = 1, 2, 3, 4) according to the following quadrature formulas: p

As a matter of fact, the value of the output CDF at a given value of y is equal to the probability of the output smaller than y , i.e.

FY (y ) = P (Y ≤ y )

p

∫ (g (x) − α1g )2fX (x) dx]1/2 ≈ [ ∑ wl (g (xl ) − α1g )2]1/2 l= 1

(13)

α3g = (14)

which means the output CDF is equal to the failure probability of Z (X , y) at a given value of y . Over the recent years, researchers have been trying to evaluate the failure probability by directly utilizing the probabilistic moments of the output [26–28]. When the central moments of the output are obtained, the failure probability can be formulated as a function of these central moments. As the first two moments are generally inadequate, higher order moments are often used. The moment theory has no shortcomings associated with the design point and needs neither iteration nor complicated computation of derivatives [26]. According to the fourth moment method [26], the following equation holds

FY (y) = P f {Z (X , y)} ≈ Φ ( − βFM )

(20)

l= 1

α2g = [

Then the following equation holds

FY (y) = P f {Z (X , y)} = P {Z (X , y) ≤ 0}

∫ g (x) fX (x) dx ≈ ∑ wl g (xl )

(12)

Now, let us define a new performance function as

Z (X , y ) = g (X ) − y

α1g =

(15)

where β FM is the reliability index in the context of fourth moment method, and Φ (⋅) is the CDF of the standard normal distribution. β FM can be obtained as [26]

α4g =

1 α23g

1 α24g

(21)

p

∫ (g (x) − α1g )3fX (x) dx ≈ α13 ∑ wl (g (xl ) − α1g )3 2g l = 1

(22)

p

∫ (g (x) − α1g )4fX (x) dx ≈ α14 ∑ wl (g (xl ) − α1g )4 2g l = 1

(23)

In Eqs. (20)–(23), the d-dimension quadrature point x l and the associated weight wl (l = 1, 2, … , p, p is the number of quadrature points) are obtained by the Smolyak algorithm described in Section 3. Now, as αlg (l = 1, 2, 3, 4) have been conveniently estimated by SGI, the unconditional CDF of the model output, FY (y), can be readily obtained at any given values of y , without any further evaluation of the model g (X). In fact, the estimate of the conditional CDF of the model output, FY | Xi (y) , can be completed in almost the same manner, which will be discussed later. In addition, it is admitted that the moment method cannot always produce a better approximation of the CDF than those derived from the classical or other competing methods. Indeed, the moment method has its own limitations, such as the possibility that two stochastic responses sharing the same four moments

50

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

can have different CDFs. However, the moment method is competitive here because it can be more easily combined with SGI to reduce the computational cost.

dy

From Section 2, it can be seen that a key problem to obtain δ i and Si(CDF ) is to estimate the expectations of s (Xi ) and A (Xi ), respectively. In fact, the heavy computational burden by traditional ways to compute δ i and Si(CDF ) is partly due to the estimate of these expectations. One has to first generate a large number of samples according to the distribution of Xi , and then compute the corresponding values of s (Xi ) and A (Xi ) with Xi fixed at these samples, and then compute the expectations of s (Xi ) and A (Xi ) to estimate

δ i and Si(CDF ) . In other words, the computational burden is comprised by two parts, one is the inner calculating part to estimate the values of s (Xi ) and A (Xi ) with Xi fixed, and the other is the outer part to estimate the expectations of s (Xi ) and A (Xi ). Now, combined with the work in Section 4.1, we implement the SGI to deal with it. The values of s (Xi ) and A (Xi ) are changing with the varying values of the input Xi , thus they can be seen as univariate functions of Xi , respectively. Based on the Gauss–Hermite quadrature rule, firstly one-dimension quadrature points x ij and weights wij (j = 1, 2, … , m) will be generated for Xi , where m is the number of the quadrature points. Then the expectations of s (Xi ) and A (Xi ) can be estimated as follows: m

∫ fX (xi ) s (Xi ) dxi ≈ ∑ wi j s (xi j ) i

j=1

(26)

which can be further written as

d(FY (y) − FY | Xi (y))

4.2. Use of SGI to estimate δ i and Si(CDF )

E Xi [s (Xi )] =

fY (y) − fY | Xi (y) = 0

=0

(27)

Eq. (27) means that the intersection points between fY (y) and fY | Xi (y) are the extreme points of FY (y) − FY | Xi (y). Liu and Homma [16] utilized the numerical differential techniques to find the solutions of Eq. (27), which is constrained by the balance between the numerical accuracy and computational burden. In our method, however, the extreme points can be easily obtained as the values of FY (y) − FY | Xi (y) at any given values of y can be immediately obtained. Suppose N values of y in its range are considered, the corresponding values of FY (y) − FY | Xi (y) can be easily calculated. Then we can directly locate the extreme points by using the MATLAB functions, imregionalmax and imregionalmin. A larger value of N often leads to more accurate results, and the increase of computational burden is negligible, as no evaluation of the model is involved here. In this work, we take the value of N as 10,000. After we find the intersection points, s (Xi ) can be easily calculated by use of Eqs. (7) and (8). Until now, the values of s (Xi ) and A (Xi ) are obtained with Xi fixed at the value of x ij . Repeat the process for m times, and the expectations of s (Xi ) and A (Xi ) can be estimated by Eqs. (24) and (25). Then the two moment-independent importance measures, δ i and Si(CDF ) , can be finally obtained by Eqs. (3) and (6). Also reminded here is that the value of E (Y ) in Eq. (6) is already obtained by Eq. (20), i.e. E (Y ) = α1g .

(24) 4.3. Review on the proposed method

m

E Xi [A (Xi )] =

∫ fX (xi ) A (Xi ) dxi ≈ ∑ wi j A (xi j ) i

j=1

(25)

In principle, the value of m should be decided based upon the order of s (Xi ) and A (Xi ), however, this may be impractical since the explicit expressions of s (Xi ) and A (Xi ) with Xi cannot be obtained. What we are sure about is that both s (Xi ) and A (Xi ) are highly nonlinear univariate functions of Xi , as they denote the differences between PDFs and CDFs, and high order polynomials are often needed to approximate a probabilistic distribution. To ensure a convergent result, we assume the value of m as 15, which will theoretically reach 29-order polynomial accuracy. In fact, it can be seen later that a higher value of m contributes little to improving the accuracy of the proposed method. Secondly, obtain the conditional CDF of the model output, FY | Xi (y) , with Xi fixed at the value of x ij . This can be done in the same way of obtaining FY (y) in Section 4.1, except that the dimension now is d − 1 instead of d as the value of Xi is fixed. After FY | Xi (y) is obtained, we are ready to estimate s (Xi ) and A (Xi ) with Xi fixed at the value of x ij . Note that now the values of FY (y) and FY | Xi (y) at any given value of y can be immediately obtained with no further evaluations of the model, which has already been shown in Section 4.1. This guarantees that we can choose any available numerical method to estimate s (Xi ) and A (Xi ). Being measured by the area closed by FY | Xi (y) and FY (y), A (Xi ) can be easily estimated as the two CDFs have been already obtained. Now let us focus on the computation of s (Xi ), which is the area between fY (y) and fY | Xi (y). By the discussion in Section 2.2,

s (Xi ) can be converted to a distance metric between FY | Xi (y) and FY (y). Thus the job now is to find the intersection points between fY (y) and fY | Xi (y) so as to compute the values of FY (y) − FY | Xi (y) at these points. The intersection points can be determined by

Let us have a review of the proposed method to get a clearer understanding of the process. In the proposed method, the SGI technique is firstly employed to estimate the first four central moments of the model output, which are further used by the moment method to establish the unconditional CDF of the output. The method can be seen as a nested-loop procedure. In the outer loop, we view s (Xi ) and A (Xi ), which are the core parts of the importance measures, as univariate functions of Xi , and generate the one-dimension quadrature points and weights to estimate the expectations of s (Xi ) and A (Xi ). In the inner loop, at each of these quadrature points of Xi , the conditional CDF of the output is obtained. The unconditional and conditional CDFs of the output are then utilized to estimate s (Xi ) and A (Xi ), where the former should be firstly converted to a distance metric between the CDFs. After the values of s (Xi ) and A (Xi ) are all obtained with Xi fixed at all the quadrature points, the expectations of s (Xi ) and A (Xi ) can be easily estimated by Eqs. (24) and (25). Then the two importance measures, δ i and Si(CDF ) , can be calculated by their definitions. For the purpose of a better illustration, the main procedure of the proposed method can be summarized in two major steps and outlined in Fig. 2. In fact, the focus of the proposed method is computation of δ i . For a random quantity, the PDF is much more difficult to estimate accurately than the CDF. In the proposed method, the shift between PDFs has been successfully converted to formulations only involving CDFs, which can be efficiently estimated by SGI. As shown in the procedure of the proposed method, computation of Si(CDF ) is accomplished during that of δ i . In other words, estimate of

Si(CDF ) can be seen as a byproduct or bonus when estimating δ i with the proposed method. Another advantage of the proposed method over the previous methods used by Borgonovo [10], Liu and Homma [16] is the high computational efficiency. The model is only evaluated when the

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

51

Fig. 2. Flowchart of the proposed method.

first four moments of the output need to be estimated. As we have discussed in Section 2.2, given a sample size N ¼10,000, which Liu and Homma [16] adopted in their work, then a total of 6.0001 × 108 model runs are needed for a 6-dimenion problem. Obviously, the computational burden is unacceptable for engineering problems with implicit models. In our method, the relationship between the needed number of model evaluations, i.e.

Ncall , and the problem dimension, i.e. d , is depicted by Fig. 3 with the accuracy level k taken as 1 and 2. From the discussion in Section 3, 1-level accuracy and 2-level accuracy will yield 3-order and 5-order polynomial exactness, respectively, which should be adequate for most problems in real applications. From Fig. 3 it can be seen that the computational burden increases with the dimension in a gentle manner, demonstrating its applicability to highly dimensional problems.

5. Examples In this section, a numerical example with analytical results is discussed at first, to prove the correctness of the proposed method. Then, three engineering examples are studied to demonstrate the efficiency and applicability of the proposed method. For the engineering examples there are no analytical results due to the complexity of functions involved, we use the method proposed by Borgonovo [10], which is named as PDF-based method and shown at the beginning of Section 2.2 [16]. The computational costs of the proposed method as well as the PDF-based method are listed and compared when possible. In the PDF-based method, an adaptive kernel density estimator (KDE) proposed by Botev et al. [32] is adopted to estimate fY (y) and fY | Xi (y). In this work, proFig. 3. Computational efforts needed by the proposed method.

blems involving only independent variables are studied.

52

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

5.1. Example 1: analytical test case

Table 2 Distribution information of the inputs for the reinforced section example.

Borgonovo et al. [12] have derived the analytical results of δ i for specific cases with different model structures and input distributions. Here we discuss a linear model which is represented by sum of normal random variables. The relationship between model output and inputs is given as

Y = g (X ) = 1.5X1 + 1.6X2 + 1.7X 3 + 1.8X 4 + 1.9X5 + 2X6

(28)

where all the input variables are independent and follow the standard normal distribution, i.e., Xi ~N (0, 1) (i = 1, 2, … , 6). This example has also been discussed by Wei et al. [33]. Notice that the performance function in this case is simply a linear combination of independent and identically distributed inputs, thus qualitative information can be deduced that the relative importance of the inputs is decided by the coefficients with each variable in the function. In this way, we can predict the importance ranking of the inputs in descending order as X6 , X5, X 4 , X3, X2 , X1. Now we perform the moment-independent importance analysis with the proposed method. Meanwhile, it can be easily found that E (Y ) = 0 in this case, which means Si(CDF ) cannot be defined in this case. Thus we only study the results of δ i in this example. Analytical solutions as well as the results obtained by the singleloop MCS method proposed by Wei et al. [33] are listed in Table 1, for comparison with those obtained by the proposed method. Results in Table 1 show that the importance ranking of the inputs is the same with that predicted beforehand. Moreover, we notice that the proposed method obtains almost accurate results compared with the analytical solutions, and behaves better than the single-loop MCS method. The comparison has demonstrated that the proposed method is correct and viable. In this example there exist analytical solutions which barely cause any computational burden, thus it is not very meaningful to show the computational efficiency of the proposed method. However, in most of the real applications with complex input-output relationship, analytical results of the moment-independent importance measures are unavailable. In these cases, the proposed method will show its advantage both in accuracy and efficiency. More discussions can be found in the following engineering examples. 5.2. Example 2: a reinforced section In this example, the ultimate bending moment of a reinforced section is considered [34]. The ultimate bending moment can be obtained as

Y = g (X ) = X1X2 X 3 − X 4

X12 X22 X5 X6

Input X1 X2 X3 X4 X5 X6

(mm2) (N/mm2) (mm) (N/mm2) (mm)

Mean

Standard deviation

Distribution

1260 250 770 0.55 30 250

63 17.5 10 0.055 4.5 5

Lognormal Lognormal Lognormal Lognormal Lognormal Lognormal

Table 3 Importance measures of the inputs for the reinforced section example. Input δi

X1 X2 X3 X4 X5 X6

Si(CDF )

Proposed method

PDF-based method

Proposed method PDF-based method

0.222(2) 0.387(1) 0.0503(3) 0.0116(5) 0.0175(4) 0.00228(6)

0.225 0.378 0.0512 0.0123 0.0186 0.00250

0.0404(2) 0.0604(1) 0.0104(3) 0.00245(5) 0.00369(4) 0.000483(6)

maximum compressive strength of the concrete, and X6 is the width of the beam. The distribution information of the inputs is listed in Table 2. Both the proposed method and the PDF-based method are applied to obtain the two moment-independent measures, and the results are listed in Table 3. It can be seen from Table 3 that the results obtained by the proposed method are in almost total agreement with those obtained by the PDF-based method. Meanwhile, as we have discussed in Section 4.3, the computational efficiency of the proposed method is much higher than the PDF-based method. The proposed method needs 5575 model evaluations with 5-order polynomial exactness, which is tremendously smaller than the number of 6.0001 × 108 needed by the PDF-based method. It is necessary to mention that the value of m is taken as 15 for Eqs. (24) and (25) in this example, as we have discussed in Section 4.2. To take a look at the effect of m on the accuracy of the proposed method, In Fig. 4 we present the relationship between the value of m and the results of δ i obtained by the proposed method. It can be seen from Fig. 4 that a comparatively small value of m (e.g. 10) will be adequate for the proposed method to obtain an acceptable result. Similar conclusions exist for the other three examples in this work. In fact, the optimum value of m can be found in an iteration process by

(29)

where X1 is the area of reinforcement, X2 the yield stress of the reinforcement, X3 the effective depth of the reinforcement, X4 a factor related to the stress–strain relationship of concrete, X5 the Table 1 Importance measures of the inputs for the analytical test case. Input

δi Analytical results

X1 X2 X3 X4 X5 X6 a

a

0.119(6 ) 0.129(5) 0.138(4) 0.148(3) 0.158(2) 0.168(1)

Single-loop MCS

Proposed method

0.105(6) 0.113(5) 0.121(4) 0.127(3) 0.135(2) 0.151(1)

0.118(6) 0.128(5) 0.137(4) 0.147(3) 0.157(2) 0.168(1)

The number in the brackets denotes the ranking of the inputs.

0.0406 0.0592 0.0107 0.00257 0.00383 0.000512

Fig. 4. Relationship between m and δi .

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

judging the convergence of the results, which, however, will consume extra computational time. It is necessary to point out that the proposed method based on SGI is applicable unless there are very unconventional function expressions which is difficult to be approximated with polynomials. In the applicability scope of the proposed method, 15 would be an appropriate value for m in the proposed method. More conclusions can be drawn from Table 3. The ranking of inputs based on δ i is the same as that based on Si(CDF ) , and the inputs are ranked in decreasing order as X2, X1, X3, X5, X4, X6. The impact of X6 is nearly negligible compared with that of X2, X1, X3, no matter the CDF or PDF of the model output is regarded. 5.3. Example 3: a riveting process In the aircraft industry, sheet metal parts are widely used, and riveting is the most common method to assemble them [35,36]. Many factors have to be considered during the riveting process, which will directly affect the riveting quality. One main factor is the squeeze force. If the squeeze force is too low, the rivets may be too loose to efficiently join different parts, while a too high squeeze force may induce excessive residual stresses which will result in stress concentration (around the hole and rivet) and initial cracks. Both cases are dangerous for aircrafts. Therefore, it is of great significance to determine the squeeze force properly for the flight safety. The real riveting process is very complex. In this work, we take the headless rivet as an example [36], and generalize the riveting process into two stages as shown in Fig. 5. In stage I the rivet is punched from state A (the initial state of the rivet before impact, without any deformation) to state B (an intermediate state of the rivet when there is no clearance between rivet and hole), then in stage II the rivet is further punched from state B to state C (the final state of the rivet after impact, with rivet heads formed). Throughout the riveting process we assume that the hole diameter keeps unchanged. To establish the relationship between the geometric dimensions of a rivet and the squeeze force, the following ideal conditions should be assumed: 1) 2) 3) 4)

The hole is not enlarged in the riveting process. The change of the rivet volume in the process can be neglected. After impact, the rivet driven head has a cylindrical shape. The material of the rivet is isotropic.

The initial volume Vol0 of the rivet before the impact (in state A) is obtained as

Vol0 =

π 2 D h 4

(30)

where D and h are the rivet diameter and length in state A, respectively. At the end of stage I, in state B the volume of the rivet, Vol1, can be obtained as

Vol1 =

π 2 D0 h1 4

53

(31)

where D0 and h1 are the rivet diameter and length in state B, respectively. After stage II, we assume the top and bottom heads of the formed rivet is state C have the same dimension, then the volume of the rivet in state C, Vol2 , can be obtained as

Vol2 =

π 2 π D0 t + 2 D12 H 4 4

(32)

where t is the whole thickness of the two sheets, D1 and H are the diameter and height of the rivet head in state C, respectively. The maximum squeeze force needed in the riveting process can be expressed as

Fmax =

π 2 D1 σmax 4

(33)

where σmax is the maximum squeeze stress of the rivet head in its formation. According to power hardening theory, the true maximum squeeze stress in the y-direction can be obtained as

σmax = K (ε y )n SHE

(34)

where K is the strength coefficient, n SHE is the strain hardening exponent of the rivet material, and ε y is the true strain in the ydirection of the rivet head in its formation. And in our model, the true strain ε y is composed of two parts: the strain ε y1 caused in stage I and the strain ε y2 caused in stage II, then the true strain ε y can be expressed as follows:

ε y = ε y1 + ε y2

(35)

where ε y1 = ln (h/h1) and ε y2 = ln ((h1 − t)/2H). Combining Eqs. (30)–(35), under the ideal assumed conditions one can obtain the maximum squeeze force needed for a certain riveting process as follows:

Fmax =

n SHE D2h − D02 t ⎞ π D2h − D02 t ⎛ ⎟⎟ K ⎜⎜ ln 4 2H 2HD2 ⎠ ⎝

(36)

In this work, the material of the rivet is 2017-T4, and we already know n SHE = 0.15, H = 2.2 mm , and the other inputs are all assumed as normal random variables, of which the distribution parameters are listed in Table 4. The importance measures are obtained by the proposed method and the PDF-based method respectively, and listed in Table 5 for comparison. For this 5-dimension example, the number of model evaluations needed by the proposed method is 3136 with 5-order polynomial exactness, while that needed by the PDF-based method is 5.0001 × 108. Obviously, the proposed method has an overwhelming advantage in efficiency over the PDF-based method. The importance results can provide helpful interpretations of the riveting process. It can be seen that the rivet length in state A, h, contributes the most to the output distribution, no matter whether PDF or CDF is considered. The next is the strength coefficient, K . The contributions of the three inputs, D , t , D0 , are considerably

Fig. 5. Simplified riveting process.

54

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

Table 4 Distribution information of the inputs for the riveting example.

Table 6 Distribution information of the inputs for the ten-bar example.

Input

Mean

Standard deviation

Distribution

Variable

Mean

Standard deviation

Distribution

D (mm) h (mm) D0 (mm) t (mm) K (MPa)

5 20 5.1 5 547.2

0.01 0.3 0.0102 0.05 5.472

Normal Normal Normal Normal Normal

Ai (m2) L (m) P1 (kN) P2 (kN) P3 (kN) E (GPa)

1
10  3 1 80 10 10 100

5
10  5 0.05 4 0.5 0.5 5

Normal Normal Normal Normal Normal Normal

Table 5 Importance measures of the inputs for the riveting example. Input δi

D h D0 t K

Si(CDF )

Proposed method

PDF-based method

Proposed method

PDF-based method

0.0695(3) 0.485(1) 0.0187(5) 0.048(4) 0.134(2)

0.0688 0.482 0.0194 0.0483 0.134

0.00436(3) 0.0212(1) 0.00120(5) 0.00306(4) 0.00805(2)

0.00444 0.0210 0.00126 0.00313 0.00821

smaller than the other inputs. To improve the robustness of the riveting process, the inputs h and K should be paid extra attention to. Fig. 7. Importance measure results for the ten-bar example.

5.4. Example 4: a ten-bar structure Consider a ten-bar structure with 15 input variables [37], which is depicted in Fig. 6. The length and sectional area of the horizontal and vertical bars are denoted as L and Ai (i = 1, 2, … , 6), the length and sectional area of diagonal bars are 2 L and Ai (i = 7, 8, … , 10). The elastic modulus of all bars is denoted as E . P1, P2 and P3 are the external loads. In this problem, the 15 input variables are denoted as X = [Ai (i = 1, 2, … , 10), L, P1, P2, P3, E]T , with the distribution information listed in Table 6. Taking the perpendicular displacement of Node 2, D , not exceeding 0.004 m as the constraint condition, the model can be established as

Y = g (X ) = 0.004 − |D|

(37)

where D is an implicit function of the inputs. As the model in this example needs to be evaluated by calling ANSYS codes, the computing process will be very time-consuming for the PDF-based method. This issue will turn increasingly severe for engineering cases where most real models are time-consuming to evaluate. Thus in this example we only apply the proposed

method based on SGI to perform importance analysis for the tenbar structure, and results are presented in Fig. 7. The results obtained by the proposed method in Fig. 7 can reflect the contributions of the inputs to the output uncertainty from the viewpoint of probabilistic distribution. The inputs can be ranked according to the importance measures in descending order as L, E, P1, P2, A1, A3, A7, A8, P3, A9, A5, A10, A4, A6, A2. In fact, the contributions by A2, A4, A5, A6 and A10 are so tiny that they can almost be neglected compared to the other inputs. The importance analysis can thus provide helpful advice for the improvement of the model, and dig out some underlying information associated with the inputs. Additionally, an interesting phenomenon can be noticed from the results of the four examples, i.e. the importance rankings of the inputs based on δ i and Si(CDF ) are the same. This phenomenon seems rather reasonable. If one input has important effects on the PDF of the model output, then there is no obvious improperness to anticipate that this input may as well greatly affects the CDF of the model output, since PDF and CDF are closely related with each other. However, we should point out that this phenomenon is just possible rather than universal, as no theoretical evidence has been found to justify it. Besides, this work is dedicated to how to calculate the moment-independent importance measures more efficiently, maybe further work will consider whether these two measures can be substituted by each other or not.

6. Conclusions

Fig. 6. Diagram of the ten-bar structure.

For importance analysis of engineering models, two momentindependent importance measures, δ i and Si(CDF ) , have been proposed by researchers. These measures can reflect the contribution by the inputs to the entire distribution of the model output. The information provided by the importance measures can be used to identify those important factors, thus guide the improvement of engineering models. Despite the advantages of the moment-independent importance measures, the application is constrained by the heavy computational burden if traditional methods are used.

C. Zhou et al. / Probabilistic Engineering Mechanics 39 (2015) 46–55

The traditional methods proposed are generally based on numerical simulations, which unavoidably generate lots of “unnecessary” samples which contribute little to the final results but still consume valuable computational time. To address this issue, in this work we propose a new method based on SGI to obtain δ i and Si(CDF ) simultaneously. Results of the engineering examples have demonstrated that the proposed method is of good applicability, being very economic in terms of the number of model evaluations while ensuring a controllable accuracy level. The proposed method offers not a substitution but a viable alternative for variable screening and ranking in reliability analysis. However, for some cases with complex functions, there may be a trade-off between accuracy and efficiency. This is mainly because that SGI will indeed be constrained by the smoothness or irregularity of the functions involved. Thus, it is necessary to point out that behavior of the proposed method is to a large extent decided by the accuracy of SGI. Meanwhile, in the present work, only the basic sparse gird technique is reviewed and used by the proposed method. The proposed method may meet trouble when dealing with problems with correlated inputs or incomplete distribution information, as isoprobabilistic transformations are needed during the procedure of the proposed method, which may induce extra errors. However, the proposed method is evolvable, since more advanced SGI techniques that are and will be developed by researchers can be easily introduced into the framework of the proposed method to improve its accuracy and robustness. The proposed method is a novel and effective estimator as long as it is properly used in right circumstances.

Acknowledgments The authors gratefully acknowledge the supports of the Nature Science Foundation of China (51175425), and Research Fund for the Doctoral Program of Higher Education of China (20116102110003). They would also like to thank the anonymous referees for the perceptive comments and useful suggestions.

References [1] G. Li, H. Rabitz, P.E. Yelvington, O.O. Oluwole, F. Bacon, C.E. Kolb, J. Schoendorf, Global sensitivity analysis for systems with independent and/or correlated inputs, J. Phys. Chem. A 114 (19) (2010) 6022–6032. [2] S. Rahman, Global sensitivity analysis by polynomial dimensional decomposition, Reliab. Eng. Syst. Saf. 96 (7) (2011) 825–837. [3] A. Saltelli, Sensitivity analysis for importance assessment, Risk Anal. 22 (3) (2002) 579–590. [4] R.L. Iman, M.E. Johnson, C.C. Watson Jr., Sensitivity analysis for computer model projections of hurricane losses, Risk Anal. 25 (5) (2005) 1277–1297. [5] I.M. Sobol, Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simul. 55 (1–3) (2001) 271–280. [6] H. Rabitz, O.F. Alis, General foundations of high-dimensional model representations, J. Math. Chem. 25 (2–3) (1999) 197–233. [7] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola, Variance based sensitivity analysis of model output. Design and estimator for the total

55

sensitivity index, Comput. Phys. Commun. 181 (2) (2010) 259–270. [8] A. Saltelli, S. Tarantola, F. Campolongo, Sensitivity analysis as an ingredient of modelling, Stat. Sci. 15 (4) (2000) 377–395. [9] C.H. Frey, S.R. Patil, Identification and review of sensitivity analysis methods, Risk Anal. 22 (3) (2002) 553–578. [10] E. Borgonovo, A new uncertainty importance measure, Reliab. Eng. Syst. Saf. 92 (6) (2007) 771–784. [11] E. Borgonovo, Measuring uncertainty importance: investigation and comparison of alternative approaches, Risk Anal. 26 (5) (2006) 1349–1361. [12] E. Borgonovo, W. Castaings, S. Tarantola, Moment independent importance measures: new results and analytical test cases, Risk Anal. 31 (3) (2011) 404–428. [13] E. Borgonovo, S. Tarantola, Moment independent and variance-based sensitivity analysis with correlations: an application to the stability of a chemical reactor, Int. J. Chem. Kinet. 40 (11) (2008) 687–698. [14] Q. Liu, T. Homma, A new importance measure for sensitivity analysis, J. Nucl. Sci. Technol. 47 (2010) 53–61. [15] A. Saltelli, Making best use of model evaluations to compute sensitivity indices, Comput. Phys. Commun. 145 (2) (2002) 280–297. [16] Q. Liu, T. Homma, A new computational method of a moment-independent uncertainty importance measure, Reliab. Eng. Syst. Saf. 94 (7) (2009) 1205–1211. [17] S.A. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, Sov. Math. Dokl. 4 (1963) 240–243. [18] F. Heiss, V. Winschel, Likelihood approximation by numerical integration on sparse grids, J. Econom. 144 (1) (2008) 62–80. [19] V. Barthelmann, E. Novak, K. Ritter, High dimensional polynomial interpolation on sparse grids, Adv. Comput. Math. 12 (4) (2000) 273–288. [20] K. Petras, Smolyak cubature of given polynomial degree with few nodes for increasing dimension, Numer. Math. 93 (4) (2003) 729–753. [21] T. Gerstner, M. Griebel, Numerical integration using sparse grids, Number Algorithms 18 (3–4) (1998) 209–232. [22] T. Gerstner, M. Griebel, Dimension-adaptive tensor-product quadrature, Computing 71 (1) (2003) 65–78. [23] M. Hegland, Adaptive sparse grids, ANZIAM J. 44 (2002) 335–353. [24] W. Gautsch, Orthogonal Polynomials: Computation and Approximation, Oxford University Press, London, 2004. [25] S. Rahman, Extended polynomial dimensional decomposition for arbitrary probability distributions, J. Eng. Mech. 135 (12) (2009) 1439–1451. [26] Y.G. Zhao, T. Ono, Moment methods for structural reliability, Struct. Saf. 23 (1) (2001) 47–75. [27] H.P. Hong, Point-estimate moment-based reliability analysis, Civ. Eng. Syst. 13 (4) (1996) 281–294. [28] C. Zhou, Z. Lu, L. Li, J. Feng, B. Wang, A new algorithm for variance based importance analysis of models with correlated inputs, Appl. Math. Model. 37 (3) (2013) 864–875. [29] P.L. Liu, A. Kiureghian, Multivariate distribution models with prescribed marginals and covariances, Probab. Eng. Mech. 1 (2) (1986) 105–112. [30] R. Lebrun, A. Dutfoy, An innovating analysis of the Nataf transformation from the copula viewpoint, Probab. Eng. Mech. 24 (3) (2009) 312–320. [31] R. Lebrun, A. Dutfoy, Do Rosenblatt and Nataf isoprobabilistic transformations really differ? Probab. Eng. Mech. 24 (4) (2009) 577–584. [32] Z.I. Botev, J.F. Grotowski, D.P. Kroese, Kernel density estimation via diffusion, Ann. Stat. 38 (5) (2010) 2916–2957. [33] P. Wei, Z. Lu, X. Yuan, Monte Carlo simulation for moment-independent sensitivity analysis, Reliab. Eng. Syst. Saf. 110 (2013) 60–67. [34] J.H. Zhou, A.S. Nowak, Integration formulas to evaluate functions of random variables, Struct. Saf. 5 (4) (1988) 267–284. [35] S.H. Cheraghi, Effect of variations in the riveting process on the quality of riveted joints, Int. J. Adv. Manuf. Technol. 39 (11–12) (2008) 1144–1155. [36] W. Mu, Y. Li, K. Zhang, H. Cheng, Mathematical modeling and simulation analysis of flush rivet pressing force, J. Northwest. Polytech. Univ. 28 (2010) 742–747 (in Chinese). [37] L. Li, Z. Lu, J. Feng, B. Wang, Moment-independent importance measure of basic variable and its state dependent parameter solution, Struct. Saf. 38 (2012) 40–47.