#### Abstract

The variation trend, failure trajectory, probability distribution, and other information vary with time and working conditions for rolling bearing vibration performance, which makes the evaluation and prediction of the evolution process difficult for the performance reliability. In view of this, the chaos theory, grey bootstrap method, and maximum entropy method were effectively fused to propose a mathematical model for the dynamic uncertainty evaluation of rolling bearing vibration performance. After reconstructing the phase space of the vibration performance time series, four local prediction methods were applied to predict the vibration values of bearings to verify the effectiveness and validity of chaos theory. The estimated true value and estimated interval were calculated using the grey bootstrap method (GBM) and maximum entropy method. Finally, the validity of the proposed model was verified by comparing the probability that the original data fall into the estimated interval with the given confidence level. The experimental results show that the proposed method can effectively predict the variation trend and failure trajectory of the vibration performance time series so as to realize the dynamic monitoring of the evolution process for rolling bearing vibration performance online.

#### 1. Introduction

The performance of rolling bearings mainly includes vibration, noise, tone quality, friction torque, wear, temperature rise, and kinematic accuracy, which vary with the movement process. Vibration performance affects the fatigue life, service precision, and reliability of rolling bearings, which can be used as an important indicator to evaluate the health status of bearings. However, the variation regulations and failure probability distributions of vibration performance are characterized by uncertainty, polytrope, diversity, nonstationarity, and nonlinearity for rolling bearings. Moreover, the variation trend, failure trajectory, probability distribution, and other information of the rolling bearing vibration performance change during the initial degradation stage, gradual degradation stage, rapid degradation stage, and sharp degradation stage, which results in the bizarre evolution of the performance reliability of rolling bearings. They are also affected by many factors, such as machining and assembly and working environment, which brings great difficulty to dynamic prediction and evaluation of the uncertainty of rolling bearing vibration performance. Domestic and foreign research studies have made some achievements in the field of vibration performance assessment and analysis of rolling bearings, which provide theoretical references and some suggestions for this investigation. Ghafari et al. [1] investigated the vibrations of balanced fault-free ball bearings by building a lumped mass-damper-spring model. Considering the impact of roller dynamic unbalance, Cui et al. [2] presented a differential equation of high-speed cylindrical roller bearings. Moreover, the mathematical model considered cage vibration characteristics as a function of roller dynamic unbalance, structural parameters, and working conditions. Dong and Chen [3] researched the Wigner–Ville spectrum based on cyclic spectral density, which was proved useful in extracting the bearing fault pattern from gearbox vibration signals. For an individual rolling bearing without historical lifetime data, it is difficult to carry out vibration performance state assessment by using the traditional method. So Xiao et al. [4] proposed the support evidence statistics method for operation reliability assessment for addressing the problem.

Dynamic uncertainty of rolling bearing vibration performance, also known as instantaneous uncertainty, refers to the extended uncertainty of vibration performance varying with time. It can be expressed by a function, and the specific value of the function is called dynamic uncertainty. Mean uncertainty, as a statistical evaluation index, can be used to evaluate the random fluctuation state of rolling bearing vibration performance. Xia and Meng [5] constructed the grey relational space between the mean of the fluctuant ranges, information entropy, correlation dimension, and box dimension to study the uncertainties of the stochastic process and information source. Cavalini et al. [6] used the fuzzy dynamic analysis to model the inherent uncertainties of the bearing parameters, i.e., the pad radius, the oil viscosity, and the radial clearance (bearing assembly clearance). Based on information poor system theory, Xia and Wang [7] combined the mean of the dynamic fluctuant range with the correlation dimension in chaos theory to address the problem about evaluating the dynamic uncertainty and the nonlinear characteristic of rolling bearings.

There is an early sign for the failure of rolling bearings in the degradation process of vibration performance. Dynamic evaluation is predicting the future sample data according to the historical sample data to evaluate the healthy status of vibration performance at one moment in the future. According to statistical theory, no method can predict the actual performance data accurately. Therefore, estimated interval and uncertainty are introduced to describe the accuracy of the predicted methods.

The vibration data collected in the test vary with time, which can be regarded as a nonlinear time series. Chaos theory, as an important method in the field of nonlinear analysis, plays an indispensable role in the prediction and fault diagnosis for the performance time series of rolling bearings with nonlinear dynamic characteristics. Tian et al. [8], Yau et al. [9], and Wu et al. [10] used chaos theory to carry out fault diagnosis by analyzing vibration signals, which provided an efficient method for monitoring the running state of rolling bearings. Xia and Chen [11] fused the chaos theory and fuzzy set theory to evaluate the nonlinear performance degradation process for rolling bearings. Based on the theoretical and mathematical models established, Chen et al. [12] applied chaos theory to analyze the nonlinear dynamic behavior of the cage on the radial plane, which established a foundation for studying the performance of high-speed cylindrical bearings.

Time delay and embedding dimension are two important parameters for reconstructing the phase space of the time series, which can be calculated by using different approaches. Tian et al. [13] used the C-C algorithm to obtain delay time and calculated the embedding dimension based on the G-P algorithm for a short-term wind speed time series. Based on the condition monitoring data of a methane compressor, Niu and Yang [14] gained the reconstruction parameters (delay time and embedding dimension) by utilizing the C-C method and false nearest neighbor method, respectively, to predict the degradation trend. Zhang et al. [15], Shang et al. [16], Li and Yuan [17], Bing et al. [18], and Li et al. [19] applied the C-C method to get the delay time and embedding dimension of the time series, simultaneously, provided a foundation for reconstructing the phase space. Hence, the C-C method is widely used in phase space reconstruction because of its small amount of computation and easy operation. Furthermore, it can estimate the time delay and embedding dimension simultaneously. The maximum Lyapunov exponent, applied to judge whether the time series has chaotic characteristics, is usually calculated by using the small data set method. Kang [20] resolved the largest Lyapunov exponent and compared it with the bifurcation plot, and the numerical results showed that the theoretical method was correct and effective. Caesarendra et al. [21] applied the largest Lyapunov exponent algorithm in monitoring the condition of low-speed slew bearings, which was proved with laboratory slew bearing vibration data and industrial bearing data.

For the prediction of the time series, traditional methods mainly include mathematical statistics and dynamic methods, but a subjective model has to be established beforehand in the process of prediction. The local prediction method in chaos theory can predict directly according to the chaotic characteristic parameters (such as time delay, embedding dimension, and maximum Lyapunov exponent), which avoids the influence of subjective factors, thereby improving the accuracy and credibility of prediction results. Local prediction methods include the zero-order local prediction method, weighted zero-order local prediction method, first-order local prediction method, weighted first-order local prediction method, and improved weighted first-order local prediction method. Because of the limitations of the zero-order local prediction method, the latter four prediction methods are used more commonly. The four local prediction methods have their advantages and disadvantages, of which the predicted values can be regarded as a small sample with four data. Therefore, effective information needs to be extracted to improve the accuracy of the predicted results. The grey bootstrap model (GBM) combines the grey model (GM) with the bootstrap resampling model, which does not depend on any probability distribution. Because the probability distributions and trends are unknown for the time series, Xia et al. [22] used the bootstrap method to generate a large amount of data (i.e., big data) and then wielded the grey prediction model to predict big data for the current time series with small data. Mitchell et al. [23] fused the grey bootstrap method, grey relation method, and adding-weight one-rank local-region method to predict the time series of the rolling bearing friction torque under the condition of an unknown probability distribution.

For sample data with unknown probability distribution, the maximum entropy method can make unbiased mean and interval estimates effectively, which is widely applied in solving probability density functions of the time series. Das and Zhou [24] combined the maximum entropy method and level set principle to detect the changes in the distribution of quantity. Edwin and Angel [25] provided a clustering method based on the maximum entropy method, which explored the space of all possible probability distributions of the data based on prior information about the clusters. Tang et al. [26] supplied a weight vector model based on the maximum entropy method for conflict management and fault diagnosis. The estimated true value can dynamically describe the change trend of the instantaneous information on rolling bearing vibration performance, which can be compared with the instantaneous value of the original sample data. The estimated interval is a dynamic description of the fluctuation range for the bearing vibration performance. The validity of the prediction model can be verified by observing whether the upper- and lower-bound curves can envelop the fluctuation trajectory of the vibration performance status closely.

In view of this, taking the time series of rolling bearing vibration performance as the research object, the time delay and embedding dimension were calculated based on the C-C method in chaos theory. Then, the small data set method was applied to calculate the maximum Lyapunov exponent. After phase space reconstruction, four local prediction methods were used to predict the vibration information of rolling bearings to verify the effectiveness and feasibility of chaos theory. For the four small sample data of each step, based on the grey bootstrap model, a large number of generated data were simulated. For a given significance level, the maximum entropy method was used to make truth value and interval estimations. Finally, the validity of the proposed model was validated by comparing the probability that the instantaneous value of the original sample data falls into the dynamic estimation interval with the given confidence level.

#### 2. Identification of Chaotic Characteristics

The maximum Lyapunov exponent *λ* can be used to judge whether the time series of rolling bearing vibration performance has chaotic characteristics. If *λ* > 0, it shows that the vibration performance has chaotic characteristics; if *λ* = 0, it shows that the vibration performance has periodic characteristics; if *λ* < 0, it shows that the vibration performance has a fixed point. But the premise is that the phase space is reconstructed based on the time delay *τ* and embedding dimension *m* to find the nearest point of each point and limit the temporary separation.

Time delay and embedding dimension are two important parameters for reconstructing the phase space of the time series, which can be calculated by using different approaches. The autocorrelation method [27], multiple autocorrelation method [28], and mutual information method [29] are commonly used to obtain time delay *τ*. The G-P method [13], false nearest neighbor method [14, 27], and Cao method [29, 30] are widely applied in gaining the embedding dimension. Some researchers assumed that time delay and embedding dimension are not independent of each other, so they utilized the C-C method [15–19] to acquire time delay and embedding dimension synchronously. The comparison between different methods is shown in Table 1.

##### 2.1. C-C Method

During the service period of rolling bearings, the vibration acceleration signals are periodically sampled. The time variable is defined as *t*, and the sampling time interval is Δ*t*. The vibration acceleration data are collected continuously, so the time series vector **X** can be obtained as follows:where *x*_{n} stands for the *n*th data of the original series **X** and *N* is the number of original data.

According to the theory of phase space reconstruction, the phase trajectory of the time series **X** can be obtained aswhere *t* is the *t*th phase trajectory; *x*(*i *+* *(*m *−* *1)*τ*) is the delay value; *m* is the embedding dimension and *τ* is the delay time, which can be calculated by the C-C method synchronously; and *M* is the number of phase points. Phase space reconstruction is the basis for predicting the future evolution of the vibration performance information of rolling bearings.

The C-C method can simultaneously estimate time delay *τ* and time window *τ*_{w}. Time delay *τ* can ensure that the original time series is independent of the embedding dimension *m* and the data components are interdependent on each other. Time window *τ*_{w} is a better variable to estimate the dimension than time delay *τ*, which is also the maximum time of data dependence on each other.

First, the 3 statistical variables statistical , , and *S*_{cor} are constructed as follows:withwhere *t* is the number of segments that the time series is divided into, *t* = 0, 1, 2, …, 200; *m* stands for the embedding dimension; *ε* is the serial number; and *δ* is the standard deviation of the time series.

The time variable *t* corresponding to the first minimum value of curve is the optimal time delay *τ*, and the time variable *t* corresponding to the minimum value of curve *S*_{cor}(*t*) is the required time window *τ*_{w}. Time window *τ*_{w}, optimal time delay *τ*, and embedding dimension *m* should satisfy the following relation:

Based on the solved time window *τ*_{w} and optimal time delay *τ*, embedding dimension *m* can be obtained according to equation (5).

In the process of solving, the following variables also need to be calculated:

The correlation integral of the time series is defined aswith

##### 2.2. Small Data Set Method

According to the solved optimal time delay *τ* and embedding dimension *m*, the phase space is reconstructed to find the nearest point of each point and limit the temporary separation. The distances of adjacent point pairs for each point after *i* discrete time steps are calculated to draw the regression line using the least-squares method.

For the point **X**_{i} on a given orbit, the distance *d*_{i}(0) to its nearest point can be given aswhere *T* is the mean cycle of the time series.

For the point **X**_{i} on a given orbit, the distance *d*_{i}() after discrete time steps is expressed as

The average value *y*() of ln *d*_{i}() can be calculated bywhere is the number of nonzero *d*_{i}(), and the slope of the regression line is the maximum Lyapunov exponent calculated by the least-squares method.

#### 3. Chaotic Prediction Methods

If the vibration time series has chaotic characteristics, then suppose that **X**(*M*) is the center trajectory (viz., the trajectory of prediction started or the phase space trajectory at the end), *L* the reference trajectories similar to the center trajectory, and **X**(*M*_{l}) the *l*th reference trajectory. The weighted zero-order local prediction method, first-order local prediction method, weighted first-order local prediction method, and improved weighted first-order local prediction method are used to reconstruct the phase space as follows [31].

##### 3.1. Weighted Zero-Order Local Prediction Method

Based on the weighted zero-order local method, the evolution rule of the phase trajectory iswithwhere **X**(*M* + 1) is the prediction result; *L* is the number of reference trajectories; *d*_{l} is the Euclidean distance between **X**(*M*) and **X**(*M*_{l}); *d*_{min} is the minimum value of *d*_{l}; and *k* is the prediction parameter, and usually *k* ≥ 1.

##### 3.2. First-Order Local Prediction Method

The most proximal points **X**(*M*)_{reference} around the center trajectory **X**(*M*) are fitted by the linear model **X**(*M* + 1) =*a***R** + *b***X**(*M*) as follows:where the points **X**(*M*_{1}), **X**(*M*_{2}), … , **X**(*M*_{L}) are the proximal points of the center trajectory **X**(*M*).

The coefficients *a* and *b* can be solved by the least-squares method, so **X**(*M* + 1) is obtained according to the formula **X**(*M* + 1) = *a***R** + *b***X**(*M*), and then the prediction value is efficiently separated.

##### 3.3. Weighted First-Order Local Prediction Method

Compared with the first-order local method, the weighted first-order local method considers the influence of weight between each proximal point and the center point; namely, the weight term is added, which can be given as follows:where *k* is the forecasting parameter, and generally *k* ≥ 1.

The linear fitting of the first-order local method can be given bywith **R** = [1, 1, ..., 1]^{T}.

Using the least-weighted-squares method to solve coefficients *a* and *b*, we obtainwithviz.,

The coefficients *a* and *b* are acquired by solving equations, and the prediction value is obtained by equation (16).

##### 3.4. Improved Weighted First-Order Local Prediction Method

The forecasting method for the improved weighted first-order local method is proposed based on the weighted first-order local method. The difference between the two is the definition of correlation between the center trajectory **X**(*M*) and the proximal points or reference trajectories **X**(*M*_{l}): the correlation of the proximal points of the weighted first-order local method is defined by the Euclidean distance, and the correlation of the improved method is defined by the cosine value of the angle as follows:where cos(*l*) is the cosine value of the angle between phase points **X**(*M*) and **X**(*M*_{l}); thus, the calculation process of the improved weighted first-order local method is similar to that of the weighted first-order local method; namely, only the Euclidean distance *d*_{l} is changed into cos(*l*).

#### 4. Grey Bootstrap Method

Using the above four prediction models: weighted zero-order local method, first-order local method, weighted first-order local method, and improved weighted first-order local method, four vibration performance information points for the *ξ*th step predicted backward can be predicted for rolling bearings, with the vector **Y**(*ξ*) expressed aswhere **Y**(*ξ*) is the predicted vibration performance data of the rolling bearing for the *ξ*th step predicted backward, *y*_{ξ}(*u*) is the *u*th data in **Y** for the *ξ*th step, and *ψ* is the maximum value of step number predicted backward for four prediction models.

If *y*_{ξ}(*u*) < 0, then a constant *c* should be selected, making *y*_{ξ}(*u*) + *c* ≥ 0. Therefore, in the actual analysis, **Y** is expressed as

Using the bootstrap method, *B* bootstrap resampling samples of size *z*, namely, the bootstrap resampling samples **V**_{bootstrap}, can be obtained by an equiprobable sampling, as follows:where **V**_{β} is the *β*th bootstrap resampling sample, *β* = 1, 2, …, *B*, and *B* is the times of bootstrap resampling and also the number of bootstrap samples, withwhere (Θ) is the Θth data in the *β*th bootstrap resampling sample.

According to the grey model GM (1, 1) [31], suppose the first-order accumulated generating operator (1-AGO) of **V**_{β} is given by

The grey generated model can be described as the differential equation as follows:where *u* is the time variable and *c*_{1} and *c*_{2} are the coefficients to be estimated.

The increment is used to replace the differential, viz.,where Δ*u* is equal to the unit interval (1). Furthermore, assume the generated vector of the mean series as follows:

The least-squares solution with the initial condition *y*_{ξβ}(1) = (1) is given bywhere the coefficients *c*_{1} and *c*_{2} are as follows:with

According to the inverse AGO, the *β*th generated data are expressed as

Therefore, *B* generated data for the vibration performance can be obtained aswhere is the *β*th generated data.

#### 5. Maximum Entropy Method

Transmuting the generated data **Y**_{B} for the vibration performance into continuous information, the expression of maximum entropy is defined aswhere *p*() is the probability density function of the data series **Y**_{B}. Through the maximum entropy principle [31], the optimal estimation of the density function based on sample information can be obtained, and the main idea of maximum entropy is that the solution is the most “unbiased” among all feasible solutions, as follows:where *S* represents the integral interval, namely, the feasible region, for the random variable .

It satisfies the constraint conditions:where *m*_{q} stands for the *q*th order origin moment, with *γ* for the highest origin moment order.

The entropy can reach its maximum by adjusting *p*(), and the probability density function *p*() can be obtained by using the Lagrange multiplier method, as follows:where *λ*_{0}, *λ*_{1}, …, *λ*_{γ} are the Lagrange multipliers and is a random variable for service accuracy.

With

Equation (37) is the probability density function of generated data **Y**_{B}, which is constructed by the maximum entropy principle. With the help of function *p*(), the true value and the confidence interval of the generated series can be implemented as follows.

According to the probability density function *p*() of the random variable , the estimated true value **X**_{0} of the series **Y**_{B} is given by

If the real number *α *∈* *(0,1), then the probability of is given bywhere is the *α* quantile of the density function *p*() and *α* is the significance level.

For the bilateral quantile, the probability is as follows:where **X**_{U} and **X**_{L} are the upper and lower boundaries of the generated series **Y**_{B}, respectively, and [**X**_{L}, **X**_{U}] is the confidence interval under the *α* level.

Therefore, four points for vibration performance information of rolling bearings are fused by effectively combining the grey bootstrap method with the maximum entropy method, and then the true value **X**_{0} and the interval [**X**_{L}, **X**_{U}] for their vibration performance information can be predicted for each moment in the future.

#### 6. Dynamic Uncertainty Evaluation Method

A subsequence vector **X**_{h} is constructed by intercepting the first *h* adjacent data from the data sequence vector **X** of rolling bearing vibration performance:where *x*_{h}(*u*) is the vibration data in the subsequence vector **X**_{h} at the moment *u*.

Dynamic evaluation is using the data sample before the moment *u* to evaluate the status of vibration performance at the moment *u*.

The fluctuation range of rolling bearing vibration performance can be expressed aswhere *U* is the estimated uncertainty, namely, the instantaneous uncertainty at the confidence level and the moment *t*. Generally, the larger the confidence level , the larger the uncertainty *U*. If = 100%, the uncertainty *U* will reach the maximum, and the estimated result is the most credible. However, the larger the uncertainty *U*, the more the estimated interval [**X**_{L}, **X**_{U}] deviates from the true value, so the more distorted the estimated result.

The number *η* of vibration data falling within the maximum entropy-estimated interval [**X**_{L}, **X**_{U}] is calculated based on the Poisson counting process, and the reliability of the prediction result is defined aswhere is the reliability of chaotic prediction. In general, does not equal the confidence level . The greater the , the higher the reliability of the evaluation results according to the definition of . In statistics and practice, it is better if > .

Thus,where *U*_{mean} is the dynamic average uncertainty and stands for that the calculation process is under the condition = 100%.

Taking into account the minimum uncertainty, the confidence level should satisfy

The statistical variable *U*_{mean} can be used as an effective index for evaluating the uncertainty of the random fluctuation status for rolling bearing vibration performance. In practical analysis, the uncertainty of vibration performance is described by *U*_{mean}, which can also be called dynamic mean uncertainty.

#### 7. Basic Modeling Ideas

A variety of mathematical methods including chaotic prediction methods, the grey bootstrap method, and the maximum entropy method, are used during the process of theoretical modeling. Each model does not exist or function here as an isolated single model but rather by complementing and interlocking with each other, which breaks through the limitations whereby one method can only solve a certain type of problem. The idea works as follows:* *Step 1: based on a time series **X** of rolling bearing vibration signals, the time delay and embedding dimension are calculated based on the C-C method in chaos theory. Step 2: then the small data set method is applied to calculate the maximum Lyapunov exponent to verify whether rolling bearing vibration performance has a chaotic characteristic. Step 3: if the rolling bearing vibration performance has a chaotic characteristic, four local prediction methods are used to predict the vibration information of rolling bearings to verify the effectiveness and feasibility of chaos theory. Step 4: four predicted values for each step in the future are obtained, comprising the small sample **Y** of size four. With the help of the grey bootstrap method, a large sample **Y**_{B} is generated from the small sample **Y**. Step 5: let the large sample **Y**_{B} be continuous data, according to the maximum entropy principle, and then solve the sample moment of each order to obtain the probability density function. The true value **X**_{0} and the interval [**X**_{L}, **X**_{U}] are calculated under the significance level *α*. Step 6: the dynamic uncertainty of rolling bearing vibration performance is predicted and evaluated by four parameters, i.e., the estimated true value, estimated interval, dynamic uncertainty, and average uncertainty. Step 7: the reliability of chaotic prediction is compared with the given confidence level to prove that the chaotic grey bootstrap model is effective for the dynamic prediction and evaluation of rolling bearing vibration performance.

#### 8. Experimental Investigations

##### 8.1. Case 1

The experimental data are collected from Hangzhou Bearing Test & Research Center. This is a strength lifetime test on the vibration performance of rolling bearings, and the test machine model is ABLT-1A, which mainly includes a test head seat, a test head, a transmission system, a loading system, a lubrication system, and a computer control system. The test bearings are angular contact ball bearings with grade P2 and of type 7008AC, which are provided by Luoyang Bearing Science &Technology Co., Ltd. The parameters of the bearings are shown in Table 2.

The research is conducted at a motor speed of 4000 r/min, an axial load of 3.5 kN, and a radial load of 2 kN. Test time and bearing vibration information are automatically recorded by the computer control system: the RMS values of vibration amplitudes are obtained at each 1 min, and then computer collects vibration data once with unit m·s^{−2}; namely, the bearing vibration signal is the RMS value of vibration amplitudes in 1 min. From the beginning of the test, if significant variation occurs in the bearing ring or the roller, or even surface fatigue spalling, the vibration value of the test machine will obviously increase and vibration performance will reduce. If the vibration value reaches a certain value, the motor will stop running, and the experiment will be over.

There are 573 vibration signals of bearing performance in the data sample, which are shown in Figure 1.

As shown in Figure 1, the vibration signals of rolling bearings have obvious characteristics of uncertainty, nonlinearity, and randomness. The first 100 vibration data belong to the initial wear stage, and the vibration signals enhance gradually and fluctuate sharply; the 100–150 vibration data are in the optimal vibration performance stage, and the vibration data are small, which shows the vibration performance is stable; the 150–270 vibration data are in the normal wear stage, and the vibration signals fluctuate sharply; the 270–573 vibration data fluctuate slightly because of smaller wear diameters of bearing components or changes in contact angles between steel balls and rings, improving the vibration performance of bearings. Therefore, the uncertainty of rolling bearing vibration performance should be predicted dynamically according to the complex and changeable vibration information.

Twenty-three steps are predicted to verify the effectiveness of the proposed model based on the previous 550 data. To estimate time delay *τ* and embedding dimension *m* using the C-C method, the three statistical variables , , and *S*_{cor} should be obtained beforehand, as shown in Figure 2.

The time variable *t* corresponding to the first minimum value of curve is the optimal time delay *τ*, and the time variable *t* corresponding to the minimum value of curve *S*_{cor} is the required time window *τ*_{w}. As shown in Figure 2, the first minimum value of curve is the optimal time delay *τ* = 3 min and the optimal estimated value of the mean orbit period is the required time window *τ*_{w} = 40. According to equation (5), the embedding dimension *m* is considered to be 14.

According to the solved optimal time delay *τ* and embedding dimension *m*, the phase space is reconstructed to find the nearest point of each point and limit the temporary separation. The distances of adjacent point pairs for each point after *i* discrete time steps are calculated to draw the regression line using the least-squares method. The average value *y*() of ln *d*_{i}() can be calculated according to equation (11), as shown in Figure 3. The slope of the regression line is the maximum Lyapunov exponent, as shown in Figure 4.

The slope of the regression line in Figure 4 is the maximum Lyapunov exponent *λ* of 0.0028, which shows that the vibration performance of the angular contact ball bearing has chaotic characteristics. Moreover, the chaos method can predict 357 steps accurately at most according to 1/*λ* = 357.14.

Four prediction models: weighted zero-order local method, first-order local method, weighted first-order local method, and improved weighted first-order local method, are used to reconstruct the phase space. According to the previous 550 data, four vibration performance information points for 23 steps backward can be predicted for rolling bearings, which are used for comparison with the actual values to valid the effectiveness of the chaos theory, as shown in Figure 5.

As shown in Figure 5, the curve predicted using the weighted zero-order local method has a linear increasing trend, which is basically consistent with the actual values in the first, third, and ninth steps; the curve predicted using the first-order local method has a linear decreasing trend, which is basically consistent with the actual values in the first, third, eighth, twelfth, fourteenth, and eighteenth steps; the curves predicted using the weighted first-order local method and improved weighted first-order local method are basically kept in a horizontal state, which are basically consistent with the actual values in the second, seventh, tenth, eleventh, fourteenth, sixteenth, seventeenth, eighteenth, nineteenth, twenty-first, and twenty-second steps. In view of this, the four prediction methods have their own advantages, so the effective information on these prediction methods should be fully integrated to accurately predict the vibration performance information for rolling bearings.

Taking the first step as an example, let the times of bootstrap resampling *B* = 20000 and the significant level *α* = 0.05, and the predicted values using four local prediction methods are sampled to obtain the grey bootstrap sample **Y**_{Bootstrap}.

In the process of grey bootstrap generation, the probability density function of the grey bootstrap sample **Y**_{Bootstrap} is calculated using the maximum entropy method, as shown in Figure 6.

Thus, the estimated true value **X**_{0} = 2.41925 ms^{−2} and the estimated interval [**X**_{L}, **X**_{U}] = [2.18166, 2.66972] ms^{−2} are obtained in the first step. Then, the estimated true value and the estimated interval can be acquired in the second step. Similarly, the estimated true value and the estimated interval can be obtained for the other 21 steps. The results are shown in Figure 7.

The number of vibration data falling outside the maximum entropy-estimated interval [**X**_{L}, **X**_{U}] is calculated based on the Poisson counting process. According to equation (44), the reliability of the evaluated results = 100%, indicating that the reliability of chaotic prediction is 100%, which satisfies * *＞* *.

The estimated uncertainties of the rolling bearing vibration performance for the 23 steps in the future are shown in Table 3 according to equation (43).

According to equation (45), the dynamic average uncertainty of rolling bearing vibration performance is calculated as *U*_{mean} = 17.18187/23 = 0.747 ms^{−2}.

##### 8.2. Case 2

This is a strength lifetime test on the vibration performance of rolling bearings, and the test machine model is ABLT-2, which mainly includes a test head seat, a test head, a transmission system, a loading system, a lubrication system, and a computer control system. The test bearings are tapered roller bearings with grade P4 and of type 33022/YA, which are processed using a new material. The parameters of the bearings are shown in Table 4.

The research is conducted at a motor speed of 1000 r/min, an axial load of 45 kN, and a radial load of 93 kN, lubricated by machinery oil no. 32. The load is added three times, and the time interval is one hour. It is added to 30% of the ordinance load for the first time, 60% of the ordinance load for the second time, and 100% of the ordinance load for the third time. The test speed is added two times, and the time interval is two hours. It is added to 50% of the ordinance speed for the first time and 100% of the ordinance speed for the second time. Test time and bearing vibration information are automatically recorded by the computer control system. If the vibration value reaches a certain value, the motor will stop running, and the experiment will be over.

There are 1092 vibration signals of bearing performance in the data sample, which are shown in Figure 8.

As shown in Figure 8, the vibration signals of rolling bearings have obvious characteristics of uncertainty, nonlinearity, and randomness. The uncertainty of rolling bearing vibration performance should be predicted dynamically according to the complex and changeable vibration information.

Twenty-two steps are predicted to verify the effectiveness of the proposed model based on the previous 1070 data. To estimate time delay *τ* and embedding dimension *m* using the C-C method, the three statistics , , and *S*_{cor} should be obtained beforehand, as shown in Figure 9.

The time variable *t* corresponding to the first minimum value of curve is the optimal time delay *τ*, and the time variable *t* corresponding to the minimum value of curve *S*_{cor} is the required time window *τ*_{w}. As shown in Figure 9, the first minimum value of curve is the optimal time delay *τ* = 3 min and the optimal estimated value of the mean orbit period is the required time window *τ*_{w} = 25. According to equation (5), the embedding dimension *m* is considered to be 9.

According to the solved optimal time delay *τ* and embedding dimension *m*, the phase space is reconstructed to find the nearest point of each point and limit the temporary separation. The distances of adjacent point pairs for each point after *i* discrete time steps are calculated to draw the regression line using the least-squares method. The average value *y*() of ln *d*_{i}() can be calculated according to equation (11), as shown in Figure 10. The slope of the regression line is the maximum Lyapunov exponent, as shown in Figure 11.

The slope of the regression line in Figure 11 is the maximum Lyapunov exponent *λ* = 0.0052, which shows that the vibration performance of the tapered roller bearings has chaotic characteristics. Moreover, the chaos method can predict 192 steps accurately at most according to 1/*λ* = 192.31.

Four prediction models: weighted zero-order local method, first-order local method, weighted first-order local method, and improved weighted first-order local method, are used to reconstruct the phase space. According to the previous 1170 data, four vibration performance information points for 22 steps backward can be predicted for rolling bearings, which are used for comparison with the actual values to valid the effectiveness of the chaos theory, as shown in Figure 12.

As shown in Figure 12, the curve predicted using the weighted zero-order local method has a linear increasing trend, which is basically consistent with the actual values in the second, third, fourth, sixth, eighth, and eleventh steps; the curve predicted using the first-order local method is basically consistent with the actual values in the second, fifth, seventh, tenth, and thirteenth steps; the curves predicted using the weighted first-order local method and improved weighted first-order local method are basically kept in a horizontal state, which are basically consistent with the actual values in the fifth, tenth, fourteenth, fifteenth, seventeenth, and twenty-second steps. In view of this, the four prediction methods have their own advantages, so the effective information on these prediction methods should be fully integrated to accurately predict the vibration performance information for rolling bearings.

Taking the first step as an example, let the times of bootstrap resampling *B* = 20000 and the significant level *α* = 0.05, and the predicted values using four local prediction methods are sampled to obtain the grey bootstrap sample **Y**_{Bootstrap}.

In the process of grey bootstrap generation, the probability density function of the grey bootstrap sample **Y**_{Bootstrap} is calculated using the maximum entropy method, as shown in Figure 13.

Thus, the estimated true value **X**_{0} = 3.51486 ms^{−2} and the estimated interval [**X**_{L}, **X**_{U}] = [3.46541, 3.55483] ms^{−2} are obtained in the first step. Then, the estimated true value and the estimated interval can be acquired in the second step. Similarly, the estimated true value and the estimated interval can be obtained for the other 20 steps. The results are shown in Figure 14.

The number of vibration data falling outside the maximum entropy-estimated interval [**X**_{L}, **X**_{U}] is calculated based on the Poisson counting process. According to equation (44), the reliability of the evaluated results = 95.45%, indicating that the reliability of chaotic prediction is 95.45%, which satisfies ≥ .

The estimated uncertainties of the rolling bearing vibration performance for the 22 steps in the future are shown in Table 5 according to equation (43).

According to equation (45), the dynamic average uncertainty of rolling bearing vibration performance is calculated as *U*_{mean} = 20.60374/22 = 0.9365 ms^{−2}.

##### 8.3. Case 3

As shown in Figures 1 and 8, the vibration time series of rolling bearings has obvious increasing, random, and uncertain trends with time. Therefore, the RMS values of vibration data can be seen as the superposition of a linear function (corresponding to the increasing trend) and a normal distribution function (corresponding to the random and uncertain trends). A simulation experiment is used to investigate the fluctuating trend of rolling bearing vibration performance to validate the effectiveness of the theory described above.

To simulate a linear function with intercept 2 and slope 0.01, a total of 100 data points are generated. To simulate a normal distribution function (white Gaussian noise) with mathematical expectation 0 and standard deviation 0.1, a total of 100 data points are generated. Hence, 100 simulated vibration signals are shown in Figure 15.

Twenty steps are predicted to verify the effectiveness of the proposed model based on the previous 80 data. To estimate time delay *τ* and embedding dimension *m* using the C-C method, the three statistical indicators , , and *S*_{cor} can be obtained beforehand, as shown in Figure 16.

The time variable *t* corresponding to the first minimum value of curve is the optimal time delay *τ*, and the time variable *t* corresponding to the minimum value of curve *S*_{cor} is the required time window *τ*_{w}. As shown in Figure 16, the first minimum value of curve is the optimal time delay *τ* = 2 and the optimal estimated value of the mean orbit period is the required time window *τ*_{w} = 2. According to equation (5), the embedding dimension *m* is considered to be 2.

According to the solved optimal time delay *τ* and embedding dimension *m*, the phase space is reconstructed to find the nearest point of each point and limit the temporary separation. The distances of adjacent point pairs for each point after *i* discrete time steps are calculated to draw the regression line using the least-squares method. The average value *y*() of ln *d*_{i}() can be calculated according to equation (11), as shown in Figure 17. The slope of the regression line is the maximum Lyapunov exponent, as shown in Figure 18.

The slope of the regression line in Figure 18 is the maximum Lyapunov exponent *λ* = 0.0007, which shows that the vibration series has chaotic characteristics. Moreover, the chaos method can predict 1428 steps accurately at most according to 1/*λ* = 1428.57.

Four prediction models: weighted zero-order local method, first-order local method, weighted first-order local method, and improved weighted first-order local method, are used to reconstruct the phase space. According to the previous 80 data, four vibration performance information points for 20 steps backward can be predicted for rolling bearings, which are used for comparison with the actual values to valid the effectiveness of the chaos theory, as shown in Figure 19.

As shown in Figure 19, the curve predicted using the weighted zero-order local method has a linear increasing trend, which is basically consistent with the actual values in the second, eighth, tenth, thirteenth, sixteenth, and twentieth steps; the curve predicted using the first-order local method is basically consistent with the actual values in the first, third, eighth, tenth, and twelfth steps; the curve predicted using the improved weighted first-order local method is basically consistent with the actual values in the seventh, eighth, tenth, fourteenth, and seventeenth steps; the curve predicted using the weighted first-order local method is basically consistent with the actual values in the first, third, seventh, eighth, tenth, thirteenth, fourteenth, fifteenth, sixteenth, seventeenth, eighteenth, nineteenth, and twelfth steps. In view of this, the four prediction methods have their own advantages, so the effective information on these prediction methods should be fully integrated to accurately predict the fluctuating trends of vibration information.

Taking the first step as an example, let the times of bootstrap resampling *B* = 20000 and the significant level *α* = 0.05, and the predicted values using four local prediction methods are sampled to obtain the grey bootstrap sample **Y**_{Bootstrap}.

In the process of grey bootstrap generation, the probability density function of the grey bootstrap sample **Y**_{Bootstrap} is calculated using the maximum entropy method, as shown in Figure 20.

Thus, the estimated true value **X**_{0} = 2.81726 and the estimated interval [**X**_{L}, **X**_{U}] = [2.56064, 3.03086] are obtained in the first step. Then, the estimated true value and the estimated interval can be acquired in the second step. Similarly, the estimated true value and the estimated interval can be obtained for the other 18 steps. The results are shown in Figure 21.

The number of vibration data falling outside the maximum entropy-estimated interval [**X**_{L}, **X**_{U}] is calculated based on the Poisson counting process. According to equation (44), the reliability of the evaluated results = 95%, indicating that the reliability of chaotic prediction is 95%, which satisfies ≥ *P*.

The estimated uncertainties of the rolling bearing vibration performance for the 20 steps in the future are shown in Table 6 according to equation (43).

According to equation (45), the dynamic average uncertainty of rolling bearing vibration performance is calculated as *U*_{mean} = 6.02505/20 = 0.30125.

In summary, the four prediction models: weighted zero-order local method, first-order local method, weighted first-order local method, and improved weighted first-order local method, can be applied to forecast the vibration information for rolling bearings. However, it is difficult to find out which method is the best or the worst of the four forecasting models for rolling bearing vibration performance because the models’ advantages or disadvantages are different under different conditions. A single prediction method can only reflect one facet of vibration information in the future, and a prediction value for each step is a character of the true value. Only by fusing and eliciting multiple aspects of the information, can a forecast of the true value be realized. The grey bootstrap method and maximum entropy principle are fused to be applied to the four forecasting results in each step: using the grey bootstrap method for sample processing the four forecast results, a large number of generated data are simulated; through the maximum entropy principle, the probability density of the generated data is calculated, the predicted true value is gained, and the prediction interval is then obtained under a given confidence level. Combining the grey bootstrap method with the maximum entropy principle to fuse the four prediction values of each step effectively, the future true value and estimated interval are obtained for rolling bearing vibration performance. According to the Poisson counting process and grey bootstrap generated data, a dynamic prediction of the uncertainty of rolling bearing vibration performance is achieved for each step in the future, and the prediction results can effectively monitor the variation information on bearing vibration performance.

#### 9. Discussion

The maximum Lyapunov exponent was calculated using chaos theory to verify the chaotic characteristic for rolling bearing vibration performance. For the angular contact ball bearings with grade P2, tapered roller bearings with grade P4, and simulated vibration sequence, the maximum Lyapunov exponents were 0.0028, 0.0052, and 0.0007, respectively, which are more than 0. So the vibration performance of the two kinds of bearings and simulated vibration sequence had a chaotic characteristic.

Each of the four local prediction methods has its own advantages, so the effective information on these prediction methods should be fully integrated to accurately predict the vibration performance information for rolling bearings. Therefore, the grey bootstrap method was fused into the maximum entropy method to estimate the true value and interval for the vibration performance of rolling bearings at each point in the future so as to dynamically predict and evaluate the uncertainty of rolling bearing vibration performance.

The maximum entropy-estimated interval described the fluctuation range of random variables of vibration performance. As shown in Figures 9 and 16, the upper- and lower-bound curves closely enveloped the random fluctuation trajectory of vibration performance. Moreover, the more acutely the vibration performance fluctuates, the greater the dynamic uncertainty, and the dynamic uncertainty had no relationship with the size of the estimated true value.

With the increase of predicted steps, the fluctuation range of the estimated intervals increased gradually, and the larger the values of estimated uncertainty, the more unsatisfactory the prediction results. It also showed the short-term effectiveness of chaotic prediction; that is, the predicted effect was very good in a certain range of predicted steps. As shown in Tables 3 and 4, the dynamic average uncertainties were 0.747 and 0.9365 for the angular contact ball bearings with grade P2 and tapered roller bearings with grade P4, respectively. The dynamic average uncertainty was smaller, and the upper- and lower-bound curves enveloped the curve of the actual vibration value more closely for the angular contact ball bearings with grade P2. It showed that the vibration performance was more stable and fluctuated less acutely for the angular contact ball bearings with grade P2 because their performance accuracy was higher.

Finally, the validity of the proposed model was verified by calculating the probability that the instantaneous values of the original data sample fall within the estimated interval and comparing it with the given confidence level. For the two bearings in the experiment and simulated vibration sequence, the reliability of the evaluated results was 100% and 95.45%, respectively, which were greater than the confidence level 95%. Therefore, the combination of chaos theory and grey bootstrap method was a correct choice.

#### 10. Conclusions

By fusing the chaotic prediction model into the grey bootstrap maximum entropy method, the estimated true value and fluctuation range are predicted for rolling bearing vibration performance at each time point in the future. The dynamic uncertainty of rolling bearing vibration performance is predicted and evaluated by four parameters, i.e., the estimated true value, estimated interval, dynamic uncertainty, and average uncertainty.

The maximum Lyapunov exponent is calculated for the vibration performance time series of rolling bearings by using the small data set method, which verifies whether the vibration performance has chaotic characteristics to lay a foundation for predicting the fluctuation trend and degradation process of rolling bearing vibration performance.

By calculating the maximum entropy-estimated interval of the vibration performance time series, the probability that the instantaneous value of the original data sample falls into the estimated interval is compared with the given confidence level to prove that the chaotic grey bootstrap model is effective for the dynamic prediction and evaluation of the fluctuation trend and degradation process of vibration performance.

For the vibration performance data obtained in experiments, the final analysis results are different for uncertainty because of different analysis models, processing methods, and parameter selection. Therefore, further studies are made to study the influence of different sample sizes and confidence levels on the uncertainty evaluation of vibration performance.

#### Data Availability

The data in case 1 can be accessed by a strength lifetime test on the vibration performance of angular contact ball bearings with grade P2 and of type 7008AC, which are provided by Luoyang Bearing Science &Technology Co., Ltd. The test is performed in Hangzhou Bearing Test & Research Center, Hangzhou, China. The original data are shown in Figure 1, and the data on dynamic analysis are shown in Figures 2–7 within the article. The data in case 2 can be accessed by a strength lifetime test on the vibration performance of tapered roller bearings with grade P4 and of type 33022/YA, which are processed using a new material. The original data are shown in Figure 8, and the data on dynamic analysis are shown in Figures 9–14 within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This project was supported by the National Natural Science Foundation of China (Grant no. 51475144) and Natural Science Foundation of Henan Province of China (Grant no. 162300410065).