Uncertainty propagation analysis by an extended sparse grid technique

Abstract

In this paper, an uncertainty propagation analysis method is developed based on an extended sparse grid technique and maximum entropy principle, aiming at improving the solving accuracy of the high-order moments and hence the fitting accuracy of the probability density function (PDF) of the system response. The proposed method incorporates the extended Gauss integration into the uncertainty propagation analysis. Moreover, assisted by the Rosenblatt transformation, the various types of extended integration points are transformed into the extended Gauss-Hermite integration points, which makes the method suitable for any type of continuous distribution. Subsequently, within the sparse grid numerical integration framework, the statistical moments of the system response are obtained based on the transformed points. Furthermore, based on the maximum entropy principle, the obtained first four-order statistical moments are used to fit the PDF of the system response. Finally, three numerical examples are investigated to demonstrate the effectiveness of the proposed method, which includes two mathematical problems with explicit expressions and an engineering application with a black-box model.

Key words： uncertainty propagation analysis; extended sparse grid; maximum entropy principle; extended Gauss integration; Rosenblatt transformation; high-order moments analysis

Introduction

In practical engineering problems, many inevitable uncertainty factors exist due to the complexity of structures, e.g., dispersion of materials, and errors in manufacturing, installation, and measurement. The coupling of multiple uncertainties often results in an insignificant deviation of the system response. The influence of the parameter uncertainty on the system response is called uncertainty propagation (UP) [¹], which is essential in designing structures or products [²–⁴]. The statistical moments are important characteristics of a system response’s probability distribution. Thus, solving the statistical moments of the system response has become a key issue in UP [⁵–¹⁰].

Some existing methods, namely, Monte Carlo simulation (MCS) [¹¹–¹³], Taylor expansion-based method [¹⁴,¹⁵], reliability-based method [¹⁶–¹⁸], polynomial chaos expansion (PCE) [¹⁹–²¹] and numerical integration-based method [¹,²²,²³] have already contributed to UP. The MCS method was first used in simulating a neutron chain reaction. Its idea is relatively simple, with simulation results inclined to theoretical solutions with the increasing samples. Therefore, MCS results are usually employed as a reference to verify the accuracy of other methods, but it requires a large number of samples to ensure acceptable accuracy, thereby causing an extremely low computational efficiency. The Taylor expansion-based method executes a Taylor expansion on the system response function by input random variables. Based on the Taylor expansion, the system response can be approximately obtained and its probability distribution conveniently calculated. However, when the response function has a relatively large nonlinearity with respect to the random variables, its accuracy may not be guaranteed. Therefore, a standard reliability-based method must utilize the first-order reliability method or the second-order reliability method, through which the failure probability under different response thresholds can be obtained and the cumulative distribution function (CDF) of the system response solved. However, this standard reliability-based method must also solve the most probable point for each response threshold, which significantly increases the computational cost. The PCE method adopts various types of orthogonal polynomials to approximate the system response according to the probability distributions of the input variables and is applicable to problems with relatively large variances of input random variables. Furthermore, this method can be improved as the order of the orthogonal polynomial increases but may have a low computational efficiency at the same time. The numerical integration-based method for the UP analysis is becoming increasingly popular. This method aims to solve the statistical moments of the system response by numerical integration and then estimate the probability distribution of the system response based on the calculated moments. The univariate dimension reduction method (UDRM) [²²] is one of the important methods in this category, which breaks down the original system response function into a series of univariate subsystem response functions and then obtains the moments of the original system response by solving the moments of the subsystem response function. UDRM has high computational efficiency for UP problems. However, its computational accuracy may not be satisfactory when the response function comprises strong interactions among the input random variables.

In recent years, another type of numerical integration technique, that is, the sparse grid numerical integration (SGNI), was introduced into the UP field [²⁴–²⁶], which was originally proposed by Smolyak [²⁷] and revised by Novak and Ritter [²⁸,²⁹]. The main concept of this method is to combine one-dimensional integration nodes into multidimensional nodes according to a specific tensor product rule and then solve the statistical moments with the multidimensional nodes. Compared with the direct integration method [³⁰], SGNI can prevent the occurrence of major dimension problems in high-dimensional systems. Furthermore, Ref. [²⁵] has demonstrated that the SGNI method can solve problems with strong interactions in the system response, which gives the method a significant application potential for UP analysis. Therefore, the SGNI method has been applied in fields such as aerospace, electronics, and others [³¹–³³]. Nevertheless, for practical cases, such as problems regarding high-dimensional variables and strong nonlinearity or a black-box model, the accuracy of the traditional SGNI (e.g., the SGNI with traditional univariate integration points) for UP analysis needs improvement, especially in solving high-order moments such as skewness and kurtosis of the system response. Although studies [³⁴] on the Kalman filter have provided univariate integration perspectives to improve the solving accuracy, it is limited in arbitrary distributions of input random variables, thereby affecting the practicability of the method.

Based on the traditional SGNI framework, this study aims to develop an efficient method for uncertainty propagation analysis to improve the accuracy of the high-order moments and fitting accuracy of the probability density function (PDF) of the system response regarding problems with arbitrary continuous distribution of input variables. In this method, all integration nodes of the input variables are transformed into the extended Gauss-Hermite integration (EGHI) nodes. Furthermore, the sparse grid technique with transformed nodes is extended to estimate the statistical moments. The maximum entropy principle is employed to calculate the PDF of the system response.

The rest of this paper is organized as follows: Section 2 provides a brief introduction to UP and SGNI; Section 3 presents a formulation of the proposed method; Section 4 analyzes the numerical examples; and Section 5 summarizes the study.

Uncertainty propagation and sparse grid method

For a UP problem, a response function or a performance function must first be established as follows:

(1)Y=g(\mathit{X}),

where \mathit{X}=(X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{...},X_n) denotes an n-dimensional vector of random variables and Y is the system response. The uncertainty propagation analysis obtains the PDF f_Y(y) of the system response through the joint PDF f_X(x) of the input random variables. In general, the statistical moments are important digital characteristics for UP analysis. Therefore, the t-order origin moments {\mu }_t and t-order central moments c_t are often used to represent the distribution of the system response:

(2){\mu }_t={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t(\mathit{X})}{f}_X(x)\mathrm{d}\mathit{X},\hspace{0.17em}\hspace{0.17em}t=\mathrm{1},\mathrm{2},\dots ,

(3)c_t={\underset{\mathrm{\Omega }}{\int }(g(\mathit{X})-{\mu }_{\mathrm{1}})}^t{f}_X(x)\mathrm{d}\mathit{X},\hspace{0.17em}\hspace{0.17em}t=\mathrm{2},\mathrm{3},\dots ,

where \mathbf{\Omega } is the integration space of the random variables. In addition, the mean \mathit{\mu }, standard deviation \mathit{\sigma }, skewness \mathit{\tau }, and kurtosis \mathit{\kappa } can be expressed as follows [³⁵]:

(4)\mu ={\mu }_{\mathrm{1}},\sigma =\sqrt{c_{\mathrm{2}}},\tau ={\displaystyle \frac{c_{\mathrm{3}}}{{\sigma }^{\mathrm{3}}}},\kappa ={\displaystyle \frac{c_{\mathrm{4}}}{{\sigma }^{\mathrm{4}}}}.

To obtain the preceding statistical moments of the system response, the numerical integration-based method can be used. When X=\xi is a one-dimensional random variable, the statistical moments can be solved by using the Gaussian integration rule as follows [³⁶]:

(5){\mu }_t={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t(X)}{f}_X(x)\mathrm{d}X={\displaystyle \sum _{j=\mathrm{1}}^n{w}_j}{g}^t({\xi }_j),\hspace{0.17em}\hspace{0.17em}t=\mathrm{1},\mathrm{2},\dots ,

(6)c_t={\displaystyle \underset{\mathrm{\Omega }}{\int }(g(X)}-{\mu }_{\mathrm{1}}{)}^t{f}_X(x)\mathrm{d}X={\displaystyle \sum _{j=\mathrm{1}}^n{w}_j}{(g({\xi }_j)-{\mu }_{\mathrm{1}})}^t,\hspace{0.17em}\hspace{0.17em}t=\mathrm{2},\mathrm{3},\dots ,

where {\xi }_j and w_j represent the j-th element in the i-th order nodes

U_{\mathrm{1}}^i

and the i-th order weights W_{\mathrm{1}}^i, respectively, and n is the total number of nodes U_{\mathrm{1}}^i. The nodes and weights can be obtained by using the quadrature rule. For example, the Gauss-Hermite quadrature rule can be used if \mathit{X} follows the normal distribution, whereas the Gauss-Legendre quadrature rule can be considered if \mathit{X} follows the uniform distribution.

When \mathit{X} is a multidimensional random vector, the multidimensional nodes {\mathit{U}}_n^k with the n-dimensional variables and k-level accuracy should be defined first. Furthermore, the SGNI provides the solution of the multidimensional nodes U_n^k [²⁵]:

(7){\mathit{U}}_n^k={\displaystyle \underset{k+\mathrm{1}\le \left|i\right|\le k+n}{\cup }{U}_{\mathrm{1}}^{i_{\mathrm{1}}}}\otimes U_{\mathrm{1}}^{i_{\mathrm{2}}}\otimes \mathrm{...}\otimes U_{\mathrm{1}}^{i_n},

where \left|i\right| (\left|i\right|=i_{\mathrm{1}}+i_{\mathrm{2}}+\mathrm{....}+i_n) denotes the summation of the multi-indices, i_j(j=\mathrm{1},\mathrm{2},\mathrm{...},n) is the univariate index in the quadrature rule, and \otimes is the tensor symbol, which is defined as follows: If R=\{a,b\}, S=\{\mathrm{0},\mathrm{1},\mathrm{2}\}, then R\otimes S=\{(a,\mathrm{0}),(a,\mathrm{1}),(a,\mathrm{2}),(b,\mathrm{0}),(b,\mathrm{1}),(b,\mathrm{2})\}. Therefore, the nodes based on Eq. (7) can be obtained by using the SGNI. For each obtained node, the corresponding weight is given as follows:

(8)w_l={(-\mathrm{1})}^{k+n-\left|i\right|}\left(\begin{array}{l}n-\mathrm{1}\\ k+n-\left|i\right|\end{array}\right)w_{\mathrm{1},j_{\mathrm{1}}}^{i_{\mathrm{1}}}{w}_{\mathrm{1},j_{\mathrm{2}}}^{i_{\mathrm{2}}}\mathrm{...}{w}_{\mathrm{1},j_n}^{i_n},

where w_{\mathrm{1},j}^i denotes the weight corresponding to the node {\xi }_{\mathrm{1},j}^i, {\xi }_{\mathrm{1},j}^i is the j-th element of the nodes U_{\mathrm{1}}^i, w_l denotes the weight corresponding to the l-th multidimensional node {\mathit{\xi }}_l, and {\mathit{\xi }}_l={[{\xi }_{\mathrm{1},j_{{}_{\mathrm{1}}}}^{i_{\mathrm{1}}},{\xi }_{\mathrm{1},j_{{}_{\mathrm{2}}}}^{i_{\mathrm{1}}},\dots ,{\xi }_{\mathrm{1},j_n}^{i_n}]}^T\in {\mathit{U}}_n^k. Thus, the statistical moments can be expressed as follows:

(9){\mu }_t={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t(\mathit{X})}{f}_X(x)\mathrm{d}\mathit{X}={\displaystyle \sum _{l=\mathrm{1}}^N{w}_l}{g}^t({\xi }_l),\hspace{0.17em}\hspace{0.17em}t=\mathrm{1},\mathrm{2},\dots ,

(10)c_t={\underset{\mathrm{\Omega }}{\int }(g(\mathit{X})-{\mu }_{\mathrm{1}})}^t{f}_X(x)\mathrm{d}\mathit{X}={\displaystyle \sum _{l=\mathrm{1}}^N{w}_l}{(g({\xi }_l)-{\mu }_{\mathrm{1}})}^t,\hspace{0.17em}\hspace{0.17em}t=\mathrm{2},\mathrm{3},\dots ,

where N stands for the number of all multidimensional nodes, which depends on the level k. When the level is specified as k=1 or 2, the relationship between the number N and the dimension n can be deduced as follows:

(11)\{\begin{array}{l}N=\mathrm{2}n+\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\mathrm{1}\\ N=\mathrm{2}{n}^{\mathrm{2}}+\mathrm{2}n+\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\mathrm{2}\end{array}.

A total of 221 multidimensional nodes exist when the dimension and level are set as n=2 and k=2, respectively. However, if the direct integration method is used to solve the system response statistical moments, then the number of multidimensional nodes will be {\mathrm{3}}^{\mathrm{10}}=\mathrm{59049} despite selecting only three nodes for each random variable.

The proposed method for UP analysis

Existing studies [²⁵,³⁷,³⁸] have reported that the SGNI is suitable for uncertainty analysis due to its good performance, especially in estimating the low order moments of the system response. However, most researchers are using the SGNI with the traditional univariate integration points, such as the Gauss-Hermite or Gauss-Legendre integration nodes. Therefore, guaranteeing the accuracy of the response function high-order moments is a challenging task. Certainly, increasing the level k could improve the accuracy but lead to a dramatic increase in computation cost. In addition, although some extended Gauss-Hermite univariate integration points can enhance the solving accuracy, it may be unsuitable for some input random variable distributions. Therefore, in this paper, the extended Gauss integration and Rosenblatt transformation are introduced into the UP analysis and the SGNI with the transformed nodes extended to estimate the system response higher-order moments. Thus, the proposed method not only guarantees the accuracy of the low-order moments but also ensures the solving accuracy of the high-order moments, based on which a fine accuracy thus can be ensured for fitting the final PDF of the response function.

Extended Gauss integration

To obtain a higher precision of the numerical integration, the extended Gauss integration is derived from the traditional Gauss integration, which can be expressed as follows [³⁹]:

(12){\displaystyle \underset{a}{\overset{b}{\int }}w(X)}y(X)dX={\displaystyle \sum _{i=\mathrm{1}}^n{A}_i}y(X_i)+{\displaystyle \sum _{j=\mathrm{1}}^m{A}_{j+n}y(X_{j+n})},

where y(X) is the integrand with integration interval [a,b], w(X) is the weight function, and X_i and A_i denote the integration nodes and weights, respectively.

Considering the various extension approaches, some researchers have studied the extended integration nodes. Kronrod [⁴⁰]and Patterson [⁴¹] first proposed a new method to handle the case with weight function w(X)=\mathrm{1} and the interval[a,b]=[-\mathrm{1},\mathrm{1}], which significantly improved the solving accuracy and computational efficiency. Additionally, Genz and Keister [⁴²] presented the fully symmetric interpolatory rules for EGHI, which considered the weight function and interval specified as w(V)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi }}}{e}^{-\frac{V^{\mathrm{2}}}{\mathrm{2}}} and [a,b]=[-\infty ,\infty ] in the standard normal space V, respectively. In fact, several EGHI nodes have already been explored by using different fully symmetric interpolatory rules. To ensure computational efficiency, the EGHI nodes can be selected in Table 1 in this paper. The integration nodes in each level contain the previous nodes in the EGHI, which effectively prevent repetitive computation. Therefore, its computational efficiency is comparable to that of the traditional Gauss-Hermite integration but with higher algebraic precision. Consequently, the solving accuracy for the following UP problems can be improved theoretically.

Tab.1 EGHI nodes and weights

Level	Integration node	Integration weight	Algebraic precision
1	V_{\mathrm{1}}^{\mathrm{1}}=\mathrm{0}	A_{\mathrm{1}}^{\mathrm{1}}=\mathrm{1}	1
2	V_{\mathrm{1}}^{\mathrm{2}}=\{-\mathrm{1.7321},\mathrm{0},\mathrm{1.7321}\}	A_{\mathrm{1}}^{\mathrm{2}}=\{\mathrm{0.1667},\mathrm{0.6667},\mathrm{0.1667}\}	5
3	V_{\mathrm{1}}^{\mathrm{3}}=\left\{\begin{array}{l}-\mathrm{2.8613},-\mathrm{0.7411},-\mathrm{4.1850},\\ -\mathrm{1.7321},\mathrm{0},\hspace{0.17em}\mathrm{1.7321},\\ \mathrm{4.1850},\mathrm{0.7411},\mathrm{2.8613}\end{array}\right\}	A_{\mathrm{1}}^{\mathrm{3}}=\left\{\begin{array}{l}\mathrm{0.0080},\mathrm{0.2701},\mathrm{0.0001},\\ \mathrm{0.0949},\mathrm{0.2540},\hspace{0.17em}\mathrm{0.0949},\\ \mathrm{0.0001},\mathrm{0.2701},\mathrm{0.0080}\end{array}\right\}	15
4	V_{\mathrm{1}}^{\mathrm{3}}=\left\{\begin{array}{l}-\mathrm{3.2053},-\mathrm{2.5961},-\mathrm{5.1870},-\mathrm{1.2304},\\ -\mathrm{6.3634},-\mathrm{2.8613},-\mathrm{0.7411},-\mathrm{4.1850},\\ -\mathrm{1.7321},\mathrm{0},\hspace{0.17em}\mathrm{1.7321},\mathrm{4.1850},\mathrm{0.7411},\\ \mathrm{2.8613},\mathrm{6.3634},\mathrm{1.2304},\mathrm{5.1870},\\ \mathrm{2.5961},\mathrm{3.2053}\end{array}\right\}	A_{\mathrm{1}}^{\mathrm{4}}=\left\{\begin{array}{l}\mathrm{0.0029},\mathrm{0.0181},\mathrm{0},\mathrm{0.0612}\\ \mathrm{0},\mathrm{0.0063},\mathrm{0.2083},\mathrm{0.0001}\\ \mathrm{0.0641},\mathrm{0.3035},\mathrm{0.0641},\mathrm{0.0001}\\ \mathrm{0.2083},\mathrm{0.0063},\mathrm{0},\mathrm{0.0612},\\ \mathrm{0},\mathrm{0.0181},\mathrm{0.0029}\end{array}\right\}	29
⋮	⋮	⋮	⋮

Transformation of extended integration nodes

In practice, although the extended Gauss integration rule has excellent computational precision, it is not suitable for cases with input random variable arbitrary distributions. For example, in the case of Weibull distribution or the student’s t-distribution, obtaining the extended integration nodes directly is challenging.

To solve these types of problems, corresponding variable transformation methods have been developed, such as the Rosenblatt transformation [⁴³,⁴⁴] and the work of Huang and Du [⁴⁵]. However, these variable transformations are based on the traditional Gauss integration. Therefore, the computational precision cannot be increased considerably.

Thus, the Rosenblatt transformation is employed into the EGHI in this paper because it has the same weight function as that of that traditional Gauss-Hermite integration. Therefore, this transformation can be applicable to any continuous distribution by using the extended integration nodes. The proposed transformation between a standard normal variable V and an original variable X can be expressed as follows:

(13)\mathbf{\Phi }(v_i)=F_{X_i}(x_i),

where \mathbf{\Phi }(\cdot ) and F_{X_i}(\cdot ) denote the CDF of the variable V and X, and v_i and x_i represent a realization of the random variable V and the standard normal variable X, respectively. Unlike the traditional transformation, v_i(i=\mathrm{1},\mathrm{2},\dots ,n+m) denotes an EGHI node, which can be selected in Table 1, and x_i denotes a node corresponding to the variable X. Therefore, Eq. (13) can be rewritten as

(14)v_i={\mathbf{\Phi }}^{-\mathrm{1}}(F_{X_i}(x_i))=T(x_i),

where {\mathbf{\Phi }}^{-\mathrm{1}}(\cdot ) is the inverse function of the standard normal variable CDF and T(\cdot ) is the function related to the variable X. Thus, the node x_i can be solved by using the following formula:

(15)x_i=T^{-\mathrm{1}}(v_i).

Furthermore, the system response statistical moments in Eqs. (2) and (3) can be replaced as follows:

(16){\mu }_t={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t(\mathit{X})}{f}_X(x)\mathrm{d}\mathit{X}={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t[T^{-\mathrm{1}}(\mathit{V})]}{f}_V(v)\mathrm{d}\mathit{V},\hspace{0.17em}\hspace{0.17em}t=\mathrm{1},\mathrm{2},\dots ,

(17)c_t={\underset{\mathrm{\Omega }}{\int }(g(\mathit{X})-{\mu }_{\mathrm{1}})}^t{f}_X(x)\mathrm{d}\mathit{X}=={\underset{\mathrm{\Omega }}{\int }\{g[T^{-\mathrm{1}}(\mathit{V})]-{\mu }_{\mathrm{1}}\}}^t{f}_V(v)\mathrm{d}\mathit{V},\hspace{0.17em}\hspace{0.17em}t=\mathrm{2},\mathrm{3},\dots ,

where T^{-\mathrm{1}}(\mathit{V})=[T^{-\mathrm{1}}(V_{\mathrm{1}}),T^{-\mathrm{1}}(V_{\mathrm{2}}),\dots ,T^{-\mathrm{1}}(V_n)]. Therefore, the moments in the UP problems can be calculated by using Eqs. (16) and (17).

Extended sparse grid technique

In this section, SGNI with the transformed extended nodes will be extended for the UP problems, denoting X_{\mathrm{1}}^i and A_{\mathrm{1}}^i with i-th level as the univariate nodes and weights of original variable \mathit{X}, respectively. The multidimensional nodes {\mathit{X}}_n^k with k-level accuracy and n-dimensional variables can be solved as follows:

(18){\mathit{X}}_n^k={\displaystyle \underset{k+\mathrm{1}\le \left|i\right|\le k+n}{\cup }{X}_{\mathrm{1}}^{i_{\mathrm{1}}}}\otimes X_{\mathrm{1}}^{i_{\mathrm{2}}}\otimes \dots \otimes X_{\mathrm{1}}^{i_n},

where \left|i\right|=i_{\mathrm{1}}+i_{\mathrm{2}}+\mathrm{....}+i_n. Based on Eq. (15), Eq. (18) can be rewritten as follows:

(19){\mathit{X}}_n^k={\displaystyle \underset{k+\mathrm{1}\le \left|i\right|\le k+n}{\cup }{T}^{-\mathrm{1}}(V_{\mathrm{1}}^{i_{\mathrm{1}}})}\otimes T^{-\mathrm{1}}(V_{\mathrm{1}}^{i_{\mathrm{2}}})\otimes \dots \otimes T^{-\mathrm{1}}(V_{\mathrm{1}}^{i_n}),

where i_j(j=\mathrm{1},\mathrm{2},\dots ,n) is the integration level of univariate nodes and V_{\mathrm{1}}^{i_j} is the nodes in the EGHI. Therefore, the multidimensional nodes corresponding to the arbitrary distributions of the input random variables can be obtained, assisted by the univariate nodes in EGHI.

To easily understand the multidimensional node formation, we take k=\mathrm{2} and n=\mathrm{2} as examples, which corresponds to the standard normal distribution. The generation of the multidimensional nodes X_{\mathrm{2}}^{\mathrm{2}} is shown in Table 2. Furthermore, Fig. 1 illustrates the construction of the nodes X_{\mathrm{2}}^{\mathrm{2}}. Under the condition of the standard normal distribution, if the accuracy level is k=\mathrm{1} and variable dimension is n, then \left|i\right|=n or \left|i\right|=n+\mathrm{1} satisfies the equation \mathrm{2}\le \left|i\right|\le \mathrm{1}+n. The former case \left|i\right|=n corresponds to a multidimensional node, namely, (\mathrm{0},\mathrm{0},\dots ,\mathrm{0}). In the latter case \left|i\right|=n+\mathrm{1}, the number of multidimensional nodes is 2n. For example, the multi-indices (i_{\mathrm{1}},i_{\mathrm{2}},\dots ,i_{n-\mathrm{1}},i_n)=(\mathrm{1},\mathrm{1},\dots ,\mathrm{1},\mathrm{2}) will generate two new nodes, namely, (\mathrm{0},\mathrm{0},\dots ,\mathrm{0},-\mathrm{1.732}) and (\mathrm{0},\mathrm{0},\dots ,\mathrm{0},\mathrm{1.732}). Similarly, when the accuracy level is k=2, we can also deduce the total number of generated nodes after removing the duplicated nodes. Therefore, the relationship between the number N_e and the dimension n can be defined as follows:

(20)\{\begin{array}{l}{N}_e=\mathrm{2}n+\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\mathrm{1}\\ N_e=\mathrm{2}{n}^{\mathrm{2}}+\mathrm{6}n+\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\mathrm{2}\end{array}.

For a generated multidimensional node, the corresponding weights can be obtained as follows [⁴⁶]:

(21)a_l={(-\mathrm{1})}^{k+n-\left|i\right|}\left(\begin{array}{l}n-\mathrm{1}\\ k+n-\left|i\right|\end{array}\right)a_{\mathrm{1},j_{\mathrm{1}}}^{i_{\mathrm{1}}}{a}_{\mathrm{1},j_{\mathrm{2}}}^{i_{\mathrm{2}}}\dots a_{\mathrm{1},j_n}^{i_n},

where a_l denotes the weight corresponding to the l-th multidimensional node {\mathit{\eta }}_l and {\mathit{\eta }}_l={[{\eta }_{\mathrm{1},j_{{}_{\mathrm{1}}}}^{i_{\mathrm{1}}},{\eta }_{\mathrm{1},j_{{}_{\mathrm{2}}}}^{i_{\mathrm{1}}},\mathrm{...},{\eta }_{\mathrm{1},j_n}^{i_n}]}^T\in {\mathit{X}}_n^k, and a_{\mathrm{1},j}^i denotes the weight corresponding to the one-dimensional node {\eta }_{\mathrm{1},j}^i.

For the extended sparse grid technique, the same multidimensional nodes will be combined multiple times. Instead of evaluating the response function several times, the same nodes are calculated only once and the respective weights should be summed up beforehand. Further details are available in the work of Heiss and Winschel [⁴⁷]. Therefore, the statistical moments of the system response in Eqs. (16) and (17) can be addressed by the following formulas:

(22){\mu }_t={\displaystyle \underset{\mathrm{\Omega }}{\int }{g}^t[T^{-\mathrm{1}}(\mathit{V})]}{f}_V(v)\mathrm{d}\mathit{V}=={\displaystyle \sum _{l=\mathrm{1}}^{N_e}{a}_l}{g}^t({\mathit{\eta }}_l),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t=\mathrm{1},\mathrm{2},\dots ,

(23)c_t={\underset{\mathrm{\Omega }}{\int }\{g[T^{-\mathrm{1}}(\mathit{V})]-{\mu }_{\mathrm{1}}\}}^t{f}_V(v)\mathrm{d}\mathit{V}=={\displaystyle \sum _{l=\mathrm{1}}^{N_e}{a}_l}{(g({\mathit{\eta }}_l)-{\mu }_{\mathrm{1}})}^t,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t=2,3,\dots .

Thereafter, the aforementioned moments can be transformed into those in Eq. (4).

Tab.2 Information on multidimensional nodes {\mathit{X}}_{\mathrm{2}}^{\mathrm{2}}

\left|i\right|	i_{\mathrm{1}}	i_{\mathrm{2}}	X_{\mathrm{1}}^{i_{\mathrm{1}}}\otimes X_{\mathrm{1}}^{i_{\mathrm{2}}}	{\mathit{X}}_{\mathrm{2}}^{\mathrm{2}}
3	1	2	\{(\mathrm{0},-\mathrm{1.732}),(\mathrm{0},\mathrm{0}),(\mathrm{0},\mathrm{1.732})\}	\left\{\begin{array}{l}(\mathrm{0},-\mathrm{1.732}),(\mathrm{0},\mathrm{0}),(\mathrm{0},\mathrm{1.732})\\ (-\mathrm{1.732},\mathrm{0}),(\mathrm{1.732},\mathrm{0}),(\mathrm{0},-\mathrm{2.861}),\\ (\mathrm{0},-\mathrm{0.741}),(\mathrm{0},-\mathrm{4.185}),(\mathrm{0},\mathrm{2.861}),\\ (\mathrm{0},\mathrm{0.741}),(\mathrm{0},\mathrm{4.185}),\\ (-\mathrm{1.732},-\mathrm{1.732}),(-\mathrm{1.732},\mathrm{1.732}),\\ (\mathrm{1.732},-\mathrm{1.732}),(\mathrm{1.732},\mathrm{1.732}),\\ (-\mathrm{2.861},\mathrm{0}),(-\mathrm{0.741},\mathrm{0}),(-\mathrm{4.185},\mathrm{0})\\ (\mathrm{2.861},\mathrm{0}),(\mathrm{0.741},\mathrm{0}),(\mathrm{4.185},\mathrm{0})\end{array}\right\}
	2	1	\{(-\mathrm{1.732},\mathrm{0}),(\mathrm{0},\mathrm{0}),(\mathrm{1.732},\mathrm{0})\}
4	1	3	\left\{\begin{array}{l}(\mathrm{0},-\mathrm{2.861}),(\mathrm{0},-\mathrm{0.741}),(\mathrm{0},-\mathrm{4.185}),\\ (\mathrm{0},-\mathrm{1.732}),(\mathrm{0},\mathrm{0}),(\mathrm{0},\mathrm{1.732}),\\ (\mathrm{0},\mathrm{2.861}),(\mathrm{0},\mathrm{0.741}),(\mathrm{0},\mathrm{4.185})\end{array}\right\}
	2	2	\left\{\begin{array}{l}(-\mathrm{1.732},-\mathrm{1.732}),(-\mathrm{1.732},\mathrm{0}),(\mathrm{0},\mathrm{0}),\\ (\mathrm{0},-\mathrm{1.732}),(-\mathrm{1.732},\mathrm{1.732}),(\mathrm{0},\mathrm{1.732}),\\ (\mathrm{1.732},-\mathrm{1.732}),(\mathrm{1.732},\mathrm{0}),(\mathrm{1.732},\mathrm{1.732})\end{array}\right\}
	3	1	\left\{\begin{array}{l}(-\mathrm{2.861},\mathrm{0}),(-\mathrm{0.741},\mathrm{0}),(-\mathrm{4.185},\mathrm{0}),\\ (-\mathrm{1.732},\mathrm{0}),(\mathrm{0},\mathrm{0}),(\mathrm{1.732},\mathrm{0}),\\ (\mathrm{2.861},\mathrm{0}),(\mathrm{0.741},\mathrm{0}),(\mathrm{4.185},\mathrm{0})\end{array}\right\}

Fig.1 Construction of multidimensional nodes {\mathit{X}}_{\mathrm{2}}^{\mathrm{2}}

Full size|PPT slide

Maximum entropy principle

After the statistical moments are obtained, the maximum entropy principle [⁷,⁴⁸,⁴⁹] can be used to fit the PDF f_Y(y) of the system response. Given the constraints of the central moments, the formulation can be expressed as follows:

(24)\begin{array}{l}\max \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}H=-{\displaystyle \underset{\Omega }{\int }{f}_Y(y)}\ln f_Y(y)\mathrm{d}Y\\ \hspace{0.17em}\hspace{0.17em}s\mathrm{.t}.\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\{\begin{array}{l}\hspace{0.17em}{\displaystyle \underset{\mathrm{\Omega }}{\int }{f}_Y(y)}\mathrm{d}Y=\mathrm{1}\\ \hspace{0.17em}{\displaystyle \underset{\mathrm{\Omega }}{\int }{(Y-\mu )}^i{f}_Y(y)}\mathrm{d}Y=c_i,i=\mathrm{1},\mathrm{2},\dots ,r\end{array}\end{array}

where H denotes the entropy, \mu is the mean of the system response Y, and c_i stands for the i-th order central moment of Y when i=\mathrm{1}, c_{\mathrm{1}}=\mathrm{1}.

The optimization problem in Eq. (24) can be solved by using the Lagrangian multiplier method with a function formulated as

(25)L=-{\displaystyle \underset{\mathrm{\Omega }}{\int }{f}_Y(y)}\ln f_Y(y)\mathrm{d}Y+({\lambda }_{\mathrm{0}}+\mathrm{1})\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{f}_Y(y)\mathrm{d}Y-\mathrm{1}}\right)+{\displaystyle \sum _{i=\mathrm{1}}^r{\lambda }_i\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(Y-\mu )}^i{f}_Y(y)\mathrm{d}Y-c_i}\right)},

where {\lambda }_i(i=\mathrm{0},\mathrm{1},\dots ,r) is the Lagrangian multiplier to be determined.

Let the derivations of L with respect to f_Y(y)dY be equal to zero. Therefore, the maximum entropy PDF can be obtained as follows:

(26)f_Y(y)=\exp \left(-{\displaystyle \sum _{i=\mathrm{0}}^r{\lambda }_i{(Y-\mu )}^i}\right).

By substituting Eq. (26) with the constraint conditions in Eq. (24), the determined Lagrangian multiplier can be obtained according to the gradient iteration. In general, the first four-order central moments are the constraint conditions. To obtain a higher fitting accuracy, a high order of the central moments can be specified. Further details are provided by Mohammad-Djafari [⁵⁰].

Computational procedure

In summary, the flowchart of the proposed method is provided in Fig. 2 and summarized in the following steps:

Step 1: The performance function g is established and the distribution function of the input random variables X is determined.

Step 2: Based on the extended Gauss-Hermite nodes in Table 1, the transformed nodes are calculated according to Eq. (15).

Step 3: The accuracy level k is specified and all multi-indices i=(i_{\mathrm{1}},i_{\mathrm{2}},\mathrm{...},i_n) are calculated, and then the multidimensional nodes {\mathit{X}}_n^k is constructed through Eq. (19).

Step 4: The weight v_l corresponding to each multi-dimensional node is calculated through Eq. (21).

Step 5: The first four-order statistical moments of the system response are calculated through Eqs. (22) and (23), and the moments are transformed through Eq. (4).

Step 6: The PDF of the system response is fit with the maximum entropy principle.

Fig.2 Computational flowchart of proposed method

Full size|PPT slide

Numerical examples and discussions

Three examples are investigated in this section. The first example considers a mathematical problem with all normal input random variable distributions, and the second example extends the distribution to the five types with 10 dimensional variables. Some cases for the first two numerical examples are discussed, which aims to show a higher accuracy and an acceptable efficiency of the proposed method. The third example is an engineering application that concerns a signal transmission problem in the power amplifier link. The UP problems for each example are analyzed using the MCS, UDRM, SGNI, and proposed methods. Furthermore, the MCS solution is used as a reference to verify the accuracy of the three other methods.

Numerical Example 1

The first example considers a performance function with n-dimensional independent input random variables:

(27)Y=g(\mathit{X})={\displaystyle \sum _{i=\mathrm{1}}^n{X}_i}+\mathrm{20}{X}_{\mathrm{1}}{}^{\mathrm{2}}{X}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \sum _{i=\mathrm{1}}^{n-\mathrm{2}}{X}_{i+\mathrm{1}}^{\mathrm{2}}{X}_{i+\mathrm{2}}^{\mathrm{2}}}-{\displaystyle \sum _{i=\mathrm{1}}^n\sin (X_i)}\exp (X_i-\mathrm{2})-\mathrm{10},

where X_i(i=\mathrm{1},\mathrm{2},\mathrm{...},n) denotes the input random variable, which follows the normal distribution N(\mathrm{1},\mathrm{0.1}). Three cases will be considered in this example. The first case discusses the accuracy of the aforementioned methods with the same accuracy level k in the SGNI and the proposed methods. The second case mainly analyzes the method efficiency when the SGNI and the proposed method achieve nearly the same accuracy. The third case considers the UP problems with different dimensions n.

Case 1: Accuracy comparison

In this case, the input random variable dimension is defined as n=\mathrm{6}. Moreover, the aforementioned methods are adopted to estimate the first four-order moments of the system response. The analysis results are listed in Table 3. First, the errors of the first two-order moments are relatively small regardless of which method is adopted, but the accuracy with the proposed method is slightly better than that with either UDRM or SGNI. The maximum error of the low-order moments is only 1.83%, which occurs in the standard deviation \sigma calculation by UDRM. Second, in terms of the higher-order moments, the accuracy for each method is different. The proposed method can ensure much better accuracy compared to that of UDRM and SGNI. For example, the errors in predicting the skewness \tau by the UDRM, SGNI, and proposed methods are 66.45%, 5.53%, and 0.24%, respectively. The respective errors in predicting the kurtosis \kappa are 14.92%, 19.16%, and 0.99%, respectively. Furthermore, the moments obtained are used to fit the PDF of the system response with the maximum entropy principle. As depicted in Fig. 3(a), the system response PDFs by using the UDRM and SGNI are relatively inaccurate, whereas the proposed method can fit the PDF of MCS effectively, which can further demonstrate the accuracy of the proposed method.

Tab.3 Results of first four-order moments of Case 1 in Example 1

Moments	MCS	UDRM (error)	SGNI (error)	Proposed method (error)
\mathit{\mu }	18.6192	18.6108 (0.05%)	18.6132 (0.03%)	18.6132 (0.030%)
\mathit{\sigma }	6.1305	6.0184 (1.83%)	6.1295 (0.02%)	6.1303 (0.002%)
\mathit{\tau }	0.5976	0.2005 (66.45%)	0.5645 (5.53%)	0.5961 (0.240%)
\mathit{\kappa }	3.5904	3.0546 (14.92%)	2.9026 (19.16%)	3.5548 (0.990%)

Fig.3 PDF of system response in Example 1. (a) PDF of Case 1; (b) PDF of Case 2; (c) PDF of Case 3

Full size|PPT slide

In addition, for the computational efficiency, the number of function calls of the MCS is 1000000, which is more than those of the other three methods. Compared with SGNI and the proposed method, the computational efficiency of the UDRM is relatively high, whereas the accuracy of the high-order moments obtained by using the UDRM is unsatisfactory, which results in a relatively poor fitting accuracy of the system response PDF. The SGNI needs 85 (k=2) function calls, which is less than that of the proposed method (109). Despite a slightly lower computational efficiency, the proposed method can ensure better high-order moment accuracy. Overall, the proposed method provides a comprehensive result between the solving accuracy and computational efficiency in solving the UP problems, which is significant for the uncertainty analysis.

Case 2: Efficiency comparison

As mentioned in Case 1, the proposed method can have higher accuracy than that of SGNI under the same level k. In this case, the efficiencies of the SGNI and proposed methods are discussed under a nearly similar accuracy. For comparison, the moment estimate solving accuracy is improved by increasing the accuracy level k until the results of the SGNI and proposed methods are substantially close to that of the MCS. The UDRM is not used for comparison because it cannot obtain a high-order moment accurate solution for this problem.

The results are presented in Table 4. In this case, for the various levels k, the accuracy of all statistical moments by using the MCS, SGNI, and proposed method is nearly the same. However, MCS requires 1000000 function calls to achieve an exact solution. The SGNI can obtain good accuracy with 4541 function evaluations. At this point, it corresponds to the level of k=\mathrm{5}, whereas the proposed method with the level k=\mathrm{3} only needs 689 function calls to obtain the same accuracy, which is nearly 6.5 times less than the number of the SGNI. The proposed method can have the same accuracy but with a higher efficiency because the univariate nodes has a higher algebraic precision. For example, the Level 3 algebraic precision in Table 1 achieves a value of 15, whereas the Level 3 traditional integration only achieves a relatively smaller algebraic precision, that is, 7 (\mathrm{2}\times \mathrm{3}+\mathrm{1}). Therefore, the proposed method is more efficient than the SGNI when they possess the same accuracy, which is important in solving vague practical problems. Moreover, the system response PDF is fitted based on the calculated moments in Fig. 3(b). The PDF, by using the SGNI and proposed methods, can be fitted well with that of the MCS due to their similar accuracy.

Tab.4 Results of first four-order moments of Case 2 in Example 1

Moments	MCS	SGNI (error)	Proposed method (error)
\mathit{\mu }	18.6192	18.6132 (0.030%)	18.6132 (0.030%)
\mathit{\sigma }	6.1305	6.1304 (0.001%)	6.1304 (0.001%)
\mathit{\tau }	0.5976	0.5978 (0.040%)	0.5978 (0.040%)
\mathit{\kappa }	3.5904	3.5964 (0.170%)	3.5962 (0.160%)

Case 3: Dimension comparison

To analyze the accuracy tendency of each order moment with the increasing dimension, the different n values are selected in this case. The minimum dimension is set as n=\mathrm{6} because the proposed method can handle the problem with high dimensional variables. Furthermore, the maximum size is set as n=\mathrm{30} due to its suitable efficiency.

The moments obtained from the UDRM, SGNI, and proposed methods are compared with the results from that of the MCS, which shows that the relative errors varied with the dimensions n plotted in Fig. 4. The relative errors for each order moments fluctuate slightly as the dimension n increases. Similar to Case 1, the aforementioned methods can guarantee the accuracy of the first two-order moments, with maximal relative errors of only 0.17% and 2.09%, respectively, which all occur in the UDRM calculation. However, neither the UDRM nor the SGNI will obtain satisfactory accuracy for the UP problems due to larger relative errors in predicting the high-order moments. For example, the maximal error of skewness by using UDRM is 67.66%, whereas the maximal error of kurtosis by using SGNI reaches 22.3%. By contrast, the high-order moment errors obtained by using the proposed method are only 4.84% and 4.93%, respectively, which are less than those of the UDRM and SGNI methods. As a result, the proposed method can guarantee the accuracy not only of the low-order moments but also of the high-order moments, which can ensure the fitting accuracy of the system response PDF. Furthermore, the obtained moments, which correspond to dimension n=\mathrm{27}, are adopted to fit the system response PDF, as shown in Fig. 3(c). The proposed method provides a more accurate system response PDF than either that of the UDRM or the SGNI, which is similar to Case 1.

Fig.4 The relative error of the first-order moments varies with the dimension n. The relative error of (a) the mean, (b) the standard deviation, (c) the skewness, and (d) the kurtosis varies with the dimension n

Full size|PPT slide

Numerical Example 2

This example is adopted to analyze the UP problems with several different distributions of input random variables with the response function created as follows:

(28)Y=g(\mathit{X})={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{10}}{X}_i}+X_{\mathrm{1}}{X}_{\mathrm{2}}+{\displaystyle \frac{\exp (-X_{\mathrm{3}})}{X_{\mathrm{4}}}}\mathrm{+sin}{X}_{\mathrm{5}}\cos X_{\mathrm{6}}+X_{\mathrm{7}}{X}_{\mathrm{8}}+\mathrm{20}{X}_{\mathrm{9}}^{\mathrm{2}}{X}_{\mathrm{10}}^{\mathrm{2}}-\mathrm{30.}

The random variable X_i distribution is given in Table 5. Two cases will be tested for the problem: One with the same level k and the other with the same accuracy, which are both considered in the SGNI and the proposed method.

Tab.5 Distributions of input random variable in Example 2

Variables	Distribution	Parameter 1	Parameter 2
X_{\mathrm{1}}	Normal	1	0.12
X_{\mathrm{2}}	Normal	5	0.5
X_{\mathrm{3}}	Weibull	1	5
X_{\mathrm{4}}	Uniform	2	6
X_{\mathrm{5}}	Lognormal	2	0.2
X_{\mathrm{6}}	Beta	2	5
X_{\mathrm{7}}	Normal	1	0.12
X_{\mathrm{8}}	Normal	1	0.12
X_{\mathrm{9}}	Normal	1	0.12
X_{\mathrm{10}}	Normal	1	0.12

Note: Parameter i (i=1, 2) denotes the i-th parameter of input random variable

Case 1: Accuracy comparison

Different from Case 1 in the first example, this case mainly discusses the problem with five different input distributions. Therefore, the transformation is involved in the UP analysis. However, without additional loss of accuracy caused by the non-normal distribution of the input random variable, the trends are almost similar to those in Example 1, as shown in Table 6. Each method can accurately estimate the system response low-order moments. For example, the standard deviation errors by using the UDRM, SGNI, and proposed methods are only 3.86%, 0.69% and 0.04%, respectively. However, with respect to a higher-order moment, such as skewness \tau, the results of these three methods are 64.04%, 27.75%, and 1.65%. Therefore, the proposed method is more accurate than the other methods. This phenomenon also occurred in the kurtosis \kappa calculation. The PDFs are then solved based on the obtained moments, as shown in Fig. 5(a). Thus, the moment estimation errors result in probability calculation errors. In general, the smaller the estimated statistical moment errors, the higher the PDF fitting accuracy. Therefore, the PDF by using the proposed method is close to that of the MCS.

In Table 6, the MCS has relatively low efficiency compared with the other methods. The numbers of the function calls of the MCS, UDRM, SGNI, and proposed methods are 1000000, 51, 221 (k=2), and 261 (k=2), respectively. Similar to Example 1, although the UDRM and SGNI have a relatively good computational efficiency, the accuracy of the high-order moments by using the two methods is unsatisfactory. On the contrary, the proposed method with a few function calls can ensure better accuracy than that of UDRM and SGNI.

Tab.6 Results of first four-order moments of Case 1 in Example 2

Moments	MCS	UDRM (error)	SGNI (error)	Proposed method (error)
\mathit{\mu }	19.5144	20.5589 (5.35%)	19.5616 (0.24%)	19.5133 (0.01%)
\mathit{\sigma }	7.6417	7.3465 (3.86%)	7.5890 (0.69%)	7.6449 (0.04%)
\mathit{\tau }	0.6202	0.2231 (64.04%)	0.4481 (27.75%)	0.6305 (1.65%)
\mathit{\kappa }	3.7232	3.0687(17.58%)	3.2924 (11.57%)	3.7960 (1.96%)

Fig.5 PDF of system response in Example 2. (a) PDF of Case 1; (b) PDF of Case 2

Full size|PPT slide

Case 2: Efficiency comparison

Similar to Example 1, the objective of this case is to compare the efficiency of the SGNI and proposed methods when they can obtain good accuracy of the moment estimates. As shown in Table 7, the errors for each moment by using the SGNI are 0.01%, 0.004%, 2.24%, and 0.04%, respectively, whereas the errors by using the proposed method are 0.01%, 0.002%, 1.73%, and 0.06%, respectively. Therefore, improving the level k, SGNI, and proposed method can result in a relatively accurate solution with k=\mathrm{7} and k=\mathrm{3}. However, the function calls of the two methods are different. For example, the SGNI requires up to 581385 function calls, but the proposed method with 2401 samples can also achieve the same accuracy. The proposed method function calls with the same k=\mathrm{3} in the first example are less than that in this case. The trend is normal due to a higher dimension in this example. The PDF results are shown in Fig. 5(b). All the PDF by using the proposed method are accurate compared with that of the MCS due to the minimal moment calculation error.

Tab.7 Results of first four-order moments of Case 2 in Example 2

Moments	MCS	SGNI (error)	Proposed method (error)
\mathit{\mu }	19.5144	19.5161 (0.01%)	19.5129 (0.01%)
\mathit{\sigma }	7.6417	7.6420 (0.004%)	7.6418 (0.002%)
\mathit{\tau }	0.6202	0.6063 (2.24%)	0.6095 (1.73%)
\mathit{\kappa }	3.7232	3.7219 (0.04%)	3.7256 (0.06%)

An engineering application

The transmit/receive (T/R) modules have attracted significant attention in recent years because of its importance in the practical target detection [⁵¹,⁵²]. The power amplifier link is an important part of the T/R modules, which mainly consists of two parts: Microstrip and printed board. In general, the printed board is divided into the FR4 printed board and microwave printed board. This example mainly considers a signal transmission problem in the power amplifier link. Furthermore, the finite element model is established as shown in Fig. 6.

Fig.6 Finite element model of power amplifier link

Full size|PPT slide

Let the input electromagnetic signal in the amplifier link be I_{\mathrm{in}}=A(\cos \mathrm{2}\pi f+{\phi }_{\mathrm{in}}) and the output be I_{\mathrm{out}}=n_{e}A(\cos \mathrm{2}\pi f+{\phi }_{\mathrm{out}}), where A denotes the amplitude of the input signal, n_e is the amplification factor, f represents the operating frequency, and {\phi }_{\mathrm{in}} and {\phi }_{\mathrm{out}} stand for the input and output phases, respectively. In practical application, we are more concerned about the uncertainty of the phase difference \Delta \phi (\Delta \phi ={\phi }_{\mathrm{out}}-{\phi }_{\mathrm{in}}). Therefore, the performance function can be expressed as follows:

(29)\Delta \phi =g(X_{\mathrm{1}},X_{\mathrm{2}},\dots ,X_n),

where X_i(i=\mathrm{1},\mathrm{2},\dots ,n) is the input variable. Fourteen variables will be considered in this example. The variables X_{\mathrm{1}}, X_{\mathrm{2}}, and X_{\mathrm{3}} denote the dielectric constant of the FR4 printed board, thickness of the FR4 printed board, and width of the microstrip in the FR4 printed board, respectively. The variables X_{\mathrm{4}}, X_{\mathrm{5}}, and X_{\mathrm{6}} denote the dielectric constant, thickness, and width of the microstrip in the microwave printed board, respectively. The variables from X_{\mathrm{7}} to X_{\mathrm{14}} denote the position sizes, which are shown in Fig. 6. The aforementioned variables are assumed to follow the different distributions due to existing uncertainty factors. The specific information is provided in Table 8.

The first four-order moments of the phase difference are obtained by using the aforementioned methods, as shown in Table 9. The errors in solving the low-order moments are completely accepted. For example, all methods provide a good estimate of the mean \mu, but the proposed method provides improved accuracy of the standard deviation \sigma. With respect to the higher-order moments, the errors in predicting the skewness \tau by using the UDRM, SGNI, and proposed methods are 74.71%, 10.97%, and 1.17%, respectively, whereas the errors in predicting the kurtosis \kappa are 24.17%, 26.27%, and 4.06%, respectively. Therefore, the proposed method provides a reasonably good estimate of the high-order moments. Figure 7 shows the phase difference PDF based on the obtained moments. The PDF from the proposed method is highly similar to that of the MCS, which can verify the accuracy of the proposed method.

Tab.8 Distributions of input random variables in Example 3

Variables	Distribution	Parameter 1/mm	Parameter 2
X_{\mathrm{1}}	Normal	4.80	0.033
X_{\mathrm{2}}	Normal	0.70	0.001
X_{\mathrm{3}}	Normal	0.90	0.017
X_{\mathrm{4}}	Normal	3.50	0.017
X_{\mathrm{5}}	Normal	1.00	0.001
X_{\mathrm{6}}	Normal	2.20	0.033
X_{\mathrm{7}}	Normal	12.90	0.017
X_{\mathrm{8}}	Normal	45.70	0.017
X_{\mathrm{9}}	Normal	96.15	0.017
X_{\mathrm{10}}	Normal	73.70	0.017
X_{\mathrm{11}}	Uniform	114.00	115 mm
X_{\mathrm{12}}	Uniform	129.10	129.5 mm
X_{\mathrm{13}}	Normal	133.40	0.017
X_{\mathrm{14}}	Normal	147.10	0.017

Note: Parameter i (i=1, 2) denotes the i-th parameter of input random variable

Tab.9 First four-order moments of phase difference in Example 3

Moments	MCS	UDRM (error)	SGNI (error)	Proposed method (error)
\mathit{\mu }	−45.0458	−45.3185 (0.60%)	−45.0591 (0.03%)	−45.0591 (0.03%)
\mathit{\sigma }	2.6849	2.5567 (4.77%)	2.6809 (0.15%)	2.6856 (0.03%)
\mathit{\tau }	0.8739	0.2210 (74.71%)	0.7780 (10.97%)	0.8637 (1.17%)
\mathit{\kappa }	4.6578	3.5320(24.17%)	3.4343 (26.27%)	4.4685 (4.06%)

Fig.7 PDF of phase difference in Example 3

Full size|PPT slide

On the one hand, through the study of the engineering application, the proposed method has higher accuracy in solving the UP problems. On the other hand, for the computational efficiency, the total number of function calls of the proposed method is 477 (k=2), which is acceptable for a practical engineering problem. The UDRM and SGNI have relatively good computational efficiency but neither methods can ensure the accuracy of the high-order moments, which will influence the practicability of the methods. Thus, the method proposed in this paper has a good comprehensive effect in practical engineering problems due to high accuracy and accepted efficiency.

Conclusions

This paper proposes an approach for uncertainty propagation analysis based on an extended sparse grid technique and the maximum entropy method, which satisfies both good accuracy and acceptable computational efficiency. Different from the traditional Rosenblatt transformation, the method proposed in this paper is extended to the extended Gauss integration. Therefore, all arbitrary variable integration nodes can be transformed into the EGHI nodes. Furthermore, the SGNI with the transformed nodes is employed to predict the first four-order moments. The system response PDF is obtained based on maximum entropy. Without loss of accuracy in the proposed transformation, our method provides relatively accurate estimates of the statistical moments, thereby ensuring a good fitting accuracy of the system response PDF. The results of the numerical examples also indicate that the proposed method is more accurate than the UDRM and SGNI methods. In terms of the first two-order moments, the accuracy of the proposed method is slightly better than that of UDRM and SGNI, but all methods are capable of ensuring the accuracy of the low-order moments. However, with respect to the high-order moments, only the proposed method has a good effect. In addition, the system response PDF has higher accuracy. When the performance function contains stronger interactions among the input random variables, the effectiveness of the proposed method improves. In addition, for the computational cost, the proposed method has relatively high efficiency compared with the MCS. All numerical examples are analyzed with the dependent variables, and the correlations among the input random variables can be considered in the future. Moreover, the proposed method extends to other interesting topics such as structural reliability analysis and optimization design.

Acknowledgments

This work was supported by the National Science Fund for Distinguished Young Scholars (Grant No. 51725502), the major program of the National Natural Science Foundation of China (Grant No. 51490662), and the National Key Research and Development Project of China (Grant No. 2016YFD0701105).

Open Access

This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the appropriate credit is given to the original author(s) and the source, and a link is provided to the Creative Commons license, indicating if changes were made.

References

1 Lee S H, Chen W. A comparative study of uncertainty propagation methods for black-box-type problems. Structural and Multidisciplinary Optimization, 2009, 37(3): 239–253

2 Wang X, Wang L, Qiu Z. Response analysis based on smallest interval-set of parameters for structures with uncertainty. Applied Mathematics and Mechanics, 2012, 33(9): 1153–1166

3 Wang X, Wang L, Qiu Z. A feasible implementation procedure for interval analysis method from measurement data. Applied Mathematical Modelling, 2014, 38(9–10): 2377–2397

4 Qiu Z P, Wang L. The need for introduction of non-probabilistic interval conceptions into structural analysis and design. Science China. Physics, Mechanics & Astronomy, 2016, 59(11): 114632

5 Gu X, Renaud J E, Batill S M, Worst case propagated uncertainty of multidisciplinary systems in robust design optimization. Structural and Multidisciplinary Optimization, 2000, 20(3): 190–213

6 Li M, Azarm S. Multiobjective collaborative robust optimization with interval uncertainty and interdisciplinary uncertainty propagation. Journal of Mechanical Design, 2008, 130(8): 081402

7 Li G, Zhang K. A combined reliability analysis approach with dimension reduction method and maximum entropy method. Structural and Multidisciplinary Optimization, 2011, 43(1): 121–134

8 Jiang Z, Li W, Apley D W, A spatial-random-process based multidisciplinary system uncertainty propagation approach with model uncertainty. Journal of Mechanical Design, 2015, 137(10): 101402

9 Mazo J, El Badry A T, Carreras J, Uncertainty propagation and sensitivity analysis of thermo-physical properties of phase change materials (PCM) in the energy demand calculations of a test cell with passive latent thermal storage. Applied Thermal Engineering, 2015, 90: 596–608

10 Li M, Mahadevan S, Missoum S, Special issue: Simulation-based design under uncertainty. Journal of Mechanical Design, 2016, 138(11): 110301

11 Madsen H O, Krenk S, Lind N C. Methods of structural safety. Mineola: Dover Publications, 2006

12 Wilson B M, Smith B L. Taylor-series and Monte-Carlo-method uncertainty estimation of the width of a probability distribution based on varying bias and random error. Measurement Science & Technology, 2013, 24(3): 035301

13 Rochman D, Zwermann W, van der Marck S C, Efficient use of Monte Carlo: Uncertainty propagation. Nuclear Science and Engineering, 2014, 177(3): 337–349

14 Hong J, Shaked S, Rosenbaum R K, Analytical uncertainty propagation in life cycle inventory and impact assessment: Application to an automobile front panel. International Journal of Life Cycle Assessment, 2010, 15(5): 499–510

15 Xu L. A proportional differential control method for a time-delay system using the Taylor expansion approximation. Applied Mathematics and Computation, 2014, 236: 391–399

16 Low B. FORM, SORM, and spatial modeling in geotechnical engineering. Structural Safety, 2014, 49: 56–64

17 Lim J, Lee B, Lee I. Post optimization for accurate and efficient reliability-based design optimization using second-order reliability method based on importance sampling and its stochastic sensitivity analysis. International Journal for Numerical Methods in Engineering, 2015, 107(2): 93–108

18 Lee I, Choi K K, Gorsich D. Sensitivity analyses of FORM-based and DRM-based performance measure approach (PMA) for reliability-based design optimization (RBDO). International Journal for Numerical Methods in Engineering, 2010, 82(1): 26–46

19 Sudret B. Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, 2008, 93(7): 964–979

20 Kersaudy P, Sudret B, Varsier N, A new surrogate modeling technique combining Kriging and polynomial chaos expansions—Application to uncertainty analysis in computational dosimetry. Journal of Computational Physics, 2015, 286: 103–117

21 Rajabi M M, Ataie-Ashtiani B, Simmons C T. Polynomial chaos expansions for uncertainty propagation and moment independent sensitivity analysis of seawater intrusion simulations. Journal of Hydrology (Amsterdam), 2015, 520: 101–122

22 Rahman S, Xu H. A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics. Probabilistic Engineering Mechanics, 2004, 19(4): 393–408

23 Xu H, Rahman S. A generalized dimension-reduction method for multidimensional integration in stochastic mechanics. International Journal for Numerical Methods in Engineering, 2004, 61(12): 1992–2019

24 Nobile F, Tempone R, Webster C G. A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM Journal on Numerical Analysis, 2008, 46(5): 2309–2345

25 Xiong F, Greene S, Chen W, A new sparse grid based method for uncertainty propagation. Structural and Multidisciplinary Optimization, 2010, 41(3): 335–349

26 He J, Gao S, Gong J. A sparse grid stochastic collocation method for structural reliability analysis. Structural Safety, 2014, 51: 29–34

27 Smolyak S A. Quadrature and interpolation formulas for tensor products of certain classes of functions. Doklady Akademii Nauk SSSR, 1963, 4: 240–243

28 Novak E, Ritter K. High dimensional integration of smooth functions over cubes. Numerische Mathematik, 1996, 75(1): 79–97

29 Novak E, Ritter K. Simple cubature formulas with high polynomial exactness. Constructive Approximation, 1999, 15(4): 499–522

30 Bathe K, Wilson E. Stability and accuracy analysis of direct integration methods. Earthquake Engineering & Structural Dynamics, 1972, 1(3): 283–291

31 Tao J, Zeng X, Cai W, Stochastic sparse-grid collocation algorithm (SSCA) for periodic steady-state analysis of nonlinear system with process variations. In: Proceedings of the 2007 Asia and South Pacific Design Automation Conference. IEEE, 2007, 474–479

32 Jia B, Xin M, Cheng Y. Sparse Gauss-Hermite quadrature filter for spacecraft attitude estimation. In: Proceedings of the 2010 American Control Conference. Baltimore: IEEE, 2010, 2873–2878

33 Petvipusit K R, Elsheikh A H, Laforce T C, Robust optimisation of CO₂ sequestration strategies under geological uncertainty using adaptive sparse grid surrogates. Computational Geosciences, 2014, 18(5): 763–778

34 Chen H, Cheng X, Dai C, Accuracy, efficiency and stability analysis of sparse-grid quadrature Kalman filter in near space hypersonic vehicles. In: Proceedings of Position, Location and Navigation Symposium-PLANS 2014, 2014 IEEE/ION. Monterey: IEEE, 2014, 27–36

35 KendallM G, Stuart A. The Advanced Theory of Statistics Volume 1: Distribution Theory. London: Charles Griffin & Company, 1958

36 Press W H, Teukolsky S A, Vetterling W T, Numerical recipes in C. Cambridge: Cambridge University Press, 1996

37 Ghosh D D, Olewnik A. Computationally efficient imprecise uncertainty propagation. Journal of Mechanical Design, 2013, 135(5): 051002

38 Ahlfeld R, Belkouchi B, Montomoli F. SAMBA: Sparse approximation of moment-based arbitrary polynomial chaos. Journal of Computational Physics, 2016, 320: 1–16

39 Patterson T. Modified optimal quadrature extensions. Numerische Mathematik, 1993, 64(1): 511–520

40 Kronrod A S.Nodes and Weights of Quadrature Formulas: Sixteen-place Tables. New York: Consultants Bureau, 1965

41 Patterson T. The optimum addition of points to quadrature formulae. Mathematics of Computation, 1968, 22(104): 847–856

42 Genz A, Keister B D. Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight. Journal of Computational and Applied Mathematics, 1996, 71(2): 299–309

43 Scarth C, Cooper J E, Weaver P M, Uncertainty quantification of aeroelastic stability of composite plate wings using lamination parameters. Composite Structures, 2014, 116: 84–93

44 Feinberg J, Langtangen H P. Chaospy: An open source tool for designing methods of uncertainty quantification. Journal of Computational Science, 2015, 11: 46–57

45 Huang B, Du X. Uncertainty analysis by dimension reduction integration and saddlepoint approximations. Journal of Mechanical Design, 2006, 128(1): 26–33

46 Gerstner T, Griebel M. Numerical integration using sparse grids. Numerical Algorithms, 1998, 18(3–4): 209–232

47 Heiss F, Winschel V. Likelihood approximation by numerical integration on sparse grids. Journal of Econometrics, 2008, 144(1): 62–80

48 Jaynes E T. Information theory and statistical mechanics. Physical Review, 1957, 106(4): 620–630

49 Phillips S J, Anderson R P, Schapire R E. Maximum entropy modeling of species geographic distributions. Ecological Modelling, 2006, 190(3–4): 231–259

50 Mohammad-Djafari A. A Matlab program to calculate the maximum entropy distributions. In: Smith C R, Erickson G J, Neudorfer P O, eds. Maximum Entropy and Bayesian Methods. Dordrecht: Springer, 1992, 221–233

51 Yeo S K, Chun J H, Kwon Y S. A 3-D X-band T/R module package with an anodized aluminum multilayer substrate for phased array radar applications. IEEE Transactions on Advanced Packaging, 2010, 33(4): 883–891

52 Pamies Porras M J, Bertuch T, Loecker C, An AESA antenna comprising an RF feeding network with strongly coupled antenna ports. IEEE Transactions on Antennas and Propagation, 2015, 63(1): 182–194