- Email: [email protected]
Numerical method for Weibull generalized renewal process and its applications in reliability analysis of NC machine tools
Computers & Industrial Engineering 63 (2012) 1128–1134
Contents lists available at SciVerse ScienceDirect
Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie
Numerical method for Weibull generalized renewal process and its applications in reliability analysis of NC machine tools Zhi-Ming Wang a,⇑, Jian-Guo Yang b a b
School of Mechanical Engineering, Huaihai Institute of Technology, 59 Cang Wu Rd., Lianyungang 222005, PR China School of Mechanical Engineering, Shanghai Jiaotong University, 800 Dongchuan Rd., Shanghai 200240, PR China
a r t i c l e
i n f o
Article history: Received 25 November 2011 Received in revised form 18 June 2012 Accepted 20 June 2012 Available online 31 August 2012 Keywords: Weibull generalized renewal process Nonlinear programming Repairable system Restoration factor Reliability
a b s t r a c t In generalized renewal process (GRP) reliability analysis for repairable systems, Monte Carlo (MC) simulation method instead of numerical method is often used to estimate model parameters because of the complexity and the difficulty of developing a mathematically tractable probabilistic model. In this paper, based on the conditional Weibull distribution for repairable systems, using negative log-likelihood as an objective function and adding inequality constraints to model parameters, a nonlinear programming approach is proposed to estimate restoration factor for the Kijima type GRP model I, as well as the model II. This method minimizes the negative log-likelihood directly, and avoids solving the complex system of equations. Three real and different types of field failure data sets with time truncation for NC machine tools are analyzed by the proposed numerical method. The sampling formulas of failure times for the GRP models I and II are derived and the effectiveness of the proposed method is validated with MC simulation method. The results show that the GRP model is superior to the ordinary renewal process (ORP) and the power law non-homogeneous Poisson process (PL-NHPP) model. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Any system that can be restored to operating condition after experiencing a failure can be defined as a repairable system. The two widely used models for analyzing failure data of repairable systems are the ordinary renewal process (ORP) model, corresponding to perfect repair (see Blischke & Murphy, 1994), and the non-homogeneous Poisson process (NHPP) model, corresponding to minimal repair, respectively. The details can be found in Saldanha, Simone, and Frutoso (2001), Majeske (2007), and Weckman, Shell, and Marvel (2001). The ORP assumes that after a repair the system returns to an ‘‘as good as new’’ condition, and the NHPP assumes that the system returns to an ‘‘as bad as old’’ condition. However, the most repair activities may realistically not result in either of these two extreme situations but in a complicated intermediate one (general repair or imperfect repair). This imperfect repair restores the system operating state to somewhere between as good as new and as bad as old. It is possible, therefore, to model other repair assumptions such as the intermediate state of ‘‘better than old but worse than new’’.
⇑ Corresponding author. Present address: School of Mechanical Engineering, Huaihai Institute of Technology, 59 Cang Wu Rd., Lianyungang 222005, PR China. Tel.: +86 0518 85895000. E-mail addresses: [email protected] (Z.-M. Wang), Jgyang@sjtu. edu.cn (J.-G. Yang). 0360-8352/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cie.2012.06.019
Kijima and Sumita (1986) and Kijima (1989) proposed two virtual age models called general renewal process (GRP) model I and II. The Kijima type model I assumes that repairs only compensate for the damage created in the last period of operation; The Kijima type model II assumes that repairs could remove all of the damage accumulated up to the last failure time. Because of the complexities and the difficulties of developing a mathematically tractable probabilistic model, based on Kijima and Sumita’s work, Kaminskiy and Krivtsov (1998, 2000), and Krivtsov (2000) presented an approximate solution to the Kijima type model I by using Monte Carlo (MC) method for certain application areas. Yanez, Joglar, and Modarres (2002) studied the application of the maximum likelihood estimation (MLE) and Bayesian method in parameter estimation for GRP model I. Jacopino, Groen, and Mosleh (2006) used MCMC sample method to model imperfect inspection and maintenance in defense aviation only for GRP model I. Veber, Nagode, and Fajdiga (2008) used EM algorithm to estimate parameters of GRP model I based on the Weibull mixture distribution of the time to first failure. However, the distribution of the time to first failure is not always known because of lack of a large failure data set. On the other hand, as mentioned above, little research has been done on the GRP model II. In this paper, based on a nonlinear constrained programming numerical method, we not only explore using GRP I model but also GRP II model to model and analyze complex repairable systems with the underlying conditional Weibull distribution.
1129
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
2. Numerical methods of stochastic process in reliability analysis of repairable systems
ln L ¼
2.1. Ordinary renewal process
i ¼ 1; 2; . . . ; n
ð1Þ
The corresponding probability density function (pdf) is
f ðxi Þ ¼
kbxb1 i
expðkxbi Þ;
i ¼ 1; 2; . . . ; n
j¼1
i¼1
ð3Þ
When the observation stops right after the last event in time truncation or failure truncation, T j tnj ;j ¼ 0. Therefore, the log-likelihood function can be expressed as follows:
ðT j Sj Þ
Estimator of b can also be got by an iteration method, thus substituting it into the second equation of system of Eqs. (7), one can get another estimator for k. 2.3. Generalized renewal process Kijima and Sumita (1986) and Kijima (1989) developed an imperfect repair model by using the idea of the virtual age process for repairable systems. If a system has the virtual age Vi1 immediately after the (i 1)th repair, then the ith TBF Xi has the following conditional cumulative distribution function (the so-called underlying distribution):
Therefore, the conditional Weibull cdf is
Fðxi jv i1 Þ ¼ 1 expfk½ðxi þ v i1 Þb v bi1 g;
½ðT j t n ;j Þ þ j
j¼1
ð9Þ
and the corresponding pdf is
i ¼ 1; 2; . . . ; n
ð10Þ
There are two Kijima type GRP models. 2.3.1. Generalized renewal process I In GRP I, it is assumed that the ith repair would only compensate for the damage created during the time between the (i 1)th and the ith failure. qi is the degree of the ith repair quality. For simplicity, we define q1 = q2 = = qn = q. With this assumption the virtual age of the system after the ith repair is: n X xi ¼ qti
v i ¼ v i1 þ qxi ¼ q
ð11Þ
i¼1
ð12Þ
and the corresponding likelihood function of the TBF for k repairable systems with time truncation can be expressed as follows:
The MLE for b and k are:
j¼1
i ¼ 1; 2; . . . ; n
i ¼ 1; 2; . . . ; n
ð4Þ
j¼1
8 Pk Pk Pnj b Pk Pnj ½ðT j t n ;j Þb lnðT j tn ;j Þþ x ln xi:j ln xi;j > i¼1 i;j > j¼1 j¼1 j j j¼1 1 > ¼ P P P Pk i¼1 nj >b k k < ðT t n ;j Þb þ xb n i¼1 i;j j¼1 j j¼1 j¼1 j j Pk > n > j¼1 j > > k ¼ Pk Pk Pnj b : b
ð8Þ
n h io f ðxi ; q; k; bjv i1 Þ ¼ kbðxi þ qti1 Þb1 exp k ðxi þ qti1 Þb ðqt i1 Þb ;
j¼1 i¼1
" # nj k X k X X k xbi;j þ ðT j tnj ;j Þb
Fðxi þ v i1 Þ Fðv i1 Þ 1 Fðv i1 Þ
The conditional pdf of GRP I model is
nj k k X X X ln L ¼ nj ðln k þ ln bÞ þ ðb 1Þ ln xi;j
j¼1 i¼1
j¼1
ð7Þ
f ðxi jv i1 Þ ¼ kbðxi þ v i1 Þb1 expfk½ðxi þ v i1 Þb v bi1 g;
where k and b are the scale and shape parameters, respectively. Therefore, the likelihood function of the TBF for k repairable systems with time truncation is ( ) nj k h iY Y b Lðxi;j ; T j tnj ;j ; k; bÞ ¼ knj bnj exp kðT j tnj ;j Þb ½xb1 expðkx Þ i;j i;j
j¼1
8 Pk b Pk Pnj ðT ln T j Sbj ln Sj Þ ln t i;j > j¼1 j j¼1 1 > ¼ P Pk i¼1 > k nj > j¼1 > : k ¼ Pk b b
ð2Þ
ði ¼ 1; 2; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ
ð6Þ
j¼1
The MLE for b and k are:
Fðxi Þ ¼ PrfX i 6 xi jV i1 ¼ v i1 g ¼
For ORP, in the case of the underlying distribution is the twoparameter Weibull distribution, the cumulative distribution function (cdf) of x for a single system is
Fðxi Þ ¼ 1
j¼1 i¼1
j¼1
We consider a component/system which under goes instantaneous general repairs after a failure. The repair time is assumed to be negligible so that the failure can be viewed as point processes. Weibull distribution is one of the most commonly used distribution in reliability analysis for repairable systems because of its flexibility and applicability to various failure processes. Therefore, we select Weibull distribution as the underlying distribution in this study. Let t be the successive failure times, and x be the time between failures (TBFs). Now, consider k repairable systems, observed from the starting time t0 to the end time Tj (j = 1, 2, . . . , k). The successive failure times of the jth system are denoted by t0 ; t1;j ; t 2;j ; . . . ; t nj ;j ; T j , and the TBF of the jth system are denoted by x0 ; x1;j ; x2;j ; . . . ; xnj ;j ; T j tnj ;j : nj is the number of failures and Tj is the truncated time for the jth system. ti,j is the failure time of the jth system at the ith failure, xi,j is the TBF of the jth system between the (i 1)th and the ith failure. For convenience, we define t0 = 0, x0 = 0.
expðkxbi Þ;
nj k k X k X X X nj ðln k þ ln bÞ þ ðb 1Þ ln t i;j k ðT bj Sbj Þ
Lðxi;j ; T j tnj ;j ; k; b; qÞ ¼
k Y
n h io knj bnj exp k ðT j tnj ;j þ qt nj ;j Þb ðqt nj ;j Þb
j¼1
ð5Þ
x
i¼1 i;j
Estimator of b can be found from the first equation of the system of Eqs. (5) by an iteration method, then substituting it into the second equation, one can get another estimator of k. 2.2. Power law non-homogeneous Poisson process
nj n h io Y ðxi;j þ qt i1;j Þb1 exp k ðxi;j þ qt i1;j Þb ðqt i1;j Þb
! ð13Þ
i¼1
Thus, the log-likelihood function is
ln L ¼
k k h i X X nj ðln k þ ln bÞ k ðT j t nj ;j þ qtnj ;j Þb ðqtnj ;j Þb j¼1
j¼1 nj k X X
k
½ðxi;j þ qti1;j Þb ðqt i1;j Þb
j¼1 i¼1
Similarly, consider k systems. S and T denote the starting time and the end time of observation accounting for censoring. The log-likelihood function of PL-NHPP model is
þ ðb 1Þ
nj k X X j¼1 i¼1
lnðxi;j þ qt i1;j Þ
ð14Þ
1130
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
Taking the derivatives of Eq. (14) with respect to k, b and q, and setting them equal to zero, we can get k k X 1X nj ½ðT j tnj ;j þ qt nj ;j Þb ðqt nj ;j Þb k j¼1 j¼1
nj k X X ½ðqti1;j þ xi;j Þb ðqti1;j Þb ¼ 0
ð15Þ
j¼1 i¼1
k X k ½ðT t nj ;j þ qtnj ;j Þb lnðT tnj ;j þ qtnj ;j Þ ðqtnj ;j Þb lnðqt nj ;j Þ j¼1 nj k X X ½ðqti1;j þ xi;j Þb lnðqti1;j þ xi;j Þ ðqti1;j Þb lnðqti1;j Þ þk j¼1 i¼1 n
j k k X X 1X nj lnðqt i1;j þ xi;j Þ ¼ 0 b j¼1 j¼1 i¼1
ð16Þ
greater than zero. According to the GRP model, there are five cases for q. For q = 0, a system is restored to an ORP (perfect maintenance) condition, while for q = 1, the system corresponds to a NHPP (minimal maintenance) state. The case of 0 < q < 1 corresponds to the intermediate ‘‘better than old but worse than new’’ repair assumption. For the cases where q > 1, the system is in a condition of worse than old. Similarly, the case where q < 0 suggests the system is restored to a state of better than new. Therefore, no constraints are added to variable q. So, it is possible to define the problem of solving the complex system of three equations as a nonlinear constrained programming problem. In this research, the objective function is negative loglikelihood function that is related to model parameters, and the constrained variables are model parameters. They can be expressed by
minð ln LÞ
ð20Þ
s:t: k > 0; 0 < b < 10 ðb 1Þ
nj k X X j¼1 i¼1
nj k X X t i1;j t i1;j ½ðqti1;j þ xi;j Þb1 kb ðqti1;j þ xi;j Þ j¼1 i¼1
ðqt i1;j Þb1 kb
k X
t nj ;j ½ðT j t nj ;j þ qtnj ;j Þb1 ðqtnj ;j Þb1 ¼ 0
j¼1
2.3.2. Generalized renewal process II The GRP II model assumes that repairs could remove all of the damage accumulated up to the last failure time. With this assumption the virtual age of a system after the ith repair is:
ð17Þ We can see that parameter estimation for GRP I model based on the MLE method leads to a complex system of three equations that should be solved simultaneously. Unfortunately, the system of equations mentioned above does not have a closed-form mathematical expression. It is difficult to solve it directly. Therefore, we use a nonlinear constrained programming approach to estimate parameters of GRP model. Mathematical model of finding minimum value of nonlinear constrained programming is
8 min f ðxÞ > > < hi ðxÞ ¼ 0; ði ¼ 1; 2; . . . ; uÞ > g ðxÞ P 0; ðj ¼ 1; 2; . . . ; v Þ > : j x2D
ðkÞ
pðx ; M k Þ ¼ f ðx Þ þ M k
( u X i¼1
i X qimþ1 xm
Therefore, the conditional pdf of GRP model II is
f ðxi ; q; k; bjv i1 Þ ¼ kbðxi þ
i1 X
qim xm Þb1 expfk½ðxi
m¼1
þ
i1 X
qim xm Þb ð
m¼1
ð18Þ
ðkÞ
hi ðx Þ
2
) v X 2 ðkÞ ; þ minð0; g j ðx ÞÞ
i1 X qim xm Þb g; i
m¼1
ð22Þ
The corresponding likelihood function of the TBF for k repairable systems with time truncation is
Lðxi;j ; T j t nj ;j ; k; b; qÞ 8 2 0 !b nj k < Y X n n nj mþ1 @k j b j exp k4 T j t n þ ¼ q x m;j j;j : m¼1 j¼1 0 !b 3 9 !b1 nj nj i1 = Y X X nj mþ1 im 5 @ q xm;j xi;j þ q xm;j ; i¼1 m¼1 m¼1 8 2 !b !b 3911 = < i1 i1 X X exp k4 xi;j þ qim xm;j qim xm;j 5 AA ; : m¼1 m¼1 ð23Þ
j¼1
ðk ¼ 1; 2; . . .Þ
ð21Þ
m¼1
¼ 1; 2; . . . ; n
where f(x) is objective function, x is variable, hi(x) is equality constrained function, gj(x) is inequality constrained function, and D is feasible region. Exterior penalty method can convert constraint problem into a series unconstraint problem by using reconstruction penalty function, and then solve it with unconstraint optimization method. So, this method is also called sequential unconstrained minimizing method. The penalty function is ðkÞ
v i ¼ qðv i1 þ xi Þ ¼ qi x1 þ qi1 x2 þ þ qxi ¼
ð19Þ
where Mk is the penalty factor of the kth calculating. The detailed procedures are given as follows:
The log-likelihood function is
ln L ¼
k X nj ðln k þ ln bÞ j¼1
(1) Set k = 1, give an initial point x(1), the first penalty factor M1 > 0, amplifying coefficient C > 1, allowable error e > 0. (2) Solve unconstraint minimal value of the penalty function p(x(k), Mk). (3) Estimate the calculating accuracy. At point x(k), if the penalty P term Mk vj¼1 ½minð0; g j ðxðkÞ ÞÞ2 < e, then stop calculation, and get an approximate minimal solution x(k) to problem; Otherwise, set Mk+1 = CMk, k = k + 1, go back step (2). To constrain variables, set b be between the values of 0 and 10, since almost all expected solutions are in this range. k can be set
2 !b !b 3 nj nj k X X X n mþ1 n mþ1 j j q xm;j q xm;j 5 k 4 T j t nj ;j þ m¼1
j¼1
m¼1
2 !b !b 3 nj k X i1 i1 X X X im im 4 xi;j þ q xm;j q xm;j 5 k m¼1
j¼1 i¼1
þ ðb 1Þ
nj k X X j¼1 i¼1
ln xi;j þ
m¼1 i1 X m¼1
! im
q
xm;j
ð24Þ
1131
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
We also use nonlinear constrained programming approach to estimate parameters of the GRP model II. 2.3.3. Validation of the generalized renewal process model I and II We use Monte Carlo simulation method to validate the GRP models I and II. For multiple systems, from Eq. (9), we can get the sampling formula of the Weibull distribution for the TBF
n o exp k½ðxi;j þ v i1;j Þb v bi1;j ¼ U i;j
ð25Þ
where Ui,j is random number of uniform distribution on the interval [0, 1]. For GRP I model, the TBF can be got by
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b v bi1;j ln Ui;j k rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b ¼ qti1;j þ ðqti1;j Þb ln U i;j k
xi;j ¼ v i1;j þ
4. Results and discussion
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 1 t i;j ¼ t i1;j þ xi;j ¼ ð1 qÞti1;j þ ðqt i1;j Þb ln U i;j ði k ¼ 1; 2; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ
Table 1 Failure data of single machine tool with time truncated for case 1.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b ¼ ln U i;j ðt0 ¼ 0Þ k
ð28Þ
For GRP II model, the TBF can be got by i1 X
qim xm;j þ
m¼1
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xi1 b 1 b im x q ln U i;j m;j m¼1 k
ð29Þ
where i1 X
qim xm;j ¼
m¼1
¼
i1 X
qim t m;j
m¼1
i1 X
i2 X
q
t m;j
m¼1
¼
i2 X
No. of failures
Failure times
No. of failures
Failure times
1 2 3 4 5 6 7 8
50.29 531.42 914.70 963.27 1029.19 1488.78 1953.16 2380.42
9 10 11 12 13 14 15
2527.76 2839.82 3371.14 3562.01 3649.22 3801.43 4140.00
i1 X qim tm1;j
m¼1 im
Using our proposed method, the results of ORP, PL-NHPP, GRP I and GRP II for three real failure data sets of NC machine tools are all given in Table 4.
ð27Þ
where
xi;j ¼
We illustrate the applications of the proposed approach with three real cases for reliability analysis of NC machine tools. Failure times of NC machine tools are always time truncated data coming from field because of the limitation of time and money. Therefore, three real field failure data sets of NC machine tools are all time truncated data (in cumulative operating hours). Case 1: These data, given in Table 1, contain 15 failure times with time truncation from a NC machine tool and the last one is right truncated time. Case 2: 29 time truncation field failure data of another NC machine tool were collected, as shown in Table 2. Case 3: Table 3 gives 28 failure time data that are from four identical machine tools, and the last one of every machine tool is right truncated time taken from Zhang, Jia, and Shen (2005). Unlike cases 1 and 2, these failure data are multiple system failure data.
ð26Þ
So the sampling formula of failure times can be expressed as follows:
t 1;j
3. Numerical examples
im1
q
Table 2 Failure data of single machine tool with time truncated for case 2.
t m;j
m¼0
qim t m;j þ qt i1;j
m¼1
i2 X
qim1 t m;j qi1 t 0;j
m¼1
¼ qti1;j þ ðq 1Þ
i2 X
qim1 t m;j
ð30Þ
m¼1
So the sampling formula of failure times can be expressed as follows:
No.
Failure times
No.
Failure times
No.
Failure times
1 2 3 4 5 6 7 8 9 10
27.51 367.52 394.52 395.64 406.75 432.49 514.17 855.57 864.85 953.02
11 12 13 14 15 16 17 18 19 20
1039.36 1357.80 1680.92 1850.55 1926.98 2398.21 2430.61 2517.04 2600.22 2796.49
21 22 23 24 25 26 27 28 29
2867.40 3299.82 3387.57 3468.58 3688.63 3780.33 3862.50 3955.48 4152.00
t i;j ¼ t i1;j þ xi;j ¼ ð1 qÞðt i1;j þ
i2 X
qim1 tm;j Þ
Table 3 Failure data of multiple machine tools with time truncated for case 3.
m¼1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u i2 u X b 1 þ t½qti1;j þ ðq 1Þ qim1 t m;j b ln U i;j k m¼1 ði ¼ 3; 4; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ where we define t 0 ¼ 0; Thus
Pi2
im1 t m;j m¼1 q
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 b 1 > < t 1;j ¼ k ln U i;j > :
t 2;j ¼ ð1 qÞt 1;j þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ðqt1;j Þb 1k ln U i;j
ð31Þ ¼ 0ði ¼ 1; 2Þ.
ð32Þ
1 2 3 4 5 6 7 8 9 10 11
No. 1
No. 2
No. 3
No. 4
16 56 116 196 356 676 1028 1428 1828 2488 2522
60 154 304 469 539
96 192 400 800 872
10 74 274 506 1226 2394 2444
1132
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
Table 4 Estimated parameters and –ln L values with the different models for three cases. Model
ORP GRP I GRP II PL-NHPP
No. 1
No. 2
No. 3
b
k
q
ln L
b
k
q
ln L
b
k
q
ln L
1.523 1.523 1.523 0.960
0.0002 0.0002 0.0002 0.0047
0 0 0 1
92.171 92.171 92.171 93.640
1.013 0.766 0.870 0.870
0.0063 0.0327 0.0199 0.0199
0 0.109 1 1
167.972 167.494 167.691 167.691
1.049 0.569 0.569 0.569
0.0028 0.0947 0.0947 0.0947
0 1 1 1
157.927 152.579 152.579 152.579
From Table 4, we can find that in case 1, a q value of 0 suggests that the machine tool’s mending leaves it in the condition of as good as new. In this case we have an ordinary renewal process; In case 2, based on the value of negative log-likelihood function, the GRP I model is superior to the GRP II model, because model I has a smaller –ln L value than that of model II. The best GRP model for the failure data, as determined by the –ln L, is the model with the lowest –ln L value. A q value of 0.109 suggests that the machine tools’ maintenance leave them in the intermediate condition of better than old but worse than new. In this case we have a GRP I; In case 3, a q value of 1 suggests that four machine tools’ repair leaves them in the condition of as bad as old, and we can use a PL-NHPP model to describe the failure process of these four NC machine tools. On the other hand, it is noticed that the results of the GRP I and II model are the same in cases 1 and 3, respectively. Through observations and simulations, Jacopino, Groen, and Mosleh (2004) and Jacopino (2005) also pointed out that at a low number of renewals there is little difference between the GRP I and II model. However, as the number of renewals increases, the difference between the two models becomes significant because of the variation in the underlying virtual age equations of both models. In case 1, we also notice that the results of the ORP model and the GRP model are the same for this particular data set; In case 3, the results of the PL-NHPP model and the GRP model are also the same. In fact, from Eqs. (4), (6), (14), and (24), we can see that when q = 0, the models for ORP, GRP I and GRP II all have the same form of log-likelihood function corresponding to perfect repair; when q = 1, PL-NHPP, GRP I and GRP II also have the same form of log-likelihood function corresponding to minimal repair. This means that the ORP and the PLNHPP are special cases of the GRP.
There are several factors which affect the speed of convergence, these factors include penalty factor Mk, amplifying coefficient C, allowable error e, and initial point x(1) as well. The proposed method is especially sensitive to the initial point x(1), we can use the estimators of the ORP model or the PL-NHPP model as starting point to speed up the convergence of the solution and get a global optimization values rapidly. The computation time using our method in a 1000 MHz Personal Computer is less than 30 s. Compared with one minute computation time of MC method used by Yanez et al. (2002), our method is more efficient. In order to validate the effectiveness of proposed method for GRP I and II model, Monte Carlo simulation method is used. Note that a large sample size generated by a computer according to the sample formulas of failure times is needed when one uses MC method to obtain close and stable estimations. For GRP model I, the details are shown as follows: (1) In the sampling formula (27), set k, b and q be the results of case 2 for GRP I model; (2) generate different sample size data with truncated time is 4152, respectively; (3) estimate model parameters and calculate –ln L values. The results are shown in Fig. 1, we can see that simulation results including –ln L values are close to real values when sample size is up to 2000. Base on the –ln L values, therefore, we say that our method is effective. The simulation value of q is slightly larger than real value, this is because that in this model, for convenience, we assume restoration factor q are all equal to each other. However, in reality, the degree of repair after a failure is different due
Fig. 1. Simulation results of GRP I for case 2.
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
to the different failure modes and positions, etc. Uematsu and Nishida (1987) and Kijima (1989) extended restoration factor q to the case that it is a random variable taking a value between 0 and 1. In this study, we can regard q as an average value of repair quality. In GRP II, Mettas and Zhao (2005) analyzed 56 failure truncated data and given the estimators of model parameters, we select them as the known values in the sampling formula (31) to validate our proposed method for GRP model II. The results are shown in Fig. 2. It is found that the simulation results are close to real values. Based on a MC simulation, Fig. 3 shows the cumulative number of failures against the time in case 2 for GRP I. It is noticed that the number of failures are different in three models. When q = 0.109, the GRP I model fits the actual data better than the other two models corresponding to q = 0 and q = 1, respectively.
1133
As mentioned in Section 1, there are several methods for GRP I in literature. Among them, the MLE method discussed by Yanez et al. (2002) is based on MC simulation with numerical method, so the computation of this method is quite time-consuming and slow. Other two methods including MCMC sampling method used by Jacopino et al. (2006) and EM algorithm used by Veber et al. (2008) are both have very limited application areas, for the time to first failure must be known and estimated from a large set of data. However, when one uses our proposed method to analyze GRP, the time to first failure distribution and a large set of data are all unnecessary. In addition, three methods mentioned above are all only for the GRP I. In this study, we not only explore using GRP I model but also GRP II model to model and analyze complex repairable systems. Therefore, compared with the other approaches in literature, our GRP numerical approach is simple and efficient.
Fig. 2. Simulation results of GRP II.
Fig. 3. Cumulative number of failures versus time plot.
1134
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
5. Conclusions In this paper, a nonlinear constrained programming method is proposed based on the conditional Weibull distribution to estimate the parameters of the GRP model with different types of data. To avoid solving complex system of equations, the proposed method uses negative log-likelihood function as objective optimal function and adds inequality constraints to model parameters. The sampling formulas of failure times for GRP models I and II are given and the effectiveness of the proposed method is validated with simulation method. The GRP model considers restoration factor in reliability analysis of repairable systems, therefore, its evaluation results are close to real maintenance circumstance than the ORP model and the PLNHPP model, the latter are only two special cases where restoration factor is 0 and 1 respectively. In case where the amount of failure data is extremely limited, other methods such as a Bayesian approach based on prior knowledge may be used to remedy the uncertainty of estimations. Acknowledgements We would like to thank the anonymous reviewer and the editor for their valuable comments and suggestions, which have greatly enhanced the clarity of the paper. The authors also wish to acknowledge the partial financial support of this study by ‘‘Top Grade NC Machine Tools and Fundamental Manufacturing Equipment’’ National Science and Technology Key Special Project of China (No. 2009ZX04014-022)’’. References Blischke, W., & Murphy, D. (1994). Warranty cost analysis. New York: Marcel Dekker. Jacopino, A. (2005). Generalization and Bayesian solution of the general renewal process for modeling the reliability effects of imperfect inspection and maintenance based on imprecise data. Dissertation for the degree of Doctor of Philosophy, University of Maryland.
Jacopino, A., Groen, F., & Mosleh, A. (2004). Behavioural study of the general renewal process. In Proceedings of the annual reliability and maintainability symposium, 26–29 January 2004 (pp. 237–242). Los Angeles, CA, USA. Jacopino, A., Groen, F., & Mosleh, A. (2006). Modeling imperfect inspection and maintenance in defense aviation through Bayesian analysis of the KIJIMA type I general renewal process (GRP). In Proceedings of annual reliability and maintainability symposium, 23–26 January 2006 (pp. 470–475). Newport Beach, California, USA. Kaminskiy, M. & Krivtsov, V. (2000). G-renewal process as a model for statistical warranty claim prediction. In Proceedings of the annual reliability and maintainability symposium, 24–27 January 2000 (pp. 276–280). Los Angeles, CA, USA. Kaminskiy, M., & Krivtsov, V. (1998). A Monte Carlo approach to repairable system reliability analysis. Probabilistic safety assessment and management. New York: Springer, pp. 1063–1068. Kijima, M. (1989). Some results for repairable systems with general repair. Journal of Applied Probability, 26, 89–102. Kijima, M., & Sumita, N. (1986). A useful generalization of renewal theory: counting process governed by non-negative Markovian increments. Journal of Applied Probability, 23, 71–88. Krivtsov, V. (2000). A Monte Carlo approach to modeling and estimation of the generalized renewal process in repairable system reliability analysis. Dissertation for the Degree of Doctor of Philosophy, University of Maryland. Majeske, K. D. (2007). A non-homogeneous Poisson process predictive model for automobile warranty claims. Reliability Engineering and System Safety, 92, 243–251. Mettas, A. & Zhao, W. (2005). Modeling and analysis of repairable systems with general repair. In Proceedings of the annual reliability and maintainability symposium, 24–27 January 2005 (pp. 176–182). Alexandria, Virginia, USA. Saldanha, P. L. C., Simone, E. A., & Frutoso, P. F. (2001). An application of nonhomogeneous Poisson point processes to the reliability analysis of service water pumps. Nuclear Engineering and Design, 210, 125–133. Uematsu, K., & Nishida, T. (1987). Branching non-homogeneous Poisson process and its application to a replacement model. Microelectronics and Reliability, 27, 685–691. Veber, B., Nagode, M., & Fajdiga, M. (2008). Generalized renewal process for repairable systems based on finite Weibull mixture. Reliability Engineering and System Safety, 93, 1147–1461. Weckman, G. R., Shell, R. L., & Marvel, J. H. (2001). Modeling the reliability of repairable systems in the aviation industry. Computers and Industrial Engineering, 40, 51–63. Yanez, M., Joglar, F., & Modarres, M. (2002). Generalized renewal process for analysis of repairable systems with limited failure experience. Reliability Engineering and System Safety, 77, 167–180. Zhang, Y., Jia, Y., & Shen, G. (2005). Research on model of failure distribution for numerical machine with random ending method. Systems Engineering Theory and Practice, 2, 134–138. in Chinese.
Contents lists available at SciVerse ScienceDirect
Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie
Numerical method for Weibull generalized renewal process and its applications in reliability analysis of NC machine tools Zhi-Ming Wang a,⇑, Jian-Guo Yang b a b
School of Mechanical Engineering, Huaihai Institute of Technology, 59 Cang Wu Rd., Lianyungang 222005, PR China School of Mechanical Engineering, Shanghai Jiaotong University, 800 Dongchuan Rd., Shanghai 200240, PR China
a r t i c l e
i n f o
Article history: Received 25 November 2011 Received in revised form 18 June 2012 Accepted 20 June 2012 Available online 31 August 2012 Keywords: Weibull generalized renewal process Nonlinear programming Repairable system Restoration factor Reliability
a b s t r a c t In generalized renewal process (GRP) reliability analysis for repairable systems, Monte Carlo (MC) simulation method instead of numerical method is often used to estimate model parameters because of the complexity and the difficulty of developing a mathematically tractable probabilistic model. In this paper, based on the conditional Weibull distribution for repairable systems, using negative log-likelihood as an objective function and adding inequality constraints to model parameters, a nonlinear programming approach is proposed to estimate restoration factor for the Kijima type GRP model I, as well as the model II. This method minimizes the negative log-likelihood directly, and avoids solving the complex system of equations. Three real and different types of field failure data sets with time truncation for NC machine tools are analyzed by the proposed numerical method. The sampling formulas of failure times for the GRP models I and II are derived and the effectiveness of the proposed method is validated with MC simulation method. The results show that the GRP model is superior to the ordinary renewal process (ORP) and the power law non-homogeneous Poisson process (PL-NHPP) model. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Any system that can be restored to operating condition after experiencing a failure can be defined as a repairable system. The two widely used models for analyzing failure data of repairable systems are the ordinary renewal process (ORP) model, corresponding to perfect repair (see Blischke & Murphy, 1994), and the non-homogeneous Poisson process (NHPP) model, corresponding to minimal repair, respectively. The details can be found in Saldanha, Simone, and Frutoso (2001), Majeske (2007), and Weckman, Shell, and Marvel (2001). The ORP assumes that after a repair the system returns to an ‘‘as good as new’’ condition, and the NHPP assumes that the system returns to an ‘‘as bad as old’’ condition. However, the most repair activities may realistically not result in either of these two extreme situations but in a complicated intermediate one (general repair or imperfect repair). This imperfect repair restores the system operating state to somewhere between as good as new and as bad as old. It is possible, therefore, to model other repair assumptions such as the intermediate state of ‘‘better than old but worse than new’’.
⇑ Corresponding author. Present address: School of Mechanical Engineering, Huaihai Institute of Technology, 59 Cang Wu Rd., Lianyungang 222005, PR China. Tel.: +86 0518 85895000. E-mail addresses: [email protected] (Z.-M. Wang), Jgyang@sjtu. edu.cn (J.-G. Yang). 0360-8352/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cie.2012.06.019
Kijima and Sumita (1986) and Kijima (1989) proposed two virtual age models called general renewal process (GRP) model I and II. The Kijima type model I assumes that repairs only compensate for the damage created in the last period of operation; The Kijima type model II assumes that repairs could remove all of the damage accumulated up to the last failure time. Because of the complexities and the difficulties of developing a mathematically tractable probabilistic model, based on Kijima and Sumita’s work, Kaminskiy and Krivtsov (1998, 2000), and Krivtsov (2000) presented an approximate solution to the Kijima type model I by using Monte Carlo (MC) method for certain application areas. Yanez, Joglar, and Modarres (2002) studied the application of the maximum likelihood estimation (MLE) and Bayesian method in parameter estimation for GRP model I. Jacopino, Groen, and Mosleh (2006) used MCMC sample method to model imperfect inspection and maintenance in defense aviation only for GRP model I. Veber, Nagode, and Fajdiga (2008) used EM algorithm to estimate parameters of GRP model I based on the Weibull mixture distribution of the time to first failure. However, the distribution of the time to first failure is not always known because of lack of a large failure data set. On the other hand, as mentioned above, little research has been done on the GRP model II. In this paper, based on a nonlinear constrained programming numerical method, we not only explore using GRP I model but also GRP II model to model and analyze complex repairable systems with the underlying conditional Weibull distribution.
1129
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
2. Numerical methods of stochastic process in reliability analysis of repairable systems
ln L ¼
2.1. Ordinary renewal process
i ¼ 1; 2; . . . ; n
ð1Þ
The corresponding probability density function (pdf) is
f ðxi Þ ¼
kbxb1 i
expðkxbi Þ;
i ¼ 1; 2; . . . ; n
j¼1
i¼1
ð3Þ
When the observation stops right after the last event in time truncation or failure truncation, T j tnj ;j ¼ 0. Therefore, the log-likelihood function can be expressed as follows:
ðT j Sj Þ
Estimator of b can also be got by an iteration method, thus substituting it into the second equation of system of Eqs. (7), one can get another estimator for k. 2.3. Generalized renewal process Kijima and Sumita (1986) and Kijima (1989) developed an imperfect repair model by using the idea of the virtual age process for repairable systems. If a system has the virtual age Vi1 immediately after the (i 1)th repair, then the ith TBF Xi has the following conditional cumulative distribution function (the so-called underlying distribution):
Therefore, the conditional Weibull cdf is
Fðxi jv i1 Þ ¼ 1 expfk½ðxi þ v i1 Þb v bi1 g;
½ðT j t n ;j Þ þ j
j¼1
ð9Þ
and the corresponding pdf is
i ¼ 1; 2; . . . ; n
ð10Þ
There are two Kijima type GRP models. 2.3.1. Generalized renewal process I In GRP I, it is assumed that the ith repair would only compensate for the damage created during the time between the (i 1)th and the ith failure. qi is the degree of the ith repair quality. For simplicity, we define q1 = q2 = = qn = q. With this assumption the virtual age of the system after the ith repair is: n X xi ¼ qti
v i ¼ v i1 þ qxi ¼ q
ð11Þ
i¼1
ð12Þ
and the corresponding likelihood function of the TBF for k repairable systems with time truncation can be expressed as follows:
The MLE for b and k are:
j¼1
i ¼ 1; 2; . . . ; n
i ¼ 1; 2; . . . ; n
ð4Þ
j¼1
8 Pk Pk Pnj b Pk Pnj ½ðT j t n ;j Þb lnðT j tn ;j Þþ x ln xi:j ln xi;j > i¼1 i;j > j¼1 j¼1 j j j¼1 1 > ¼ P P P Pk i¼1 nj >b k k < ðT t n ;j Þb þ xb n i¼1 i;j j¼1 j j¼1 j¼1 j j Pk > n > j¼1 j > > k ¼ Pk Pk Pnj b : b
ð8Þ
n h io f ðxi ; q; k; bjv i1 Þ ¼ kbðxi þ qti1 Þb1 exp k ðxi þ qti1 Þb ðqt i1 Þb ;
j¼1 i¼1
" # nj k X k X X k xbi;j þ ðT j tnj ;j Þb
Fðxi þ v i1 Þ Fðv i1 Þ 1 Fðv i1 Þ
The conditional pdf of GRP I model is
nj k k X X X ln L ¼ nj ðln k þ ln bÞ þ ðb 1Þ ln xi;j
j¼1 i¼1
j¼1
ð7Þ
f ðxi jv i1 Þ ¼ kbðxi þ v i1 Þb1 expfk½ðxi þ v i1 Þb v bi1 g;
where k and b are the scale and shape parameters, respectively. Therefore, the likelihood function of the TBF for k repairable systems with time truncation is ( ) nj k h iY Y b Lðxi;j ; T j tnj ;j ; k; bÞ ¼ knj bnj exp kðT j tnj ;j Þb ½xb1 expðkx Þ i;j i;j
j¼1
8 Pk b Pk Pnj ðT ln T j Sbj ln Sj Þ ln t i;j > j¼1 j j¼1 1 > ¼ P Pk i¼1 > k nj > j¼1 > : k ¼ Pk b b
ð2Þ
ði ¼ 1; 2; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ
ð6Þ
j¼1
The MLE for b and k are:
Fðxi Þ ¼ PrfX i 6 xi jV i1 ¼ v i1 g ¼
For ORP, in the case of the underlying distribution is the twoparameter Weibull distribution, the cumulative distribution function (cdf) of x for a single system is
Fðxi Þ ¼ 1
j¼1 i¼1
j¼1
We consider a component/system which under goes instantaneous general repairs after a failure. The repair time is assumed to be negligible so that the failure can be viewed as point processes. Weibull distribution is one of the most commonly used distribution in reliability analysis for repairable systems because of its flexibility and applicability to various failure processes. Therefore, we select Weibull distribution as the underlying distribution in this study. Let t be the successive failure times, and x be the time between failures (TBFs). Now, consider k repairable systems, observed from the starting time t0 to the end time Tj (j = 1, 2, . . . , k). The successive failure times of the jth system are denoted by t0 ; t1;j ; t 2;j ; . . . ; t nj ;j ; T j , and the TBF of the jth system are denoted by x0 ; x1;j ; x2;j ; . . . ; xnj ;j ; T j tnj ;j : nj is the number of failures and Tj is the truncated time for the jth system. ti,j is the failure time of the jth system at the ith failure, xi,j is the TBF of the jth system between the (i 1)th and the ith failure. For convenience, we define t0 = 0, x0 = 0.
expðkxbi Þ;
nj k k X k X X X nj ðln k þ ln bÞ þ ðb 1Þ ln t i;j k ðT bj Sbj Þ
Lðxi;j ; T j tnj ;j ; k; b; qÞ ¼
k Y
n h io knj bnj exp k ðT j tnj ;j þ qt nj ;j Þb ðqt nj ;j Þb
j¼1
ð5Þ
x
i¼1 i;j
Estimator of b can be found from the first equation of the system of Eqs. (5) by an iteration method, then substituting it into the second equation, one can get another estimator of k. 2.2. Power law non-homogeneous Poisson process
nj n h io Y ðxi;j þ qt i1;j Þb1 exp k ðxi;j þ qt i1;j Þb ðqt i1;j Þb
! ð13Þ
i¼1
Thus, the log-likelihood function is
ln L ¼
k k h i X X nj ðln k þ ln bÞ k ðT j t nj ;j þ qtnj ;j Þb ðqtnj ;j Þb j¼1
j¼1 nj k X X
k
½ðxi;j þ qti1;j Þb ðqt i1;j Þb
j¼1 i¼1
Similarly, consider k systems. S and T denote the starting time and the end time of observation accounting for censoring. The log-likelihood function of PL-NHPP model is
þ ðb 1Þ
nj k X X j¼1 i¼1
lnðxi;j þ qt i1;j Þ
ð14Þ
1130
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
Taking the derivatives of Eq. (14) with respect to k, b and q, and setting them equal to zero, we can get k k X 1X nj ½ðT j tnj ;j þ qt nj ;j Þb ðqt nj ;j Þb k j¼1 j¼1
nj k X X ½ðqti1;j þ xi;j Þb ðqti1;j Þb ¼ 0
ð15Þ
j¼1 i¼1
k X k ½ðT t nj ;j þ qtnj ;j Þb lnðT tnj ;j þ qtnj ;j Þ ðqtnj ;j Þb lnðqt nj ;j Þ j¼1 nj k X X ½ðqti1;j þ xi;j Þb lnðqti1;j þ xi;j Þ ðqti1;j Þb lnðqti1;j Þ þk j¼1 i¼1 n
j k k X X 1X nj lnðqt i1;j þ xi;j Þ ¼ 0 b j¼1 j¼1 i¼1
ð16Þ
greater than zero. According to the GRP model, there are five cases for q. For q = 0, a system is restored to an ORP (perfect maintenance) condition, while for q = 1, the system corresponds to a NHPP (minimal maintenance) state. The case of 0 < q < 1 corresponds to the intermediate ‘‘better than old but worse than new’’ repair assumption. For the cases where q > 1, the system is in a condition of worse than old. Similarly, the case where q < 0 suggests the system is restored to a state of better than new. Therefore, no constraints are added to variable q. So, it is possible to define the problem of solving the complex system of three equations as a nonlinear constrained programming problem. In this research, the objective function is negative loglikelihood function that is related to model parameters, and the constrained variables are model parameters. They can be expressed by
minð ln LÞ
ð20Þ
s:t: k > 0; 0 < b < 10 ðb 1Þ
nj k X X j¼1 i¼1
nj k X X t i1;j t i1;j ½ðqti1;j þ xi;j Þb1 kb ðqti1;j þ xi;j Þ j¼1 i¼1
ðqt i1;j Þb1 kb
k X
t nj ;j ½ðT j t nj ;j þ qtnj ;j Þb1 ðqtnj ;j Þb1 ¼ 0
j¼1
2.3.2. Generalized renewal process II The GRP II model assumes that repairs could remove all of the damage accumulated up to the last failure time. With this assumption the virtual age of a system after the ith repair is:
ð17Þ We can see that parameter estimation for GRP I model based on the MLE method leads to a complex system of three equations that should be solved simultaneously. Unfortunately, the system of equations mentioned above does not have a closed-form mathematical expression. It is difficult to solve it directly. Therefore, we use a nonlinear constrained programming approach to estimate parameters of GRP model. Mathematical model of finding minimum value of nonlinear constrained programming is
8 min f ðxÞ > > < hi ðxÞ ¼ 0; ði ¼ 1; 2; . . . ; uÞ > g ðxÞ P 0; ðj ¼ 1; 2; . . . ; v Þ > : j x2D
ðkÞ
pðx ; M k Þ ¼ f ðx Þ þ M k
( u X i¼1
i X qimþ1 xm
Therefore, the conditional pdf of GRP model II is
f ðxi ; q; k; bjv i1 Þ ¼ kbðxi þ
i1 X
qim xm Þb1 expfk½ðxi
m¼1
þ
i1 X
qim xm Þb ð
m¼1
ð18Þ
ðkÞ
hi ðx Þ
2
) v X 2 ðkÞ ; þ minð0; g j ðx ÞÞ
i1 X qim xm Þb g; i
m¼1
ð22Þ
The corresponding likelihood function of the TBF for k repairable systems with time truncation is
Lðxi;j ; T j t nj ;j ; k; b; qÞ 8 2 0 !b nj k < Y X n n nj mþ1 @k j b j exp k4 T j t n þ ¼ q x m;j j;j : m¼1 j¼1 0 !b 3 9 !b1 nj nj i1 = Y X X nj mþ1 im 5 @ q xm;j xi;j þ q xm;j ; i¼1 m¼1 m¼1 8 2 !b !b 3911 = < i1 i1 X X exp k4 xi;j þ qim xm;j qim xm;j 5 AA ; : m¼1 m¼1 ð23Þ
j¼1
ðk ¼ 1; 2; . . .Þ
ð21Þ
m¼1
¼ 1; 2; . . . ; n
where f(x) is objective function, x is variable, hi(x) is equality constrained function, gj(x) is inequality constrained function, and D is feasible region. Exterior penalty method can convert constraint problem into a series unconstraint problem by using reconstruction penalty function, and then solve it with unconstraint optimization method. So, this method is also called sequential unconstrained minimizing method. The penalty function is ðkÞ
v i ¼ qðv i1 þ xi Þ ¼ qi x1 þ qi1 x2 þ þ qxi ¼
ð19Þ
where Mk is the penalty factor of the kth calculating. The detailed procedures are given as follows:
The log-likelihood function is
ln L ¼
k X nj ðln k þ ln bÞ j¼1
(1) Set k = 1, give an initial point x(1), the first penalty factor M1 > 0, amplifying coefficient C > 1, allowable error e > 0. (2) Solve unconstraint minimal value of the penalty function p(x(k), Mk). (3) Estimate the calculating accuracy. At point x(k), if the penalty P term Mk vj¼1 ½minð0; g j ðxðkÞ ÞÞ2 < e, then stop calculation, and get an approximate minimal solution x(k) to problem; Otherwise, set Mk+1 = CMk, k = k + 1, go back step (2). To constrain variables, set b be between the values of 0 and 10, since almost all expected solutions are in this range. k can be set
2 !b !b 3 nj nj k X X X n mþ1 n mþ1 j j q xm;j q xm;j 5 k 4 T j t nj ;j þ m¼1
j¼1
m¼1
2 !b !b 3 nj k X i1 i1 X X X im im 4 xi;j þ q xm;j q xm;j 5 k m¼1
j¼1 i¼1
þ ðb 1Þ
nj k X X j¼1 i¼1
ln xi;j þ
m¼1 i1 X m¼1
! im
q
xm;j
ð24Þ
1131
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
We also use nonlinear constrained programming approach to estimate parameters of the GRP model II. 2.3.3. Validation of the generalized renewal process model I and II We use Monte Carlo simulation method to validate the GRP models I and II. For multiple systems, from Eq. (9), we can get the sampling formula of the Weibull distribution for the TBF
n o exp k½ðxi;j þ v i1;j Þb v bi1;j ¼ U i;j
ð25Þ
where Ui,j is random number of uniform distribution on the interval [0, 1]. For GRP I model, the TBF can be got by
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b v bi1;j ln Ui;j k rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b ¼ qti1;j þ ðqti1;j Þb ln U i;j k
xi;j ¼ v i1;j þ
4. Results and discussion
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 1 t i;j ¼ t i1;j þ xi;j ¼ ð1 qÞti1;j þ ðqt i1;j Þb ln U i;j ði k ¼ 1; 2; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ
Table 1 Failure data of single machine tool with time truncated for case 1.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b ¼ ln U i;j ðt0 ¼ 0Þ k
ð28Þ
For GRP II model, the TBF can be got by i1 X
qim xm;j þ
m¼1
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xi1 b 1 b im x q ln U i;j m;j m¼1 k
ð29Þ
where i1 X
qim xm;j ¼
m¼1
¼
i1 X
qim t m;j
m¼1
i1 X
i2 X
q
t m;j
m¼1
¼
i2 X
No. of failures
Failure times
No. of failures
Failure times
1 2 3 4 5 6 7 8
50.29 531.42 914.70 963.27 1029.19 1488.78 1953.16 2380.42
9 10 11 12 13 14 15
2527.76 2839.82 3371.14 3562.01 3649.22 3801.43 4140.00
i1 X qim tm1;j
m¼1 im
Using our proposed method, the results of ORP, PL-NHPP, GRP I and GRP II for three real failure data sets of NC machine tools are all given in Table 4.
ð27Þ
where
xi;j ¼
We illustrate the applications of the proposed approach with three real cases for reliability analysis of NC machine tools. Failure times of NC machine tools are always time truncated data coming from field because of the limitation of time and money. Therefore, three real field failure data sets of NC machine tools are all time truncated data (in cumulative operating hours). Case 1: These data, given in Table 1, contain 15 failure times with time truncation from a NC machine tool and the last one is right truncated time. Case 2: 29 time truncation field failure data of another NC machine tool were collected, as shown in Table 2. Case 3: Table 3 gives 28 failure time data that are from four identical machine tools, and the last one of every machine tool is right truncated time taken from Zhang, Jia, and Shen (2005). Unlike cases 1 and 2, these failure data are multiple system failure data.
ð26Þ
So the sampling formula of failure times can be expressed as follows:
t 1;j
3. Numerical examples
im1
q
Table 2 Failure data of single machine tool with time truncated for case 2.
t m;j
m¼0
qim t m;j þ qt i1;j
m¼1
i2 X
qim1 t m;j qi1 t 0;j
m¼1
¼ qti1;j þ ðq 1Þ
i2 X
qim1 t m;j
ð30Þ
m¼1
So the sampling formula of failure times can be expressed as follows:
No.
Failure times
No.
Failure times
No.
Failure times
1 2 3 4 5 6 7 8 9 10
27.51 367.52 394.52 395.64 406.75 432.49 514.17 855.57 864.85 953.02
11 12 13 14 15 16 17 18 19 20
1039.36 1357.80 1680.92 1850.55 1926.98 2398.21 2430.61 2517.04 2600.22 2796.49
21 22 23 24 25 26 27 28 29
2867.40 3299.82 3387.57 3468.58 3688.63 3780.33 3862.50 3955.48 4152.00
t i;j ¼ t i1;j þ xi;j ¼ ð1 qÞðt i1;j þ
i2 X
qim1 tm;j Þ
Table 3 Failure data of multiple machine tools with time truncated for case 3.
m¼1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u i2 u X b 1 þ t½qti1;j þ ðq 1Þ qim1 t m;j b ln U i;j k m¼1 ði ¼ 3; 4; . . . ; nj ; j ¼ 1; 2; . . . ; kÞ where we define t 0 ¼ 0; Thus
Pi2
im1 t m;j m¼1 q
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 b 1 > < t 1;j ¼ k ln U i;j > :
t 2;j ¼ ð1 qÞt 1;j þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ðqt1;j Þb 1k ln U i;j
ð31Þ ¼ 0ði ¼ 1; 2Þ.
ð32Þ
1 2 3 4 5 6 7 8 9 10 11
No. 1
No. 2
No. 3
No. 4
16 56 116 196 356 676 1028 1428 1828 2488 2522
60 154 304 469 539
96 192 400 800 872
10 74 274 506 1226 2394 2444
1132
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
Table 4 Estimated parameters and –ln L values with the different models for three cases. Model
ORP GRP I GRP II PL-NHPP
No. 1
No. 2
No. 3
b
k
q
ln L
b
k
q
ln L
b
k
q
ln L
1.523 1.523 1.523 0.960
0.0002 0.0002 0.0002 0.0047
0 0 0 1
92.171 92.171 92.171 93.640
1.013 0.766 0.870 0.870
0.0063 0.0327 0.0199 0.0199
0 0.109 1 1
167.972 167.494 167.691 167.691
1.049 0.569 0.569 0.569
0.0028 0.0947 0.0947 0.0947
0 1 1 1
157.927 152.579 152.579 152.579
From Table 4, we can find that in case 1, a q value of 0 suggests that the machine tool’s mending leaves it in the condition of as good as new. In this case we have an ordinary renewal process; In case 2, based on the value of negative log-likelihood function, the GRP I model is superior to the GRP II model, because model I has a smaller –ln L value than that of model II. The best GRP model for the failure data, as determined by the –ln L, is the model with the lowest –ln L value. A q value of 0.109 suggests that the machine tools’ maintenance leave them in the intermediate condition of better than old but worse than new. In this case we have a GRP I; In case 3, a q value of 1 suggests that four machine tools’ repair leaves them in the condition of as bad as old, and we can use a PL-NHPP model to describe the failure process of these four NC machine tools. On the other hand, it is noticed that the results of the GRP I and II model are the same in cases 1 and 3, respectively. Through observations and simulations, Jacopino, Groen, and Mosleh (2004) and Jacopino (2005) also pointed out that at a low number of renewals there is little difference between the GRP I and II model. However, as the number of renewals increases, the difference between the two models becomes significant because of the variation in the underlying virtual age equations of both models. In case 1, we also notice that the results of the ORP model and the GRP model are the same for this particular data set; In case 3, the results of the PL-NHPP model and the GRP model are also the same. In fact, from Eqs. (4), (6), (14), and (24), we can see that when q = 0, the models for ORP, GRP I and GRP II all have the same form of log-likelihood function corresponding to perfect repair; when q = 1, PL-NHPP, GRP I and GRP II also have the same form of log-likelihood function corresponding to minimal repair. This means that the ORP and the PLNHPP are special cases of the GRP.
There are several factors which affect the speed of convergence, these factors include penalty factor Mk, amplifying coefficient C, allowable error e, and initial point x(1) as well. The proposed method is especially sensitive to the initial point x(1), we can use the estimators of the ORP model or the PL-NHPP model as starting point to speed up the convergence of the solution and get a global optimization values rapidly. The computation time using our method in a 1000 MHz Personal Computer is less than 30 s. Compared with one minute computation time of MC method used by Yanez et al. (2002), our method is more efficient. In order to validate the effectiveness of proposed method for GRP I and II model, Monte Carlo simulation method is used. Note that a large sample size generated by a computer according to the sample formulas of failure times is needed when one uses MC method to obtain close and stable estimations. For GRP model I, the details are shown as follows: (1) In the sampling formula (27), set k, b and q be the results of case 2 for GRP I model; (2) generate different sample size data with truncated time is 4152, respectively; (3) estimate model parameters and calculate –ln L values. The results are shown in Fig. 1, we can see that simulation results including –ln L values are close to real values when sample size is up to 2000. Base on the –ln L values, therefore, we say that our method is effective. The simulation value of q is slightly larger than real value, this is because that in this model, for convenience, we assume restoration factor q are all equal to each other. However, in reality, the degree of repair after a failure is different due
Fig. 1. Simulation results of GRP I for case 2.
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
to the different failure modes and positions, etc. Uematsu and Nishida (1987) and Kijima (1989) extended restoration factor q to the case that it is a random variable taking a value between 0 and 1. In this study, we can regard q as an average value of repair quality. In GRP II, Mettas and Zhao (2005) analyzed 56 failure truncated data and given the estimators of model parameters, we select them as the known values in the sampling formula (31) to validate our proposed method for GRP model II. The results are shown in Fig. 2. It is found that the simulation results are close to real values. Based on a MC simulation, Fig. 3 shows the cumulative number of failures against the time in case 2 for GRP I. It is noticed that the number of failures are different in three models. When q = 0.109, the GRP I model fits the actual data better than the other two models corresponding to q = 0 and q = 1, respectively.
1133
As mentioned in Section 1, there are several methods for GRP I in literature. Among them, the MLE method discussed by Yanez et al. (2002) is based on MC simulation with numerical method, so the computation of this method is quite time-consuming and slow. Other two methods including MCMC sampling method used by Jacopino et al. (2006) and EM algorithm used by Veber et al. (2008) are both have very limited application areas, for the time to first failure must be known and estimated from a large set of data. However, when one uses our proposed method to analyze GRP, the time to first failure distribution and a large set of data are all unnecessary. In addition, three methods mentioned above are all only for the GRP I. In this study, we not only explore using GRP I model but also GRP II model to model and analyze complex repairable systems. Therefore, compared with the other approaches in literature, our GRP numerical approach is simple and efficient.
Fig. 2. Simulation results of GRP II.
Fig. 3. Cumulative number of failures versus time plot.
1134
Z.-M. Wang, J.-G. Yang / Computers & Industrial Engineering 63 (2012) 1128–1134
5. Conclusions In this paper, a nonlinear constrained programming method is proposed based on the conditional Weibull distribution to estimate the parameters of the GRP model with different types of data. To avoid solving complex system of equations, the proposed method uses negative log-likelihood function as objective optimal function and adds inequality constraints to model parameters. The sampling formulas of failure times for GRP models I and II are given and the effectiveness of the proposed method is validated with simulation method. The GRP model considers restoration factor in reliability analysis of repairable systems, therefore, its evaluation results are close to real maintenance circumstance than the ORP model and the PLNHPP model, the latter are only two special cases where restoration factor is 0 and 1 respectively. In case where the amount of failure data is extremely limited, other methods such as a Bayesian approach based on prior knowledge may be used to remedy the uncertainty of estimations. Acknowledgements We would like to thank the anonymous reviewer and the editor for their valuable comments and suggestions, which have greatly enhanced the clarity of the paper. The authors also wish to acknowledge the partial financial support of this study by ‘‘Top Grade NC Machine Tools and Fundamental Manufacturing Equipment’’ National Science and Technology Key Special Project of China (No. 2009ZX04014-022)’’. References Blischke, W., & Murphy, D. (1994). Warranty cost analysis. New York: Marcel Dekker. Jacopino, A. (2005). Generalization and Bayesian solution of the general renewal process for modeling the reliability effects of imperfect inspection and maintenance based on imprecise data. Dissertation for the degree of Doctor of Philosophy, University of Maryland.
Jacopino, A., Groen, F., & Mosleh, A. (2004). Behavioural study of the general renewal process. In Proceedings of the annual reliability and maintainability symposium, 26–29 January 2004 (pp. 237–242). Los Angeles, CA, USA. Jacopino, A., Groen, F., & Mosleh, A. (2006). Modeling imperfect inspection and maintenance in defense aviation through Bayesian analysis of the KIJIMA type I general renewal process (GRP). In Proceedings of annual reliability and maintainability symposium, 23–26 January 2006 (pp. 470–475). Newport Beach, California, USA. Kaminskiy, M. & Krivtsov, V. (2000). G-renewal process as a model for statistical warranty claim prediction. In Proceedings of the annual reliability and maintainability symposium, 24–27 January 2000 (pp. 276–280). Los Angeles, CA, USA. Kaminskiy, M., & Krivtsov, V. (1998). A Monte Carlo approach to repairable system reliability analysis. Probabilistic safety assessment and management. New York: Springer, pp. 1063–1068. Kijima, M. (1989). Some results for repairable systems with general repair. Journal of Applied Probability, 26, 89–102. Kijima, M., & Sumita, N. (1986). A useful generalization of renewal theory: counting process governed by non-negative Markovian increments. Journal of Applied Probability, 23, 71–88. Krivtsov, V. (2000). A Monte Carlo approach to modeling and estimation of the generalized renewal process in repairable system reliability analysis. Dissertation for the Degree of Doctor of Philosophy, University of Maryland. Majeske, K. D. (2007). A non-homogeneous Poisson process predictive model for automobile warranty claims. Reliability Engineering and System Safety, 92, 243–251. Mettas, A. & Zhao, W. (2005). Modeling and analysis of repairable systems with general repair. In Proceedings of the annual reliability and maintainability symposium, 24–27 January 2005 (pp. 176–182). Alexandria, Virginia, USA. Saldanha, P. L. C., Simone, E. A., & Frutoso, P. F. (2001). An application of nonhomogeneous Poisson point processes to the reliability analysis of service water pumps. Nuclear Engineering and Design, 210, 125–133. Uematsu, K., & Nishida, T. (1987). Branching non-homogeneous Poisson process and its application to a replacement model. Microelectronics and Reliability, 27, 685–691. Veber, B., Nagode, M., & Fajdiga, M. (2008). Generalized renewal process for repairable systems based on finite Weibull mixture. Reliability Engineering and System Safety, 93, 1147–1461. Weckman, G. R., Shell, R. L., & Marvel, J. H. (2001). Modeling the reliability of repairable systems in the aviation industry. Computers and Industrial Engineering, 40, 51–63. Yanez, M., Joglar, F., & Modarres, M. (2002). Generalized renewal process for analysis of repairable systems with limited failure experience. Reliability Engineering and System Safety, 77, 167–180. Zhang, Y., Jia, Y., & Shen, G. (2005). Research on model of failure distribution for numerical machine with random ending method. Systems Engineering Theory and Practice, 2, 134–138. in Chinese.