*Article* **Modified Dimension Reduction-Based Polynomial Chaos Expansion for Nonstandard Uncertainty Propagation and Its Application in Reliability Analysis**

**Jeongeun Son and Yuncheng Du \***

Department of Chemical & Biomolecular Engineering, Clarkson University, Potsdam, NY 13699, USA; son@clarkson.edu

**\*** Correspondence: ydu@clarkson.edu; Tel.: +1-315-268-2284

**Abstract:** This paper presents an algorithm for efficient uncertainty quantification (UQ) in the presence of many uncertainties that follow a nonstandard distribution (e.g., lognormal). Using the polynomial chaos expansion (PCE), the algorithm builds surrogate models of uncertainty as functions of a standard distribution (e.g., Gaussian variables). The key to build these surrogate models is to calculate PCE coefficients of model outputs, which is computationally challenging, especially when dealing with models defined by complex functions (e.g., nonpolynomial terms) under many uncertainties. To address this issue, an algorithm that integrates the PCE with the generalized dimension reduction method (gDRM) is utilized to convert the high-dimensional integrals, required to calculate the PCE coefficients of model predictions, into several lower-dimensional ones that can be rapidly solved with quadrature rules. The accuracy of the algorithm is validated with four examples in structural reliability analysis and compared to other existing techniques, such as Monte Carlo simulations and the least angle regression-based PCE. Our results show our algorithm provides accurate UQ results and is computationally efficient when dealing with many uncertainties, thus laying the foundation to address UQ in complex control systems.

**Keywords:** dimension reduction; polynomial chaos expansion; uncertainty analysis; nonstandard distribution; statistical moments

#### **1. Introduction**

Uncertainty, originating from the inherent randomness of a complex system, is common in first-principle models widely used in various engineering problems. Uncertainty quantification (UQ), which quantitatively studies the impact of uncertainty on a system's performance, is important since uncertainty in model parameters (e.g., external loads, geometry, and material properties in a complex system) can degrade the accuracy of model predictions, thus affecting prediction-based control design and analysis. Many UQ approaches have been developed and applied to engineering problems to improve a system's performance [1–8]. For example, reliability assessment has been explored for the controlled structures of viscoelastic dampers [3], and the stochastic responses of bridge structures under uncertainty have been studied [4]. Further, reliability-based design optimization (RBDO), used to consider uncertainty in process design, has been developed for structural analysis, including the sequential approximate optimization to deal with multimodal random variables [6] and the adaptive surrogate-based algorithm to replace expensive-toevaluate models [7]. A hybrid reliability-based optimization method was also proposed for thermal structure design [9]. However, several challenges remain unsolved. For example, a major challenge of existing UQ tools is how to quantify the effect of a large number of uncertainties on model responses in an accurate and computationally efficient manner.

There are several popular UQ methods, including Monte Carlo (MC) simulations [10,11], Taylor expansion-based methods [12,13], and polynomial chaos expansion (PCE) [14,15].

**Citation:** Son, J.; Du, Y. Modified Dimension Reduction-Based Polynomial Chaos Expansion for Nonstandard Uncertainty Propagation and Its Application in Reliability Analysis. *Processes* **2021**, *9*, 1856. https://doi.org/10.3390/ pr9101856

Academic Editors: Francisco Ronay López-Estrada and Guillermo Valencia-Palomo

Received: 23 September 2021 Accepted: 15 October 2021 Published: 19 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Among these, MC has been widely used since it is easy to implement. It only requires random generation of samples from the distribution of uncertainty and repetitive execution of the first principle models with each sample. Despite the apparent simplicity, MC and its associated approaches are insufficient because they often require a large number of samples to ensure UQ accuracy, thus increasing the computational cost. To reduce the computational cost, different variants of MC were developed, which include an active learning approach to combine Kriging and MC simulation (i.e., AK-MCS) [8]. Further, the Taylor expansion-based method is another way for UQ, which is effective for less complicated models [16,17]. This includes the first-order reliability method (FORM) and second-order reliability method (SORM) [12,18] in structural reliability. However, when models involve complex nonpolynomial terms, the UQ accuracy of the Taylor expansionbased method cannot be guaranteed [12,13,18].

Recently, PCE has been used in different areas including the communities of structural reliability and control theory [4,19]. This approach was proposed by Wiener [14] to quantify uncertainty that follows a normal distribution and later extended to uncertainty with several standard distributions (e.g., uniform) by Xiu [20,21]. The PCE-based methods approximate uncertainty and its effect on model responses with a spectral expansion as a function of PCE coefficients. The PCE coefficients of uncertainty in models are given by a prior probability density function (PDF), but coefficients of model outputs are unknown. Thus, the key of PCE is to solve the PCE coefficients of model responses. Once the coefficients of model responses (e.g., performance function) are obtained, they can be used to approximate the PDF of performance function [19,20] and to evaluate failure probability in reliability analysis [18,22,23]. Since each model response can be represented explicitly with random variables of input uncertainties, it is useful to obtain the statistics of the response. Depending on how the PCE coefficients of model responses are calculated, it is categorized into intrusive and nonintrusive methods [20,24,25].

For the intrusive method, the coefficients of model outputs are obtained by projecting the governing equations into PCE, which requires calculating the inner product between a basis polynomial and the equations described with a truncated PCE in a stochastic Galerkin approach. The inner product generates a family of nested deterministic models to represent the original stochastic models. Calculating the coefficients with the inner product can be computationally prohibitive when models have complex functions (e.g., nonpolynomial functions) and many uncertainties. Nonintrusive PCE is another way, which has been used in the reliability analysis of bridge, truss, and frame structures [4,19,26,27]. Similar to MC, the convergence and accuracy of nonintrusive PCE depend on samples used for UQ. These samples are often referred to as collocation points. While practical, the convergence and accuracy of the nonintrusive PCE-based UQ can be affected by the total number of collocation points and by how the collocation points are generated [28,29]. The nonintrusive methods can be further categorized as a regression-based approach and projection-based approach, depending on how the PCE expressions of output responses are calculated. The former estimates the PCE coefficients by minimizing the mean square error of the approximation of model responses, while each PCE coefficient for the latter is defined with a multidimensional integral. It should be noted that as the number of uncertainties increase, both the size of dimension in the integral and the total number of PCE terms increase. Thus, the projection-based approach can be computationally demanding.

To address these issues, our previous works [29,30] have presented algorithms to combine the PCE with the generalized dimension reduction method (gDRM) or *M*-variate DRM [31]. This approach approximates high-dimensional integrals in the spectral projection (SP) [28] with a few low-dimensional ones. These low-dimensional integrals involve at most *M* integration variables [29] that are much smaller than the number of uncertainties. To further reduce the computational cost, the lower-dimensional integrals resulting from the gDRM were calculated with a sampling-based method (i.e., Gaussian quadrature rules) [30]. However, our previous works only deal with uncertainty that follows a standard distribution as done in the classic PCE theory.

Uncertainty may have nonstandard distributions. For example, lognormal distribution describes several natural phenomena, e.g., the bubble size distributions in the gas-liquid and gas-liquid-solid systems [32]. Because there is an explicit relationship between the lognormal and normal variables, we capitalize on the relationship following the work by Ghanem in [33] to efficiently deal with the lognormal uncertainty. Several tools were developed to approximate lognormal random variables [20,24,34–39], which include the isoprobabilistic transformation, such as the Rosenblatt or Nataf transformation [37,38], that transforms the original random variable into independent normal variables. Besides, the lognormal distribution was described in [39] by constructing orthogonal polynomials using the Stieltjes–Wigert polynomials. Compared to existing techniques, our algorithm approximates a lognormal distribution with a normal random variable based on the explicit relationship [33] and couples the PCE-based approximation with the gDRM to quantify the effect of many uncertainties on model responses in a computationally efficient way.

In summary, this work expands our previous work to couple the gDRM with the PCE to quantify the effect of nonstandard uncertainties on model outputs. For algorithm validation, lognormal and Weibull distributions are chosen as the testbed, which are approximated with standard random variables [20,33]. In addition, four examples in structural reliability analysis were chosen to show the efficiency of the algorithm and to discuss the UQ accuracy by comparing the results with recent works in the literature.

This paper is organized as follows. A brief description of the PCE-based UQ is given in Section 2. Section 3 shows procedures to approximate nonstandard uncertainty, followed by the gDRM-based PCE to deal with high-dimensional UQ problems. For algorithm verification, four numerical examples in structural reliability analysis are presented in Section 4, and a brief conclusion is given in Section 5.

#### **2. Background of Polynomial Chaos Expansion (PCE)**

Suppose that a stochastic process can be defined as *u* = *K*(*g*), where *g* = (*g*1, ..., *gn*) ∈ <sup>R</sup>*<sup>n</sup>* is a Gaussian vector involving *<sup>n</sup>* uncertain parameters (*<sup>n</sup>* <sup>≥</sup> 1) and *<sup>K</sup>* is the function to describe the relationship between a model response *u* and *g*. It is assumed each parameter in *g* is independent, i.e., any correlation among uncertain parameters is not considered. In this work, uncertainty in model parameters, i.e., input random variable of the models, is hereafter referred to as parametric uncertainty.

To obtain a PCE expression of *u*, the first step is to rewrite each parametric uncertainty *gi* as a function of a Gaussian variable *xi* as in [20]:

$$g\_i = g\_i(\mathbf{x}\_i) = \sum\_{k=0}^p \overline{\mathbb{g}}\_{i,k} H\_k(\mathbf{x}\_i) \tag{1}$$

where *gi*,*<sup>k</sup>* are the PCE coefficients to estimate the *i* th parametric uncertainty, which is defined by the *i* th standard normal distribution, i.e., *xi*~<sup>N</sup> (0, 1). These coefficients *gi*,*<sup>k</sup>* are often assumed to be a given a priori or can be estimated with parameter estimation techniques [40]. In (1), *Hk* is the one-dimensional Hermite polynomial basis function, and the polynomial order *p* defines the number of terms to approximate *gi*. Uncertainty in model parameters will introduce uncertainty in model responses. Thus, once the PCE coefficients of uncertainty in (1) are available, the PCE expansion of model response *u* can be subsequently defined with random variables as [20]:

$$\mu(\mathbf{x}) = K(\mathbf{g}) = \sum\_{j=0}^{m-1} \overline{\pi}\_{\mathbf{j}} H\_{j,n}(\mathbf{x}) \tag{2}$$

where *x* = (*x*1,..., *xi*,..., *xn*) and *Hj*,*n*(*x*) are the *j* th *n*-dimensional orthogonal polynomials defined by basis functions {*Hk*(*xi*)} for each uncertainty. Explicit expressions of *Hj*,*n*(*x*) can be found in [15]. In (2), *uj* are the deterministic PCE coefficients of model response *u*, which are unknown and have to be calculated. Further, *m* is the total number of terms to approximate *u*(*x*), which can be calculated with the number of parametric uncertainty *n* and the polynomial order *p* as in [20]:

$$m = \frac{(n+p)!}{p!n!} \tag{3}$$

Unlike *gi*,*<sup>k</sup>* , PCE coefficients of model response, *uj* , are unknown and have to be calculated by projecting (2) onto each polynomial basis function *Hj*,*<sup>n</sup>* and by using the spectral projection (SP), which can be defined as in [20]:

$$\overline{u}\_{\rangle} = \frac{\langle \mu(\mathbf{x}) H\_{\mathbf{j},\mathbf{n}}(\mathbf{x}) \rangle}{\langle H\_{\mathbf{j},\mathbf{n}}^2(\mathbf{x}) \rangle} \left( \forall j \in \{0, \dots, m - 1\} \right) \tag{4}$$

where the angle brackets · represent the inner product operator for the *n*-dimensional random space R*<sup>n</sup>* defined by *x*. For example, the inner product of a function *f*(*x*) is defined as:

$$
\langle f(\mathbf{x}) \rangle = \int\_{\mathbb{R}^n} f(\mathbf{x}) \mathcal{W}(\mathbf{x}) \, d\mathbf{x} \tag{5}
$$

where W(*x*) is the weighted function defined by the joint PDFs of random variables *x*. Note that the inner product of *f*(*x*) in (5) can be further specified by the expectation operator *E*[·], i.e., *E*[ *f*(*x*)] = *f*(*x*). Accordingly, the PCE coefficient *uj* in (4) can be further written as in [20]:

$$\overline{u}\_{\dot{j}} = \frac{1}{\gamma\_{\dot{j}}} \int\_{\mathbb{R}^n} u(\mathbf{x}) H\_{\dot{j},n}(\mathbf{x}) \mathcal{W}(\mathbf{x}) \, d\mathbf{x} \tag{6}$$

where *γ<sup>j</sup>* = *E H*<sup>2</sup> *j*,*n* is the normalization factor. It is worth mentioning that the *j* th normalization factor of the Hermite polynomial basis can be easily calculated as *γ<sup>j</sup>* = *j*! [20]. In addition, the multidimensional integral in (6) is solved over the entire random domain R*<sup>n</sup>* of *x*, and W(*x*) here is the Gaussian weighted function. By substituting the PCE coefficients in (6) into (2), uncertainty in model responses can be rapidly estimated. For example, the statistical central moments *μ<sup>u</sup>* and *σ*<sup>2</sup> *<sup>u</sup>*, i.e., mean and variance of model response *u*, can be analytically calculated as in [20]:

$$\mu\_u = E[u] = E\left[\sum\_{j=0}^{m-1} \overline{u}\_j H\_{j,n}\right] = \overline{u}\_0 E\left[H\_{0,n}\right] + E\left[\sum\_{j=1}^{m-1} \overline{u}\_j H\_{j,n}\right] = \overline{u}\_0 \tag{7}$$

$$\begin{split} \sigma\_{u}^{2} &= E\left[ \left\{ u - E[u] \right\}^{2} \right] = E\left[ \left\{ \sum\_{j=0}^{m-1} \overline{u}\_{j} H\_{j,u} - \overline{u}\_{j,0} \right\}^{2} \right] \\ &= E\left[ \left\{ \sum\_{j=1}^{m-1} \overline{u}\_{j} H\_{j,u} \right\}^{2} \right] = \sum\_{j=1}^{m-1} \left\{ \overline{u}\_{j} \right\}^{2} E\left[ H\_{j,u}^{2} \right] \end{split} \tag{8}$$

where *H*0,*<sup>n</sup>* ≡ 1 is the constant polynomial. As seen in (7) and (8), the first PCE coefficient *u*<sup>0</sup> can be used to estimate the mean, while the rest of the coefficients (*uj* =0) and the expectation of the squared basis functions (*E H*<sup>2</sup> *j* =0,*n* ) can be used to estimate the uncertainty around the mean of *u* [20]. Similarly, other higher-order statistical central moments, such as skewness *su* and kurtosis *κu*, can be rapidly calculated with the high-order PCE coefficients as in [41].

$$s\_{\boldsymbol{u}} = \mathrm{E}\left[\left\{\boldsymbol{u} - \mathrm{E}[\boldsymbol{u}]\right\}^3\right] \frac{1}{\sigma\_{\boldsymbol{u}}^3} = \mathrm{E}\left[\left\{\sum\_{j=0}^{m-1} \overline{\boldsymbol{u}}\_j H\_{j,\boldsymbol{n}} - \overline{\boldsymbol{u}}\_{j,0}\right\}^3\right] \frac{1}{\sigma\_{\boldsymbol{u}}^3} = \mathrm{E}\left[\left\{\sum\_{j=1}^{m-1} \overline{\boldsymbol{u}}\_j H\_{j,\boldsymbol{n}}\right\}^3\right] \frac{1}{\sigma\_{\boldsymbol{u}}^3} \tag{9}$$

$$\kappa\_{\boldsymbol{u}} = \mathbb{E}\left[\left\{\boldsymbol{u} - \mathbb{E}[\boldsymbol{u}]\right\}^4\right] \frac{1}{\sigma\_{\boldsymbol{u}}^4} = \mathbb{E}\left[\left\{\sum\_{j=0}^{m-1} \overline{u}\_j H\_{j,\boldsymbol{n}} - \overline{u}\_{j,0} \right\}^4\right] \frac{1}{\sigma\_{\boldsymbol{u}}^4} = \mathbb{E}\left[\left\{\sum\_{j=1}^{m-1} \overline{u}\_j H\_{j,\boldsymbol{n}} \right\}^4\right] \frac{1}{\sigma\_{\boldsymbol{u}}^4} \tag{10}$$

Details to calculate these statistical central moments can be found in [41–43]. These central moments are used to evaluate the reliability index and failure probability for structural reliability analysis in moment methods [18,22]. Thus, they are chosen in this work to discuss the performance of different UQ algorithms in terms of UQ accuracy and computational efficiency in Section 4.

#### **3. Efficient UQ for Nonstandard Uncertainties**

#### *3.1. Approximation of Lognormal Uncertainties*

Since lognormal uncertainty has an explicit relationship with Gaussian uncertainty, it will be transformed and estimated with the Gaussian variable as in [33]. For clarity, only a parametric uncertainty, which refers to an input random variable of the model, is considered, thus the subscript *i* in (1) will not be used.

Suppose a normally distributed uncertainty *g*(*x*) is approximated with a random variable *x*; and a lognormal uncertainty *l*(*x*) can be calculated by taking the exponential operator *exp*[·] of *g*(*x*), which results in a relationship as in [33]:

$$\lg(\mathbf{x}) = \ln l(\mathbf{x}), \ l(\mathbf{x}) = \exp[\mathbf{g}(\mathbf{x})] \tag{11}$$

In (11), *g*(*x*) will be referred to as the normalized Gaussian distribution, and the mean and the variance of *g*(*x*) are defined as *μ<sup>g</sup>* and *σ*<sup>2</sup> *<sup>g</sup>* , respectively. As done for the Gaussian uncertainty in (1), let one assume that *l*(*x*) can be approximated with a PCE expansion as in:

$$d(\mathfrak{x}) = \sum\_{k=0}^{p} \tilde{l}\_k H\_k(\mathfrak{x}) \tag{12}$$

where *lk* are the PCE coefficients, describing the original lognormal uncertainty, which will be derived based on normalized Gaussian variables. As defined in (1), *Hk* is the Hermite polynomial basis function. Since the first coefficient of PCE in (12) represents the mean value as described in (7), *l*<sup>0</sup> can be easily given as *μ<sup>l</sup>* = *e μg*+(*σ*<sup>2</sup> *<sup>g</sup>*/2) [33]. For other higher-order terms, we apply the SP to calculate PCE coefficients as a function of the normalized Gaussian random variables. Following the similar procedures to calculate the coefficients of model responses in (4), the resulting PCE expansion of the lognormal uncertainty in (12) can be represented as in [20,33]:

$$\begin{split} I(\mathbf{x}) &= \sum\_{k=0}^{p} \tilde{l}\_{k} H\_{k}(\mathbf{x}) = \tilde{l}\_{0} \sum\_{k=0}^{p} \frac{\sigma\_{\mathcal{S}}^{k}}{\mathcal{S}!} H\_{k}(\mathbf{x}) \\ &= \tilde{l}\_{0} \left( 1 + \sigma\_{\mathcal{S}} \mathbf{x} + \frac{1}{2!} \sigma\_{\mathcal{S}}^{2} (\mathbf{x}^{2} - 1) + \frac{1}{3!} \sigma\_{\mathcal{S}}^{3} (\mathbf{x}^{3} - 3\mathbf{x}) + \dots \right) \end{split} \tag{13}$$

Once the PCE coefficients of a lognormal uncertainty are available, the resulting uncertainty in model response as in (2) can be approximated following the similar procedures as in Section 2.

Compared to the PDF of the Gaussian random variable that has a symmetrical shape, the required polynomial order *p* in (13) for approximating the original lognormal uncertainty needs to be carefully selected to capture its asymmetric property, such as the long tail. As an example, Figure 1 shows the simulation results of a lognormal uncertainty, where three different polynomial orders (i.e., *p* = 2, 3, and 4) were used. Note that the PDFs shown in Figure 1 are estimated with the kernel density estimation function in MATLAB. The accuracy to approximate a lognormal distribution with different polynomial orders is compared to MC with 107 samples.

As seen in Figure 1, the accuracy to approximate a lognormal uncertainty can be affected by *p* and the magnitude of uncertainty. When uncertainty is small, the approximation of a lognormal uncertainty with different *p* values exhibits almost identical results as in Figure 1a. Whereas, when uncertainty is large, there is a noticeable difference between the PCE-based approximation and MC. However, when *p* was increased (>2), as seen in Figure 1b, the difference between the PCE-based approximation and MC is insignificant. Thus, we focus on two different values of *p* (2 vs. 3) in this work to show the performance of the UQ method, which will be discussed in Section 4.

**Figure 1.** Illustration of the effect of polynomial order on the approximation accuracy of a lognormal uncertainty. The PCE-based approximation of a lognormal uncertainty is built with the normalized Gaussian variable with different polynomial orders and compared to the MC-based method. The mean value of *l* was set to 1 in both graphs, but two different standard deviations were used: 0.1 was used in (**a**) and 0.3 was used in (**b**).

#### *3.2. Approximation of Other Nonstandard Uncertainties*

Compared to approximating a lognormal uncertainty with a PCE expression that has an explicit correlation to a Gaussian random variable, the approximation of other nonstandard uncertainties can be made with its probability distribution function [20]. To estimate a nonstandard uncertainty in PCE, suppose *Y* is an uncertainty with nonstandard distributions (e.g., a Weibull distribution), which has a PCE expression as:

$$Y(\mathbf{x}) = \sum\_{k=0}^{p} \overline{Y}\_k \Psi\_k(\mathbf{x}) \tag{14}$$

where *Yk* are PCE coefficients, and {Ψ*k*} are polynomial basis functions. Considering *Y* follows a nonstandard distribution, the calculation of *Yk* is difficult since the explicit correlation between *Y* and *x* is unknown [20]. In this work, a technique in [20] is used to estimate an uncertainty that follows a nonstandard distribution. The resulting PCE coefficients in (14) can be calculated as:

$$\overline{Y}\_{k} = \frac{1}{\gamma\_{k}} E\left[F\_{Y}^{-1}(F\_{\mathbf{x}}(\mathbf{x})) \Psi\_{k}(\mathbf{x})\right] = \frac{1}{\gamma\_{k}} \int\_{\mathbb{R}^{n}} F\_{Y}^{-1}(F\_{\mathbf{x}}(\mathbf{x})) \Psi\_{k}(\mathbf{x}) \mathcal{W}(\mathbf{x}) d\mathbf{x} \tag{15}$$

where *Fx* is the cumulative distribution function (CDF) defined as *Fx*(*x*) = *<sup>x</sup>* <sup>−</sup><sup>∞</sup> <sup>W</sup>(*t*)*dt*, and *F*−<sup>1</sup> *<sup>Y</sup>* is the inverse of a cumulative distribution function of *Y*, i.e., *FY*. For example, let one suppose that a Gaussian random variable *x* is used to approximate a Weibulldistributed uncertainty *Y*. Then, W(*x*) is the probability distribution function of *x*; Ψ*<sup>k</sup>* is the Hermite polynomial basis function; *Fx* is the CDF of the Gaussian random variable *x*; and *F*−<sup>1</sup> *<sup>Y</sup>* is the inverse of the CDF of Weibull-distributed uncertainty *Y*. Details about the derivation of (15) and its approximation are not given for brevity and can be found in [20].

#### *3.3. Modified gDRM-Based PCE Using Quadrature Rules*

When the PCE coefficients of uncertainty are available, such as (1) for the normal and (12) for lognormal uncertainty, SP is used to calculate the PCE coefficients of model response in (2), which requires calculating multivariate integrals as in (6). As in (3), the total number of terms (*m*) to accurately approximate uncertainty in model response in (2) is a function of the polynomial order (*p*) and the total number of uncertainties (*n*). When *p* is large (e.g., >2) and when the number of uncertainties increases, the number of terms in (2) can be greatly increased. This can increase the computational cost to calculate the PCE coefficients of model outputs. To reduce the computational cost, the gDRM will be used with quadrature rules to quickly calculate the PCE coefficients as in [30]. This approach is hereafter referred to as the modified gDRM-based PCE (or mgDRM-PCE).

Suppose that a continuous and differentiable function can be defined as *f*(*x*) = *u*(*x*)*Hj*,*n*(*x*), which is a portion of the integrand in (6). Then, the multivariate integral in (6) is rewritten as:

$$E[f(\mathbf{x})] = \int\_{\mathbb{R}^n} f(\mathbf{x}) \mathcal{W}(\mathbf{x}) \, d\mathbf{x} \tag{16}$$

where W(*x*) represents the weighted function described by the joint PDFs of *x*. Note that the function *f*(*x*) is defined to easily calculate a specific PCE coefficient *uj* in (6), and this expression can also be reused for different *j* = 0, 1, . . . , *m* − 1.

Following the definition of gDRM [31] and the representation in (16), each of the unknown PCE coefficients of the model response in (6) can be calculated as in [29]:

$$\overline{\pi\_j}(\mathbf{x}) = \frac{1}{\gamma\_j} E[f(\mathbf{x})] \cong \frac{1}{\gamma\_j} \sum\_{r=0}^{M} (-1)^r \binom{n-M+r-1}{r} E[f\_{M-r}] \tag{17}$$

where *fM*−*<sup>r</sup>* in the expectation operator *<sup>E</sup>*[·] is a function of (*<sup>M</sup>* − *<sup>r</sup>*) random variables selected from *x*. In this way, several lower-dimensional functions, involving up to *M* random variables of *x*, can be used to approximate the original function *f*, which is written as in [31]:

$$f\_{M-r} = \sum\_{d\_1 < d\_2 < \dots < d\_{M-r}} f\left(0, \dots, 0, \mathbf{x}\_{d\_1}, 0, \dots, 0, \mathbf{x}\_{d\_2}, 0, \dots, 0, \mathbf{x}\_{d\_{M-r}}, 0\right) \tag{18}$$

where *d*1, *d*2, ... , *dM*−*<sup>r</sup>* ∈ {1, 2, . . . , *n*}. Since the gDRM considers combinations of the random variables in *<sup>x</sup>*, the total number of terms in (18) can be calculated as *<sup>n</sup> M* − *r* . Notably, the last term in (18), when *M* = *r*, can be calculated as *f*<sup>0</sup> = *f*(0), which can be solved by setting the mean values of all random variables *x* to 0. For example, when a bi-variate dimension reduction method (BiDRM) is used, i.e., *M* = 2, the *n*-variate function *<sup>f</sup>*(*x*) will be approximated with *<sup>n</sup>* 2 two-variate functions, *<sup>n</sup>* 1 one-variate functions, and a constant term *<sup>f</sup>*0. Thus, the expectation of *fM*−*r*, i.e., *<sup>E</sup>*[ *fM*−*r*] in (17), results in a (*M* − *r*)-dimensional integral for each combination of (*M* − *r*)-random variables selected from *<sup>x</sup>*, which gives *<sup>n</sup> M* − *r* (*M* − *r*)-dimensional integrals in total.

As noted in [29,30], the accuracy to approximate a high-dimensional integral with several lower-dimensional ones increases as *M* increases. For example, compared to the BiDRM, a tri-variate DRM (TriDRM) approximates a high-dimensional integral in (16) with *<sup>n</sup>* 3 three-dimensional integrals, *<sup>n</sup>* 2 two-dimensional integrals, *<sup>n</sup>* 1 one-dimensional integrals, and a constant term *f*0. As such, extra effort is required to calculate *<sup>n</sup>* 3 additional three-dimensional integrals. Since the total number of lowerdimensional integrals increases, the computational time to approximate (16) may increase. To address this issue, we will estimate the resulting lower-dimensional integrals with quadrature rules [44] to reduce the computational cost [30,31]. It should be noted that this approach reduces the computational difficulty to evaluate high-dimensional integrals using sampling methods, which appears to be methodologically similar to the anchored analysisof-variance (ANOVA) expansion in [45]. However, the anchored ANOVA expansion is dependent on the choice of the anchor point, which impacts the accuracy of the expansion and the truncation dimension [45], while the mgDRM does not require any strategy for the decomposition procedure in (18). A detailed description on the ANOVA expansion can be found in [45,46].

To calculate the expectation of one-variate function *E*[ *f*1]—an integral only involves a single random variable, e.g., *xd*<sup>1</sup> of *f*<sup>1</sup> in (18), a one-dimensional quadrature rule can be defined as [30]:

$$\int\_{-\infty}^{\infty} f\left(0, \dots, 0, \mathbf{x}\_{d\_1}, 0, \dots, 0\right) \mathcal{W}(\mathbf{x}\_{d\_1}) d\mathbf{x}\_{d\_1} \cong \sum\_{q\_1=1}^{\theta\_1} f\left(0, \dots, 0, \mathbf{x}\_{d\_1}^{q\_1}, 0, \dots, 0\right) \cdot a\_{d\_1}^{q\_1} \tag{19}$$

where *x q*1 *d*1 , *α q*1 *d*1 *θ*<sup>1</sup> *q*1=1 is a set of quadrature points (*x q*1 *d*1 ) and their corresponding weights (*α q*1 *d*1 ) used to compute the one-dimensional integral resulting from the gDRM. In this work, Gauss–Hermite quadrature rules [44] are used to estimate low-dimensional integrals in (17) in the gDRM step.

To generate quadrature points to approximate multidimensional integrals in (17), a full tensor product grid can be constructed based on the one-dimensional quadrature rules. Note that these lower-dimensional integrals resulting from the gDRM only involve a few integration variables. Thus, the approximation with the full tensor product grid can be finished in real time.

To better illustrate the quadrature points-based approximation, suppose the PCE coefficients of model outputs in (17) are calculated with the TriDRM. This implies that these lower-dimensional integrals to approximate the multidimensional integrals in SP involve at most three random variables in *x* (or *M* = 3). That is, a multidimensional integral (*E*[ *<sup>f</sup>*(*x*)]) is estimated with *<sup>n</sup>* 1 one-dimensional, *<sup>n</sup>* 2 two-dimensional, and *<sup>n</sup>* 3 three-dimensional integrals in total, which can be computed with the quadrature rules as in [30]:

$$E[f\_1] \cong \sum\_{d\_1} \left\{ \sum\_{q\_1=1}^{\theta\_1} f\left(0, \dots, 0, \mathbf{x}\_{d\_1}^{q\_1}, 0, \dots, 0\right) \cdot \mathbf{a}\_{d\_1}^{q\_1} \right\} \tag{20}$$

$$E[f\_2] \cong \sum\_{d\_1 < d\_2} \left\{ \sum\_{q\_1=1}^{\theta\_1} \sum\_{q\_2=1}^{\theta\_2} f\left(0, \dots, 0, \mathbf{x}\_{d\_1}^{q\_1}, 0, \dots, 0, \mathbf{x}\_{d\_2}^{q\_2}, 0, \dots, 0\right) \cdot \left(\mathbf{a}\_{d\_1}^{q\_1} \otimes \mathbf{a}\_{d\_2}^{q\_2}\right) \right\} \tag{21}$$

$$E[f\_3] \cong \sum\_{d\_1 < d\_2 < d\_3} \left\{ \sum\_{q\_1 = 1}^{\theta\_1} \sum\_{q\_2 = 1}^{\theta\_2} \sum\_{q\_3 = 1}^{\theta\_3} f\left(0, \dots, 0, x\_{d\_1}^{q\_1}, 0, \dots, 0, x\_{d\_2}^{q\_2}, 0, \dots, 0, x\_{d\_3}^{q\_3}, 0\right) \right. \tag{22}$$
 
$$\left\{ a\_{d\_1}^{q\_1} \otimes a\_{d\_2}^{q\_2} \otimes a\_{d\_3}^{q\_3} \right\}$$

where *θ<sup>i</sup>* (*i* = 1, 2, 3) is the number of quadrature points, *x qi di* , for each integration variable; *α qi di* is the weight of each quadrature point; and ⊗ means the tensor product. In summary, the total number of quadrature points for low-dimensional integrals can be defined as *<sup>Q</sup>* <sup>=</sup> *<sup>M</sup>*−*<sup>r</sup>* ∏ *i*=1 *θ<sup>i</sup>* = *θM*−*<sup>r</sup> <sup>i</sup>* [44]. It is important to note that as the dimension of a random space increases, the number of integrals resulting from gDRM increases. In this case, the total number of evaluations of these integrals with quadrature points is also increased, which may increase the computational cost. This issue will be discussed with benchmark examples in Section 4.

#### *3.4. Summary of the UQ Algorithm*

Our algorithm quantifies the effect of nonstandard uncertainty on model responses by coupling PCE with the gDRM. A flowchart summarizing the procedures of the modified gDRM-based PCE to deal with nonstandard random variables is shown in Figure 2.

As in Figure 2, nonstandard uncertainties are firstly approximated with standard random variables. Once the PCE-based surrogate model of parametric uncertainty, i.e., an input random variable of the models is available, model responses can be approximated with coupled PCE coefficients as in (2). The PCE coefficients of model responses are calculated with (4). This produces high-dimensional integrals, which are not trivial to solve, especially when the number of uncertainties is large and when the model involves nonpolynomial terms. To address this, the gDRM, such as the BiDRM and TriDRM noted in Section 3.3, is used to convert the high-dimensional integral into few low-dimensional ones. To further decrease the computational cost, quadrature points generated with Gaussianquadrature rules are used to numerically solve these low-dimensional integrals. The approximation of lower-dimensional integrals is highlighted in green boxes in Figure 2. Once the PCE coefficients of model responses are available, the statistical moments of

model outputs (e.g., mean, variance, skewness, and kurtosis) can be quickly calculated with (7)~(10).

**Figure 2.** Schematic of the modified gDRM-based PCE for UQ under many uncertainties that follow nonstandard distributions.

#### **4. Benchmark Examples in Structural Reliability Analysis**

Four examples in structural reliability analysis are chosen to illustrate the performance of the UQ algorithm to deal with nonstandard uncertainties. The correlation among uncertainties is not considered, since our objective is to study the applicability of the algorithm to tackle many uncertainties. When uncertainties are correlated, the number of random variables to approximate uncertainties will be decreased, thus reducing the number of dimensions. To evaluate the UQ accuracy, statistical moments (mean, standard deviation, skewness, and kurtosis) are calculated and compared with other techniques, which include MC simulations and the least angle regression-based PCE.

#### *4.1. Example 1: Linear Performance Function*

As given in Figure 3, a linear performance function is used to study the plastic collapse mechanisms of a one-bay frame, which is mathematically described as in [47,48]:

$$Z = \mathfrak{u}(X) = X\_1 + 2X\_2 + 2X\_3 + X\_4 - 5X\_5 - 5X\_6 \tag{23}$$

Following [47,48], *X*1~*X*<sup>6</sup> in (23) are assumed to be independent and lognormally distributed; the mean and standard deviation of each variable are listed in Table 1. In this case study, all parameters are considered as parametric uncertainties, resulting in a six-dimensional random space. Since uncertainty follows a lognormal distribution, Hermite polynomials were used as the basis functions to build the PCE-based surrogate model of the response *Z* in (23).

The first step to build a surrogate model of *Z* in (23) is to formulate the PCE expressions of lognormally distributed uncertainties, *X*1~*X*6. In this case study, (13) was used to build a PCE model for each uncertainty with normalized Gaussian variables that are related to the lognormal uncertainty as in Section 3.1. The second coefficients of the normalized Gaussian variables (i.e., *σ<sup>g</sup>* for *X*1~*X*6) were calculated as 0.0998 for *X*1~*X*<sup>4</sup> and 0.2936 for *X*<sup>5</sup> and *X*6, respectively. These were also used to derive the explicit PCE expressions of lognormal uncertainties as in (13). Once the PCE expressions of uncertainties (*X*1~*X*6) were available, the mgDRM-PCE in Section 3.3 was used to compute the PCE coefficients of *Z* in (23).

**Figure 3.** Illustration of the one-bay frame for a plastic collapse mechanism [48]. Reprinted from [48], Copyright (2019), with permission from Elsevier.

**Table 1.** Details of the uncertain variables in Example 1 [47,48]. Reprinted from [48], Copyright (2019), with permission from Elsevier.


In this case study, different polynomial orders of *X*1~*X*<sup>6</sup> (*p* = 2 vs. 3) and different numbers of random variables in the gDRM (*M* = 2 vs. 3) were used. The UQ accuracy with different combinations of *p* and *M* was compared to the results of MC by calculating the relative error (*R*), which is defined as in:

$$
\epsilon\_R = \left| \frac{\varepsilon\_{\rm MC} - \varepsilon}{\varepsilon\_{\rm MC}} \right| \tag{24}
$$

where *eMC* is a reference of a specific statistical moment calculated with MC, and *e* represents the corresponding statistical moment calculated with other UQ methods. The simulation results are summarized in Table 2 and compared with existing algorithms: a recent work in [48] that integrates the BiDRM with a high-order unscented transformation [49] and the least angle regression-based PCE (LAR-based PCE) [50,51]. For MC, 107 samples were used for each uncertainty to ensure UQ accuracy.

In Table 2, mBiDRM-PCE means that the quadrature rules were used to estimate integrals from the BiDRM, which involve at most two integration variables, while mTriDRM-PCE means that the quadrature rules were applied to the TriDRM, which has at most three integration variables. The BiDRM converted a six-dimensional integral into 15 two- and 6 one-dimensional ones, while 20 three-, 15 two-, and 6 one-dimensional integrals were used in the TriDRM. In addition, five quadrature points of each dimension were generated following the Gauss–Hermite quadrature rules (i.e., *θ<sup>i</sup>* = 5), for which a full-tensor product grid was built to calculate the low-dimensional integrals from the BiDRM and the TriDRM.

In addition, for the LAR-based PCE, the LAR procedure was implemented to select a set of basis polynomials to build a sparse PCE and estimate the unknown PCE coefficients of the response *Z* in (23) with the least square regression. It is important to note the LARbased PCE in this work does not involve the adaptive solver as in [52] because our objective is to assess the accuracy of the surrogate model derived by the gDRM-based PCE. To obtain model evaluations of *Z*, 1000 pairs of random samples of parametric uncertainties were

used. Thus, the total number of model runs is 1000. Further, the LAR method was used to identify the sensitive basis polynomials of output *Z* until the maximum correlation between the residual and the basis polynomials is less than a preselected effective value, which was set to 1.0 × <sup>10</sup>−6. In this case study, two different polynomial orders (*<sup>p</sup>* = 2 and 3) were studied. Thus, the number of PCE terms of LAR-based PCE was 13 and 19, respectively. In contrast, the number of PCE terms for our algorithm was 28 and 84, respectively.


**Table 2.** Summary of statistical moments and their relative errors in Example 1.

As seen in Table 2, the difference in the mean value of *Z* (*μZ*) calculated with the mBiDRM-PCE and mTriDRM-PCE is insignificant, and their results are almost identical to MC. For the other higher-order statistical moments, it was found that the relative error is relatively larger when *p* was set to 2, as compared to the results when *p* was set to 3. This indicates that more polynomial terms to approximate a lognormal uncertainty can improve the UQ accuracy. Besides, compared to the LAR-based PCE, our algorithm shows identical results for the four statistical moments, which demonstrates the accuracy of the gDRM-based PCE. It was also found that, when *p* was set to 3, the mBiDRM-PCE outperforms the results in a recent work [48], for which only the BiDRM was used. This shows the advantage of combining the PCE with the gDRM, because PCE can not only provide an analytical expression for model responses as a function of random variables, but also explicitly quantify the impact of approximation on UQ accuracy.

For comparison, the simulation results of mBiDRM-PCE and mTriDRM-PCE with different polynomial orders and MC are shown in Figure 4, where the cumulative density function (CDF) of each UQ method was approximated with the built-in kernel density estimation function in MATLAB. As seen, the CDF of *Z* in (23) is almost identical to the results obtained with MC, when *p* was 3. This shows our algorithm can accurately quantify the effect of lognormal uncertainties on model responses.

#### *4.2. Example 2: Roof Structure*

The performance function of a roof truss structure as illustrated in Figure 5 was used to study the UQ accuracy of our algorithm to deal with different nonstandard uncertainties, when the system involves nonlinear polynomial functions. The performance function is described as in [19,53]:

$$Z = \mu(X) = 0.03 - 0.5ql^2 \left(\frac{3.81}{A\_2 E\_2} + \frac{1.13}{A\_1 E\_1}\right) \tag{25}$$

where *q* and *l* represent a vertical load and roof span that are used to calculate the nodal load *P* = *ql*/4, and *A* and *E* are the cross-sectional area and elastic modulus, respectively. In this case study, all parameters in the performance function in (25) are uncertain with

predefined PDFs. This yields a six-dimensional random space (*n* = 6). Details about the statistical description of uncertainties are given in Table 3.

**Figure 4.** Comparison of the cumulative density functions of *Z* calculated with the algorithm in this work and MC in Example 1. (**a**) mBiDRM-PCE and (**b**) mTriDRM-PCE.

**Figure 5.** Illustration of a roof truss structure [19,53]. Reprinted from [19], Copyright (2018), with permission from Elsevier.

**Table 3.** Details of the uncertain variables in Example 2 [19,53]. Reprinted from [19], Copyright (2018), with permission from Elsevier.


Like Example 1, Hermite polynomials were used as the basis functions, and the PCE expressions of lognormal uncertainties were constructed with normalized Gaussian variables. In addition, it is assumed that *l* follows a Weibull distribution, for which the PCE expression was approximated with standard random variables using its PDF as described in Section 3.2. Additionally, to compute the unknown PCE coefficients of the performance function *Z* in (25), the mBiDRM and mTriDRM were used to reduce the computational cost by converting the high-dimensional integrals into low-dimensional ones, which were further approximated with quadrature points constructed by the Gaussian–Hermite quadrature rules. Note that the number of quadrature points for each dimension *θ<sup>i</sup>* was set to 5 as done in Example 1. To quantify the UQ accuracy, the relative errors *<sup>R</sup>* defined in (24)

were calculated by referring to the results of MC. For MC, 107 samples were used to ensure UQ accuracy.

The results of Example 2 are summarized in Table 4. It is worth mentioning that the results from previous works (e.g., [19,53]) were not given. This is because the uncertainty in [19,53] follows a normal distribution and we intentionally introduced lognormal and Weibull distributions to test the performance of our method to deal with different types of nonstandard uncertainties. Following Example 1 in the previous section, however, the LAR-based PCE was chosen to compare the UQ accuracy with two different polynomial orders (*p* = 2 and 3). The simulation results are given in Table 4. The total number of PCE terms for our algorithm is 28 and 84, when the polynomial order for each uncertain parameter was set to 2 and 3, respectively, while the number of PCE terms of the LAR-based PCE is 28 and 77 for each polynomial order, respectively.


**Table 4.** Summary of statistical moments and their relative errors in Example 2.

Similar to Example 1, it was found that the UQ accuracy can be affected by the polynomial order (*p*) and the number of random variables in the gDRM (*M*). It is important to note that the mean values (*μZ*) of different PCE-based gDRM methods in Table 4 reached the same result (0.0064), since only two significant figures were shown. Specifically, as the number of *p* and *M* increases, the UQ accuracy can be improved. For example, when *M* was set to 3 for the mTriDRM and when *p* was set to 3, the relative error *<sup>R</sup>* of *κ<sup>Z</sup>* is one order of magnitude smaller than other methods (~0.0202%). It was also found that when the polynomial order *p* was set to 2, the UQ accuracy of LAR-based PCE is close to our method (mBiDRM and mTriDRM) for all four statistical moments. When *p* was set to 3, the UQ accuracy of mTriDRM-PCE has the same order of magnitude as the LAR-based PCE. This clearly shows the algorithm in this work can deal with nonlinear problems with nonstandard uncertainties and provide accurate results as one of the most popular techniques in the literature.

For comparison, Figure 6 shows the CDFs of the performance function *Z* in (25), with different combinations of *p* and *M*. As seen, the UQ results of our algorithm converge to MC, as the number of *p* and *M* increases. It is important to note that the accuracy of MC is related to the total number of samples used for UQ. When the model is highly nonlinear and when the number of uncertainties is large, a large number of samples is required, thereby increasing the computational cost for MC. This is further discussed with two more complicated examples below.

**Figure 6.** Comparison of the cumulative density functions of *Z* calculated with the algorithm in this work and with MC in Example 2. (**a**) mBiDRM-PCE and (**b**) mTriDRM-PCE.

#### *4.3. Example 3: Truss Structure with 13 Members*

A 13-bar truss structure [48] was chosen to study the performance of the UQ algorithm for dealing with a combination of different types of uncertainties. A schematic of the truss structure is shown in Figure 7, which involves 8 nodes and 13 members.

**Figure 7.** Illustration of a 13-bar truss structure [48]. Reprinted from [48], Copyright (2019), with permission from Elsevier.

In this example, it is assumed that there are three external loads, i.e., *P*1, *P*2, and *P*3, imposed on nodes 6, 7, and 8, respectively, which follow a normal distribution. Additionally, the elastance *E* and the sectional area *A* for each member are assumed to be independent and lognormally distributed. The statistical description of each parameter is listed in Table 5. Finite element analysis was used to quantify uncertainty in the performance function *Z* of the 13-bar truss structure, which is mathematically defined as in [48]:

$$Z = \mathfrak{u}(X) = h\_{\max} - h \tag{26}$$

where *hmax* is the threshold describing the maximum allowable deflection, i.e., displacement, and *h* is the vertical displacement on node 3. In reliability analysis, when the performance function *Z* exceeds a limit, e.g., zero in (26), it would be considered as a failure event. Additionally, the probability of failure events is referred to as the failure probability, which is often used for structural reliability analysis [54]. Thus, it is essential to calculate the statistics of the prediction *h* in a precise and computationally efficient way. Here we focus on evaluating the displacement *h* with different UQ methods. The proposed methods in this work were integrated with finite element analysis to assess the precise prediction *h*, and the performance was validated with other UQ techniques, such as MC, LAR-based PCE, and the method in [48].


**Table 5.** Details of the uncertain variables used in Example 3 [48]. Reprinted from [48], Copyright (2019), with permission from Elsevier.

As seen in Table 5, there are two lognormal uncertainties and three normal uncertainties. For the lognormal uncertainty, the normalized Gaussian variables and Hermite basis functions were used to build the PCE surrogate models in (13), while for normal uncertainty, the formulation of PCE models is defined in (1). Once these PCE models of uncertainties are available, the resulting surrogate model of *h* in (26) is described with unknown PCE coefficients that can be solved with the mgDRM involving five quadrature points in each dimension (i.e., *θ<sup>i</sup>* = 5). Two different polynomial orders (2 vs. 3) and two different values of *M* (2 vs. 3) were considered. Table 6 summarizes the simulation results. It is important to note that the total number of samples for each parametric uncertainty in MC was set to 107 to ensure UQ accuracy.

**Table 6.** Summary of statistical moments and their relative errors in Example 3.


For the implementation of LAR-based PCE, 1000 model runs were used to evaluate the truss structure model as in Figure 7. Further, the sparse PCE terms were determined with the LAR algorithm to find the sensitive basis polynomials with a predefined value (i.e., 1.0 × <sup>10</sup><sup>−</sup>6) as done in previous examples. Note that the value was used to terminate the LAR algorithm when the maximum correlation between the residual and the basis polynomials was below the value. Details of the LAR procedures can be found in [50,51]. In addition, the resulting number of PCE terms for the LAR-based PCE was 21 and 40, when the polynomial order *p* was 2 and 3, respectively. In contrast, the total numbers of PCE coefficients of our method were 21 and 56, which were determined using (3).

As seen in Table 6, it was found that UQ accuracy can be affected by the polynomial order (*p*) and the total number of integration variables (*M*) used in the gDRM. For example, when the mTriDRM-PCE was used (*p* = 3 and *M* = 3), the relative error *<sup>R</sup>* of *κ<sup>h</sup>* is ~0.1628% in Table 6. This is smaller than the results (~0.4273%) in [48]. As compared to the LAR-based PCE, our algorithm provides comparable results. This shows the potential of the gDRM-based PCE to deal with more complicated problems, especially when estimating the higher-order statistical moments since it was previously considered challenging [55]. To graphically study the UQ accuracy with different combinations of *p* and *M*, Figure 8 shows the estimated CDFs of the displacement *h*. As seen, the mTriDRM-PCE, when *p* was 3, can provide almost identical results, as compared to MC.

**Figure 8.** Comparison of the cumulative density function of the displacement *h* calculated with the algorithm in this work and with MC in Example 3. (**a**) mBiDRM-PCE and (**b**) mTriDRM-PCE.

As compared to the first two examples, UQ in this case study requires integrating the finite element analysis with the PCE, which may increase the computational burden for UQ. We further studied the computational time to calculate the displacement *h* with an office desktop (Core i5-8400 central processing unit (CPU) at 2.80 GHz). Using MC with 107 samples for each uncertainty, it was found that ~31.03 min were required to simulate all nodes in Figure 7. For the mBiDRM-PCE, it took ~21.36 and 48.06 s to simulate all nodes, when the polynomial order (*p*) was set to 2 and 3, respectively. In addition, for the mTriDRM-PCE, it was found that ~26.28 and 67.31 s were required, when *p* was 2 and 3, respectively. The mTriDRM-PCE took longer than mBiDRM-PCE due to the larger number of integrals that need to be solved. Further, depending on the polynomial order *p*, the requisite number of terms in PCE to approximate the displacement *h* differs as shown in (3), which leads to the difference in computational times of each method (e.g., mBiDRMand mTriDRM-PCE).

#### *4.4. Example 4: Planar Truss Structure with 23 Members*

A planar truss structure with 23 members as in Figure 9 is considered in this case study. It is assumed that all parameters cannot be known with certainty, resulting in a ninedimensional random space. This allows us to validate the UQ accuracy of the proposed method for dealing with many uncertainties. Details of model parameters are given in Table 7.

**Figure 9.** Illustration of a planar truss structure with 23 members [19,56]. Reprinted from [56], Copyright (2006), with permission from Elsevier.


**Table 7.** Details of the uncertain variables in Example 4 [19,56]. Reprinted from [56], Copyright (2006), with permission from Elsevier.

In this case study, *E* is the elastic modulus; *A*<sup>1</sup> and *A*<sup>2</sup> are the cross-sectional area corresponding to the specific member; and *P*<sup>1</sup> to *P*<sup>6</sup> represent external loads imposed on nodes from 8 to 13. As done in Example 3, the finite-element analysis was used to calculate the performance function *Z* as in [19,56]:

$$Z = \mathfrak{u}(X) = v\_{\max} - v \tag{27}$$

where *vmax* represents a specific threshold, i.e., the maximum deflection, and *v* is the vertical displacement on node 4 that will be approximated with different UQ methods. Our objective in this example is to accurately approximate uncertainty in *v* such that the failure probability can be quantified for structural reliability analysis.

As seen in Table 7, three parameters (*E*, *A*1, and *A*2) follow a lognormal distribution, while the rest of the parameters follow a Weibull distribution. Thus, the resulting random space is nine-dimensional, i.e., *n* = 9. To derive a PCE expression of *v* and to compute the statistical moments, the finite element analysis was coupled with the proposed methods, i.e., mBiDRM- and mTriDRM-PCE, and two different polynomial orders, i.e., *p* = 2 and 3, were investigated. Note that in these methods, the BiDRM approximates a nine-dimensional integral with 36 two- and 9 one-dimensional integrals, while the TriDRM requires the evaluations of 84 three-, 36 two-, and 9 one-dimensional ones. These integrals were calculated with the Gauss–Hermite quadrature rule [44], where the number of quadrature points for each random variable was set to 5. In addition, the number of samples for MC was set to 10<sup>7</sup> in order to compute the relative errors *<sup>R</sup>* in (24) for comparison purposes. The UQ results and relative errors (*R*) are shown in Table 8. Since the distribution of uncertainties is different as compared to [19,56], the results in these works were not given. However, LAR-based PCE was simulated for two different polynomial orders (*p* = 2 and 3) for comparison. The simulation results are summarized in Table 8. For the LAR-based PCE, the total number of PCE terms was 55 and 212, when *p* was set to 2 and 3, respectively. In constrast, 55 and 220 PCE coefficients were used in our algorithm for the output response, when *p* was set to 2 and 3, respectively.

**Table 8.** Summary of statistical moments and their relative errors in Example 4.


Similar to previous examples, our method provides accurate results, as the number of polynomial order (*p*) and the integration variables (*M*) in gDRM increase. For example, the mTriDRM outperforms others when *p* was 3, since the relative errors are smaller as in Table 8. In addition, it was found that the mTriDRM-PCE can provide comparable results to the LAR-based PCE for all statistical moments, thus confirming the accuracy of the mTriDRM-PCE-based algorithm in this work. For comparison, the CDF of the displacement *v* on node 4 was approximated with the proposed method and is shown in Figure 10. As compared to MC, it was found that mTriDRM-PCE provided the most accurate result, when *p* was set to 3. This shows the proposed algorithm can deal with complicated problems involving many uncertainties.

**Figure 10.** Comparison of the cumulative density functions of the displacement *v* calculated with the algorithm in this work and with MC in Example 4. (**a**) mBiDRM-PCE and (**b**) mTriDRM-PCE.

We further studied the computational time of the proposed method and compared the results to MC. Similar to the previous case study, it was found that the computational time for MC is larger, as compared to the gDRM-based method. For example, ~55.54 min were required in this case study, when 107 samples were used for each parametric uncertainty for MC. However, the computational time for the mTriDRM-PCE was found to be ~11.52 min, when *p* was set to 3, which is about 80% lower than MC. It is also important to note that the gDRM-PCE based UQ method, as compared to sampling-based MC, can provide an analytical expression of the displacement *v* as a function of uncertainties. This provides mathematical explanations to gain a deep understanding of the problem for improved reliability analysis. The analytical expression, for instance, can be combined with sensitivity analysis techniques to find the most sensitive random variable that can significantly affect the system's performance.

#### **5. Conclusions**

Using the polynomial chaos expansion (PCE), an uncertainty quantification (UQ) algorithm is presented to deal with many uncertainties that follow nonstandard distributions. To build PCE-based surrogate models, standard random variables are used to identify the relationship between nonstandard and standard distributions. To reduce the computational cost for calculating the PCE coefficients of model outputs, a generalized dimension reduction method is used to transform a high-dimensional integral resulting from the spectral projection (SP) into a few lower-dimensional integrals, which can be rapidly solved with quadrature rules in real-time. To show the UQ accuracy of our algorithms, four examples of structural reliability analysis were used. As compared to Monte Carlo simulations and other works in the literature, our results show the superior performance of the algorithms in terms of UQ accuracy and computational time. This shows the potential of the algorithm to tackle UQ in more complicated engineering problems that require consideration of many uncertainties. However, when uncertainty has a nonstandard distribution with a large variance, more PCE coefficients of uncertainty might be needed in the proposed approach to ensure UQ accuracy. In this case, it can lead to intensive computational burden in the presence of many uncertainties. Future work will improve the proposed algorithm to address this issue by identifying orthogonal polynomial functions only for a given uncertainty.

**Author Contributions:** Conceptualization, Y.D.; methodology, Y.D. and J.S.; validation, J.S.; investigation, J.S.; writing—original draft preparation, J.S.; writing—review and editing, Y.D.; supervision, Y.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the National Science Foundation, Division of Civil, Mechanical and Manufacturing Innovation (CMMI), under the award No. 1727487.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


**Vicente Borja-Jaimes 1,\*, Manuel Adam-Medina 1, Betty Yolanda López-Zapata 2, Luis Gerardo Vela Valdés 1, Luisana Claudio Pachecano <sup>1</sup> and Eduardo Mael Sánchez Coronado <sup>3</sup>**


**Abstract:** A fault detection and isolation (FDI) approach based on nonlinear sliding mode observers for a wind turbine model is presented. Problems surrounding pitch and drive train system FDI are addressed. This topic has generated great interest because the early detection of faults in these components allows avoiding irreparable damage in wind turbines. A fault diagnosis strategy using nonlinear sliding mode observer banks is proposed due to its ability to handle model uncertainties and external disturbances. Unlike the reported solutions, the solution approach does not need a priori knowledge of the faults and considers system uncertainty. The robustness to disturbances, uncertainties, and measurement noise is shown in the dynamic of the generated residuals, which is sensible to only one kind of fault. To show the effectiveness of the proposed FDI approach, numerical examples based on a wind turbine benchmark model, considering closed loop applications, are presented.

**Keywords:** fault detection and isolation (FDI); sliding mode observer and wind turbine; nonlinear systems

#### **1. Introduction**

Wind energy has been one of the fastest growing renewable energy sources in recent years [1]. Recently, alternative research has emerged, known as hybrid wind and photovoltaic (PV)-based systems, which provide information on the improvements of these two activities together, as described in [2]. Therefore, development of new wind turbines and wind energy conversion systems, to increase the efficiency and reliability of the energy conversion process, has become one of the main activities (and challenges) of the power industry, as well as the design of new approaches for fault diagnosis. From the reliability analysis carried out by leading countries in the production of wind energy, it is known that most of the faults that occur in wind turbines are located in the pitch system, drive train, and the control system [3]. Moreover, other faults that commonly occur are located on input and output sensors of the different subsystems that make up the wind energy conversion system [4].

In addition to their susceptibility to faults, wind turbines present significant control challenges due to the fact that they include several subsystems and due to the non-linear nature of their input [5]. Therefore, there are different approaches that describe the wind turbine model, such as in [6], where a method is proposed to analyze the parameters in the propeller design, such as drag force, lift force, thrust coefficient, power coefficient, and lift coefficient, among others. These mathematical formulations used in the QBlade software were based on the blade element moment (BEM) theory. In the work [7], specific

**Citation:** Borja-Jaimes, V.; Adam-Medina, M.; López-Zapata, B.Y.; Vela Valdés, L.G.; Claudio Pachecano, L.; Sánchez Coronado, E.M. Sliding Mode Observer-Based Fault Detection and Isolation Approach for a Wind Turbine Benchmark. *Processes* **2022**, *10*, 54. https://doi.org/10.3390/pr10010054

Academic Editors: Francisco Ronay López-Estrada and Guillermo Valencia-Palomo

Received: 15 November 2021 Accepted: 8 December 2021 Published: 28 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

airfoil simulations were performed for wind turbines of less than 1 kW, considering that the nominal wind speed was 8.4 m/s. In the presented simulations, lift and drag coefficients for different angles of attack are compared. Furthermore, a blade design and performance analysis were published in [8], where the main objective was to optimize the number of blades and tip speed ratio. A frequency regulation model for a wind farm and a thermal power plant was established in [9], where the main variable to be controlled was the load frequency of the produced power. To obtain an equivalent model of a wind farm, different wind turbines were compared, computing the gap metrics of 11 models under different operating conditions. One way to avoid wind turbine failures is through early fault detection methods. An FDI-based approach to dynamic systems, where faults can be detected and located, is described in [10]. This is based on data-driven, intelligent, and analytical modeling methods. In the literature, data-based diagnostic systems use multi-sensors that are a set of signals with fusion technology, which allows estimating, combining, and correlating the data obtained from sensors or information sources [11]. In the sense of the use of multi-sensors, but applied to fault diagnosis in engines, there is a proposal using infrared thermal imaging data, and vibration measurements for automatic fault detection [12]. In the multisensory integration to perform FDI in rotating machines, there is work using data fusion and artificial neural networks [13], FDI was also done using learning algorithms limited to rotor faults [14], the latter later proposed another work based on the same fault detection strategy, but on the gearbox [15]. Multi-sensor techniques represent an interesting and novel alternative to perform FDI in wind turbines.

On the other hand, model-based methods for fault diagnosis in wind turbines present good results, such as in [16], where the diagnosis is carried out using Kalman filters for residual generation. The model-based approach is more frequent in the literature than the data-driven approach, since it is possible to use existing mathematical methods based on differential equations to design control approaches based on fault diagnosis for wind turbine systems [17]. The methods proposed in the literature generally use thresholds for the residual values generated by the difference between the outputs of observer banks and the output signals of the system, as shown in [18]. Fault detection in wind turbine systems is reported in [19] and fault diagnosis based on state estimation is presented in [20]; these works describe the detection and isolation of faults in sensors and actuators of reference wind turbines. An approach to detect faults in actuators and sensors of wind turbines using unknown input observers (UIOs), is presented in [21]. Fault detection is carried out using residual signals and fault localization is achieved by estimating states, output signals, and controlled input signals without considering step faults. In [22], an adaptive diagnosis based on an observer is designed to perform FDI in wind turbines. In [23], three different FDD approaches are addressed, in the power train subsystem, in the generator, and in the wind turbine converter. A robust fault estimation scheme is presented in [24], where actuator and sensor faults in the pitch and drive train system are considered. Robustness is achieved by decoupling the faulty dynamics from the system using a series of coordinate transformations. For the non-faulty dynamics, an observed reduced order input is used for state estimation. While for dynamic faults, an accommodation approach based on sliding mode observers is presented. In [25], a wind speed estimator is developed to remedy the inaccuracies in the single point measurement in the WECS gondola by using the discrete Kalman filter to perform the aerodynamic torque estimation and subsequently calculate the effective wind speed. The design procedure, calculation, and estimation procedure are repeated by implementing an extended Kalman filter (EKF). The results of both algorithms are compared by simulation.

Sliding mode designs for observer and/or control are well known for their robustness properties [26]. Therefore, in [27], an approach based on adaptive sliding mode observers in wind turbines is proposed to perform fault diagnosis as part of an adaptive fault tolerant control (AFTC). In this approach, the adaptive sliding mode observer is used to simultaneously detect both the actuator and sensor faults. In the FTC design, model uncertainty and wind speed variations are considered as unknown disturbances. This allows accurate state

estimation and reconstruction of the disturbances. However, the proposed scheme does not consider the faults in the pitch and drive train.

From the literature review, it can be observed that, to solve the problem of detecting faults in the pitch system and handling dynamics introduced by faults, most of the authors have presented methods of compensating for the dynamics introduced by faults. However, the main drawback of the methods presented is that a priori knowledge of the magnitude of the faults is required. On the other hand, the problem of the lack of precision in the measurement of the rotor torque has been addressed through estimation approaches, which do not consider the uncertainty present in the estimation and the system itself.

In this work, a benchmark model is used, which consists mainly of three parts: the blade and pitch system, the generator and converter, as well as the drive train, for which a three-bladed horizontal wind turbine coupled to an electric generator with a capacity of 4.8 MW is considered. The model allows proposing faults in the actuators, sensors, and systems in the pitch actuators, the drive train, as well as in the converter system. In addition, it offers the possibility to analyze and study the operation with different types of structures for the detection and localization of faults in the wind turbine model in a more realistic context [28]. A fault diagnosis approach is proposed on this model. The main contribution of this work is the design of a fault diagnosis for the pitch system and the drive train of a wind turbine model-based. The novelty of the proposed approach is that the problem of the dynamics introduced by the faults, and the uncertainty associated with the system, in a simple way, is addressed, and eliminates the need to know a priori the behavior of the faults. Fault detection and isolation is done using the concept of residual signals and fault signatures. For the residual generation, a bank of sliding mode observers is proposed, which are designed to be sensitive only to one type of fault and to be robust to both, the dynamics of the faults and system uncertainty. This is accomplished through proper handling of the observer's error injection term and a transformation of the model in which the error introduced by the faults and the uncertainty of the system act on the input signal. The results obtained in simulation show the effectiveness of the proposed scheme.

The document is organized as follows: in Section 2 we present the wind turbine model, sliding mode observer design, and FDI scheme based on sliding mode observers, Section 3 describes the results obtained from numeric examples; Section 4 is a discussion about the results analyzed and finally the conclusion is presented.

#### **2. Materials and Methods**

#### *2.1. Wind Turbine Model*

The nonlinear wind turbine model considered in this paper is the benchmark described in [28]. This benchmark model considers the wind turbine at the system level as shown in Figure 1. This figure shows the interconnection of the subsystems that make up the wind turbine. The input to the system is the wind speed *vω*, which causes the turbine blades to rotate. This motion propagates through the rotor at a speed *ωr*, and a torque *τr*. The drive train couples this motion with the shaft of the generator, which rotates at a speed *ωg*, and has a torque *τg*. The generator is in charge of transforming the mechanical energy into electrical energy. The electrical power generated by the system is controlled by modifying the turbine aerodynamics through pitch control or by controlling the rotational speed of the generator (torque control). In either case, the objective of the control system is to track the reference electric power *Pref* , or, if not possible at least, to minimize the tracking error. The pitch control block in Figure 1 receives as inputs the measurements of the pitch angle *βi*, the generator speed *ωg*, and the generated electrical power *Pg*. The output of the block is the pitch position reference signal *βref* , which will allow the capture of the kinetic energy needed to produce the desired electrical power output. On the other hand, the torque control block has as input the generator speed measurement and the generated electrical power. The controller generates as output the reference torque *τref* , which will allow us to reach the desired electrical power output.

**Figure 1.** Interconnection of subsystems of the wind turbine in benchmark.

The mathematical model of each component of the wind turbine shown in Figure 1 is presented in the following lines.

#### 2.1.1. Blade and Pitch Model

This model is a combination of the aerodynamic model and the wind and pitch model. The wind turbine nonlinear aerodynamic model is modeled as a torque acting on the blades. This aerodynamic torque is defined as

$$
\pi\_\mathbf{r}(t) = \frac{1}{2} \,\,\pi\rho\boldsymbol{R}^3 v\_\omega^2 \,\,(t) \,\mathbf{S} \,\,(\lambda, \beta) \tag{1}
$$

where *S*(*λ*, *β*) is the power coefficient, *R* is the rotor ratio, *ρ* is the density of the air, and *vw*(*t*) is the wind speed.

The hydraulic pitch system is modeled as a second order transfer function

$$\frac{\beta\_l(s)}{\beta\_{ref}(s)} = \frac{\psi^2}{s^2 + 2\zeta\psi\,s + \psi^2} \tag{2}$$

where *β<sup>i</sup>* is the *i th* pitch system whose position is measured with two sensors *βi*,*m*<sup>1</sup> and *βi*,*m*2; *βref*(*t*) is the reference value, *ζ* is the damping factor, and *ψ* is the natural frequency.

#### 2.1.2. Drive Train Model

The drive train system is modeled by a simple two-mass model. Therefore, a drive train model can be represented by

$$\begin{aligned} J\_r \dot{\omega}\_r(t) &= \tau\_r(t) - K\_{\text{d}t} \theta - (H\_{\text{d}t} - F\_r) \omega\_r(t) + \frac{H\_{\text{d}t}}{N\_{\text{g}}} \omega\_{\text{g}}(t) \\ J\_{\text{g}} \dot{\omega}\_{\text{g}}(t) &= \frac{(\eta\_{\text{d}t} K\_{\text{d}t})}{N\_{\text{g}}} \theta + \frac{(\eta\_{\text{d}t} H\_{\text{d}t})}{N\_{\text{g}}} \omega\_{\text{g}}(t) - \left(\frac{(\eta\_{\text{d}t} H\_{\text{d}t})}{N\_{\text{g}}^2} + F\_{\text{g}}\right) \omega\_{\text{g}}(t) - \tau\_{\text{g}}(t) \\ \dot{\theta}\ (t) &= \omega\_{r}(t) - \frac{1}{N\_{\text{g}}} \omega\_{\text{g}}(t) \end{aligned} \tag{3}$$

where *ω<sup>r</sup>* and *ω<sup>g</sup>* are the rotor and generator speed, both measured by two sensors *ωr*,*mi* and *ωg*,*mi* (with *i* = 1, 2). Moreover, *θ* denotes the torsion angle of the drive train. *τ<sup>g</sup>* and *τ<sup>r</sup>* are the converter and aerodynamic torques. *Jr* y *Jg* are the moment of inertia of the lowand high-speed shafts, respectively. *Fr* and *Fg* are the viscous frictions of the low- and high-speed shafts, respectively. *Kdt*. is the torsion stiffness of the drive train, *Hdt* is the torsion damping coefficient of the drive train, *Ng* is the gear ratio and *ηdt* is the efficiency of the drive train.

#### 2.1.3. Generator and Converter Model

The generator and converter dynamics can be modeled by a first–order transfer function:

$$\frac{\pi\_{\pounds}(s)}{\pi\_{\parrow\reflectbox{g}}(s)} = \frac{\gamma}{s+\gamma} \tag{4}$$

where *τg*,*ref* is the reference value and *γ* is the cutoff frequency. The power produced by the generator is given by

$$P\_{\mathcal{K}}(t) = \eta\_{\mathcal{K}} \omega\_{\mathcal{K}}(t) \tau\_{\mathcal{K}}(t) \tag{5}$$

The efficiency of the generator is determined by *η<sup>g</sup>* and *Pg* is the power produced by the generator.

#### *2.2. Sliding Mode Observer*

The objective of an observer is to estimate the states of a system, which are generally not measurable. In essence, an observer is a mathematical model of the system, which receives the inputs of the system and the estimation error obtained from the difference between the output of the system and the output of the observer himself.

The simplest form of an observer is the Luenberger observer, in which the estimation error is fed back linearly through a gain. However, the Luenberger observer is very sensitive to measurement noise and the success of the estimation is highly dependent on the accuracy of the system model. Therefore, if the model presents some uncertainty, or is affected by external disturbances, the estimation results will be very poor. To overcome this drawback, the design of a type of non-linear observer is proposed, known as a sliding mode observer. This observer shows to be insensitive to parametric variations and presents robustness to disturbances, which represents a great advantage over other types of observers [29]. In particular, in this work, these properties are used to obtain accurate estimates of the measurements made by the sensors, even in the presence of system faults and uncertainty.

To appreciate the properties of the sliding mode observers, namely, their robustness to modeling uncertainties and disturbances, the design considers a system with uncertainty with representation in the state space given by:

$$\begin{array}{l}\dot{\mathbf{x}}(t) = A\mathbf{x}(t) + Bu(t) + D\tilde{\xi}(t, \mathbf{x}, \mathbf{u})\\y(t) = Cx(t)\end{array} \tag{6}$$

where *<sup>A</sup>* ∈ *<sup>R</sup>n*×*n*, *<sup>B</sup>* ∈ *<sup>R</sup>n*×*m*, *<sup>C</sup>* ∈ *<sup>R</sup>p*×*n*, and *<sup>D</sup>* ∈ *<sup>R</sup>n*×*<sup>q</sup>* with *<sup>p</sup>* ≥ *<sup>q</sup>*. Assume that the matrices *<sup>B</sup>*, *<sup>C</sup>* are the full rank and the function *<sup>ξ</sup>* : *<sup>R</sup>*<sup>+</sup> × *<sup>R</sup><sup>n</sup>* × *<sup>R</sup><sup>m</sup>* → *<sup>R</sup><sup>q</sup>* is unknown, but bounded so that *<sup>ξ</sup>*(*t*, *<sup>x</sup>*, *<sup>u</sup>*) ≤ *<sup>r</sup>*<sup>1</sup> *<sup>u</sup>* +*α*(*t*, *<sup>y</sup>*) and *<sup>α</sup>* : *<sup>R</sup>*<sup>+</sup> × *<sup>R</sup><sup>p</sup>* → *<sup>R</sup>*<sup>+</sup> .

Suppose that exits a liner change of coordinates *T* so that the system can be written as

$$\begin{array}{c} \dot{\mathbf{x}}\_1(t) = A\_{11}\mathbf{x}\_1(t) + A\_{12}\mathbf{x}\_2(t) + B\_1\mathbf{u}(t) \\ \dot{\mathbf{x}}\_2(t) = A\_{21}\mathbf{x}\_1(t) + A\_{22}\mathbf{x}\_2(t) + B\_2\mathbf{u}(t) + D\_2 \not\subset \\ \mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) \end{array} \tag{7}$$

where *<sup>x</sup>*<sup>1</sup> <sup>∈</sup> *<sup>R</sup>*(*n*−*p*), *<sup>x</sup>*<sup>2</sup> <sup>∈</sup> *<sup>R</sup>p*, and the matrix *<sup>A</sup>*<sup>11</sup> has stable eigenvalues. The observer structure that will be considered can be written in form

$$\begin{aligned} \dot{\mathfrak{k}}\_1(t) &= A\_{11}\mathfrak{k}\_1(t) + A\_{12}\mathfrak{k}\_2(t) + B\_1 u(t) - A\_{12}e\_y(t) \\ \dot{\mathfrak{k}}\_2(t) &= A\_{21}\mathfrak{k}\_1(t) + A\_{22}\mathfrak{k}\_2(t) + B\_2 u(t) - (A\_{22} - A\_{22}^s)e\_y(t) + w \\ \mathfrak{j}(t) &= \mathfrak{k}\_2(t) \end{aligned} \tag{8}$$

where *A<sup>s</sup>* <sup>22</sup> is a stable design matrix, *ey*(*t*) = *y*ˆ(*t*) − *y*(*t*) and the discontinuous vector *w* is defined as

$$w = \begin{cases} \ -\sigma \parallel D\_2 \parallel \frac{P\_2 \epsilon\_y(t)}{\|P\_2 \epsilon\_y(t)\|} \text{ if } \epsilon\_{\mathcal{Y}} \neq 0\\ \ 0 \text{ otherwise} \end{cases} \tag{9}$$

where *<sup>P</sup>*<sup>2</sup> ∈ *<sup>R</sup>p*×*<sup>p</sup>* is a symmetric positive defined Lyapunov matrix for *<sup>A</sup><sup>s</sup>* 22.

The scalar *σ* is chosen so that *ξ* < *σ*. If the state estimation error is defined as *e*1(*t*) = *x*ˆ1(*t*) − *x*1(*t*), then it is straightforward to show

$$\begin{array}{c} \dot{e}\_1(t) = A\_{11}e\_1(t) \\ \dot{e}\_\mathcal{Y}(t) = A\_{21}e\_1(t) + A\_{22}^s e\_1(t) + w - D\_2 \not\subset \end{array} \tag{10}$$

It is shown in [30] that the error system Equation (10) is quadratically stable and a sliding motion takes place, forcing *ey* (*t*) = 0 in finite time.

If *x*ˆ(*t*) represents the state estimate for *x*(*t*) and *e*(*t*) = *x*ˆ(*t*) − *x*(*t*), it follows that if

$$G\_l = T^{-1} \begin{bmatrix} A\_{12} \\ A\_{22} - A\_{22}^s \end{bmatrix} \text{ and } G\_{nl} = [D\_2] T^{-1} \begin{bmatrix} 0 \\ I\_p \end{bmatrix} \tag{11}$$

then the robust observer given in Equation (8) can be written in terms of the original coordinates in the form

$$
\dot{\hat{x}}(t) = A\mathbf{x}(t) + B\boldsymbol{u}(t) - \mathbf{G}\_l \mathbf{C}e(t) + \mathbf{G}\_{nl}w \tag{12}
$$

where the nonlinear function is given by

$$w = \begin{cases} \begin{array}{c} -\sigma \frac{P\_2 \mathbb{C}c(t)}{\|P\_2 \mathbb{C}c(t)\|} \text{ if } \mathbb{C}c \neq 0\\ 0 \text{ otherwise} \end{array} \end{cases} \tag{13}$$

It is important to note that, even when the model considered for the design of the sliding mode observer has an uncertainty term, the dynamics of the observer does not depend on the uncertainty; therefore, the convergence of the observer is achieved even in the presence of uncertainty.

The methodology for the design of the sliding mode observer described in this section is applied in the following sections for the development of fault diagnosis schemes in the pitch and drive train system.

#### *2.3. FDI Scheme Based on Sliding Mode Observers*

The principal purpose of the FDI scheme is to generate an alarm when a fault develops in the system and to identify the source and location. To detect and isolate faults on the wind turbine, a model based FDI is proposed. To achieve fault detection and isolation in the model-based approach with the observer-based method, some authors have proposed particular configurations of observer banks. In this work, the proposed scheme is based on the configuration presented in [31], and is shown in Figure 2. In this scheme, a set of observers dedicated to isolate *p* sensor faults are designed, the proposed configuration allows observer *p* to use all inputs and only output *p*. In this way, observer *p* is sensitive to sensor *p* faults. The residue generator block generates a set of structured residues from the difference of the measured variables and those estimated by the observer bank. In the decision block, the decision is made to detect the fault if the residual exceeds a predefined threshold.

Since the wind turbine system is subject to modeling uncertainty and unknown disturbance, to achieve the effective fault detection, a robust fault detection scheme is needed. Since the sliding mode techniques are well known for their robustness properties, the bank of the observer is designed using these techniques.

The novelty of the proposed scheme is that it does not require a priori knowledge of the failures and considers the uncertainty of the system. This is achieved by treating the dynamics of the faults in the pitch system and the uncertainty in the estimation of the rotor torque, as coupled disturbances acting on the input signal. Under this condition, it is possible to observe the characteristics of the sliding mode observer. The faults considered in this paper cover the pitch systems and drive train, and these are independent. Therefore, an FDI scheme was designed separately for each one.

The configuration proposed for the development of the fault diagnosis scheme in the pitch system is presented below, considering the general structure described in Figure 2.

#### 2.3.1. FDI Architecture for Pitch System

In this section, the FDI scheme for each pitch position system is presented. It considers only sensor faults, which are characterized by either a fixed value or a change gain factor on the sensor measurements. Since all three pitch positions are measured with two sensors, the aim of the FDI scheme is to detect and isolate faults in one of the two sensors. In order to perform the detection and isolation of faults in the pitch system, the architecture shown in Figure 3 is proposed. The figure shows the particular configuration of the observer bank used to estimate the pitch angle measurements, which is measured with two sensors for reasons due to hardware redundancy. Since the pitch system is integrated by three identical subsystems, there are a total of six sensors. The configuration of the observer bank allows the observer in charge of estimating the output 1 of pitch *i*, i.e., observer *i*, 1, to use only the measurement of sensor 1 of pitch *i*. In this way, the observer is sensitive only to faults in that sensor. The residual generator block generates a set of structured residues from the difference between the variables measured and those estimated by the observers. In the decision block, an alarm is generated if any of the residuals exceeds a pre-established threshold, which is established from the behavior of the residuals before and after the fault.

**Figure 3.** Architecture for FDI in the pitch system.

The FDI design only considers that one fault occurs at a time and, therefore, two fault scenarios for each pitch position are considered. The first assumes that sensor 1 is faulty and sensor 2 is healthy. The second assumes sensor 2 is faulty and sensor 2 is healthy. Under this consideration, two residual signals are needed for each pitch position system. Because residual signals are generated from the difference between the estimates of the

measurement and real measurement, a set of two sliding mode observers for each pitch position system is needed.

Since sensor faults directly influence the pitch system through the feedback, in order to consider this influence, the function *ξ<sup>i</sup> pb* is used to model the position error introduced by the sensor fault. Then, the continuous-time state space model of *i* th pitch position is as follows: . *i*

$$\begin{array}{l} \dot{\boldsymbol{x}}\_{pb}^{i}(t) = \boldsymbol{A}\_{pb}^{i}\boldsymbol{x}\_{pb}^{i}(t) + \boldsymbol{B}\_{pb}^{i}\boldsymbol{\beta}\_{ref}^{i}(t) + \boldsymbol{D}\_{}^{\pi i}\boldsymbol{y}\_{pb} \\ \boldsymbol{y}\_{pb}^{i}(t) = \boldsymbol{C}\_{pb}^{i}\boldsymbol{x}\_{p}^{i}(t) \end{array} \tag{14}$$

where *x<sup>i</sup> pb*(*t*) =  . *βi* (*t*), *β<sup>i</sup>* ] *T* . *Apb*, *Bpb* and *Cpb* are the system matrices and *ξ<sup>i</sup> pb* represents the position error introduced by the sensor fault given by *ξ<sup>i</sup> pb* = *β<sup>i</sup>* − 0.5*i*(*βi*,*m*<sup>1</sup> + *βi*,*m*2) with *i* = 1, 2, 3.

For system (14), we used the sliding mode observer described in Section 2.2. From Equation (12), two sliding mode observers for the *i* th pitch position can be designed as:

$$\dot{\mathfrak{X}}\_{pb1}^{\dot{i}}(t) = A\_{pb}^{i} \mathfrak{X}\_{pb1}^{i}(t) + B\_{pb}^{i} \mathfrak{X}\_{ref}^{i}(t) - G\_{I}^{i} \mathfrak{e}\_{y1}^{i} + G\_{nl}^{i} w\_{1}^{i} \tag{15}$$

$$\mathfrak{Y}\_{pb1}^{i}(t) = \mathfrak{F}\_{i,m1}(t)$$

$$\begin{array}{c} \dot{\mathfrak{X}}\_{pb2}^{i}(t) = A\_{pb}^{i}\mathfrak{X}\_{pb2}^{i}(t) + B\_{pb}^{i}\mathfrak{X}\_{ref}^{i}(t) - G\_{I}^{i}e\_{y2}^{i} + G\_{nl}^{i}w\_{2}^{i} \\ \mathfrak{Y}\_{pb2}^{i}(t) = \mathfrak{F}\_{i,m2}^{i}(t) \end{array} \tag{16}$$

where *x*ˆ*<sup>i</sup> pb*(*t*) and *<sup>y</sup>*ˆ*<sup>i</sup> pb* denote, respectively, the estimate of state and the estimate of output. *Gi <sup>l</sup>* and *<sup>G</sup><sup>i</sup> nl* are appropriate gain matrices to guarantee the stability and convergence of the dynamics of error. The estimation error is defined as *e<sup>i</sup> <sup>y</sup>*<sup>1</sup> <sup>=</sup> *<sup>β</sup>*<sup>ˆ</sup> *<sup>i</sup>*,*m*1(*t*) − *βi*,*m*1(*t*) and *ei <sup>y</sup>*<sup>2</sup> <sup>=</sup> *<sup>β</sup>*<sup>ˆ</sup> *<sup>i</sup>*,*m*2(*t*) − *βi*,*m*2(*t*). The *w* represents a discontinuous switched component to induce a sliding motion.

The observer proposed has the capacity to handle the influences of a fault if *ξ* < *σ*, where *σ* represents the gain in the discontinues injection term *w* and is obtained from Equation (12).

Under this consideration, the sliding mode observer proposed can estimate the state accurately, even in the presence of sensor faults; therefore, it is clear that the estimated *i* th pitch position must be close to the real value before the fault occurrence and after the sensor fault. To isolate sensor faults, a structured set of residual signals is constructed to be sensitive to only one fault and is given in Table 1.


**Table 1.** Pitch system FDI logic.

The residual signal *ri*,*<sup>j</sup>* is used to indicate a fault in the sensor *j* of the *i* th-pitch position system, with *i* = 1, 2, 3 and *j* = 1, 2. Each residual is sensitive to only one fault. The details of the residual combinations, which are used to isolate different faults, are given by the fault signature.

The following subsection describes the proposed architecture for the transmission train fault diagnosis scheme. In a similar way to the scheme proposed for the pitch system, for the generation of residual signals, an observer bank configuration is proposed.

#### 2.3.2. FDI Configurations for Drive Train System

To achieve fault detection and isolation in the drive train, the architecture shown in Figure 4 is proposed, where each observer has all drive train inputs, i.e., rotor torque and generator torque, and only the corresponding sensor measurement. The drive train has hardware redundancy so that the rotor and generator speed measurement is performed with two sensors respectively. Therefore, a bank of four observers is required so that each observer is sensitive only to faults present in the sensor feeding the output signal corresponding to the observer. Since the rotor torque cannot be measured, an estimator based on an algebraic map that depends on the measurement of wind speed, pitch angle, and generator speed is used in this scheme. This algebraic map is described in more detail in [19]. In the residual generator block, the estimation errors of each of the outputs are used as residuals. Then in the decision-making block, an alarm is generated if any of the residuals exceed a preset threshold.

**Figure 4.** FDI architecture for drive train.

To ensure physical redundancy, both the generator and rotor speed are measured with two sensors. Therefore, the aim of the FDI scheme is to detect and isolate faults in one of the four sensors. Since both the rotor and generator speed are measured using encoders, the fault signals are characterized by either a fixed value or a change gain factor on the sensor measurements. The FDI design only considers that one fault occurs at a time and, therefore, two fault scenarios for both rotor and generator speed sensors are considered. The first assumes that sensor 2 is faulty and sensor 1 is healthy. The second assumes that sensor 1 is faulty and that sensor 2 is healthy. Under this consideration, four residual signals are needed to isolate a fault in the drive train system. Therefore, a set of four sliding mode observers to estimate each sensor measurement was designed.

On the other hand, since the aerodynamic torque is not possible to measure, it depends on wind speed, which can only be roughly measured, and imposes a challenge to the design of the FDI scheme. To overcome this issue, the nonlinear aerodynamics torque is obtained using an aerodynamic torque estimator based on the measured rotor speed *ωr*, pitch position *βi*, wind speed *vω*, and a mapping of the toque coefficient S (λ, β). It is important to note that due to the difficulty to accurately estimate the aerodynamics torque, an unknown bounded disturbance, *ξdt*, is added in the model. The disturbance is modeled as a bounded matching uncertainty entering the system through the input channel.

For the design of the sliding mode observers, consider the model in state space of the drive train system written as

$$\begin{array}{l} \dot{\mathbf{x}}\_{dt}(t) = A\_{dt}\mathbf{x}\_{dt}(t) + B\_{dt}\boldsymbol{u}\_{dt}(t) + D\boldsymbol{\xi}\_{dt} \\ y\_{dt}(t) = \mathbf{C}\_{dt}\mathbf{x}\_{dt}(t) \end{array} \tag{17}$$

where *xdt*(*t*) = *ωr*(*t*), *ωg*(*t*), *θ*Δ(*t*) !*T* , *<sup>u</sup>* <sup>=</sup> *<sup>τ</sup>r*(*t*) *<sup>τ</sup>g*(*t*) !*<sup>T</sup>* and *Adt*, *Bdt*, and *Cdt* are the system matrices, and *ξdt*, represents represent the system uncertainty, which are unknown but bounded.

Based on Equation (12), two sliding mode observers of FDI on the drive train are designed, in which *ωr*(*t*) and *ωg*(*t*) measurements are used separately.

The proposed sliding mode observers to estimate the rotor speed are given by:

$$\begin{aligned} \dot{\mathfrak{X}}\_{dt}^{\dot{j}}(t) &= A\_{dt}(t)\mathfrak{X}\_{dt}^{\dot{j}}(t) + B\_{dt}u\_{dt}(t) - G\_{l}^{\dot{j}}e\_{yr}^{j} + G\_{\text{nl}}^{\dot{j}}w\_{r}^{j} \\ \mathcal{Y}\_{dt}^{\dot{j}}(t) &= \hat{\omega}\_{r,mj} \end{aligned} \tag{18}$$

and the proposed sliding mode observers to estimate the generator speed are given by:

$$\begin{aligned} \dot{\mathfrak{x}}\_{dt}^{j}(t) &= A\_{dt}(t)\mathfrak{x}\_{dt}^{j}(t) + B\_{dt}u\_{dt}(t) - G\_{l}^{j}e\_{\mathcal{Y}\mathcal{S}}^{j} + G\_{n}^{j}v\_{\mathcal{S}}^{j} \\ \mathfrak{y}\_{dt}^{\text{gj}}(t) &= \hat{\omega}\_{\mathbb{S},m\text{j}} \end{aligned} \tag{19}$$

with *e j yr*(*t*) = *ω*ˆ *<sup>r</sup>*,*mj* − *ωr*,*mj*, *e j yg*(*t*) = *ω*ˆ *<sup>g</sup>*,*mj* − *ωg*,*mj*, and *j* = 1, 2.

where *ω*ˆ *<sup>r</sup>*,*mj* and *ω*ˆ *<sup>g</sup>*,*mj* are the rotor speed and generator speed estimations and are computed by observer 1 and 2, respectively. *Gl* and *Gnl* are the observer gains matrices.

The residuals of the FDI scheme are presented in Table 2. From Table 2, each residual is sensitive to only one fault. The residuals *r*4,1 and *r*4,2 indicate the rotor speed sensor faults, while the residuals *r*5,1 and *r*5,2 indicate the generator speed sensor faults.


**Table 2.** Drive train system FDI logic.

In the next section, the effectiveness of the proposed configuration to diagnose faults in both the pitch system and the transmission train is put to the test.

#### **3. Results**

In this section, the proposed FDI scheme based on sliding mode observers is tested on the benchmark model of the wind turbine.

The test only considers sensor faults in the pitch system and the drive train. The type of fault is either electrical or mechanical and may result in a fixed value or a change gain factor on the measurements. In general, the severity of these faults is low, but if these are not handled correctly, can result in a breakdown. In order to prevent this situation, early fault detection and isolation are needed.

The wind turbine is a nonlinear system, driven by a disturbance, namely the wind. The residual signals can be nonzero, even when there are no faults. Therefore, the design of the FDI scheme is expected to be robust, in order to achieve a low false alarm rate and a high fault detection rate.

To solve this challenge, the proposed FDI scheme is based on a set of sliding mode observers, robust against disturbances and model uncertainties. The effectiveness of the proposed scheme is demonstrated by a test that includes five sensor faults, three for the pitch and two for the drive train system. These faults are defined in [25] and set as:


The results of the FDI scheme for the pitch 1 system when the fault occurs are shown in Figure 5. Fault 1 occurs in sensor 1 of pitch 1, in the time interval from 2000 to 2100 s. As can be seen in the residuals plot in Figure 5, the behavior of the residuals allows identifying the presence of fault 1. It can be observed that the residual *r*1,1 presents a remarkable abrupt change during the time interval of the occurrence of fault 1, while the residual *r*1,2 remains very close to zero. Therefore, it can be concluded that a fault occurred in the pitch 1 system, specifically in sensor 1. Ideally, it would be expected that the residuals remain at zero, while there is no fault in the sensors; however, due to the presence of noise in the sensors, the residuals present small deviations around zero, even when the fault has not occurred.

**Figure 5.** Faults in the first pitch system.

For the pitch 2 system, a scaling fault is considered during a time interval of 2300–2400 s, acting on sensor 2. Figure 6 shows the behavior of the residuals generated from the observer bank in sliding modes. In the graph, it can be noticed that the residual *r*2,2 presents an erratic behavior in the fault incidence interval, while the residual *r*2,1 does not present any atypical behavior and remains at a value close to zero. Therefore, from the fault signature of the two residuals, it is shown that the fault occurred in sensor 2 of pitch 2. We should note that the magnitude of the change of the fault-sensitive residual 2 is much smaller with respect to the behavior of the fault-sensitive residual 1. This is because a scaling error fault does not have as noticeable an effect on the sensor signal as a fixed value fault does.

For the pitch 3 system, the fault occurs in sensor 1. It is a fixed value fault acting in the time interval 2600–2700 s. The results of the fault detection scheme in the pitch 3 system are presented in Figure 7. It can be observed that fault 3 can be detected by means of the residual *r*3,1, which presents a significant deviation in the time interval of the fault. While the residual *r*3,2 does not present significant changes, it remains very close to zero. From the behavior of the residuals of the pitch system, it can be observed that these were shown to be sensitive to the presence of faults, thus solving the problem of fault detection in the pitch system. Moreover, since each of the presented residuals was shown to be sensitive to only one fault, fault isolation is achieved with the help of the fault signature table.

**Figure 6.** Fault bias in the second pitch during the time 2300 to 2400 s.

**Figure 7.** Fault bias in the third pitch during the time 2600 to 2700 s.

The proposed scheme for the detection and isolation of faults in the drive train employs a bank of four observers whose function is to estimate the turbine rotor speed and the generator speed. From the difference between the estimates made by the observers and the measurements obtained by the sensors, a set of four residuals is generated. Figure 8 shows the behavior of the residuals designed to be sensitive to fault 4. Fault 4 corresponds to a scaling error in the signal of sensor 1 in charge of measuring the rotor speed. Since sensor 2 is fault free, it can be observed that the residual *r*4,2 remains practically invariant during the whole simulation interval. While the residual *r*4,1 shows an atypical behavior during the1500–1600 s interval. The magnitude of the deviation observed in the residual in the mentioned interval is high enough to infer the presence of a fault in the sensor.

Finally, the behavior of the residuals sensitive to fault 5 is shown in Figure 9. Fault 5 is a fixed value in sensor 2 responsible for measuring the generator speed, during the time interval 1000–1100 s. It is important to note that, from the point of view of fault dynamics, in a fixed value fault a more significant deviation of the residual can be observed, this compared to the case of a scaling error type fault. This is evident in the magnitude of the deviation of the residual *r*5,2. From which it can be concluded that a fault has occurred

in sensor 2 responsible for the measurement generator speed. While the behavior of the residual *r*5,1 confirms that sensor 1 is free of fault.

**Figure 8.** The abrupt fault of rotor speed sensor 1 in the time period of 1500 to 1600 s.

**Figure 9.** The abrupt fault in the speed of the generator in the period time of 1000 to 1100 s.

#### **4. Discussion**

The faults considered in this paper only cover the pitch system and drive train, and since these are modeled independent, an FDI scheme was designed separately for each one. The FDI scheme proposed for the pitch system considers the error introduced by sensor faults as an uncertainty acting on the input channel. Under this condition, a sliding mode observer bank was capable of estimating the system outputs, even in the presence of faults is proposed. The differences between the estimates made by the observers and the sensor measurements were used to regenerate a structured set of residuals for each pitch system. The structure of the FDI scheme for the pitch system is shown in Figure 3. It considers a set of six observers. To detect sensor faults, each sliding mode observer is dedicated to detect one sensor fault and generate one residual. The residual signals are generated from the measurements and estimates made by the bank of observers. After successful fault detection, the next task is to isolate the faults, which is done using a signature Table 1.

Secondly, to overcome the problem of coupled dynamics and the absence of aerodynamic torque measurement, the proposed FDI approach for the drive train employed an

aerodynamic torque estimator based on an algebraic map. The uncertainty associated with the estimation was modeled as a disturbance in the input channel of the system model. The aerodynamic torque estimate was fed to a bank of observers in sliding modes each of which employed different sensor measurements in order to be sensitive to only one fault. Each observer was dedicated to one sensor fault; consequently, the generated residual was sensitive to only one sensor fault. For the task of isolating the faults, a signature Table 2 was used.

The effectiveness of the FDI schemes is demonstrated from the evaluation of five sensor fault scenarios as described in the results section. Due to the measurement noises and the uncertainties in the system, the residual signals are not zero, even for the fault-free case. To overcome this issue, residual evaluation and threshold are employed. The thresholds are designed based on the expressions of the residuals and uncertainties and measurement noise in the system. From the results presented in the last section, it can be easily seen that the residual remains around zero before and after the occurrence of fault, but during the fault interval, the residual rises significantly. To achieve a low false alarm, the thresholds are set at the maximum absolute value of the residuals. Therefore, a fault is presented in the system if it causes significant deviations from zero on some of the residuals.

#### **5. Conclusions**

In this article, a fault detection and isolation approach is proposed for the sensors of the pitch system and drive train in a benchmark wind turbine model. This approach employs a bank of sliding mode observers, designed to receive all of the inputs of the system and only one measurement of the set of available sensors. This allows each observer to be sensitive only to faults present in the sensor input signal. The dynamics introduced by the fault in a sensor of the pitch system, as well as the uncertainty associated with the estimation of the rotor torque in the transmission train, are modeled as disturbances acting on the input channel of their respective subsystem. Under this condition, the convergence of the observers turns out to be independent of the fault dynamics and the uncertainty present in the system. The difference between observer bank outputs and the measurements of the sensors is used to generate a structured set of residues for the pitch system and the drive train. It is shown that the residuals are capable of detecting a sensor fault occurrence, and, on the other hand, they remain very close to zero when the system is fault-free. This allows decision-making on fault detection using the concept of a fixed threshold. The value of the thresholds was selected based on the behavior of the residues before and after the appearance of the fault. To achieve fault isolation, fault signature matrix constructed from the residuals was used.

**Author Contributions:** Conceptualization, V.B.-J. and M.A.-M.; methodology, V.B.-J. and B.Y.L.-Z.; software, V.B.-J. and E.M.S.C.; validation, V.B.-J. and M.A.-M.; formal analysis, M.A.-M. and V.B.-J.; investigation, V.B.-J., M.A.-M. and B.Y.L.-Z.; resources, M.A.-M., B.Y.L.-Z., E.M.S.C., L.C.P.; writing—original draft preparation, V.B.-J., B.Y.L.-Z., L.G.V.V.; writing—review and editing, M.A.-M., E.M.S.C., V.B.-J., L.G.V.V., L.C.P.; supervision, M.A.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The development of this project is the product of the support provided by the National Council of Science and Technology (CONACYT), the National Technological Institute of Mexico (TecNM)/CENIDET.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**

