1 Introduction

Spindle is a key unit of NC machine tools. The reliability of spindle plays an important role in reliability of NC machine tools, so how to evaluate spindle reliability of accurately is great significance^[1]. Reliability evaluation of spindle is the basis of reliability design, test, operation and maintenance of NC machine tools.

It's well known that Weibull distribution is the most widely used in reliability studies such as software^[2-3], machine tool^[4-7], wind^[8-10] and human^[11-12]. The two-parameter Weibull distribution^{[6, 13-14]}, three-parameter Weibull distribution^[15] and mixed Weibull distribution^[16-17] are commonly used in reliability engineering.

Two-parameter and three-parameter Weibull distribution are only adapted to the single failure mode. In the case of failure modes are multiple for complex electromechanical system, mixed Weibull distribution is usually used. Different methods are used to estimate the parameters of mixed Weibull distribution. Maximum likelihood estimation (MLE) is widely used to estimate the parameters of mixed Weibull^[18-20]. Tan^[16] proposed a new approach to MLE of Weibull distribution with interval data. Mixed Weibull distributions was used for fitting failure times data by Razali^[17], and Maximum likelihood estimation was used to estimate the parameters.

Xia^[21] analyzed the reliability-based approach to bridge condition assessments with mixed Weibull distribution for finding the multi-modal peak stresses. Zhang^[22] established a mixed Weibull proportional hazard model for mechanical system failure prediction utilising lifetime and monitoring data. Castet^[23] provided an advanced parametric fit based on mixed Weibull distribution for nonparametric satellite reliability. From references, it will be found that most researchers used the MLE for parameters evaluation of mixed Weibull distribution. But MLE is more complex and may not exist, or it is difficult to get the final solution even if they exit.

The graphical method is another method for the estimation of mixed Weibull distribution^[24-26]. It is introduced by Jiang et al.^[24], and has been used extensively to estimate the parameters of mixed Weibull distribution. The method has been applied by Paulo^[25] to analyze the wind speed. Although graphical method can get the estimated value, the subjectivity is very strong and the result is not accuracy.

Because the defect of parameters estimation with graphical method, we have proposed a new method to estimate the parameters of mixed Weibull. The method is combined the graphical method and genetic algorithm (GA). It can identify the failure data obeying mixed Weibull distribution with WPP, and then obtain the observed value of the parameters (shape and scale) with graphical method. The following can get the search interval of shape and scale parameters with loose coefficient, finally optimized the parameters with GA. The method for the estimation of parameters is quick and accurate.

2 Failure Model Identification 2.1 Failure Data Collection

For evaluating the reliability level of spindle,42 spindles have been used in the reliability field test and 49 failures are collected. The time between failures is shown in Table 1.

2.2 Weibull Probability Paper

The time between failures of spindle may obey several types of distribution, such as normal distribution, two or three parameter Weibull distribution and mixed Weibull distribution.To selecting the best distribution, Weibull probability paper (WPP) is drawn to identify which distribution is fit^[24]. Weibull probability paper (WPP) is shown in Fig. 1.

If there is no obvious turning point on the WPP, presented in a straight line, two-parameter or three-parameter Weibull distribution is used to fit the data. If there is an obvious turning point on the WPP, the spindle will meet the three-parameter Weibull distribution or the mixed Weibull. As shown in Fig. 1, there is not a straight line and with an obvious turning point. So this paper proposes using mixed Weibull distribution for reliability modeling^[27].

3 Reliability Modeling 3.1 Mixed Weibull Model

The mixed Weibull model function is shown as follows^[28]:

R(t)=_pR₁(t)+_qR₂(t)=

$p\exp \left[{ - {{\left( {\frac{t}{{{a_1}}}} \right)}^{\beta 1}}} \right] + q\exp \left[{ - {{\left( {\frac{t}{{{a_2}}}} \right)}^{\beta 2}}} \right]$

(1)

where β₁ and β₂ are shape parameters; α₁ and α₂ are scale parameters; p and q are the proportion of the two-parameter Weibull distribution in mixed Weibull model.

For the spindle, there are two types of failure.One is sudden failure, it means failure by lost strength and the shape parameter β₁＜1, another is abrasion failure and the shape parameter β₂＞1^[27]. If α₁=α₂ and β₁=β₂, the mixed Weibull has transformed the two-parameter Weibull.

3.2 Parameter Estimation

Because the problems of graphical method and maximum likelihood estimation, this paper proposes a method that combined graphical method and genetic algorithm(GA) for parameter estimation.

GA is a kind of adaptive and intelligent search technology with global optimization. The most successful applications are complex nonlinear optimization problem and it will find a satisfactory solution. The search interval of shape and scale for GA is got with graphical method. The steps of parameter estimation are shown as follows:

(1) Calculate x_i and y_i(i=1,2,…,n).

Let all the time between failures in Table 1 in ascending order, and record as t_i,(i=1,2,…,n)^[27-28]. So x_i and y_i are got with Eq.(2).

$\left\{ \begin{align} & {{x}_{i}}=\ln \left( {{t}_{i}} \right) \\ & {{y}_{i}}=\ln \left[ -\ln \left( 1-\frac{i-0.3}{n+0.4} \right) \right] \\ \end{align} \right.$

(2)

(2) Plot asymptote in WPP.

Drawing (x_i,y_i),i=1,2,…,n, in the WPP, then searching two asymptote line y=ax+b as shown in Fig. 2. The asymptote record as L₁,L₂, and the expressions as follows:

$\left\{ \begin{align} & {{L}_{1}}:y=0.799lx-5.9315 \\ & {{L}_{2}}:y=3.290lx-25.5467 \\ \end{align} \right.$

(3)

(3) Define the search interval.

Through Eq.(3) to Eq.(9), it can obtain the observed value slope a^ and ordinate b^ via least square. Meanwhile, from Eq.(10) it can be got shape parameters β₁,β₂, and scale parameters α₁,α₂.

$\left\{ \begin{align} & \hat{b}={{l}_{xy}}/{{l}_{xx}} \\ & \hat{a}=\bar{y}-{{{\hat{b}}}_{xx}} \\ \end{align} \right.$

(4)

${{l}_{xx}}=\sum\limits_{i=1}^{n}{{{\left( {{x}_{i}}-\overset{-}{\mathop{\text{ }x}}\, \right)}^{2}}}=\sum\limits_{i=1}^{n}{x}\overset{2}{\mathop{\text{ }i}}\,-n{{\overset{-}{\mathop{\text{ }x}}\,}^{2}}$

(5)

${{l}_{xx}}=\sum\limits_{i=1}^{n}{{{\left( {{y}_{i}}-\overset{-}{\mathop{\text{ }y}}\, \right)}^{2}}}=\sum\limits_{i=1}^{n}{y}\overset{2}{\mathop{\text{ }i}}\,-n{{\overset{-}{\mathop{\text{ }y}}\,}^{2}}$

(6)

${l_{xx}} = \mathop \sum \limits_{i = 1}^n \left( {{x_i} - \mathop x\limits^ - } \right)\left( {{y_i} - \mathop y\limits^ - } \right) = \mathop \sum \limits_{i = 1}^n {x_i}{y_i} - n\mathop x\limits^ - \mathop y\limits^ - $

(7)

$\mathop x\limits^ - = \frac{1}{n}\mathop \sum \limits_{i = 1}^n {x_i}$

(8)

$\mathop y\limits^ - = \frac{1}{n}\mathop \sum \limits_{i = 1}^n {y_i}$

(9)

$\left\{ \begin{align} & \hat{\beta }=\hat{b} \\ & \hat{a}=\exp \left( -\hat{a}/\hat{b} \right) \\ \end{align} \right.$

(10)

According to the given loose coefficient k(0≤k≤1), and the expression of L₁ and L₂, we can obtained the four parameters observation through the least square α₁₀=1 668,β₁₀=0.779 2,α₂₀=2 355 and β₂₀=3.290 1. The search interval of each parameter is as follows:

$\left\{ \begin{align} & {{a}_{1}}\in \left[ \left( 1-k \right)\bullet {{a}_{10}},\left( 1+k \right)\bullet {{a}_{10}} \right] \\ & {{\beta }_{1}}\in \left[ \left( 1-k \right)\bullet {{\beta }_{10}},\left( 1+k \right)\bullet {{\beta }_{10}} \right] \\ & {{a}_{2}}\in \left[ \left( 1-k \right)\bullet {{a}_{20}},\left( 1+k \right)\bullet {{a}_{20}} \right] \\ & {{\beta }_{2}}\in \left[ \left( 1-k \right)\bullet {{\beta }_{20}},\left( 1+k \right)\bullet {{\beta }_{20}} \right] \\ \end{align} \right.$

(11)

For expanding the scope of parameters' search interval from Eq.(11), we have chosen loose coefficient k=0.8. Because shape parameter β₁＜1 and β₂＞1, so the max value of β₁ and the min value of β₂ is 1 . The interval of parameters is shown in Table 2.

(4) Parameters optimization.

Fitness function is the key of genetic algorithm. It has important influence on the speed of convergence and reconciliation precision. Because spindle failure data has the feature of small sample, so we extend d test as a fitness function.

Let all the failure data in ascending order, then computing the difference of hypothesis distribution (mixed Weibull distribution) function F₀(t_i) and the empirical distribution function F_n(t_i). The smaller of the difference, the better of it fits the hypothesis distribution. The largest value of the difference's absolute is regarded as test statistics D_n. Its expression is shown as follows:

${D_n} = \mathop {\sup }\limits_{ - \infty ＜ x ＜ + \infty } |{F_n}\left( {{t_i}} \right) - {F_0}\left( {{t_i}} \right)| = \max \left\{ {{d_i}} \right\} ＜ {D_n}{,_a}$

(12)

where F₀(t_i) is mixed Weibull distribution function; F_n(t_i) is empirical distribution function and the expression as shown in Eq.(13); d_i is distribution of deviation as shown in Eq.(14); D_n,α is critical value and the expression as shown in Eq.(15) when confidence level α=0.1.

${{F}_{n}}\left( {{t}_{i}} \right)=\left\{ \begin{align} & 0, t＜{{t}_{i}} \\ & \frac{i-0.3}{n+{{0.4}^{,}}}{{t}_{i}}＜t＜{{t}_{i}}_{+1} \\ & 1, t\ge {{t}_{n}} \\ \end{align} \right.$

(13)

${d_i}| = |{F_n}\left( {{t_i}} \right) - {F_0}\left( {{t_i}} \right)|$

(14)

${D_n}{,_a} = \frac{{1.22}}{{\sqrt n }}$

(15)

The test statistic D_n and the critical value D_n,α are analyzed. If the critical value is greater than the test statistics, the hypothesis is accepted, otherwise rejected. So the fitness is shown as follows:(x)=D_n-D_n,_a=

$\max \left\{ {|\mathop {{F_n}}\limits^{n, a} \left( {{x_i}} \right)|} \right\} - {D_{n, a}}$

(16)

Fitness function value is much smaller, which indicates that the hypothesis distribution is close to the empirical distribution, it means the hypothesis distribution is more accurate. At the same time, if Fitness(x)＜0, we can consider that the parameters have passed d test, and got the parameters estimation. We can get the parameters estimation and the fitness with Matlab. There are three kinds of operation in GA. Selection-reproduction, crossover as well as mutation. In the process of GA, the initial population is 200, and the generation gap is 0.9, and the crossover rate is 0.7, and the mutation rate is 0.01. After 500 iterations, the parameter estimates are obtained.

Because fitness is less than zero, so the estimated value of parameters is effective, and the time between failure obey mixed Weibull distribution. The estimated value is shown in Table 3. The reliability model graphic-GA is shown in Eq.(17). Eqs.(18)-(20) are the functions of graphical method, piecewise Weibull and two-Weibull.

$R\left( t \right) = 0.3254 \times \exp - \left[{{{\left( {\frac{t}{{361}}} \right)}^{0.8676}}} \right] + 0.6746 \times \exp \left[{{{\left( {\frac{t}{{2712}}} \right)}^{3.4895}}} \right]$

(17)

${R_1}\left( t \right) = 0.5814 \times \exp - \left[{{{\left( {\frac{t}{{1668}}} \right)}^{0.7992}}} \right] + 0.4186 \times \exp \left[{{{\left( {\frac{t}{{2355}}} \right)}^{3.2901}}} \right]$

(18)

${{R}_{2}}\left( t \right)=\left\{ \begin{align} & \exp \left[ {{\left( \frac{t}{5422} \right)}^{0.6681}} \right],0\le t\le 940 \\ & 0.786\times \exp \left[ {{\left( \frac{t}{2311} \right)}^{29547}} \right],940＜t \\ \end{align} \right.$

(19)

${R_3}\left( t \right) = \exp \left[{{{\left( {\frac{t}{{2126}}} \right)}^{0.8359}}} \right]$

(20)

3.3 Model Optimization 3.3.1 Average cumulative error test

The average cumulative error is the average of absolute value of difference between empirical distribution function and hypothesis distribution function. The smallest of the average cumulative error of the fitting function is the best fitting function. The average cumulative error is usually calculated with Eq.(21).

$\delta = \frac{1}{n}\mathop \sum \limits_{i = 1}^n |\mathop R\limits^ \wedge \left( {{t_i}} \right) - R\left( {{t_i}} \right)|$

(21)

Where R^(t_i) is reliability function of empirical distribution function; R(t_i) is reliability function of hypothesis distribution function.

3.3.2 Mean square error test

Mean square error (MSE) test refers to the mean value of the square root of the error of each measurement error. The smaller of MSE, the better of data fitting. MSE usually calculate with Eq.(22).

$\varepsilon = \frac{1}{n}\sqrt {\mathop \sum \limits_{i = 1}^n \left[{\mathop R\limits^ \wedge \left( {{t_i}} \right) - R\left( {{t_i}} \right)} \right]} $

(22)

The four reliability model has been optimized by d test, the average cumulative error test and MSE respectively. All the tests indicate that the mixed Weibull distribution is best. The test value is shown in Table 4.

Fig. 3 shows reliability function curve of spindle with different methods. From Fig. 3, the function with graphical-GA is the best to fit the scatter diagram. The spindle failure includes sudden failure and abrasion failure. The first sub section describes the sudden failure, which reflectes that the spindle is on early failure period. Because all the spindles are minted, it does not running-in well. The second sub section describes that the spindle is on abrasion failure period. After a period of use, the early malfunction has been eliminated.

4 Reliability Evaluation 4.1 MTBF Point Estimation

MTBF(Mean Time between Failure) point estimation is shown as follows:

$MTBF = p \bullet {a_1} \bullet \Gamma \left( {1 + \frac{1}{{{\beta _1}}}} \right) + q \bullet {a_2}\Gamma \left( {1 + \frac{1}{{{\beta _2}}}} \right)$

(23)

So MTBF of spindle can calculate with Eq.(21).

$MTBF = p \bullet {a_1} \bullet \Gamma \left( {1 + \frac{1}{{{\beta _1}}}} \right) + q \bullet {a_2}\Gamma \left( {1 + \frac{1}{{{\beta _2}}}} \right) = 1744\left( h \right)$

4.2 MTBF Interval Estimation

The hypothesis that θ^* is estimate value of the unknown parameter θ. P is the probability of θ less than positive number δ, which is shown in Eq.(24).

$P|\theta * - \theta | ＜ \delta = 1 - a$

(24)

where α is significance level; 1－α is probability of θ in (θ^*－δ,θ^*+δ), and (θ^－δ,θ^*+δ) is confidence interval.

The reliability of spindle is time terminated testing, and confidence interval 1－α is 90%. It will get quantiles $x\frac{{^2a}}{2}$ and $x{}_1^2{ - _{\frac{a}{2}}}$ as follows:

$p\left[{x\frac{{{{^2}_a}}}{2}\left( {2r} \right) \leqslant \frac{{2T*}}{\theta } \leqslant x\mathop 1\limits^2 - \frac{a}{2}\left( {2r + 2} \right)} \right] = 1 - a$

(25)

Hence Eq.(25) can be re-written as

$p\left[{\frac{{2T*}}{{x\mathop 2\limits^1 - \frac{a}{2}\left( {2r + 2} \right)}} \leqslant \theta \leqslant \frac{{2T*}}{{x\frac{{{{^2}_a}}}{2}\left( {2r} \right)}}} \right] = 1 - a$

(26)

The two-sided confidence interval of MTBF can be intended as follows:

$\frac{{2T*}}{{{x^2}\mathop {_{0.95}\left( {2r + 2} \right)}\limits^{} }} \leqslant \theta \leqslant \frac{{2T*}}{{{x^2}_{0.05}\left( {2r} \right)}}$

(27)

where r is failure number, and T^ is total testing time. Hence the two-sided confidence interval of spindle MTBF is (1 387.1,2 264.5), calculated with Eq.(27).

$\eqalign{ & {\theta _{\min }} = \frac{{2T*}}{{{x^2}\mathop {_{0.95}\left( {2r + 2} \right)}\limits^{} }} = \frac{{2 \times 86236}}{{{x^2}\mathop {_{0.95}\left( {2 \times 49 + 2} \right)}\limits^{} }} = \cr & {\theta _{\min }} = \frac{{2T*}}{{{x^2}\mathop {_{0.05}\left( {2r + 2} \right)}\limits^{} }} = \frac{{2 \times 86236}}{{{x^2}\mathop {_{0.05}\left( {2 \times 49} \right)}\limits^{} }} = 2264\left( h \right) \cr} $

5 Conclusions

In this paper, a new method has introduced to evaluating reliability of spindle. From the preceding analysis, the conclusions can be drawn as follows:

1) According to time terminated reliability testing of spindle in manufacture, failure data were collected, then making model identification through WPP. The WPP has the obvious turning point, so the data of spindle may be obey mixed Weibull distribution.

2) In view of graphical method with less precision and the complexity in calculating of MLE, it has combined graphical method with genetic algorithm for parameters estimation, and processed d test as a fitness function. With the comparative analysis, it can verify that graphical-GA methods has the estimation more objective and higher precision.

3) MTBF point estimation and interval estimation of spindle with mixed Weibull distributionis evaluated. The result can provide the basis for design and reliability growth of spindle.