A Novel Model for Evaluating the Operation Performance Status of Rolling Bearings Based on Hierarchical Maximum Entropy Bayesian Method

Abstract

Information such as probability distribution, performance degradation trajectory, and performance reliability function varies with the service status of rolling bearings, which is difficult to analyze and evaluate using traditional reliability theory. Adding equipment operation status to evaluate the bearing operation performance status has become the focus of current research to ensure the effective maintenance of the system, reduce faults, and improve quality under the condition of traditional probability statistics. So, a mathematical model is established by proposing the hierarchical maximum entropy Bayesian method (HMEBM), which is used to evaluate the operation performance status of rolling bearings. When calculating the posterior probability density function (PPDF), the similarities between time series regarded as a weighting coefficient are calculated using overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method. The experiment investigation shows that the variation degree of the optimal vibration performance status can be calculated more accurately for each time series relative to the intrinsic series.

1. Introduction

As the key components of rotating machinery, whether rolling bearings are in normal working condition directly affects the running state of the host. As the running environment of rolling bearings becomes more and more complex, changeable, and harsh, the performance analysis and evaluation of rolling bearings are facing serious challenges. Therefore, performance degradation analysis, evaluation, and fault diagnosis are very urgent and necessary for rolling bearings, which will directly affect the safety and stability of the whole host system [1,2,3,4,5,6]. The nonlinear contact and collision between the components (inner race, outer race, rolling element, and cage) of rolling bearings lead to nonlinear and complex dynamic characteristics in the performance degradation progress. Most attention focused on rolling bearing concerns its performance degradation index, so it is generally necessary to evaluate the degradation state of rolling bearing vibration performance effectively.

The performance maintaining relative reliability (PMRR) is used to characterize the performance degradation degree of rolling bearings. Performance maintaining reliability (PMR) is the probability of a rolling bearing running at the optimum performance status, which can be expressed as a function [7]. During the period of optimal vibration performance, there is almost no possibility of performance failure. Not only is the vibration data value small, but also the data fluctuation is not violent. Owing to the interaction and coupling among internal and external factors, it is difficult to solve PMR using the existing dynamic equations and statistical theory.

Adding equipment operation status to evaluate the performance reliability has become the focus of current research to ensure the effective maintenance of the system, reduce faults, and improve quality under the condition of traditional probability statistics [8]. The similarity between time series plays an important role in the analysis and evaluation of performance degradation, which can be calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method [9,10,11]. Therefore, this factor should be considered in the dynamic evaluation process of PMR bearing vibration performance. When calculating the PPDF, the similarity between time series can be regarded as a weighting coefficient, so the variation degree of the optimal vibration performance can be calculated more accurately.

Traditional performance degradation evaluation methods often need complete performance data and assume that the data are deterministic. Weibull distribution plays an important role in the research of bearing reliability and bearing life in the traditional theoretical system. Many scholars build models using Weibull distribution or compare the test results with the calculation results using Weibull distribution, which greatly promotes the development of bearing performance reliability [12,13,14,15]. However, owing to the high reliability and long performance degradation period of rolling bearings, it is difficult to obtain complete vibration performance data in the experiment, and the collected test data have the typical characteristics of poor information. Poor information occurs when the characteristic information of objects researched is incomplete or inadequate [16,17,18,19]. Moreover, the vibration data of rolling bearings are uncertain in the testing process because of the uncertainty factors such as the system error of the testing device, the changeability of the service environment, and so on, which has great limitations for the analysis and evaluation of bearing performance reliability [20,21].

In the existing data-analysis methods, the analysis result of using a single method often has one-sidedness in some aspects, so the evaluation results of PMR can be obtained more comprehensively by fusing several different methods. The bootstrap method can simulate the probability distribution of data samples through re-sampling [22,23,24,25], which can separate the systematic errors in dynamic evaluation process by using the nuclear concept, but it needs to take advantage of the prior information of some rules of data sampling. The maximum entropy principle is used to calculate the probability distribution of data samples while making the subjective estimation error minimum [17,24,26,27]. Bayesian theory fully combines the prior information with the current sample information to obtain the posterior sample information [28,29,30]. The current real-time updating method is based on Bayesian theory, which only uses the observed data sample to update the prior probability distribution, but cannot apply other available information such as the moments of parameters or functions of moments in the probability distribution. Therefore, a variety of methods can be selected for fusion according to the research needs to make up for the limitations of a single method. Ye et al. [7] initially put forward a new concept—accuracy maintaining reliability (AMR) of super-precision rolling bearings, effectively fusing gray bootstrap method and maximum entropy method to predict the failure degree of a bearing successfully maintaining its optimum service accuracy status in the future. However, in the process of calculation, the accuracy threshold was given in advance, which leads to artificial subjective error. So, it is necessary to study a more general real-time performance reliability evaluation method.

In view of this, this paper proposes a HMEBM to establish the reliability evaluation model for evaluating the operation performance status of rolling bearings under the condition that the prior sample information of vibration performance time-series is unknown. The ideas are as follows: based on the vibration acceleration data collected during the service period of rolling bearings, the maximum entropy method was used to calculate the PDF of different time series. Then, the Bayesian method was applied to obtain the PPDF of different time series. Similarities between time series were calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method, which were regarded as weighting coefficients to calculate the PMR of different time series more accurately. Bootstrap method and maximum entropy method were used to calculate the estimated true value and estimated interval of PMRR. Finally, the failure probability of the rolling bearings maintaining the optimal vibration performance status was analyzed. The flow diagram of the proposed method is shown in Figure 1.

2. Mathematical Models

During the service period of rolling bearings, the vibration acceleration data are periodically recorded. The time variable is defined as t, and w time series are obtained as X_w.

$$X_w=\left(x_w\left(1),x_w\right(2),\cdots ,x_w(k),\cdots ,x_w(N)\right);k=1,2,\cdots ,N;w=1,2,\dots ,r$$

(1)

where w stands for the order number of time series X_w; r is the number of time series; x_w(k) is the kth performance data in time series X_w; N is the number of original data.

The intrinsic series is the time series where the performance data are recorded during the optimum vibration performance status for rolling bearings, which is recorded as the first time series and expressed by X₁.

2.1. Solving PDF

Transmuting the performance data into continuous information, the PDF f_w(x) of the time series X_w with maximum entropy is defined as

$$f_w(x)=\exp (c_{0w}+{\displaystyle \sum _{i=1}^j{c}_{iw}}{x}^i)$$

(2)

where c_0w is the first Lagrange multiplier and c_iw is the (i + 1)th Lagrange multiplier for the time series X_w; i is the order number of origin moment, i = 1, 2,…,j; j stands for the highest origin moment order; generally, j = 5.

According to the maximum entropy principle, 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:

$$H_w(x)=-\underset{{\Omega }_w}{\int }{f}_w\left(x\right)\ln f_w(x)\mathrm{d}x\to \max$$

(3)

where H_w(x) is the information entropy and Ω_w represents the feasible domain for the data sample of the time series X_w, and Ω_w = [x_minw, x_maxw]; x_minw and x_maxw are the lower-bound value and upper-bound value in the time series X_w; lnf_w(x) is the logarithmic value of f_w(x).

Equation (3) satisfies the constraint conditions

$${\displaystyle \underset{\Omega w}{\int }{f}_w(x)\mathrm{d}x}=1$$

(4)

$${\displaystyle \underset{{\Omega }_w}{\int }{x}^i{f}_w(x)\mathrm{d}x}=m_{iw}$$

(5)

where m_iw stands for the ith order origin moment for the data sample of the time series X_w, and m_0w = 1.

The entropy can reach its maximum by adjusting f_w(x), and the PDF f_w(x) can be obtained by using the Lagrange multiplier method.

The Lagrange function L_w(x) can be expressed as

$$L_w(x)=H_w(x)+(c_{0w}+1)[{\displaystyle \underset{{\Omega }_w}{\int }{f}_w(x)}\mathrm{d}x-m_{0w}]+{\displaystyle \sum _{i=1}^j{c}_{iw}[{\displaystyle \underset{{\Omega }_w}{\int }{x}^i{f}_w\left(x\right)}\mathrm{d}x-m_{iw}}]$$

(6)

The first Lagrange multiplier can be given by

$$c_{0w}=-\ln [{\displaystyle \underset{{\Omega }_w}{\int }\exp ({\displaystyle \sum _{i=1}^j{c}_{iw}{x}^i})}\mathrm{d}x]$$

(7)

Other j Lagrange multipliers should satisfy the constraint condition

$$1-\frac{{\displaystyle \underset{{\Omega }_w}{\int }{x}^i\exp ({\displaystyle \sum _{i=1}^j{c}_{iw}{x}^i})}\mathrm{d}x}{m_{iw}{\displaystyle \underset{{\Omega }_w}{\int }\exp ({\displaystyle \sum _{i=1}^j{c}_{iw}{x}^i})}\mathrm{d}x}=0$$

(8)

To ensure solution convergence, the original data interval is mapped to interval [-e, e] by the substitution of the variable. Let

$$x=a_{w}t+b_w$$

(9)

where a_w and b_w are mapping parameters for the data sample of the time series X_w; t∈[−e, e], e has a value of 2.71828.

Based on the mapping of the original data in the interval [-e, e], the PDF in Equation (2) can be obtained as

$$f_w(t)=\exp [c_{0w}+{\displaystyle \sum _{i=1}^j{c}_{iw}{(a_{w}t+b_w)}^i}]$$

(10)

In order to be more adaptable to researchers’ habits, the variable t is replaced by the variable x.

2.2. Parameter Estimation

Set a significant level and let α∈(0,1). The confidence level P is given by

$$P=(1-\alpha )\times 100\%$$

(11)

Assume that the maximum entropy estimated interval is [x_Lw, x_Uw] under confidence level P. The lower boundary value x_Lw should satisfy

$$\frac{1}{2}\alpha ={\displaystyle {\int }_{x_{\min w}}^{x_{\mathrm{L}w}}{f}_w\left(x\right)}\mathrm{d}x$$

(12)

The upper boundary value x_Uw should satisfy

$$1-\frac{1}{2}\alpha ={\displaystyle {\int }_{x_{\min w}}^{x_{\mathrm{U}w}}{f}_w\left(x\right)}\mathrm{d}x$$

(13)

2.3. Calculating PPDF

Consider the data sample in the intrinsic series X₁ as the prior sample, which is obtained during the period of the optimal vibration performance status for rolling bearings. In order to study the variation process of the vibration performance reliability for rolling bearings, the PPDF for each time series is constructed according to Bayesian statistics as

$$hy{f}_w(x)=\frac{f_1\left(x\right)f_w\left(x\right)}{{\displaystyle \underset{{\Omega }_{1w}}{\int }{f}_1\left(x\right)f_w(x)\mathrm{d}x}}$$

(14)

where hyf_w(x) is the PPDF of the wth time series X_w; Ω_1w is the intersection of the feasible regions of the data samples for the intrinsic series X₁ and the wth time series X_w.

The similarity between time series plays an important role in the analysis and evaluation of bearing vibration performance variation. Therefore, this factor should be considered in the dynamic evaluation process of PMR of bearing vibration performance. When calculating the PPDF, the similarity between time series can be regarded as a weighting coefficient so that the variation degree of the optimal vibration performance can be calculated more accurately for each time series relative to the intrinsic series. Moreover, the result of analysis using a single method often has one-sidedness in some aspects, so the evaluation results of PMR can be obtained more comprehensively by fusing several different methods.

2.4. Overlapping Area Method

The PMR for vibration performance of rolling bearings is used to express the possibility that rolling bearings run at the optimal vibration performance status. The PMR for the intrinsic series is defined as R₁ = 1. The PMR R_w (1) for the wth time series is expressed using the overlapping area method by

$$R_w\left(1\right)=S_{1w}={\displaystyle \underset{x=x_{1\mathrm{L}}}{\overset{x_{1w}}{\int }}hy{f}_1(x)\mathrm{d}x+}{\eta }_w{\displaystyle \underset{x=x_{1w}}{\overset{x_{2w}}{\int }}hy{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x_{2w}}{\overset{x_{1\mathrm{U}}}{\int }}hy{f}_1(x)\mathrm{d}x}$$

(15)

where S_1w is the overlapping area for the PPDF of the wth time series and the intrinsic series; x_1L and x_1U are the lower and upper bound values of confidence intervals for the PPDF of intrinsic series, respectively; x_1w and x_2w are the abscissa values of the intersections for the PPDF of the wth time series and the intrinsic series, and x_1w < x_2w; η_w is the overlapping area for the PDF of the wth time series and the intrinsic series.

$${\eta }_w={\displaystyle \underset{x=x^{\ast }{}_{{}_{1\mathrm{L}}}}{\overset{x^{\ast }{}_{{}_{1w}}}{\int }}{f}_1(x)\mathrm{d}x+}{\displaystyle \underset{x=x^{\ast }{}_{{}_{1w}}}{\overset{x^{\ast }{}_{{}_{2w}}}{\int }}{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x^{\ast }{}_{{}_{2w}}}{\overset{x^{\ast }{}_{{}_{\mathrm{U}1}}}{\int }}{f}_1(x)\mathrm{d}x}$$

(16)

where x^*_1L and x^*_1U are the lower- and upper-bound values of confidence intervals for the PDF of intrinsic series, respectively; x^*_1w and x^*_2w are the abscissa values of the intersections for the PDF of the wth time series and the intrinsic series, and x^*_1w < x^*_2w.

2.5. Membership Degree Method

X_w is mapped into [0,1] by mapping formula to ensure that all data are fuzzy numbers.

$$z_w\left(k\right)=\frac{x_w\left(k\right)-x_w{}_{\min }}{x_w{}_{\max }-x_w{}_{\min }}$$

(17)

Let Z_w be a fuzzy subset on the finite field Q. The elements between columns have different attributes, and the elements in the same column have the same attributes. Z_w is described as

$$Z_w=\left(z_w\left(1\right),z_w\left(2\right),\dots ,z_w\left(k\right),\dots ,z_w\left(N\right)\right)$$

(18)

When studying the conformity degree of Z_w (w = 1,2,…,r) relative to Z₁, define the absolute difference as

$${\mathsf{\Delta }}_w\left(k\right)=\left|z_w\left(k\right)-z_1\left(k\right)\right|$$

(19)

$$\mathsf{\Delta }{k}_{\max }=\underset{w}{\max } {\mathsf{\Delta }}_w\left(k\right)$$

(20)

The membership functions of elements with the same attributes are established as

$${\mu }_{wk}={\mu }_{wk}\left(z_w\left(k\right),z_1\left(k\right)\right)=1-\frac{{\mathsf{\Delta }}_w\left(k\right)}{\mathsf{\Delta }{k}_{\max }}$$

(21)

The average membership degree is given as

$${\mu }_w=\frac{1}{N}{\displaystyle \sum _{k=1}^N{\mu }_{wk}},w=1,2,\dots ,r$$

(22)

The membership degree can reflect the similarity degree between data series. The larger the average membership degree μ_w is, the more significant the relationship between Z_w and Z₁ is; and on the contrary, the less significant the relationship between Z_w and Z₁ is.

The PMR R_w(2) for the wth time series is expressed using the membership degree method by

$$R_w\left(2\right)={\displaystyle \underset{x=x_{1\mathrm{L}}}{\overset{x_{1w}}{\int }}hy{f}_1(x)\mathrm{d}x+}{\mu }_w{\displaystyle \underset{x=x_{1w}}{\overset{x_{2w}}{\int }}hy{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x_{2w}}{\overset{x_{1\mathrm{U}}}{\int }}hy{f}_1(x)\mathrm{d}x}$$

(23)

2.6. Approach Degree Method

Let X₁ and X_w be fuzzy subsets on the finite field Q. μ_1k∈[0,1] and μ_wk∈[0,1] are the membership degrees of X₁ and X_w, respectively. Minkowski distance d_pk is defined as

$$d_{pk}=d_p\left(X_1,X_w\right)={\left(\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left|{\mu }_{1k}-{\mu }_{wk}\right|}^p}\right)}^{\frac{1}{p}}$$

(24)

where N is the number of elements in sets X₁ and X_w; p is a constant, and generally is considered to be a positive integer.

When p = 1, Minkowski distance becomes Hamming distance; when p = 2, Minkowski distance becomes Euclidean distance.

The approach degree between two fuzzy subsets X₀ and X_i is defined as

$$\beta ={\beta }_p\left(X_1,X_w\right)=\frac{1}{2}\left. [X_1\circ X_w+\left(1-X_1\stackrel{⌢}{\circ }{X}_w\right)\right]$$

(25)

where

$$X_1\circ X_w=\underset{k=1}{\overset{N}{\vee }}\left({\mu }_{1k}\wedge {\mu }_{wk}\right)$$

(26)

$$X_1\stackrel{⌢}{\circ }{X}_w=\underset{k=1}{\overset{N}{\wedge }}\left({\mu }_{1k}\vee {\mu }_{wk}\right)$$

(27)

The Hamming approach degree β_1w is given as

$${\beta }_{1w}={\beta }_1\left(X_1,X_w\right)=1-\frac{1}{N}{\displaystyle \sum _{k=1}^N\left|{\mu }_{1k}-{\mu }_{wk}\right|}$$

(28)

The PMR R_w(3) for the wth time series is expressed using the Hamming approach degree method by

$$R_w\left(3\right)={\displaystyle \underset{x=x_{1\mathrm{L}}}{\overset{x_{1w}}{\int }}hy{f}_1(x)\mathrm{d}x+}{\beta }_{1w}{\displaystyle \underset{x=x_{1w}}{\overset{x_{2w}}{\int }}hy{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x_{2w}}{\overset{x_{1\mathrm{U}}}{\int }}hy{f}_1(x)\mathrm{d}x}$$

(29)

Euclidean approach degree β_2w is defined as

$${\beta }_{2w}={\beta }_2\left(X_1,X_w\right)=1-{\left(\frac{1}{N}{{\displaystyle \sum _{k=1}^N\left|{\mu }_{1k}-{\mu }_{wk}\right|}}^2\right)}^{1/2}$$

(30)

The PMR R_w(4) for the wth time series is expressed using the Euclidean approach degree method by

$$R_w\left(4\right)={\displaystyle \underset{x=x_{1\mathrm{L}}}{\overset{x_{1w}}{\int }}hy{f}_1(x)\mathrm{d}x+}{\beta }_{2w}{\displaystyle \underset{x=x_{1w}}{\overset{x_{2w}}{\int }}hy{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x_{2w}}{\overset{x_{1\mathrm{U}}}{\int }}hy{f}_1(x)\mathrm{d}x}$$

(31)

Cardinal approach degree β_3w is given as

$${\beta }_{3w}={\beta }_3\left(X_1,X_w\right)=\frac{2{\displaystyle \sum _{k=1}^N\left({\mu }_{1k}\wedge {\mu }_{wk}\right)}}{{\displaystyle \sum _{k=1}^N{\mu }_{1k}}+{\displaystyle \sum _{k=1}^N{\mu }_{wk}}}$$

(32)

The PMR R_w(5) for the wth time series is expressed using the Cardinal approach degree method by

$$R_w\left(5\right)={\displaystyle \underset{x=x_{1\mathrm{L}}}{\overset{x_{1w}}{\int }}hy{f}_1(x)\mathrm{d}x+}{\beta }_{3w}{\displaystyle \underset{x=x_{1w}}{\overset{x_{2w}}{\int }}hy{f}_w(x)\mathrm{d}x+}{\displaystyle \underset{x=x_{2w}}{\overset{x_{1\mathrm{U}}}{\int }}hy{f}_1(x)\mathrm{d}x}$$

(33)

The approach degree describes the similarity degree between X₁ and X_w. The larger the approach degree is, the more significant the relationship between X₁ and X_w is; on the contrary, the less significant the relationship between X₁ and X_w is.

2.7. Dynamic Evaluation of PMR

According to the above five values of R_w, small data samples of PMR of rolling bearing vibration performance are obtained for each time series.

$$R_w=\left(R_w\left(1\right),R_w\left(2\right),\dots ,R_w\left(5\right)\right)=\left(R_w\left(\gamma \right)\right);\gamma =1, 2,\dots ,5;w=1,2,\dots ,r;$$

(34)

where R_w is the data sample of PMR for the wth time series; R_w(γ) is the γth data in the PMR data sample for the wth time series.

Using the bootstrap method, B bootstrap re-sampling samples of size z, namely the bootstrap re-sampling samples R_wbootstrap, can be obtained by an equiprobable sampling as

$$R_w{}_{\mathrm{Bootstrap}}=\left(R_w{}_1,R_w{}_2,\dots ,R_w{}_{\theta },\dots ,R_w{}_B\right)$$

(35)

where R_wθ is the θth bootstrap re-sampling sample, θ = 1, 2, …, B; B is the times of the bootstrap re-sampling, and also the number of bootstrap samples, with

$$R_w{}_{\mathsf{\theta }}=\left. [R_w{}_{\mathsf{\theta }}\left(\mathit{\Theta }\right)\right];\mathit{\Theta }=1,2,\dots ,z$$

(36)

where R_wθ(Θ) is the Θth data in the θth bootstrap re-sampling sample of PMR for the wth time series.

The maximum entropy method is used to calculate the PDF of the generated sample R_wbootstrap. According to the PDF, the true value and upper and lower bound values are estimated for the data sample of PMR of each time series.

The PDF of the PMR data sample can be calculated as

$$f(R_w)=\exp [c^{\ast }{}_{0w}+{\displaystyle \sum _{i=1}^j{c}^{\ast }{}_{iw}{(a^{\ast }{}_w{R}_w+b^{\ast }{}_w)}^i}\left)\right]$$

(37)

where c^*_0w is the first Lagrange multiplier and c^*_iw is the (i + 1)th Lagrange multiplier for the time series R_w; i is the order number of origin moment, i = 1,2,…,j; j stands for the highest origin moment order; generally, j =5.where a^*_w and b^*_w are mapping parameters for the data sample of the time series R_w.

The estimated true value of the PMR is obtained as

$$R_{w0}={\displaystyle \underset{S}{\int }{R}_{w}f\left(R_w\right)}\mathrm{d}{R}_w$$

(38)

Set a significant level and let α∈(0,1). The maximum entropy estimated interval is given as

$$[R_{w\mathrm{L}},R_w{}_{\mathrm{U}}]=[R_w{}_{\frac{\alpha }{2}},R_w{}_{{{}_{1-}}_{\frac{\alpha }{2}}}]$$

(39)

with

$$\frac{\alpha }{2}={\displaystyle {\int }_{R_{w\min }}^{R_{w_{\frac{\alpha }{2}}}}f\left(R_w\right)}\mathrm{d}{R}_w$$

(40)

$$1-\frac{\alpha }{2}={\displaystyle {\int }_{R_{w\min }}^{R_{w_{1-\frac{\alpha }{2}}}}f\left(R_w\right)}\mathrm{d}{R}_w$$

(41)

where R_wL is the lower bound value and R_wU is the upper bound value of the PMR for the wth time series.

The variation probability of the PMR for different time series relative to the intrinsic series is defined by

$$P_w\left(\gamma \right)=1-R_w\left(\gamma \right)$$

(42)

According to the concept of relative error in measurement theory, the failure degree of rolling bearings running at the optimal vibration performance status—that is, the PMRR for the vibration performance of rolling bearings—is expressed by

$$d_w\left(\gamma \right)=\frac{R_w\left(\gamma \right)-R_{10}}{R_{10}}\times 100\%$$

(43)

where R₁ is the PMR for the intrinsic series of rolling bearings.

d_w(γ) < 0 indicates that the PMR for the time interval corresponding with the wth time series is less than the PMR for the time interval corresponding with the intrinsic series. The smaller d_w(γ) is, the greater the failure probability of the rolling bearings maintaining the optimal vibration performance status to work for the time interval corresponding with the wth time series.

$$d_{w0}=\frac{R_{w0}-R_{10}}{R_{10}}\times 100\%$$

(44)

$$d_{w\mathrm{L}}=\frac{R_{w\mathrm{L}}-R_{10}}{R_{10}}\times 100\%$$

(45)

$$d_{w\mathrm{U}}=\frac{R_{w\mathrm{U}}-R_{10}}{R_{10}}\times 100\%$$

(46)

where d_w0 is the estimated true value, d_wL is the lower-bound value and d_wU is the upper-bound value of the PMRR for the wth time series.

The flow diagram of the proposed method is shown in Figure 1.

3. Experimental Verification

3.1. Case 1

This is a strength lifetime test on the vibration performance of rolling bearings. Experimental data are collected in the whole life-cycle of rolling bearings in Hangzhou Bearing Test & Research Center. The test machine model is an ABLT-1A. This machine mainly consists of a test head seat, test head, transmission system, loading system, lubrication system, and a computer control system. The physical drawings of the testing machine and test head are shown in Figure 2 and Figure 3a–c. The test bearings and support bearings are angular contact ball bearings with grade P2 and type 7008AC provided by Luoyang Bearing Science &Technology Co., Ltd. (Luoyang, China). The inner ring, outer ring, and rolling elements are made of high carbon chromium bearing steel, which has the following characteristics: density 7.8 g/mm³, elastic modulus 2.08 × 10⁵ N/mm², Poisson’s ratio 0.3 and Hardness 700 HV10. The bearing parameters are shown in

Table 1. Parameters of angular contact ball bearings with type 7008AC.

Parameters	Values	Parameters	Values
Inner diameter d/mm	40	Number of balls	19
Outer diameter D/mm	68	Ball diameter/mm	7.144
Width B/mm	15	Contact angle/°	25

.

The research is conducted at a motor speed 6000 r/min, an axial load of 4.17 kN, and a radial load of 4.58 kN. The DAQ board type is PCI-1711U, and the data-acquisition rate is 20 KHz. The vibration acceleration sensor used is a YD-1 piezoelectric sensor produced by Far East Vibration System Engineering Technology Co., Ltd. (Beijing, China), with a measuring range of ± 2000 g and a resolution of 0.0001 g. The RMS values of vibration amplitudes are obtained by a 1 min interval. The computer collects the vibration data in the RMS value at this interval with a unit of m·s⁻². If significant variation occurs in the bearing ring or the roller, or even if surface fatigue spalling is noticed, the vibration value of the test machine will obviously increase, and the vibration performance will reduce. If the vibration value reaches 25 m·s⁻², the motor will stop running and the experiment will be terminated. (The vibration threshold limit is set according to the service condition of bearings and experience of experimental operators.) The vibration data are automatically collected by the computer control system, as shown in Figure 4.

The vibration time series of rolling bearings have obvious nonlinearity, randomness, and uncertainty in Figure 4. Therefore, the vibration performance reliability of bearings should be predicted and evaluated dynamically based on the complex and changeable vibration information. The vibration data from the 1st to the 212th point range from 0.3 to 2.5 m·s⁻². The vibration performance of the bearings appears stable during this period. The vibration data from the 213th to the 472nd point range from 2.5 to 5.5 m·s⁻². During this period, the vibration signals are stronger and highly fluctuating. The vibration data from the 473rd to the 2574th point range from 2.5 to 5.1 m·s⁻². During this period, the vibration performance of the bearings is stable. The values of vibration data from the 2575th to the 6659th point range from 3 to 5.5 m·s⁻². During this period, the vibration signals of bearings are strengthened. The values from the 6660th to the 7446th data points are around 3.4–7.7 m·s⁻². During this period, the vibration signals are stronger and highly fluctuating. The values of the 7447th to 7793rd data points are in the range of 8.2–22.9 m·s⁻². During this period, the vibration signals increase linearly. Therefore, it can be considered that the vibration performance of the bearings is at the initial wear stage during the period corresponding to the 1st to the 472nd data point. The vibration performance is at the optimal vibration performance state during the period corresponding to the 473rd to the 2574th data point, which is regarded as the first time series (intrinsic series). The vibration performance is at the normal wear stage during the period corresponding to the 2575th to the 6659th data point, which is regarded as the second time series. The vibration performance is at the degeneration stage during the period corresponding from the 6660th to the 7446th data point, which is regarded as the third time series. The vibration performance is at the deterioration stage during the period corresponding to the 7447th to the 7793rd data point, which is regarded as the fourth time series.

3.1.1. PDF of Data Samples of Time Series

For the first time series, various order origin moments can be obtained using the maximum entropy method as [m₁₁, m₂₁, m₃₁, m₄₁, m₅₁] = [−0.7710, 1.1829, −1.6809, 3.0058, −5.0425]; the Lagrange multipliers [c₀₁, c₁₁, c₂₁, c₃₁, c₄₁, c₅₁] = [−0.526176, −1.0192, −0.8837, −0.3003, 0.0325, 0.0628]; the mapping parameters a₁ = 1.8819 and b₁ = −7.1511. Set the significance level α is 0.01; that is, the confidence level P = 99%. The maximum entropy estimated interval is [2.5507, 4.4738] m·s⁻² for the first time series. According to Equation (10), the probability density estimated function f₁(x) is calculated as shown in Figure 5.

For the second time series, various order origin moments can be obtained using the maximum entropy method as [m₁₂, m₂₂, m₃₂, m₄₂, m₅₂] = [−0.4015, 0.8709, −0.7686, 1.6745, −1.9542]; the Lagrange multipliers [c₀₂, c₁₂, c₂₂, c₃₂, c₄₂, c₅₂] = [−0.3629, −0.3569, −0.2062, −0.1813, −0.1369, 0.0291]; the mapping parameters a₂ = 1.9572 and b₂ = −8.3179. The significance level α is set as 0.01; that is, the confidence level P = 99%. The maximum entropy estimated interval is [3.1823, 5.1008] m·s⁻² for the second time series.

For the third time series, various order origin moments can be obtained using the maximum entropy method as [m₁₃, m₂₃, m₃₃, m₄₃, m₅₃] = [−0.8840, 1.5477, −2.2030, 4.5372, −7.4059]; the Lagrange multipliers [c₀₃, c₁₃, c₂₃, c₃₃, c₄₃, c₅₃] = [−1.4292, −1.1189, −0.2384, 0.0274, −0.0268, 0.0145]; the mapping parameters a₃ = 1.1379 and b₃ = −6.3153. The significance level α is set as 0.01; that is, the confidence level P = 99%. The maximum entropy estimated interval is [3.3405, 7.3435] m·s⁻² for the third time series.

For the fourth time series, various order origin moments can be obtained using the maximum entropy method as [m₁₄, m₂₄, m₃₄, m₄₄, m₅₄] = [−0.1634, 1.2181, −0.2573, 3.2716, −0.7843]; the Lagrange multipliers [c₀₄, c₁₄, c₂₄, c₃₄, c₄₄, c₅₄] = [−2.3395, −0.5070, −0.0421, 0.2472, −0.0660, −0.0299]; the mapping parameters a₄ = 0.3328 and b₄ = −5.1758. The significance level α is set as 0.01; that is, the confidence level P = 99%. The maximum entropy estimated interval is [8.1599, 22.5150]m·s⁻² for the fourth time series.

The PDF f_w(x) of four time series are shown in Figure 6.

3.1.2. PPDF of Data Samples of Time Series

Intersection interval [x_w1, x_w2] and intersection area η_w are calculated for the PDF of the wth time series and the intrinsic series, as shown in

Table 2. Intersection intervals and overlapped areas of PDF.

Sequence Number w	Intersection Interval [xw1, xw2]/(m·s−2)	Intersection Area ηw
1	[2.5000, 5.1000]	1
2	[3.6312, /]	0.3590
3	[3.9596, /]	0.0730
4	[/, /]	0

. The similarity degrees between time series η_w, μ_w, β_1w, β_2w and β_3w are calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method, as shown in

Table 3. Similarity degrees calculated using different methods.

Sequence Number w	Similarity Degrees
	Membership Degree μw	Euclidean Approach Degree β1w	Hamming Approach Degree β2w	Cardinal Approach Degree β3w
1	1	1	1	1
2	0.7897	0.7538	0.7897	0.8825
3	0.6094	0.6036	0.6094	0.7573
4	0.0023	0.0017	0.0023	0.0045

.

The experimental data in Figure 3 show that the vibration performance has obvious different uncertainty and nonlinearity for different fault diameters, which belongs to the poor-information problem with an unknown trend. In order to study the variation trend of PMR for the vibration performance, the product functions f₁(x)f_w(x) and posterior probability density functions hyf_w(x) of the four time series are constructed using the Bayesian principle, as shown in Figure 7 and Figure 8.

3.1.3. PMR of Time Series

The intersection interval [x_1w, x_2w] and intersection area S_w are calculated for the PPDF of the wth time series and the intrinsic series, as shown in

Table 4. Intersection intervals and overlapped areas of PPDF.

Sequence Number w	Intersection Interval [xw1, xw2]/(m·s−2)	Intersection Area ηw
1	[2.5000, 5.1000]	1
2	[3.3579, /]	0.6341
3	[3.5821, /]	0.4107
4	[/, /]	0

. The values of PMR, R_w(1), R_w(2), R_w(3), R_w(4), and R_w(5) are calculated using overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method for different time series, as shown in

Table 5. PMR calculated using different methods for different time series.

Sequence Number w	Values of PMR/%
	Overlapping Area Method Rw(1)	Membership Degree Method Rw(2)	Hamming Approach Degree Method Rw(3)	Euclidean Approach Degree Method Rw(4)	Cardinal Approach Degree Method Rw(5)
1	100.00	100.00	100.00	100.00	100.00
2	59.82	61.35	60.99	61.35	62.26
3	36.71	38.18	38.14	38.18	39.27
4	0.00	0.00	0.00	0.00	0.00

. The PMR of each time series decreases gradually, and the slope of the curve decreases gradually, as shown in Figure 9.

3.1.4. PMRR of Time Series

The values of PMRR, d_w(1), d_w(2), d_w(3), d_w(4), and d_w(5) are calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method, as shown in

Table 6. PMRR calculated using different methods.

Sequence Number w	Values of PMRR/%
	Overlapping Area Method dw(1)	Membership Degree Method dw(2)	Hamming Approach Degree Method dw(3)	Euclidean Approach Degree Method dw(4)	Cardinal Approach Degree Method dw(5)
1	0	0	0	0	0
2	−40.18	−38.65	−39.01	−38.65	−37.74
3	−63.29	−61.82	−61.86	−61.82	−60.73
4	−100	−100	−100	−100	−100

.

3.1.5. Fusion Results of Multiple Weighting Methods

In the process of bootstrap generation, take the sample data of PMR of the second time series as an example. Let the sampling number q = 5, the times for the bootstrap re-sampling B = 10,000, and the significance level α = 0. R₂ = (R₂(1), R₂(2), R₂(3), R₂(4), R₂(5)) = (59.82%, 61.35%, 60.99%, 61.35%, 62.26%). The bootstrap re-sampling samples R_2bootstrap can be obtained by an equiprobable sampling, as shown in Figure 10.

The PDF of the PCR data sample can be calculated for the second time series, as shown in Figure 11. Thus, the estimated true value of PMR R₂₀ = 61.18% and the maximum entropy estimated interval [R_2L, R_2U] = [59.68%, 62.39%] are obtained for the sample data of PMR of the second time series. Similarly, the estimated true values and the estimated intervals can be obtained for the sample data of PMR of the other three time series. The results are shown in

Table 7. Estimated true values and estimated intervals for different time series.

Sequence Number w	Estimated True Value Rw0/%	Estimated Intervals [R2L, R2U]/%
1	100	/
2	61.18	[59.68, 62.39]
3	38.13	[36.71, 39.55]
4	0	/

.

The calculated results of the estimated truth curve and the upper- and lower-bound curve are shown for different time series in Figure 12. The estimated true value, lower-bound value and upper-bound value of PMR decrease gradually and have trends of nonlinear reduction over time. From the first time series to the second time series, the estimated true value of PMR is reduced from 100% to 61.18% very rapidly; from the second time series to the third time series, the estimated true value of PMR is reduced from 61.18% to 38.13% slowly; from the third time series to fourth time series, the estimated true value of PMR is reduced from 38.13% to 0 rapidly, relatively.

The PMRR for vibration performance of rolling bearings are calculated for different time series, as shown in

Table 8. Estimated true values and estimated intervals of PMRR.

Sequence Number w	Estimated True Value dw0/%	Estimated Intervals [d2L, d2U]/%
1	0	[/, /]
2	−38.82	[−40.32, −37.61]
3	−61.87	[−63.29, −60.45]
4	−100	[/, /]

and Figure 13. The estimated true value, lower-bound value and upper-bound value of PMRR decrease gradually and have trends of nonlinear reduction. From the first time series to second time series, the estimated true value of PCRR is very rapidly reduced from 0% to −38.82%; from the second time series to the third time series, the estimated true value of PCRR is slowly reduced from −38.82% to −61.87%; from the third time series to the fourth time series, the estimated true value of PCRR is reduced from −61.87% to −100% relatively rapidly.

This analysis result is more in line with engineering practice, because the prior information and current sample information of the vibration performance degradation information of bearings in service were comprehensively considered by the proposed method. In addition, the HMEBM also considers the similarity between time series, which contains the hidden information of the degradation process of optimal vibration performance state.

3.2. Case 2

The test machine and the bearing used in this case are exactly the same as those of Case 1. The vibration data are shown in Figure 14 by changing the test conditions of the motor to a speed of 4000 r/min, an axial load of 4.17 kN, and a radial load of 4.58 kN.

3.2.1. PDF of Data Samples of Time Series (Case 2)

For the first time series, various order origin moments can be obtained using the maximum entropy method as [m₁₁, m₂₁, m₃₁, m₄₁, m₅₁] = [−0.0463, 0.8915, −0.2456, 2.1047, −0.7743]; the Lagrange multipliers [c₀₁, c₁₁, c₂₁, c₃₁, c₄₁, c₅₁] = [0.6246, 0.3796, −0.4400, −0.3470, −0.0294, 0.0463]; the mapping parameters a₁ = 4.4481 and b₁ = −9.1186. The significance level α is set as 0.01; that is, the confidence level P = 99%. The maximum entropy estimated interval is [1.5354, 2.5003] m·s⁻² for the first time series.

The PDF f_w(x) of four time series are shown in Figure 15.

3.2.2. PPDF of Data Samples of Time Series (Case 2)

The intersection interval [x_w1, x_w2] and intersection area η_w are calculated for the PDF of the wth time series and the intrinsic series as shown in

Table 9. Intersection intervals and overlapped areas of PDF (Case 2).

Sequence Number w	Intersection Interval [xw1, xw2]/(m·s−2)	Intersection Area ηw
1	[1.5000, 2.6000]	1
2	[2.2776, /]	0.6272
3	[2.3655, /]	0.4404
4	[/, /]	0

. The similarity degrees between time series η_w, μ_w, β_1w, β_2w, and β_3w are calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method, as shown in

Table 10. Similarity degrees calculated using different methods (Case 2).

Sequence Number w	Similarity degrees
	Membership Degree μw	Euclidean Approach Degree β1w	Hamming Approach Degree β2w	Cardinal Approach Degree β3w
1	1	1	1	1
2	0.9053	0.8718	0.9053	0.9503
3	0.7138	0.6722	0.7138	0.8329
4	0.0066	0.0049	0.0066	0.0131

.

The experimental data in Figure 14 show that the vibration performance has obvious different uncertainty and nonlinearity for different fault diameters, which belongs to the poor-information problem with an unknown trend. In order to study the variation trend of PMR for the vibration performance, the product functions f₁(x)f_w(x) and posterior probability density functions hyf_w(x) of the four time series are constructed using the Bayesian principle, as shown in Figure 16 and Figure 17.

3.2.3. PMR of Time Series (Case 2)

The intersection interval [x_1w, x_2w] and intersection area S_w are calculated for the PPDF of the wth time series and the intrinsic series, as shown in

Table 11. Intersection intervals and overlapped areas of PPDF (Case 2).

Sequence Number w	Intersection Interval [xw1, xw2]/(m·s−2)	Intersection Area ηw
1	[1.5000, 2.6000]	1
2	[1.7124, 1.9161, 2.2126]	0.7473
3	[1.7431, 1.9981, 2.2740]	0.6751
4	[/, /]	0

. The values of PMR, R_w(1), R_w(2), R_w(3), R_w(4), and R_w(5) are calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method for different time series, as shown in

Table 12. PMR calculated using different methods for different time series (Case 2).

Sequence Number w	Values of PMR/%
	Overlapping Area Method Rw(1)	Membership Degree Method Rw(2)	Hamming Approach Degree Method Rw(3)	Euclidean Approach Degree Method Rw(4)	Cardinal Approach Degree Method Rw(5)
1	100.00	100.00	100.00	100.00	100.00
2	64.30	72.08	71.14	72.08	73.34
3	57.54	62.41	61.67	62.41	64.53
4	0.00	0.00	0.00	0.00	0.00

. The PMR of each time series decreases gradually, and the slope of the curve decreases gradually, as shown in Figure 18.

3.2.4. PMRR of Time Series (Case 2)

The values of PMRR, d_w(1), d_w(2), d_w(3), d_w(4), and d_w(5) are calculated using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method, as shown in

Table 13. PMRR calculated using different methods (Case 2).

Sequence Number w	Values of PMRR/%
	Overlapping Area Method dw(1)	Membership Degree Method dw(2)	Hamming Approach Degree Method dw(3)	Euclidean Approach Degree Method dw(4)	Cardinal Approach Degree Method dw(5)
1	0	0	0	0	0
2	−35.7	−27.92	−28.86	−27.92	−26.66
3	−42.46	−37.59	−38.33	−37.59	−35.47
4	−100	−100	−100	−100	−100

.

3.2.5. Fusion Results of Multiple Weighting Methods (Case 2)

In the process of bootstrap generation, take the sample data of PMR of the second time series as an example. Let the sampling number q = 5, the times for the bootstrap re-sampling B = 10,000, and the significant level α = 0. R₂ = (R₂(1), R₂(2), R₂(3), R₂(4), R₂(5)) = (59.82%, 61.35%, 60.99%, 61.35%, 62.26%). The bootstrap re-sampling samples R_2bootstrap can be obtained by an equiprobable sampling, as shown in Figure 19.

The PDF of the PCR data sample can be calculated for the second time series, as shown in Figure 20. Thus, the estimated true value of PMR R₂₀ = 61.18% and the maximum entropy estimated interval [R_2L, R_2U] = [59.68%, 62.39%] are obtained for the sample data of PMR of the second time series. Similarly, the estimated true values and the estimated intervals can be obtained for the sample data of PMR of the other three time series. The results are shown in

Table 14. Estimated true values and estimated intervals for different time series (Case 2).

Sequence Number w	Estimated True Value Rw0/%	Estimated Intervals [R2L, R2U]/%
1	100	/
2	72.16	[71.02, 73.46]
3	62.72	[61.51, 64.69]
4	0	/

.

The calculated results of the estimated truth curve and the upper and lower bound curve are shown for different time series in Figure 21. The estimated true value, lower-bound value and upper-bound value of PMR decrease gradually and have trends of nonlinear reduction over time. From the first time series to the second time series, the estimated true value of PMR is rapidly reduced from 100% to 72.16%; from the second time series to the third time series, the estimated true value of PMR is slowly reduced from 72.16% to 62.72%; from the third time series to the fourth time series, the estimated true value of PMR is reduced from 62.72% to 0 very rapidly.

The PMRR for vibration performance of rolling bearings are calculated for different time series, as shown in

Table 15. Estimated true values and estimated intervals of PMRR (Case 2).

Sequence Number w	Estimated True Value dw0/%	Estimated Intervals [d2L, d2U]/%
1	0	[/, /]
2	−27.84	[−28.98, −26.54]
3	−37.28	[−38.49, −35.31]
4	−100	[/, /]

and Figure 22. The estimated true value, lower-bound value and upper-bound value of PMRR decrease gradually and have trends of nonlinear reduction. From the first time series to the second time series, the estimated true value of PCRR is rapidly reduced from 0% to −27.84%; from the second time series to the third time series, the estimated true value of PCRR is slowly reduced from −27.84% to −37.28%; from the third time series to the fourth time series, the estimated true value of PCRR is reduced from −37.28% to −100% very rapidly.

This analysis result is more in line with engineering practice because the prior information and current sample information of the vibration performance degradation information of bearings in service were comprehensively considered by the proposed method. In addition, the HMEBM also considers the similarity between time series, which contains hidden information regarding the degradation process of optimal vibration performance state.

4. Conclusions

Based on the HMEBM, the reliability model is established to evaluate the operation performance status of rolling bearings, which has no requirements for the priori information of data samples, types, and components of bearings.

1.

Considering the operation status of rolling bearing, the variation degree of the optimal vibration performance status can be calculated more accurately to ensure effective maintenance of the system, reduce faults, and improve quality under the condition of traditional probability statistics.

2.

The similarities between time series are obtained using the overlapping area method, membership degree method, Hamming approach degree method, Euclidean approach degree method, and cardinal approach degree method.

3.

The maximum entropy method, Bayesian theory, and bootstrap method are fused fully to discover more information in time series of bearing vibration performance. The estimated true values and maximum entropy estimated intervals of PMR and PMRR are calculated to dynamically monitor the health status of rolling bearings online.

The performance degradation law of rolling bearing will be further studied under a different load and speed. In order to verify the universality of the proposed method, HMEBM will be used to evaluate the performance degradation process of other type of bearings.