A hybrid prognostics approach to estimate the residual useful life of a planetary gearbox with a local defect
Zhe Cheng^{1}
^{1}Science and Technology on Integrated Logistics Support Laboratory, National University of Defense Technology, Changsha, P. R. China
^{1}Corresponding author
Journal of Vibroengineering, Vol. 17, Issue 2, 2015, p. 682694.
Received 9 June 2014; received in revised form 16 July 2014; accepted 7 September 2014; published 31 March 2015
Copyright © 2015 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
JVE Conferences
A hybrid prognostics approach for the monioring of a planetary gearbox with the local defect is presented. This hybrid method can predict the remaining useful life (RUL) of planetary gearbox with a fatigue crack. The method consists of a dynamical model for simulation data generation, a statistical algorithm for feature selection and weighting, and a modified grey model for RUL prediction. Experimental studies are conducted to validate and demonstrate the feasibility of the proposed method for RUL prediction of a cracked sun gear in planetary gearbox. And the validation has a promising result.
Keywords: planetary gearbox, defect, dynamical model, feature selection and weighing, residual useful life, prognosis.
1. Introduction
Engine and drive train failures are the leading causes of mishaps in rotorcrafts. As a critical component in main transmission gearbox of rotorcraft, the defect in planetary gear set is widely researched. One common defect is gear tooth crack caused by cyclic stress fatigue, which may lead to gear tooth breakage or complete failure of main transmission gearbox. However it is particularly challenging to detect and predict the remaining useful life (RUL) of these components, as the gear tooth fault of planetary gear set is located deep inside the main transmission gearbox.
As failure prognosis for planetary gearbox/rotorcraft transmission system has been studied widely, there are a lot of papers published in several key journals and at conferences. The common research papers in this area are datadriven approaches. Heimes [1] utilized advanced recurrent neural network architecture to estimate the RUL of the system. Peel [2] used a data driven approach of multiple neural network model to estimate the RUL left of an unspecified complex system. Orchard [3, 4] introduced a methodology for accurate and precise prediction of a failing component based on particle filtering and learning strategies. Swanson [5] used Kalman filters to track changes in relevant features like vibration levels, mode frequencies, or other waveform signature features, and then used the predicted states to estimate the future failure hazard, probability of survival, and RUL in an automated and objective methodology. Tse [6] used recurrent neural networks to forecast the rate of machine deterioration.
Physicsbased methods have also been used to predict the failure time in mechanical systems. Ray [7] presents a stochastic model of fatigue crack damage in metallic materials to predict failure, and this method was validated by KL analysis of test data. Li [8] presented a hybrid modelbased method to predict the RUL of a gear with a fatigue crack. Experimental studies were conducted to validate and demonstrate the feasibility of the proposed method for prognosis of a cracked gear. Kacprzynski [9, 10] utilized stochastic, physicsof failure models controlled failures of single spur gear teeth and then on a helical gear contained within a drive train system. The method was validated with transitional runtofailure data developed at Penn State ARL. Gu [11] presented a physicsoffailure based prognostics and health management approach for effective reliability prediction. The implementation procedure of this approach included failure modes, mechanisms, and effects analysis, data reduction and feature extraction from the lifecycle loads, damage accumulation, and assessment of uncertainty.
Generally, the methods aforementioned have their own shortages: 1) datadriven methods cannot provide sufficient guarantee about their prediction accuracy and definite physical meaning on the relationship between the index and damage; 2) physicsbased approaches need high computational costs associated with the physical models and these methods are only suitable for offline applications; 3) the data sets used for the prognosis of planetary gearbox include strong environmental noise components and many frequency components of other moving parts.
In view of the above analysis, a hybrid prognostics approach is presented for planetary gearbox. This method is composed of a dynamical model for vibration response analysis, a twostep statistical algorithm for feature selection and a modified grey model for RUL prediction. In comparison with the existing methods reported above, the proposed method harnesses the merits that the computational costs of analytical model are lower and the relationship between indicator and fault is easy to understand. This method has better prediction accuracy compared to the data driven method. Damage seeded tests and runtofailure test are conducted to test the performance of the proposed method.
2. The dynamical model of planetary gearbox with defect
The model of the gear tooth failure can be used to analyze the dynamic features in order to provide new signal source and suitable tools for predicting these failures. In the models, gear tooth defect is commonly assessed by the variation of dynamic parameters due to defect, such as stiffness and damping. The finite element analysis (FEA) and the analytical model are the general approach to simulate the defect seeded in part [1220]. Considering that FEA requires much computations time, this research will focus on the tooth stiffness reduction due to defect based on analytical model of planetary gear set.
2.1. The dynamical model of the healthy planetary gearbox
Due to the multiple components and its orbital motion of the planets, the epicyclic stage of transmission system is more complex than the ordinary gear train. As the internal ring gear is fixed to the top of the transmission casing and the planet gear support is very rigid, the central displacement of the ring gear and the planet gears in the planetary gearbox can be ignored to simplify the model. Then a lumped parameter dynamical formulation is employed to develop the dynamical model of a common planetary gearbox. Deduced from Lagrange equations, the dynamical differential equations of this planetary gearbox could be listed as:
where ${\theta}_{s}$, ${\theta}_{pi}$ and ${\theta}_{c}$ denote the rotation angle of sungear, planet gear and carrier respectively. ${T}_{D}$ and ${T}_{L}$ denote the driving torque and the loading torque respectively. $r$ denotes radius and ${r}_{b}$ is the radius of basic circle. $m$ is mass, $I$ is the rotational inertia of subcomponents, for planet gear, sun gear and carrier, $I=1/2m{{r}_{b}}^{2}$. $N$ is the number of planet gear and $\alpha $ is the mesh angle. The subscripts $r$, $s$, ${p}_{i}$ and $c$ denote ring gear, sun gear, the $i$th planet gear and carrier. $D$ and $P$ is the adhesive engaging force and the elastic engaging force respectively.
Eq. (1) is positive semidefinite, nonlinear equations, which have $N+$2 degreesoffreedom (dof’s). To translate the angle displacements of rigid bodies into the relative linear displacement, the relative displacement between the $i$th planet gear ${x}_{spi}$ and the sun gear, and the relative displacement between the carrier and the sun gear are defined as:
To simplify the model, the manufacturing errors and the mesh errors in the model are ignored as can be referred in [15]. After that, the internal meshing side clearance and the external meshing side clearance are defined as 2${b}_{spi}$ and 2${b}_{rpi}$, the adhesive engaging force $D$ and the elastic engaging force $P$ can be presented as:
where ${K}_{spi}$ denotes the mesh stiffness between the planet gear and the sun gear, ${K}_{rpi}$ denotes the mesh stiffness between the internal ring gear and the planet gear; ${C}_{spi}$ denotes the mesh damping between the planet gear and the sun gear; ${C}_{rpi}$ denotes the mesh damping between the ring gear and the planet gear. $K\left(t\right)$ is time varying mesh stiffness, and $f(x,b)$ is nonlinear clearance function, which is defined as:
where $b$ is the clearance constant.
After that, Eq. (2) and Eq. (3) are substituted into Eq. (1) and simplifying the equations, the dynamic model of 2KH planetary gear set can be obtained. The equivalent mass is denoted as $M$, for the planet gear and the sun gear, $M=I/{r}_{b}^{2}\text{,}$ for the carrier, ${M}_{c}={I}_{c}/{{r}_{bc}}^{2}+\sum _{i=1}^{N}{(m}_{pi}/\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\alpha )$.
2.2. The dynamical model of planetary gearbox with fatigue crack seeded
Fatigue crack is one of the common damages in planetary gearbox, and the cause of this damage is cyclic stressing of the gear tooth beyond its fatigue limit. Generally, the bending fatigue crack starts in the root section and it progresses until the tooth or part of it breaks off. The common fatigue crack occurs in the usual tooth root fillet section, so the tooth fatigue crack of sun gear is seeded in the model of planetary gear set in this section, as can be seen in Fig. 1.
Fig. 1. Gearmesh stiffness variation of the gear pair with tooth crack
a)
b)
To simplify the model of damage, the shape of tooth root crack is approximated with straight lines, which are defined by the length l and the direction angle $\varsigma $ of the crack, as shown in Fig. 1. The level of crack is defined as $s$, which is determined by l and $\varsigma $, $s=s({l}_{j},\varsigma )$. For the mathematical tractability to solve the dynamical model, the crack size and its threshold value are defined as ${c}_{j}$ and ${c}_{f}\text{,}$${c}_{j}=2{m}^{\mathrm{\text{'}}}\mathrm{c}\mathrm{o}\mathrm{t}\varsigma {l}_{j}\mathrm{s}\mathrm{i}\mathrm{n}\varsigma \text{,}$${\stackrel{~}{c}}_{j}={c}_{j}/{c}_{f}\text{,}$$j=$ 1, 2, 3,…. And then ${s}_{j}=s({l}_{j},\varsigma )=s\left({\stackrel{~}{c}}_{j}\right)$, ${s}_{j}\in \mathrm{}$[0 %, 100 %]. Referring [16, 17], the gear mesh stiffness of the gear pair is calculated by taking into account the geometric change due to the tooth crack. The details can be referred to [16].
The variation amplitude of gear mesh stiffness, which is caused by sun gear tooth crack, is defined as $\mathrm{\Delta}K({f}_{d},t)$, thus the time varying mesh stiffness of gear pair with tooth crack can be expressed as:
where ${K}_{0}\left({f}_{0},t\right)$ is the mesh stiffness of the healthy gear pair, ${f}_{0}$ is the mesh frequency, and ${f}_{d}$ is the meshing frequency of the gear pair with defect. ${f}_{d}$ is the function of the rotary frequency of sun gear ${f}_{s}$, the rotary frequency of carrier ${f}_{c}$ and the number of planet gears $N$, ${f}_{d}=N({f}_{s}{f}_{c})$. The time vary mesh stiffness is illustrated in Fig. 1(b). The dynamical model of planetary gearbox with sun gear tooth crack is acquired by substituting Eq. (5) into the model deduced above.
2.3. Simulation of the dynamical models
The parameters in the dynamical model above should be set up according to Table 1. Parameter s is varied as [0 %:5 %:100 %] for the defect seeded model. Matlab 7.1 is used to execute this simulation, and fourorder RungeKutta method is selected as the solution for the models. The step size and the duration for the solution procedure are set up to 0.0001 second and 1 second respectively.
Table 1. Parameter value in the models
Parameter (Unit)

Value

Parameter (Unit)

Value

Modulus (mm)

2.5

${r}_{c}$ (m)

0.108

Tooth number of sun gear

28

Driving torque (N⋅m)

100

Tooth number of planet gear

32

Loading torque (N⋅m)

220

Tooth number of ring gear

92

${K}_{spi}$ (N/m)

1.806×10^{9}

Number of planet gear

4

${K}_{rpi}$ (N/m)

2.212×10^{9}

Tooth width (mm)

12

${C}_{spi}$ (N⋅s/m)

4.474×10^{3}

Pressure angle (Deg)

20

${C}_{rpi}$ (N⋅s/m)

5.764×10^{3}

${m}_{s}$ (kg)

1.144

${b}_{spi}$, ${b}_{rpi}$ (m)

0.0001

${m}_{pi}$ (kg)

1.757

${I}_{pi}$ (m)

0.0016

${m}_{r}$ (kg)

3.016

${I}_{s}$ (m)

0.0011

${r}_{bs}$ (m)

0.046

${I}_{c}$ (m)

0.032

${r}_{bpi}$ (m)

0.057

Material

40Cr

The dynamic responses of the dynamical model in both the time domain and the frequency domain are shown in Fig. 2, which is for a healthy planetary gearbox ($s=$ 0 %) and 3 different levels of crack ($s=$ 20 %, 40 %, 80 %). The healthy model is characterized by the dominance of the gear mesh frequency, denoted by 1$X$ and its harmonics of 2$X$6$X$, and etc. For the defect seeded models, the amplitude modulation in time domain of the gear mesh signal can be clearly observed, as can be seen in Fig. 2(a). As a consequence, many new frequency components appear at the left side of the dominant frequency (near 1/2$X$). The amplitude of the new frequency components increases with the growth of crack level, but the amplitude of dominant frequency and its harmonics decreases at the same time, as can be seen in Fig. 2(b). These features are very close to the experimental observations in the previous study literatures [16, 17, 21].
Fig. 2. Dynamical response of the models
a) Time domain
b) Frequency domain
2.4. Validation of the dynamical models
The test data sets in health condition and damage seeded condition are used to validate the dynamical model deduced in this paper. Above all, Root Mean Square (RMS) of each data set is calculated to reduce the effect of environment noise as the aid of the common algorithm TSA, and then the curves of error percentage of RMS between test data and simulation data of the dynamical models are depicted in Fig. 4. The real line is the error percentage between test data set and simulation data set in health condition, which is about 5 %, and the imaginal line is for damage seeded condition, which is about 15 %, as can be seen in Fig. 3. As the reason of that the damage seeded in the dynamical model above is simplified, which is different from the damage seeded in test. As a result, the imaginal line is above the real line obviously. Generally, the error percentage is low and the accuracy of dynamical model would be improved by considering more nonlinear factors of the damage seeded in the model.
Fig. 3. Validation results of dynamical models
3. Feature selection and weighting
From our literature review, 27 features have been used in different cases of gearbox condition monitoring. The features are organised into four groups and assigned with serial numbers.
1) Features in time domain are used most frequently in gearbox diagnosis [21]. They include root mean squared (RMS), crest factor (CF), energy ratio (ER), kurtosis, standard deviation, energy operator, absolute mean value, clearance factor and impulse factor, which are assigned with the serial numbers 1 to 9 respectively for the ease of identification in the process of feature selection.
2) There are many other traditional statistical feature parameters, such as FM0, FM4, NA4, M6A, M8A, NB4, NA4*, NB4*, M6A* and M8A*, generally used in planetary gearbox condition monitoring [22]. The serial numbers of these features are 10 to 19.
3) Other kinds of feature parameters based on frequency spectrum of vibration signal are widely used to detect and diagnose faults in rotorcraft power trains, such as mean frequency (MF), frequency centre (FC), root mean square frequency (RMSF) and standard deviation frequency (STDF) [23]. The serial numbers of these features are 20 to 23.
4) Other than the features presented above, there are some important features which have been validated in literatures, which are named IntraRevolution Energy Variance (IREV) [24], spectrum kurtosis (SK) [25, 26], local spectrum kurtosis (LSK) [27], NSR [28]. The serial numbers of these features are 24 to 27.
The simulation data sets generated by the dynamical model above are utilized to calculate all the feature parameters listed above.
3.1. Feature selection
An optimal feature suitable for prognostics should have three merits: 1) Sensitivity, it means that the feature has a optimal classification distance; 2) Stability, which means the feature has the stable classification performance in different operation conditions (including loading and rotational speed). The target features will be selected from the 27 feature parameters above.
A statistic algrithm named twosample Ztest is commonly used to measure the distance of two class samples [24, 27]. In this research, this algrithm is applied as sensitivity analysis algorithm, and then it is modified to analyze the stability and the relationality of feature parameters. The feature selection procedure can be described as follows:
1) Sensitivity analysis.
For the $k$th feature parameter, calculating the classification distance ${Z}_{k}$ of the healthy samples and the fault seeded samples under the identical operating condition by two samples Ztest procedure:
where ${\mathbf{X}}_{1}$ and ${\mathbf{X}}_{2}$ are the healthy sample set and the fault seeded sample set, ${S}_{\mathbf{X}}$ and $\stackrel{}{\mathbf{X}}$ are the standard deviation and the mean of $\mathbf{X}$, respectively. And $n$ is the number of samples for each sample set. $k=$1, 2,...,$K$, $K$ is the number of feature parameters above. As the value of ${Z}_{k}$ increases, the sensitivity of the kth feature parameter will grow, that means the value of the $k$th feature parameter will have an obvious change while there is damage.
2) Stability analysis.
The classification distance of the kth feature parameter under the jth condition is defined as ${Z}_{k,j}$, $j=$1, 2,..., $J$. Then, the similarity ratio of classification distance for different condition, which is named stability in this research, is calculated:
where ${\mathbf{Z}}_{k}$ is the classification distance vector of the $k$th feature parameter for different condition, ${Z}_{k,m}$ is the mth element of ${\mathbf{Z}}_{k}$, ${S}_{{\mathbf{Z}}_{k}}$ is the standard deviation of ${\mathbf{Z}}_{k}$, $J$ is the number of operation conditions, $n$ is the element number of ${\mathbf{Z}}_{k}$. The numerator in Eq. (7) shows the general statistical characteristic of ${\mathbf{Z}}_{k}\text{,}$ and the denominator in Eq. (7) denotes the statistical diversity of classification distance of the $k$th feature parameter for different operation conditions. So ${S}_{k}$ indicates the stability of the kth feature parameter to different operation condition. As the value of ${S}_{k}$ increases, the stability of the kth feature parameter will rise, which means the $k$th feature parameter will keep stable while the operation condition (e.g. rotary speed and loading) changes.
According to the feature selection algorithms in this section, the feature parameter suitable for prognosis can be selected from the potential feature set when their evaluation criterion is predefined.
3.2. Feature weighting
As the feature vector includes more than one feature parameters, the feature weighting is necessary for fusing the analysis results of all the parameters in feature vector to achieve more reliable results. In this research, the weighting of feature parameter is determined based on the sensitivity and stability performance of the feature parameter calculated above. And also note that the weighting sum of all the feature parameters is equal to 1. Thus the parameters in feature vector can be weighted as follows:
where $i$ is the serial number of the feature parameter, $N$ is the number of feature parameter in the target vector.
4. RUL estimation based on modified grey model
4.1. Prognosis approach based on grey model
The simplest and most elementary form of grey model is GM(1,1) [31]. It is proposed as a forecasting model to extrapolate the original data sequence. The concept of the grey theory was first derived by Deng in 1986 [31]. It describes GM(1,1) model as approximations via leastsquare estimations at the firstorder differential equation level. GM(1,1) has been widely used in the prediction of mechanical systems [3138].
Vibration signal is an important information resource for monitoring rotary mechanicals. With the propagation of crack, the value of vibration indicators will increase generally. This process has two characteristics: 1) $X\left({t}_{i}\right)\ge $ 0 and 2) $X\left({t}_{i}\right)\ge X\left({t}_{i1}\right)$, which show the data set has strong tendency. These two characteristics meet the need to build grey model (GM), so GM(1,1) can be used to predict the future condition of planetary gearbox. The prediction procedure is as bellows:
Define ${x}^{\left(0\right)}=\left({x}^{\left(0\right)}\left(1\right),{x}^{\left(0\right)}\left(2\right),\cdot \cdot \cdot ,{x}^{\left(0\right)}\left(n\right)\right)$ as the feature condition series, ${x}^{\left(1\right)}$ is the accumulated generating operation (AGO) of ${x}^{\left(0\right)}$, ${x}^{\left(1\right)}\left(k\right)=\sum _{i=1}^{k}{x}^{\left(0\right)}\left(i\right)$, $k\in \left\{1,\cdot \cdot \cdot ,n\right\}$; ${z}^{\left(1\right)}$is the mean generating operation of ${x}^{\left(1\right)}\text{,}$${z}^{\left(1\right)}\left(k\right)=0.5{x}^{\left(1\right)}\left(k\right)+0.5{x}^{\left(1\right)}\left(k1\right)\text{,}$$k\in \left\{2,\cdot \cdot \cdot ,n\right\}\text{.}$ GM(1,1) is:
The discrete model for prediction is:
where:
4.2. RUL estimation based on the modified grey model
To improve the prediction accuracy of the grey model in Eq. (10), a residual efficient $\epsilon $ is added to the grey input $u$ of the grey model for the addition of novel information from the physicsbased models:
Then the discrete form is:
The residual efficient of feature ${F}_{i}$ is defined as ${\epsilon}_{i}$. After the feature parameters used for RUL prediction are determined, $\epsilon $ could be calculated based on the given information of simulation signal from dynamical model developed above:
where ${\widehat{x}}^{\left(1\right)}\left(1\right)$ and ${\widehat{x}}^{\left(1\right)}\left(k+1\right)$ are AGO of the time series of feature parameter $F$, and the elements of the time series are calculated based on simulation signal.
The onestep and multistep prediction result of ${F}_{i}$ is ${x}_{i}^{\left(1\right)}\left(k+1\right)$, which could be obtained based on test data or real data according to Eq. (13). After that the value of prediction results ${x}_{i}^{\left(1\right)}\left(k+1\right)$ are fused utilizing the feature weight vector ${W}_{i}$ according to Eq. (15):
where ${x}_{i}^{\left(1\right)}\left(k+1\right)$ is the prediction result of ${F}_{i}$, and ${W}_{i}$ is the weight of ${F}_{i}$. Then ${x}_{}^{\left(1\right)}\left(k+1\right)$ in Eq. (15) is the expected value of future condition prediction.
In this research, RUL ${N}_{f}$ is expressed as the stress cycle of gear crack tooth $n$. ${N}_{f}$ is corresponding to ${x}_{}^{\left(1\right)}\left(k+1\right)$ according to Eq. (16):
where $\xi $ is determined based on the historical data of damage evolution test.
4.3. Description of the RUL estimation result
Logarithmic normal distribution (LND) is often used to describe the prognosis result of mechanical parts. In this research, the RUL probability distribution function of tooth fatigue crack in planetary gearbox is considered to belong to LND. The common form of LND is twoparameterdistribution, the density function is:
where $t$ is time, $\sigma $ and $\mu $ are the standard deviation and the mean, respectively.
As there are a lot of uncertainty issues in equipment operation, the accuracy of prediction result will decrease quickly while the prediction time extends, and then the variance of RUL distribution will grow. By iteration of modified grey model, multistep prediction will be carried. Variance $\sigma $ could be expressed as the function of prediction life step $n$, $\sigma =\sigma \left(n\right)=\sigma {e}^{n1}$. The RUL distribution is obtained by substituting the expressions of variance and expectation into Eq. (17):
Eq. (15) and Eq. (18) are combined, and then the prognosis result of tooth root fatigue crack in planetary gear sets is generated.
The relationship of RUL $N$ and failure time $T$ is correlated with operation profile in test/field, such as the rotary speed variation.
To evaluate the prognosis results, the prediction error $\eta $ is defined as the difference between the actual failure time ${T}_{af}$ and the predicted failure time ${T}_{pf}$:
where ${T}_{af}$ is the actual failure time, ${T}_{pf}$ is the predicted failure time and ${T}_{c}$ is the current time.
5. Validation with test data
5.1. Test rig and the tooth crack propagation test
As can be seen in Fig. 4(a), the test rig provides torque and speed reduction through three stages, which simulates the main transmission system of general rotorcraft. The second stage is a planetary gearbox, of which the internal gear rig is splined to the top of gearbox casing.
To validate the prognostics approach in this research, a number of experiments have been conducted. The vibration signal generated by the planetary gearbox is picked up by an accelerometer bolted on the top of planetary gearbox casing and the electrical signal is transferred to the data acquisition system. The sampling frequency $fs$ is set up to 10 kHz. The signal is lowpass filtered at 2 kHz through a 4th order filter to limit aliasing distortion and retain waveform integrity of the effective signal. Data is stored in a PC for post processing. A number of 204,800 data points has been acquired in tests corresponding to a timehistory length of 20 s.
In this research, an incipient crack is seeded in the sun gear, as can be seen in Fig. 4(b), and a runtofailure test is carried out. In the test, the interval of sampling is 10 min. For every 50 minutes, the test is paused and the actual size of fatigue crack in propagation is checked. When the fatigue crack level reaches the threshold value, the test is ended. The threshold value of crack level is set up according to the experience of the prior test. The total run time of this crack seeded test is 736 min. From the test, the runtofailure transitional data of fatigue crack propagation in sun gear of planetary gearbox is acquired.
Fig. 4. Configuration of the test rig and the sungear with tooth root crack
a) The test rig
b) The sungear with tooth root crack
5.2. Validation results with test data
Based on the simulation data of the dynamical models, feature selection and weighting has been carried out before the test data validation. Only a few valueble features are kept and used, as can be seen in Fig. 5. As there is no prior knowledge of setting the threshold, a median value is used as selection criteria. In this research, [20, 20] is selected to be the threshold vector to illustrate the proposed methodology.
Fig. 5. Feature selection results
a)
b)
In the feature selection procedure, the feature parameters are calculated based on the simulation data sets of the dynamical models with diferent level damge seeded. The procedure of the feature selection could be clearly dipicted in Fig. 5. The first step is the sensetivity analysis of feature parameters, the threshold value for this step is set up to 20, and the feature parameters #10, #11, #12, #13, #14, #16, #17, #18, #19, #22, #24, #25, #26, #27 are in the second step for stability analysis, as shown in Fig. 5(a). In the second step, the threshold value is set up to 20, and the selection result of this step is the feature parameters #25, #26 and #27, as shown in Fig. 5(a).
The feature parameters #25, #26 and #27 have been selected to be the remaining features, as shown in Fig. 5(b). These features are SK, LSK and NSR, and they are labeled ${F}_{1}$, ${F}_{2}$ and ${F}_{3}$ in this paper for convenience. Following the procedure in Eq. (9), the weights of these selected features are 0.2752, 0.2013 and 0.5235 respectively.
To validate the proposed RUL prediction approach in this research, the test data sets ${X}_{1}$~${X}_{4}$ at time 360 min and 460 min are used for analysis. 50 feature samples of each test data set are extracted. The RUL for the test data at $t=$ 360 min and $t=$ 460 min is predicted. The prediction result of test data is recorded in Table 2, where ${F}_{0}$ denotes the healthy case, ${F}_{t}$ denotes the failure case, ${F}_{A}$ is the prediction result at $t=$ 360 min, ${F}_{B}$ is the prediction result at $t=$ 460 min. The RUL at $t=$ 360 min is 681 min with the confidence of 95 %. And the RUL at $t=$ 460 min is 705 min with the confidence of 95 %, while the real RUL of this test is 736 min.
According to Eq. (17), the prediction accuracy at $t=$ 360 min is 85.37 % and at $t=$ 460 min is 88.77 %. It could be seen from these prediction results that, as the prediction time reaches failure time, the accuracy of prediction result will be improved, in other words, the expected RUL of prediction result will approximate the real RUL, while the variance of prediction result decreases quickly. Compared to the traditional GM(1,1), the modified GM(1,1) has better accuracy and lower prediction error, which can be seen in Table 2.
Table 2. Comparation of the prediction results for the modified GM and the traditional model
Model type

Prediction start time (min)

Failure time for 95 % confidence (min)

Real failure time (min)

Prediction error (min)

Accuracy (%)

Traditional GM(1,1)

360

573

736

163

56.65

460

621

115

58.33


Modified GM(1,1)

360

681

54

85.37


460

705

31

88.77

6. Conclusions
A hybrid RUL esitimation methodology for planetary gearbox is devised by integrating dynamical models for damage simulation signal generation, a twostep statistical algorithm for feature selection and weighting, and a modified grey model for prediction, to forecast the RUL of a cracked sungear in planetary gearbox. The proposed method is validated with the sun gear runtofailure test. The experimental study results indicate that this prognostics approach has high accuracy to predict RUL of sun gear tooth fatigue crack in planetary gearbox. However, the RUL prediction of this approach is a little conservative. It is suggested that the method should improve the model parameter of dynamical model and GM(1,1) for certain applications.
Acknowledgements
This research is supported by the National Natural Sciences Foundation of China (Grant No. 51205401) and Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20124307120010) and Foundation for the Author of Provincial Excellent Doctoral Dissertation of Hunan Province of China (Grant No. YB2013B002). Valuable comments on the paper from anonymous reviewers are very much appreciated.
References
 Heimes F. Recurrent neural networks for remaining useful life estimation. International Conference on Prognostics and Health Management, CA, USA, 2008. [CrossRef]
 Peel L. Data driven prognostics using a kalman filter ensemble of neural network models. International Conference on Prognostics and Health Management, CA, USA, 2008. [CrossRef]
 Orchard M. A Particle filtering framework for failure prognosis. Proceedings of World Tribology Congress, Washington, D.C., USA, 2005. [CrossRef]
 Orchard M. Advances in uncertainty representation and management for particle filtering applied to prognostics. International Conference on Prognostics and Health Management, CA, USA, 2008. [CrossRef]
 Swanson D. A general prognostic tracking algorithm for predictive maintenance. IEEE Aerospace Conference, Big Sky, 2001, p. 29712977. [CrossRef]
 Tse P., Atherton D. Prediction of machine deterioration using vibration based fault trends and recurrent neural networks. Journal of Vibration and Acoustics, Vol. 121, Issue 3, 1999, p. 355362. [CrossRef]
 Ray A. Stochastic modeling of fatigue crack damage for risk analysis and remaining life prediction. Journal of Dynamic System, Measurement and Control, Vol. 121, Issue 3, 1999, p. 386393. [CrossRef]
 Li C., Lee H. Gear fatigue crack prognosis using embedded model, gear dynamic model and fracture mechanics. Mechanical Systems and Signal Processing, Vol. 19, Issue 4, 2005, p. 836846. [CrossRef]
 Kacprzynski G., Sarlashkar A., Roemer M., etc. Predicting remaining life by fusing the physics of failure modeling with diagnostics. Journal of Materials, Vol. 56, Issue 3, 2004, p. 2935. [CrossRef]
 Kacprzynski G., Maynard K., Roemer M. Enhancement of physicsoffailure prognostic models with system level features. Aerospace Conference Proceedings, 2002. [CrossRef]
 Gu J., Pecht M. Prognostics and health management using physicsoffailure. Proceedings of the Annual Reliability and Maintainability Symposium, 2008. [CrossRef]
 Pimsarn M., Kazerounian K. Efficient evaluation of spur gear tooth mesh load using pseudointerference stiffness estimation method. Mechanism and Machine Theory, Vol. 37, Issue 8, 2002, p. 769786. [CrossRef]
 Wang J. Numerical and experimental analysis of spur gears in mesh. Ph.D., Curtin University of Technology, 2003. [CrossRef]
 Ambarisha V., Parker R. Nonlinear dynamics of planetary gears using analytical and finite element models. Journal of Sound and Vibration, Vol. 302, Issue 3, 2007, p. 577595. [CrossRef]
 Chaari F., Fakhfakh T., Haddar M. Dynamic analysis of a planetary gear failure caused by tooth pitting and cracking. Journal of Failure Analysis and Prevention, Vol. 6, Issue 2, 2006, p. 3944. [CrossRef]
 Chaari F., Baccar W. Effect of spalling or tooth breakage on gearmesh stiffness and dynamic response of a onestage spur gear transmission. European Journal of Mechanics A/Solids, Vol. 27, Issue 4, 2008, p. 691705. [CrossRef]
 Hbaieb R., Chaari F., Fakhfakh T., Haddar M. Influence of eccentricity, profile error and tooth pitting on helical planetary gear vibration. Machine Dynamics Problems, Vol. 29, Issue 1, 2005, p. 532. [CrossRef]
 Walha L., Louati J., Fakhfakh T., Haddar M. Dynamic response of two stages gear system damaged by teeth defects. Machine Dynamics Problems, Vol. 29, Issue 1, 2005, p. 107124. [CrossRef]
 Choy F., Polychuk V., Zakarajsek J., Handschuh R., Townsend D. Analysis of the effects of surface pitting and wear on the vibrations of a gear transmission system. Tribology, Vol. 29, Issue 1, 1996, p. 7783. [CrossRef]
 Yuksel C., Kahraman A. Dynamic tooth loads of planetary gear sets having tooth profile wear. Mechanism and Machine Theory, Vol. 39, Issue 7, 2004, p. 695715. [CrossRef]
 Dempsey P., Lewicki D., Le D. Investigation of current methods to identify helicopter gear health. Aerospace Conference, Big Sky, Montana, 2007. [CrossRef]
 Samuel P., Pines D. A review of vibrationbased techniques for helicopter transmission diagnostics. Journal of Sound and Vibration, Vol. 282, Issue 1, 2005, p. 475508. [CrossRef]
 Yang B., Kim K. Application of dempstershafer theory in fault diagnosis of induction motors using vibration and current signals. Mechanical Systems and Signal Processing, Vol. 20, Issue 2, 2006, p. 403420. [CrossRef]
 Wu B., Saxena A. An approach to fault diagnosis of helicopter planetary gears. IEEE Autotestcon, 2004, p. 475481. [CrossRef]
 Barszcz T., Randall R. Application of spectral kurtosis for detection of a tooth crack in the planetary gear of a wind turbine. Mechanical Systems and Signal Processing, Vol. 23, Issue 4, 2009, p. 13521365. [CrossRef]
 Antoni J., Randall R. The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines. Mechanical Systems and Signal Processing, Vol. 20, Issue 2, 2006, p. 308331. [CrossRef]
 Cheng Z., Hu N., Gao J. An approach to detect damage quantitatively for planetary gear sets based on physical models. Advanced Science Letters, Vol. 4, Issue 4, 2011, p. 16951701. [CrossRef]
 Cheng Z., Hu N. Application of GRA in damage detection and mode identification of planetary gear set. Journal of Grey System, Vol. 23, Issue 4, 2011, p. 315326. [CrossRef]
 Wen K. Grey systems: Modeling and prediction. Yang’s Scientific Research Institute, Tucson, 2004. [CrossRef]
 Wen K., Changchien S., Yeh C., Wang C., Lin H. Apply MATLAB in Grey System Theory. CHWA Publisher, Taipei, 2006. [CrossRef]
 Deng J. Introduction to grey system theory. Journal of Grey System, Vol. 1, Issue 1, 1989, p. 124. [CrossRef]
 Deng J. Control problems of grey system. System Control Letters, Vol. 1, Issue 5, 1982, p. 285294. [CrossRef]
 Deng J. A Novel grey model GM(1,1τ,r): Generalizing GM(1,1). Journal of Grey System, Vol. 13, Issue 1, 2001, p. 18. [CrossRef]
 Lin S., Wu S. An Intelligent Webbased GRA/Cointegration Analysis for Systematic Risk. International Journal of Computers, Vol. 4, Issue 4, 2010, p. 223234. [CrossRef]
 Wang Y., Liao R., Chen W., Du L. A online forecast method for predicting the gasinoil concentrations in oillfilled transformer. International Conference on Condition Monitoring and Diagnosis, Changwon, Korea, 2006. [CrossRef]
 Anjali C., Nirmal K. Applying grey theory prediction model on the dga data of the transformer oil and using it for fault diagnosis. Transactions on Power Systems, Vol. 4, Issue 1, 2009, p. 4352. [CrossRef]
 Cempel C. Forecasting the global and partial system condition by means of multidimensional condition monitoring methods. Journal of Theoretical and Applied Mechanics, Vol. 46, Issue 3, 2008, p. 777797. [CrossRef]
 Cempel C. Decomposition of the symptom observation matrix and grey forecasting in vibration condition monitoring of machines. International Journal of Applied Mathematics and Computation Science, Vol. 18, Issue 4, 2008, p. 569579. [CrossRef]
 Chen Z., Yang Y., Hu Z. Detecting and predicting early faults of complex rotating machinery based on cyclostationary time series model. Journal of Vibration and Acoustics, Vol. 128, Issue 5, 2006, p. 666671. [CrossRef]
 Cheng Z., Hu N., Zhang X. Crack level estimation approach for planetary gearbox based on simulation signal and GRA. Journal of Sound and Vibration, Vol. 331, Issue 26, 2012, p. 58535863. [CrossRef]