A multivariate CBM model with a random and time-dependent failure threshold

Reliability Engineering and System Safety 119 (2013) 178–185

Contents lists available at SciVerse ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

A multivariate CBM model with a random and time-dependent failure threshold R. Jiang n Faculty of Automotive and Mechanical Engineering, Changsha University of Science and Technology, Changsha, Hunan 410114, China

art ic l e i nf o

a b s t r a c t

Article history: Received 25 January 2013 Received in revised form 14 May 2013 Accepted 14 May 2013 Available online 29 May 2013

In a condition-based maintenance setting, the degradation of an item is usually represented by several condition variables, and they can be combined into a composite condition variable. In this case, the functional failure threshold associated with the composite condition variable is usually not a fixed and known constant. It is an open issue to model the failure threshold and accordingly determine a threshold of preventive maintenance (PM). This paper addresses this issue. The condition variables are combined using a weighted power model, the failure threshold is represented by the Gaussian process model, and the PM threshold is determined by two approaches. Based on the gamma process and stress–strength interference models, the distributions of time to failure and to the PM threshold are derived, respectively. The appropriateness of the approach is illustrated by a real-world example. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Condition-based maintenance Composite condition variable Weighted power model Failure threshold PM threshold

1. Introduction Condition-based maintenance (CBM) is a widely used maintenance strategy for critical and deteriorating items (e.g., see [1]). In simple cases, the degradation process of an item is represented by a single condition variable, and the item fails when the condition variable reaches a certain critical level, which is referred as functional failure threshold. In this case, the failure threshold is often fixed and known (e.g., see [2], [3] and the literature cited therein). The use of a fixed failure threshold is reasonable if the monitored condition variable directly relates to the state of the item (e.g., wear amount or drop amount in a certain performance). However, in most of practical problems, the degradation process is represented by several condition variables, which indirectly relate to the state of the item (e.g., see [4] and [5]). In this case, the condition variables are usually combined into a composite condition variable and its failure threshold is no longer known and fixed rather than random and time-dependent due to the following reasons [6]: (a) the critical threshold value can vary appreciably among users of products, and (b) the level of degradation that will cause a failure is not explicitly known. Jiang and Jardine [2] present a random failure threshold in a proportional hazard model setting. Wang and Coit [6] consider n

Tel.: +86 731 82309122; fax: +86 731 82309621. E-mail addresses: [email protected], [email protected]

0951-8320/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ress.2013.05.023

a random and time-independent failure threshold to assess reliability based on degradation modeling. Usynin et al. [7] examine the issue how variability in the failure threshold affects the reliability prediction in conjunction with cumulative damage models. Wang and Carr [8] and Wang et al. [9] use the exponential distribution with a location parameter and the normal distribution to represent the failure threshold. As a whole, the studies in CBM involving a random failure threshold are limited. When the failure threshold is fixed, an alarm threshold is usually set to trigger a preventive maintenance (PM) action (e.g., see [3]). In this case, the item is definitely at working state when the composite condition variable is at or before the alarm threshold. However, this is not true if the failure threshold is random so that the item can fail at or before the alarm threshold. This raises a new problem, i.e., how to set a PM threshold when the failure threshold is random. In summary, a systematic approach is needed to address the following issues: (a) Fusing multiple condition variables into a composite condition variable in an appropriate manner, (b) Building the failure threshold model of the composite condition variable, which is random and possibly time-dependent, and (c) Determining the PM threshold of the composite condition variable. This paper presents such an approach to address these issues. The proposed approach uses a weighted power model to represent the composite condition variable, and uses the Gaussian

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

process model to represent the failure threshold of the composite condition variable. Based on these two models, the PM threshold is determined based on a cost rate model or the tradeoff BX life approach proposed by Jiang [10]. The resulting multivariate CBM model will be applied to derive the distributions of time to failure and to the PM threshold. The approach is illustrated by a realworld example. The paper is organized as follows. The proposed approach is presented in Section 2. The applications of the CBM model are discussed in Section 3. The results are illustrated in Section 4. The paper is concluded in Section 5. 2. Proposed approach The approach involves the following four steps:


Step 1: Preliminary treatment of data,
Step 2: Development of the composite condition variable

179

expectation approach. Suppose that zij ¼ zþ 0 is a censored observation of Zi. We replace the censored observation zij with an equivalent failure observation z′ij ¼ z0 þ Mðz0 Þ, where Mðz0 Þ is the mean residual value given by Z ∞ 1 ðz−z0 ÞdF i ðzÞ ð3Þ Mðz0 Þ ¼ 1−F i ðz0 Þ z0 For the Weibull distribution, we have (e.g., see [12]) zij0 ¼ ηi Γð1 þ 1=βi Þ

Ga ððz0 =ηi Þβi ; 1 þ 1=βi ; 1Þ 1−F i ðz0 Þ

ð4Þ

where Ga ðx; u; vÞ is the gamma distribution function with the shape parameter u and scale parameter v and evaluated at x. In such a way, the censored observations are transformed into equivalently complete observations. In the following discussion, we assume that the data have been transformed into the complete data and still use zij to denote zij0 .

model,


Step 3: Development of the failure threshold model, and
Step 4: Determination of the PM threshold. Specific details are presented in the following four subsections. 2.1. Preliminary treatment of data The starting point is a set of observations of times to failure T and condition variables Zi (1 ≤ i ≤m) at failure, given by ðt j ; zij ; 1 ≤i ≤m; 1 ≤j ≤nÞ

ð1Þ

2.1.3. Dimensionless transformation of data Assume Z i ≥0; otherwise, a transformation can be performed to meet this inequality. The condition variable Zi is measured using a specific scale and can be increasing or decreasing with time. As such, a transformation is needed to make all the condition variables increase (or decrease) and to delete the effect of the scale. In this paper, we transform condition variables to the increasing case. If the dataset is complete, the common normalization transformation for the case of Zi being increasing is given by

where j indicates the j-th failure event; zij is the realization of Zi at tj and can be either complete or censored. The preliminary treatment of data deals with the following three issues:

X i ¼ ðZ i −Z i;min Þ=ðZ i;max −Z i;min Þ


to evaluate the failure prediction capability of each condition

X i ¼ ðZ i;max −Z i Þ=ðZ i;max −Z i;min Þ





variable based on their variability, to transform censored observations into equivalently complete observations so as to facilitate further analysis, and to transform the data into a dimensionless form so as to delete the effect of scales of the condition variables. Specific details are outlined as follows.

2.1.1. Failure prediction capability of a condition variable Let Zi,f denote the value of Zi at failure, which is a random variable. If the coefficient of variation (CV) of Zi,f is small, then Zi has a good failure prediction capability (e.g., see [2,5] and [11]). To make the notation simple, we will no longer differentiate Zi from Zi,f in the remainder of the paper. Let ρi denote the CV of Zi. To evaluate ρi, we fit the data subset (zij ; 1 ≤j ≤n) to an appropriate distribution Fi(z) using the maximum likelihood method. Without the loss of generality, we consider the Weibull distribution with the shape parameter βi and scale parameter ηi. The CV of Zi is given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Γð1 þ 2=βi Þ ρi ¼ −1 ð2Þ Γ 2 ð1 þ 1=βi Þ where Γð:Þ is the gamma function. As such, the failure prediction capabilities of Zi (1 ≤ i ≤m) can be evaluated based on the values of ρi's. If Zk has the smallest CV, then it is the principal condition variable. 2.1.2. Transformation of censored data To facilitate further analysis, we need to transform a censored observation into an equivalent failure observation using an

ð5Þ

and the transformation for the case of Zi being decreasing is given by ð6Þ

where Zi,min and Zi,max are the maximum and minimum values of the dataset. It is noted that Xi is dimensionless. In this paper, we will not use (5) or (6) due to the following two reasons: (a) The weighted power model to be presented in (11) requires X i 40; and (b) If the dataset is incomplete, Zi,max or Zi,min is uncertain. We revise (5) and (6) to obtain the transformation with desired properties. For the case that Zi is increasing, the lower limit of Zi (denoted as ziL) is usually known. To revise (5), we use ziL to replace Zi,min, and ziα to replace Zi,max, where ziα ( ¼ F −1 i ðαÞ) is the α-fractile of Zi. The α should be sufficiently large (e.g., 0.99) so that zia 4 Z i for all or most of the observations. However, zia o Z i;max is allowed. Usually, we have ziL ¼ 0. In this case, (5) is revised as X i ¼ Z i =ziα

ð7Þ

If the upper limit of Zi is known, we use the upper limit to replace ziα. For the case that Zi is decreasing, the range of Zi is usually known. To revise (6), we use the lower [upper] limit to replace Z i;min [Z i;max ]. For the lifetime data, the lower limit is usually equal to zero. As such, the lifetime data are transformed by τ ¼ t=t α . As such, the data given by (1) can be written into the following form: ðτj ; xij ; 1 ≤i ≤m; 1 ≤j ≤nÞ

ð8Þ

180

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

2.2. Composite condition variable model

2.3. Failure threshold associated with composite condition variable

2.2.1. Traditional approach There are various approaches to fuse multiple condition variables into a single composite condition variable (e.g., see [4,5,13] and the literature cited therein); and traditional approach is through a linear combination of condition variables given by

The parameters of the composite condition variable model are determined through minimizing the CV of Y so that Y (at failures) varies in small range and is not sensitive to the item age t. As such, for a given t, it is reasonable to assume that Y(t) (at failures) follows the normal distribution with mean μ(t) and standard deviation s(t). The mean and standard deviation can change with age. The change rates are expected to be small so that we can assume:

m

Y ¼ ∑ ai X i ; ai ≥0; ak ¼ 1 i¼1

ð9Þ

where k is associated with the principal condition variable Xk (corresponding to Zk). The linear model given by (9) has only m−1 unknown parameters. Using the data given by (8) to (9), we can obtain a sample of Y (at failures) given by ðτj ; yj ; 1 ≤j ≤nÞ

ð10Þ

μðtÞ ¼ μ0 þ s1 t; sðtÞ ¼ s0 þ s2 t; s0 ≥0

ð15Þ

Eq. (15) implies that the failure threshold is represented by the Gaussian process with linear mean and standard deviation. In other words, for a given t, the composite condition variable Y(t) at failure is a normal random variable with the distribution function: FðyjtÞ ¼ Φðy; μðtÞ; sðtÞÞ;

ð16Þ

Let μi and sy denote the sample mean and sample standard deviation of Y, respectively. The CV of Y is given by ρy ¼ sy =μy . To make Y have good failure prediction capability, the parameters are determined so that ρy achieves its minimum or 1=ρy achieves its maximum.

where Φðy; μ; sÞ is the normal cdf with the location parameter μ and scale parameter s. As a result, the failure threshold model is fully specified by (15) and (16). The parameter set θ ¼ ðμ0 ; s1 ; s0 ; s2 Þ can be estimated using the maximum likelihood method for the data given by (10) with the log-likelihood function:

2.2.2. Weighted power model It is noted that (9) can be written as

lnðLÞ ¼ ∑ ln½ϕðyj ; μðt j Þ; sðt j ÞÞ

m

Y ¼ K ∑ wi X i ¼ KM 1 i¼1

n

j¼1

ð11Þ

m where K ¼ Σ m i ¼ 1 ai , wi ¼ ai =K, and M 1 ¼ Σ i ¼ 1 wi X i , which is the weighted arithmetic mean. Eq. (11) implies that the composite condition variable is proportional to the weighted arithmetic mean of the constituting condition variables. Usually, the condition variables contribute to the degradation of an item in a more complex way. For example, if a certain condition variable has a large value and the other variables have small values, the item can fail though their linear combination is possibly not too big. To represent this kind of complex relation, we assume that the composite condition variable is proportional to the generalized mean of the constituting condition variables. This yields the weighted power model given by

Y ¼ KM p ðX 1 ; :::; X m Þ

ð12Þ

where M p ðX 1 ; :::; X m Þ is the generalized mean (e.g., see [14]) given by m

M p ðX 1 ; :::; X m Þ ¼ ð ∑ wi X pi Þ1=p ; p∈ð−∞; ∞Þ i¼1

ð13Þ

where X i 4 0, wi ≥0 and Σ m i ¼ 1 wi ¼ 1. When p ¼ −1, 0, 1 and 2, Mp is the weighted harmonic, geometric, arithmetic and quadratic means, respectively. Specially, M −∞ ¼ minðX 1 ; :::; X m Þ when p-−∞, and M ∞ ¼ maxðX 1 ; :::; X m Þ when p-∞. In other words, the generalized mean with a big [small] p emphasizes big [small] values among (X 1 ; :::; X m ). Let ai ¼ K p wi , (12) can be written as m

Y ¼ ð ∑ ai X pi Þ1=p i¼1

ð14Þ

Letting ak ¼ 1, the weighted power model given by (14) has m unknown parameters (i.e., p, ai, 1 ≤i ≤ m, i≠k). The approach outlined in Section 2.2.1 can be used to estimate these parameters. The weighted power model reduces into the linear model (11) when p ¼1. In a CBM setting, it can be reasonably expected to have p⪢1. In this sense, the weighted power model outperforms the linear model, as illustrated later.

ð17Þ

where ϕð:Þ is the normal pdf, and t j ¼ τj t α . In this paper, we use the Solver of Microsoft Office Excel to perform the maximization (or minimization) operation. 2.4. PM threshold Let Yf(t) denote the random failure threshold and YPM ¼ω(t) denote the PM threshold. As mentioned earlier, the item can fail at or before the PM threshold. For a given t, a preventive replacement is performed if ωðtÞ ≤YðtÞ o Y f ðtÞ, where the event “YðtÞ o Y f ðtÞ” means “the item is working” at t. Clearly, a large PM threshold value results in poor reliability and a small value implies large loss of useful life. The PM threshold must be appropriately chosen so that a good tradeoff between the reliability and life can be achieved. If the preventive and failure replacement costs can be appropriately specified, the PM threshold can be obtained through minimizing the cost rate associated with the classic age replacement policy [15]. If it is difficult to accurately specify the cost parameters, the tradeoff BX life approach proposed by Jiang [10] can be used to obtain the PM threshold. The basic idea and main results of the tradeoff BX life approach are briefly outlined as follows. 2.4.1. Tradeoff BX life Let f ðtÞ, FðtÞ, RðtÞ and rðtÞ denote the density, distribution, reliability and failure rate functions of a random life T, respectively. The BX life is defined by BX ¼ F −1 ðxÞ; x ¼ X=100∈ð0; 1Þ

ð18Þ

An item will be preventively replaced at age BX if its life is larger than BX; otherwise it is correctively replaced at the time when it fails. The probability of a preventive replacement is 1−x and the probability of a corrective replacement is x. The expected useful life of the item is given by Z BX μðxÞ ¼ RðtÞdt ð19Þ 0

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

181

3.1. Prediction of degradation level

0.3 0.25

For a given t, assume that the composite condition variable Y(t) follows a gamma distribution with the cdf Ga ðy; u; vÞ and the pdf g a ðy; u; vÞ. The stationary gamma process is defined by Ga ðy; u; vÞ with u ¼ λt, where λ ( 4 0) represents the degradation rate. The model has two parameters: λ and v. Using the data given by (10), the parameters can be estimated using the maximum likelihood method with the log-likelihood function given by

x*

0.2 0.15 0.1 0.05 0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

CV

n

lnðLÞ ¼ ∑ ln½g a ðyi ; λt i ; vÞ

Fig. 1. Plot of xn versus CV for the normal distribution.

j¼1

The useful life contributed by the preventively replaced item equals μp ðxÞ ¼ ð1−xÞBX ; and the useful life contributed by the failed item equals μf ðxÞ ¼ μðxÞ−μp ðxÞ. It is desired that μp ðxÞ is large and μf ðxÞ is small (implying high reliability). As a result, the tradeoff BX life is defined by the value of BX when μp ðxÞ achieves its maximum. Letting dμp ðxÞ=dx ¼ 0 and after some simplifications, the tradeoff BX life is specified by the solution of rðBX ÞBX ¼ 1

ð20Þ

Let xn or BnX denote the solution. For the Weibull, gamma and lognormal life distributions, Jiang [10] shows that xn increases with CV. This is also true for the normal life distribution as shown in Fig. 1. It clearly shows that the tradeoff BX life corresponds to a high reliability if CV is small.

2.4.2. PM threshold of Y The CBM is applicable for the situations where the item is gradually deteriorating with age, implying the item age contributes to the item failure to some extent. To maintain the scale of the age unchanged and further improve the failure prediction capability, we combine the item age and composite condition variable into an equivalent age V by a linear combination: T V ¼ þ γY; γ 4 0 tα

ð21Þ

For a given value of γ, from (10), we can obtain a sample of V: (vj ; 1 ≤ j ≤n). Once more, the parameter γ is determined by minimizing the sample CV of V. When CV is small (this is the case in the current context), the normal distribution with the parameters μv and sv is an appropriate model for fitting the sample of V. Define a critical value (denoted as v0) of V. When V ¼v0, a PM action is triggered. The value of v0 is determined from the classic age replacement policy or the tradeoff BX life approach associated with the distribution of V. Once v0 is specified, from (21) the PM threshold is given by v0 −τ Y PM ðtÞ ¼ ωðtÞ ¼ γ

ð22Þ

This is a decreasing straight line in the Y–t plane. As such, the item age has an upper limit t u ¼ v0 t α , which is usually very large.

ð23Þ

Suppose that we have an observation Y(t)¼ y0 ( o ωðt 0 Þ) at t0. For t 4 t 0 , the stationary gamma process has independent increment ΔY ¼ YðtÞ−y0 , whose pdf and cdf are given respectively by gðΔyÞ ¼ g a ðΔy; λðt−t 0 Þ; vÞ; GðΔyÞ ¼ Ga ðΔy; λðt−t 0 Þ; vÞ

ð24Þ

As such, the degradation level at t is given by YðtÞ ¼ y0 þ ΔY. 3.2. Prediction of residual life The distribution of time to failure is given by F f ðtÞ ¼ PrðYðtÞ≥Y f Þ, where Yf is the failure threshold given by (16) and the event “YðtÞ≥Y f ” means “the item is at failure state” at t. In terms of the stress–strength interference model, Y(t) can be viewed as the stress and Yf can be viewed as the strength. As a result, the distribution of time to failure is given by Z ∞  ϕðy; μðtÞ; sðtÞÞ  1−Ga ðy−y0 ; λðt−t 0 Þ; vÞ dy F f ðtÞ ¼ ð25Þ 1−Φðy ; μðtÞ; sðtÞÞ y0 0 For a set of values of t, a numerical integration method can be used to evaluate Ff(t). This generates a set of data pairs (t k ; F f ;k ; 1 ≤k ≤N). The data can be fitted to a standard distribution with the location parameter t0 using the least square method, and the parameters can be estimated by minimizing the following: N

SSE ¼ ∑ ½F f ðt k −t 0 ; θÞ−F f ;k 2 k¼1

ð26Þ

The life distribution and statistic characteristics (e.g., mean and variance) of the population without condition information can be derived by letting t0 ¼0 and y0 ¼ 0. 3.3. Prediction of time to the PM threshold In a similar argument, the distribution of time to the PM threshold is given by F p ðtÞ ¼ PrðYðtÞ≥Y PM Þ ¼ 1−Ga ½ωðtÞ−y0 ; λðt−t 0 Þ; vÞ

ð27Þ

The distribution of time to the PM threshold for the population without condition information is directly given by (27) through letting t0 ¼0 and y0 ¼0. It can be also approximated by a standard distribution with the location parameter t0 using the least square method outlined above.

4. Illustration 3. Applications 4.1. Background and data Given an observation yt ( o ωðtÞ) of Y at time t, we need to infer the residual life of the item and time to the PM threshold. This involves inferring the future degradation level of a surviving item. The stationary gamma process model (e.g., see [16]) can be used for this purpose. We first look at the prediction of degradation level in Section 3.1, and then address the predictions of residual life and time to the PM threshold in Sections 3.2 and 3.3, respectively.

This illustration comes from [4] and deals with the degradation of reduction cells for aluminum production. Three condition variables of characterizing the degradation process of reduction cells are the iron contamination of the molten aluminum from cracking of the carbon lining Z1 in percent, the horizontal distortion of the steel box Z2 in inches and the amount of cathode drop Z3 in inches. Cells usually fail when the deterioration reaches

182

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

Table 1 Data on failure age and condition variables.

Table 3 Correlation coefficients between failure age and condition variables.

j

tj

z1j

z2j

z3j

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

1540 1415 660 999 1193 1006 869 1035 797 296 775 1424 1169 1500 728 670 841

0.61 1.95 1.20 1.48 0.44 1.22 1.33 2.69 2.51 43.00 2.15 3.60 0.46 0.33 1.15 2.83 4.30

0.60 4.80 3.00 3.00 4.00 4.00 4.70 1.30 2.80 1.50 4.30 3.00 4.30 3.80 3.30 3.75 2.80

0.66 0.67 0.74 0.74 0.72 0.88 0.66 0.68 0.76 0.72 0.74 0.97 0.88 0.76 0.79 0.76 0.70

Table 2 Distributional parameters and characteristics of failure age and condition variables.

β η μ s ρ t0.99, z0.99

T

Z1

Z2

Z3

3.3198 1109.7 995.8 330.5 0.3319 1758.0

1.5726 2.1139 1.8984 1.2342 0.6501 5.5826

3.2057 3.5972 3.2221 1.1034 0.3425 5.7924

8.4081 0.7939 0.7495 0.1061 0.1416 0.9520

a certain level. Table 1 shows the failure age tj in days of service and the corresponding values of the condition variables for 17 nominally identical reduction cells that were operated to failure under uniform conditions. Whitmore et al. [4] use a bivariate Wiener process model to fit the data. The model requires that each condition variable value at t¼0 should equal zero. As such, Z3 is excluded from the analysis since it is considered to have a non-zero initial value. They consider the linear model given by (9), and estimate the model parameters by maximizing a squared correlation coefficient between the composite condition variable and the degradation process. The “true” degradation process is unobservable so that the correlation coefficient has to be estimated from the failure times and condition variable values at failure. Their approach is applicable only for a complete sample, and hence the censored observation for item 10 is approximately viewed as a complete observation.

T Z1 Z2 Z3

T

Z1

Z2

Z3

1 −0.4415 0.0782 0.1319

−0.4415 1 −0.3030 0.0277

0.0782 −0.3030 1 0.2187

0.1319 0.0277 0.2187 1

Table 4 Dimensionless data on failure age and condition variables. j

τj

x1j

x2j

x3j

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0.8760 0.8049 0.3754 0.5683 0.6786 0.5723 0.4943 0.5887 0.4534 0.1684 0.4409 0.8100 0.6650 0.8533 0.4141 0.3811 0.4784

0.1093 0.3493 0.2150 0.2651 0.0788 0.2185 0.2382 0.4819 0.4496 0.7093 0.3851 0.6449 0.0824 0.0591 0.2060 0.5069 0.7702

0.1036 0.8287 0.5179 0.5179 0.6906 0.6906 0.8114 0.2244 0.4834 0.2590 0.7423 0.5179 0.7423 0.6560 0.5697 0.6474 0.4834

0.6932 0.7037 0.7773 0.7773 0.7563 0.9243 0.6932 0.7142 0.7983 0.7563 0.7773 1.0189 0.9243 0.7983 0.8298 0.7983 0.7353

This implies that the approach of minimizing CV is more appropriate than the approach of maximizing the correlation coefficient.

4.2.2. Transformation of censored observation The value of Z1 of Item 10 is a censored observation. From (4), we have the mean residual value M(3) ¼0.96. As such, the observation “43.00” is replaced with 3.96.

4.2.3. Dimensionless transformation of data All the three condition variables are positive and increasing with time, and hence (7) can be used to make the data dimensionless. The values of ziα with α ¼ 0.99 are shown in the last row of Table 2. The dimensionless transformed data are shown in Table 4. It is noted that only one datum is slightly larger than 1. As mentioned earlier, this is allowed.

4.2. Preliminary treatment of data 4.3. Composite condition variable model 4.2.1. Failure prediction capabilities of age and condition variables Since the dataset for Z1 is incomplete, the sample CV cannot be evaluated directly. We first fit a specific dataset to a distribution and then evaluate the CV from the fitted model. We fit the data of T and Zi (1 ≤ i ≤3) in Table 1 to the Weibull distributions, and the results are shown in Table 2. As seen from the sixth row of Table 2, among the three condition variables, Z3 has the smallest CV (i.e., the best failure prediction capability), and hence is chosen as the principal condition variable (i.e., k ¼3). For comparison purpose, we calculate the correlation coefficients between the failure ages and the condition variables and the results are shown in Table 3. As seen, the absolute value of the correlation coefficient between T and Z1 is relatively large, but the failure prediction capability of Z1 is poor since its CV is the largest.

Using the approach outlined in Section 2.2, we obtained the parameters of the composite condition variable for the linear and weighted power models, which are shown in the 2nd and 4th row of Table 5, respectively. As expected, p⪢1. Namely, the weighted power model emphasizes big values among the condition variables. The values of 1=ρy for the two models are 9.05 and 10.51. The last row of Table 5 shows the inverse CV values of individual condition variables, which are obtained from the transformed data. Compared with 1=ρ3 ¼ 8.8703, the linear model has a minor improvement to the failure prediction capability, and the weighted power model has a significant improvement in failure prediction capability due to 1=ρy ⪢1=ρ3 .

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

Table 5 Parameters of Y and contributions of Zi to Y. Z2

Z1 p¼1 p=9.8842

ai ci, %

0.0785 3.11

0.0685 4.46

ai ci, % 1/ρi

1.8932 6.10 1.5106

1.0325 12.74 2.7039

183

Table 6 Parameters of the mean and envelope lines. Z3

p

1 92.44

9.05

1 81.16 8.8703

10.51

cc

0.9904

b d

Mean

Lower

Upper

0.7985 0.0483

0.8037 −0.1472

0.7585 0.3242

0.9974

Table 7 Values of v0 obtained from different approaches. 1.2 1

Approach

Cost ratio

v0

x

v0tα

Cost model

cf/cp ¼4 cf/cp ¼6

10.7877 10.5279

0.0154 0.0087

18,964 18,508

11.3667

0.0476

19,982

y

0.8 0.6

Tradeoff BX

0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

τ Fig. 2. Plots of y(τ), mean and envelope lines.

μ0 ¼ 0:7997; s1 t α ¼ 0:0451; s0 ¼ 0; s2 t α ¼ 0:1091

The contribution and contribution rate of Zi to Yf can be evaluated respectively by n

m

j¼1

l¼1

C i ¼ ai ∑ xpij ; ci ¼ C i = ∑ C l

ð28Þ

The contribution rates are shown in the 3rd and 5th row of Table 5 for the two models. As seen, Z3 makes a major contribution to failure. As a result, a great effort in condition monitoring should be given to Z3. The second largest contribution comes from Z2, and Z1 has the smallest contribution. The last column of Table 5 shows the correlation coefficients (denoted as cc) between 1=ρi and ci, which are very close to 1. This implies that a small CV indicates a large contribution or good failure prediction capability. This further confirms that the approach of minimizing CV is appropriate and the weighted power model significantly outperforms the linear model in terms of the failure prediction capability.

4.4. Failure threshold model 4.4.1. Checking linear assumptions of μðtÞ and sðtÞ Fig. 2 shows the plot of y versus τ (i.e., those dots), and three straight lines (i.e., mean, lower and upper envelope lines). The straight line model is lðτÞ ¼ b þ dτ

ð29Þ

The mean line can be easily obtained by regression and the parameters are shown in the 2nd column of Table 6. The lower envelope line can be obtained by minimizing n

SSD ¼ ∑ ½yðτj Þ−lðτj Þ2 j¼1

4.4.2. Determination of μðtÞ and sðtÞ Using the maximum likelihood method, we have the parameter estimates of the Gaussian process model:

ð30Þ

subject to the constraint lðτj Þ ≤yðτj Þ. Similarly, the upper envelope line can be obtained by minimizing SSD given by (30) subject to the constraint lðτj Þ≥yðτj Þ. Their parameters are shown in the last two columns of Table 6. As seen from the table and figure, the mean line slightly increases and the range between the upper and lower envelope lines linearly increases with τ. This validates the assumptions about μðtÞ and sðtÞ.

ð31Þ

Since tα is large, the values of s1 and s2 are actually small. 4.5. PM threshold of Y 4.5.1. Distribution of equivalent age V Using the approach outlined in Section 2.4.2, we obtained the parameter of the equivalent age model, γ ¼ 15:4685. According to (21), the contribution of the composite condition variable to the item failure can be calculated as below n

γ ∑ yj cy ¼

j¼1

n

ð32Þ

∑ ðt j =t α þ γyj Þ

j¼1

For the current example, cy ¼ 95.76%, implying that the main contribution to the item failure comes from the composite condition variable. Fitting the sample of V to the normal distribution yields μv ¼ 13:3402 and sv ¼ 1:2167 with 1=ρv ¼ 10:96. Since 1=ρv 4 1=ρy (¼10.51), the failure prediction capability gets further improved by the equivalent age model. 4.5.2. Value of v0 from the cost model Let cf [cp] denote the failure [preventive] replacement cost. For the distribution of V and a given value of cf/cp, the optimal value of v0 can be obtained from the age replacement policy. For cf/cp ¼4 and 6, the second and third rows of Table 7 shows the optimal values of v0 and the corresponding failure probabilities at v0. The last column of Table 7 shows the age upper bound t u ¼ v0 t α . As mentioned earlier, tu is very large. Fig. 3 shows the PM thresholds that correspond to different cost ratios. The lowermost PM threshold corresponds to cf/cp ¼6. 4.5.3. Value of v0 from the tradeoff BX life approach If it is hard to specify the cost ratio, the tradeoff BX life approach can be used. For the current example, the results obtained from this approach are shown in the last row of Table 7. The solution is equivalent to the solution of the cost model with cf =cp ≈2:14. The corresponding PM threshold (i.e., the uppermost one) is also shown in Fig. 3.

184

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

1.2 1

y(t)

0.8 0.6

PM threshold

0.4 0.2 0 0

500

1000

1500

2000

t

4.6.2. Distribution of time to the PM threshold for population For t 0 ¼ 0, y0 ¼ 0 and v0 ¼ 11:36, from (27) we obtained a set of values of F p ðtÞ, which are shown in Fig. 4 (i.e., those dots). Fitting these values to the Weibull distribution, we obtained the parameters shown in the last row of Table 8. As seen from Fig. 4, the fitted Weibull distribution is in good agreement with the data points. The last column of Table 8 shows that 75.06% of the useful life can be utilized if the item is preventively replaced at the PM threshold obtained from the trade-off BX life approach.

Fig. 3. PM threshold of Y.

5. Conclusions

Table 8 Distributional parameters and means of time to failure and PM thresholds.

F(t) Ff(t) Fp(t)

β

η

μ

1-m/mE

3.3198 3.2134 3.3104

1109.7 1172.9 879.1

995.8 1050.7 788.7

0.0523 0 0.2494

1

CDFs

0.8

Fp(t)

0.6

Ff(t)

0.4 0.2 0

0

500

1000

1500

2000

t Fig. 4. Plots of Ff(t) and Fp(t).

4.6. Applications of the CBM model 4.6.1. Population life distribution without condition information We can directly fit the data of time to failure, (t j ; 1 ≤j ≤17), to an appropriate model F(t) without considering the values of yðt j Þ. Such a life distribution is sometimes termed as the marginal distribution of failure time [4]. Considering the Weibull, gamma, normal, and lognormal distributions as candidate models, it is found that the Weibull distribution provides the best fit to the data in terms of log maximum likelihood value, and the estimated parameters are shown in the 2nd row of Table 8. On the other hand, the marginal distribution can be obtained from (25) with t 0 ¼ 0 and y0 ¼ 0. This needs to estimate the parameters of the gamma process model. Using the approach outlined in Section 3.1, we obtained the following parameter estimates:

In this paper, we have considered the CBM problem involving multiple condition variables. In this context, basic models are the composite condition variable model and its failure and PM thresholds. We proposed using the weighted power model to fuse the condition variables into a composite condition variable and using the Gaussian process model to represent the random failure threshold of the composite condition variable. We combined the item age and the composite condition variable into an equivalent age, and defined the PM threshold using the classic age replacement policy or the tradeoff BX life associated with the distribution of the equivalent age. We applied these models to derive the distributions of time to failure and PM thresholds. The approach has been illustrated by a real-world example. Main conclusions and findings have been the following: (a) The approach proposed in this paper allows censored observations. (b) The coefficient of variation is an appropriate measure for evaluating the failure prediction capability of a condition variable and identifying key condition indicators. (c) The composite condition variable generated from the weighted power model provides better failure prediction capability than individual condition variables. In the CBM context, a large power parameter p is desired, and hence the weighted power model significantly outperforms the linear model. (d) In the context of multiple condition variables, the failure threshold is usually random and time-dependent. The Gaussian process model is appropriate for representing such a failure threshold. (e) When the failure threshold is random, the item can fail at or before the PM threshold. The cost model and the tradeoff BX life approach can be used to specify the PM threshold. Topics for the future study include optimization of the condition monitoring scheme, and optimization of an alarm limit that provides the maintenance decision maker with a sufficient lead time (i.e., the time interval between the decision instant and the time to failure).

Acknowledgments λ ¼ 8:4386  10−3 ; v ¼ 9:8342  10−2

ð33Þ

A set of values of Ff(t) obtained from (25) is shown in Fig. 4 (i.e., those dots). Fitting these values to the Weibull distribution, we obtained the parameters shown in the 3rd row of Table 8. The last column of Table 8 shows that there is the relative error of 5.23% between the means obtained from F(t) and Ff(t). This confirms the appropriateness of the prediction method presented in Section 3. The fitted Weibull distribution is also shown in Fig. 4. As seen, it is in good agreement with the data points. This confirms the appropriateness of the approximation method.

The author would like to thank the two referees for their helpful comments and suggestions which have greatly enhanced the clarity of the paper. The research was supported by the National Natural Science Foundation (No. 71071026). References [1] Jardine AKS, Lin D, Banjevic D. A review on machinery diagnostics and prognostics implementing condition-based maintenance. Mechanical Systems and Signal Processing 2006;20(7):1483–510.

R. Jiang / Reliability Engineering and System Safety 119 (2013) 178–185

[2] Jiang R, Jardine AKS. Health state evaluation of an item: a general framework and graphical representation. Reliability Engineering and System Safety 2008;93(1):89–99. [3] Jiang R. Optimization of alarm threshold and sequential inspection scheme. Reliability Engineering and System Safety 2010;95(3):208–15. [4] Whitmore GA, Crowder MJ, Lawless JF. Failure inference from a marker process based on a bivariate Wiener model. Lifetime Data Analysis 1998;4:229–51. [5] Jiang R, Jardine AKS. Composite scale modeling in the presence of censored data. Reliability Engineering and System Safety 2006;91(7):756–64. [6] Wang P, Coit D. Reliability assessment based on degradation modeling with random or uncertain failure threshold. In: Proceedings of the 2007 annual reliability and maintainability symposium. p. 392–7. [7] Usynin A, Hines JW, Urmanov A. Uncertain failure thresholds in cumulative damage models. In: Proceedings of the 2008 annual reliability and maintainability symposium. p. 334–40. [8] Wang W, Carr M. An adapted Brownion motion model for plant residual life prediction. 2010 Prognostics and System Health Management Conference. Macau, China, MU3048.

185

[9] Wang W, Carr M, Xu W, Kobbacy KAH. A model for residual life prediction based on Brownian motion with an adaptive drift. Microelectronics Reliability 2011;51:285–93. [10] Jiang R. A tradeoff BX life and its applications. Reliability Engineering and System Safety 2013;113:1–6. [11] Jiang R. A general proportional model and modelling procedure. Quality and Reliability Engineering International 2012;28(6):634–47. [12] Jiang R, Murthy DNP. A study of Weibull shape parameter: properties and significance. Reliability Engineering and System Safety 2011;96(12):1619–26. [13] Wang W. A simulation-based multivariate Bayesian control chart for real time condition-based maintenance of complex systems. European Journal of Operational Research 2012;218:726–34. [14] Bullen PS. Handbook of means and their inequalities. Dordrecht: Kluwer Academic Publishers; 2003. [15] Barlow RE, Proschan F. Mathematical theory of reliability. New York: John Wiley & Sons; 1965. [16] van Noortwijk JM. A survey of the application of gamma processes in maintenance. Reliability Engineering and System Safety 2009;94:2–21.