Polynomial Chaos Expansion for Probabilistic Uncertainty Propagation

Abstract

Uncertainty propagation (UP) methods are of great importance to design optimization under uncertainty. As a well-known and rigorous probabilistic UP approach, the polynomial chaos expansion (PCE) technique has been widely studied and applied. However, there is a lack of comprehensive overviews and studies of the latest advances of the PCE methods, and there is still a large gap between the academic research and engineering application for PCE due to its high computational cost. In this chapter, latest advances of the PCE theory and method are elaborated, in which the newly developed data-driven PCE method that does not depend on the complete information of input probabilistic distribution as the common PCE approaches is introduced and improved. Meanwhile, the least angle regression technique and the trust region scenario are, respectively, extended to reduce the computational cost of data-driven PCE to accommodate it to practical engineering design applications. In addition, comprehensive comparisons are made to explore the relative merits of the most commonly used PCE approaches in the literature to help designers to choose more suitable PCE techniques in probabilistic design optimization.

1. Introduction

Uncertainties are ubiquitous in engineering problems, which can roughly be categorized as aleatory and epistemic uncertainty [1, 2]. The former represents natural or physical randomness that cannot be controlled or reduced by designers or experimentalists, while the latter refers to reducible uncertainty resulting from a lack of data or knowledge. In systems design, all sources of uncertainties need to be propagated to assess the uncertainty of system quantities of interest, i.e., uncertainty propagation (UP). As is well known, UP is of great importance to design under uncertainty, which greatly determines the efficiency of the design. Since generally sufficient data are available for aleatory uncertainties, probabilistic methods are commonly employed for computing response distribution statistics based on the probability distribution specifications of input [3, 4]. Conversely, for epistemic uncertainties, data are generally sparse, making the use of probability distribution assertions questionable and typically leading to nonprobabilistic approaches, such as the fuzzy, evidence, and interval theories [5–7]. This chapter mainly focuses on propagating the aleatory uncertainties to assess the uncertainty of system quantities of interest using probabilistic methods, which is shown in Figure 1.

Figure 1.

Illustration of uncertainty propagation.

A wide variety of probabilistic UP approaches for the analysis of aleatory uncertainties have been developed [8], among which the polynomial chaos expansion (PCE) technique is a rigorous approach due to its strong mathematical basis and ability to produce functional representations of stochastic quantities. With PCE, the function with random inputs can be represented as a stochastic metamodel, based on which lower-order statistical moments as well as reliability of the function output can be derived efficiently to facilitate the implementation of design optimization under uncertainty scenarios such as robust design [9] and reliability-based design [10]. The original PCE method is an intrusive approach in the sense that it requires extensive modifications in existing deterministic codes of the analysis model, which is generally limited to research where the specialist has full control of all model equations as well as detailed knowledge of the software. Alternatively, nonintrusive approaches have been developed without modifying the original analysis model, gaining increasing attention, thus is the focus of this chapter. As a well-known PCE approach, the generalized PCE (gPCE) method based on the Askey scheme [11, 12] has been widely applied to UP for its higher accuracy and better convergence [13, 14] compared to the classic Wiener PCE [15]. Generally, the random input does not necessarily follow the five types of probabilistic distributions (i.e., normal, uniform, exponential, beta, and gamma) in the Askey scheme. In this case, the transformation should be made to transfer each random input variable to one of the five distributions. It would induce substantially lower convergence rate, which makes the nonoptimal application of Askey polynomial chaos computationally inefficient [8]. Therefore, the Gram-Schmidt PCE (GS-PCE) [16] and multielement PCE (ME-PCE) [17] methods have been developed to accommodate arbitrary distributions through constructing their own orthogonal polynomials rather than referring to the Askey scheme.

All the PCE methods discussed above are constructed based on the assumption that the exact knowledge of the involved joint multivariate probability density function (PDF) of all random input variables exists. Generally, by assumption of independence of the random variables, the joint PDF is factorized into univariate PDFs of each random variable in introducing PCE in the literature. However, the random input could exist as some raw data with a complicated cumulative histogram, such as bi-modal or multi-modal type, for which it is often difficult to obtain the analytical expression of its PDF accurately. Under these scenarios, all the above PCE approaches become ineffective since they all have to assume the PDFs to be complete. To address this issue, the data-driven PCE (DD-PCE) method has been proposed [18], in which its accuracy and convergence with diverse statistical distributions and raw data are tested and well demonstrated. With this PCE method, the one-dimensional orthogonal polynomial basis is constructed directly based on a set of data of the random input variables by matching certain order of their statistic moments, rather than the complete distributions as in the existing PCE methods, including gPCE, GS-PCE, and ME-PCE.

At present, great research achievements about PCE have been made in the literature, which have also been applied to practical engineering problems to save the computational cost in UP. However, there is still a large gap between the academic study and engineering application for the PCE theory due to the following reasons: (1) the complete information of input PDF often is not known in engineering, which cannot be solved by most PCE methods presented in the literature; (2) the computational cost of existing PCE approaches is still very high, which cannot be afforded in practical problems, especially when applied to design optimization; and (3) there is a lack of comprehensive exploration of the relative merits of all the PCE approaches to help designers to choose more suitable PCE techniques in design under uncertainty.

2. Data-driven polynomial chaos expansion method

Most PCE methods presented in the literature are constructed based on the assumption that the exact knowledge of the involved PDF of each random input variable exists. However, the PDF of a random parameter could exist as some raw data or numerically as a complicated cumulative histogram, such as bimodal or multimodal type, which is often difficult to obtain the analytical expression of its PDF accurately. To address this issue, the data-driven PCE method (DD-PCE for short in this chapter) has been proposed. DD-PCE follows the similar general procedure as that of the well-known gPCE method. For gPCE, the one-dimensional orthogonal polynomial basis simply comes from the Askey scheme in Table 1 and is a function of standard random variables. While for DD-PCE, the one-dimensional orthogonal polynomial basis is constructed directly based on the data of random input by matching certain order of statistic moments of the random inputs and is a function of the original random variables.

2.1. Procedure of data-driven PCE method

Step 1. Represent the output y as a PCE model of order p.

$y\approx {\displaystyle {\displaystyle \sum _{i=0}^P}{b}_i}{\text{Φ}}_i(X)={\displaystyle {\displaystyle \sum _{i=0}^P}{b}_i}{\displaystyle {\displaystyle \prod _{j=1}^d}{P}_j^{\left({\alpha }_j^i\right)}(X_j)}$ E1

where P+1 ( $1+P=\left(d+p\right)!/\left(d!p!\right)$ ) is the number of PCE coefficients b_i that is the same as gPCE; ${\text{Φ}}_i\left(X\right)$ is the d-dimensional orthogonal polynomial produced by the full tensor product of one-dimensional orthogonal polynomials $P_j^{\left({\alpha }_j^i\right)}(X_j)$ ; and ${\alpha }_j^i$ represents the order of $P_j^{\left({\alpha }_j^i\right)}(X_j)$ and clearly satisfies $0\le {\displaystyle {\displaystyle \sum _{j=1}^d}{\alpha }_j^i}\le p$ .

$P_j^{\left({\alpha }_j^i\right)}(X_j)$ corresponding to the jth dimensional random input variable x_j in Eq. (1) is defined as below, in which the index ${\alpha }_j^i$ is replaced by k_j for simplicity below:

$P_j{}^{\left(k_j\right)}(X_j)={\displaystyle {\displaystyle \sum _{s=0}^{k_j}}{p}_{s,j}^{\left(k_j\right)}\text{*}{\left(X_j\right)}^s,}j=1,2,\dots ,d$ E2

where $p_{s,j}^{\left(k_j\right)}$ is the unknown polynomial coefficient to be solved.

Step 2. Solve the unknown polynomial coefficient $p_{s,j}^{\left(k_j\right)}$ to construct the one-dimensional orthogonal polynomial basis.

Since the construction of $P_j^{\left({\alpha }_j^i\right)}(X_j)$ on each dimension is the same, the subscript j denoting the dimension number is omitted thereafter for simplicity. Based on the property of orthogonality, one clearly has

${\displaystyle {\int }_{x\in \text{Ω}}{P}^{\left(k\right)}\left(X\right)P^{\left(l\right)}\left(X\right)}d\text{Γ}\left(X\right)={\delta }_{kl},\forall k,l=0,1,\dots ,p$ E3

where δ_kl is the Kronecker delta, Ω is the original stochastic span, and Γ(X) represents the cumulative distribution function of the random variable X.

It is assumed that all the coefficients $p_s^{\left(k\right)}$ in Eq. (2) are not equal to 0, and then $P^{\left(0\right)}=p_0^{\left(0\right)}$ . For simplicity, the coefficient of the highest degree term in each P^(k) is set as $p_k^{\left(k\right)}=1,\forall k$ . According to Eq. (3), one has

${\displaystyle {\int }_{x\in \text{Ω}}{p}_0^{\left(0\right)}}\left[{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^s}\right]d\text{Γ}\left(X\right)=0$ E4

In the same way as above, one has

$\begin{array}{l}{\displaystyle {\int }_{x\in \text{Ω}}\left[{\displaystyle {\displaystyle \sum _{s=0}^1}{p}_s^{\left(1\right)}{X}^s}\right]}\left[{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^s}\right]d\text{Γ}\left(X\right)=0\\ ⋮⋮\\ {\displaystyle {\int }_{x\in \text{Ω}}\left[{\displaystyle {\displaystyle \sum _{s=0}^{k-1}}{p}_s^{\left(k-1\right)}{X}^s}\right]}\left[{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^s}\right]d\text{Γ}\left(X\right)=0\end{array}$ E5

There are totally k equations in Eqs. (4) and (5). Through substituting Eq. (4) into the first equation in Eq. (5), and then substituting Eq. (4) and the first equation in Eq. (5) to the second equation in Eq. (5), and so on, one set of new equations can be derived:

$\begin{array}{l}{\displaystyle {\int }_{x\in \text{Ω}}{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^s}}d\text{Γ}\left(X\right)=0\\ {\displaystyle {\int }_{x\in \text{Ω}}{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^{s+1}}}d\text{Γ}\left(X\right)=0\\ ⋮\\ {\displaystyle {\int }_{x\in \text{Ω}}{\displaystyle {\displaystyle \sum _{s=0}^k}{p}_s^{\left(k\right)}{X}^{s+k-1}}}d\text{Γ}\left(X\right)=0\end{array}$ E6

It is observed that ${\displaystyle {\int }_{\xi \in \text{Ω}}{X}^k}d\text{Γ}\left(X\right)$ is actually the kth order statistic moment of x, i.e., ${\displaystyle {\int }_{x\in \text{Ω}}{X}^k}d\text{Γ}\left(X\right)={\mu }_k$ . Therefore, Eq. (6) can be rewritten as

$\left[\begin{array}{cccc}{\mu }_0& {\mu }_1& \cdots & {\mu }_k\\ {\mu }_1& {\mu }_2& \cdots & {\mu }_{k+1}\\ ⋮& ⋮& ⋮& ⋮\\ {\mu }_{k-1}& {\mu }_k& \cdots & {\mu }_{2k-1}\\ 0& 0& \cdots & 1\end{array}\right]\left[\begin{array}{c}{p}_0^{\left(k\right)}\\ p_1^{\left(k\right)}\\ ⋮\\ p_{k-1}^{\left(k\right)}\\ p_k^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]$ E7

where ${\mu }_i(i=0,1,\dots ,2k-1)$ is the ith order statistic moment of x, which can be easily calculated from the given input data statistically or the PDFs of random inputs by integral. Of course, when the number of given data is not large enough, errors would be induced in the moment calculation.

Clearly, to obtain a k-order one-dimensional orthogonal polynomial basis, 0 to (2k − 1)-order statistic moments of x should be matched, which can be calculated based on the PDF or statistically based on the data set. Of course, when the number of data is not large enough, errors would be induced in the moment calculation. The polynomial coefficients for the one-dimensional orthogonal polynomial basis can be obtained by solving Eq. (7) with the Cramer’s Rule.

Step 3. Calculate the PCE coefficients b_i by the least square regression technique.

Step 4. Once the PCE coefficients are obtained, a stochastic metamodel (i.e., PCE model) that is much cheaper than the original model is provided. Evaluate on the PCE model by Monte Carlo simulation (MCS) to obtain the probabilistic characteristics of y. Since the PCE model is cheap, a large amount of sample points can be used. For the statistic moments, the analytical expressions can also be conveniently derived based on the PCE coefficients:

$\{\begin{array}{l}E\left[y\right]=E\left[{\displaystyle {\displaystyle \sum _{i=0}^P}{b}_i{\psi }_i\left(X\right)}\right]=b_0\\ {\sigma }^2\left[y\right]=E\left[y^2\right]-E^2\left[y\right]={\displaystyle {\displaystyle \sum _{i=0}^P}{b}_i^{2}E\left[{\psi }_i^2\left(X\right)\right]}-E^2\left[y\right]\\ Skew\left[y\right]=E\left[{\left(\frac{y-E\left[y\right]}{\sigma \left[y\right]}\right)}^3\right]=\frac{E\left[y^3\right]-3E\left[y\right]{\sigma }^2\left[y\right]-E^3\left[y\right]}{{\sigma }^3\left[y\right]}\\ Kur\left[y\right]=\frac{E\left[{\left(y-E\left[y\right]\right)}^4\right]}{{\sigma }^4\left[y\right]}=\frac{E\left[y^4\right]-4E\left[y\right]E\left[y^3\right]+6E^2\left[y\right]{\sigma }^2\left[y\right]+3E^4\left[y\right]}{{\sigma }^4\left[y\right]}\end{array}$ E8

2.2. Extension of Galerkin projection to DD-PCE

In the existing work about DD-PCE, only the regression method is employed to calculate the PCE coefficients. To the experience of the authors, the matrix during regression may become ill-conditioned during regression for higher-dimensional problems since the sample points required for regression that is often set as two times of the number of PCE coefficients P + 1 [19] is increased greatly causing a large-scale matrix during regression. To solve higher-dimensional problems, the Galerkin projection method in conjunction with the sparse grid technique has been widely used in gPCE due to its high accuracy, robustness, and convergence [20], which has also been observed and demonstrated during our earlier studies on PCE in recent years. In this section, the Galerkin projection method for PCE coefficients calculation is extended to the DD-PCE approach to address higher-dimensional UP problems. Figure 2 shows the general procedure of the improved DD-PCE method.

With the projection method, the Galerkin projection is conducted on each side of Eq. (1):

Figure 2.

Procedure of the improved DD-PCE.

$⟨y{\text{Φ}}_j\left(X\right)⟩=⟨{\displaystyle {\displaystyle \sum _{i=0}^P}{b}_i{\text{Φ}}_i\left(X\right)}{\text{Φ}}_j\left(X\right)⟩,\left(j=0,1,\cdots ,P\right)$ E9

where ⟨•⟩ represents the operation of inner product as below

$⟨g,f⟩={\displaystyle \int g}fdH\left(X\right)$ E10

where H(X) is the joint cumulative distribution function of random input variables X.

Based on the orthogonality property of orthogonal polynomials, the PCE coefficient can be calculated as

$b_i=E\left[y{\text{Φ}}_i\left(X\right)\right]/E\left[{\text{Φ}}_i\left(X\right){\text{Φ}}_i\left(X\right)\right],\left(i=0,1,\cdots ,P\right)$ E11

Similar to gPCE, the key point is the computation of the numerator in Eq. (11), which can be expressed as

$E\left[y{\text{Φ}}_i\left(X\right)\right]={\displaystyle {\int }_{\xi \in \text{Ω}}y{\text{Φ}}_i\left(X\right)d{\rm H}\left(X\right)}$ E12

The Gaussian quadrature technique, such as full factorial numerical integration (FFNI) and spare grid numerical integration, has been widely used to calculate the numerator in the existing gPCE approaches, with which the one-dimensional Gaussian quadrature nodes and weighs are directly derived by multiplying some scaling factors on the nodes and weights from the existing Gaussian quadrature formulae and then the tensor product is employed to obtain the multidimensional nodes. For some common type of probability distributions, for example, normal, uniform, and exponential distributions, their PDFs have the similar formulations as the weighting functions of the Gaussian-Hermite, Gaussian-Legendre, and Gaussian-Laguerre quadrature formula. Therefore, l_i and w_i can be conveniently obtained based on the tabulated nodes and weights of Gaussian quadrature formula [21], which are shown in Table 2, where $l_i^{G-H}$ and ${\omega }_i^{G-H}$ , $l_i^{G-La}$ and ${\omega }_i^{G-La}$ , $l_i^{G-Le}$ and ${\omega }_i^{G-Le}$ , respectively, represent the quadrature nodes and weights of Gaussian-Hermite, Gaussian-Laguerre, and Gaussian-Legendre quadrature formula; λ is the parameter of exponential distribution; and μ₁ and μ₀ denote the lower and upper bounds of uniform distribution.

Distribution types	PDFs	Polynomials	Weights	Intervals
Normal	$\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}$	Hermite Hn(x)	$e^{-x^2/2}$	[−∞, +∞]
Uniform	1/2	Legendre Pn(x)	1	[−1, 1]
Beta	$\frac{{\left(1-x\right)}^{\alpha }{\left(1+x\right)}^{\beta }}{2^{\alpha +\beta +1}B\left(\alpha +1,\beta +1\right)}$	Jacobi $P_n^{\left(\alpha ,\beta \right)}\left(x\right)$	${\left(1-x\right)}^{\alpha }{\left(1+x\right)}^{\beta }$	[−1, 1]
Exponential	e^−x	Laguerre Ln(x)	e^−x	[0, +∞]
Gamma	$\frac{x^{\alpha }{e}^{-x}}{\text{Γ}\left(\alpha +1\right)}$	General Laguerre $L_n^{\left(\alpha ,\beta \right)}$	x^αe^−x	[0, +∞]

Table 1.

Random variable types and the corresponding orthogonal polynomials.

Normal		Exponential		Uniform
li	ωi	li	ωi	li	ωi
$\sqrt{2}\sigma l_i^{G-H}+\mu$	$\frac{{\omega }_i^{G-H}}{\sqrt{\pi }}$	$\frac{l_i^{G-La}}{\lambda }$	${\omega }_i^{G-La}$	$\frac{{\mu }_1-{\mu }_0}{2}{l}_i^{G-Le}+\frac{{\mu }_1+{\mu }_0}{2}$	$\frac{{\omega }_i^{G-Le}}{2}$

Table 2.

l_i and ω_i calculated based on Gaussian quadrature.

However, the distributions of random inputs may not follow the Askey scheme, or are even nontrivial, or even exist in some raw data with a cumulative histogram of complicated shapes. Thus, such way to derive these nodes and weighs is not applicable in this case. In this work, a simple method is proposed based on the moment-matching equations below to obtain the one-dimensional quadrature nodes and weights.

$\begin{array}{l}{\omega }_0+{\omega }_1+\cdots +{\omega }_n={\displaystyle {\int }_{x\in \text{Ω}}1}d\text{Γ}\left(x\right)\\ {\omega }_0{l}_0+{\omega }_1{l}_1+\cdots +{\omega }_n{l}_n={\displaystyle {\int }_{x\in \text{Ω}}x}d\text{Γ}\left(x\right)\\ ⋮\\ {\omega }_0(l_0{)}^n+{\omega }_1(l_1{)}^n+\cdots +{\omega }_n(l_n{)}^n={\displaystyle {\int }_{x\in \text{Ω}}{x}^r}d\text{Γ}\left(x\right)\end{array}$ E13

where l_i and ω_i (i = 0, 1, …, n) are respectively the ith one-dimensional Gaussian quadrature nodes and weights, which theoretically can be obtained by solving Eq. (13).

However, Eq. (13) are multivariate nonlinear equations, which are difficult to solve when the number of equations is large (n + 1 > 7). It is noted that the one-dimensional polynomial basis P^(k) corresponding to each dimension constructed above is orthogonal. Therefore, its zeros are just the Gaussian quadrature nodes l_i, which can be easily obtained by solving P^(k) = 0. Through substituting l_i into Eq. (13), the n + 1 weights ω_i can be conveniently calculated. To calculate Eq. (13) of PCE order p, generally at least p + 1 one-dimensional nodes should be generated to ensure the accuracy, i.e., n ≥ p, which means that 0 to at least pth statistic moments of the random variable X should be matched. In this work, n is set as n = p.

In the same way, the nodes and weights in other dimensions are obtained conveniently. Then, the numerator can be calculated by the full factorial numerical integration (FFNI) method [8] for lower-dimensional problems (d < 4) as

$E[y{\text{Φ}}_i(X)]=\text{E}[Z(X)]\approx {\displaystyle {\displaystyle \sum _{i_1=1}^{m_1}}{\omega }_{i_1}}\cdots {\displaystyle {\displaystyle \sum _{i_j=1}^{m_j}}{\omega }_{i_j}}\cdots {\displaystyle {\displaystyle \sum _{i_d=1}^{m_d}}{\omega }_{i_d}}Z\left(l_{i_1},\cdots ,l_{i_j},\cdots ,l_{i_d}\right)\text{=}{\displaystyle {\displaystyle \sum _{j=1}^N}{W}_{j}Z(L_j)}$ E14

where $l_{i_j}$ and ${\omega }_{i_j}$ , respectively, represent the one-dimensional nodes and weights of the jth random input variable, which can be obtained using the way introduced above; L_i and W_i (i = 1, …, N) are the d-dimensional nodes and weights, respectively.

Generally, m is set as m ≥ p + 1 (p is the order of the PCE model). If the number of nodes N for calculating E[yΦ_i(X)] is too small, which is not matched with the PCE order, large error would be induced. Therefore, the conclusion that the higher the PCE order, the more accurate the UP results is based on the fact that E[yΦ_i(X)] has been calculated accurately enough. Clearly, the number of nodes N is increased exponentially with the increase of dimension d, causing curse of dimensionality. Therefore, FFNI is only suitable for lower-dimensional problems (d < 4).

For higher-dimensional problems (d ≥ 4), the sparse grid numerical integration method [22] can be used to calculate E[yΦ_i(X)] to reduce the computational cost:

$E\left[y{\text{Φ}}_i\left(X\right)\right]=E\left[Z\left(X\right)\right]\approx {\displaystyle {\displaystyle \sum _{q-d+1\le \left|i\right|\le q}^{}}{(-1)}^{q-\left|i\right|}\left(\begin{array}{l}d-1\\ q-\left|i\right|\end{array}\right)({\omega }_{i_1}\dots {\omega }_{i_j}\dots {\omega }_{i_{\text{d}}})}Z\left(l_{i_1},\cdots ,l_{i_j},\cdots ,l_{i_d}\right)$ E15

where $\left|i\right|=i_1+,\dots ,+i_d$ and i₁, …, i_d are the accuracy index corresponding to each dimension.

For the FFNI-based method, if m nodes are selected on each dimension (m₁ =…= m_d = m), 2m − 1 accuracy level can be obtained. For the sparse grid-based method, 2k + 1 accuracy level can be obtained with the accuracy level k = q − d. For example, if k = 2 and d = 8, for the sparse grid-based method, 17 nodes are required yielding 5th (2×2 + 1)-order accuracy level. For the FFNI-based method, to obtain the same accuracy level 5 (2×3 − 1), m should be m = 3 requiring 3⁸ nodes. Clearly, to obtain the same accuracy level, the number of nodes of the sparse grid-based method is much smaller than that of the FFNI-based method if d is relatively large.

In this chapter, we focus on extending the Galerkin projection to the DD-PCE method to address higher-dimensional UP problems and then exploring the relative merits of these PCE approaches. For the case with only small data sets, both DD-PCE and the existing distribution-based method (gPCE) may produce large errors for UP, and the estimation of PDF for the existing PCE methods is problem dependent and very subjective. It is difficult to make a comparison effectively between DD-PCE and the existing PCE methods. Therefore, during the comparison, it is assumed that there are enough data of the random input to ensure the accuracy of the moments.

2.3. Comparative study of various PCE methods

In this section, the enhanced DD-PCE method, the recognized gPCE method, and the GS-PCE method that can address arbitrary random distributions are applied to uncertainty propagation to calculate the first four statistic moments (mean μ, standard deviation σ, skewness β₁, kurtosis β₂) and probability of failure (P_f), of which the results are compared to help designers to choose the most suitable PCE method for UP. To comprehensively compare the three PCE approaches, four cases are respectively tested on four mathematical functions with varying nonlinearity and dimension shown in Table 3 and N, U, Exp, Wbl, Rayl, and Logn denote normal, uniform, exponential, Weibull, Rayleigh, and lognormal distribution, respectively. P_f is defined as P_f = probability (y ≤ 0).

Function 1: y = x1+x2+x3

Case 1: x1 ~U(1,2), x2 ~N(1,0.2), x3 ~ Exp(0.5)

Case 2: x1 ~Wbl(2,6), x2 ~Rayl(3), x3 ~ Logn(0,0.25)

Case 3: x1~BD, x2~ BD, x3~N(0,0.2)

Case 4: 500 and 10^7 sample points x1~BM, x2~ BM, x3~N(-0.8,0.2)

Function 2: y = sin(x1) − cos^2(x2) + x3sin(x1) + 0.9

Case 1: x1 ~N(0.5,0.2), x2 ~U(0,1.5), x3 ~ Exp(0.1)

Case 2: x1 ~Wbl(2,3), x2 ~Rayl(0.2), x3 ~ Logn(0,0.25)

Case 3: x1~ BD, x2~ BD, x3~U(0,1)

Case 4: 500 and 10^7 sample points x1~BM, x2~ BM, x3~U(0.4,2)

Function 3: y = e^−x1cos(x2) + x3e^−x4x5 − e^−x6

Case 1: x1 ~N(1,0.2), x2 ~U(−1,1), x3 ~ N(1,0.2), x4 ~ U(−1,1), x5 ~ N(0,0.2), x6 ~ U(0,2)

Case 2: x1 ~Wbl(1,5), x2 ~Rayl(0.5), x3 ~ Logn(0.5,0.25) ,x4 ~ Rayl (0.3), x5 ~ Wbl(1,5), x6 ~ Rayl(1)

Case 3: x1~ BD, x2~ BD, x3~ N(2,0.2), x4 ~ U(-1,0), x5 ~ N(1,0.2), x6 ~ U(−1,4)

Case 4: 500&10^7sample points x1~BM, x2~Rayl (0.3), x3~BM ,x4 ~Rayl (0.3), x5 ~BM, x6 ~Rayl (1)

Function 4: y = x1^2x2^2 − x3^2x4^2 + x5^2x6^2 − x7^2x8^2 + x9^2x10^2

Case 1: x1 ~ N(1,0.2), x2 ~ U(0,2), x3 ~ N(0,0.2), x4 ~ U(0,2), x5 ~ N(1,0.2), x6 ~ U(0,2), x7 ~ N(0,0.2), x8 ~ U(0,2), x9 ~ N(1,0.2), x10 ~ U(0,2)

Case 2: x1 ~Wbl(1,5), x2 ~Rayl(1), x3 ~Wbl(1,5), x4 ~Rayl(0.3), x5 ~Wbl(1,5), x6 ~Rayl(1), x7 ~Wbl(1,5), x8 ~Rayl(0.3), x9 ~Wbl(1,5), x10 ~Rayl(1)

Case 3: x1 ~ N(1,0.2), x2 ~ N(1,0.2), x3 ~ BD, x4 ~ BD, x5 ~ N(1,0.2), x6 ~ N(1,0.2), x7 ~ BD, x8 ~ BD, x9 ~ N(1,0.2), x10 ~ N(1,0.2)

Case 4: 500 and 10^7 sample points x1 ~ N(1.5,0.2), x2 ~ N(1,0.2), x3 ~ BM, x4 ~ BM, x5 ~ N(1,0.2), x6 ~ N(1,0.2), x7 ~ N(0,0.2), x8 ~ N(0,0.2), x9 ~ N(1,0.2), x10 ~ N(1,0.2)

Table 3.

Test functions and random input information of four cases.

The PCE order is set as p = 5 for all the functions for comparison, which means that 0–9th statistic moments of the random inputs should be matched to construct the one-dimensional orthogonal polynomials for the DD-PCE approach. For the first and second functions, FFNI-based Galerkin projection is used to calculate the PCE coefficients, while for the latter two, the sparse grid-based method with accuracy level k = 4 is used since the dimension is higher. The results of MCS with 10⁷ runs are used to benchmark the effectiveness of the three methods.

In Case 1, all the random input distributions are known and belong to the Askey scheme. The test results are shown in Tables 4–7, where the bold numbers with underline are the relatively best results and e represents the relative errors of the first four moments (μ, σ, β₁, β₂) with respect to MCS. P_f estimated by MCS is presented with 95% confidence interval. The results marked with * are from the sparse grid-based method. From these tables, it is found that with the same number of function calls (denoted as Ns), DD-PCE, gPCE, and GS-PCE produce almost the same results of the statistic moments, which are very similar to those of MCS (with the largest error as 2.6927%). The estimation of P_f for all the methods is within the 95% confidence interval with respect to MCS, indicating the high accuracy of UP. Although the orthogonal polynomial basis for DD-PCE is constructed by matching only 0–9th statistic moments of the random input variable instead of the complete PDFs for gPCE and GS-PCE, the results are accurate enough in this case. Moreover, the application of sparse grid technique to DD-PCE can greatly reduce the function calls for higher-dimensional problems (see Tables 5 and 6), while exhibiting good accuracy. Especially for the fourth function, with FFNI, the computational cost is very large (N_s = 976,562).

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0.0050	0	0.0100
eσ (%)	–	0.0164	(0.0164)	0.0164
$e_{{\beta }_1}$ (%)	–	0.1367	0.1367	0.1094
$e_{{\beta }_2}$ (%)	–	0.4877	0.3032	0.2199
Pf (1e^−3)	[8.5185,8.6328]	8.5472	8.5901	8.5688
Ns	10^7	125	125	125

Table 4.

Results of function 1 (Case 1).

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0.0115	0	0.0231
eσ (%)	–	0.0516	0	0.0258
$e_{{\beta }_1}$ (%)	–	0	0.4202	5.4852
$e_{{\beta }_2}$ (%)	–	0.1725	0.0814	0.0958
Pf ( 1e^−3)	[3.1403,3.2101]	3.1713	3.2017	3.1936
Ns	10^7	125	125	125

Table 5.

Results of function 2 (Case 1).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0	0.0112	0.0225
eσ (%)	–	0.0288	0.0288	(0.0288
$e_{{\beta }_1}$ (%)	–	2.2284	2.6927	1.7642
$e_{{\beta }_2}$ (%)	–	0.6040	0.6074	0.5028
Pf (1e^−3)	[4.8454,4.9318]	4.8993	4.8669	4.9074
Ns	10^7	1820	1820	1820

Table 6.

Results of function 3 (Case 1).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0.0123	0.0296	0.0074
eσ (%)	–	0.0402	0.0723	0.0522
$e_{{\beta }_1}$ (%)	–	0.1018	0.0890	0.1399
$e_{{\beta }_2}$ (%)	–	0.1050	0	0.1326
Pf (1e^−3)	[4.2476,4.3286]	4.2627	4.2881	4.2562
Ns	10^7	10,626	10,626	10,626

Table 7.

Results of function 4 (Case 1).

In Case 2, all the random input distributions are known but do not belong to the Askey scheme. In this case, the Rosenblatt transformation is employed for the gPCE method first. However, DD-PCE and GS-PCE can be directly used. The results are shown in Tables 8–11. It is observed that overall DD-PCE and GS-PCE perform better than gPCE, yielding results that are close to those of MCS. The reason is that the transformation in gPCE would induce error. Specifically, in Tables 9 and 10, the gPCE method causes relatively large errors due to the transformation. In addition, note the numbers with shadow, they are clearly larger than those of DD-PCE and GS-PCE, and P_f is outside the range of the 95% confidence interval of MCS. The interpretation is that since the first function is linear, the impact of transformation employed in gPCE on the accuracy of UP is small; while, for the second and third functions, they are more complicated and nonlinear (including trigonometric and exponential terms), the error induced by the transformation employed in gPCE is amplified more. The fourth function is a nonlinear polynomial one, which is easier to be handled than functions 2 and 3 in doing UP. Therefore, the results are generally accurate except P_f that is still outside the range of the 95% confidence interval of MCS. Moreover, the application of sparse grid greatly reduces N_s, exhibiting good potential applications for higher-dimensional problems.

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0.0196	0.0087	0.0175
eσ (%)	–	0.0298	0.0099	0.0199
$e_{{\beta }_1}$ (%)	–	0.2573	0.2059	0.2059
$e_{{\beta }_2}$ (%)	–	0.2170	0.2263	0.0899
Pf (1e^−4)	[1.9818,2.1602]	2.0360	2.1490	2.0480
Ns	10^7	125	125	125

Table 8.

Results of function 1 (Case 2).

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0.0243	0.0182	0.0061
eσ (%)	–	0.0467	0.2101	0
$e_{{\beta }_1}$ (%)	–	1.8877	8.0227	2.5956
$e_{{\beta }_2}$ (%)	–	0.0307	1.1659	0.0279
Pf (1e^−4)	[9.0052,9.3808]	9.0130	7.9720	9.0250
Ns	10^7	125	125	125

Table 9.

Results of function 2 (Case2).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0	0.0084	0
eσ (%)	–	0.0443	0.0887	0.0443
$e_{{\beta }_1}$ (%)	–	0.3471	0.6480	0.4397
$e_{{\beta }_2}$ (%)	–	0.0419	0.1927	0.1368
Pf (1e^−3)	[1.0859,1.1271]	1.0963	1.2291	1.1188
Ns	10^7	1820	1820	1820

Table 10.

Results of function 3 (Case2).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0.0240	0.0180	0.0320
eσ (%)	–	0.0111	0.0722	0.0250
$e_{{\beta }_1}$ (%)	–	0.2170	0.1979	0.2362
$e_{{\beta }_2}$ (%)	–	0.4229	1.9117	0.4582
Pf (1e^−3)	[4.4019,4.4843]	4.4635	4.6942	4.4200
Ns	10^7	10,626	10,626	10,626

Table 11.

Results of function 4 (Case 2).

In Case 3, the PDFs of some variables is bounded (BD) as below,

$f\left(x\right)=\{\begin{array}{c}2x,0<x<1\\ 0,\mathrm{otherwise}\end{array}$ E16

and the rest of the variables follow typical distributions. In this case, the Rosenblatt transformation is also employed for the gPCE method first.

From the results in Tables 12–15, it is found that generally large errors are induced by gPCE, especially the numbers with shadow in the tables. Since the first two variables follow the distribution bounded in an interval, the error induced by the transformation is large and all values of P_f are outside the confidence intervals for gPCE. While, the results of DD-PCE and GS-PCE are generally accurate and comparable, which are still very close to those of MCS. It should be noted that although the error of gPCE is the largest, all P_f by the three methods are outside the confidence intervals for function 3 (italic numbers) since this function is the most nonlinear and complicated. Hence, we increase the PCE order p and accuracy level k of the sparse grid to p = 6 and k = 5, and the results of P_f for DD-PCE, gPCE, and GS-PCE are 3.1263, 3.1446, and 3.1350, exhibiting evident improvement. Clearly with the same Ns, DD-PCE and GS-PCE are much more accurate than gPCE when nontrivial distribution is involved. These results further demonstrates the effectiveness and advantage of the enhanced DD-PCE for UP.

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0	0.0150	0
eσ (%)	–	0.0195	24.1063	0
$e_{{\beta }_1}$ (%)	–	0.1359	36.9565	0.1132
$e_{{\beta }_2}$ (%)	–	0.0545	12.3239	0.0545
Pf (1e^−3)	[4.9841,5.0717]	5.0038	5.2620	5.0333
Ns	10^7	125	125	125

Table 12.

Results of function 1 (Case 3).

Methods	MCS	DD-PCE	gPCE	GS-PCE
eμ (%)	–	0.0083	0.0914	0.0083
eσ (%)	–	0.0213	19.7662	0.0213
$e_{{\beta }_1}$ (%)	–	0.4186	123.2093	0.3256
$e_{{\beta }_2}$ (%)	–	0.0555	12.7841	0.0476
Pf (1^e−3)	[1.4429,1.4903]	1.4449	1.7890	1.4452
Ns	10^7	125	125	125

Table 13.

Results of function 2 (Case 3).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0.0359	0.7473	0.0598
eσ (%)	–	0.3983	(4.2798	0.3693
$e_{{\beta }_1}$ (%)	–	0.1221	22.5570	0.2036
$e_{{\beta }_2}$ (%)	–	0.6186	77.1134	0.6321
Pf (1e^−3)	[3.1972,3.2676]	2.6222	8.9269	2.6071
Ns	10^7	1820	1820	1820

Table 14.

Results of function 3 (Case 3).

Methods	MCS	DD-PCE*	gPCE*	GS-PCE*
eμ (%)	–	0.0039	6.4980	0.0194
eσ (%)	–	0.0409	8.2618	0.0164
$e_{{\beta }_1}$ (%)	–	0.1187	50.3681	0.0475
$e_{{\beta }_2}$ (%)	–	0.1720	11.8984	0.1949
Pf (1e^−3)	[8.6089,8.7237]	8.6559	0.8227	8.6728
Ns	10^7	10,626	10,626	10,626

Table 15.

Results of function 4 (Case 3).

In Case 4, the distributions of the random input variables are unknown and only some data exist. Although, based on the data, the analytical PDF can be obtained through some experience systems, such as Johnson or Pearson system [8], if the distribution of the data is very complicated, such as with a complicated cumulative histogram of bi- or multimodes, it is often very difficult to obtain the analytical PDF accurately. As is well-known that the Pearson system based on the first four statistic moments of the random variable would produce large errors for bimode (BM) or multimode PDFs. Evidently, the existing PCE approaches, including gPCE and GS-PCE, may produce large errors since they all depend on the exact PDFs of the random inputs in this case. However, DD-PCE can still work since it is a data-driven approach. To explore the effectiveness and advantage of DD-PCE over the other two approaches, it is assumed that the input data for some random input variables have a complicated bimode (BM) histogram shown in Figure 3 and the data for the rest from the typical distributions. Therefore, for the convenience and effectiveness of test, all the input data are generated based on the PDFs, of which the PDF of BM distribution is shown in Eq. (17). It should be pointed out that the PDFs actually are unknown and only some data exist in practice.

Figure 3.

PDF plot of the bimodal distribution.

$f_{\text{PDF}}=\frac{0.647}{0.1\sqrt{2\pi }}\exp \left(-\frac{x^2}{2\times {0.1}^2}\right)+\frac{0.353}{0.2\sqrt{2\pi }}\exp \left(-\frac{{\left(x-1\right)}^2}{2\times {0.2}^2}\right),x\in [-\infty ,+\infty ]$ E17

We tested small (500) and large (10⁷) numbers of input data to investigate the impact of number of data on the accuracy of UP. The results are shown in Tables 16–19, from which it is noticed that the results of DD-PCE are generally very close to those of MCS when the number of sample points of the random input variables is large (10⁷). When only 500 sample points are used, the errors are much larger. It means that the accuracy of DD-PCE is improved with the increase of the number of sample points. The reason is very simple that with the increase of the number of sample points, the statistic moments of random input variables calculated are more accurate, which would undoubtedly increase the accuracy of UP. The observation exhibits great agreements to what has been reported in work of Oladyshkin and Nowak. Similar to Case 3, the estimated P_f is outside the confidence intervals for function 3 since this function is the most nonlinear and the random distribution is more irregular, which can be improved by increasing N_s. This means that the generally the more nonlinear the function and the more irregular the random input distribution, the more difficult it is to achieve accurate UP results. These results further demonstrate the effectiveness and advantage of the enhanced DD-PCE method for UP.

Methods	MCS	DD-PCE (10^7)	DD-PCE (500)
eμ (%)	–	0.0066	1.4873
eσ (%)	–	0.0196	0.0688
$e_{{\beta }_1}$ (%)	–	0.0150	0.0451
$e_{{\beta }_2}$ (%)	–	0.0052	3.2327
Pf (1e^−3)	[1.4772,1.5252]	1.5069	0
Ns	10^7	125	125

Table 16.

Results of function 1 (Case 4).

Methods	MCS	DD-PCE(10^7)	DD-PCE(500)
eμ (%)	–	0.0132	0.4350
eσ (%)	–	0.0109	0.1957
$e_{{\beta }_1}$ (%)	–	0.1159	13.4783
$e_{{\beta }_2}$ (%)	–	0.0131	0.8956
Pf (1e^−3)	[6.4478,6.5474]	6.4703	8.000
Ns	10^7	125	125

Table 17.

Results of function 2 (Case 4).

Methods	MCS	DD-PCE(10^7)	DD-PCE(500)
eμ (%)	–	0.0327	0.6047
eσ (%)	–	2.7503	5.3717
$e_{{\beta }_1}$ (%)	–	3.8373	9.5932
$e_{{\beta }_2}$ (%)	–	0.5563	1.3573
Pf (1e^−3)	[7.7830,7.8924]	6.6667	6.0000
Ns	10^7	1820	1820

Table 18.

Results of function 3 (Case 4).

Methods	MCS	DD-PCE(10^7)	DD-PCE(500)
eμ (%)	–	0.0024	0.1925
eσ (%)	–	0.0241	3.5156
$e_{{\beta }_1}$ (%)	–	0.4149	214.4537
$e_{{\beta }_2}$ (%)	–	0.0170	11.9346
Pf (1e^−3)	[9.2650,9.3842]	9.2937	0
Ns	10^7	10626	10626

Table 19.

Results of function 4 (Case 4).

Function 1: f1 = X1X2

Case 1: 10^5 raw data

Case 2: X1~ $\mathcal{N}$ (1,0.2^2), X2~ $\mathcal{U}$ (0.4,1.6)

Case 3: X1 and X2 ~ $\mathcal{B}\mathcal{D}$

Function 2: $f_2=-X_1^2{X}_2^2-2X_3^4+3X_4^2-0.5X_5+4.5$

Case 1: 10^5 raw data

Case 2: X1~ $\mathcal{N}$ (1,0.2^2), X2~ $\mathcal{U}$ (0.4,1.6), X3~ $\mathcal{E}$ (0.1), X4~ $\mathcal{U}$ (−0.5,1),X5~ $\mathcal{U}$ (0.5,1).

Case 3: X1~ $\mathcal{B}\mathcal{D}$ , X2~ $\mathcal{U}$ (0.4,1.6).^2, X3~ $\mathcal{U}$ (0.5,1) .^2, X4~ $\mathcal{U}$ (−0.5,1), X5~ $\mathcal{U}$ (0.5,1).

Function 3: $f_3=-20\exp \left(-0.2\sqrt{\frac{1}{10}{\displaystyle {\displaystyle \sum _{i=1}^{10}}{x}_i^2}}\right)-\exp \left(\frac{1}{10}{\displaystyle {\displaystyle \sum _{i=1}^{10}}\cos \left(2\pi x_i\right)}\right)$

Case 1: 10^5 raw data

Case 2: X1~ $\mathcal{N}$ (1,0.2^2), X2~ $\mathcal{U}$ (0.4,1.6), X3~ $\mathcal{U}$ (−1.5,15), X4~ $\mathcal{U}$ (−1,2), X5~ $\mathcal{U}$ (−15, 1), X6~ $\mathcal{N}$ (2,0.2^2), X7~ $\mathcal{U}$ (−3,3), X8~ $\mathcal{U}$ (−15,1.5), X9~ $\mathcal{U}$ (−2,15), X10~ $\mathcal{U}$ (−2,15).

Case 3: X1 and X2 ~ $\mathcal{B}\mathcal{D}$ , X3~ $\mathcal{U}$ (−1.5,15), X4~ $\mathcal{U}$ (−1,2), X5~ $\mathcal{U}$ (−15,1), X6~ $\mathcal{N}$ (2,0.2^2), X7~ $\mathcal{U}$ (−3,3), X8~ $\mathcal{U}$ (−15,1.5), X9~ $\mathcal{U}$ (−2,15), X10~ $\mathcal{U}$ (−2,15).

Table 20.

Test functions.

To study the convergence property of the enhanced DD-PCE method, the errors (e) of moments and P_f with different PCE orders obtained by the proposed one as well as gPCE and GS-PCE are shown in Figures 4–7, taking Function 2, for example. Clearly, similar to the existing two methods, with the increase of the PCE order, the errors decrease significantly, exhibiting an approximate exponential convergence rate. Meanwhile, it is observed that the speed of convergence in Case 1 (Askey scheme) is the fastest. Generally, the more irregular the input distribution and the more nonlinear the function, the slower is the convergence process. In addition, it is also noticed that for Case 3, since x₁ and x₂ follow the nontrivial distribution, the convergence rate is very slow for gPCE (see left in Figure 6) due to the error induced by the transformation.

Figure 4.

Errors with respect to different PCE orders (Case 1).

Figure 5.

Errors with respect to different PCE orders (Case 2).

Figure 6.

Errors with respect to different PCE orders (Case 3).

Figure 7.

Errors with respect to different PCE orders (Case 4).

2.4. Summary

Overall, the three approaches produce comparably good results when the random inputs follow the Askey scheme. However, gPCE is the most mature and convenient to be implemented since there is no need to construct the orthogonal polynomials. When the PDFs of random inputs are unknown but do not follow the Askey scheme, large errors would be induced by the transformation for gPCE and the rest two PCE methods are comparable in accuracy and implementation complexity. It should also be pointed out that for DD-PCE, when constructing one-dimensional polynomials, the statistic moments (often 0–10 order) should be calculated first. If large gap exists between the high-order and low-order moments, the matrix singularity would happen in solving the linear equations (Eq. (7)). Therefore, in this case, GS-PCE is preferable especially when the function is highly nonlinear. When the PDF is unknown and cannot be obtained accurately, such as when random inputs exist as some raw data with a complicated cumulative histogram, only the DD-PCE method can still perform well since it is a data-driven method instead of the probabilistic-distribution-driven, while large errors would be produced if GS-PCE and gPCE are employed. However, more efforts should be made to solve the numerical problems in the DD-PCE method to make it more robust and applicable in constructing the one-dimensional orthogonal polynomials.

3. A sparse data-driven PCE method

The size of the truncated polynomial terms in the full PCE model is increased with the increase of the dimension of random inputs d and the order of PCE model p (see Eq. (1)), resulting in a significant growth of the computational cost. Therefore, attempts are made in this section on the full DD-PCE method introduced in Section 2 to reduce the computational cost. Accordingly, a sparse PCE approximation, which only contains a small number of polynomial terms compared to a classical full representation, is eventually provided by using the least angle regression (LAR) theory [23] and the sequential sampling method. The original LAR method is used for variables selection, aiming to find the most important variables with respect to a function response. In this work, LAR is extended to select some polynomial terms Ф_i(x) from the full PCE model that have the greatest impact on the model response $y\approx M(x)={\displaystyle {\displaystyle {\sum }_{i=0}^P}{b}_i{\text{Φ}}_i(x)}$ in a similar way.

Although the computational cost and accuracy are dependent on the PCE order, how to determine a suitable order that compromises between accuracy and efficiency is not within the scope of this chapter. In common situations, PCE of order p = 2 or 3 can produce results with good agreement to MCS for the output PDF estimation [24]. For more rigorous approaches of adaptively determining the order of the PCE model rather than specifying it in advance, readers can refer to references [25, 26].

3.1. Procedure of data-driven PCE method

A step-by-step description of the proposed method is given in detail as below with a side-by-side flowchart in Figure 8.

Figure 8.

The flowchart of the sparse DD-PCE method.

Step 1. Given the information of the random inputs (raw data or probabilistic distributions), specify the PCE order p, and then construct the full DD-PCE model without computing the PCE coefficients.

Step 2. Generate the initial input sample points X = [x₁,…,x_m,…,x_N]^T according to the distributions of the random inputs or select the sample points from the given raw data, where x_m = [x_m1,…,x_md]. Meanwhile, calculate the corresponding real function responses y = [y₁,…,y_m,…,y_N]^T.

X is standardized to have mean 0 and unit length, and that the response y has mean 0.

$\frac{1}{N}{\displaystyle {\displaystyle \sum _{m=1}^N}{x}_{\text{mn}}}=0,\sqrt{{\displaystyle {\displaystyle \sum _{m=1}^N}{x}_{\text{mn}}^2}}=1(n=1,\dots ,d),\frac{1}{N}{\displaystyle {\displaystyle \sum _{m=1}^N}{y}_m}=0$ E18

Then one has all the standardized data as

${\displaystyle \overline{X}}=\left[\begin{array}{l}{{\displaystyle \overline{x}}}_{11},{{\displaystyle \overline{x}}}_{12},\dots ,{{\displaystyle \overline{x}}}_{1d}\\ {{\displaystyle \overline{x}}}_{21},{{\displaystyle \overline{x}}}_{22},\dots ,{{\displaystyle \overline{x}}}_{2d}\\ \cdots \cdots \\ {{\displaystyle \overline{x}}}_{N1},{{\displaystyle \overline{x}}}_{N2},\dots ,{{\displaystyle \overline{x}}}_{Nd}\end{array}\right],{\displaystyle \overline{y}}={\left({{\displaystyle \overline{y}}}_1,\dots ,{{\displaystyle \overline{y}}}_N\right)}^T$ E19

Step 3. Set the iteration number as K = 0 and compute the values of all polynomial terms Ф_i(x) (i = 0, 1, …, P) of the full PCE model in Eq. (1) by, respectively, substituting each input sample point x_m into them. Then one obtains the information matrix as

$\text{Φ}=\left[\begin{array}{l}{\text{Φ}}_0\left(x_1\right){\text{Φ}}_1\left(x_1\right),\dots ,{\text{Φ}}_P\left(x_1\right)\\ ⋮\dots ⋮\\ {\text{Φ}}_0\left(x_N\right){\text{Φ}}_1\left(x_N\right),\dots ,{\text{Φ}}_P\left(x_N\right)\end{array}\right]$ E20

Step 4. The LAR algorithm is employed to automatically detect some number (often K + 1) of significant orthogonal polynomial terms from the first K + 1 terms Ф_i(x) (i = 0, 1, …, K) in Eq. (1), which will be retained to construct a sparse candidate PCE model that has a smaller scale than the full PCE model. For the introduction of the original LAR algorithm, readers can refer to reference [23] for more details.

Step 5. To save the computational cost, the leave-one-out cross-validation method [27] is employed to evaluate the accuracy of the candidate sparse PCE model constructed above, which is represented as the leave-one-out error analytically as below:

$Er{r}_{\text{LOO}}=\frac{1}{N}{{\displaystyle {\displaystyle \sum _{j=1}^N}\left(\frac{g\left(x_j\right)-{{\displaystyle \widehat{g}}}_I^{(-j)}\left(x_j\right)}{1-h_j}\right)}}^2$ E21

where g(x_j) is the response value from the original model at the sample point x_j; ${{\displaystyle \widehat{g}}}_I^{(-j)}$ represents the candidate sparse PCE model comprised of all the selected polynomial terms, of which the indices are stored in I; the PCE coefficients of ${{\displaystyle \widehat{g}}}_I^{(-j)}$ are computed through using the ordinary least-square regression method based on the leave-one-out approach, i.e., the sample points for regression are X^(-j) = [x₁,…,x_j−1, x_j+1,…,x_N]^T and y^(−j) = [y₁,…,y_j-1, y_j+1,…,y_N]^T.

Once the PCE coefficients are calculated, the predicted value by the candidate sparse PCE model at the sample point x_j is calculated as ${{\displaystyle \widehat{g}}}_I^{(-j)}\left(x_j\right)$ ; h_j is the jth diagonal element of the matrix Φ_A(Φ^T_AΦ_A)⁻¹Φ^T_A, where Φ_A is a N × k matrix comprised of all the selected column vectors Ф_i = [Ф_i(x₁), …, Ф_i(x_N)]^T (i ∈ I) and k is the number of selected polynomial terms.

To evaluate the accuracy more effectively, the relative error is employed based on Err_LOO as

${\epsilon }_{\text{LOO}}=Er{r}_{\text{LOO}}/{\displaystyle \widehat{V}}\left(y\right)$ E22

where ${\displaystyle \widehat{V}}\left(y\right)$ denotes the empirical variance of the response sample points y, which is calculated by

${\displaystyle \widehat{V}}\left(y\right)=\frac{1}{N-1}{{\displaystyle {\displaystyle \sum _{j=1}^N}\left(y_j-{\displaystyle \overline{y}}\right)}}^2,{\displaystyle \overline{y}}=\frac{1}{N}{\displaystyle {\displaystyle \sum _{j=1}^N}{y}_j}$ E23

Step 6. Check the stop criterion:

If the accuracy ε_LOO satisfies the target threshold ε, i.e., ε_Loo ≤ ε, the procedure will be stopped, the PCE model obtained by LAR in Step 4 will be considered as the final one, and all the sample points will be used for regression to calculate the PCE coefficients of the current sparse PCE model;

If ε_Loo > ε and K < P, set K = K + 1 and go to Step 4 to find another candidate sparse PCE model by LAR;

If ε_Loo > ε and K = P, generate some new sample points X_new with the sequential sampling technique and calculate the corresponding responses y_new, and add the new sample points into the old ones as X = [X; X_new] and y = [y; y_new], then go to Step 3 to find another candidate sparse PCE model.

In this work, if the PDF of random input is known, a large number of sample points are generated as the database according to the PDF beforehand; if the PDF of random input is unknown, the raw data are considered as the database. Each sample point in the database has its own index. The initial sample points are selected from the database through randomly and uniformly generating their indices. Then these sample points will be removed from the database and the rest will be indexed again. Similarly, by randomly and uniformly generating the indices, the sequential sample points will be selected from the reduced database. By using this sampling strategy, the sample points are distributed uniformly as far as possible, which is helpful to improve the accuracy of the PCE coefficient calculation.

Step 7. Based on the final sparse PCE model, the probabilistic properties of y can be obtained by running MCS or analytically.

3.2. Comparative study

In this section, the proposed sparse DD-PCE method (shortened as sDD-PCE hereafter) is applied to three mathematical examples to calculate the mean and variance of the output responses. The full DD-PCE (shortened as fDD-PCE hereafter) method that adopts a full PCE structure and one-stage sampling with the size of one times the number of PCE coefficients is also applied to UP, of which the results are compared to those of sDD-PCE to demonstrate its effectiveness and advantage.

The test examples of varying dimensions including their input information are shown in Table 20, in which the symbols $\mathcal{N},\mathcal{U}$ and $\mathcal{E}$ respectively, denote normal, uniform, and exponential distribution. To fully explore the applicability of sDD-PCE, three different cases of the random input information that almost cover all the situations in practice (Case 1: raw data; Case 2: common distribution; Case 3: nontrivial distribution) are considered. The nontrivial bimodal distribution (denoted as $\mathcal{B}\mathcal{D}$ ) used in Section 2.3 (Eq. (16)) is considered.

Another type of nontrivial distribution considered here is invented by conducting square operation on the sample points from some common distributions (see Case 3 in Function 2). The target accuracy ε of sDD-PCE is set as 10⁻⁵. Meanwhile, to ensure the effectiveness of comparison between sDD-PCE and fDD-PCE, the order of the PCE model p is set as the same (p = 3, 4, 5) for both methods. MCS with 10⁸ runs is conducted to benchmark the accuracy of both methods. In Case 1, the probabilistic distributions of all the random input variables are unknown and only a number of raw data (10⁵) exist, which cannot be solved by the traditional PCE methods, such as gPCE. Clearly, the more the raw data, the more reliable the results will be. Considering that the main objective of this paper is to investigate the effectiveness and capability of sDD-PCE in reducing the computational cost, it is assumed that a large number of raw data (10⁵) exist of the random inputs.

The results are listed in Tables 21–23, in which e_m and e_v, respectively, denote the errors (%) of mean and variance relative to the results of MCS, N denotes the number of total sample points (function evaluations) used for PCE coefficients estimation during regression, and Na represents that the result cannot be obtained.

	fDD-PCE			sDD-PCE
em	0.321	0.099	0.044	0.330	0.201	0.181
ev	0.232	0.813	0.173	0.203	0.099	0.068
p	3	4	5	3	4	5
N	10	15	21	20	30	20

Table 21.

Results of function 1 (Case 1).

	fDD-PCE			sDD-PCE
em	6.162	Na	Na	8.803	7.263	2.402
ev	10.182	Na	Na	16.670	5.026	8.882
p	3	4	5	3	4	5
N	56	126	252	20	20	30

Table 22.

Results of function 2 (Case 1).

	fDD-PCE			sDD-PCE
em	Na	Na	Na	0.045	0.739	0.239
ev	Na	Na	Na	18.134	12.882	2.479
p	3	4	5	3	4	5
N	286	1001	3003	30	30	30

Table 23.

Results of function 3 (Case 1).

From the results some noteworthy observations are made. First, generally with high PCE order (p = 5), the results of sDD-PCE are accurate. Second, for low-dimensional problem (d = 2, Function 1), the efficiency and accuracy of sDD-PCE and fDD-PCE are almost comparable. Specially, for lower orders p = 3 and 4, sDD-PCE is even less efficient. The interpretation is that in addition to the regression process, the sample points are also required during the construction of the sparse PCE model for sDD-PCE, while for fDD-PCE, the sample points are only used during regression. Moreover, for em with p = 3 (lower order) of Function 1, fDD-PCE is even more accurate with higher efficiency (see underlined numbers). The reason may be that for low-dimensional problems with low-order PCE models, the size of the total polynomial terms is already small and the sparse structure of sDD-PCE is of little help in reducing the number of sample points since additional sample points are required during the selection of important polynomial terms. Therefore, fDD-PCE may produce more accurate results than sDD-PCE since it maintains more information. This will be verified later. Third, with the increase of dimension (from d = 2, 5 to d = 10), N is increased significantly with the increase of p for fDD-PCE, causing matrix ill-conditioned problem. So some results (p = 4 and 5) even cannot be obtained by fDD-PCE. Specially, for Function 3, the dimension is high (d = 10), fDD-PCE cannot produce results for any order p. However, for sDD-PCE, no remarkable increase in N is noticed since it adopts a sparse PCE model that can adaptively remove the insignificant polynomial terms, while its accuracy is generally improved clearly exhibiting small error relative to MCS. When p = 5, only 13 polynomial terms are selected from 3003 total terms for Function 3; while for Function 1, 4 are selected from 21 total terms. Therefore, the larger the dimension, the more obvious the advantage of sDD-PCE over fDD-PCE in efficiency.

In Case 2, the PDFs of all the random inputs are known and assumed to follow common distributions. This is a general case that can be solved by the traditional probabilistic distribution-based PCE methods. The results are shown in Tables 24–26. Generally with high PCE order (p = 5), the results of sDD-PCE are accurate, demonstrating its effectiveness in dealing with random inputs with known PDFs. Meanwhile, for low-dimensional problem (Function 1), generally sDD-PCE is more accurate with the similar N as fDD-PCE. However, for lower order (p = 2) of Function 1, fDD-PCE is even more accurate than sDD-PCE, but with much smaller N. This observation is consistent with what has been noticed in Case 1 and the reason is that additional sample points are required to selecting important polynomial terms. With the increase of dimension, N is increased significantly with the increase of p for fDD-PCE. However, for sDD-PCE, remarkable improvement in the accuracy is noticed without a remarkable increase in N. These results show great agreements to what has been noticed in Case 1.

	fDD-PCE			sDD-PCE
em	0.083	0.044	0.060	0.710	0.010	0.100
ev	0.468	0.758	0.211	0.975	0.061	0.061
p	3	4	5	3	4	5
N	10	15	21	30	15	20

Table 24.

Results of function 1 (Case 2).

	fDD-PCE			sDD-PCE
em	24.401	Na	Na	1.244	0.490	0.216
ev	39.578	Na	Na	4.380	3.271	2.837
p	3	4	5	3	4	5
N	56	126	252	20	20	30

Table 25.

Results of function 2 (Case 2).

	fDD-PCE			sDD-PCE
em	Na	Na	Na	3.461	4.432	0.317
ev	Na	Na	Na	20.155	6.217	4.223
p	3	4	5	3	4	5
N	286	1001	3003	30	30	30

Table 26.

Results of function 3 (Case 2).

In Case 3, the PDFs of all the random inputs are known; however, some of them follow nontrivial distributions. In this case, the traditional gPCE method cannot work well since large errors would be induced in transforming such nontrivial distributions to certain ones in the Askey scheme. The results are shown in Tables 27–29, which exhibit great agreements to what has been observed in Case 1 and Case 2. The proposed sDD-PCE method can significantly reduce the number of sample points while with high accuracy. The higher the dimension, the more advantageous the adaptive sparse structure of sDD-PCE can be. In this case, only 11 polynomial terms are selected from 3003 total terms for d = 10 with sDD-PCE. Moreover, sDD-PCE can produce accurate and efficient results for nontrivial distributed random inputs.

	fDD-PCE			sDD-PCE
em	1.210	0.854	0.302	1.366	1.044	0.161
ev	2.321	0.748	0.815	0.805	0.161	0.000
p	3	4	5	3	4	5
N	10	15	21	10	10	10

Table 27.

Results of function 1 (Case 3).

	fDD-PCE			sDD-PCE
em	3.324	Na	Na	5.718	1.383	0.680
ev	7.855	Na	Na	7.634	7.322	2.290
p	3	4	5	3	4	5
N	56	126	252	20	30	30

Table 28.

Results of function 2 (Case 3).

	fDD-PCE			sDD-PCE
em	Na	Na	Na	4.114	2.212	0.112
ev	Na	Na	Na	48.894	15.817	3.101
p	3	4	5	3	4	5
N	286	1001	3003	30	30	30

Table 29.

Results of function 3 (Case 3).

To verify the guess that for low-dimensional problems with low-order PCE models, fDD-PCE may produce more accurate results than sDD-PCE since it maintains more information. Another test is conducted for Function 1 with lower order p = 2 with all the three cases, of which the results are shown in Table 30. Just as expected, fDD-PCE is clearly more accurate than sDD-PCE while generally with less sample points. For Function 1 with p = 2, the total number of polynomial terms is 6, which is very small. With sDD-PCE, only the last polynomial term is removed, while more points are required in removing the insignificant polynomials. So the sparse scheme does not have obvious impact under this circumstance. Therefore, it is concluded that the developed sDD-PCE method is particularly applicable to high-dimensional problems, especially those requiring a high order PCE model.

	Case 1		Case 2		Case 3
	fDD	sDD	fDD	sDD	fDD	sDD
em	0.2801	0.146	0.0366	0.244	0.414	0.807
ev	0.6344	0.367	0.3577	0.431	0.552	0.477
N	6	7	6	10	6	18

Table 30.

Results of function 1 (p = 2).

3.3. Summary

The developed sDD-PCE can reduce the number of polynomial terms in the PCE model, thus reducing the computational cost. Generally, the larger the random input dimension, the more obvious the advantage of the developed sDD-PCE over fDD-PCE in efficiency. The sDD-PCE method is much more efficient than fDD-PCE in solving high-dimensional problems, especially those requiring a high order PCE model.

4. Sparse DD-PCE-based robust optimization using trust region

In Section 3, to reduce the computational cost of DD-PCE, a sparse DD-PCE method has been developed by removing some insignificant polynomial terms from the full PCE model, thus decreasing the number of samples for regression in computing PCE coefficients. However, when the sparse DD-PCE is applied to robust optimization, it is conventionally a triple-loop process (see Figure 9): the inner one tries to identify the insignificant polynomial terms of the PCE model (the dash box); the middle is UP; the outer is the search for optima, which clearly is still very time-consuming for problems with expensive simulation models.

Figure 9.

The triple-loop formulation of sDD-PCE-based robust optimization.

As has been mentioned in Section 3, during each optimization iteration, although the sample points required for regression during UP of sDD-PCE are greatly reduced, certain additional number of sample points are required to identify the insignificant polynomial terms by the inner loop. If at some iteration design points, almost the same sparse polynomial terms are retained, the inner loop can clearly be avoided, thus saving the computational cost. To address this issue, the trust region technique widely used in nonlinear optimization is extended in this section. During optimizing, a trust region is dynamically defined. If the updated design point lies in the current trust region, it is considered that the insignificant terms of its PCE model remain unchanged compared to those of the last design point, i.e., the inner loop is eliminated at the updated design point. Meanwhile, to further save the computational cost, the sample points lying in the overlapping area of two adjacent sampling regions are reused for the PCE coefficient regression for the updated design point. The proposed robust optimization procedure employing sparse DD-PCE in conjunction with the trust region scenario is applied to several examples of robust optimization, of which the results are compared to those obtained by the robust optimization without the trust region method, to demonstrate its effectiveness and advantage.

4.1. The trust region scenario

The trust region method is a traditional approach that has been widely used in nonlinear numerical optimization [28]. The basic idea of the trust region method is that in the trust region of the current iteration design point, the second-order Taylor expansion is used to approximate the original objective function. If the accuracy of the current second-order Taylor expansion is satisfied, the size of the trust region is increased to speed up the convergence, and if not it is reduced to improve the accuracy of approximation. To reduce the computational cost of design optimization, the idea of the trust region technique has been extended and applied to reliability-based wing design optimization [29], multifidelity wing aero-structural optimization [30], and multifidelity surrogate-based wing optimization [31], which has been widely believed as an efficient strategy in design optimization. For example, when the trust region technique is applied to meta-model-based design optimization, during optimization, the sample points are sequentially generated in the trust region and the radius of the trust region is dynamically adjusted based on the accuracy of the meta-model in the local region.

4.2. Robust design using sparse data-driven PCE and trust region

The scenario of trust region is extended here to reduce the computational cost of sDD-PCE-based robust optimization. The basic idea is that the radius of a trust region is determined by the distance between two successive design points and the variation of the corresponding objective function values. If the updated design point ${\mu }_x^{k+1}$ lies in the current trust region, it is considered that the insignificant terms of its PCE model remain unchanged compared to those of the last design point ${\mu }_X^k$ , i.e., the inner loop is eliminated at the updated design point. Meanwhile, the sample points lying in the overlapping area of two adjacent sampling regions are reused for the PCE coefficient regression for the updated design point to further save the computational cost. Generally, for a practical engineering optimization problem, there is only one performance function that is computationally expensive. Therefore, only one PCE model is required to be constructed and the UP for the rest of the functions can be conveniently implemented by MCS. In this study, it is assumed that the PCE model is only constructed for the objective function and the general steps of the proposed method is as below.

Step 0: Set the iteration number as k = 1 and the initial staring design point ${\mu }_x^0$ , do robust optimization with sDD-PCE without trust region and obtain a new design variable ${\mu }_x^k$ , where the Latin Hypercube sample points are generated around ${\mu }_x^0$ to calculate the PCE coefficients.

Step 1: After the kth optimization iteration, define/update the trust region at the current obtained new design point ${\mu }_x^k$ as a rectangle with each length as

$r_1=max\left\{{\zeta }_1{\left|{\mu }_x^k\right|}_2,{\zeta }_2\right\},r_2=max\left\{{\zeta }_1\left|Y_{}^k\right|,{\zeta }_2\right\}$ E24

where ${\left|{\mu }_x^k\right|}_2=\sqrt{{\displaystyle {\displaystyle {\sum }_{i=1}^d}{\left({\mu }_{x_i}^k\right)}^2}}$ and |Y^k| is the absolute value of corresponding objective function at ${\mu }_x^k$ , i.e., $\left|Y_{}^k\right|=abs\left(Y({\mu }_x^k)\right)$ ; ζ₁ and ζ₂ are user-defined parameters, which can be constants or functions with respect to the iteration number k.

Step 2: During the (k + 1)th optimization iteration, the obtained new design point is ${\mu }_x^{k\text{+}1}$ . Before conducting UP, calculate the variation between two successive design points ${\mu }_x^k$ and ${\mu }_x^{k+1}$ as $\Delta x={\left|{\mu }_x^{k+1}-{\mu }_x^k\right|}_2=\sqrt{{\displaystyle {\displaystyle {\sum }_{i=1}^d}{\left({\mu }_{x_i}^{k+1}-{\mu }_{x_i}^k\right)}^2}}$ and the variation of the objective function $\Delta Y=\left|Y_{}^{k+1}\left({\mu }_x^{k\text{+1}}\right)-Y^k\left({\mu }_x^k\right)\right|$ .

Step 3: If ΔX ≤ r₂ and ΔY ≤ r₂ both are satisfied, ${\mu }_x^{k+1}$ is considered to be located in the trust region of ${\mu }_x^k$ defined in Eq. (45), and go to Step 4; if either ΔX ≤ r₁ or ΔY ≤ r₂ cannot be satisfied, ${\mu }_x^{k+1}$ is considered to be not located in the trust region of ${\mu }_x^k$ defined in Eq. (45), and go to Step 5.

Step 4: The retained polynomial terms Ф_i(x) at the updated new design point ${\mu }_x^{k+1}$ are kept as the same as those for the last design point ${\mu }_x^k$ , indicating that the inner loop of UP conducted on ${\mu }_x^{k+1}$ is removed. The Latin Hypercube sample points are generated around ${\mu }_x^{k+1}$ according to the distribution type and parameters of X with the same number of sample points as that used at the last design point ${\mu }_x^k$ to calculate the PCE coefficients. Meanwhile, the sample points located in the overlapping area of the two successive sampling regions are identified and reused for PCE coefficients calculation to improve the accuracy.

Step 5: The inner loop is conducted on the updated design point ${\mu }_x^{k+1}$ to detect the significant polynomial terms. Similarly, the sample points located in the overlapping area of the two successive sampling regions are also reused at the updated design point ${\mu }_x^{k+1}$ in detecting the significant polynomial terms and calculating the PCE coefficients to save the computational cost.

Step 6: Set k = k + 1, based on the results of UP, search for the next updated new design point ${\mu }_x^{k+1}$ and go to Step 1.

The above procedure will continue until the convergent criterion is satisfied. Figure 10 shows the case that the sample points in the previous optimization iteration are reused in the two successive iterations. As is seen that two points are located in the overlapping area of two successive sampling regions, thus are reused in the next iteration for regression to identify the significant polynomials/calculate the PCE coefficients. In this way, the computational cost can be further reduced.

Figure 10.

Illustration of the reuse of sample points.

4.3. Comparative studies

The first example is the Ackley Function:

$f\left(X\right)=-20\exp \left(-0.2\sqrt{\frac{1}{d}{\displaystyle {\displaystyle \sum _{j=1}^d}{X}_j^2}}\right)-\exp \left(\frac{1}{d}{\displaystyle {\displaystyle \sum _{j=1}^d}\cos \left(2\pi X_j\right)}\right)+22.71282,d=10$ E25

The robust design optimization of this example is:

$\begin{array}{l}\min F={\mu }_f+k{\sigma }_f\\ -10\le {\mu }_{X_j}\le 10,j=1,2,\dots ,10\end{array}$ E26

All the design variables are considered to follow uniform distribution with variation of ±0.2 around their mean values and k in Eq. (26) is set as k = 20. In this study, the fmincon function in Matlab is used for optimization, and ζ₁ and ζ₂ in Eq. (45) are set as ζ₁ = 0.5 and ζ₂ = 0.5. Meanwhile, the obtained optimal design variables of sDD-PCE-based robust design with and without trust region scenario as well as the deterministic design without considering any uncertainties are respectively substituted into Eq. (26) through MCS (with 1e⁶ runs) to calculate the mean μ_f and standard deviation σ_f of the objective function.

The results are shown in Table 31, from which it is found that compared to the robust optimization without the trust region scenario (denoted by without), the obtained performance results (μ_f, σ_f, and F) of the robust optimization with the trust region scenario (denoted by with) are comparable. However, the number of function calls (denoted as Funcall) is clearly reduced. The decrease in computational cost is attributed to the application of trust region scenario and the reusing of sample points. Meanwhile, the optimal designs of the two robust designs are both less sensitive to uncertainties (smaller σ_f) compared to the results of deterministic design (denoted by DD). These results demonstrate the effective and advantage of the proposed method.

	${\mu }_X^{\text{*}}$	μf	σf	F	Funcall
DD	[0,0,0,0,0,0,0,0,0,0]	1.8839	0.4390	10.6639	---
with	[0.6246,0.7066,0.6687,0.7796,0.5744,; 0.6784,0.7470,0.6333,0.6578,0.6904]	4.5014	0.1377	7.2554	12,735
without	[0.6564,0.6935,0.6984,0.7036,0.6691,; 0.0299,0.0141,0.6407,0.0205,0.0038]	3.7457	0.2003	7.7517	16,840

Table 31.

Results of the Ackley Function.

The second example is the robust design optimization of an automobile torque arm, shown in Figure 11.

Figure 11.

Automobile torque arm.

In this problem, the four geometrical parameters (a, d₁, d₂, and l) are considered as design variables, and the yielding strength S_y, Young’s modulus E, and the applied force Q are deterministic parameters.

$\begin{array}{l}\min f\left(a,d_1,d_2,l\right)=\frac{\pi a{d}_2^2}{4}+2\left(l-\frac{d_1}{2}-\frac{d_2}{2}\right)a^2\\ s.t.g_1\left(a,d_2,l\right)=\frac{Q\left(2l-d_2\right)d_2}{4I{S}_y}-1.0\le 0\\ g_2\left(a,d_1,d_2,l\right)=1.0-\frac{{\pi }^{2}E{a}^4}{3{\left(2l-d_1-d_2\right)}^2}\frac{d_2-d_1}{Ql}\le 0\\ 5\le a\le 15,45\le d_1\le 55,55\le d_2\le 65,110\le l\le 210\end{array}$ E27

where the objective function f represents the volume of the arm, the first constraint g₁ denotes the yielding failure at section A-A, the second constraint g₂ denotes the buckling failure at the two connecting rods, and $I=a^2{\left(d_2-a\right)}^2/2+a^4/6$ .

The distribution parameters of the four design variables and design parameters are shown in Table 32.

Random variables	Distribution	Lower bound	Upper bound
a	Uniform	μa −0.5 mm	μa +0.5 mm
d1	Uniform	μd1 −0.5 mm	μd1 +0.5 mm
d2	Uniform	μd2 −0.5 mm	μd2 +0.5 mm
l	Uniform	μl −0.5 mm	μl +0.5 mm
Parameters		Values
Q	Deterministic	5500 N
Sy	Deterministic	170 N/mm^2
E	Deterministic	2.1 × 10^10 N/mm^2

Table 32.

Distribution parameters for design variables and design parameters.

The corresponding robust design optimization model is formulated as

$\begin{array}{l}\min F={\omega }_1\frac{{\mu }_f}{{\mu }_f^{*}}+{\omega }_2\frac{{\sigma }_f}{{\sigma }_f^{*}},{\omega }_1=0.5,{\omega }_2=0.5\\ s.t.G_1\left(a,d_2,l\right)={\mu }_{g_1}+k{\sigma }_{g_1}\le 0\\ G_2\left(a,d_1,d_2,l\right)={\mu }_{g_2}+k{\sigma }_{g_2}\le 0\\ 5\le {\mu }_a\le 15,45\le {\mu }_{d_1}\le 55,55\le {\mu }_{d_2}\le 65,110\le {\mu }_l\le 210\end{array}$ E28

As has been mentioned above, the PCE model is only constructed for the objective function and the results are shown in Table 33. It is noticed that the robust optimization designs with and without the trust region scenario yields comparable results, while the function calls (objective function calls) required by design with trust region is evidently smaller. The deterministic design cannot even obtain a feasible optimal solution with both constraint violated (>0), since it does not consider uncertainties during design. These results further demonstrate the effectiveness and advantage of the proposed method.

	${\mu }_X^{\text{*}}$	μf	σf	F	G1	G2	Funcall
DD	[8.13, 55.00, 55.00, 110.00]	2.6616e^4	1.2171e^3	1	0.1848	5.1509e^4	82
with	[8.53, 54.10, 58.67, 111.03]	3.1027e^4	1.3355e^3	1.1315	−0.0123	−1.1833e^5	658
without	[8.57, 52.68, 57.50, 110.00]	3.0332e^4	1.3093e^3	1.1077	−4.0000e^−4	−1.2913e^2	1283

Table 33.

Results of automobile torque arm.

4.4. Summary

The employment of the trust region in sDD-PCE-based robust optimization can evidently reduce the computational cost. However, the determination of the trust region in this chapter is still very subjective and a more rigorous method should be explored. In this section as well as Section 3, the scenarios of sparse PCE and trust region are only employed to DD-PCE to save the computational cost. However, the methods proposed here are also applicable to other PCE approaches, such as gPCE and GS-PCE.

In this chapter, the latest advances in PCE theory and approach for probabilistic UP are comprehensively presented in detail. However, it does not limit the application of PCE to nonprobabilistic UP to address epistemic uncertainties. Sudret and Schöbi have proposed a two-level metamodeling approach using nonintrusive sparse PCE to surrogate the exact computational model to facilitate the uncertainty quantification analysis, in which the input variables are modeled by probability-boxes (p-boxes), accounting for both aleatory and epistemic uncertainty [32]. The Fuzzy uncertainty propagation in composites has been implemented using Gram-Schmidt polynomial chaos expansion, in which the parameter uncertainties are represented by fuzzy membership functions [5]. A general framework has been proposed for a dynamical uncertain system to deal with both aleatory and epistemic uncertainty using PCE, where the uncertain parameters are described through random variables and/or fuzzy variables [33]. The mix UP approach is proposed, in which the inner loop PDFs are calculated using the PCE, and outer loop bounds can be computed with optimization-based interval estimation [34]. PCE has also been applied for solving Bayesian inverse problem as “surrogate posterior.” However, it has been indicated that the accuracy cannot always be ensured by PCE since a sufficiently accurate PCE for this problem requires a high order, making PCE impractical compared to directly sampling the posterior [35].