Bearing remain life prediction based on weighted complex SVM models

3.1. SVM prediction model

SVM is a machine learning tool that uses statistical learning theory to solve multi-dimensional functions. It is based on structural risk minimization principles, which overcomes the extra-learning problem of ANN.

The learning process of regression model LS-SVM is essentially a problem in quadratic programming. Given a set of data points ($x_1$,$y_1$)…($x_k$,$y_k$) such that $x_i\in R^n$ as input and $y_i\in R$ as target output, the regression problem is to find a function such as Eq. (3):

where $\varphi \left(x\right)$ is the high dimensional feature space, which is nonlinear mapped from the input space $x$. $w$ is the weight vector, $b$ the bias [24].

After training, the corresponding$y$can be found through $f\left(x\right)$ for the $x$ outside the sample. The$\mathrm{\epsilon }$-support vector regression ($\epsilon$-SVR) by Vapnik controls the precision of the algorithm through a specified tolerance error$\mathrm{\epsilon }$. The error of the sample is $\xi$, regardless of the loss, when $\left|\xi \right|\le \epsilon$; else consider the loss as $\left|\xi \right|-\epsilon$. First, map the sample into a high dimensional feature space by a nonlinear mapping function and convert the problem of the nonlinear function estimates into a linear regression problem in a high dimensional feature space. If we let $\phi \left(x\right)$ be the conversion function from the sample-space into the high dimension feature space, then the problem of solving the parameters of $f\left(x\right)$ is converted to solving an optimization problem Eq. (4) with the constraints in Eq. (5):

The feature space is one of high dimensionality and the target function is non-differentiable. In general, the LS-SVM regression problem is solved by establishing a Lagrange function, convert this problem to a dual optimization, i.e., problem Eq. (6) with constraint of Eq. (7) in order to determine the Lagrange multipliers ${\alpha }_i$, ${\alpha }_i^{*}$:

where ${\alpha }_i$, ${\alpha }_i^{*}$ are Lagrange multipliers and${\alpha }_i$, ${\alpha }_i^{*}\ge 0$. ${\alpha }_i\times {\alpha }_i^{*}=0$. $C$ evaluates the tradeoff between the empirical risk and the smoothness of the model.

The LS-SVM regression problem has therefore been transformed into a quadratic programming problem. The regression equation can be obtained by solving this problem. With the kernel function $K(x_i,x_j)$, the corresponding regression function is provided by Eq. (8):

where the kernel function $K(x_i,x_j)$ is an internal product of vectors $x_i$ and $x_j$ in feature spaces $\phi \left(x_i\right)$ and $\phi \left(x_j\right)$.

3.2. The LS-SVM model parameters selection

The particle swarm optimization algorithm was first proposed in 1995. It is an optimization method based on a set of particles whose coordinates are potential solutions in the search space. Particles in PSO will change their coordinates (their solutions) by migrating. During migrating, each particle adjusts its own coordinates based on its own past experience and other particles’ past experiences.

The PSO was chosen to optimize selected the SVM parameters through the following formula:

where the subscript $i$ represents the $i$th particle. $j$ represents the $j$-dimensional.

The subscript $t$ represents the$t$generation. $v_{ij}\left(t\right)$ is the velocity of the $i$th particle in the $t$th iteration; $x_{ij}\left(t\right)$ is the position of the $i$th particle; $p_{ij}\left(t\right)$ is the pbest position of the $i$th particle; $p_{g_j}$ is the gbest position (pbest represents the local optimum of the particles, gbest represents the overall situation optimum of the particles); The $w$ represents the inertia weight. $c_1$, $c_2$ are learning factors. $r_1~U\left(\mathrm{0,1}\right)$, $r_2~U\left(\mathrm{0,1}\right)$ represent two independent random functions.

The root-mean square error (RMSE) is used as the criteria to evaluate the proposed algorithm accuracy:

where $N$ represents the total number of data points in the test set, $y_i$ is actual value in training set or test set, and ${\widehat{y}}_i$ represents the predicted value of the model.

The process of optimizing the parameters $\gamma$, $\sigma$ based on the PSO is given below:

(1) At the beginning of the optimization process, randomly initialize population size, $c_1$, $c_2$, $\omega$, $rand\left(1\right)$, $rand\left(2\right)$, determine the termination condition, positions and velocities of the particle, mapping the LS-SVM parameters $\gamma$, $\sigma$ into a group of particles, initialize the initial position of each particle, pbest, gbest of the particles;

(2) When training the SVM, using the Eq. (11) as the PSO fitness function;

(3) Use the target parameters $\gamma$, $\sigma$ as the particles, use their initial values as the SVM parameters in step 2), the corresponding value of the fitness function as the optimal solution of the $\gamma$, $\sigma$;

(4) Use the initial error value of step 2) as the particle’s initial fitness value, search the optimal value as the global fitness value among the initial fitness value, and the corresponding particles as the current global optimal solution;

(5) Updating the velocity and position vector;

(6) Re-substituted the updated parameters $\gamma$, $\sigma$ into the SVM model, and re-training the SVM model according to the step 2), save the output value, calculate the fitness value of the particles again;

(7) Comparing the saved global fitness value gotten in step (6) with the current particle’s fitness value, if the global fitness value superior to the current particle’s fitness value, update the current particle’s fitness value according to the step (5), and update the current particle’s optimal value equal to the corresponding particle’s optimal value gotten in step (6);

(8) While the termination conditions are not met, return to step (5);

(9) End loop.

3.3. The prediction strategy and the structure of the SVM model

Traditional forecasting methods mainly achieve single-step prediction, when those methods are used for multi-steps prediction, they can’t get an overall development trend of the series. Multi-steps prediction method has the ability to obtain overall information of the series which provides the possibility for long-term prediction. There are two typical alternatives to build multi-steps life prediction model. One is iterated prediction and the other is direct prediction. The comparison of the two strategies can be found in a number of literatures [25]. Marcellino [26] presented a large-scale empirical comparison of iterated versus direct prediction. The results show that iterated prediction typically outperforms the direct prediction. So, the iterated multi-steps prediction strategy has numerous advantages and will be adopted in this paper.

In order to determine the structure of the SVM, we constructed a three layers SVM prediction model. But, to achieve the multi-steps time series life prediction there is a basic problem should be suppressed. That is how many essential observations (inputs) are used for forecasting the future value (The output node number is 1), so-called embedding dimension $d$. In order to suppressed the problem, the CAO method [27] which is particularly efficient to determine the minimum embedding dimension through the expansion of neighbor point in the embedding space, is employed to select an appropriate embedding dimension $d$. Then, the SVM input node number is determined.

To effective select an appropriate embedding dimension based on the CAO method, the phase space reconstruction method should be mentioned. The fundamental theorem of phase space reconstruction is pioneered by Takens [28]. For an N-point time series $X=\{x_1,x_2,\dots ,x_N\}$, a sequence of vectors $y_i$ in a new space can be generated as: $y_{i\left(d\right)}=\{x_i,x_{i+\tau },\cdots ,x_{i+(d-1)\tau }\}$, where $i=\mathrm{1,2},\dots ,N_d$, $N_d=N-(d-1)\tau$ is the length of the reconstructed vector $y_i$; $d$ is the embedding dimension of the reconstructed state space; and $\tau$ is embedding delay time. The time delay $\tau$ is chosen through the autocorrelation function [29]:

where $x_i^{\mathrm{\text{'}}}=x_i-\overline{x}$, $\overline{x}$ is the average value of the time series. The optimal time delay $\tau$ is determined when the first minimum value of $C\left(\tau \right)$ occurs.

The embedding dimension $d$ is chosen through CAO method, defining the quantity as follows:

where $‖•‖$ is the Euclidian distance and is given by the maximum norm. $y_i\left(d\right)$ means the $i$th reconstructed vector and $n(i,d)$ is an integer, so that $y_{n(i,d)}\left(d\right)$ is the nearest neighbor of $y_i\left(d\right)$ in the embedding dimension $d$. A new quantity is defined as the mean value of all $a(i,d)\mathrm{\text{'}}s$:

$E\left(d\right)$ is only dependent on the dimension $d$ and time delay $\tau$. To investigate its variation from $d$ to $d+$ 1, the parameter $E_1$ is given by:

By increasing the value of $d$, the value $E_1\left(d\right)$ is also increased and it stops increasing when the time-series comes from a deterministic process. If a plateau is observed for $d\ge d_0$, then $d_0+1$ is the minimum embedding dimension. But $E_1\left(d\right)$ has the problem of slowly increasing or has stopped changing if $d$ is sufficiently large. CAO introduced another quantity $E_2\left(d\right)$ to overcome the problem:

where:

Through CAO method, the embedding dimension$d$of the SVM prediction model is chosen. The structure of the SVM model is determined.

5.1. Case 1

The vibration signals used in this paper are provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati [30]. The experimental data sets are generated from bearing run-to-failure tests under constant load conditions on a specially designed test rig as shown in Fig. 2. The bearing test rig hosts four test Rexnord ZA-2115 double row bearings on one shaft. The shaft is driven by an AC motor and coupled by rub belts. The rotation speed is kept constant at 2000 rpm. A radial load of 6000 lbs is added to the shaft and bearing by a spring mechanism. A PCB 353B33 high sensitivity quartz ICP accelerometer is installed on each bearing housing. The NI 6062E card is used to collected the vibration data every 10 minutes. The data sampling rate is 20 kHz and the data length is 20,480 points as shown in Fig. 3. Data collection is done using a National Instruments LabVIEW program. Four bearings are used for full time accelerated life test, and 3 bearings recession data were used for training the SVM, getting 3 forecasting models and 1 bearing recession data for prediction.

The Time-domain, frequency domain, time-frequency domain methods are used to deal with the collected vibration data as described in Section 2. The extract features are showed in Figs. 4-6.

Fig. 2The test rig

Fig. 3The collected vibration signals

Fig. 4Sixteen indicators of the time domain: a) mean value; b) RMS; c) mean-square-root amplitude; d) absolute value; e) skewness; f) kurtosis; g) variance; h) the maximum value; i) the minimum value; j) the peak-to-peak value; k) index indictor; l) peak-to-peak indictor; m) pulse indictor; n) margin index; o) skewness index; p) kurtosis indicator

Fig. 5Twelve indicators of the frequency domain: a) the mean frequency; b) frequency standard deviation; c) the frequency characteristics 1; d) the frequency characteristics 2; e) frequency center; f) the frequency characteristics 3; g) RMS frequency; h) the frequency characteristics 4; i) the frequency characteristics 5; j) the frequency characteristics 6; k) the frequency characteristics 7; l) the frequency characteristics 8

Fig. 6Six IMF energy indicators of the time-frequency domain: a) IMF1 energy indicator; b) IMF2 energy indicator; c) IMF3 energy indicator; d) IMF4 energy indicator; e) IMF5 energy indicator; f) IMF6 energy indicator

From Figs. 4-6 we can see that, some features are not appropriate to reflect the bearing degradation process, for example, the feature (a), (e) in Fig. 4 still show as a straight line, where (l), (m), (n) features in Fig. 4 contain a lot of noise information, feature (d) in Fig. 5 and feature (e), (f), (g) in Fig. 6 didn’t sensitive to the early fault information. So, we need to extract the typical feature from these original features to better reflect the bearing degradation process,

The PCA is used to reduce the dimensionality of calculated features, in this research, the first principal component is chosen, in order to compare the effect of the recession indictor of the Kurtosis and the PCA, the result is shown in Fig. 7 and Fig. 8.

Fig. 7The kurtosis performance of the life-cycle stage indictor about the measured bearing

Fig. 8The PCA performance of the life-cycle stage indictor about the measured bearing

From Fig. 7 we can see that, the Kurtosis is not sensitive to the bearing running fault, especially the early fault, it will on duty until the end of the test bearing life, the Kurtosis showing a fluctuation. From the Fig. 8 we can see that, the machine is in normal condition during the time correlated with the first 700 points. After that time, the condition of the bearing suddenly changes. It indicates that there are some faults occurring in this bearing. Then, though the autocorrelation function method the$\mathrm{\tau }$is set as 36 for PC1 value as shown in Fig. 9.

After determine the delay time, the embedding dimension is selected by the CAO method. The results are shown in Fig. 10, the optimal embedding dimension $d$ for PC1, is chosen as 7.

Fig. 9The τ is set as 36 for PC1 value though the autocorrelation function

Fig. 10The optimal embedding dimension d is selected by the CAO method

Based on the selected optimal embedding dimension and the delay time, the LS-SVM model is used to achieve the multi-steps prediction. The essential observations (inputs) are used for forecasting the future value is set to 7, the output node number is 1.

The LS-SVM model is used to achieve the multi-steps prediction, the feature of PC1 is used (fitted into 984 points of 4 bearings). The iterated multi-steps prediction strategy is used. The PSO is used to obtain the main parameters of the model, the particle swarm population size is set to 100, the number of the particles is set to 20. The fitness function is set to get the minimum prediction error with the optimized parameters. The prediction error is set to 0.001. The PSO particle’s dimension is set to 2, the$w$is set to 0.5, the $c_1$ is set to 1, the $c_2$ is set to 1, based on the CAO method to determine the embedding dimension number, the LS-SVM model input number is set to 7. The optimized obtain parameters $\gamma$ is 90.3, the ${\sigma }^2$ is 20. Then the two parameters are used to build the LS-SVM model to train and predict the value. The training data length is set to 984 points and the test data point is set to 840 points. The degradation curves and the fitted curve (the fitting number is 5) of the 3 training bearings are showing in Fig. 11. The test beating recession curve and the fitted curve is showing in Fig. 12.

Fig. 11The degradation trend carve of three training bearing: a) the training bearing 1 and its fitted curve; b) the training bearing 1 and its fitted curve; c) the training bearing 1 and its fitted curve

Fig. 12The degradation trend carve of the test bearing

From the Fig. 11 and Fig. 12 we can see that, the recession curve and its fitted curve is similar to the recession curve and fitted curve of training bearing 2, and far away with bearing1, so the SVM prediction model should have a heavy weigh factor than the bearing1. From the Fig. 12 we can see that there are 164 points, the 140 point is chosen for predicting (The remaining life should be 24 point, one point represents the actual work hour). The prediction results are:

The weighted factors of the 3 training bearing are: $w=$ model 1: 0.0516; model 2: 0.8900; model 3: 0.0584.

The predicted remaining life is 18.3608. The prediction value is closer to the actual value 24, get a good prediction result.

In order to compare the different remaining life prediction methods that are proposed by other lecture, this research achieve a comparison among the proposed methods. 1) The vibration signals is proposed by the EMD energy entropy, and the SOM neural network is used to processed the features, the confidence value model is used to achieve the remaining life prediction [6]. 2) The time series data is processed by the phase space reconstructing method and the weighted Normalized cross correlation model is constructed to achieve the remaining life [7]. 3) The accurate bearing remaining useful life prediction based on Weibull distribution and artificial neural network [17]. 4) The bearing remaining life prediction model is achieved by the Hybrid PSO–SVM-based method [22]. 5) The proposed method in this research. The remaining life prediction results of different methods are showing in Table 3.

From the result we can see that the neural network and the SVM model are used for the remaining life prediction, however the neural network model are not stable, the result may changed even the parameters are same when constructing the model, the SVM model can work better than the ANN model, however, if only one model is used in remaining life prediction, the result will be affected by the constructed model, so the weighted models are necessary for remaining life prediction, and based on the results, we can see the proposed model work better than the proposed model in other lectures.

5.2. Case 2

In Case 1, the 4 bearing are working in the same state, they all work in the unified operating conditions. In order to further validate the proposed method, the other data provide also by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati [28] are used (8 bearing full time recession data) together with the 4 bearing data used in Case1 for validation. In this research the bearing 1 recession indictor of case 1 is used for test, and the remaining 11 bearings recession data are used for training. As can be seen from the Case 1, sine the work condition of the 3 bearing used in Case 1 are very different with the 8 bearings, so the weight factors of the 3 bearing may greater than the 8 bearings. The iterated multi-steps prediction strategy is used. The PSO is used to obtain the main parameters of the model, the particle swarm population size is set to 300, the number of the particles is set to 20. The fitness function is set to get the minimum prediction error with the optimized parameters. The prediction error is set to 0.005. The PSO particle’s dimension is set to 2, the $w$ is set to 0.1, the $c_1$ is set to 1, the $c_2$is set to 1, based on the CAO method to determine the embedding dimension number, the LS-SVM model input number is set to 7. The optimized obtain parameters $\gamma$ is 620, the ${\sigma }^2$ is 511. Then the two parameters are used to build the LS-SVM model to train and predict the value. The recession indictor (Unified interpolated to 984 points) and the fitted curve (the fit number is 5) of 11 training bearing are showing in Fig. 13. There are 984 points, the 440 point is chosen for predicting (The remaining life should be 544 point, one point represents the actual work hour).

The weight factor of every bearing prediction SVM models is: $w=$ model 1: 0.0048; model 2: 0.0043; model 3: 0.0044; model 4: 0.0034; model 5: 0.0053; model 6: 0.0056; model 7: 0.0039; model 8: 0.0041; model 9: 0.0004; model 10:0.0000; model 11: 0.9638.

The remaining life of the test bearing 1 is 559.2252. From the got result we can see that the proposed method can get the remaining life 544 of the bearing precisely, which validate the useful of the proposed method.

However, from the results of Case 1 and Case 2 we can also get that, in order to get a more accurate prediction result, the collected data must be representative, and the working condition of the test bearing must be similar with the training bearing, the prediction result can be precisely.

Fig. 13The degradation trend curves of eleven training bearing: a) the bearing degradation trend curve of bearing 1 (new added); b) the bearing degradation trend curve of bearing 2 (new added); c) the bearing degradation trend curve of bearing 3 (new added); d) the bearing degradation trend curve of bearing 4 (new added); e) the bearing degradation trend curve of bearing 5 (new added); f) the bearing degradation trend curve of bearing 6 (new added); g) the bearing degradation trend curve of bearing 7 (new added); h) the bearing degradation trend curve of bearing 8 (new added); i) the bearing degradation trend curve of bearing 1 (go in case 1); j) the bearing degradation trend curve of bearing 2 (go in case 1); k) the bearing degradation trend curve of bearing 3 (got in case 1)

Fig. 14The degradation trend carve of the test bearing