# Remaining Useful Life Prediction of High-Frequency Swing Self-Lubricating Liner

Authors: Xiuhong Hao; Shuqiang Wang; Mengfan Chen; Deng Pan
Journal: Shock and Vibration
 (2021)

---

## Abstract

The remaining useful life (RUL) prediction of self-lubricating spherical plain bearings is essential for replacement decision-making and the reliability of high-end equipment. The high-frequency swing self-lubricating liner (HSLL) is the key component of self-lubricating spherical plain bearings under high-frequency oscillation conditions. In this study, a RUL prediction method was proposed based on the Wiener process and grey system theory. First, the predictive processing of the wear depth was carried out using the grey model GM(1,1) to reduce the randomness and enhance the inherent regularity of the life test data. A degradation process model was established and the RUL was predicted online with the model parameter estimates based on the Bayesian updating strategy. Finally, examples were provided to elaborate the RUL prediction of the HSLL. The results show that the prediction accuracy of the proposed RUL prediction model is higher than that of the simple Wiener process during the entire residual life cycle of the HSLL. Based on the original wear data, the prediction accuracy of the RUL exhibited a strong dependence on prior samples and was relatively low owing to the larger deviation of the wear rate between the test sample and prior samples.


---

## Body

## 1. Introduction

As key activity connectors that are applied extensively to aerospace, engineering machinery, and water projects, self-lubricating spherical plain bearings (SSPB) offer numerous advantages, such as a compact structure, no need for supplementary lubricant, and long life. Moreover, it directly affects the safety and reliability of high-end equipment [1, 2]. A woven self-lubricating liner is critical for ensuring the service performance index of SSPB. The remaining useful life (RUL) prediction of the self-lubricating liner serves as the core and foundation of the service life assessment of SSPB, as well as the fault prediction and health management of high-end equipment [3–5].With the rapid development of machine learning theory, RUL predictions based on support vector machines have been developed [6, 7]. The performance degradation parameters of the mechanical parts are fitted and the degradation laws are obtained based on the artificial intelligence algorithm. The RUL is calculated with the failure threshold. However, the probability distribution function (PDF), which can embody the random uncertainty of the RUL, cannot be obtained. That is, it is impossible to achieve the probabilistic prediction of the mechanical parts and the evaluation of the RUL prediction rapidly for batch products [8].As a superior degradation modeling of the random process, RUL predictions based on the Wiener process have been used extensively in the key components of aerospace, railway, and electronic applications. Deng et al. carried out a turbofan engine degradation simulation and improved the near-failure prediction accuracy using the surrogate Wiener propagation model and a long short-term memory network [9]. Guan et al. presented a combined model with a Wiener process-based wear model and unequal interval weighted grey linear regression. The RUL of the Guangzhou Metro was analyzed with high precision [10]. Pan et al. proposed a RUL prediction method based on the Wiener degradation model by considering temporal uncertainty, measurement uncertainty, and unit-to-unit heterogeneity. The effectiveness of the Wiener degradation model was verified by the degradation dataset of a light-emitting diode [11]. Kong et al. constructed a random process model by inserting the jump at the change point in the degradation process, which was described by the linear Wiener process. The degradation processes of the bearings were simulated [12].Because numerous interference factors exist during life testing or equipment running, excessive randomness appears in the degradation data and the accuracy of the RUL prediction is significantly reduced. The grey data can be whitened by the accumulated generating operation (AGO) based on the grey system theory [13]. The randomness of the original data is subsequently significantly reduced, and the inherent regularity is excavated. The grey system theory has been used extensively in areas such as the RUL prediction of mechanical parts and trend predictions of social phenomenon development [14, 15]. Sun et al. predicted the nitrous oxide emissions for 2030 in six countries based on three advanced mathematical grey prediction models [14]. Huang et al. proposed a grey online modeling surface roughness monitoring system to predict the surface roughness in end milling operations [16]. Li et al. constructed a grey model for concrete acidification prediction; they analyzed the effects of the pH value, concrete cover thickness, and surface coating on the service life [17]. Ene and Öztürk established a forecasting system for discarded end-of-life vehicles based on the grey system theory [18]. The prediction accuracy was improved by parameter optimization, Fourier series, and Markov chain correction.In this study, life tests were carried out on a high-frequency swinging self-lubricating liner (HSLL). A RUL model was proposed based on the grey system theory and Wiener process. To improve the prediction accuracy, the original wear depth of the HSLL was conducted, and the prediction data of the wear depth were obtained. Thereafter, the degradation process model of the wear depth was established and the RUL was predicted at any moment with the model parameter estimates based on the Bayesian updating strategy. Examples were provided to illustrate the effectiveness of the RUL model.This paper introduced grey system theory into the simple Wiener process. The inherent regularity of the original wear data was clearly reflected through predictive processing based on the GM(1,1) model. Compared to the simple Wiener process, the prediction accuracy of the proposed method for the RUL was improved significantly in the entire residual life cycle. And the uncertainty of the model in this paper is smaller than the simple Wiener process.


## 2. Life Test of HSLL

The HSLL used in the life test was developed by Yanshan University. Brown polytetrafluoroethylene fibers and Kevlar49 fibers were woven into a satin structure, impregnated, semicured, and cured using a modified phenolic resin to obtain the test material. The HSLL had a maximum swing frequency of 15 Hz, a maximum load of 2.8 kN, and a maximum swing angle of ±10°. According to SAE AS 81819A, the maximum wear depth of the HSLL after 25,000 oscillations should be ≤ 0.127 mm [19].The test samples were in the form of a cylinder contact, as illustrated in Figure1. Two half-rings and a test shaft made of no. 45 steel with chrome plating on the surface were placed in a group. The failure threshold of the HSLL was Xf = 0.254 mm. The life tester was developed by Yanshan University and is shown in Figure 2. The maximum swing frequency of the tester was 30 Hz, the maximum load was 80 kN, and the maximum swing angle was ±15°. The operating parameters for the life test were the most stringent conditions of the self-lubricating liner; that is, f = 15 Hz, F = 2.8 kN, and θ = ±10°. The wear depth data was measured through the displacement sensor and transmitted to a computer. The online measurement and collection of the data can be realized. The life tester had the function of self-compensation of the wear depth; that is, it could compensate for the measurement error caused by the change in the environmental temperature and deformation of the parts under loading.Figure 1
Friction pair for life test.Figure 2
Testing machine and self-lubricating fabric liner under test. (a) Testing machine. (b) Self-lubricating fabric liner under test.
(a)(b)Based on the HSLL, all samples were prepared under the same bonding process and were tested under the same conditions using the same tester. The wear depth, which is the most critical characteristic parameter, was measured and evaluated as the degradation parameter. Considering that the chrome face on no. 45 steel has a higher degree of hardness and excellent wear resistance, the wear depth between the friction pairs was mainly located on the liner. Therefore, the measured wear depth was taken as the HSLL wear depth. The life tester had an automatic alarm function for failure and excessive wear. The experimenter periodically judged the running state of the life tester and the increasing tendency of the wear depth to avoid the occurrence of abnormal test results. Life tests were conducted on the six groups of samples and the curves of the wear depth were measured, as illustrated in Figure3.Figure 3
Wear curves of the self-lubricating liner.


## 3. Original Data Processing Based on Grey System Theory

In HSLL life tests, vibration and other random factors still existed, and random errors of the wear depth were inevitable, although displacement compensation with the environmental temperature and loading were considered. These will have a certain influence on the prediction accuracy of the RUL. Grey systems with several uncertain factors are the main objects of studies in the grey system theory. The correct description and effective monitoring of the system operation behavior and evolution law can be realized when partially known information is generated and developed; for example, by using the AGO. Once the test data of the wear depth are processed based on the grey system theory, the regularity of the wear depth can be enhanced, and the prediction accuracy of the RUL can be improved.


### 3.1. Grey Model GM(1,1)

The grey model GM(1,1) is the core model of the grey system theory, and the corresponding differential equations can be obtained by generating the discrete grey data into the regular white data [20, 21]. The original series of wear depths for the HSLL can be defined as (1)X0=x01,x02,…,x0n,where the superscript (0) of X(0) is the original series of the wear depth.The first-order AGO can be expressed asX(1) and is generated from X(0): (2)X1=x11,x12,…,x1n,where x1k=∑i=1kx0i,k=1,2,…,n.The background value of GM(1,1) is obtained fromX(1):(3)Z1=z12,z13,…,z1n,where z1k=μx1k+x1k−1,k=2,3,…,n, and μ is generally 0.5.Thereafter, the grey differential equation of GM(1,1) can be obtained as(4)x0k+az1k=b,where a and b are the developing coefficient and grey controlled variable, respectively.The grey derivative ofX(1) is defined in backward difference quotient [22]:(5)dx1kdk=Δx1kΔk=x0k=x1k−x1k−1,k=2,3,…,n.For the GM(1,1) presented in equation (4), if the time k = 2, 3, …, n, of x(0) (k) is regarded as a continuous variable t, X(1) can be regarded as a function of t, and it can be expressed as x(1) = x(1)(t). Meanwhile, let x(0)(k) correspond to dx(1)/dt, and let the background value z(1)(k) correspond to x(1)(t). Then, the whitening differential equation of GM(1,1) can be obtained as follows:(6)dx1tdt+ax1t=b.


### 3.2. Original Data Processing

First, inspection is necessary for the original wear data to ensure the feasibility of the modeling method. Taking sample 1 as an example, the class ratio of the original seriesX(0) can be calculated as follows:(7)λk=x0k−1x0k,k=2,3,…,n.If all class ratiosλ(k) are in the capacitable coverage range Θ=e−2/n+1,e2/n+1, the original series X(0) can be solved by establishing GM(1,1). Otherwise, it is necessary to transform the original series X(0) to make it fall into the capacitable coverage range, that is, by taking the appropriate constant c for the translation transformation. The transformation can be defined as(8)y0k=x0k+c,k=1,2,…,n.The minimum and maximum values can be calculated, which are 0.9841 and 1.0173, respectively, by substituting the original seriesX(0) of sample 1 into equation (7). Because the capacitable coverage range is Θ = (0.9744, 1.0260), the translation transformation need not be carried out.Moreover,a and b in GM(1,1) can be estimated by the least-squares method and the grey differential equation of the original series X(0):(9)ab=DTD−1DTY,where(10)D=−z121−z131⋮⋮−z1n1,Y=x02x03⋮x0n.By substituting the original seriesX(0) of sample 1 into equations (9) and (10), a and b can be calculated as 5.17 × 10−4 and 98.31, respectively. The first-order AGO can be obtained by substituting a and b into equation (6):(11)x^1k+1=x01−bae−ak+ba,k=0,2,…,n−1.Then, the predicted value of the original degradation series can be obtained as follows:(12)x^0k+1=x^1k+1−x^1k,k=0,2,…,n−1.Finally, the average relative residual can be determined as follows:(13)ε¯=∑k=2n1n−1x0k−x^0kx0k.If the average relative residualε¯ is less than 0.2, the fitting precision of the grey model GM(1,1) to the original series X(0) meets the general requirements. If ε¯ is less than 0.1, the GM(1,1) model exhibits excellent fitting precision on the original series X(0).The predicted values and corresponding average relative residuals can be obtained by substituting the wear data of samples 1, 2, and 3 into the above equations. In this case,ε¯ of samples 1, 2, and 3 are 0.0577, 0.0377, and 0.1402, respectively. Therefore, the GM(1,1) models exhibit high fitting precision, and the prediction curves of samples 1, 2, and 3 are illustrated in Figure 4. The figure indicates that the influence of the external disturbance factors was significantly weakened and the inherent laws of the degenerate variables were extracted effectively.Figure 4
Prediction results of samples 1, 2, and 3.
(a)(b)(c)


## 3.1. Grey Model GM(1,1)

The grey model GM(1,1) is the core model of the grey system theory, and the corresponding differential equations can be obtained by generating the discrete grey data into the regular white data [20, 21]. The original series of wear depths for the HSLL can be defined as (1)X0=x01,x02,…,x0n,where the superscript (0) of X(0) is the original series of the wear depth.The first-order AGO can be expressed asX(1) and is generated from X(0): (2)X1=x11,x12,…,x1n,where x1k=∑i=1kx0i,k=1,2,…,n.The background value of GM(1,1) is obtained fromX(1):(3)Z1=z12,z13,…,z1n,where z1k=μx1k+x1k−1,k=2,3,…,n, and μ is generally 0.5.Thereafter, the grey differential equation of GM(1,1) can be obtained as(4)x0k+az1k=b,where a and b are the developing coefficient and grey controlled variable, respectively.The grey derivative ofX(1) is defined in backward difference quotient [22]:(5)dx1kdk=Δx1kΔk=x0k=x1k−x1k−1,k=2,3,…,n.For the GM(1,1) presented in equation (4), if the time k = 2, 3, …, n, of x(0) (k) is regarded as a continuous variable t, X(1) can be regarded as a function of t, and it can be expressed as x(1) = x(1)(t). Meanwhile, let x(0)(k) correspond to dx(1)/dt, and let the background value z(1)(k) correspond to x(1)(t). Then, the whitening differential equation of GM(1,1) can be obtained as follows:(6)dx1tdt+ax1t=b.


## 3.2. Original Data Processing

First, inspection is necessary for the original wear data to ensure the feasibility of the modeling method. Taking sample 1 as an example, the class ratio of the original seriesX(0) can be calculated as follows:(7)λk=x0k−1x0k,k=2,3,…,n.If all class ratiosλ(k) are in the capacitable coverage range Θ=e−2/n+1,e2/n+1, the original series X(0) can be solved by establishing GM(1,1). Otherwise, it is necessary to transform the original series X(0) to make it fall into the capacitable coverage range, that is, by taking the appropriate constant c for the translation transformation. The transformation can be defined as(8)y0k=x0k+c,k=1,2,…,n.The minimum and maximum values can be calculated, which are 0.9841 and 1.0173, respectively, by substituting the original seriesX(0) of sample 1 into equation (7). Because the capacitable coverage range is Θ = (0.9744, 1.0260), the translation transformation need not be carried out.Moreover,a and b in GM(1,1) can be estimated by the least-squares method and the grey differential equation of the original series X(0):(9)ab=DTD−1DTY,where(10)D=−z121−z131⋮⋮−z1n1,Y=x02x03⋮x0n.By substituting the original seriesX(0) of sample 1 into equations (9) and (10), a and b can be calculated as 5.17 × 10−4 and 98.31, respectively. The first-order AGO can be obtained by substituting a and b into equation (6):(11)x^1k+1=x01−bae−ak+ba,k=0,2,…,n−1.Then, the predicted value of the original degradation series can be obtained as follows:(12)x^0k+1=x^1k+1−x^1k,k=0,2,…,n−1.Finally, the average relative residual can be determined as follows:(13)ε¯=∑k=2n1n−1x0k−x^0kx0k.If the average relative residualε¯ is less than 0.2, the fitting precision of the grey model GM(1,1) to the original series X(0) meets the general requirements. If ε¯ is less than 0.1, the GM(1,1) model exhibits excellent fitting precision on the original series X(0).The predicted values and corresponding average relative residuals can be obtained by substituting the wear data of samples 1, 2, and 3 into the above equations. In this case,ε¯ of samples 1, 2, and 3 are 0.0577, 0.0377, and 0.1402, respectively. Therefore, the GM(1,1) models exhibit high fitting precision, and the prediction curves of samples 1, 2, and 3 are illustrated in Figure 4. The figure indicates that the influence of the external disturbance factors was significantly weakened and the inherent laws of the degenerate variables were extracted effectively.Figure 4
Prediction results of samples 1, 2, and 3.
(a)(b)(c)


## 4. RUL Prediction Based on Wiener Process

### 4.1. Degradation Process Model of HSLL

#### 4.1.1. Model of Wear Process

The degradation process model of the wear depth for the HSLL could be constructed using the Wiener process. The wear depth at timet is defined as X(t) and is expressed as(14)Xt=X0+αt+σBBt,where X(0) is the wear depth at the initial time and X(0) is generally equal to 0. However, for the HSLL, X(0) is not equal to 0 and data conversion must be conducted by X(t)-X(0). Moreover, α is the drift coefficient and follows a normal distribution, that is, α∼Nμα,σα2 [23, 24]; σB is the diffusion coefficient, which can describe the inconsistency and instability of the HSLL life test as well as the effects of certain random factors such as vibration; B(t) is the standard Brownian motion and follows a normal distribution, that is, B(t)∼N(0, t); α is the random parameter; and σB is the common parameter.During the online RUL prediction of the single sample, the parametric estimated value must be accurately obtained, and it is difficult to meet the requirements depending only on the wear depth of the selected sample. The distributed hyperparameter ofα and the estimated value of σB should be obtained from the wear data of historical prior samples. Then, the random parameter α of the selected sample can be realized by online updating using the Bayesian method.


#### 4.1.2. Parameter Estimation of Wear Process Model

(1) Prior Parameter Estimation. The parameter estimation of the random parameter α and common parameter σB in the wear process model were executed with samples 1, 2, and 3 considered as historical prior samples. The sampling time and sampling number were the same for all samples. The sampling numbers were defined as m. The wear depth of the historical prior samples was fitted by the GM(1, 1) model, and the corresponding predicted values are given as follows: (15)X=X1,1X1,2⋯X1,mX2,1X2,2⋯X2,mX3,1X3,2⋯X3,m,where Xi,j is the prediction value of sample i at time ti,j, 1 ≤ j ≤ m.Following the grey data processing, the degradation process model of the wear depth was established:(16)Xi,j=α0t+σBBti.j,where α0 is the prior value of the random parameter α and follows a normal distribution; that is, π0α0∼Nμα,0,σα,02. The hyperparameters μα,0 and σα,02 could be calculated from historical prior samples.The wear depth of samplei, Xi = (Xi,1, Xi,2,…, Xi,m)T, follows a multivariate Gaussian distribution. The mean and covariance can be obtained as follows: (17)μ=μα,0t,∑=σα,02ttT+Ω,where t = (t1, t2,…, tm)T, Ω=σB2Q, and(18)Q=t1t1⋯t1t1t2⋯t2⋮⋮⋱⋮t1t2⋯tm.The log-likelihood function can be generated based on the wear data of the historical prior samples:(19)Lθ|X1:3=−3m2ln2π−32ln∑−12∑i=13Xi−μα,0T∑−1Xi−μα,0,where θ is the estimated parameter, θ = (μα, 0, σα, 0, σB), and X1:3 denotes the wear data of the historical prior samples.To simplify the terms in the log-likelihood, we use the results(20)∑=Ω1+σα,02tTΩ−1t,∑−1=Ω−1−σα,021+σα,02tTΩ−1tΩ−1ttTΩ−1.The first partial derivatives ofL(θ|X1:3) for μα,0 and σα,0 can be expressed as(21)∂∂μα,0Lθ|X1:3=∑i=13tTΩ−1Xi−3μα,0tTΩ−1t1+σα,02tTΩ−1t,(22)∂∂σα,0Lθ|X1:3=−3σα,0tTΩ−1t1+σα,02tTΩ−1t+σα,0∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t1+σα,02tTΩ−1t2.The expressions ofμα,0 and σα,0 can be obtained using equations (21) and (22) equal to 0(23)μ^α,0=∑i=13tTΩ−1Xi3tTΩ−1t,(24)σ^α,0=∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t3tTΩ−1t2−1tTΩ−1t1/2.By substituting equations (23) and (24) into equation (19), the contour likelihood function of σB can be given by [25](25)LσB|X1:3=−3m2ln2π−32−32lnΩ−12∑i=13XiTΩ−1Xi−∑i=13tTΩ−1Xi2tTΩ−1t−32ln∑i=13tTΩ−1Xi23tTΩ−1t−∑i=13tTΩ−1X232tTΩ−1t.The maximum of equation (25) can be calculated by the optimization toolbox of MATLAB, and the maximum likelihood estimate of the common parameter σB can be determined as follows:(26)σ^B=5.7095×10−5.By substituting equation (26) into equations (23) and (24), the estimated values of μα,0 and σα,0 can be obtained as follows:(27)μ^α,0=4.9750×10−5,σ^α,0=4.1708×10−6.(2) Real-Time Parameter Updating. The wear data of samples 4, 5, and 6 were selected as the test sets and the random parameter α was updated in real time. Following the grey data processing, the wear process model can be expressed as (28)Xk=αktk+σBBtk,where Xk is the wear depth at time tk of the single sample online and αk is the posterior value of the random parameter α and follows a normal distribution; that is, πkαk∼Nμα,k,σα,k2.It was assumed that the wear depth beforetk is X1:k = (X1, X2,…, Xk)T of the single sample online, and the time interval is Δt. Based on the Bayesian method, the posterior distribution of the random parameter α can be determined by [26, 27](29)παk|X1:k∝LX1:k|αkπ0α0,where the prior information follows a normal distribution, that is, π0α0∼Nμα,0,σα,02, and L(X1:k|αk) is a likelihood function with given random parameters.Based on equation (14) and the basic nature of the Brownian motion, L(X1:k|αk) can be obtained as follows:(30)LX1:k|αk=exp−∑q=1kXq−Xq−1−αkΔt2/2σB2Δt∏q=1k2πσB2Δt.Because bothL(X1:k|αk) and π0(α0) follow a normal distribution, π(αk|X1:k) also follows a normal distribution. Based on equations (29) and (30), the mean and variance of π(αk|X1:k) can be obtained as follows, respectively:(31)μα,k=μα,0/σα,02+∑q=1kXq−Xq−1/σB2∑q=1kΔt2/σB2Δt+1/σα,02,σα,k=1∑q=1kΔt2/σB2Δt+1/σα,02.


### 4.2. RUL Prediction

The RUL distribution could be obtained using the concept of the first arrival times during the wear process. When the random process {X(t), t ≥ 0} reached the failure threshold, the time was defined as the sample life:(32)T=inft:Xt≥Xf|X0<Xf,where Xf is the failure threshold, X(0) is the wear depth of the initial moment, and X(0) = 0.The residual useful lifeLk in the current runtime tk can be defined as(33)Lk=lk:T−tk|T>tk,where lk is the online RUL value.Because the RUL of the Wiener process with linear drift follows an inverse Gaussian distribution [28], the PDF of Lk can be obtained as follows:(34)fLklk=Xf−Xk2πlk3σα,k2lk+σB2exp−Xf−Xk−μα,klk22lkσα,k2lk+σB2.Based on the observed data, the PDF of the RUL for sample 4 could be calculated and is illustrated in Figure5.Figure 5
PDF of RUL.Figure5 indicates that the PDF curves become narrower with wear data information increases, which indicates the uncertainty of RUL prediction is increasingly small, and the results of RUL prediction will be increasingly accurate.The RUL of a single sample can be written as(35)ELk=EELk|μα,k,σα,k=EXf−Xkμα,k,σα,k=Xf−Xkσα,k2exp−μα,k22σα,k2∫0μα,kexpx22σα,k2dx=2Xf−Xkσα,kDμα,k2σα,k,where D(z) is the Dawson function with respect to z:(36)Dz=exp−z2∫0zexpx2dx,z=μα,k2σα,k.Whenz is sufficiently large, according to the approximate property of the integral, D(z) ≈ 1/(2z). The RUL can be obtained as(37)ELk≈Xf−Xkμα,k.600 min of the wear data was selected as the initial prediction point. The random parameters were updated every 300 min. The RUL of samples 4, 5, and 6 could be predicted with the origin wear depth and grey processing data, respectively. The curves of the actual RUL and the prediction RUL were depicted in Figure6.Figure 6
RUL and confidence interval for HSLL based on original wear depth and grey processing data. (a) Origin wear depth. (b) Grey processing data. (c) Confidence interval of samples 4. (d) Confidence interval of samples 5. (e) Confidence interval of samples 6.
(a)(b)(c)(d)(e)Figure6(a) indicates that the predicted values of the three samples were close at the initial prediction point. The actual and predicted values of sample 5 exhibited better consistency. For samples 4 and 6, which had lower/higher wear rates, the prediction deviation was relatively large. However, the predicted life values gradually became closer to the actual life values and were almost equal to the actual life values at the end of the wear process. Following the wear depth being processed by the GM(1,1) model, as shown in Figure 6(b), the initial prediction point of predicted values were also basically close, but the prediction accuracy of the RUL was improved significantly in the entire residual life cycle.Simultaneously, Figures6(c)–6(e) display the 95% pointwise confidence intervals for the RUL based on these prediction points. Basically, all confidence intervals can cover the actual values of the RUL. The actual RUL line of sample 5 was positioned centrally over the confidence interval, for the origin wear depth, and the confidence intervals of samples 4 and 6 were relatively biased. And compared to the simple Wiener process, the 95% confidence interval of the proposed method in the paper is narrower, which can help reduce uncertainty in residual life prediction results and can better adapt to the individual needs of residual life prediction. These further indicate that the errors of the proposed method are lower.Meanwhile, the RULs of samples 2, 4, and 6 when samples 1, 3, and 5 were selected as the prior samples are shown in Figure7.Figure 7
RUL and confidence interval for HSLL based on original wear depth and grey processing data. (a) Origin wear depth. (b) Grey processing data. (c) Confidence interval of samples 2. (d) Confidence interval of samples 4. (e) Confidence interval of samples 6.
(a)(b)(c)(d)(e)Figure7 had the RUL prediction curves and the 95% pointwise confidence intervals similar to Figure 6. When the differences of the wear rates for samples 2, 4, and 6 were larger, the prediction deviations were relatively large at the initial phase. But the predicted life values gradually became closer to the actual life values with the wear process. Similarly, the overall prediction accuracy of the proposed method is higher than the original wear data.Synthesizing Figures5–7, the PDF of RUL became narrower with the continuous accumulation of data. That is, the uncertainty of RUL prediction becomes smaller and smaller. And the relative error between the prediction of the RUL and the actual RUL is smaller. These happened because the results of model parameter estimation were gradually closer to the reality of sample degradation with the test data increasing.To analyze the prediction precision of the proposed RUL method, multiple evaluation indicators, namely, the root mean square error (RMSE), mean absolute error (MAE), and symmetric mean absolute percentage error (SMAPE), were calculated and can be given as follows:(38)RMSE=1m∑i=1mhti−yi2,(39)MAE=1m∑i=1mhti−yi,(40)SMAPE=1m∑i=1mhti−yihti+yi/2×100%,where m is a positive integer, h(ti) is the predicted value at point i, and yi is the actual RUL at point i.The evaluation indicators of samples 4, 5, and 6 can be calculated using equations (38)–(40), as shown in Table 1.Table 1
The RMSE, MAE, and SMAPE of the predicted values.
Indicator456Original wear dataGrey processing dataOriginal wear dataGrey processing dataOriginal wear dataGrey processing dataRMSE0.62930.25160.20540.23310.93600.4874MAE0.54980.15810.18460.18860.81920.3597SMAPE34.60%9.30%9.60%8.64%26.59%13.33%Table1 shows that the assessment results were similar for RMSE and MAE. For sample 5, the prediction accuracy based on the original wear data was slightly higher than that based on the grey processing data. But for samples 4 and 6, the prediction accuracy of the proposed method was much higher. Based on the SMAPE, the predicted trends of the proposed method were all more accurate.The broken line graph of RMSE, MAE, and SMAPE shown in Figure8 reflects the changes in the RUL prediction accuracy of the predicted values at point i. All three indicators decreased gradually, and the prediction accuracy all gradually increased. Because of the great difference in the wear rates among samples, the values of the three indicators based on the original wear data varied greatly. However, all three indicators based on the grey processing data remained at relatively low levels, and the effectiveness of the proposed method in this paper was demonstrated.Figure 8
RMSE MAE and SMAPE of predicted RUL values. (a) RMSE. (b) MAE. (c) SMAPE.
(a)(b)(c)


## 4.1. Degradation Process Model of HSLL

### 4.1.1. Model of Wear Process

The degradation process model of the wear depth for the HSLL could be constructed using the Wiener process. The wear depth at timet is defined as X(t) and is expressed as(14)Xt=X0+αt+σBBt,where X(0) is the wear depth at the initial time and X(0) is generally equal to 0. However, for the HSLL, X(0) is not equal to 0 and data conversion must be conducted by X(t)-X(0). Moreover, α is the drift coefficient and follows a normal distribution, that is, α∼Nμα,σα2 [23, 24]; σB is the diffusion coefficient, which can describe the inconsistency and instability of the HSLL life test as well as the effects of certain random factors such as vibration; B(t) is the standard Brownian motion and follows a normal distribution, that is, B(t)∼N(0, t); α is the random parameter; and σB is the common parameter.During the online RUL prediction of the single sample, the parametric estimated value must be accurately obtained, and it is difficult to meet the requirements depending only on the wear depth of the selected sample. The distributed hyperparameter ofα and the estimated value of σB should be obtained from the wear data of historical prior samples. Then, the random parameter α of the selected sample can be realized by online updating using the Bayesian method.


### 4.1.2. Parameter Estimation of Wear Process Model

(1) Prior Parameter Estimation. The parameter estimation of the random parameter α and common parameter σB in the wear process model were executed with samples 1, 2, and 3 considered as historical prior samples. The sampling time and sampling number were the same for all samples. The sampling numbers were defined as m. The wear depth of the historical prior samples was fitted by the GM(1, 1) model, and the corresponding predicted values are given as follows: (15)X=X1,1X1,2⋯X1,mX2,1X2,2⋯X2,mX3,1X3,2⋯X3,m,where Xi,j is the prediction value of sample i at time ti,j, 1 ≤ j ≤ m.Following the grey data processing, the degradation process model of the wear depth was established:(16)Xi,j=α0t+σBBti.j,where α0 is the prior value of the random parameter α and follows a normal distribution; that is, π0α0∼Nμα,0,σα,02. The hyperparameters μα,0 and σα,02 could be calculated from historical prior samples.The wear depth of samplei, Xi = (Xi,1, Xi,2,…, Xi,m)T, follows a multivariate Gaussian distribution. The mean and covariance can be obtained as follows: (17)μ=μα,0t,∑=σα,02ttT+Ω,where t = (t1, t2,…, tm)T, Ω=σB2Q, and(18)Q=t1t1⋯t1t1t2⋯t2⋮⋮⋱⋮t1t2⋯tm.The log-likelihood function can be generated based on the wear data of the historical prior samples:(19)Lθ|X1:3=−3m2ln2π−32ln∑−12∑i=13Xi−μα,0T∑−1Xi−μα,0,where θ is the estimated parameter, θ = (μα, 0, σα, 0, σB), and X1:3 denotes the wear data of the historical prior samples.To simplify the terms in the log-likelihood, we use the results(20)∑=Ω1+σα,02tTΩ−1t,∑−1=Ω−1−σα,021+σα,02tTΩ−1tΩ−1ttTΩ−1.The first partial derivatives ofL(θ|X1:3) for μα,0 and σα,0 can be expressed as(21)∂∂μα,0Lθ|X1:3=∑i=13tTΩ−1Xi−3μα,0tTΩ−1t1+σα,02tTΩ−1t,(22)∂∂σα,0Lθ|X1:3=−3σα,0tTΩ−1t1+σα,02tTΩ−1t+σα,0∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t1+σα,02tTΩ−1t2.The expressions ofμα,0 and σα,0 can be obtained using equations (21) and (22) equal to 0(23)μ^α,0=∑i=13tTΩ−1Xi3tTΩ−1t,(24)σ^α,0=∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t3tTΩ−1t2−1tTΩ−1t1/2.By substituting equations (23) and (24) into equation (19), the contour likelihood function of σB can be given by [25](25)LσB|X1:3=−3m2ln2π−32−32lnΩ−12∑i=13XiTΩ−1Xi−∑i=13tTΩ−1Xi2tTΩ−1t−32ln∑i=13tTΩ−1Xi23tTΩ−1t−∑i=13tTΩ−1X232tTΩ−1t.The maximum of equation (25) can be calculated by the optimization toolbox of MATLAB, and the maximum likelihood estimate of the common parameter σB can be determined as follows:(26)σ^B=5.7095×10−5.By substituting equation (26) into equations (23) and (24), the estimated values of μα,0 and σα,0 can be obtained as follows:(27)μ^α,0=4.9750×10−5,σ^α,0=4.1708×10−6.(2) Real-Time Parameter Updating. The wear data of samples 4, 5, and 6 were selected as the test sets and the random parameter α was updated in real time. Following the grey data processing, the wear process model can be expressed as (28)Xk=αktk+σBBtk,where Xk is the wear depth at time tk of the single sample online and αk is the posterior value of the random parameter α and follows a normal distribution; that is, πkαk∼Nμα,k,σα,k2.It was assumed that the wear depth beforetk is X1:k = (X1, X2,…, Xk)T of the single sample online, and the time interval is Δt. Based on the Bayesian method, the posterior distribution of the random parameter α can be determined by [26, 27](29)παk|X1:k∝LX1:k|αkπ0α0,where the prior information follows a normal distribution, that is, π0α0∼Nμα,0,σα,02, and L(X1:k|αk) is a likelihood function with given random parameters.Based on equation (14) and the basic nature of the Brownian motion, L(X1:k|αk) can be obtained as follows:(30)LX1:k|αk=exp−∑q=1kXq−Xq−1−αkΔt2/2σB2Δt∏q=1k2πσB2Δt.Because bothL(X1:k|αk) and π0(α0) follow a normal distribution, π(αk|X1:k) also follows a normal distribution. Based on equations (29) and (30), the mean and variance of π(αk|X1:k) can be obtained as follows, respectively:(31)μα,k=μα,0/σα,02+∑q=1kXq−Xq−1/σB2∑q=1kΔt2/σB2Δt+1/σα,02,σα,k=1∑q=1kΔt2/σB2Δt+1/σα,02.


## 4.1.1. Model of Wear Process

The degradation process model of the wear depth for the HSLL could be constructed using the Wiener process. The wear depth at timet is defined as X(t) and is expressed as(14)Xt=X0+αt+σBBt,where X(0) is the wear depth at the initial time and X(0) is generally equal to 0. However, for the HSLL, X(0) is not equal to 0 and data conversion must be conducted by X(t)-X(0). Moreover, α is the drift coefficient and follows a normal distribution, that is, α∼Nμα,σα2 [23, 24]; σB is the diffusion coefficient, which can describe the inconsistency and instability of the HSLL life test as well as the effects of certain random factors such as vibration; B(t) is the standard Brownian motion and follows a normal distribution, that is, B(t)∼N(0, t); α is the random parameter; and σB is the common parameter.During the online RUL prediction of the single sample, the parametric estimated value must be accurately obtained, and it is difficult to meet the requirements depending only on the wear depth of the selected sample. The distributed hyperparameter ofα and the estimated value of σB should be obtained from the wear data of historical prior samples. Then, the random parameter α of the selected sample can be realized by online updating using the Bayesian method.


## 4.1.2. Parameter Estimation of Wear Process Model

(1) Prior Parameter Estimation. The parameter estimation of the random parameter α and common parameter σB in the wear process model were executed with samples 1, 2, and 3 considered as historical prior samples. The sampling time and sampling number were the same for all samples. The sampling numbers were defined as m. The wear depth of the historical prior samples was fitted by the GM(1, 1) model, and the corresponding predicted values are given as follows: (15)X=X1,1X1,2⋯X1,mX2,1X2,2⋯X2,mX3,1X3,2⋯X3,m,where Xi,j is the prediction value of sample i at time ti,j, 1 ≤ j ≤ m.Following the grey data processing, the degradation process model of the wear depth was established:(16)Xi,j=α0t+σBBti.j,where α0 is the prior value of the random parameter α and follows a normal distribution; that is, π0α0∼Nμα,0,σα,02. The hyperparameters μα,0 and σα,02 could be calculated from historical prior samples.The wear depth of samplei, Xi = (Xi,1, Xi,2,…, Xi,m)T, follows a multivariate Gaussian distribution. The mean and covariance can be obtained as follows: (17)μ=μα,0t,∑=σα,02ttT+Ω,where t = (t1, t2,…, tm)T, Ω=σB2Q, and(18)Q=t1t1⋯t1t1t2⋯t2⋮⋮⋱⋮t1t2⋯tm.The log-likelihood function can be generated based on the wear data of the historical prior samples:(19)Lθ|X1:3=−3m2ln2π−32ln∑−12∑i=13Xi−μα,0T∑−1Xi−μα,0,where θ is the estimated parameter, θ = (μα, 0, σα, 0, σB), and X1:3 denotes the wear data of the historical prior samples.To simplify the terms in the log-likelihood, we use the results(20)∑=Ω1+σα,02tTΩ−1t,∑−1=Ω−1−σα,021+σα,02tTΩ−1tΩ−1ttTΩ−1.The first partial derivatives ofL(θ|X1:3) for μα,0 and σα,0 can be expressed as(21)∂∂μα,0Lθ|X1:3=∑i=13tTΩ−1Xi−3μα,0tTΩ−1t1+σα,02tTΩ−1t,(22)∂∂σα,0Lθ|X1:3=−3σα,0tTΩ−1t1+σα,02tTΩ−1t+σα,0∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t1+σα,02tTΩ−1t2.The expressions ofμα,0 and σα,0 can be obtained using equations (21) and (22) equal to 0(23)μ^α,0=∑i=13tTΩ−1Xi3tTΩ−1t,(24)σ^α,0=∑i=13Xi−μα,0tTΩ−1ttTΩ−1Xi−μα,0t3tTΩ−1t2−1tTΩ−1t1/2.By substituting equations (23) and (24) into equation (19), the contour likelihood function of σB can be given by [25](25)LσB|X1:3=−3m2ln2π−32−32lnΩ−12∑i=13XiTΩ−1Xi−∑i=13tTΩ−1Xi2tTΩ−1t−32ln∑i=13tTΩ−1Xi23tTΩ−1t−∑i=13tTΩ−1X232tTΩ−1t.The maximum of equation (25) can be calculated by the optimization toolbox of MATLAB, and the maximum likelihood estimate of the common parameter σB can be determined as follows:(26)σ^B=5.7095×10−5.By substituting equation (26) into equations (23) and (24), the estimated values of μα,0 and σα,0 can be obtained as follows:(27)μ^α,0=4.9750×10−5,σ^α,0=4.1708×10−6.(2) Real-Time Parameter Updating. The wear data of samples 4, 5, and 6 were selected as the test sets and the random parameter α was updated in real time. Following the grey data processing, the wear process model can be expressed as (28)Xk=αktk+σBBtk,where Xk is the wear depth at time tk of the single sample online and αk is the posterior value of the random parameter α and follows a normal distribution; that is, πkαk∼Nμα,k,σα,k2.It was assumed that the wear depth beforetk is X1:k = (X1, X2,…, Xk)T of the single sample online, and the time interval is Δt. Based on the Bayesian method, the posterior distribution of the random parameter α can be determined by [26, 27](29)παk|X1:k∝LX1:k|αkπ0α0,where the prior information follows a normal distribution, that is, π0α0∼Nμα,0,σα,02, and L(X1:k|αk) is a likelihood function with given random parameters.Based on equation (14) and the basic nature of the Brownian motion, L(X1:k|αk) can be obtained as follows:(30)LX1:k|αk=exp−∑q=1kXq−Xq−1−αkΔt2/2σB2Δt∏q=1k2πσB2Δt.Because bothL(X1:k|αk) and π0(α0) follow a normal distribution, π(αk|X1:k) also follows a normal distribution. Based on equations (29) and (30), the mean and variance of π(αk|X1:k) can be obtained as follows, respectively:(31)μα,k=μα,0/σα,02+∑q=1kXq−Xq−1/σB2∑q=1kΔt2/σB2Δt+1/σα,02,σα,k=1∑q=1kΔt2/σB2Δt+1/σα,02.


## 4.2. RUL Prediction

The RUL distribution could be obtained using the concept of the first arrival times during the wear process. When the random process {X(t), t ≥ 0} reached the failure threshold, the time was defined as the sample life:(32)T=inft:Xt≥Xf|X0<Xf,where Xf is the failure threshold, X(0) is the wear depth of the initial moment, and X(0) = 0.The residual useful lifeLk in the current runtime tk can be defined as(33)Lk=lk:T−tk|T>tk,where lk is the online RUL value.Because the RUL of the Wiener process with linear drift follows an inverse Gaussian distribution [28], the PDF of Lk can be obtained as follows:(34)fLklk=Xf−Xk2πlk3σα,k2lk+σB2exp−Xf−Xk−μα,klk22lkσα,k2lk+σB2.Based on the observed data, the PDF of the RUL for sample 4 could be calculated and is illustrated in Figure5.Figure 5
PDF of RUL.Figure5 indicates that the PDF curves become narrower with wear data information increases, which indicates the uncertainty of RUL prediction is increasingly small, and the results of RUL prediction will be increasingly accurate.The RUL of a single sample can be written as(35)ELk=EELk|μα,k,σα,k=EXf−Xkμα,k,σα,k=Xf−Xkσα,k2exp−μα,k22σα,k2∫0μα,kexpx22σα,k2dx=2Xf−Xkσα,kDμα,k2σα,k,where D(z) is the Dawson function with respect to z:(36)Dz=exp−z2∫0zexpx2dx,z=μα,k2σα,k.Whenz is sufficiently large, according to the approximate property of the integral, D(z) ≈ 1/(2z). The RUL can be obtained as(37)ELk≈Xf−Xkμα,k.600 min of the wear data was selected as the initial prediction point. The random parameters were updated every 300 min. The RUL of samples 4, 5, and 6 could be predicted with the origin wear depth and grey processing data, respectively. The curves of the actual RUL and the prediction RUL were depicted in Figure6.Figure 6
RUL and confidence interval for HSLL based on original wear depth and grey processing data. (a) Origin wear depth. (b) Grey processing data. (c) Confidence interval of samples 4. (d) Confidence interval of samples 5. (e) Confidence interval of samples 6.
(a)(b)(c)(d)(e)Figure6(a) indicates that the predicted values of the three samples were close at the initial prediction point. The actual and predicted values of sample 5 exhibited better consistency. For samples 4 and 6, which had lower/higher wear rates, the prediction deviation was relatively large. However, the predicted life values gradually became closer to the actual life values and were almost equal to the actual life values at the end of the wear process. Following the wear depth being processed by the GM(1,1) model, as shown in Figure 6(b), the initial prediction point of predicted values were also basically close, but the prediction accuracy of the RUL was improved significantly in the entire residual life cycle.Simultaneously, Figures6(c)–6(e) display the 95% pointwise confidence intervals for the RUL based on these prediction points. Basically, all confidence intervals can cover the actual values of the RUL. The actual RUL line of sample 5 was positioned centrally over the confidence interval, for the origin wear depth, and the confidence intervals of samples 4 and 6 were relatively biased. And compared to the simple Wiener process, the 95% confidence interval of the proposed method in the paper is narrower, which can help reduce uncertainty in residual life prediction results and can better adapt to the individual needs of residual life prediction. These further indicate that the errors of the proposed method are lower.Meanwhile, the RULs of samples 2, 4, and 6 when samples 1, 3, and 5 were selected as the prior samples are shown in Figure7.Figure 7
RUL and confidence interval for HSLL based on original wear depth and grey processing data. (a) Origin wear depth. (b) Grey processing data. (c) Confidence interval of samples 2. (d) Confidence interval of samples 4. (e) Confidence interval of samples 6.
(a)(b)(c)(d)(e)Figure7 had the RUL prediction curves and the 95% pointwise confidence intervals similar to Figure 6. When the differences of the wear rates for samples 2, 4, and 6 were larger, the prediction deviations were relatively large at the initial phase. But the predicted life values gradually became closer to the actual life values with the wear process. Similarly, the overall prediction accuracy of the proposed method is higher than the original wear data.Synthesizing Figures5–7, the PDF of RUL became narrower with the continuous accumulation of data. That is, the uncertainty of RUL prediction becomes smaller and smaller. And the relative error between the prediction of the RUL and the actual RUL is smaller. These happened because the results of model parameter estimation were gradually closer to the reality of sample degradation with the test data increasing.To analyze the prediction precision of the proposed RUL method, multiple evaluation indicators, namely, the root mean square error (RMSE), mean absolute error (MAE), and symmetric mean absolute percentage error (SMAPE), were calculated and can be given as follows:(38)RMSE=1m∑i=1mhti−yi2,(39)MAE=1m∑i=1mhti−yi,(40)SMAPE=1m∑i=1mhti−yihti+yi/2×100%,where m is a positive integer, h(ti) is the predicted value at point i, and yi is the actual RUL at point i.The evaluation indicators of samples 4, 5, and 6 can be calculated using equations (38)–(40), as shown in Table 1.Table 1
The RMSE, MAE, and SMAPE of the predicted values.
Indicator456Original wear dataGrey processing dataOriginal wear dataGrey processing dataOriginal wear dataGrey processing dataRMSE0.62930.25160.20540.23310.93600.4874MAE0.54980.15810.18460.18860.81920.3597SMAPE34.60%9.30%9.60%8.64%26.59%13.33%Table1 shows that the assessment results were similar for RMSE and MAE. For sample 5, the prediction accuracy based on the original wear data was slightly higher than that based on the grey processing data. But for samples 4 and 6, the prediction accuracy of the proposed method was much higher. Based on the SMAPE, the predicted trends of the proposed method were all more accurate.The broken line graph of RMSE, MAE, and SMAPE shown in Figure8 reflects the changes in the RUL prediction accuracy of the predicted values at point i. All three indicators decreased gradually, and the prediction accuracy all gradually increased. Because of the great difference in the wear rates among samples, the values of the three indicators based on the original wear data varied greatly. However, all three indicators based on the grey processing data remained at relatively low levels, and the effectiveness of the proposed method in this paper was demonstrated.Figure 8
RMSE MAE and SMAPE of predicted RUL values. (a) RMSE. (b) MAE. (c) SMAPE.
(a)(b)(c)


## 5. Conclusions

The predictive processing of the wear depth was executed based on the GM(1,1) model. The randomness of the original wear data was reduced and the inherent regularity was clearly reflected. The RUL prediction for the HSLL was accomplished using the Wiener process model and data whitening operation on the original wear depth. Based on the original wear data and simple Wiener process, the prediction accuracy of the RUL exhibited a strong dependence on the prior samples, and the prediction accuracy was relatively low owing to the deviation of the wear rate between the test sample and prior samples. Compared to the simple Wiener process, the prediction accuracy of the proposed method for the RUL was improved significantly in the entire residual life cycle. For 95% confidence interval, the width of an interval for the improved method is narrower than the simple method, which indicates the smaller uncertainty of the model proposed.Although the RUL prediction of self-lubricating spherical plain bearings based on the Wiener process and grey system theory has exhibited a relatively high prediction accuracy, it is important that the origin data is nonnegative and performs well on the smooth performance. And depending on the gradual degradation data rather than the saltation data, the application of the proposed method may be limited; hence, further study is required.



---

*Source: 8843374-2021-01-30.xml*