Stochastic Identification and Analysis of Long-Term Degradation Through Health Index Data

1. Introduction

Diagnostics and prognostics using long-term monitoring data have increased with the development of condition monitoring systems. Efficient use of long-term condition monitoring data, collected over months or even years, is a crucial element for diagnosis and prognosis. In recent years, many methods have been developed and published in the area, which can be categorised into three main groups: Data-driven approach (including machine learning (ML) based approach [1,2,3,4], and statistic based model approach [5,6,7,8] ), physics-based approach [9,10], and hybrid approach [11,12]. Physics model-based approaches try to explain the degradation process by taking advantage of the physics of the degradation process based on damage and fracture mechanics [13]. This family of approaches could provide accurate results. However, the exact physics knowledge of complex systems is not always available or too expensive to extract [14]. Therefore, employing the physics-based approaches has significant restrictions to use in real applications. The data-driven approach attempts to build a model of the degradation process using historical data. The data-driven approach is usually more effective than the physic-based approach when we are faced with complex engineering systems such as wind turbines, aircraft, and mining machines.This class of approach is divided into two subclasses: machine learning-based model and statistical-based model. The machine learning-based approach is a powerful tool for modelling, segmenting, and predicting complex systems where the degradation processes are hard to interconnect with physics or statistical-based models. However, these methods require a large amount of historical degradation process data for training, which are not accessible for most practical equipment. Compared to previous methods, statistical-based approaches do not require much degradation data and mechanical knowledge of the equipment. Furthermore, the statistical-based model has enormous potential to consider degradation uncertainties. The hybrid approach wants to use the advantages of physical-based and data-driven approach through their integration. More detailed information on the hybrid approach can be found in literature [11,12].

Selecting the appropriate statistical model (stochastic or random variable), which has closely matched the degradation process in an actual application, is crucial to address diagnostic and prognostic problems. Therefore, different statistical model-based approaches have been developed in the literature to consider the problem. Some of these stochastic model-based approaches are developed based on Gaussian processes, such as the Brownian motion (BM), also known as the Wiener process [15]. Wang et al. [16] developed a degradation model that employs linear BM that calculated RUL assuming that its distribution is inverse Gaussian. Bian et al. [17] introduced a covariance-dependent degradation model. However, for models developed on the basis of BM, it is difficult to describe the long-range dependence. To solve this problem, Xi et al. [18] used a degradation model based on fractional Brownian motion (FBM), which provided an approximate explicit solution to estimate the remaining useful life (RUL).

Another part of the stochastic model has been developed based on the non-Gaussian process to model a degradation process such as the Gamma process [19], the generalised Cauchy process (GC) [20], and fractional Levy stable motion (FLSM) [21].

However, the mentioned stochastic models are practical and suitable only under particular consideration. For example, both BM and FBM are developed on the basis of stationary increments. Also, the increment of both mentioned methods has followed the same distribution. However, BM and FBM can describe non-stationary degradation processes when only the drift term has a nonlinear trend. The degradation model with this non-stationary characteristic is recalled only in the increment expectations, as the drift terms can express the deterministic part of the degradation process [22]. Generally, the heterogeneity of actual degradation processes is represented in deterministic and random parts. Therefore, it is necessary to consider non-stationary characteristics in a random part. Some methods describe the non-stationary part with a specific function, such as power law. However, this assumption is too simple for complex systems [23].

The autoregressive model (AR) is one of the influential and well-known approaches used to model and analyse time series [5,24,25]. Nevertheless, classical AR modelling is suitable only for stationary time series, and in numerous cases assumptions about stationarity are too restrictive or inappropriate [22]. Nonetheless, Żuławiński et al. l [23] introduced a new framework for the long-term HI model, which in part of this framework uses a robust AR approach to identify the characteristics of random parts. They used a robust estimator to identify the scale (variance in the Gaussian case) and normalised a random component. Then they fitted the AR model to the stationary part of the signal. This framework can also model both the deterministic and random parts of the degradation process. However, it is difficult to use this model for online applications.

In the case of non-stationary time series, time-varying parameters or adaptive modelling seem to be potential options. TVC-AR models have been developed since the early 1980s, particularly in the Bayesian framework; see [26] for an excellent review. TVC-AR models are powerful tools for expressing non-stationary time series with quasi-periodic latent features [27]. These types of time series appear in several applications that involve, for example, structure health monitoring [28,29], biomedical [30,31,32], speech signal processing [33], and financial time series [34,35].

To address the above-mentioned problems, we introduced a framework for identifying and modelling complex long-term condition monitoring data with non-stationary characteristics based on TVC-AR. We describe how to identify all components (including deterministic and random parts) and how finally to build a model for prognosis and remaining useful life when we have a signal with non-stationary characteristics in both deterministic and random components. Also, this model contains all dependent characteristics such as location and variance (scale). The contributions in this paper are summarised as follows.

• A model of HI data was proposed as a three-segment sequence with non-stationary characteristics in both random and deterministic components to describe the degradation process, which can be used to simulate the artificial data set.
• Online identification of the time-varying random characteristics component like mean (location) and variance (scale), and also the dependency between them, is described for the TVC-AR model.
• A long-term data model based on TVC-AR is proposed for identification and modelling, and extensive experiments are carried out on the simulated data set and FEMTO and wind turbine datasets to verify its effectiveness.

The paper is formed as follows: After the introduction, in Section 2, the critical parts of the proposed approach are defined in theory. Then, in Section 3, the suggested model is simulated, and the results are presented. The results of using the proposed approach to two benchmark data sets are shown in Section 4 with an indication of all intermediate steps. Ultimately, in Section 5, the discussion made for the results of the previous section and the conclusions are formed in Section 6.

2. Methodology and Theory

2.1. Degradation Model

In the PHM community, there are several stochastic models that are used to describe the degradation process. These models are usually composed of the deterministic part (which is trying to qualify the global trend of the degradation process ) and the random part (which is trying to consider the uncertainty of the degradation process). These models are usually selected based on the physic of the degradation process with different deterministic and random trends [23]. In this work, due to the complex trend of growing crack, we selected a time-varying model with 3 different regimes. The first regime has a constant deterministic trend which refers to the healthy state of the system (the machine is working stable ). The second regime has a linear growing trend which assigns to the degradation of the machine (when the length of the crack is growing up gradually). The last regime has an exponential trend which is attributed to the critical state of the degradation machine (when the length of the crack is growing up dramatically). Also, as we discussed modeling the random part has very important in reflecting the uncertainty. Therefore, we proposed a model with non-stationary characteristics. The scale’s trend (variance) of the random part is changed during the degradation process based on the regime we are in. It may also include the dependency on the random component during the process. In Table 1, we show the main characteristics of the data corresponding to different stages for more details on the model.

Also in Figure 1, it can be seen as an example of the degradation curve that generated based on the Table 1.

2.2. Methodology

The flowchart of this methodology procedure is presented in Figure 2. First, we express the general approach for the proposed framework. Then, we present the algorithms used in this framework in detail.

2.3. Theory

2.3.1. Deterministic Component

First, we identify the deterministic component in the long-term data. The deterministic component indicates a global trend of the degradation process, which plays a crucial role in describing, analysing, and predicting the degradation process. Most failures in machinery failure are related to the initiation and growth of cracks until the failure. Therefore, according to the nature of the growing crack in material, our preliminary research on the real data sets indicates that a single deterministic function cannot adequately describe this component [23]. More precisely, we observe a much more complex situation when the type of the deterministic component (called here trend) changes depending on the regime corresponding to the good, warning, and alarm regimes. Therefore, in this paper, we compute the empirical location measure for overlapping segments from windows of a given length w to identify the deterministic component without considering long-term data segmentation problems. This step is essential for identifying the deterministic component that is further removed from the raw data. After removing the deterministic component, the series is denoted as $\left\{y\right(n\left)\right\}$

2.3.2. Separating the Random and Deterministic Component

In the general framework, the first crucial step is identifying the deterministic trend. The classical statistic assumed as the location measure is just a sample mean. In this case, we compute the moving average (MA) for overlapping segments (with number of overlapping samples o) of length w and assume it as the deterministic component of the signal.

2.3.3. Random Component

In this subsection, we tried to identify and model random components. Identifying, analysing, and modeling random components are essential because they help us consider the degradation process’s uncertainty, which may be related to the time-varying load, changing environmental conditions, etc.

• Time varying autoregressive model (TVC-AR)
• Model

The TVC-AR model of order p is an extension of the autoregressive (AR) model stated by assuming time varying (instead of constant) coefficients $\{{\phi }_l\left(n\right);l=1,...,p\}$. The evolution of the model is then described with the equation [36]:

$$y_n=\sum _{l=1}^p{\phi }_l\left(n\right)y_{n-l}+{\epsilon }_n.$$

(1)

As usual we assume that ${\epsilon }_n$ is a Gaussian white noise series with zero mean and unknown variance ${\sigma }^2$. ${\epsilon }_n$ is an innovation to the model, thus it is independent from $y_{l-1}$ for $l>0$. Bayesian framework is used to model the time-varying coefficients of the AR model, by assuming $\{{\phi }_l\left(n\right);l=1,...,p\}$ are also random variables.

Stochastic constraints (also called smoothness priors) are put on the AR coefficients, which is expressed by stochastic difference equations:

$${\Delta }^q{\phi }_l\left(n\right)={\nu }_{nl}$$

(2)

where ${\Delta }^q$ is the qth order difference operator defined by $\Delta {\phi }_l\left(n\right)={\phi }_l\left(n\right)-{\phi }_l(n-1)$, ${\Delta }^q{\phi }_l\left(n\right)={\Delta }^{q-1}\left(\Delta {\phi }_l\left(n\right)\right)$. The innovation term ${\nu }_{nl}$ is a Gaussian white noise sequence (indexed with n) with zero mean and unknown variance ${\tau }^2$. It is also supposed that ${\nu }_{nl}$ and ${\epsilon }_n$ are independent of each other. The behavior of the model depends on the value of the difference order q, this hyperparameter is usually choosed considering lower order cases such as q equal to 1 or 2.

• State space representations

Each TVC-AR model can be expressed by a state space representation:

$$\begin{array}{c}{x}_n=F{x}_{n-1}+G{w}_n,\hfill \\ y_n=H_n{x}_n+{\epsilon }_n.\hfill \end{array}$$

(3)

The matrix components of the model depend on specific order q. We present here the case of $q=1$ and $q=2$. If $q=1$, the state space model is parameterised as follows:

$$x_n={({\phi }_1\left(n\right),...,{\phi }_p\left(n\right))}^T,w_n={({\nu }_{n1},...,{\nu }_{np})}^T,F=I_p,G=I_p,H_n=(y_{n-1},...,y_{n-p}),$$

(4)

where ${(·)}^T$ denotes the matrix transposition, $I_p$ the $p\times p$ identity matrix . For $q=2$, the state space model is given by:

$$\begin{array}{c}{x}_n={({\phi }_1\left(n\right),...,{\phi }_p\left(n\right),{\phi }_1(n-1),...,{\phi }_p(n-1))}^T,\hfill \\ w_n={({\nu }_{n1},...,{\nu }_{np})}^T,\hfill \\ F=\left[\begin{array}{cc}2I_p& -I_p\\ I_p& 0\end{array}\right],\hfill \\ G=\left[\begin{array}{c}{I}_p\\ 0\end{array}\right],\hfill \\ H_n=(y_{n-1},...,y_{n-p},0,...,0).\hfill \end{array}$$

(5)

The covariance matrices of ${\omega }_n$ and ${\epsilon }_n$ are respectively given by $Q=E\left\{w_n{w_n}^T\right\}={\tau }^2{I}_p,R=E\left\{{\epsilon }_n^2\right\}={\sigma }^2$. Given observations $y_1,...,y_N$ and the initial conditions $x_{0|0}$ and $V_{0|0}$, we use Kalman filter algorithm to obtain the conditional mean and covariance matrix of the state vector at each time $n=1,...,N$. This Bayesian procedure consists of repeated two steps of prediction and filtering:

• [Prediction]

  $$\begin{array}{c}{x}_{n|n-1}=F{x}_{n-1|n-1},\hfill \\ V_{n|n-1}=F{V}_{n-1|n-1}{F}^T+{\tau }^{2}G{G}^T.\hfill \end{array}$$

  (6)
• [Update]

  $$\begin{array}{c}{k}_n=V_{n|n-1}{H}_n^T{(H_n{V}_{n|n-1}{H}_n^T+{\sigma }^2)}^{-1},\hfill \\ x_{n|n}=x_{n|n-1}+k_n(y_n-H_n{x}_{n|n-1}),\hfill \\ V_{n|n}=(I_m-k_n{H}_n)V_{n|n-1}.\hfill \end{array}$$

  (7)

  Above, $x_{i|j}$ and $V_{i|j}$ denote respectively the mean and covariance matrix of the state vector (estimated from the Kalman filter) $x_i$ given the data $y_1,...,y_j$. The term $k_n$ is called Kalman gain.

After performing Kalman filter procedure, by the following backward iteration we can obtain smoothed estimation $x_{n|N}$ of the state vector and corresponding covariance matrix $V_{n|N}$:

• [Smoothing]

  $$\begin{array}{c}{A}_n=V_{n|n}{F}^T{V}_{n+1|n}^{-1},\hfill \\ x_{n|N}=x_{n|n}+A_n(x_{n+1|N}-x_{n+1|n}),\hfill \\ V_{n|N}=V_{n|n}+A_n(V_{n+1|N}-V_{n+1|n})A_n^T.\hfill \end{array}$$

  (8)

  Note that the state vector $x_n$ consists of the coefficients ${\phi }_l\left(n\right)$, and thus their estimates ${\widehat{a}}_l\left(n\right)$ are directly obtained via $x_{n|N}$.
• Estimation and identification of the model

The conditional density function of $y_n$, conditioned on the data $y_1,...,y_{n-1}$, is given by the following formula:

$$f\left(y_n\right|y_1,...,y_{n-1};{\sigma }^2,{\tau }^2)={\left(2\pi {\nu }^2\left(n\right)\right)}^{-{\displaystyle \frac{1}{2}}}exp\left\{-{\displaystyle \frac{{(y_n-H_n{x}_{n|n-1})}^2}{2{\nu }^2\left(n\right)}}\right\},$$

(9)

where ${\nu }^2\left(n\right)=H_n{V}_{n|n-1}{H}_n^T+{\sigma }^2$ denotes the variance of the difference between measurement $y_n$ and prediction $H_n{x}_{n|n-1}$, based on observations $y_1,...,y_{n-1}$. Since, given p and q (the order of the model and the differentiation order), the joint density function of the random vector $y_1,...,y_N$ is following:

$$f(y_1,...,y_N|{\sigma }^2,{\tau }^2)=\prod _{n=1}^{N}f\left(y_n\right|y_1,...,y_{n-1};{\sigma }^2,{\tau }^2),$$

(10)

the log-likelihood for the hyperparameters ${\sigma }^2$ and ${\tau }^2$ can be approximated with:

$$\ell ({\sigma }^2,{\tau }^2;p,q)=-{\displaystyle \frac{1}{2}}\left\{Nlog2\pi +\sum _{n=1}^N(log{\nu }^2\left(n\right)+{\displaystyle \frac{{(y_n-H_n{x}_{n|n-1})}^2}{{\nu }^2\left(n\right)}})\right\}.$$

(11)

In order to estimate the hyperparameters ${\sigma }^2$ and ${\tau }^2$, the method of maximization of the log-likelihood is applied. For more in-depth details about TVC-AR and estimation of its hyperparameters, the reader can look into [36].

3. Simulation

3.1. Generating the Degradation Data

In the following section based on Section 2.1, the health index is generated. As we discussed in Section 2.1 the health index (HI) is constructed from two time-varying parts (deterministic and random part) as follows:

$$\begin{array}{c}\hfill HI\left(t\right)=D\left(t\right)+R\left(t\right),\end{array}$$

(12)

where $D\left(t\right)$ and $R\left(t\right)$ are represented as deterministic and random parts of the degradation process respectively. Both of these parts are following the assumption of three regimes based on Table 1. The behavior of the deterministic part in Equation (12) has the following form:

$$\begin{array}{c}\hfill D\left(t\right)=\left\{\begin{array}{cc}{c}_1\hfill & 0<t\le {\tau }_1,\hfill \\ a_{1}t+c_2\hfill & {\tau }_1<t\le {\tau }_2,\hfill \\ a_{2}exp\left(b_{1}t\right)+c_3\hfill & {\tau }_2<t\le N,\hfill \end{array}\right.\end{array}$$

(13)

where $\tau$ is denote changing point between regimes. $c_1$ is a constant value that is used to represent a healthy regime. Also, the constants $a_1$, $a_2$, and $b_1$ are used to construct the linear and exponential trend which corresponds to the degradation and critical regimes. Furthermore, the parameter of $c_2$ and $c_3$ are tuned in such a way that keeps the continuity of the degradation curve ( health index) in changing points.

In the simulation part, we consider Gaussian distributions for the random component of the degradation process. The time-varying random component corresponding to $R\left(t\right)$ is generated in the following way:

$$\begin{array}{c}\hfill R\left(t\right)=SC\left(t\right)\widehat{R}\left(t\right),\end{array}$$

(14)

where $SC\left(t\right)$ corresponds to the scale (variance ) of the random part which is constructed in the following way:

$$\begin{array}{c}\hfill SC\left(t\right)=\left\{\begin{array}{cc}{a}_{3}t+b_2\hfill & 0<t\le {\tau }_1,\hfill \\ a_{4}t+b_3\hfill & {\tau }_1<t\le {\tau }_2,\hfill \\ a_{5}exp\left(b_{4}t\right)\hfill & {\tau }_2<t\le N,\hfill \end{array}\right.\end{array}$$

(15)

where values $a_3,b_2,a_4,b_3,a_5,b_4$ are constant and tuned in this way which $SC\left(1\right)={\sigma }_1$, $SC\left({\tau }_1\right)={\sigma }_2$, $SC\left({\tau }_2\right)={\sigma }_3$ and $SC\left(N\right)={\sigma }_4$.

And $\widehat{R}\left(t\right)$ is AR1 with time-varying coefficient which is constructed as follow:

$$\begin{array}{c}\hfill \widehat{R}(t+1)=a_{6}t\widehat{R}\left(t\right)+ϵ\left(t\right),\end{array}$$

(16)

where $ϵ\left(t\right)$ is independent identically distributed random variables (iid) and it arises from Gaussian distribution $ϵ\left(t\right)\sim \mathcal{N}(0,1)$. For simplicity, we consider that the distribution in each regime is the same and time-varying AR coefficient is increased linearly; however, as was mentioned in Section 2.1, in practice, it may be different for different regimes.

In panel (a) in Figure 3 we present the deterministic component $D\left(t\right)$ and the scale function $SC\left(t\right)$ for the following values of the parameters: ${\tau }_1=1500$, ${\tau }_2=2500$, $N=3000$, ${\sigma }_1=0.25$, ${\sigma }_2=1$, ${\sigma }_3=3$, ${\sigma }_4=12$ and $c_1=0$.

3.2. Results of Proposed Approach

This subsection applies the proposed methodology to data generated by the suggested model, see Figure 3. The proposed method’s implementation results are presented in Figure 4. Panel (a) shows the simulated health index’s identified trend (deterministic component). As can be seen in panel (a), the proposed method could identify the deterministic component properly. Also, the panel (b) of Figure 4 illustrated the random component of the generated health index after removing the identified deterministic part from the health index. Panel (c) of Figure 4 demonstrated the identified time-varying scale (variance) of the random component. By comparing of the real scale (variance) and identified scale proved the efficacy of the proposed approach to the identified scale. In the end, the identified time-varying AR coefficient is shown in Figure 4, as can be seen in panel (d); the proposed method detected the AR coefficient as proper.

The proposed procedure was repeated for 100 simulations from the same model to examine the performance of the proposed methodology, and the results are presented in Figure 5. The left column is illustrated the results for trend, variance, and AR(1) coefficient for 100 simulated data. The real value is presented in blue, and the identified values are shown in cyan. It can be seen that the proposed method could detect and follow the actual values with acceptable performance. Also, in the right column, the RMSE is calculated for the left column plot and is shown with the boxplot.

4. Real Data Analysis

This section will apply and evaluate the proposed methodology for available real datasets. These data are typically employed as benchmark datasets for different papers and competitions and have specific behavior corresponding to noise properties and deterministic trends. In the following, essential information about objects, experiments, and data will be recalled, and appropriate references are provided.

4.1. FEMTO Dataset

The FEMTO data set is collected from the PRONOSTIA test rig; see Figure 6 by the Franche-Comté Electronics Mechanics Thermal Science and Optics–Sciences and Technologies Institute (FEMTO). The data were collected by using two acceleration sensors and one temperature sensor during the test. The sampling frequency for collecting acceleration data is 25600 Hz, which is recorded data for 0.1 seconds every 10 seconds [37]. This data set has been used in many research in recent years for health index construction [38,39,40,41,42,43,44,45], health index evaluation, or segmentation of health index [46,47,48,49,50,51], and the prediction of RUL [52,53,54,55,56,57].

4.2. Wind Turbine Data Set

This data set is collected data from the sensor mounted in the high-speed bearing shaft of the wind turbine, see Figure 8. In this data set the bearing inner race energy is calculated every 10 min for +50 days, see Figure 9. More information about the procedure of constructing a health index could be found in the following reference [58]. Finally, the inner race-bearing fault has occurred, which was proved by inspection, see Figure 8. Also, this data set (we call a wind turbine data set) has been used to predict RUL by several papers in recent years; see, e.g., [58,59,60,61].

4.3. Result for FEMTO Data Set

The results for the FEMTO data are illustrated in Figure 10. To detect the trend component (deterministic) (panel (a)), we selected the length of the window 51. As can be seen, the detected trend changes over time; i.e., for the first stage, it is nearly constant; then, it transforms its nature and acts as a linear function, while in the last stage it behaves as an exponential function. Also, this identified trend confirmed the correctness of the assumption of the three stages (constant, linear, and exponential) degradation model, which is discussed in Section 2.1. After separating the deterministic trend from the health index, we obtain the random component (panel (b)). As can be seen, it conforms to a non-homogeneous sequence with a non-constant scale (time-varying scale). Panel (c) presents the identified scale of the random component, and we selected $[0.995,0.9]$ as the discount factors for AR and scale (variance) discount factors, respectively. Also, these results confirm our assumption about increasing the random component over time, and one can detect that the scale of a random part grows in a non-linear manner. It can be seen in the last part of the curves that the amplitude of the scaled noise increases dramatically, which can have a significant effect on the quality of the prognosis. The first four time-varying AR coefficients are presented in panel (d). It should be noted that here we selected four coefficients for AR as hyperparameters. As can be seen in panel (d), these coefficients change over time with a fluctuating trend. However, the values of these mentioned coefficients are not significant. The residual of the TVC-AR model is illustrated in panel (e), while its empirical ACF is presented in panel (f). The plot of the empirical autocorrelation function shows that the data can be considered independent observations. Eventually, to confirm our first assumption about Gaussian noise, we plotted the empirical and theoretical tails (panel (g)) plus the well-known non-Gaussian distribution, including the stable and student t distribution. For more information on these mentioned non-Gaussian distributions, please see Appendix A. We conclude that the residual series corresponds to the alpha distribution with $\alpha =1.902$. Additionally, the probability density functions (PDF) of the theoretical distributions of the residual signal are presented in panel (h). This result also rejected our assumption about Gaussian noise. However, the level of impulsivity is relatively low.

We analyse quantile lines constructed based on the identified parameters and proposed a fitted model for FEMTO data sets to confirm our results. We should note that the presentation of quantile lines constructed on the basis of the fitted model is one of the most typical techniques for validating that an offered model is properly fitted to the data. This approach is frequently employed in both the literature and practical applications. The procedure for building the quantile lines is as follows: first, the model is fitted to the real data, then we simulate the number of trajectories by using the fitted model, and for every time point we estimate the quantiles (at appropriate levels). They are named quantile lines. If the real data fall into the constructed intervals (with reasonable probability), then we can confirm that the fitted model is proper. The fitted models have been synthesised and the simulation results are illustrated in Figure 11. In this figure, the real datasets are shown by purple lines, and by blue lines, the constructed quantile lines on the levels of $5\%$ and $95\%$ obtained based on 400 simulated trajectories corresponding to the fitted models. Furthermore, based on the results presented in Figure 11, it can be concluded that the fitted models maintain the specific nature of the data sets and can be considered the optimal ones, for example, for the prognosis.

4.4. Result for Wind Turbine

The results for the wind turbine data are illustrated in Figure 12. To detect the trend (deterministic) component (panel (a)), we selected the window of length 51. As can be seen, the detected trend (deterministic component) changes over time as the FEMTO data set. After separating the deterministic trend from the health index, we obtain the random component (panel (b)). As can be seen, it contains a non-homogeneous sequence with a non-constant scale (time-varying scale). Panel (c) presents the identified scale of the random component, and we selected $[0.995,0.99]$ as the AR and scale (variance) discount factors, respectively. Also, these results confirm our assumption about increasing the random component over time, and one can detect that the scale of a random part grows in a non-linear manner. The scale increases differently for each regime ($0-2000$ and $2000-5000$), after which the situation becomes more complicated. The first ten time-varying AR coefficients are presented in panel (d). It should be noted that here we selected ten coefficients for AR as hyperparameters. As seen in panel (d), these coefficients converge to a constant value over time. The values of these mentioned coefficients, particularly the first coefficient of the AR model, are valuable. The residual of the TVC-AR model is illustrated in panel (e), while its empirical ACF is presented in panel (f). The plot of the empirical autocorrelation function shows that the data can be considered independent observations. Eventually, we plotted the empirical and theoretical tails (panel (g)) plus the well-known non-Gaussian distribution, including the stable and student t distribution. We conclude that the residual series corresponds to a non-Gaussian distribution. Additionally, the probability density functions (PDF) of the theoretical distributions of the residual signal are presented in panel (h). This result also rejected our assumption about Gaussian noise. However, the level of impulsivity is relatively low.

In addition, we analyse quantile lines constructed based on the fitted model for the wind turbine data set to confirm our results as the previous one. The fitted models have been synthesised and the simulation results are as illustrated in Figure 13. In this figure, the real datasets are shown by purple and blue lines, the constructed quantile lines on the levels of $5\%$ and $95\%$ obtained based on 400 simulated trajectories corresponding to the fitted models. Also, based on the results presented in Figure 13, it can be concluded that fitted models maintain the specific nature of the data sets and can be considered the optimal ones, for instance, for the prognosis.

5. Discussion

The results presented for both real data sets confirm the efficiency of the proposed approach. The deterministic trend, scale (variance), and AR coefficients are detected as time-varying functions for both cases. Also, based on the FEMTO results, the deterministic part of the FEMTO data set follows the assumption of a constant, linear and exponential function trend for the degradation process. On the other hand, for the wind turbine data set, the exponential trend can not be clearly seen that may arise from this fact; this data set is not entirely run to failure data [58]. In addition, in the deterministic trends of the wind turbine data set, a few fluctuations can be seen that may originate from different reasons and phenomena, such as self-healing. Likewise, the variation of the scale (variance) during the time for FEMTO and wind turbine data sets ultimately confirms our first assumption about the time-varying scale (variance), and it can be seen that the scale (variance) FEMTO data set follows the assumption of constant, linear and exponential function, while for the scale of the wind turbine data set, we can see some complex trend, including fluctuations after t = 2000, which can cause problems in the efficiency of the model to be used to applications of prediction by over and underestimate. We should consider the fact that we assumed that the observed noise is Gaussian, while it may not be a correct assumption for such data with non-Gaussian characteristics. Also, the Kalman filter used in our proposed model is not robust against intense non-Gaussian noise, and it may be better to use a robust version. We will try to deal with this in our future work. Correspondingly, the time-varying AR coefficients for both cases; however, the AR coefficients for wind turbines are not significantly giant and, after some time, have approximately constant trends that are predictable based on the fact that this data set is gathered under constant speed and condition. Also, we additionally presented the residual of the TVC-AR model and fit different distributions (Gaussian and non-Gaussian) to investigate whether the noise with Gaussian characteristics assumption for these data sets is valid. As seen in the results section, the results for FEMTO and wind turbine data sets are more compatible with non-Gaussian distributions like stable noise and student t noise, while it is not so far from Gaussian noise.

Finally, we also compared the results to the method introduced by Żuławiński et al. [23]; this method considers non-Gaussian characteristics, non-homogeneous manners, and autodependence in the time series. However, this method considers constant AR coefficients. for more information on this model, please find this reference [23]. For this work, we compared quantile lines constructed based on fitted models according to the identified variables for FEMTO and the wind turbine data set. As can be seen in panels (a) and (b) in Figure 14, both models could keep the specific nature of the datasets as well, however, as can be seen in panel (b), the 5 percent quantile detected by the proposed methods is a little bit more acquired rather than Żuławiński et al. model. Also, for the bold difference between the results of these two models, we presented 95 percent of the quantile line constructed by these two models plus the difference; see panels (c) and (d) in Figure 14. So, we can figure out that the quantiles for the first and the second regime (healthy stage and degradation regime ) take the same values, while for the critical stage, we can see the Żuławiński et al. model. take a higher value than our proposed model. And we should consider this; however, our proposed method is not robust as Żuławiński et al. model, and it can be affected by strong non-Gaussian characteristics, due to this model developed by the Bayesian theorem can tolerate soft non-Gaussian characteristics. Also, our proposed method can identify and model the time-varying AR coefficients, which is a crucial point, particularly for applications that work under non-stationary conditions. In contrast, Żuławiński et al. model considers constant AR coefficients. Furthermore, the proposed method can be used for online applications, which is a critical feature for prognosis and diagnosis.

6. Conclusions

This paper discusses modelling and identifying long-term data with non-stationary characteristics from condition monitoring systems. The primary purpose of this paper was to introduce a method for modelled such data. In this research, we proposed another way to model of long-term data to generalise the classical model.This model has the potential to be used in online applications. Moreover, we considered that the characteristics of the random part might be time-varying, such as the growing scale ( variance). The last novelty was modelling the random components by time-varying autoregressive time series, which can give us the ability to describe the time-varying dependency of a random component. Next, we suggest reliable algorithms to identify the components mentioned above. Ultimately, we presented an approach to identify each component and synthesise all of them to simulate the exemplary trajectory of the suggested model. The proposed procedure has been applied to the simulated and two real data sets named FEMTO and the wind turbine data set. The simulated data set was designed based on the three-stage model see Figure 1 with changing scale (variance) and coloured noise to be closer and more realistic to the real data set. The results of the proposed approach to the simulated data set confirmed the method’s efficacy in identifying and modelling the time-varying deterministic part and the random part’s time-varying scale (variance) and AR coefficients.

Also, the proposed procedure has been applied to two real data (FEMTO and wind turbine datasets). The results of these two data sets prove our first assumptions about non-homogeneous manners and time-varying autodependence in the time series. Furthermore, we have seen non-Gaussian characteristics in residues of the TVC-AR model, which can influence the results of this proposed approach. For these two real cases, the non-Gaussian characteristics were not so strong; the proposed approach could cover the effects of it, but certainly, if the level of the non-Gaussianty is increased, the efficiency of the method will be decreased. So in the future, we will try to deal with this problem.

Appendix A.1. Stable Distribution

The stable distribution (also known as Lévy alpha-stable distribution) is defined by its characteristic function and is characterised by four parameters: $\alpha$ (stability), $\beta$ (skewness), $\sigma$ (scale) and $\mu$ (location). However, by considering the symmetric case with a standardized scale, we can assume $\beta =\mu =0$, $c=1$, and the corresponding characteristic function is reduced to the following equation:

$$\begin{array}{c}\hfill E\left[e^{itX}\right]=e^{-{\left|t\right|}^{\alpha }}.\end{array}$$

(A1)

The parameter $\alpha$ is called the stability index and takes the value from $(0,2]$ interval. It should be noted that the stable distribution reduces to the Gaussian distribution when $\alpha =2$. In case $\alpha$ decreases, the distribution becomes significantly non-Gaussian and heavy-tailed, [62].

Appendix A.2. Student T Distribution

The Student’s t distribution is known to be one of a family of curves of one parameter. This distribution is usually employed to test a hypothesis regarding the population means when the population standard deviation is unknown. The probability density function (PDF) of the Student’s t distribution is the following:

$$\begin{array}{c}y=F\left(x\right|\nu )={\displaystyle \frac{\Gamma \left(\frac{\nu +1}{\nu }\right)}{\frac{\nu }{2}}}{\displaystyle \frac{1}{\sqrt{\nu \pi }}}{\displaystyle \frac{1}{{(1+\frac{x^2}{\nu })}^{\frac{\nu +1}{2}}}}\hfill \end{array}$$

(A2)

where $\nu$ is the degree of freedom and $\gamma (.)$ is the Gamma function. The result y is the probability of observing a particular value of x from the Student’s t distribution with $\nu$ degrees of freedom.

It should be nothe that: The mean of the Student’s t distribution is $\mu =0$ (mean) for degrees of freedom $\nu$ greater than 1. If $\nu$ equals 1, then the mean is undefined.

The scale (variance) of the Student’s t distribution is ${\displaystyle \frac{\nu }{\nu -2}}$ for degrees of freedom $\nu$ greater than 2. If $\nu$ is less than or equal to 2, then the scale (variance) is undefined.