Fuzzy reliability prediction of rotating machinery product with accelerated testing data

3.1. Feature extraction and principal component analysis

Feature extraction is to understand the performance of machinery product from the vibration signal $s\left(t\right)$ at different aspect, which mainly includes time-domain, frequency-domain, time-frequency domain and information entropy as follows [2, 3].

3.2. Acquiring degradation index through GMM method

Based on the new $k$ dimensional feature set, the next step is to acquire the degradation index. For bearings, a gradual deterioration may happen due to ball wearing or inner impurity particles, which means that the extracted features should present the changing process. For instance, during the first $m$ monitor intervals, the tested product is in health status. With the time goes on, the most recent $m$ intervals should be in a partially degradation status. Following this idea, a GMM modeling method is introduced to obtain the degradation index through the measure of overlap rate [10, 20, 21].

Supposed that $\mathbf{y}=\left[y_1,y_2,\dots ,y_m\right]$ present one status with the $m$ columns of feature set $\mathbf{Y}$. The probability density function (PDF) of$\mathrm{}{y}_q$, $q=\text{1,}\text{2,…,}m$, follows a mixture distribution with $N$ Gaussian components that is:

where, ${\omega }_p$ is the weight of the $p$th mixture and $Norm\left({\mathbf{\mu }}_p,{\mathbf{\Sigma }}_p\right)$ is:

Through the expectation maximization (EM) algorithm, the values of unknown parameters, i.e.$\mathrm{}{\omega }_p$, ${\mu }_p$, ${\mathrm{\Sigma }}_p$, will be obtained when the likelihood function value $L\left(\mathbf{y}|\mathbf{\mu },\mathbf{\Sigma }\right)={\prod }_{q=1}^{m}f\left(y_q\right)$ is converged or the number of iteration reaches the maximum $i_{\mathrm{m}\mathrm{a}\mathrm{x}}$ [10, 22].

3.3. Basic EM algorithm for the estimation of GMM model parameters

1) $i\leftarrow 0$, initialize the unknown parameters for all Gaussian components, ${\omega }_{p\left(0\right)}$, ${\mu }_{p\left(0\right)}$, ${\mathrm{\Sigma }}_{p\left(0\right)}$.

2) do $i\leftarrow i+1$.

$E$ step: compute the posterior probability of sample of $y_q$ in the $p$th component:

$M$ step: update ${\omega }_p$, ${\mu }_p$ and ${\mathrm{\Sigma }}_p$:

3) end if $i=i_{\mathrm{m}\mathrm{a}\mathrm{x}}$ or $\left|\left(L_i-L_{\left(i-1\right)}\right)/L_{\left(i-1\right)}\right|<\epsilon$.

Therefore, the GMM model for each $m$-interval (i.e. one health state) can be computed through abovementioned algorithm. The confidence value (CV) between the baseline model (i.e. normal state) $f_0$ and the most recent $h$th model $f_h$ can then be calculated by the overlap between two GMM models [10, 23]:

If the $h$th state is still in normal, the CV value will be near 1. Otherwise, it will be near 0 if the $h$th state is barely normal which indicates the corresponding maintenance should be given to avoid the system failure. Hence, the CV results can be used as the degradation index for machinery product health assessment.

4.1. Basic fuzzy arithmetic

Zadeh [24] introduced fuzzy set theory to solve the problem of uncertainty due to imprecise or vague through applying a membership function to data rather than crisp one. Membership functions, e.g. triangular, rectangular or trapezoidal, are widely used. Without loss of generality, triangular membership function is selected to simply the calculations. Let $\tilde{X}=tfn\left(\stackrel{-}{x},l,r\right)$ be a fuzzy number with triangular membership function, see Fig. 2. Herein, $\stackrel{-}{x}$, $l$ and $r$ donate the center, left boundary and right boundary, respectively. The ${\mu }_{\tilde{X}}\left(x\right)\in \left[0,1\right]$ presents the membership value of $x$ in$\mathrm{}\tilde{X}$, i.e.:

Fig. 2Fuzzy number and its alpha-cuts

The realization of fuzzy numbers are always conducted with alpha-cuts, which decomposes the membership function at the vertical axis with membership value $\alpha$ (see Fig. 2), i.e. $0\le \alpha \le 1$. The alpha-cuts of $\tilde{X}$ is defined as:

where, $c\left(\alpha \right)$ and $d\left(\alpha \right)$ are increasing and decreasing function of $\alpha$ with $c\left(\alpha \right)\le d\left(\alpha \right)$ [25]. When $\alpha =0$, [$c\left(0\right),d\left(0\right)]$ is called the support of fuzzy number $\tilde{X}$.

Assuming two triangular numbers $\tilde{A}$ and $\tilde{B}$ with their alpha-cuts, namely $\tilde{A}\left[\alpha \right]=\left[c_1\left(\alpha \right),d_1\left(\alpha \right)\right]$ and $\tilde{B}\left[\alpha \right]=\left[c_2\left(\alpha \right),d_2\left(\alpha \right)\right]$, the operator can be calculated with basic interval arithmetic and alpha-cuts. The results of addition, subtraction, multiplication and division between two fuzzy numbers are $\tilde{C}\left[\alpha \right]=\tilde{A}\left[\alpha \right]+\tilde{B}\left[\alpha \right]$, $\tilde{A}\left[\alpha \right]-\tilde{B}\left[\alpha \right]$, $\tilde{A}\left[\alpha \right]\tilde{B}\left[\alpha \right]$ and $\tilde{A}\left[\alpha \right]/\tilde{B}\left[\alpha \right]$, respectively. Provided that zero does not belong to $\tilde{B}\left[\alpha \right]$ for all $\alpha$ in division [18, 25]. Assuming a function is $h:\left[a,b\right]\to \mathcal{R}$, $\tilde{Z}=h\left(\tilde{X}\right)$ is a fuzzy function and the membership function of $\tilde{Z}$ can be calculated through Zadeh’s extension principle, i.e.

If $h$ is continuous, then the alpha-cuts for $\tilde{Z}$, namely $\tilde{Z}\left[\alpha \right]=[z_1\left(\alpha \right),z_2\left(\alpha \right)]$, are [26]:

4.2. Fuzzy regression method and q-TTFs

For the analysis of imprecise degradation data, fuzzy regression method is a powerful tool to model the degradation path containing the epistemic uncertainty by means of fuzzy theory. There are two options when using this method to model the degradation path. The first one utilizes a membership function to each degradation data. Then, the regression method with fuzzy parameters is used to model the degradation path [17, 27]. The second one directly uses the regression method to model the degradation path and select the (1-$\alpha$) 100 % confidence intervals of the parameters as their fuzzy estimators [18]. Taking the CV results in Section 3.2 as an example, the data are imprecise because of the inherent uncertainties. Since no prior information is given to establish the membership function for CV results, the second method is selected in this paper.

In general, the degradation paths can be classified into three types: linear, convex and concave. Thus, the linear, exponential and power law relationships are widely used for degradation modeling. And, the basic regression model is defined as:

where, $D(·)$ is the expected degradation path, $\theta =\left[a,b\right]$ is the constant parameters, $\epsilon$ is the error term which follows normal distribution. Herein, one may use the statistical fitting indexes to select the most appropriate degradation model among the three candidate models in Eq. (12), e.g. $R^2$ and $R_{adj}^2$. Meanwhile, the unknown parameters of the selected degradation model can be estimated through least square method.

Then, the (1-$\alpha$) 100 % confidence interval for vector parameters are regarded as the fuzzy estimators $\tilde{\theta }$ for parameters $\widehat{\theta }$, which is [18]:

where, $n_{cv}$ is the number of CV results, $p$ is the number of unknown parameters which is equal to 2. For the sake of simplicity, the symmetric triangular membership function is selected for fuzzy parameters, i.e. $l=r$ in Eq. (8), and their alpha-cuts are defined as $\tilde{\theta }\left[\alpha \right]$ for all $\alpha \in \left[\text{0, 1}\right]$.

Thus, the fuzzy regression model is:

When the failure threshold $C$ is available, the corresponding fuzzy failure time at each operating condition $\tilde{T}\left[\alpha \right]=\left[\left(t_1|C\right)\left(\alpha \right),\left(t_2|C\right)\left(\alpha \right)\right]$ will be computed at each alpha-cut using interval arithmetic.

4.3. Acceleration model

The machinery products are experienced with hasher than normal conditions. In order to evaluate the lifetime and reliability results at normal operating condition, the relationship between q-TTFs and accelerated conditions should be obtained through acceleration model.

Generally, acceleration models are given through physical mechanism analysis or empirical experiences. Examples include Arrhenius model, Eyring model or inverse power model [28]. When no prior knowledge is acquired for machinery products, like the bearing under different load and speed conditions, the following log-linear relationship is recommended. Similar application can be found in [29]:

where, ${\lambda }_k$ is the coefficient, $k=$ 1, 2,…,5, $S$ and $L$ donate speed and load for bearing operating conditions, $T$ is the failure time. The optimal parameters for Eq. (15) can be obtained through least square method by “fitlm” function in Matlab. Similarly, the fuzzy acceleration model can be given to describe the relationship between the environmental variables, i.e. load, speed, and the fuzzy q-TTF $\tilde{T}\left[\alpha \right]$ under each alpha-cut for all $\alpha \in \left[\text{0, 1}\right]$.

Then, under each alpha-cut, the averaged fuzzy q-TTFs at both normal operating condition and the $j$th stress level, i.e. ${\widetilde{\stackrel{-}{T}}}_0\left[\alpha \right]$ and ${\widetilde{\stackrel{-}{T}}}_j\left[\alpha \right]$, can be computed by substituting the corresponding speed and load into the optimal function. Therefore, the fuzzy acceleration factor (AF) is defined as:

With the results from Eq. (16), fuzzy q-TTFs $\tilde{T}\left[\alpha \right]$ can be extrapolated to normal operating condition ${\tilde{T}}_0\left[\alpha \right]$.

Besides the lifetime and reliability evaluation, another advantage of accelerated testing is that it can provide information about the stresses or their combinations which are sensitive to the failure of machinery products. Such information can provide reliable basis for reliability growth and design improvement for products. The following algorithm is proposed for such purpose through the analysis of fuzzy coefficients ${\tilde{\lambda }}_k\left[\alpha \right]$.

4.4. Analysis of sensitive stresses

Set $i_{\mathrm{m}\mathrm{a}\mathrm{x}}=M$, which denotes the number of intervals for membership value in [0, 1].

1) $i\leftarrow 0$, compute the fuzzy coefficient at each alpha-cut $\alpha \leftarrow i/M$.

Substitute all fuzzy failure time ${\tilde{T}}_j\left[\alpha \right]$ and its corresponding $j$th condition $S_j$ and $L_j$:

$model\leftarrow fitlm\left({\tilde{T}}_j\left[\alpha \right]{,}^{\text{'}}T~1+S+L+S*L+S^2+L^2\text{'}\right)$;

${\lambda }_k\left[\alpha \right]\leftarrow model.Coefficients.Estimate\left(k\right)$;

${\widetilde{\stackrel{-}{T}}}_0\left[\alpha \right]\leftarrow predict\left(model,\left[S_0,L_0\right]\right)$;

${\widetilde{\stackrel{-}{T}}}_j\left[\alpha \right]\leftarrow predict\left(model,\left[S_j,L_j\right]\right)$;

Compute fuzzy acceleration factor through Eq. (16).

2) End if $i\leftarrow i_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

3) Construct membership function for all fuzzy coefficient ${\tilde{\lambda }}_k\left[\alpha \right]$ and analyze the relationship between them, compute the main stress.

5.1. Lifetime distribution for q-TTFs

In order to select the lifetime distribution under fuzzy q-TTFs, Anderson-Darling test [30] is used to test the goodness of fit of the candidate distribution to the obtained q-TTFs at each alpha-cut. The following procedure is proposed to rank the candidate distributions which are normal, exponential, extreme value, Log-normal and Weibull distributions.

1) Set$\mathrm{}{i}_{\mathrm{m}\mathrm{a}\mathrm{x}}=M$, $h_j=0$ which denotes the total score of the $j$th distribution at all alpha-cuts in [0, 1], and the significance level (default is 5 %).

2) $i\leftarrow \text{0}$, compute the fuzzy coefficient at each alpha-cut $\alpha \leftarrow i/M$.

Calculate the score whether A-D test reject the null hypothesis that failure times ${\tilde{T}}_0\left[\alpha \right]$ are from the $j$th distribution.

$h\leftarrow adtest\left({\tilde{T}}_0\left[\alpha \right]{,}^{\mathrm{\text{'}}}{j}^{th}\mathrm{}distributio{n}^{\mathrm{\text{'}}}\right)$;

If yes,$\mathrm{}h=0$, $h_j=h_j$; else $h=1$, $h_j=h_j+1$.

3) End if $i\leftarrow i_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

4) Rank the $h_j$ in an ascending order and select the lowest (or the first) one as the most appropriate lifetime distribution.

For some specific applications, like bearings, studies have shown that the failure times follow Weibull distribution [31] which can be directly used for fuzzy reliability analysis:

where, $\eta$ is scale parameter, and $\beta$ is shape parameter. In classical reliability engineering, $\beta >$1 means that the product experiences a wear-out phase from the perspective of the bathtub curve. For fuzzy q-TTFs, the fuzzy parameters, i.e. $\tilde{\eta }\left[\alpha \right]$ and $\tilde{\beta }\left[\alpha \right]$, are estimated through lifetime distribution fitting for ${\tilde{T}}_0\left[\alpha \right]$ at each alpha-cut.

5.2. Fuzzy reliability analysis

The selected lifetime distribution $f\left(\bullet \right)$ can then be used for reliability prediction. The fuzzy reliability at time interval [$t_1$, $t_2$] is defined as (Weibull distribution as an example):

Other reliability indexes, e.g. Mean Time to Failure (MTTF) or hazard function, can compute through substituting fuzzy parameters into classical analytical results at each alpha-cut:

Furthermore, if other lifetime distribution is satisfied, e.g. Log-normal distribution, Eqs. (18) and (19) can be rewritten accordingly.

6.1. Experiment platform

The analyzed data is from IEEE PHM 2012 data challenge and the bearing experiment platform PRONOSTA is shown in Fig. 3, which is consist of three parts: the rotating, degradation generation and measurement parts [15]. The sampling frequencies are 25.6 kHz for the embedded acceleration measures and 10 Hz for temperature ones. The platform can perform bearing accelerated degradation testing, while monitoring vibration and temperature signals for fault diagnosis and remaining useful life (RUL) prediction.

Fig. 3The bearing experiment platform [17]

The tested 17 bearings are experienced three operating conditions with different rotating speed and load, i.e. $S=$1800 rpm and $L=$ 4000 N (Condition 1), $S=$1650 rpm and $L=$4200 N (Condition 2) and $S=$1500 rpm and $L=$5000 N (Condition 3) and the bearing failure happens when the accelerometer exceeds 20 g. Table 1 gives the details about the tested bearings which are divided into two sets for training and testing (algorithms or models). Due to the inner non-linear property for bearing features, it is quite difficult to accurately estimate the RULs [15, 32]. Hence, the proposed framework is used to conduct bearing reliability analysis with accelerated data. Since the bearing reliability is an index for the population but not individual, the bearings whose lifetimes are significant lower or higher than others are omitted, i.e. 1_2, 1_4, 2_1, 2_7. In addition, the vibration signals from vertical axis is selected since it has been verified to be better reflect the health state [33].

6.2. Degradation modeling and fuzzy q-TTFs

According to the procedure in Section 2, the original signal is at first de-noised by the wavelet method at level 4 using the principle of Stein’s unbiased risk with soft-thresholding [34]. Then, the $p=$ 25 features in four domains are extracted and then reduced through PCA method where $CP=$ 0.8. After that, the degradation index CV can be obtained using GMM method in Section 3.2. Herein, the choice for the number of features as one state needs more attention, especially for the baseline normal state. We used different combination of $m$ values in three operating conditions and compared their CV results with the announced bearing failure in the learning set, where the acceleration level exceeds 20 g. Thus, $m=$ 1100, 300 and 250 are selected as one health state for three conditions. In addition, the new failure threshold for the degradation index is set to 0.2, i.e. $C=$ 0.2. The GMM results are shown in Fig. 4.

Fig. 4The degradation processes for bearings under three accelerated levels with their fitted paths

Then, three degradation models given in Eq. (12) are used to model the path of CV results at each operating condition. The fitting results are listed in Table 2. The linear model gives bad fitting results for bearing 2_3, 2_5, 2_6 although it performs well for other bearings. Meanwhile, Fig. 4 shows that the degradation processes are non-linear. For the other two models, the exponential model is selected with larger $R^2$ and $R_{adj}^2$ values than power model for most of the bearings, which will be used for the following q-TTF extrapolation.

The statistical results for the selected exponential degradation model are given in Table 3, where Anderson-Darling test is used to check the normal assumption of the error term. Hence, the degradation paths for bearings No. 1_3, 2_3, 2_5 and 3_2 are not compromised with the normal assumption although the $R^2$ and $R_{adj}^2$ suggest good fitting. As seen from Fig. 4, the bearing 1_3 experienced a sudden decrease and it will fail soon which is not accordance with others in operating condition 1. Thus, the bearing 1_3 is omitted for the following analysis. In terms of bearings No. 2_3, 2_5 and 3_2, the CV values are already lower than the failure threshold 0.2 around 2 (h). Considered that the tails of the degradation paths present a state that is much deviate from the normal state. The CV results are ignored from time 1.5, 2 and 2.111(h) from the three bearings. As to bearing 2_2, it seems to perform better after 1.5833 (h) which, however, has already exceeded the threshold. Thus, the CV values after 1.5833 (h) are ignored for bearing 2_2. The new results for this four bearings are given in Table 3 by adding a bracket, e.g. [2_3].

The fuzzy parameters of the exponential degradation model are constructed as illustrated in Section 4.2. The 95 % confidence values are shown at the right of Table 3. For simplicity, the symmetrical triangle membership function is selected to present the belief of certain parameter through the membership values. Thus, the fuzzy q-TTFs can be computed by substituting fuzzy parameters, i.e. $\tilde{a}\left[\alpha \right]$ and $\tilde{b}\left[\alpha \right]$, and the failure threshold $C=$ 0.2 into Eq. (14) at each alpha-cut, which are interval values as $\tilde{T}\left[\alpha \right]=\left[\left(t_1|C\right)\left(\alpha \right),\left(t_2|C\right)\left(\alpha \right)\right]$. The results are given at the left side of Table 4.

6.3. Acceleration model and analysis of sensitive stresses

After obtaining the fuzzy q-TTFs under three operating conditions, Eqs. (15) and (16) in Section 4.3 can then compute the fuzzy acceleration factors (see Fig. 5) with their fuzzy coefficients, and convert the failure times into normal operating conditions. Assumed that the normal condition is 3000 (rpm) in Speed and 1200 (N) in load, the results are given at the right side of Table 4.

From Fig. 5, it can conclude that the support of fuzzy acceleration factors under three acceleration operating conditions are [1.53, 2.56], [4.96, 14.90] and [3.95, 11.12], while the peak values are 2.44, 10.57 and 8.58, respectively. Intuitively, both condition 2 & 3 are harsher than condition 1 with larger values of acceleration factor under each alpha-cut for all $\alpha \in \left[\text{0, 1}\right]$. And condition 2 is severer than condition 3 when the belief of the results $\alpha >\mathrm{}$0.8, otherwise their difference are undistinguishable. Furthermore, the fuzzy coefficients of the acceleration model can be used to analyze the sensitivity of speed and load to the lifetimes of bearings through the procedure in Section 4.3. The results show that$\mathrm{}{\lambda }_0$, ${\lambda }_1$ and$\mathrm{}{\lambda }_3$ are zero for all $\alpha \in \left[\text{0,}\text{1}\right]$ which means that no contribution to the failure of bearings has been obtained for the single stress of speed and the interactive stress of speed and load. Thus, the sensitive stresses are load, self-correlation of speed and load, which are related to coefficient$\mathrm{}{\lambda }_2$, ${\lambda }_4$ and$\mathrm{}{\lambda }_5$, see Fig. 6. The results can help understand the hazard rating for each stress or their combination to the failure of bearings, and improve the design or operating conditions to increase the lifespan of bearings.

Fig. 5Fuzzy acceleration factors for three accelerated operating conditions

Fig. 6Fuzzy coefficients in acceleration models

a)

b)

c)

6.4. Lifetime distribution selection and fuzzy reliability analysis for bearings

Once the fuzzy q-TTFs are given, the selection of lifetime distribution is essential for reliability analysis. In this case, the procedure of distribution selection proposed in Section 5.1 is used to compute the score of each candidate distribution for fuzzy q-TTFs ${\tilde{T}}_0\left[\alpha \right]$ at both sides. The membership value is equally divided into 10 intervals, which means that $M=$ 10 and the maximum score will be$\mathrm{}{h}_{\mathrm{m}\mathrm{a}\mathrm{x}}=2M+1=\text{21}$. Fig. 7 presents the results and indicates that normal, Log-normal and Weibull distributions are appropriate for describing bearing lifetime data, while exponential distribution is not suitable at all and extreme value distribution seems well with only one score. Although it has been demonstrated that the failure times for bearing follows Weibull distribution [31], the normal and Log-normal distributions also perform well in this study. The reason may be from the small sample (12 bearing lifetime data), which leads to the similar shape that those three distributions can describe. Thus, special attention should be given for other applications if there have no such prior knowledge. For bearings, the Weibull distribution is used for computing fuzzy reliability and MTTF.

Fig. 7The score of each distribution for bearing failure times

The fuzzy parameter for Weibull distribution is calculated at each alpha-cut, see Fig. 8. For all $\alpha \in \left[\text{0, 1}\right]$, $\tilde{\beta }\left[\alpha \right]$ is greater than 1 which means that the bearing experienced a wear-out phase. Interested that with different belief of the results, the fuzzy value $\tilde{\beta }\left[\alpha \right]$ varied from [2.80, 3.25] to the peaking 9.30. That is to say, with different belief, the recognition on the severity of bearing deterioration is significantly different. While in traditional analysis, only one crisp value can be obtained and such phenomenon cannot be identified.

Fig. 8Fuzzy parameters for the Weibull lifetime distribution

As described in the challenge details [15], the lifetimes for the majority of the bearings are from 1 to 7 (h). It may be interesting to analyze the reliability of this kind of bearing during the next hour in normal operating condition. Thus, the time interval [7, 8] (h) is substituted to Eq. (18) with fuzzy Weibull parameters. Fig. 9 shows the results, where the reliability results is in [0.8573, 0.9969] and have the peak value at 0.9999. The results demonstrate that the bearings have a higher possibility that is safety enough for operating under normal operating condition.

For users, the maintenance time is that of interest for making schedule to avoid the occurrence of failure and reduce their related costs. Hence, the fuzzy MTTF values is then computed using Eq. (19), see Fig. 10. The MTTF goes from 7.83 to 36.82 (h) and have the highest peak at 19.83 (h). In this case, it has the lowest possibility that maintenance schedule should be given before 7.83 (h), while the highest at 19.83 (h). Compared with traditional crisp value, the related risk can be directly expressed by fuzzy MTTF. For instance, if the traditional value is at the left side of fuzzy MTTF, the schedule will be conservative. Otherwise, it will be radical.

Fig. 9Fuzzy reliability for bearings at time interval [7, 8] (h)

Fig. 10Fuzzy MTTF for bearings