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

**Liang Ye <sup>1</sup> , Yusheng Hu 2,3, \*, Sier Deng 1 , Wenhu Zhang 1 , Yongcun Cui <sup>1</sup> and Jia Xu 2,3**

> <sup>1</sup> School of Mechatronics Engineering, Henan University of Science and Technology, Luoyang 471003, China; 9945070@haust.edu.cn (L.Y.); dse@haust.edu.cn (S.D.); zwh@haust.edu.cn (W.Z.); 9945020@haust.edu.cn (Y.C.)

<sup>2</sup> State Key Laboratory of Air-Conditioning Equipment and System Energy Conservation, Gree Electric Appliances Co., Ltd., Zhuhai 519070, China; jidianyihao@163.com


**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.

**Keywords:** rolling bearing; performance degradation; variation degree; probability density function; similarities between time 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–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

**Citation:** Ye, L.; Hu, Y.; Deng, S.; Zhang, W.; Cui, Y.; Xu, J. A Novel Model for Evaluating the Operation Performance Status of Rolling Bearings Based on Hierarchical Maximum Entropy Bayesian Method. *Lubricants* **2022**, *10*, 97. https:// doi.org/10.3390/lubricants10050097

Received: 10 April 2022 Accepted: 5 May 2022 Published: 13 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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–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–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–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–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.

**Figure 1.** Flow diagram of proposed method.

#### **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 *Xw*.

$$\mathbf{X}\_{\overline{w}} = (\mathbf{x}\_{\overline{w}}(1), \mathbf{x}\_{\overline{w}}(2), \dots, \mathbf{x}\_{\overline{w}}(k), \dots, \mathbf{x}\_{\overline{w}}(N) \text{ )}; \mathbf{k} = 1, 2, \dots, N; \mathbf{w} = 1, 2, \dots, r \tag{1}$$

where *w* stands for the order number of time series *Xw*; *r* is the number of time series; *xw*(*k*) is the *k*th performance data in time series *Xw*; *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*1.

#### *2.1. Solving PDF*

Transmuting the performance data into continuous information, the PDF *fw*(*x*) of the time series *X<sup>w</sup>* with maximum entropy is defined as

$$f\_w(\mathbf{x}) = \exp(c\_{0w} + \sum\_{i=1}^{j} c\_{iw} \mathbf{x}^i) \tag{2}$$

where *c*0*<sup>w</sup>* is the first Lagrange multiplier and *c*i*<sup>w</sup>* is the (*i* + 1)th Lagrange multiplier for the time series *Xw*; *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(\mathbf{x}) = -\int\_{\tilde{\Omega}\_w} f\_w(\mathbf{x}) \ln f\_w(\mathbf{x}) d\mathbf{x} \to \max \tag{3}$$

where *Hw*(*x*) is the information entropy and Ω*<sup>w</sup>* represents the feasible domain for the data sample of the time series *Xw*, and Ω*<sup>w</sup>* = [*x*min*w*, *x*max*w*]; *x*min*<sup>w</sup>* and *x*max*<sup>w</sup>* are the lower-bound value and upper-bound value in the time series *Xw*; ln*fw*(*x*) is the logarithmic value of *fw*(*x*).

Equation (3) satisfies the constraint conditions

$$\int\_{\Omega^{w}} f\_{w}(\mathbf{x})d\mathbf{x} = 1 \tag{4}$$

$$\int\_{\Omega\_w} \mathbf{x}^i f\_w(\mathbf{x}) d\mathbf{x} = m\_{iw} \tag{5}$$

where *miw* stands for the *i*th order origin moment for the data sample of the time series *Xw*, and *m*0*<sup>w</sup>* = 1.

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

The Lagrange function *Lw*(*x*) can be expressed as

$$L\_{\mathfrak{w}}(\mathbf{x}) = H\_{\mathfrak{w}}(\mathbf{x}) + (c\_{0w} + 1)[\int\_{\Omega\_{\mathfrak{w}}} f\_{\mathfrak{w}}(\mathbf{x}) \mathbf{dx} - m\_{0\mathfrak{w}}] + \sum\_{i=1}^{j} c\_{i\mathfrak{w}} [\int\_{\Omega\_{\mathfrak{w}}} \mathbf{x}^{\mathfrak{i}} f\_{\mathfrak{w}}(\mathbf{x}) \mathbf{dx} - m\_{l\mathfrak{w}}] \tag{6}$$

The first Lagrange multiplier can be given by

$$\mathcal{L}\_{0w} = -\ln\left[\int \exp\{\sum\_{i=1}^{j} c\_{iw} \mathbf{x}^{i}\} \mathbf{dx}\right] \tag{7}$$

Other *j* Lagrange multipliers should satisfy the constraint condition

$$1 - \frac{\int \mathbf{x}^i \exp(\sum\_{i=1}^j c\_{iw}\mathbf{x}^i) \mathbf{dx}}{m\_{iw} \int \exp(\sum\_{i=1}^j c\_{iw}\mathbf{x}^i) \mathbf{dx}} = 0 \tag{8}$$

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

$$\mathbf{x} = a\_w \mathbf{t} + b\_w \tag{9}$$

where *a<sup>w</sup>* and *b<sup>w</sup>* are mapping parameters for the data sample of the time series *Xw*; *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\_{\boldsymbol{w}}(t) = \exp[c\_{0\boldsymbol{w}} + \sum\_{i=1}^{j} c\_{i\boldsymbol{w}} (a\_{\boldsymbol{w}}t + b\_{\boldsymbol{w}})^i] \tag{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 - \mathfrak{a}) \times 100\% \tag{11}$$

Assume that the maximum entropy estimated interval is [*x*L*w*, *x*U*w*] under confidence level *P*. The lower boundary value *x*L*<sup>w</sup>* should satisfy

$$\frac{1}{2}\mathfrak{a} = \int\_{\chi\_{\text{min}}}^{\chi\_{\text{L}w}} f\_w(\mathfrak{x})d\mathfrak{x} \tag{12}$$

The upper boundary value *x*U*<sup>w</sup>* should satisfy

$$1 - \frac{1}{2}\mathfrak{a} = \int\_{\mathcal{X}\_{\text{min}w}}^{\mathcal{X}\_{\text{U}w}} f\_w(\mathfrak{x})d\mathfrak{x} \tag{13}$$

#### *2.3. Calculating PPDF*

Consider the data sample in the intrinsic series *X*<sup>1</sup> 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

$$hyf\_w(\mathbf{x}) = \frac{f\_1(\mathbf{x})f\_w(\mathbf{x})}{\int\_{\Omega\_{1w}} f\_1(\mathbf{x})f\_w(\mathbf{x})d\mathbf{x}} \tag{14}$$

where *hyfw*(*x*) is the PPDF of the *w*th time series *Xw*; Ω1*<sup>w</sup>* is the intersection of the feasible regions of the data samples for the intrinsic series *X*<sup>1</sup> and the *w*th time series *Xw*.

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*<sup>1</sup> = 1. The PMR *R<sup>w</sup>* (1) for the *w*th time series is expressed using the overlapping area method by

$$R\_{w}(1) = S\_{1w} = \int\_{x=x\_{1L}}^{x\_{1w}} hyf\_{1}(x)dx + \eta\_{w} \int\_{x=x\_{1w}}^{x\_{2w}} hyf\_{w}(x)dx + \int\_{x=x\_{2w}}^{x\_{1U}} hyf\_{1}(x)dx\tag{15}$$

where *S*1*<sup>w</sup>* is the overlapping area for the PPDF of the *w*th 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*1*<sup>w</sup>* and *x*2*<sup>w</sup>* are the abscissa values of the intersections for the PPDF of the *w*th time series and the intrinsic series, and *x*1*<sup>w</sup>* < *x*2*w*; *η<sup>w</sup>* is the overlapping area for the PDF of the *w*th time series and the intrinsic series.

$$\eta\_{\rm av} = \int\_{x=\mathbf{x}^\* \atop \mathbf{x} = \mathbf{x}^\* \mathbf{1}}^{\mathbf{x}^\*} f\_1(\mathbf{x}) \mathbf{dx} + \int\_{\mathbf{x} = \mathbf{x}^\* \atop \mathbf{x} = \mathbf{x}^\* \mathbf{1}}^{\mathbf{x}^\*} f\_{\rm av}(\mathbf{x}) \mathbf{dx} + \int\_{\mathbf{x} = \mathbf{x}^\* \atop \mathbf{x} = \mathbf{x}^\* \atop \mathbf{x} \neq \mathbf{x}}^{\mathbf{x}^\*} f\_1(\mathbf{x}) \mathbf{dx} \tag{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 \** <sup>1</sup>*<sup>w</sup>* and *x \** <sup>2</sup>*<sup>w</sup>* are the abscissa values of the intersections for the PDF of the *w*th time series and the intrinsic series, and *x \** <sup>1</sup>*<sup>w</sup>* < *x \** <sup>2</sup>*w*.

#### *2.5. Membership Degree Method*

*X<sup>w</sup>* is mapped into [0,1] by mapping formula to ensure that all data are fuzzy numbers.

$$z\_w(k) = \frac{\varkappa\_w(k) - \varkappa\_{w \text{min}}}{\varkappa\_{w \text{max}} - \varkappa\_{w \text{min}}} \tag{17}$$

Let *Z<sup>w</sup>* 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<sup>w</sup>* is described as

$$\mathbf{Z}\_w = (z\_w(1), z\_w(2), \dots, z\_w(k), \dots, z\_w(N)) \tag{18}$$

When studying the conformity degree of *Z<sup>w</sup>* (*w* = 1,2, . . . ,*r*) relative to *Z*1, define the absolute difference as

$$
\Delta\_w(k) = |z\_w(k) - z\_1(k)|\tag{19}
$$

$$
\Delta k\_{\text{max}} = \max\_{w} \Delta\_{w}(k) \tag{20}
$$

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

$$
\mu\_{wk} = \mu\_{wk}(z\_w(k), z\_1(k)) = 1 - \frac{\Delta\_w(k)}{\Delta k\_{\text{max}}} \tag{21}
$$

The average membership degree is given as

$$\mu\_w = \frac{1}{N} \sum\_{k=1}^{N} \mu\_{wk}, w = 1, 2, \dots, r \tag{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<sup>w</sup>* and *Z*<sup>1</sup> is; and on the contrary, the less significant the relationship between *Z<sup>w</sup>* and *Z*<sup>1</sup> is.

The PMR *Rw*(2) for the *w*th time series is expressed using the membership degree method by

$$R\_w(\mathbf{2}) = \int\_{\mathbf{x} = \mathbf{x}\_{\text{IL}}}^{\mathbf{x}\_{\text{1w}}} hyf\_1(\mathbf{x}) \mathbf{dx} + \mu\_w \int\_{\mathbf{x} = \mathbf{x}\_{\text{1w}}}^{\mathbf{x}\_{\text{2w}}} hyf\_w(\mathbf{x}) \mathbf{dx} + \int\_{\mathbf{x} = \mathbf{x}\_{2w}}^{\mathbf{x}\_{\text{1U}}} hyf\_1(\mathbf{x}) \mathbf{dx} \tag{23}$$

#### *2.6. Approach Degree Method*

Let *X*<sup>1</sup> and *X<sup>w</sup>* be fuzzy subsets on the finite field *Q*. µ1k∈[0,1] and µwk∈[0,1] are the membership degrees of *X*<sup>1</sup> and *Xw*, respectively. Minkowski distance *dpk* is defined as

$$d\_{pk} = d\_p(\mathbf{X}\_1 \mathbf{X}\_w) = \left(\frac{1}{N} \sum\_{k=1}^{N} |\mu\_{1k} - \mu\_{wk}|^p\right)^{\frac{1}{p}} \tag{24}$$

where *N* is the number of elements in sets *X*<sup>1</sup> and *Xw*; *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*<sup>0</sup> and *X<sup>i</sup>* is defined as

$$\beta = \beta\_p(\mathbf{X}\_1, \mathbf{X}\_w) = \frac{1}{2} \left[ \mathbf{X}\_1 \circ \mathbf{X}\_w + \left( 1 - \mathbf{X}\_1 \xrightarrow{\frown} \mathbf{X}\_w \right) \right] \tag{25}$$

where

$$\mathbf{X}\_1 \circ \mathbf{X}\_w = \bigvee\_{k=1}^N (\mu\_{1k} \wedge \mu\_{wk}) \tag{26}$$

$$\mathbf{X}\_1 \stackrel{\frown}{\circ} \mathbf{X}\_w = \bigwedge\_{k=1}^N (\mu\_{1k} \lor \mu\_{wk}) \tag{27}$$

The Hamming approach degree *β*1*<sup>w</sup>* is given as

$$\beta\_{1w} = \beta\_1(\mathbf{X}\_1, \mathbf{X}\_w) = 1 - \frac{1}{N} \sum\_{k=1}^{N} |\mu\_{1k} - \mu\_{wk}| \tag{28}$$

The PMR *Rw*(3) for the *w*th time series is expressed using the Hamming approach degree method by

$$\mathcal{R}\_{\mathbf{w}}(\mathfrak{Z}) = \int\_{\mathbf{x} = \mathbf{x}\_{\mathrm{IL}}}^{\mathbf{x}\_{\mathrm{1w}}} h y f\_{1}(\mathbf{x}) \mathbf{dx} + \beta\_{\mathbf{1}w} \int\_{\mathbf{x} = \mathbf{x}\_{\mathrm{1w}}}^{\mathbf{x}\_{\mathrm{2w}}} h y f\_{\mathbf{w}}(\mathbf{x}) \mathbf{dx} + \int\_{\mathbf{x} = \mathbf{x}\_{\mathrm{2w}}}^{\mathbf{x}\_{\mathrm{1U}}} h y f\_{1}(\mathbf{x}) \mathbf{dx} \tag{29}$$

Euclidean approach degree *β*2*<sup>w</sup>* is defined as

$$\beta\_{2w} = \beta\_2(\mathbf{X}\_1, \mathbf{X}\_w) = 1 - \left(\frac{1}{N} \sum\_{k=1}^{N} |\mu\_{1k} - \mu\_{wk}|\right)^{1/2} \tag{30}$$

The PMR *Rw*(4) for the *w*th time series is expressed using the Euclidean approach degree method by

$$R\_w(4) = \int\_{x=x\_{1L}}^{x\_{1w}} hyf\_1(x)dx + \beta\_{2w} \int\_{x=x\_{1w}}^{x\_{2w}} hyf\_w(x)dx + \int\_{x=x\_{2w}}^{x\_{1U}} hyf\_1(x)dx\tag{31}$$

Cardinal approach degree *β*3*<sup>w</sup>* is given as

$$\beta\_{3w} = \beta\_3(\mathbf{X}\_1, \mathbf{X}\_w) = \frac{2\sum\_{k=1}^{N} (\mu\_{1k} \wedge \mu\_{wk})}{\sum\_{k=1}^{N} \mu\_{1k} + \sum\_{k=1}^{N} \mu\_{wk}} \tag{32}$$

The PMR *Rw*(5) for the *w*th time series is expressed using the Cardinal approach degree method by

$$R\_w(\mathbf{5}) = \int\_{x=\mathbf{x}\_{\mathrm{IL}}}^{\mathbf{x}\_{\mathrm{1w}}} hyf\_1(\mathbf{x}) \mathrm{d}x + \beta\_{\mathrm{3w}} \int\_{x=\mathbf{x}\_{\mathrm{1w}}}^{\mathbf{x}\_{\mathrm{2w}}} hyf\_w(\mathbf{x}) \mathrm{d}x + \int\_{x=\mathbf{x}\_{\mathrm{2w}}}^{\mathbf{x}\_{\mathrm{1U}}} hyf\_1(\mathbf{x}) \mathrm{d}x \tag{33}$$

The approach degree describes the similarity degree between *X*<sup>1</sup> and *Xw*. The larger the approach degree is, the more significant the relationship between *X*<sup>1</sup> and *X<sup>w</sup>* is; on the contrary, the less significant the relationship between *X*<sup>1</sup> and *X<sup>w</sup>* is.

#### *2.7. Dynamic Evaluation of PMR*

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

$$R\_w = (R\_w(1), R\_w(2), \dots, R\_w(5)) = (R\_w(\gamma)); \gamma = 1, 2, \dots, 5; w = 1, 2, \dots, r; \tag{34}$$

where *R<sup>w</sup>* is the data sample of PMR for the *w*th time series; *Rw*(*γ*) is the *γ*th data in the PMR data sample for the *w*th time series.

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

$$R\_{w\text{Boosttrap}} = (R\_{w1}, R\_{w2}, \dots, R\_{w\theta\prime}, \dots, R\_{wB}) \tag{35}$$

where *Rw<sup>θ</sup>* 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\Theta} = [R\_{w\Theta}(\Theta)]; \Theta = 1, 2, \dots, z \tag{36}$$

where *Rw<sup>θ</sup>* (Θ) is the Θth data in the *θ*th bootstrap re-sampling sample of PMR for the *w*th time series.

The maximum entropy method is used to calculate the PDF of the generated sample *Rw*bootstrap. 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(\mathcal{R}\_w) = \exp[\mathcal{c}^\*\_{0w} + \sum\_{i=1}^j \mathcal{c}^\*\_{iw} (a^\* \!\_w \mathcal{R}\_w + b^\* \!\_w)^i] \tag{37}$$

where *c \** <sup>0</sup>*<sup>w</sup>* is the first Lagrange multiplier and *c \* iw* is the (*i* + 1)th Lagrange multiplier for the time series *Rw*; *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 *Rw*.

The estimated true value of the PMR is obtained as

$$R\_{w0} = \int\_{\mathcal{S}} R\_w f(R\_w) \,\mathrm{d}R\_w \tag{38}$$

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

$$\begin{bmatrix} \left[ R\_{w\mathbf{L}\prime} R\_{w\mathbf{U}} \right] \end{bmatrix} = \left[ R\_{w\frac{\mathbf{g}}{2}\prime} R\_{w\_{1-\frac{\mathbf{g}}{2}}} \right] \tag{39}$$

with

$$\frac{\alpha}{2} = \int\_{R\_{\text{wmin}}}^{R\_{w\_{\frac{\alpha}{2}}}} f(R\_w) \, \text{d}R\_w \tag{40}$$

$$1 - \frac{a}{2} = \int\_{R\_{w\text{min}}}^{R\_{w\_1 - \frac{a}{2}}} f(R\_w) \mathbf{d}R\_w \tag{41}$$

where *Rw*<sup>L</sup> is the lower bound value and *Rw*<sup>U</sup> is the upper bound value of the PMR for the *w*th time series.

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

$$P\_w(\gamma) = 1 - R\_w(\gamma) \tag{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(\gamma) = \frac{R\_w(\gamma) - R\_{10}}{R\_{10}} \times 100\% \tag{43}$$

where *R*<sup>1</sup> is the PMR for the intrinsic series of rolling bearings.

*dw*(*γ*) < 0 indicates that the PMR for the time interval corresponding with the *w*th time series is less than the PMR for the time interval corresponding with the intrinsic series. The smaller *dw*(*γ*) 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 *w*th time series.

$$d\_{w0} = \frac{R\_{w0} - R\_{10}}{R\_{10}} \times 100\% \tag{44}$$

$$d\_{\rm wL} = \frac{R\_{\rm wL} - R\_{10}}{R\_{10}} \times 100\% \tag{45}$$

$$d\_{w\mathbf{U}} = \frac{R\_{w\mathbf{U}} - R\_{10}}{R\_{10}} \times 100\% \tag{46}$$

where *dw*<sup>0</sup> is the estimated true value, *dw*<sup>L</sup> is the lower-bound value and *dw*<sup>U</sup> is the upperbound value of the PMRR for the *w*th 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 Figures 2 and 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<sup>3</sup> , elastic modulus 2.08 <sup>×</sup> <sup>10</sup><sup>5</sup> N/mm<sup>2</sup> , Poisson's ratio 0.3 and Hardness 700 HV10. The bearing parameters are shown in Table 1.

**Figure 2.** Physical drawing of testing machine.

**Figure 3.** Testing head. (**a**) Physical drawing, (**b**) Main view, (**c**) Section view.


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

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 −2 . 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 −2 , 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.

**Figure 4.** Vibration signals of bearing.

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 −2 . 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 −2 . 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 −2 . 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 −2 . 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 −2 . 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 −2 . 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*11, *m*21, *m*31, *m*41, *m*51] = [−0.7710, 1.1829, −1.6809, 3.0058, −5.0425]; the Lagrange multipliers [*c*01, *c*11, *c*21, *c*31, *c*41, *c*51] = [−0.526176, −1.0192, −0.8837, −0.3003, 0.0325, 0.0628]; the mapping parameters *a*<sup>1</sup> = 1.8819 and *b*<sup>1</sup> = −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 −2 for the first time series. According to Equation (10), the probability density estimated function *f* <sup>1</sup>(*x*) is calculated as shown in Figure 5.

**Figure 5.** PDF of data sample of the intrinsic sequence.

For the second time series, various order origin moments can be obtained using the maximum entropy method as [*m*12, *m*22, *m*32, *m*42, *m*52] = [−0.4015, 0.8709, −0.7686, 1.6745, −1.9542]; the Lagrange multipliers [*c*02, *c*12, *c*22, *c*32, *c*42, *c*52] = [−0.3629, −0.3569, −0.2062, −0.1813, −0.1369, 0.0291]; the mapping parameters *a*<sup>2</sup> = 1.9572 and *b*<sup>2</sup> = −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 −2 for the second time series.

For the third time series, various order origin moments can be obtained using the maximum entropy method as [*m*13, *m*23, *m*33, *m*43, *m*53] = [−0.8840, 1.5477, −2.2030, 4.5372, −7.4059]; the Lagrange multipliers [*c*03, *c*13, *c*23, *c*33, *c*43, *c*53] = [−1.4292, −1.1189, −0.2384, 0.0274, −0.0268, 0.0145]; the mapping parameters *a*<sup>3</sup> = 1.1379 and *b*<sup>3</sup> = −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 −2 for the third time series.

For the fourth time series, various order origin moments can be obtained using the maximum entropy method as [*m*14, *m*24, *m*34, *m*44, *m*54] = [−0.1634, 1.2181, −0.2573, 3.2716, −0.7843]; the Lagrange multipliers [*c*04, *c*14, *c*24, *c*34, *c*44, *c*54] = [−2.3395, −0.5070, −0.0421, 0.2472, −0.0660, −0.0299]; the mapping parameters *a*<sup>4</sup> = 0.3328 and *b*<sup>4</sup> = −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 −2 for the fourth time series.

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

**Figure 6.** PDF of data samples of time series.

#### 3.1.2. PPDF of Data Samples of Time Series

Intersection interval [*xw*1, *xw*2] and intersection area *η<sup>w</sup>* are calculated for the PDF of the *w*th time series and the intrinsic series, as shown in Table 2. The similarity degrees between time series *ηw*, *µw*, *β*1*w*, *β*2*<sup>w</sup>* and *β*3*<sup>w</sup>* 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.

**Table 2.** Intersection intervals and overlapped areas of PDF.


**Table 3.** Similarity degrees calculated using different methods.


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* <sup>1</sup>(*x*)*fw*(*x*) and posterior probability density functions *hyfw*(*x*) of the four time series are constructed using the Bayesian principle, as shown in Figures 7 and 8.

**Figure 7.** Product *f* <sup>1</sup> (x)*fw*(x) of PDF.

**Figure 8.** PPDF of 4 time series.

#### 3.1.3. PMR of Time Series

The intersection interval [*x*1*w*, *x*2*w*] and intersection area *S<sup>w</sup>* are calculated for the PPDF of the *w*th time series and the intrinsic series, as shown in Table 4. The values of PMR, *Rw*(1), *Rw*(2), *Rw*(3), *Rw*(4), and *Rw*(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. The PMR of each time series decreases gradually, and the slope of the curve decreases gradually, as shown in Figure 9.




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

**Figure 9.** PMR of time series.

#### 3.1.4. PMRR of Time Series

The values of PMRR, *dw*(1), *dw*(2), *dw*(3), *dw*(4), and *dw*(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.

**Table 6.** PMRR calculated using different methods.


#### 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*<sup>2</sup> = (*R*2(1), *R*2(2), *R*2(3), *R*2(4), *R*2(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.

**Figure 10.** Generated data of PMR of second time series.

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*<sup>20</sup> = 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.

**Figure 11.** PDF of PMR of the second time series.



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, lowerbound 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.

**Figure 12.** Maximum entropy estimated results of PMR.

The PMRR for vibration performance of rolling bearings are calculated for different time series, as shown in Table 8 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.

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

**Figure 13.** Maximum entropy estimated results of PMRR.

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.

**Figure 14.** Vibration signals of bearing (Case 2).

#### 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*11, *m*21, *m*31, *m*41, *m*51] = [−0.0463, 0.8915, −0.2456, 2.1047, −0.7743]; the Lagrange multipliers [*c*01, *c*11, *c*21, *c*31, *c*41, *c*51] = [0.6246, 0.3796, −0.4400, −0.3470, −0.0294, 0.0463]; the mapping parameters *a*<sup>1</sup> = 4.4481 and *b*<sup>1</sup> = −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 −2 for the first time series.

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

**Figure 15.** PDF of data samples of time series (Case 2).

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

The intersection interval [*xw*1, *xw*2] and intersection area *η<sup>w</sup>* are calculated for the PDF of the *w*th time series and the intrinsic series as shown in Table 9. The similarity degrees between time series *ηw*, *µw*, *β*1*w*, *β*2*w,* and *β*3*<sup>w</sup>* 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.

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



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

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* <sup>1</sup>(*x*)*fw*(*x*) and posterior probability density functions *hyfw*(*x*) of the four time series are constructed using the Bayesian principle, as shown in Figures 16 and 17.

**Figure 16.** Product *f* <sup>1</sup> (x)*fw*(x) of PDF (Case 2).

**Figure 17.** PPDF of 4 time series (Case 2).

#### 3.2.3. PMR of Time Series (Case 2)

The intersection interval [*x*1*w*, *x*2*w*] and intersection area *S<sup>w</sup>* are calculated for the PPDF of the *w*th time series and the intrinsic series, as shown in Table 11. The values of PMR, *Rw*(1), *Rw*(2), *Rw*(3), *Rw*(4), and *Rw*(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. The PMR of each time series decreases gradually, and the slope of the curve decreases gradually, as shown in Figure 18.

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


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


**Figure 18.** PMR of time series (Case 2).

#### 3.2.4. PMRR of Time Series (Case 2)

The values of PMRR, *dw*(1), *dw*(2), *dw*(3), *dw*(4), and *dw*(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.


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

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*<sup>2</sup> = (*R*2(1), *R*2(2), *R*2(3), *R*2(4), *R*2(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.

**Figure 19.** Generated data of PMR of second time series (Case 2).

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*<sup>20</sup> = 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.

**Figure 20.** PDF of PMR of the second time series (Case 2).


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

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, lowerbound 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.

**Figure 21.** Maximum entropy estimated results of PMR (Case 2).

The PMRR for vibration performance of rolling bearings are calculated for different time series, as shown in Table 15 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.

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


**Figure 22.** Maximum entropy estimated results of PMRR (Case 2).

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.


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.

**Author Contributions:** Conceptualization, L.Y.; methodology, L.Y.; validation, L.Y. and Y.H.; formal analysis, W.Z.; investigation, L.Y.; resources, L.Y.; data curation, L.Y.; writing—original draft preparation, L.Y.; writing—review and editing, W.Z.; visualization, Y.C.; supervision, J.X.; project administration, S.D.; funding acquisition, Y.H., W.Z. and Y.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [Opening Foundation of State Key Laboratory of Air-Conditioning Equipment and System Energy Conservation] grant number [ACSKL2019KT02], [National Natural Science Foundation of China] grant number [U1804145], And The APC was funded by [Yusheng Hu], [Youth Programs of the National Natural Science Foundation of China] grant number [51905152 and 52005158], And The APC was funded by [Yusheng Hu, Xiaoqiang Wang, Wenhu Zhang and Yongcun Cui].

**Data Availability Statement:** The data used to support the findings of this study are available from the corresponding author upon request.

**Conflicts of Interest:** The authors declare that there is no conflict of interest regarding the publication of this paper.

#### **Nomenclature**


#### **References**

