# Multifractional Brownian Motion and Quantum-Behaved Partial Swarm Optimization for Bearing Degradation Forecasting

Authors: Song Wanqing; Xiaoxian Chen; Carlo Cattani; Enrico Zio
Journal: Complexity
 (2020)

---

## Abstract

Gradual degradation of the bearing vibration signal is usually studied as a nonstationary stochastic time series. Roller bearings are working at high speed in a heavy load environment so that the combination of bearing faults gradually degraded during the rotation might lead to unpredicted catastrophic accidents. The degradation process has the property of long-range dependence (LRD), so that the fractional Brownian motion (fBm) is taken into account for a prediction model. Because of the dramatic changes in the bearing degradation process, the Hurst exponent that describes the fBm will change during the degradation process. A priori Hurst value of the conventional fBm in the prediction is fixed, thus inducing a minor accuracy of the prediction. To avoid this problem, we propose an improved prediction method. Based on the following steps, at the initial data processing, a skip-over factor is selected as the characteristics parameter of the bearing degradation process. A multifractional Brownian motion (mfBm) replaces the fBm for the degradation modeling. We will show that also our mfBm has the same property of long-range dependence as the fBm. Moreover, a time-varying Hurst exponentH(t) is taken to replace the constant H in fBm. Finally, we apply the quantum-behaved partial swarm optimization (QPSO) to optimize H(t) for a finite interval. Some tests and corresponding experimental results will show that our model QPSO + mfBm have a much better performance on the prediction effect than fBm.


---

## Body

## 1. Introduction

ROLLER bearings represent the largest application of rotating mechanisms. Since they are usually working in high speed and heavy load environment, then bearing faults might gradually lead to some catastrophic accidents [1]. Therefore, it is important to study the long-range time forecasting of the bearing life to optimize the service assistance [2].There already exist some models for the analysis of bearing faults. The grey model (GM) was applied to forecast bearing vibration signal [3]. However, under heavy load conditions, the accuracy of forecasting for the bearing vibration signal is usually influenced by many factors, such as deterioration [4], pitting, cracking, etc. Moreover, the GM model lacks a reflection on the randomness of the data; it implies that the GM model should be enhanced. Because of the lacking information, the GM approach only considers the size of the instant forecasting value. Also the hidden Markov model (HMM) has been applied to forecast the bearing fault [5, 6]. But the disadvantage of the HMM is that a state of the HMM depends only on the state itself, and this dependence is time-independent. Besides, the HMM can only predict discrete (step) catastrophe failures and cannot forecast the trend of continuous fault changes. The error back-propagation (BP) algorithm is also widely used in forecasting bearing vibration, but the selection structure of the BP algorithm is not single, its training efficiency is not high, so that overcomplicated structure may cause the overfitting phenomenon. As a result, it is impossible to achieve accurate forecasting results [7]. Compared to the artificial neural network algorithm [8], support vector machine (SVM) is based on the structural risk minimization principle, which has a better performance [9]. However, the SVM needs to set the penalty parameter that makes the prediction model more complex. Moreover, the calculation of the remaining useful life for the bearing is crucial. For this reason, some methods were developed, such as the Gaussian process (GP) model [10] and the recursive maximum likelihood estimation (RMLE) [11]. When using the RMLE to solve the unknown noise in the system caused by nonstationary operating conditions, the noise value is relatively large. The shortcoming of the GP is that the prediction results are performed on a relatively limited data set. Such studies employed the regression-based adaptive predictive models [12], the exponential model [13], and the Extended Kalman Filtering [14].Since the bearing failure is a slowly varying process, the degradation process has a long-range time memory. The past states can affect the future degradation states. This property is called long-range dependent (LRD). There are two kinds of methods to identify whether the process is of LRD. The autocorrelation function (ACF) is nonintegrable [15], or the power spectrum density (PSD) of the process is divergent at the origin [16].Fractional Brownian motion (fBm) is a nonstationary stochastic process with fractal properties of self-similarity and LRD [17] and is widely applied in many forecasting problems, like, e.g., image processing [18] and network traffic analysis [19]. FBm has also some applications in degradation modeling and remaining useful life prediction [20]. Zhang et al. [21] use fBm to describe the degradation of the blast furnace wall. Qin and Lin [22] combine fBm with Delft3D, WRF model, and GIS to predict the trajectory of harmful algal blooms. Xi et al. [23] apply fBm to predict the remaining useful life of lithium batteries. Song et al. [24] use fBm for short-term power load forecasting. Gupta et al. [25] establish the relationship between DCT and discrete fBm to provide a theoretical basis for applying DCT to fBm signal reconstruction. FBm can be considered as a generalized form of the Brownian Motion (BM). Consequently, the forecasting results of the fBm model are more flexible than those obtained from the BM model. FBm is characterized by one parameter called Hurst exponent. The expression range of Hurst exponent is 0<H<1. Except for H=0.5, fBm is an LRD model. Under the condition of H=0.5, fBm degenerates to the BM, which is a model with independent increment process.However, the conventional fBm degradation model has some restrictions. In fact, in the process of bearing degradation, the degradation process at the early stage of the fault is relatively stable, so that there will be a huge oscillation in the later stage. The dramatic change for the bearing degradation process indicates that the Hurst exponent will change at different stages of degradation, thus implying thatH usually has multiple values. The Hurst value of fBm is a constant and will not change again during the current prediction phase.In order to solve this problem, we propose an improved model called multifractional Brownian motion (mfBm). Compared with fBm, the fixed exponentH will be replaced by a time-varying exponent H(t). Moreover the quantum-behaved particle swarm optimization (QPSO) will be proposed to optimize the fundamental parameters [26]. The main idea of QPSO is to optimize the search strategy of the particle swarm optimization (PSO) to get the global optimum solution.In order to extract the data from bearing degradation, the skip-over factor is proposed to measure the tendency of the gradual collapse [27]. In our model, we propose an integrated approach for the mfBm and QPSO for forecasting a bearing skip-over series. The computed H for the skip-over series will be used to deduce parameters in the mfBm model. The QPSO will be used to get the global optimal Ht value in the mfBm model; the optimal can be used to make the accurate forecasting. Finally, we adopt mfBm and mfBm-QPSO models for the forecasting experiment. Monte Carlo method is used to show that the combination of mfBm and QPSO is superior to the fBm model of bearing degradation forecasting.This paper is organized as follows. Section2 describes the fBm model to forecast the bearing skip-over series. The optimization of Ht by QPSO for mfBm model is given in Section 3. Section 4 introduces the computation of the skip-over factor. Some examples and experimental tests are given in Section 5, together with a discussion on our model compared with previous ones. Conclusion is given in Section 6.


## 2. The Forecasting Model of fBm

Let the fBm model be denoted byBHt. Thus, BHt is defined as follows [28]:(1)BH0=b0,BHt−BH0=1ΓH+1/2×∫−∞0t−sH−1/2−−sH−1/2dBs+∫0tt−sH−1/2dBs,where Γ. is the gamma function, τ is the time lag. H is the Hurst exponent. Reported methods for estimating H are variance method, absolute value method, curve fitting method, rescaled range analysis (R/S), periodic graph method, and wavelet method [28–31]. Here R/S is used to estimate the Hurst exponent. The time series of length T is divided into k subinterval of length n, where T=kn. The average value of each subinterval composes a new series Yi. The difference between the maximum and minimum values of the new series and the standard of deviation S2n are calculated. Then the ratio between the two values is as follows:(2)RSn=maxYi−minYiS2n,Yi=∑j=1nXj−Xn¯,where Xn¯ and S2n are the mean and variance of the subinterval, respectively. The value of R/S looks like cnH as n⟶∞, where c is a constant independent of n.The increment process of the fBm called fractional Gaussian noise (fGn) is stationary. The ACF of the fGn is expressed as(3)RHτ=σ2ε2H−22×τε+12H−2τε2H+τε−12H,where ε>0 is used by smoothing fBm so that the smoothed fBm is differentiable. σ2 is the intensity of the fGn, and it is written by(4)σ2=Γ1−2HcosπHπH.Then the fBm can generate a specific stochastic time series as the prediction for next power load series. The stochastic differential equations (SDEs) [32, 33] of the prediction model is written by(5)dXt=μXtdt+σXtdBHt,where μ is the drift parameter, σ is the diffusion coefficient. The Hurst exponent is the most significant parameter in fBm; once it is obtained, the two parameters μ and σ of the fBm model can be calculated by maximum likelihood estimation (MLE) [34, 35].Discretize (5), the increments of the series are in the form(6)ΔXt=μXtΔt+σXtwtΔtH,where wt is unit white noise. The prediction model of the fBm is rewritten as(7)Xt+1=Xt+μXtΔt+σXtwtΔtH.In fBm, the Hurst exponentH is a constant, and it will not change again during the current prediction phase. However, as the degradation progresses, the Hurst exponent is time-varying; thus, the series generated by (6) cannot characterize the nonlinear trend. Therefore, fBm is replaced by the mfBm. For mfBm, Ht will change at different stages of degradation, and it is a time-varying parameter. To further improve the prediction accuracy, the QPSO is used to optimize the Ht of the mfBm.


## 3. Quantum-Behaved Particle Swarm for the Hurst Exponent Optimization

The QPSO introduces quantum computing into the PSO, so that each particle in the space has the quantum behavior. Since each particle of the PSO moves in the search space with a velocity vector to find its local optimum and the global optimum, the disadvantage of the PSO is that it will be trapped in its optimal local solution and cannot reach the optimal global one. To avoid this problem, the particles in QPSO have a probability to appear anywhere in space; thus, QPSO has the ability to avoid falling into a locally optimal solution by searching a globally optimal solution. It has also fewer parameters to be adjusted. The optimization computation ofHt by QPSO will be exhibited as follows.Clerc and Kennedy [36] have proved that each particle must converge to the local attractor pi,j to provide the convergence for the PSO algorithm. They defined it as(8)pi,jt=φPi,jt+1−φPg,jt,φ∼U0,1.Comparing with the convergence of the PSO, the QPSO defines one-dimensional delta potential well for each dimension at the pointpit, and each particle has a quantum behavior, which can be described by the Schrödinger [37]. The distribution function of the particle position is given as(9)FHijt+1=e−2pi,jt−Hi,jt+1/Li,jt.where Li,jt is the standard deviation of the distribution. The position of each particle can be obtained by [38](10)Hi,jt+1=pi,jt±Li,jt2ln1u,where u∼U0,1, mt is defined as the average of a global mean best position:(11)mt=m1t,m2t,…,m2t=1N∑i=1NPi,1t,1N∑i=1NPi,2t,…,1N∑i=1NPi,Dt,,where N is the population size, and pi is the local best position of ith particle. The values of Li,jt and the position Hi,jt+1 are evaluated as follows:(12)Li,jt=2β⋅mjt−Hi,jt,Hi,jt+1=pi,jt±β⋅mjt−Hi,jt⋅ln1u,where β<1.781, which is the only parameter that can guarantee the convergence of the particle. Continuing with the previous example, we use the previous 30 point to forecast next 30 point.(1)
InitializeH, let the number of particles N=20, and each Hi,jt=2∗randint1,1,0,2 with dimension equal to 3, and the number of iterations t=20. The 20 particles of H is shown in Table 1. Initialize the current position, individual optical positions, and the population’s global optimal position. Then each particle is used to find the global optimal solution.(2)
Initialize the local optimal positionPi,j0 for the current position Hi,j0. The mean optimal position mjt is computed by (11). Select a suitable value Î2 and evaluate the objective function value fHi,jt.(3)
The local attractorpi,jt is computed by (8). If rand0,1>0.5, then(13)Hi,jt+1=pi,jt+β⋅mjt−Hi,jt⋅ln1u,else(14)Hi,jt+1=pi,jt−β⋅mjt−Hi,jt⋅ln1u.(4)
The optimal positionpi,jt of each particle Hi,jt is updated by Pg,j=argminPiPi.(5)
The global optimal position is given byHgbest=argminPgPg.(6)
End of iteration.Table 1
The 20 particles ofH.
0.7513
0.8909
0.1493
0.8143
0.1966
0.2551
0.9593
0.2575
0.2435
0.2511
0.5060
0.5472
0.8407
0.9293
0.6160
0.6991
0.1386
0.2543
0.3500
0.4733During the process of updating the particle’s position, it is necessary to have a comparison between the fitness value of the current individual position of each particle and the fitness value of the global optimal position. The loop iteration is used to find the global optimal solution. Once the iteration number or minimum error is reached in QPSO, the optimal approximation power curve for theH can be found. This effective value of H parameter will be used by the mfBm model for the next data series forecasting.


## 4. Computing the Skip-Over Series

The dimensionless parameter skip-over factor is essentially an improvement of variance. Variance is an indicator related to energy. Energy-related indices will be disturbed by drastic operating load changes. The skip-over factor reduces the sensitivity to energy; thus, the gradual degradation trend can be reflected better. The skip-over factor is the average of the sum of squares of the difference between the data and the mean of the minimum and the measure of the difference between the source data and the expected value of the minimum value. Based on the whole-period sampling, the vibration waveforms of bearings often show skip-over characteristic, such as the failure of foundation loosening. With the faults deteriorating, the skip-over factor also changes. The essence of the skip-over factor is amplitude modulation of a waveform. The computation of skip-over factor is as follows:(1)
Firstly, the bearing vibration signals are sampled and normalized:(15)xi=kxixmax,k>1,i=1,2,…,n,wherek is a scale factor.(2)
The vibration signals of the bearing after normalization are divided intox11,x12,…,x1m;…; xn1,xn2,…,xnm.(3)
Letm=1024 [27], among set x11,x12,…,x1m, extract minimum value xmin=x1,p,1<p<m; then get the average value for all m sets:(16)x¯=1n∑i=1nxip.(4)
Calculated variance(17)Dx=1n∑i=1nxip−x¯2.DefineDx as the skip-over factor.Figure1 shows the prediction flowchart of the QPSO + mfBm.Figure 1
The prediction flowchart of the QPSO + mfBm.


## 5. Experiment

The actual data was generated from accelerated life test of roller bearing Center for Intelligent Maintenance Systems (IMS), University of Cincinnati. Four Rexnord ZA-2115 double-row bearings were installed on the shaft as shown in Figure2. The shaft is driven by an AC motor and coupled by rub belts. Accelerometers were installed on the bearing housing. The rotation speed of the shaft was 2000 rpm, the sampling rate was set at 20 kHz, and the data length was 20,480 points. Bearing 1 of vibration signals with an outer race failure is selected.Figure 2
Bearing test rig.Figure3 shows the full-lifecycle tendency of bearing within 7 days of the skip-over series. The skip-over series has relative stability and sensitivity to initial fault. The initial fault appears about 7000 min. Here the ACF method is used to calculate the Hurst exponent of the series. Denote the measured ACF of the skip-over series by rk, Rk be the modeled ACF. For simplicity, we use the normalized ACF. Equation (7) is rewritten as(18)Rk=12k+12H−2k2H+k−12H.Figure 3
Skip-over series for full-life tendency of bearing 7 days.LetM=ER−r2 be the mean square errors. As the M is minimized, the corresponding H value is the Hurst exponent. Figure 4 shows the fitting ACF curve, the calculated Hurst exponent is H=0.849. The high fit of the measured curve to the model curve also indicates that the skip-over series is in accordance with the fBm.Figure 4
The fitting ACF curve. Red line is the ACF of real data; blue line is the ACF of the theoretical curve.TheH of this skip-over series is 0.5<H<1. We focus on fBm, PSO + mfBm, and QPSO + mfBm models. The value of each 30 points is selected to calculate one H, the model parameters μ, and σ of the fBm change with H. Then, we can deduce X^t+1⟵X^t by (6), the fBm can generate the approximation series. The parameter estimation results are shown in Table 2. In this prediction model, the estimation of parameters is affected by two aspects: the observed value and the time interval. The feature extraction process is to integrate multiple points of the original vibration signal into one point, that is, skip-over factor. After feature extraction, the values of skip-over factor are determined, and the time interval is changed. When the time interval is small, the calculated parameter values will increase. The larger the time interval, the smaller the corresponding parameter values. Here, the time interval is set as Δt=0.1, resulting in a large calculated value. Figure 5 shows the prediction by three models. In the fBm prediction, the average relative error of fBm is 4.0%. In PSO + mfBm prediction, the average relative error of mfBm is 3.2%. In QPSO + mfBm prediction, the average relative error of the QPSO + mfBm can be reduced to 2.5%. Obviously, QPSO + mfBm has performed better on prediction accuracy.Table 2
The parameters value of model.
Model
μ
σ
FBm
0.1345
321.9634
PSO + mfBm
78.3264
4.2128e + 04
QPSO + mfBm
97.6945
2.7573e + 05Figure 5
Forecasting with different models.In order to obtain the best prediction accuracy, different combinations of observation data length and prediction step length are tested. Figures6(a)–6(d) show that 30 points to forecast next 20 points, 30 points to forecast next 30 points, 40 points to forecast next 30 points, 50 points to forecast the next 30 points, respectively. The parameter estimation results are shown in Table 3 and the average relative error is shown in Table 4.Figure 6
Forecasting skip-over series difference step long.
(a)
(b)
(c)
(d)Table 3
The parameter value of QPSO + MfBm
Figure1
μ
σ
6(a)
97.6945
2.7573e + 05
6(b)
97.6945
2.7573e + 05
6(c)
106.8199
4.6921e + 05
6(d)
80.4281
8.8612e + 07Table 4
Average relative error of different steps.
a (%)
b (%)
c (%)
d (%)
2.68
2.49
2.59
2.54Obviously, 30 points predict next 30 points has minimum average error in statistic meaning. Figure7 shows prediction results by three models from 600th to 981th points.Figure 7
Forecasting bearing skip-over series by two models.Figure8 shows a box-scatter plot describing the relative error distribution of the forecasting skip-over series. The mean relative error of the QPSO + mfBm is lower than the other two models. Meanwhile, the relative error distribution of the QPSO + mfBm is more concentrated than the other two models; it means the QPSO + mfBm has lower relative error and stable prediction performance.Figure 8
Box-scatter plot of relative error distribution in three models.


## 6. Conclusion

In this paper, a hybrid framework (mfBm and QPSO) for predicting a bearing skip-over series is proposed. We have used the ACF method to estimate the Hurst value of the entire degradation process. The fitting results show that the degradation process of the bearing has the LRD property; moreover, the fBm-based model gives the more reasonable description of the process. However, we have shown that the fBM-based model can be significantly improved by adopting the QPSO + mfBm. By using this improved method, we have shown that the predicted value is more consistent with the actual value, and the average relative error is significantly lower. In terms of relative error distribution, the QPSO + mfBm is more efficient. The experimental results have also shown that the improvement to the LRD prediction method for fBm is effective.



---

*Source: 8543131-2020-01-31.xml*