Predicting the Remaining Useful Life of Turbofan Engines Using Fractional Lévy Stable Motion with Long-Range Dependence

Abstract

Remaining useful life prediction guarantees a reliable and safe operation of turbofan engines. Long-range dependence (LRD) and heavy-tailed characteristics of degradation modeling make this method advantageous for the prediction of RUL. In this study, we propose fractional Lévy stable motion for degradation modeling. First, we define fractional Lévy stable motion simulation algorithms. Then, we demonstrate the LRD and heavy-tailed property of fLsm to provide support for the model. The proposed method is validated with the C-MAPSS dataset obtained from the turbofan engine. Principle components analysis (PCA) is conducted to extract sources of variance. Experimental data show that the predictive model based on fLsm with exponential drift exhibits superior accuracy relative to the existing methods.

1. Introduction

Developing algorithms and software modules that allow processing online data from IoT devices to predict turbofan engines’ remaining useful life (RUL), i.e., the number of work cycles after the last completed cycle, is a relevant research topic. Remaining useful life prediction can effectively eliminate potential safety hazards of passenger flights [1]. The traditional approach to equipment maintenance involves a combination of breakdown and preventive maintenance [2]. The former means repairing and replacing equipment when it breaks down. The latter indicates the repairing of equipment prior to the failure, and the best timing of maintenance is chosen by the RUL prediction [3].

The degradation process of turbofan engines is related to environmental temperature, material, load changes, and other factors, making it difficult to predict the remaining useful life (RUL). Currently, the mainstream RUL prediction methods are artificial intelligence, stochastic model, and data-driven, based on acoustic monitoring of turbofan engines. The primary artificial intelligence algorithms utilized are LSTM and deep learning approaches. In recent years, various improved LSTM [4,5] and deep learning optimization algorithms [6,7,8] have been proposed to improve the accuracy of RUL prediction. These advancements also aim to combine both approaches [9,10], which necessitates a large number of data samples. The stochastic model method establishes a stochastic model for the non-stationarity and randomness of the degradation process [11,12,13]. Data-driven research often needs as many data samples as possible combined with LSTM and deep learning [14,15].

Each of these methods has its advantages and disadvantages, but the biggest problem is that the prediction length of non-stationary time series is too short. The reason for the shortcoming is that the fundamentals of these algorithms are similar to the least squares fitting methods predictive models of RUL based on artificial intelligence, which aim to establish a mapping relationship between the monitoring dataset and the failure time. However, this relationship between the monitoring data and RUL is difficult to express clearly and is always implemented as a “black box” characteristic [16].

Additionally, a method exists that involves combining turbofan engine vibration signals and acoustic signals using a multi-sensor fusion to create a two-dimensional signal image [17,18,19]. Then, we can apply deep learning to predict the RUL. Deep learning is very mature in two-dimensional feature extraction. In practical applications, the fluctuations of aircraft engines during operation are too large, and thus, we hope that the incipient faults can be identified as early as possible. However, the current deep learning models are not sensitive and accurate enough.

Degradation is a slow and non-stationary stochastic process. The degradation sequence waveform exhibits characteristics of self-similarity and long-range dependence. This means each segment of the waveform is self-similar to the overall trend of change, and the degradation state at any given moment is closely related to its previous states; this property is known as long-range dependence. Based on this concept, the objective is to establish a non-stationary stochastic prediction model that possesses both self-similar and long-range dependent characteristics. The data series generated by this model will closely mimic the actual degradation sequences observed in the turbofan engine, resulting in high prediction accuracy and extended forecasting capability.

This paper proposes a non-stationary stochastic fractional-order Levy stable operation model, demonstrating its self-similarity and long-range dependence. By analyzing the collected time series data of turbofan engine vibrations, we utilize the Hurst exponent to determine the presence of long-range dependence. Subsequently, the model parameters are identified, ensuring that the series generated by the model highly resembles the degradation process of the turbofan engine. The model is trained with the principle components analysis (PCA) output.

A review of the Brownian motion (BM)-based RUL predictive examines several problems of this approach [20]. However, there are still some problems with the BM model. The BM model is a Markov process, but in practice, the aging data exhibit LRD [21]. The fractional Brownian motion extends the BM model to the fractional order. The fBm model introduces the LRD characteristics through the Hurst parameter. Song et al. developed a series of iterative differential predictive models based on the fBm model [22,23]. Gao et al. obtained the analytical expression of the probability density function (PDF) for the RUL prediction based on fBm [24]. In fact, Levy stable motion and fractional Brownian motion were applied to different studies [23,25,26,27,28].

This article is organized as follows: Section 2 proposes the fLsm model, and Section 3 provides LRD properties and self-similar characteristics of the fLsm model. The degradation model based on fLsm is presented in Section 4, and a case study with the C-MAPSS dataset is provided in Section 5 [29].

2. Fractional Lévy Stable Motion

2.1. Probability Density Function of the Lévy Stable Distribution

The connection between the generalized central limit theorem and Lévy stable distribution implies that the analytical expression of the PDF in the closed form does not exist [30]. Therefore, the properties of the Lévy stable distribution are expressed by its characteristic function defined by $\phi \left(\theta : \alpha , \beta , \delta , \mu \right)$, which is defined as follows:

$$\phi \left(\theta :\alpha ,\beta ,\delta ,\mu \right)=E\left[e^{j\theta x}\right]=\left\{\begin{array}{c} exp\left\{j\mu \theta -\delta {\left|\theta \right|}^{\alpha }\left[1-j\beta \frac{\theta }{\left|\theta \right|}\mathit{tan}\left(\frac{\pi \alpha }{2}\right)\right]\right\}, \alpha \ne 1\\ exp\left\{j\mu \theta -\delta {\left|\theta \right|}^{\alpha }\left[1+j\beta \frac{\theta }{\left|\theta \right|}\frac{2}{\pi }\mathit{ln}\left|\theta \right|\right]\right\}, \alpha =1\end{array}\right.$$

(1)

Usually, $X~S_{\alpha }\left(\beta ,\mu ,\delta \right)$ is used to indicate that the random variable follows the Lévy stable distribution, described by the characteristic exponent $\alpha$, which measures the thickness of the tails of the distribution. $\beta$ is the skew index, $\delta$ is the drift coefficient, and $\mu$ is the diffusion parameter. The shape of the distribution is governed by parameters $\alpha$ and $\beta$. Linear transformation of the distribution is determined by $\delta$ and $\mu$. The parameters of Equation (1) must satisfy the following restrictions: $\alpha \in \left(0,2\right]$, $\beta \in [-1,1]$, $\mu \in \mathbb{R}$, $\delta >0$. The effect of the governing parameters on the PDF is shown in Figure 1.

2.2. Fractional Lévy Stable Motion

fLsm is obtained by extending the Lévy stable motion into the fractional order based on the derivation of the fBm model. The fBm model can be simply described as:

$$B_H\left(t\right)\triangleq \frac{1}{\Gamma \left(H+\frac{1}{2}\right)}\left\{{\int }_{-\mathrm{\infty }}^0\left[{\left(t-u\right)}^{H-\frac{1}{2}}-{\left(-u\right)}^{H-\frac{1}{2}}\right]\omega \left(u\right)du+{\int }_0^t{\left(t-u\right)}^{H-\frac{1}{2}}\omega \left(u\right)du\right\},$$

(2)

where $0<H<1$ is the Hurst parameter, $\omega \left(u\right)$ is the Gaussian white noise with zero expectation, and $\Gamma \left(.\right)$ is the gamma function. Let us replace $H-\frac{1}{2}$ in Equation (2) with $H-\frac{1}{\alpha }$ and the Gaussian white noise with Lévy white noise. Then, the expression of the fLsm takes the form:

$$L_{H,\alpha }\left(t\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{a\left[{\left(t-s\right)}_{+}^{H-\frac{1}{\alpha }}-{\left(-s\right)}_{+}^{H-\frac{1}{\alpha }}\right]+b\left[{\left(t-s\right)}_{-}^{H-\frac{1}{\alpha }}-{\left(-s\right)}_{-}^{H-\frac{1}{\alpha }}\right]\right\}Mds,$$

(3)

where $\left(x_{+}\right)={\left(-x\right)}_{-}=\left\{\begin{array}{c}x, x>0\\ 0, x\le 0\end{array}\right.$, a and b are real constants, $\left|a\right|+\left|b\right|>0$, and M is the Lévy stability measure in the Lebesgue measure space. Note that the fLsm model becomes the fBm model when $\alpha =2$.

2.3. Simulation Algorithm

Fractional Lévy stable noise can be expressed as follows:

$$X_{H,\alpha }=D_t^{-a}{\omega }_a\left(t\right)=\frac{1}{\Gamma \left(a\right)}{\int }_0^t{\left(t-s\right)}^{a-1}{\omega }_a\left(s\right)ds,0<a<0.5,$$

(4)

where ${\omega }_a\left(t\right)$ is the Lévy stable white noise and $X_{H,\alpha }$ denotes the generated fractional Lévy stable noise sequence. The fractional Lévy stable noise degenerates to fractional Gaussian noise when $\alpha =2$.

The fractional Lévy stable noise is the incremental process of the fLsm model, so the simulation sequence of the fLsm model can be obtained by integrating or accumulating Equation (4) as follows:

$$L_{H,\alpha }=D_t^{-1}{X}_{H,\alpha }=\frac{1}{\Gamma \left(a+1\right)}{\int }_0^t{\left(t-s\right)}^a{\omega }_a\left(s\right)ds,0<a<0.5,$$

(5)

Figure 2 shows the simulated sequences of the fLsm for H = 0.8 and $\alpha$= 1.4, 2.

3. The Property of the fLsm

The self-similar nature of the fLsm is represented by parameter H ($H\ne 1/\alpha$) [31], and its LRD characteristic is given by the relationship between α and H:αH > 1.

3.1. LRD Property

To ensure that the fLsm model effectively fits the LRD of stochastic sequences, we limit the range of parameters α and H to α ∈ (1, 2) and H ∈ (1/2, 1), respectively.

$$L_{H,\alpha }\left(t\right)={\int }_{-\infty }^{\infty }\left\{a\left[{\left(t-s\right)}_{+}^{H-\frac{1}{\alpha }}-{\left(-s\right)}_{+}^{H-\frac{1}{\alpha }}\right]+b\left[{\left(t-s\right)}_{-}^{H-\frac{1}{\alpha }}-{\left(-s\right)}_{-}^{H-\frac{1}{\alpha }}\right]\right\}{M}_{\alpha }\left(ds\right),$$

(6)

Allowing us to derive the following equation can be derived:

① when $H=\frac{1}{\alpha }$, the value of fLsm is as follows:

$$L_{H,\alpha }\left(t\right)={\int }_{-\infty }^{\infty }\left(a-b\right)M_{\alpha }\left(ds\right),$$

(7)

if $a \mathrm{a}\mathrm{n}\mathrm{d} b$ are constants, then the $L_{H,\alpha }\left(t\right)$ series is mutually independent without any relativity.

② when $H<\frac{1}{\alpha }$, the fLsm integration kernel includes all the states of the past time, but the integration is inversely proportional to the past states; therefore, the predicted results are inversely proportional to the actual results. This is called short-range dependence (SRD).

③ when $H>\frac{1}{\alpha }$, the fLsm integration kernel contains all the states of the past time $s$, but the integration is directly proportional to the past state; this is called long-range dependence (LRD).

The LRD properties of the fLsm model are determined by $H$ and $\alpha$. The larger the product, the stronger the LRD properties. Special note $0<\alpha \le 1$: the fLsm model does not have LRD characteristics because $H$ is limited to the interval $\left(0,1\right)$, as shown in Figure 3.

3.2. Self-Similarity Property

A stochastic time series $X\left(t\right)$ is defined as having the following characteristics below:

$$X(t)\stackrel{d}{\Rightarrow }X(ct)\stackrel{d}{\Rightarrow }{c}^{H}X(t),$$

(8)

Furthermore, $X\left(t\right)$ is referred to as a self-similar process with exponent $H$.

The fLsm model is presented below:

$$L_{H,\alpha }\left(t\right)={\int }_{-\infty }^t({a(t-s)}^{H-\frac{1}{\alpha }}-b{s}^{H-\frac{1}{\alpha }})M\left(ds\right),$$

(9)

For any $n\ge 1$, the series $M\left({∆}_1\right),M\left({∆}_2\right),\cdots ,M\left({∆}_n\right)$ are mutually independent and ${∆}_1,{∆}_2,\cdots ,{∆}_n$ are mutually disjoint. Therefore, the series can be expressed using the self-similar parameter $H=\frac{1}{2}$:

$$M\left(cds\right)=c^{\frac{1}{2}}M\left(ds\right),$$

(10)

From (9) and (10), we obtain:

$$\begin{gathered} L_{H,\alpha }\left(ct\right)={\int }_{-\infty }^{ct}\left({a\left(ct-s\right)}^{H-\frac{1}{\alpha }}-b{s}^{H-\frac{1}{\alpha }}\right)M\left(ds\right) \\ =c^{H-\frac{1}{\alpha }+\frac{1}{2}}{\int }_{-\infty }^t({a(t-s)}^{H-\frac{1}{\alpha }}-b{s}^{H-\frac{1}{\alpha }})M\left(ds\right)=c^{H-\frac{1}{\alpha }+\frac{1}{2}}{L}_{H,\alpha }\left(t\right), \end{gathered}$$

(11)

Fundamentally, fLsm is a self-similar process characterized by the parameter $H-\frac{1}{\alpha }+\frac{1}{2}$.

4. Degradation Modeling of Fractional Lévy Stable Motion

4.1. Model Derivation

Lévy movement is the Brownian movement extend. In terms of the Langevin equation, the degradation model based on the fBm model is defined as follows:

$$dX\left(t\right)={\mu }_{B}dt+{\delta }_{B}d{B}_H\left(t\right),$$

(12)

where $\left\{X\left(t\right),t\ge 0\right\}$ is the degradation process and ${\mu }_B$ is the coefficient of the drift term, which describes the degradation speed rate; the term occupies a dominant position. ${\delta }_B$ is the coefficient of the diffusion term, which describes the stochastic and LRD characteristics. $\left\{B_H\left(t\right),t\ge 0\right\}$ is the standard fBm model. $d{B}_H\left(t\right)$ is an increment.

$d{L}_{H,\alpha }\left(t\right)$ replaces $d{B}_H\left(t\right)$ so that we obtain

$$dX\left(t\right)={\mu }_{B}dt+{\delta }_{B}d{L}_{H,\alpha }\left(t\right),$$

(13)

${\mu }_B$ is a time-varying function $\mu \left(t;\theta \right)$, which is called drift function; $\theta$ is the parameter vector in the drift function, and ${\int }_0^t\mu \left(t;\theta \right)dt$ can change as demonstrated below:

$${\int }_0^t\mu \left(t;\theta \right)dt=\mu \varphi \left(t\right),$$

(14)

$\varphi \left(t\right)$ usually has three defined functions: linear, nonlinear, and exponential power rate.

The iterative degradation model based on the fLsm is represented as follows:

$$X\left(t_k+∆t_k\right)=X\left(t_k\right)+{\int }_{t_k}^{t_k+∆t_k}\mu \left(s;\theta \right)ds+\delta \omega \left({∆t}_k\right),$$

(15)

$$\omega \left(∆t_k\right)=L_{H,\alpha }\left(t_k+∆t_k\right)-L_{H,\alpha }\left(t_k\right),$$

(16)

As a self-similar stochastic process, the following condition is satisfied for fLsm:

$$L_{H,\alpha }\left(at\right)-L_{H,\alpha }\left(0\right)\genfrac{}{}{0pt}{}{∆}{=}{a}^H\left[L_{H,\alpha }\left(t\right)-L_{H,\alpha }\left(0\right)\right],$$

(17)

and according to the Maruyama model:

$${\int }_0^{t}f\left(\tau \right){\left(d\tau \right)}^a=a{\int }_0^{\tau }{\left(t-\tau \right)}^{a-1}f\left(\tau \right)d\tau ,$$

(18)

$$dx=f\left(t\right){\left(dt\right)}^a,$$

(19)

Then, it can be concluded that:

$$L_{H,\alpha }\left(t_k+∆t_k\right)-L_{H,\alpha }\left(t_k\right)={\omega }_{\alpha }\left(t\right){\left(∆t_k\right)}^{H-\frac{1}{2}+\frac{1}{\alpha }},$$

(20)

where ${\omega }_{\alpha }\left(t\right)$ is the Lévy white noise conforming to $S_{\alpha }\left(0,1,0\right)$.

Therefore, the iterative form of the degradation model is transformed into:

$$X\left(t_k+∆t_k\right)-X\left(t_k\right)={\int }_{t_k}^{t_k+∆t_k}\mu \left(s;\theta \right)ds+\delta {\omega }_{\alpha }\left(t\right){\left(∆t_k\right)}^{H-\frac{1}{2}+\frac{1}{\alpha }},$$

(21)

Assume that $\alpha =1.75$, $\beta =0$, $\delta =0.195$, H = 0.75,$∆t_k=2,$ and the simulation time is 100. Figure 4 presents the simulated degradation path.

4.2. Parameters Estimation

The Hurst parameter H plays a key role in the modeling process of fLsm. The rescaled range method is the most widely used Hurst parameter estimation algorithm. Therefore, we propose to adopt the rescaled range method to calculate the Hurst parameter.

Other parameters of fLsm are estimated by the improved eigenfunction-based method. The steps are shown in Figure 5.

The proposed eigenfunction method first calculates diffusion parameter $\delta$. This procedure allows the value of the diffusion parameter $\delta$ to be estimated directly without the need to base it on the other three parameters. In addition, there are two auxiliary points to note in the estimation process of the new eigenfunction method: $\theta =1$ and ${\theta }_0$. One is the auxiliary point located at $\theta =1$. This is because the absolute value of the eigenfunction at $\theta =1$ is equal regardless of the value taken for the parameters of the eigenfunction:

$$\left|\phi \left(\theta ;\alpha ,\beta ,\mu ,\delta \right)\right|=e^{-\delta {\left|\theta \right|}^{\alpha }}{|}_{\theta =1}=e^{-\delta },$$

(22)

The second is the auxiliary parameter ${\theta }_0$, which is introduced to calculate the estimate of the stability index $\alpha$ and can be obtained with the help of the diffusion parameter $\delta$.

In order to estimate the diffusion parameter $\delta$ of the sample sequence $\left\{X\left(t\right),t=1,2,3,\cdots \right\}$, the eigenfunction of Equation (23) can be considered as a complex function associated with $\theta$.

$$\left|\phi \left(\theta ;\alpha ,\beta ,\mu ,\delta \right)\right|=\left|E\left\{e^{j\theta x}\right\}\right|=e^{-\delta {\left|\theta \right|}^{\alpha }},$$

(23)

Then, the logarithm is found for each side of Equation (23).

$$\mathit{ln}\left|\phi \left(\theta ;\alpha ,\beta ,\mu ,\delta \right)\right|=-\delta {\left|\theta \right|}^{\alpha },$$

(24)

Calculate the value of Equation (24) assuming $\theta =1$ to obtain an estimate of the diffusion parameter $\delta$.

$$\delta =-\mathit{ln}\left|\phi \left(1;\alpha ,\beta ,\mu ,\delta \right)\right|=-\mathit{ln}\left|E\left\{e^{jx}\right\}\right|,$$

(25)

The estimated value of the diffusion parameter $\delta$ is given by the following equation.

$$\widehat{\delta }=-\mathit{ln}\left|\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right|=-\mathit{ln}\frac{1}{N}\left|{\sum }_{i=1}^N{e}^{j{x}_i}\right|,$$

(26)

Before calculating the estimate of the stability index $\alpha$ and the second auxiliary parameter ${\theta }_0$ needs to be calculated in addition to the known auxiliary point $\theta =1$. The characteristic functions can degenerate to different models when the stability index $\alpha$ equals different values. Therefore, the new eigenfunction method uses the Cauchy and Gaussian eigenfunctions corresponding to $\alpha =1$ and $\alpha =2$ to calculate the auxiliary parameter ${\theta }_0$. This is carried out by calculating the maximum distance between the absolute values of the Cauchy eigenfunction $e^{-\delta \theta }$ and the Gaussian eigenfunction $e^{-\delta {\theta }^2}$.

$$\frac{d}{d\theta }\left|e^{-\delta {\theta }^2}-e^{-\delta \theta }\right|=0,$$

(27)

$$2\theta e^{-\delta {\theta }^2}=e^{-\delta \theta },$$

(28)

Taking the logarithm of Equation (28) yields the following Equation (29).

$$\frac{\mathit{ln}2\theta }{{\theta }^2-\theta }{|}_{{\theta }_0}=\widehat{\delta },$$

(29)

The value of the auxiliary parameter ${\theta }_0$ is obtained by solving Equation (29). Substituting the obtained ${\theta }_0$ into Equations (24) and (25):

$${\theta }_0^{\alpha }=\frac{\mathit{ln}\left|E\left\{e^{j{\theta }_{0}x}\right\}\right|}{\mathit{ln}\left|E\left\{e^{jx}\right\}\right|}=\frac{\mathit{ln}\left|\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)\right|}{\mathit{ln}\left|\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right|},$$

(30)

Therefore, the estimated value of the stability index $\alpha$ can be expressed as:

$$\widehat{\alpha }=\frac{\frac{\mathit{ln}\left|\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)\right|}{\mathit{ln}\left|\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right|}}{\mathit{ln}{\theta }_0}\ast {\mathit{log}}_{{\theta }_0}\left(\frac{\mathit{ln}\left|\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)\right|}{\mathit{ln}\left|\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right|}\right),$$

(31)

where $\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)=\frac{1}{N}\left|\sum _{i=1}^N{e}^{j{\theta }_0{x}_i}\right|$.

The estimates of skewness index $\beta$ and drift coefficient $\mu$ are calculated from the logarithmic form of the characteristic function.

$$\mathit{ln}\left|\phi \left(\theta ;\alpha ,\beta ,\mu ,\delta \right)\right|=-\delta {\left|\theta \right|}^{\alpha }+j\left[\delta {\left|\theta \right|}^{\alpha }\beta \frac{\theta }{\left|\theta \right|}tan\left(\frac{\pi \alpha }{2}\right)+\mu \theta \right],$$

(32)

Suppose $\theta =1$, then Equation (32) can be expressed as:

$$\ln \phi \left(1;\alpha ,\beta ,\mu ,\delta \right)=-\delta +j\left[\delta \beta tan\left(\frac{\pi \alpha }{2}\right)+\mu \right],$$

(33)

Equation (33) reveals that the estimated value of the diffusion parameter can be expressed in terms of the real part. This is expressed as follows:

$$\widehat{\delta }=-Re\left\{\ln \widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right\},$$

(34)

Furthermore, considering $\ln y=\ln \left|y\right|+jarg\left(y+2k\pi \right)$, we can obtain:

$$Re\left\{\ln \widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right\}=\ln |\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)|,$$

(35)

$$Im\left\{\ln \widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)\right\}=\arg \widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right),$$

(36)

Considering the two auxiliary parameters $\theta =1$ and ${\theta }_0$, one can solve Equations (35) and (36) to obtain the skewness index $\beta$ and the drift coefficient $\mu$.

$$\widehat{\beta }=\frac{Im\left\{\ln |\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)|-{\theta }_0\ln |\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)|\right\}}{\widehat{\delta }tan\left(\frac{\pi \widehat{\alpha }}{2}\right)\left({{\theta }_0}^{\alpha }-{\theta }_0\right)},$$

(37)

$$\widehat{\mu }=\frac{Im\left\{{{\theta }_0}^{\widehat{\alpha }}\ln |\widehat{\phi }\left(1;\alpha ,\beta ,\mu ,\delta \right)|-\ln |\widehat{\phi }\left({\theta }_0;\alpha ,\beta ,\mu ,\delta \right)|\right\}}{{{\theta }_0}^{\alpha }-{\theta }_0},$$

(38)

5. Case Study

The datasets consist of multiple multi-dimensional time series data. Each data set is divided into training and test samples. A separate time series comes from different engines of the same type and represents changes in various characteristics: pressure, temperature, equipment rotation speeds, etc. Systems start up with different degrees of initial wear and production changes. Each engine has 21 sensors that collect various measurements related to the engine’s condition during operation. The engine runs normally at the beginning of each time series, and a malfunction occurs at some point during the series. The error grows in magnitude in the training set until the system crashes. The proposed method was applied to the subset FD001 of the C-MAPSS dataset [29] for aero-engine unit prognostics. A schematic representation of the engine is depicted in Figure 6.

5.1. Feature Selection of Health Indicators

The subset FD001 varies in operating and fault conditions, and the data are contaminated with sensor noise. Therefore, we applied filtering to the training data. The filtering was carried out using the smooth filter function in MATLAB with a setting window of 10. Since there are considerable differences between the magnitudes of the data, we applied normalization to the data after filtering. In order to enhance degradation features, we applied PCA and fusion in the construction of HI (Figure 7 and Figure 8).

5.2. Prediction of Remaining Useful Life

The definition of RUL is given as

$$T=inf\left\{t_k:X\left(t_k\right)\ge \omega |X\left(t_0\right)<\omega \right\},$$

(39)

where $\omega$ indicates the degradation threshold. The non-Gaussian nature of fLsm makes it impossible to obtain the analytical formulas for the RUL prediction; therefore, we refer to the Monte Carlo method [32,33]. The point RUL prediction is the mode of a large number of RUL simulations [34,35]. The empirical threshold value $\omega$ was set to be 4.25. The fLsm model with three drift terms was chosen. The corresponding parameters are listed in

Table 1. Input parameter of the fLsm model for different degradation trends.

Model	H	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\mu }$	$\mathit{\delta }$
$M_1$	0.8669	1.7019	0.0686	1.3723	1.9623
$M_2$	0.8669	1.9473	0.0700	1.2531	2.0652
$M_3$	0.8669	1.9564	0.0715	1.3275	1.9514

. The corresponding results of RUL are shown in Figure 9.

We use four standard error evaluation metrics to compare the predictive accuracy: mean absolute error, root mean square error, mean absolute percentage error, and Health Degree [36].

$$MAE=\frac{1}{m}{\sum }_{i=1}^m\left|r_i-r_i^{\mathrm{*}}\right|,$$

(40)

$$RMSE=\sqrt{\frac{1}{m}{\sum }_{i=1}^m{\left(r_i-r_i^{\mathrm{*}}\right)}^2},$$

(41)

$$MAPE=\frac{1}{m}{\sum }_{i=1}^m\left|\frac{r_i-r_i^{\mathrm{*}}}{r_i^{\mathrm{*}}}\right|\times 100\mathrm{\%},$$

(42)

$$HD=1-\frac{{\sum }_{i=1}^m{\left(r_i-r_i^{\mathrm{*}}\right)}^2}{{\sum }_{i=1}^m{\left(r_i-\overline{r}\right)}^2},$$

(43)

where $r_i$ denotes the predicted value of the RUL, which is obtained for the i-th prediction starting point. $r_i^{*}$ is the actual value of the RUL at the i-th prediction starting point, m is the number of predictions, and $\overline{r}$ denotes the mean value of $r_i$. The results of error estimation are given in

Table 2. Prediction accuracy of RUL.

Model	MAE	RMSE	MAPE	HD
$M_1$	1.5000	1.6432	0.0147	0.9869
$M_2$	1.3000	1.4491	0.0153	0.9898
$M_3$	1.2000	1.3416	0.0054	0.9913

. As we can see, the exponential drift term of model M3 is superior to the power drift term of model M2 and the linear drift term M1.

The fLsm model with an exponential drift is compared with the fGC model [37,38], the fBm model [22,23], and the LSTM model [39]. The results are shown in Figure 10. The LSTM neural network uses a single-input–single-output model with 200 hidden units in LSTM layers, a stacking layer of 2, and an initial learning rate of 0.005. The results of error estimation are shown in

Table 3. Prediction accuracy of RUL.

Model	MAE	RMSE	MAPE	HD
fLsm	1.2000	1.3416	0.0054	0.9913
fGC	1.4000	1.5632	0.0151	0.9893
fBm	1.6000	1.8439	0.0373	0.9815
LSTM	1.9000	1.9748	0.0476	0.9811

. Thus, the proposed approach outperformed other predictive models.

6. Conclusions

This paper proposes an iterative degradation model based on fLsm with LRD and heavy-tailed characteristics. We have also demonstrated that fLsm has LRD property and a heavy probabilistic tail. The RUL prediction is based on the Monte Carlo simulation of the degradation model, which is the Langevin-type stochastic differential equation driven by fLsm. A new eigenfunction method and rescaled range method were applied for the parameter identification of the model.

The model’s training data, i.e., the HI, underwent PCA analysis, identifying degradation features and enabling us to construct HI using the dimension reduction PCA algorithm. The proposed model with an exponential drift is compared with the fGC, fBm, and the LSTM model. Real turbofan engine data were used for the case study, and the results indicate the advantage of the fLsm prediction model over other methods. We may utilize a more advanced drift term in the RUL prediction model based on fLsm.

In the future study, we will expand application fields, and the prediction accuracy and prediction length in time will be further improved by parameter optimization [40,41].