Particle Filter Based Approach for Remaining Useful Life Prediction of High-Speed Shaft Bearing in Wind Turbine Generators

Abstract

This work involves a novel data-driven procedure using vibration analysis for bearing health prognosis. In this work, we investigate the time-domain features and applying spectral kurtosis features in order to extract the damage indicators which eventually represent the degradation of the high-speed shaft bearing (HSSB). These damages were characterized by their Monotonicity, Trendability, and Prognosability. The most appropriate indicator was then used as a health index for the remaining useful life (RUL) prediction task. In this study, we used an integrated approach based on Particle Filter approach which was then developed for direct RUL prediction of HSSB. This methodology was validated using real world vibration data wind turbine gearbox. The experimental results and the prognostics metrics like fitness degree equal to 0.9941 shown that the Particle Filter approach is more feasible prediction tool.

1. Introduction

The economic development was based on low cost energy. The growth of renewable energy, especially that is produced by wind turbine generator (WTG) helps continued growth while curbing CO₂ emission. While the costs of WTG power production is low, operations and maintenance costs are higher than expected due to higher rates of electrical and mechanical failures. According to national renewable energy laboratory (NAREL), failure in HSSB accounts for 48% of all drive-train damage. Mechanical failure in (WTGs) leads to sudden downtime and electricity production cessation causes a higher than expected maintenance cost. Prognostics and health management (PHM) of WTGs aim to predict the future behavior of the generator’s health condition by estimating the RUL of HSSB. This allows requires a proactive maintenance policy that can reduce the operations and maintenance cost and provides better balance of plant.

Mechanical failures in WTG bearings drive train is often the result of moisture contamination of the oil, especially in offshore application. Contamination causes reduction in oil lubrication, reducing the life cycle [1]. Bearing failures can create long down time and costly maintenance. To maintain balance of plant, implementation of a Condition Based Maintenance (CBM) program should be considered. CBM along with PHM are able to identify the failure and estimate the RUL of the elements.

PHM is a maintenance paradigm that merges diagnostics, future loads and a damage propagation model to estimate the RUL. Diagnosis describes the current state and the damaged component in the system. Applying the current state to a propagation model, based on future estimated load, and a point where maintenance is appropriate (e.g. the threshold), the RUL [2, 3, 4, 5] can be estimated. As shown in Figure 1 the PHM cycle is composed of three principal parts ranged as follows:

• Observe contains two steps: Acquire data from sensors installed on the machine, then data processing/feature extraction to generate a condition indicator (CIs). Some other techniques, such as digital signal processing improve the signal to noise and quantify the representative statistics of component damage.

• Analysis is a three steps process. First, a condition assessment is made of the observed CIs. Second, the given diagnostic represents is given representing the component state of health. This also includes fault isolation and identification. Third, using the current state of health and the threshold, estimate the RUL and the associated confidence level of the RUL for the given component.

• Act is a two-step process. First, decision support tool evaluates the evidence, which can generate a maintenance intervention. The last step is Human Machine interfaces, which display the PHM status.

Figure 1.

PHM cycle retyped and adapted from the ISO 13374.

A PHM architecture [6, 7] leads to implement a paradigm that supports maintenance planning with the goal of eliminating unscheduled maintenance and improving operational readiness.

The activity associated with PHM is growing, with a number of scientific papers and several reviews have been published in CBM domain. Lee et al. [4] introduce different strategies for the detection of failures associated with rotating equipment.

Lee goes on to describe many of the analysis algorithm that can be used, and the analysis performance (advantages and disadvantages).

The review by Jardine et al. in [8] is an excellent overall summary of CBM. Even though this publication is old, all the proposed approaches used in data processing and maintenance decision support are still valid until now.

In [9], the authors applied a Kalman smoother approach with confidence bounds in order to predict the RUL of HSSB. The prognostic approach is hybrid using both a data driven approach and physics-based degradation model. Vibrating data collected over 50 days was used to estimate the damage associated with propagating inner race crack. Using a Kalman filter, the unknown parameters associated with Paris’ Law model were determined in order to estimate RUL.

In [10], Kramti et al. applied Elman neural network (ENN) technique. They proposed an approach using time-domain features and frequency-domain via Spectral Kurtosis (SK) as inputs. These features were extracted from raw vibration bearing signal in order to predict the RUL of HSSB. The architecture of an ENN is built with two hidden layers. The first hidden layer is composed of five neurons, and the second hidden layer is composed of three neurons. The hidden layers transfer function is Logarithmic sigmoid. The output layer is composed of one neuron using a pure linear transfer function ranging between 0 and 1. This ENN gave a prediction horizon of 20 days. While a powerful model, there were some gaps in the predicted RUL when compared to the true RUL with large fluctuations.

The authors in [11] used a support vector regression (SVR) approach based on classical, time domain features. This SVR model used spectral kurtosis (SK) to predict the RUL of HSSB. Spectral kurtosis derived indices reduce the noise by using the Short-Time Fourier Transform (STFT) and provide good trendability and monotonicity metrics.

In the SVR study, the model was trained by 60% and tested by 40% of area under SK index data. In other words, it was trained using 40% of the data and tested using 60% of the same type of data used in the first step. The experimental results have shown that the estimated RUL based on area under SK index tracks the actual RUL with a small prediction error using 60% of training data. Unfortunately, the model did not predict any SK values.

The [12], the authors used Acoustic Emission (AE) data which were delivered from four bearings. The fault feature was Root Means Square (RMS) and Signal Intensity Estimator (SIE). All bearings features were fitted. The bearing feature bearing 1, were used as inputs of training process of three learning machines using: Gaussian Process Regression (GPR), support vector machine regression (SVMR) and multilayer artificial neural network (ANN). The other three bearings features are used as input of process test. The experimental results have shown that ANN gave the lowest error compared to SVMR and GPR.

The proposed method aims to improve the reliability of the HSSB and reduce maintenance cost. Therefore, the new proposed failure prognostics method is based on the analysis of the behavior of vibration signal. The proposed failure prognostics method uses the most suitable feature, which is then mapped to a Health index (HI). The RUL is then the estimated time for the current HI, to the HI threshold were it is appropriate to do maintenance.

The remainder of this work is organized as follows. Section 2 reports the bearing characteristic and describes the data used in the experiment. In Section 3 the proposed methods and techniques used like features extraction methodology, particle filter prognosis is introduced with the model degradation. In Section 4 we detail the procedure to obtain the RUL. Section 5 A discussion and a comparison with some previous work and methodologies in the literature. Finally, we conclude this chapter and future work are synthetized in Section 6.

2. Experimental steps

The vibration signals data were obtained from an online condition monitoring system from Green Power Monitoring System (GPMS). During 50 days of measurements, data was collected at a rate of once per 10 minutes (144 acquisitions per day), where once per day raw vibrating signal data was downloaded. The sensor monitoring the high-speed bearing was sampled at 97656/second for six seconds, along with tachometer data. The data was collected from a real wind turbine generator (S88, Suzlon) with 2 MW electric power generation.

The failing bearing supports the high-speed shaft, which drives the generator. This shaft rate was approximately 1800 revolutions per minute (these are doubled feed induction machines, so that the gearbox/generator is not synchronous at 60 Hz, while the output of the generator is electronically controlled to be line synchronous).

The vibrating signal data given by the Green Power Monitoring System from USA. The vibrating sensor was made by MEMS accelerometer (analog device ADXL001, the bandwidth is a 32 kHz, the resonance at 22 kHz, with a sensitivity of +/−70 Gs. The data was sampled by a delta-sigma analog to digital converter ADS1271 (24 bit). The vibrating sensor was installed above the HSSB.

On the fiftieth day, the HSSB was inspected, where the damage to the inner race fault was verified, as indicated in Figure 2. The HSSB model is 32222-J2 [13] tapered roller bearing it is made by SKF. Bearing dimensions are 200 mm in outside diameter; the inside diameter is 110 mm, the width is 56 mm, it has a 20 rolling elements each one has 46 mm of width, the taper angle is 16°, weighs 7.10 kg and the speed-limit is 3200 r/min.

Figure 2.

Inner race fault.

The rotating speed of WTG depends on climatic condition, the pressure altitude, offshore or onshore location. Indeed, these environmental conditions directly influence the bearings radial and axial loads. When wind speed is near the cut-in wind (i.e. the lowest wind speed needed for power generation), the main shaft will be in lower speed (high speed shaft perhaps 1500 rpm). When the wind speed is higher than cut-in wind, the main shaft speed operates higher on the power curve, close to 1800 RPM, the Figure 3 shows the mean of speed shaft over 50 days. Loads are moderated by the wind turbine blade pitch. Figure 4 shows the concatenation of the 50 days of measurements data in a 6-second period each time. It is clear that the RMS increases with time and damage propagation.

Figure 3.

The mean speed variation in high speed shaft over 50 days.

Figure 4.

The historical vibration data ending with inner race fault.

1. Proposed prognostic approach

   In this chapter, we propose an effective prognostic method to estimate degradation in the bearing. The proposed method for HSSB failure is summarized in Figure 5. The proposed method consists mainly of three task features definition, feature selection and RUL prediction. During the first task of the proposed method, two types of features are used to extract information from vibration data: classical time-domain features and a SK feature (which is a frequency-domain feature). Once the features definition task has been done, feature selection was performed. In this task, we use three metrics in order to determine a suitable feature, which will inform the RUL prediction step. Finally, the RUL prediction step implies the use of the particle filter method to predict the RUL according to prognostics metrics.

2. Features

   1. The Classical features

      Classical features are presented in Table 1 [5, 14] were applied on historical vibration signal run-to-failure over 50 days as shown in Figure 6(a).

      Generally, time-domain features are well proven and historically, are the basis of many condition-monitoring systems. The classic features used in this study are; RMS, Kurtosis, Skewness, Peak to peak, Crest Factor, Mean, Standard deviation (Std), Energy, and Entropy. A detailed description of the classical features are given as follows in Table 1.

   2. Features derived from SK

Figure 5.

Flow chart of the proposed prognosis method.

Figure 6.

The trend for the 50 days (a) representing the time domain, classical features, and (b) features derived from SK.

Feature name	Mathematical expression
RMS	${\left(\frac{1}{N}{\sum }_{i=1}^N{x}_i^2\right)}^{\frac{1}{2}}$
Kurtosis	$\frac{1}{N}{\sum }_{i=1}^N\frac{{\left(x_i-\overline{x}\right)}^4}{{\rho }^4}$
Skewness	$\frac{1}{N}{\sum }_{i=1}^N\frac{{\left(x_i-\overline{x}\right)}^3}{{\rho }^3}$
Peak to peak	$x_{\max }-x_{\min }$
Crest Factor	$\frac{x_{\max }}{\mathit{RMS}}$
Mean	$\overline{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$
Standard deviation	$\sigma ={\left(\frac{1}{N}{\sum }_{i=1}^N\left(x_i-\mathit{\text{mean}}\right)\right)}^{\frac{1}{2}}$
Energy	${\sum }_{i=1}^N{x}_i^2$
Entropy	$-{\sum }_{i=1}^N{x}_{i}log\left(x_i\right)$

Table 1.

The proposed features.

In order to estimate the RUL process, it is necessary to identify features that show trendability, monotonicity and that can be used as a surrogate for component damage. For this study, we used features derived from SK. SK typically involves the band-pass filtering of the raw data to remove signal that is not associated with the fault. Also, it improves the SNR.

The approach requires a band-pass filter to find a bandwidth that emphasizes the demodulated impulsive signature (associated with bearing fault) which is hidden in the raw vibrating signal. For a full treatment on SK, pleases refer to Randall and Antoni [15].

The kurtosis is defined as:

$$K=\frac{\frac{1}{N}{\sum }_{i=1}^N{\left(x_i-\overline{x}\right)}^4}{{\left(\frac{1}{N}{\sum }_{i=1}^N{\left(x_i-\overline{x}\right)}^2\right)}^2}$$E1

Where i is the sample index, x and x̄ are the sample time index and sample mean respectively and N is the number of samples. This normalized fourth moment is defines the “peakedness” of the signal. The spectral kurtosis of signal is described as the kurtosis of its spectral elements. The SK is defined as follow [16].

$$\mathit{SK}\left(f\right)=\frac{⟨\left|X^4\left(t,f\right)\right|⟩}{{⟨X^2\left(t,f\right)⟩}^2}-2$$E2

Where X²(t,f) and X⁴(t,f) are the second-order and fourth-order cumulate respectively of a band-pass filtered signal of x(t) around f. ‹•› correspond to the time frequency averaging operator.

The most important characteristics of this description defined as:

• In case of stationary system, the SK is a permanent function of frequency.

• In case of stationary Gaussian system, the SK is the Gaussian process.

The SK has been shown to be a more powerful indicator of damage than raw signal kurtosis. SK detects the high-frequency train of impulse derived from a damaged bearing. It can be shown that the SK of a non-stationary system x(n) and a damaged stationary noise b(t) source is defined as

$${\mathit{SK}}_{\left(x+b\right)}\left(f\right)=\frac{{\mathit{SK}}_x\left(f\right)}{{\left(1+\rho \left(f\right)\right)}^2}+\frac{\rho {\left(f\right)}^2{\mathit{SK}}_b}{{\left(1+\rho \left(f\right)\right)}^2}$$E3

Where f ≠ 0, then ρ(f) is the signal-to-noise ratio (SNR) as function of frequency.

If b (t) is an additive stationary Gaussian noise independent of x(t), then the SK is transformed into.

$${\mathit{SK}}_{\left(x+b\right)}\left(f\right)=\frac{{\mathit{SK}}_x\left(f\right)}{{\left(1+\rho \left(f\right)\right)}^2}$$E4

Now, it can be seen that the SK is able to characterize and detect bearing damage that is masked by a non-stationary signal in the frequency domain. Additionally, this shows that the when using the SK, for a nominal, stationary signal, the value of SK is approximately 0. For non-Gaussian, damaged signal (e.g. transients), the value is ≠ 0, see Figure 7 for more explication.

Figure 7.

Sate of health in HSSB using spectral kurtosis.

1. 3. HSSB degradation model

The principal cause of damage and early bearing failures are overload, inadequate lubrication, lubrication contamination, or corrosion. Bearing faults can appear on inner race, outer race, rollers and cage. Bearing damage may result in large, quasi-periodic impacts, which degrade exponentially over time [17, 18].

When the model degradation with a defined level of damage, the measured data can be used to calibrate and identify the parameters of model. If the model parameters are known, they can be applied to estimate the prospective behavior of the damage.

• The model degradation is defined as

$$d=aexp{\left(\mathit{bt}\right)}^2$$E5

Where d is the magnitude degradation, a, b are the model constant parameters and t is the time index. This model allows a “best fit” trend of the health index of HSSB.

1. C. Particle filter based prognosis

In this section we present a brief review of Particle Filter (PF). More detailed study on PF can be found in [19] In this chapter we present only the basic theory.

PF has been applied in numerous engineering domain such as robotics, aerospace, automatic control, etc. and more currently in diagnosis and prognosis. PF may also be known as the sequential Monte Carlo method. The main uses of PF are to accurately model the degradation state with a set of particles. The PF has a corresponding state values, and a correlated set of particles weights, which correspond to the discrete Probability masses of the distinctive particles.

In PF, the Bayesian update is processed in sequential mode with samples (or particles) having the information probability of hidden parameters: when a new measured data is obtained, the posterior step is used as the information for the present step, and the parameters are updated by multiplying it with the likelihood function.

The particle can be created and updated recursively by the use of non-linear state-transition model, illustrating the evolution of the system under control. The Bayesian tracking task is described by two equations as follows

• The state equation, the state transition function f.

$$x_k=f\left(x_{k-1},{\theta }_k,v_k\right)$$E6

• The model observation, measurement function h.

$$z_k=h\left(x_k,{\omega }_k\right)$$E7

Where k is the time step index, x_k is the system state at the preceded step in this work x_k represents the damage state of HSSB. θ_k is a model parameters vector, z_k is measured data, ω_k and v_k are respectively measurement data and process noise. All these variables fluctuate at each time step, and the progress from the k-1 to the k step is produced by the transition function f. As mentioned before, this work is about prognostics area so the state transition function f is designed as damage model.

According to the state model in (Eq. (6)), the HSSB degradation model in (Eq. (5)) can be edit in the following shape

$$x_k=x_{k-1}exp{\left(b_k\mathrm{\Delta }t\right)}^2$$E8

Where the process noise v_k is neglected because it can be managed through the uncertainty in model parameters. In case of measurement function, it is supposed that z_k is the selected feature that reflects the HI of HSSB (see Figure 8) This feature includes measurement Gaussian noise ω_k ∼ N (0,σ) which is applied with unknown standard deviation σ. Consequently the unknown parameters are θ = [b,σ], containing the damage state x_k which is acquired based on the model parameter b_k.

Figure 8.

Health index data.

A probability density function (PDF) p(x_k|z_1:k) is needed to obtain the distribution of the possible states of x at time k. The initial state is estimated by the state distribution p(x₀|z₀) = p(x₀). It is assumed that the PDF is known. The optimal Bayesian solution is given by iterating prediction and update functions respectively:

• Prediction

$$P⟨x_k|z_{1:k-1}⟩=\int P⟨x_k\left|x_{k-1}\right.⟩P⟨x_{k-1}\left|z_{1:k-1}\right.⟩{\mathit{dx}}_{k-1}$$E9

• Update

$$P⟨x_k|z_{1:k}⟩=\frac{P⟨z_k\left|x_k\right.⟩P⟨x_k\left|z_{1:k-1}\right.⟩}{P⟨z_k\left|z_{1:k-1}\right.⟩}$$E10

Unfortunately, this solution cannot be obtained systematically, but PF is a robust approach designed to obtain an approximate solution via feedback.

Figure 9 shows an algorithm of the PF which can be recapitulated as follow [20, 21].

The primary step consists in subdividing the initial state distribution p(x₀) into n samples called also particles. The next three steps are then reiterated until the appropriate results are achieved see Figure 9.

1. Prediction: in this step the particles are generated through the state model from k-1 to the k step providing each incremental time a new PDF. The information from preceding step should be fully accessible.

2. Update: the model parameters and state degradation are updated. The measurement data are utilized to computing the likelihood function p(z₀|x_k), and provides weights to the particles. The new particles states converge to the actual one provided by the last measured data which have the higher likelihood, as they are more suitable to describe the system state.

3. Resampling: the concept of resampling is to keep all samples (i.e., particles) have the same weight. Particle that have a low weight are removed. And particles with the high weights are manifold. This operation can favorite the best filter results. Namely a higher number of particles with poor weight can degenerate the filter results in the previous stages.

Figure 9.

Illustration of PF principle estimation process and prediction process.

All these steps are used during the learning process. In case there is no available data, no measured data z_k should be used to compute the likelihood function, and the prognosis pass into prediction phase. In this prediction step, the particles are diffused via the state model. When the failure threshold is crossed the latest distribution of particles are the most appropriate state, (i.e., the state expressed via the weighted average of the particles).

3. Tracking fault degradation

1. A. Features definition

In order to define the degradation of HSSB, two types of features have been used [5]. The first is time domain indices of classical features such as: kurtosis, skewness, mean, standard deviation (std), peak to peak, root means square (RMS), energy, entropy and crest factor. The second are features based on the SK but operated on by the same classical features: kurtosis-SK, skewness-SK, mean-SK, standard deviation-SK (std-SK), peak to peak-SK, root means square-SK (RMS-SK), energy-SK, entropy-SK, crest factor-SK and area under curve-SK which is added. These features are extracted for the 50 day data set, as shown in Figure 6.

1. B. Feature selection

Feature selection performed to identify the best analysis for classification and prognosis. The principal idea is find the most relevant features that provide useful information and that can predict a future state. This was done by two methods:

• By transforming the same feature set into new domain, such as main component evaluation or linear discriminant evaluation [4]. This approach unfortunately generates new feature sets in the new domain that are dissimilar to actual features.

• By choosing a feature that allows for prognosis. As noted in recent literature on prognosis, metrics such as: trendability, prognosability and monotonicity [22, 23] can be used to determine the most appropriate feature set.

Monotonicity can define the main negative or positive trend of the feature. This is a powerful metric to detect degradation because degradation in bearing is an irreversible and a growing process. The monotonicity of a group features is affected by the mean between the number of positive and negative step for every assessment point of time. Suppose that n is the number of assessment point of time, the monotonicity will be defined as follows

$$\mathit{\text{Mono}}=\left|\frac{\mathit{no}\mathit{\text{of}}\left(\frac{d}{\mathit{dt}}<0\right)-\mathit{no}\mathit{\text{of}}\left(\frac{d}{\mathit{dt}}>0\right)}{n-1}\right|$$E11

The range of monotonicity value between 0 and 1, non-monotonic features will take the value of 0 and greatly monotonic features will take the value of 1.

Trendability quantifies the correlation of the features vs. time. If the feature is constant, the correlation with time will be 0. However, if the derivative of the feature is linear, the correlation with time will take on a non 0 value. In the same way, correlation can change with increase in non-linearity (i.e, a nonlinear feature will result in low correlation). Trendability is ranged between −1 and 1 it is defined as follows

$$\begin{array}{l}\mathit{\text{Tren}}=\left|\mathit{\text{corrcoef}}\left(\mathit{\text{time}},\mathit{\text{feature}}\right)\right|\\ =\frac{n\left(\sum \mathit{\text{timefeature}}\right)-\left(\sum \mathit{\text{time}}\right)\left(\sum \mathit{\text{feature}}\right)}{\sqrt{\left[n\sum {\mathit{\text{time}}}^2-{\left(\sum \mathit{\text{time}}\right)}^2\right]-\left[n\sum {\mathit{\text{feature}}}^2-{\left(\sum \mathit{\text{feature}}\right)}^2\right]}}\end{array}$$E12

Prognosability is the exponential of the standard deviation of degradation feature measure divided by the difference between final and first value of degradation feature measure.

$$\mathit{\text{Prog}}=exp\left(-\frac{\mathit{std}\left(deg\mathit{\text{measures}}\right)}{\mathit{\text{mean}}\left(\left|\mathit{\text{final}}deg\mathit{\text{measure}}-\mathit{\text{first}}deg\mathit{\text{measure}}\right|\right)}\right)$$E13

According the numerical values of suitability (Eq. (14)) in Tables 2 and 3, the most significant feature for the prognostic task among all ones, is the mean-SK which has an exponential growth. Note the feature is usually combined with noise. The noise can obscure the trend and reduce the power of the RUL estimation.

Features	Monotonicity	Trendability	Prognosability	Suitability
Kurtosis	0.1020	0.8081	0.7810	1.6911
Std	0.1020	0.6555	0.6944	1.4519
Peak to peak	0.1020	0.8004	0.7599	1.6623
RMS	0.1020	0.6479	0.6892	1.4391
Skewness	0.0612	0.1007	0.7913	0.9532
Energy	0.1020	0.6319	0.7084	1.4423
Crest factor	0.0204	0.7852	0.7419	1.5475
Mean	0.0204	0.2496	0.5232	0.7932
Entropy	0.0204	0.5275	0.5393	1.0872

Table 2.

Classical features results during 50 days.

SK-features	Monotonicity	Trendability	Prognosability	Suitability
Kurtosis-SK	0.1020	0.8873	0.6469	1.6362
Std-SK	0.0612	0.7702	0.7957	1.6271
Peak to peak-SK	0.0204	0.7997	0.7988	1.6189
RMS-SK	0.0204	0.7928	0.7938	1.6070
Skewness-SK	0.0204	0.9024	0.6881	1.6109
Energy-SK	0.0612	0.7061	0.8316	1.5989
Crest factor-SK	0.0612	0.8896	0.6525	1.6033
Mean-SK	0.1429	0.8852	0.7730	1.8011
Entropy-SK	0.0612	0.8332	0.5920	1.4864
Area-under curve	0.0612	0.8851	0.7730	1.7193

Table 3.

Derived features-SK results during 50 days.

$$\mathit{\text{Suitability}}=\mathit{\text{Monotonicity}}+\mathit{\text{Prognosability}}+\mathit{\text{Trendability}}$$E14

Extracted features are usually associated with noise. The noise with opposite trend can sometimes be harmful to the RUL prediction. In addition, one of the feature performance metrics, introduced above is not robust to noise. Therefore, a causal moving mean filter with a lag window of 5 steps is applied to the most suitable extracted features, where “causal” means no future value is used in the moving mean filtering. The Figure 10 showing the mean-SK feature before and after smoothing, the smoothed mean-SK is used as health indicator of HSSB.

The smoothed mean-SK feature as shown in Figure 8 is normalized in range [0 1]. 0 corresponds to 0% degradation and 1 corresponds to 100% degradation. This feature is used as Health index (HI) data in the following prognosis steps and the smoothed feature is normalized using (Eq. (15)).

Figure 10.

The selected feature.

$$\mathit{\text{mean}}-\mathit{SK}\text{'}=\frac{\mathit{\text{mean}}-\mathit{SK}-\mathit{\text{mean}}-{\mathit{SK}}_{\min }}{\mathit{\text{mean}}-{\mathit{SK}}_{\max }-\mathit{\text{mean}}-{\mathit{SK}}_{\min }}$$E15

4. RUL prognostic

1. A. RUL prognosis based particle filter approach

Measured Data is applied to predict model parameters, which are used then to estimate the RUL. The damage model equation in (Eq.(8)) is defined as follow: the time interval ∆t equal to 1 which corresponds to one inspection per day. Also, the model parameters b_k and the damage state at the previous step x_k-1, and the standard deviation of measurement error s. The initial distribution is P× Q matrix of probability parameters of the initial distribution, which P is the number of unknown parameters and Q is the probability parameters. In this case, the prior information are not available, it is supposed that the initial distribution of these unknown parameters (P = 3) are uniform, where the probability parameters are (Q = 2), lower and upper bounds:

$$x_0\sim U\text{(0.01,0.03),} b_0\sim U\text{(0.038, 0.04),} s_0\sim U\text{(0.251, 0.253).}$$

The other setting for the prognosis using particle filter are the number of samples (or particles) n and significance level for adjusting the prediction interval (PI) and the confidence interval (CI). In this work, it is used (n = 5000 particles) and 90% of significance level. For more details according to the number of samples please read the work published in [19].

The other results can be plotted such as the model parameter and the prediction of failure state can be obtained by the use of sampling results during the updating process. The sampling results can be displayed for any variable at each step. In this work the exact values of the model parameters are known, the result should be compared with the known values. The exact value of b = 0.03919 and s = 0.2520, the failure state can be calculated using the (Eq. (8)). More graphical results are shown in the following plots: Figures 11–13.

Figure 11.

Estimation of b model parameter.

Figure 12.

The s value estimation.

Figure 13.

Particle filter RUL prediction.

Once the model parameters are classified as a physisc based approach, the mathematical exponontial fucntion in (Eq. (8)) is trained using a the derived data. The particle filter uses these as imputs to predict the remaining time until the degradation propagates to a maintenance threshold.

The expremental results were made using a total of 50 measurement data points for HSSB HI (see Figure 8). The Particle filter process is runing only with 48% of measured data (24 points) and the rest of degradation trend were predicted by the proposed method.

The main goal of the particle filter based prognosis is to predict the degradation behavior by the use of the exponential degradation model. If the model accuratly represents the HI, then it can be appliyed to find the RUL using (Eq. (16)). The true RUL is achieved by substracting the current time (24 days) from the failure threshold time (50 days). The threshold time is given by the last value of HI which is 1. According the Figure 13, and (Eq. (16)) the true RUL is 26 days.

$${\mathit{\text{True}}}_{\mathit{RUL}}=t_{\mathit{\text{FT}}}-t_{\mathit{CT}}$$E16

Figure 14 shows, 5 percentile equal to 16, median equal to 24 and 95 percentile equal to 28 which are caused by 90% of significance level interval. The median value result can be compared with the true RUL in Figure 13 which is computed using the (Eq. (16)). Therefore, the median value of RUL prediction of 24 days is fairly accurate compared with the true RUL of 26 days. The RUL prediction can be more exact by decreasing the time interval after the current time.

1. B. Particle filter perfermance

In this section, the discussion is going to compare the robustness and performance of the proposed particle filter approach. Defining prognostics metrics allows comparison and evaluation of different RUL algorithm. By comparing the true and predicted RUL, we can define statistical metrics to measure performance. This comparison needs to be evaluated by the following metrics.

• Prognostic horizon (PH): The prognostic horizon [24] is described as the difference between the first time when the predicted RUL continuously resides in the accuracy zone and the last cycle time (i.e. End Of Life (EOL)). According the Figure 13 the accuracy zone has an invariable bound with a value of ±5% error with respect to the last cycle, plotted as two parallel broken lines. The first time when predicted RUL resides in the accuracy zone is at the cycle 24 and the last cycle is 50, thereby PH is 26. The prognostics approaches with a larger PH designate a greater performance, which gives earlier RUL prediction with more reliability.

• Convergence: Eventually, the convergence [25] can be defined with a non-negative error metric (Eq. (17)) of precision or prediction accuracy between the true RUL and predicted RUL. We consider the actual or linear RUL values (denoted X), the predicted RUL (denoted X̂), and the length of predicted data is n. The RMSE is defined as follows: In this chapter, the relative error is illustrated below:

Figure 14.

RUl destribution with percentiles.

$$E=\left|X{}_i-\widehat{X_i}\right|$$E17

• RUL error: presents the value computed between true RUL and predicted RUL for each day. The lower value confirms that this method has a good way to predict bearing degraded mode, mathematically the error percent defined as [12, 26].

$$\mathit{\text{Error}}\left(\%\right)=\frac{X{}_i-\widehat{X_i}}{X{}_i}\times 100$$E18

• Root Mean Square Error (RMSE) is used to measure the precision of prediction because it is able to make the error and predicted value at the same magnitude [27] We consider the actual or linear RUL values (denoted X), the predicted RUL (denoted X̂), and the length of predicted data is n. The RMSE is defined as follows

$$\mathit{\text{RMSE}}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(X{}_i-{\stackrel{⌢}{X}}_i\right)}^2}$$E19

• Mean Absolute Percentage Error (MAPE) is generally applied to determine the error size as a percentage. However, it is not advisable with small data sets. The MAPE is defined as follows

$$\mathit{\text{MAPE}}=\frac{100}{n}\sum _{i=n}^n\left|\frac{X_i-{\stackrel{⌢}{X}}_i}{X_i}\right|$$E20

• Fitness degree (sometimes called the R² coefficient), gives an indication of good prognostics when the R² value is close to 1. R² is defined as fellow

$$R^2=1-\left(\frac{{\left(X_i-\widehat{X_i}\right)}^2}{X_i^2}\right)$$E21

• Relative Error Analysis (REA) is used to measure the precision which is defined in percentage. It presents relative information between the measurement and the size of data measured. The error is proportional to the size of the RUL being measured. REA is defined as follows

$$\mathit{REA}=\frac{100}{n}\sum _{i=1}^n\frac{X_i-\widehat{X_i}}{X_i}$$E22

• Accuracy metric (A) calculates the “exactitude” between the true RUL and the predicted RUL. The result of this metric is presented as a percentage, if the accuracy is close to 100% that prove the predicted RUL is similar to the true RUL. The accuracy is defined as

$$A\left(\%\right)=100\times \left[1-\frac{\left|X_i-\widehat{X_i}\right|}{X_i}\right]$$E23

These above mentioned metrics are used to evaluate PF results as shown in Table 4.

Metrics	Particle Filter
Error %	7.69
RMSE	0.2828
MAPE	0.1538
R^2	0.9941
REA	0.1538
A %	92.30
PH	26

Table 4.

Perfermance of the proposed methods.

The PF method adopted in this chapter is called hybrid prognosis approach. This study is based on the combination of data driven approach and exponential model degradation. This combination makes a very powerful prognostics tool. This idea can be extended to combine parameters model and state prediction. After the comparison and the discussion, it is proved that RUL prediction using particle filter method provides a more accurate PH with a 26 cycles. Although we used in the training step only 48% of measured data. It should be noted, this method is able to be applied on any HSSB in wind farm, with an initial value parameter setting (eg; failure threshold).

5. Discussion and comparison with some previous works

Concerning resent work and according to our bibliographic study, this chapter presents a pedagogic implementation of a hybrid prognostic approach based on particle filter for bearing PHM. We encourage all researchers to work on this to have universal approach for bearing prognosis and complete one of the aims of industry 4.0 challenges. In [11], authors have applied the same data of this chapter the RUL prediction based on SVR. A smaller estimation error was found in 60% of training data compared to using 40% of training data. The SVR process is considered as internal RUL prediction. In addition, the SVR model parameters prediction was done after reaching 60% of degradation and that cannot be done online. It is impossible to define online the time where the degradation reaches 60%. In addition, it is hard to build the SVR model and validate it before the recording of the next raw vibration data. Some specific systems need to generate RUL prediction in little times due to the short lifetime of the used bearing and the required precision and excellence such as in robotics or nuclear application. In [5] the proposed Elman Neural Network (ENN) is motivated by a feature extraction from raw vibration data. The feature reduction is considered very important, as non-informative feature will be then discarded. Therefore, the online computational time will be reduced and ENN converge can be easily reached. Consequently, selecting suitable features is a prerequisite for accurate prognostics. ENN based prognosis is powerful but the implementation is costly and complicated. PF used in this chapter is a powerful tool for bearing failure prognosis; the implementation is very easy for and do not require a large data. The main thing to have the best RUL prediction is to build the right exponential degradation model with the true initial parameters, which can make an online PHM. As shown in Figure 13 the exponential curve trend the HI over 50 days with imperceptible fluctuations see Table 4.

6. Conclusion

This chapter introduces how to use approach based on condition monitoring data on HSSB for WTG. Classical statistical features and features derived from SK were used to elicit the bearing health state from the raw vibration data. Trendability, Monotonicity, Prognosability and Suitability are used as metric indices to obtain the corresponding feature for training step. The selected feature is used as Health index for the two proposed prognosis approach in order to predict the best RUL with higher performances.

The acquired results indicate that the Particle filter is more feasible tool for HSSB RUL prediction where the error equal to 7.69% and the degradation model with estimated parameters presents better trends for HI, compared to existing works.

As future work, the proposed method needs to be evaluated in a large amount of HSSB over a very long period. Also, the investigation of time-frequency-domain features will be considered. In addition, we invite next work to focus on external prognostic and to investigate new methodology or adopt some existing ones for dynamic feature selection. This ensures more alignment with industry 4.0 requirements.

Acknowledgments

The authors of this chapter would like to thank green power monitoring systems (GPMS) in USA for the permission to use their bearing data.