Remaining useful life prediction of the ball screw system based on weighted Mahalanobis distance and an exponential model

2.1. Bayesian theory

Most dynamic systems can be described by the following two models: a process model and an observation model:

where ${\mathbf{x}}_k$ is the actual state at time $t_k$, and ${\mathbf{y}}_k$ denotes the corresponding observation. $f(\cdot )$ and $h(\cdot )$ represent the state transition function and the observation function, respectively. ${\mathbf{n}}_k$ is the independent and identically distributed (i.i.d.) process noise, whereas ${\mathbf{v}}_k$ is the i.i.d. observation noise.

Let ${\mathbf{x}}_{0:k}=\{{\mathbf{x}}_0,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_k\}$ and ${\mathbf{y}}_{1:k}=\{{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_k\}$ denote all the available states and observations, respectively, and suppose the states follow a first order Markov process, i.e., $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{0:k-1})=p({\mathbf{x}}_k\left|{\mathbf{x}}_{k-1}\right)$. With the application of the Bayesian theory, the posterior probability density function (PDF) of the state ${\mathbf{x}}_k$ at $t_k$, namely, $p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k})$, can be calculated recursively using the following equations:

where $p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1})$ stands for the prior PDF at $t_k$, and $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ is the transition density determined by the process model. $p\left({\mathbf{y}}_k\right|{\mathbf{x}}_k)$ denotes the likelihood defined by the observation model, and $p\left({\mathbf{y}}_k\right|{\mathbf{y}}_{1:k-1})$ represents the evidence, as expressed by:

2.2. Particle filtering

Generally, it is impractical to solve Eq. (3) and Eq. (4) analytically, so the PF is utilized to approximate the posterior PDF by:

where $\{{\mathbf{x}}_k^i{\}}_{i=1}^N$ represents the particles, and $\{{\omega }_k^i{\}}_{i=1}^N$ denotes the associated weights. $N$ is the number of particles, and $\delta (\cdot )$ stands for the Dirac function.

Particles $\{{\mathbf{x}}_k^i{\}}_{i=1}^N$ are sampled from an importance PDF $q\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k})$, and their weights can be calculated by:

For easy implementation, the transition density, $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$, is often selected as the importance PDF, i.e., $q\left({\mathbf{x}}_k^i\right|{\mathbf{x}}_{k-1}^i,{\mathbf{y}}_{1:k})=p({\mathbf{x}}_k^i\left|{\mathbf{x}}_{k-1}^i\right)$, and then Eq. (7) becomes:

The procedures of the standard PF are summarized as following [34]:

1) Initialization: Set $k=$ 0, and sample particles $\{{\mathbf{x}}_0^i{\}}_{i=1}^N$ from the initial distribution $p\left({\mathbf{x}}_0\right)$. Also, initialize the weight of each particle as ${\omega }_0^i=1/N$.

2) Importance sampling: Set $k=k+$ 1, and calculate the transition density $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ with Eq. (1). Set the importance PDF $q\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k})=p({\mathbf{x}}_k\left|{\mathbf{x}}_{k-1}\right)$, and then draw particles $\{{\mathbf{x}}_k^i{\}}_{i=1}^N$ from it.

3) Weights updating: Update the weights with newly obtained measurements with Eq. (7), and then normalize the weights with the following:

4) Resampling: Remove particles with small weights and copy particles with large weights. A resampling implementation is described in detail as following [33]:

a) Set $d=$ 1: $N$, and generate a random value $u_d$ from the uniform distribution $U\left(\mathrm{0,1}\right)$.

b) Set $j=$ 1: $N$, and calculate the cumulative distribution function of weights as ${\sum }_{i=1}^j{\omega }_k^i$. If ${\sum }_{i=1}^j{\omega }_k^i\ge u_d$, duplicate ${\mathbf{x}}_k^j$ as a new particle ${\mathbf{x}}_k^d$ with the weight of $1/N$, and go back to Step a; otherwise, set $j=j+1$, and return to Step b.

5) State estimation: Estimate the current state with the resampled particles and the associated weights by:

Then, return to Step 2 and repeat Step 2–5 for the next inspection time.

3.1. HI construction

3.2. RUL prediction

Due to the complex and heterogeneous working conditions, the ball screw system usually deteriorates stochastically. Therefore, we attempt to describe the degradation path of the ball screw system with an exponential Wiener process model, as denoted by:

where $\beta$ is a random variable, and $\epsilon \left(t\right)=\sigma W\left(t\right)$ is a centered Brownian motion representing the stochastic effect.

For convenience, Eq. (20) is transformed into a linear form by using the logged observations, as represented by:

According to Eq. (21), we can obtain the following difference equation:

where ${\eta }_k=\sigma \left(W(t_k\right)-\sigma \left(W\left(t_{k-1}\right)\right)$ obeys $N(0,{\sigma }^2\mathrm{\Delta }{t}_k)$ with $\mathrm{\Delta }{t}_k=t_k-t_{k-1}$. For simplicity, let $\mathbf{\Theta }=[\beta ,\sigma ]$ represent the model parameters.

After degradation modeling, the model parameters and the health state of the system can be estimated by combining the model and the condition monitoring information. At $t_k\text{,}$ let ${\mathbf{y}}_{1:k}=\{y_1,y_2,\dots ,y_k\}$ represent the available logged measurements, and suppose the model parameters $\mathbf{\Theta }$ are known. Because the error increments ${\eta }_k$ are i.i.d. normal random variables, the conditional joint PDF of ${\mathbf{y}}_{1:k}$ can be obtained by:

In practice, the parameters $\mathbf{\Theta }$ are usually unknown. However, based on their prior PDF $p\left({\mathrm{\Theta }}_0\right)$ and Eq. (23), we can determine their joint posterior PDF $p\left({\mathbf{\Theta }}_k\right|{\mathbf{y}}_{1:k})$ at $t_k$ through the Bayesian theory. After the health state estimation and parameters update, the PDF of the RUL at $t_k$ can be predicted.

The definition of the RUL at $t_k$ can be expressed as:

where $\mathrm{i}\mathrm{n}\mathrm{f}\{\cdot \}$ represents the lower bound of a variable, $\gamma$ is a predefined failure threshold, and $z(t_k+l_k)$ is the predicted health state at $t_k+l_k$, which can be calculated by:

According to the independent increment property of the Brownian motion, the PDF of the RUL at $t_k$ can be represented as:

In this paper, the PF is utilized to combine the HI and the degradation model for updating parameters, estimating the health state, and predicting the RUL. To update the parameters along with the health status, parameters $\beta$ and $\sigma$ are regarded as latent states of the system as well as the health state $z$. Then the state space model of the system become ${\mathbf{x}}_k=[z_k,{\beta }_k,{\sigma }_k]$, and the RUL prediction process with the PF includes the following steps:

1) At initial time $t_0$, sample particles $\{{\mathbf{x}}_0^i{\}}_{i=1}^N$ from the initial distribution $p\left({\mathbf{x}}_0\right)$ and initialize the particle weights by ${\omega }_0^i=1/N$.

2) Once a new measurement $y_k$ is obtained at $t_k(k\ge 1)$, sample particles $\{{\mathbf{x}}_k^i{\}}_{i=1}^N$ according to the importance PDF defined by the following state transition functions:

3) Update the particle weights ${\omega }_k^i$ by:

4) Resample particles $\{{\mathbf{x}}_k^d{\}}_{i=1}^N$ from $\{{\mathbf{x}}_k^i{\}}_{i=1}^N$ according to $\{{\omega }_k^i{\}}_{i=1}^N$, using the resampling algorithm described in Section 2.2.

5) Update the health state and the model parameters at $t_k$ with Eq. (10).

6) With the estimated HI, ${\widehat{z}}_k$, and the model parameters, ${\widehat{\beta }}_k$ and ${\widehat{\sigma }}_k$, predict the RUL by:

7) Set $k=k+1$, return to Step 2 and repeat Step 2-6 until ${\widehat{z}}_k>\gamma$.

4.1. Experimental system and vibration data

Since the ball screw is a highly reliable product, it usually takes a long time for the ball screw system to fail. Consequently, it is difficult to acquire the condition monitoring data during the entire working life of the ball screw system. To make up for this deficiency, we designed and manufactured an experimental system to conduct accelerated life tests for the ball screw system. With the experimental platform, the degradation speed of the ball screw system is expedited by increasing the system load, which causes the time amount for the system to become invalid to be reduced. Thus, adequate degradation data can be collected to study ball screw system diagnostics and prognostics. Fig. 3 shows the overview of the experimental system, and it consists of two parts: the testbed for the accelerated life test of the ball screw system, and the control cabinet.

The testbed depicted in Fig. 4 can simulate the actual movement of the ball screw system and exert the load on the ball screw. To collect the condition monitoring data, two uniaxial accelerometers are fixed on the bearings at the drive end and the floating end, respectively, and a triaxial accelerometer is attached to the ball nut, as illustrated in Fig. 5. The main dimensions of the ball screw are listed in Table 2. During the testing process, the rotating speed was 1,000 rpm, and the screw rotation was controlled by the servo motor. The ball screw was subjected to a load of 100 N axially throughout the test, which was produced by the magnetic powder brake and transmitted by the gear and rack.

Fig. 3The experimental system for accelerated life tests of the ball screw system

Fig. 4The testbed for accelerated life tests of the ball screw system

Fig. 5Installation positions of accelerometers

In the test, vibration signals were collected every 30 min, and the sampling frequency was 5,000 Hz. We used a digital microscope to observe the wear status of the ball screw and take pictures, as described in Fig. 6. With the microscope, we can obtain the enlarged view of the ball screw, as illustrated in Fig. 6(b). By increasing the magnification times, we can see more details about the condition of the ball screw, as shown in Fig. 6(c) and Fig. 6(d), which display the photographs of the ball groove before and after the test, respectively. It is seen that pittings appeared on the surface of the screw at the end of the test. Vibration signals acquired from uniaxial accelerometer 1 throughout the test are shown in Fig. 7. In the beginning, the vibration data remained relatively stable, and then small increases occurred. In the late stage of the degradation, the signals increased rapidly until the ball screw system failed.

Fig. 6a) The screw, b) the enlarged view of the screw, c) the enlarged view of the ball groove before the test, d) the enlarged view of the ball groove after the test

Fig. 7Vibration signals of the ball screw system

4.2. HI construction of the ball screw system

To construct the HI for the ball screw system, 40 original features, as listed in Table 1, were extracted from the raw vibration data. Then we calculated the trendability of each feature, as depicted in Fig. 8. To select features that appropriately correlated to the defect development and propagation, features whose values of trendability were 0.5 or below were removed, and then 22 features were obtained. To identify the redundant information among these features, they were clustered with the correlation clustering algorithm; the PBM-index values with different cluster numbers are shown in Fig. 9. The index achieved its maximum value when $H=$ 3, so these 22 features were divided into three classes. After clustering, the feature with the highest trendability in each class was chosen as the representative feature for its cluster. Accordingly, we obtained three typical features for feature fusion. To compare the representative features, they were normalized and displayed in Fig. 10. It is observed that different typical features represent distinct deterioration trends of the ball screw system. Specifically, Y5 contains abundant information about the early stage of the degradation, and Y12 describes the fluctuation in the middle of the damage propagation. Furthermore, Y32 is sensitive to the late stage of the degradation.

Fig. 8Trendability of original features

Fig. 9PBM-index values with different cluster numbers

To take advantage of each typical feature, the MD was employed to fuse the information represented by the features. Before calculating the MD, the representative features were weighted according to their values of trendability. Then we could obtain the WDMD at any measurement time using Eq. (18), as depicted in Fig. 11. To verify the benefits of the WDMD, the MD of 40 original features were calculated for comparison, as shown in Fig. 12. It is seen that the WDMD indicates a more evident deterioration tendency during the entire lifecycle. To compare these two HIs quantitatively, the trendability values of the WDMD and the MD were calculated, as illustrated in Table 3, and it is found that the trendability value of the WDMD is higher than that of the MD. Therefore, the WDMD is more sensitive to the degradation development of the ball screw system, and it is more suitable for RUL prediction.

Fig. 10Representative features

Fig. 11The WDMD of the ball screw system

Fig. 12The MD of the ball screw system

4.3. RUL prediction of the ball screw system

In this section, the PF was used to integrate the WDMD with the presented exponential model for model parameters update, health state estimation, and RUL prediction. The number of particles was set to be 1,000. Fig. 13 shows the estimation results of the WDMD during the whole degradation process, and it is seen that the proposed model can track the deterioration path well.

Fig. 13WDMD estimation results.

The update process of the model parameters $\beta$ and $\sigma$ is displayed in Fig. 14. As more measurements became available, $\sigma$ converged to a stable value after 12 samples with fast speed. In contrast, $\beta$ experienced fluctuations during the estimation process. At each inspection time index, the PDF of the RUL can be predicted with the estimated health state and model parameters.

Fig. 14The update process of model parameters: a) β, b) σ

a)

b)

To manifest the merits of the presented model, both the linear Wiener process model [40] and the nonlinear Wiener process model [28] were employed for comparison. The means of the predicted RUL by three models at different measurement intervals are shown in Fig. 15. The proposed model converged to the real RUL after 83 points, whereas the other two models differed greatly from the actual RUL value at this stage. Among these three models, the linear model generated the largest prediction errors, because the degradation of the ball screw system did not increase linearly. The RUL estimation results of the nonlinear model were more accurate than those of the linear model, but it also could not estimate the RUL accurately until the end of the deterioration stage. For the presented exponential model, it forecasted the RUL with large errors in the beginning since there were not enough WDMD values. As more measurements became available, the prediction errors were reduced, and the predicted value converged to the real RUL gradually.

Fig. 15The means of the predicted RUL by three models at different inspection time indexes

To compare the models mentioned above quantitatively, we calculated the prediction errors at $t_k$ by:

where $L_k^{\mathrm{*}}$ is the real value of the RUL, and $L_k$ denotes the predicted RUL.

The prediction errors of the three models at all time indexes were calculated using Eq. (30). For comparison, we divided the results of each model into three intervals and calculated the percentage of points in each interval, as shown in Table 4. It is observed that the proposed model predicts the RUL for 52.08 % of points with a prediction error of no higher than 30 %, whereas the nonlinear model’s percentage of points with prediction errors of 50 % or lower is only 30.56 %. The linear model estimates only 4.86 % of points with prediction errors of 50 % or lower. Therefore, the proposed model can predict the RUL of the ball screw system with the lowest prediction error.

The superior properties of the presented model can be explained as follows. The ball screw system deterioration has an exponentially growing trend comprehensively, so that the exponential model can describe the degradation process best. In addition, the Brownian motion parameter is considered and updated with the PF in our method, which can indicate the Brownian error in real time. Hence, the proposed model can improve the accuracy of ball screw system RUL estimation.