Abstract

We develop a reliability model for systems with s-dependent degradation processes using copulas. The proposed model accommodates assumptions of s-dependence among degradation processes and allows for different marginal distributions. This flexibility makes the model more attractive compared with the multivariate distribution model, which lay on the limitation of the homogeneous marginal distribution and can only describe linear correlation. Marginal degradation process is modeled by the inverse Gaussian (IG) process with time scale transformation. Furthermore, we incorporate random drift to account for the possible heterogeneity in population. This paper also develops the statistical inference method using EM algorithm with two-stage procedure. The comparison results of the reliability estimation under both s-dependent and s-independent assumptions are illustrated in the illustrative example to demonstrate the applicability of the proposed method.

1. Introduction

For systems with high reliability, it is challenging to do reliability assessment in a timely manner because usually minor failures occur during such short time. Thus it is difficult to assess reliability based only on limited lifetime data. Compared with traditional lifetime data, degradation data contains more information, which records the accumulation of damage over time. Degradation analysis aims to characterize the underlying failure process and uses more information than lifetime data analysis. In reality, a system may have multiple degradation processes or consists of multiple degrading components. The reliability of a system is usually calculated under the assumption of s-independence. However, assuming s-independence among degradation processes may underestimate system reliability. It is more realistic to assume some sort of dependence among degradation processes.

In practice, the continuous degradation process is often modeled by a stochastic process. Wiener process and the Gamma process are two classes of stochastic processes that have been widely used in degradation modeling [1]. A distinct feature of the Wiener process is that its degradation measures are not necessarily monotone [2], which is not applicable in many cases. As an alternative, the Gamma process is often used when monotonicity is required [3]. Although the Wiener process and Gamma process have wide applications in degradation modeling, only these two classes of models cannot fit all degradation data well. Wang and Xu [4] showed an application in which both Wiener and Gamma processes do not fit the data well. They try to fit the GaAs Laser data by using Wiener and Gamma process models, but probability plots show that both of them fit the data poorly. They suggest inverse Gaussian (IG) process as another good choice for degradation modeling which provides a monotone degradation path. The IG process was proposed by Wasan [5] in 1968. Notwithstanding wide applications of the IG distribution, the IG process was scarcely used in degradation modeling. Ye and Chen [6] attributed the reason of its unclear physical meanings to reliability engineers and prove that IG process is the limit of a compound Poisson process whose arrival rate goes to infinity while the jump size goes to zero in a certain way. Their research provides a physical meaning for the IG process as a degradation model.

For multiple degradation processes, Wang and Coit [7] introduced a general modeling approach to estimating the reliability of a system using multivariate normal distribution under the condition that degradation mechanism is unknown. Pan and Balakrishnan [8] introduced a reliability model for systems with bivariate degradation processes by utilizing bivariate Birnbaum-Saunders distribution. In their model, the degradation process is described by Gamma process and distribution of the first passage time is approximated by Birnbaum-Saunders distribution. In recent years, much attention has been placed on modeling the s-dependence behavior between degradation variables by copulas, which allows us to couple degradation variables from different distribution families. Sari et al. [9] introduced a two-stage model for bivariate processes using the copula function. Experiment data from LED lighting systems are analyzed to illustrate the application of the proposed model. Pan et al. [10] proposed a bivariate Wiener degradation process to describe the dependence between degradation characteristics, where marginal distribution and multivariate dependence can be modeled separately, and marginal distribution functions of degradation measures may follow different distributions families.

Based on the bivariate Wiener degradation model, Wang et al. [11] put forward an adaptive method for residual life estimation. In their method, the dependence of degradation variables is characterized by the Frank copula function. Hao and Su [12] presented a bivariate nonlinear Wiener degradation model and use Markov Chain Monte Carlo (MCMC) method to estimate the model parameters. Wang and Pham [13] established a s-dependent competing risk model to handle the dependence among multiple degradation processes and random shocks using copulas. Two types of dependence are considered in their model: dependence among degradation processes and dependence between degradation and shocks.

Previous research has been focused on utilizing multivariate distributions. This approach has some drawbacks. One major drawback is that it limits the distribution of each degradation variable. Degradation variables should have the same distribution. Meanwhile, such approaches are applicable only when degradation variables are linearly correlated. Some researchers use copulas to model the dependence among degradation processes, which allows for different marginal distributions. A more interesting problem that has not been addressed in the literature is the effect of copulas on s-dependence modeling for systems with IG processes. This paper makes a further study on s-dependence modeling using IG process and copulas. We first use the IG process with random drift and time scale transformation to model the monotonic degradation process; then we employ the copula method to fit the joint distribution of multiple degradation processes. We also develop the statistical inference method using EM algorithm with two-stage procedure to estimate model parameters.

The remainder of the paper is organized as follows. Section 2 introduces an IG process with random drift and time scale transformation. Section 3 proposes a system reliability model under s-dependent assumption using copulas. The statistical inference for the model parameters is discussed in Section 4. In Section 5, a simulation study is conducted to demonstrate the quality of the estimators. Section 6 gives an example of crack growth problem to demonstrate the proposed method. Section 7 concludes the paper.

2. IG Process with Random Drift and Time Scale Transformation

In practice, many systems or components have their performance characteristics whose degradation over time can be measured. The degradation path over time can be modeled by a stochastic process . In our study, we use the inverse Gaussian (IG) process to describe monotonic degradation process.

Definition 1. The inverse Gaussian process is defined as the stochastic process satisfying the following.(i) follows an inverse Gaussian distribution with mean and scale .(ii) has independent increments on nonoverlapping intervals; that is, and are independent, .

Consider that the IG process with follows , where is the transformed time scale; the transformation depends on particular mechanism. General discussions about time transformation can be found in Whitmore and Schenkelberg [15]. is the drift parameter; is the diffusion parameter. Let and ; thus follows with mean and variance , where denotes the IG distribution with probability density function (PDF): and cumulative distribution function (CDF)

In many engineering applications, the failure of a system is defined when the degradation process crosses a critical threshold level for the first time. The corresponding time is

Because IG process has monotonic property, the reliability function of the first passage time can be obtained from the relation as

Suppose that a system has degradation processes and each of them is governed by an IG process. We assume that all samples have common inspection times. The inference can be easily extended to the case where inspection times are different. In a degradation test, units are tested and the degradation measures of all units are observed at ordered inspection times with observations , , , , where is the total number of degradation processes, is the total number of test units, and is the total number of inspection times. In general, the degradation data can be presented in the form

Because of the independent property of IG process, increments for nonoverlapping intervals are independent. Therefore, we focus on the degradation increment. Define , , , and for , , . We have where and are parameters for th degradation process.

The joint PDF of is given by

Considering heterogeneity within the population, random effects are often used to account for heterogeneous degradation rates across units [16, 17]. Normality assumption for has been adopted by many studies [18]. To avoid the negative values of and simplify the deviation, we assume that conditional on follows a truncated normal distribution . This means that is fixed but unknown for each sample. Its probability density function (PDF) for is given by where is the PDF of standard normal distribution, is the CDF of standard normal distribution, and . Note that, conditional on , the unconditional PDF of can be obtained by integrating out of (7), which yields

The random drift model covers simple model without random effects when . The time scale transformation can also be eliminated if we let . When both and , the random drift model with time scale transformation degenerates to the regular IG process. This means that the random drift model with time scale transformation can be viewed as an extended model which includes several IG process models as special cases. In degradation data analysis, an extended model is a much flexible choice to fit a specific dataset.

3. Modeling for Systems with Multiple Degradation Processes

Although, in many studies, multiple degradation processes are assumed to have independent property, it may be more realistic to assume some sort of dependence among degradation processes. According to Sklar’s theorem [19], if is a joint distribution function with continuous margins , then there exists a unique copula such that, for all in ,

This theorem provides the theoretical foundation for the application of copulas. Based on Sklar’s theorem, any multivariate distribution can be decomposed into a copula and its marginal. Thus, copula functions provide a much more flexible and realistic way to study multivariate distributions. Table 1 lists some well-known one-parameter Archimedean copulas, in which and are marginal distributions.

Properties of the copula have been explored in a number of studies; for example, see Jia et al. [20].

Consider a system subject to multiple s-dependent degradation processes. The degradation processes are denoted as , where is the total number of degradation processes and are degradation process for . The joint distribution of is where is the CDF of for .

The system is considered to be failed if either degradation process reaches its corresponding failure threshold, which is known as for . Therefore, the reliability of a system with s-dependent degradation processes is where means safe region, within which the system is working and is the CDF of for .

4. Statistical Inference

4.1. Time Scale Determination

In a few situations we understand well the parametric form for time scale based on physics mechanism or prior knowledge. When there is no such failure mechanism to describe the time scale, we use a statistical method to help us determine the parametric form. From (6), we can see that the mean degradation path is , where can be obtained by using the property of truncated normal distribution. We may average the samples to obtain the mean degradation curve and specify a parametric form for by fitting the mean path. Linear relationship, power law relationship, and exponential law relationship are the most used parametric forms.

4.2. Parameter Estimation

Suppose that the multivariate distribution is specified by margins with CDF and PDF , , and a copula . Let be the parameter in copula and the parameters in the th degradation process. We will remark that is the parameter in time scale . Now, consider the joint distribution of degradation increments in the interval . The joint PDF can be computed by where and are the th marginal CDF and PDF with parameter ; , is the density of copula ; and is degradation increment in the interval .

Using MLE, the full log-likelihood function of the copula with marginal function can be expressed as

The parameters can be estimated by maximizing the above log-likelihood function. However, the log-likelihood is computationally difficult to work with, or even infeasible. Here we use a two-stage procedure of firstly estimating the marginal parameters and then estimating the copula parameter from the joint likelihood with the marginal parameters estimated in the first stage. Computational difficulty is greatly reduced since each stage has a very small number of parameters. Joe [21] has studied the efficiency properties of this two-stage estimation procedure. He points out that this method has good efficiency and can sometimes be equivalent to MLE method.

Stage 1 (marginal parameters). Each marginal distribution has its own parameter set . In the first stage, the log-likelihood of the marginals is separately maximized to get estimates .

Based on (9), the log-likelihood function of each degradation process can be expressed as where is the total number of degradation processes, is the total number of test units, and is the total number of inspection times.

The log-likelihood function is very flat around . Direct maximization of the likelihood function often yields a solution far away from the MLE. We use Expectation-Maximization (EM) algorithm to solve this problem as it is free of the rounding error. We denote as the unobserved data and as the observed data. Given the complete data including and , the log-likelihood of unknown parameters , up to a constant, can be expressed as

For the -Step, we need to estimate given the observed data and current estimated parameters. Ye and Chen [6] have proven that the conditional distribution of follows . By using the property of truncated normal distribution, we can easily obtain and for each unit. After substituting them into (16), we have where

It is noted that only depends on , , and only depends on , . So, in the -Step, we can compute the updated parameter estimates separately. The estimation equations for , can be easily obtained by maximizing (18); that is,

We can also obtain , by solving

After obtaining updated parameter estimates in the -Step, we can iterate to obtain new estimates until convergence. It is worth noting that both (20) and (21) contain only two parameters. It is much easier to find the optimal solution, as opposed to the direct maximization of the original log-likelihood function.

Stage 2 (copula parameters). With given in Stage 1, the estimation of the parameter from the copula function can be performed as where is the parameter space. is the density function derived from copula .

The copula parameter can be estimated by maximizing the full log-likelihood function in (22).

4.3. Goodness-of-Fit and Model Selection

When there are several potential copula functions available, the most straightforward way to quantitatively select the model is to use the Akaike Information Criterion (AIC). The AIC is defined as where is the number of unknown parameters in the statistical model and is the maximized value of the log-likelihood function for the estimated model.

Based on this criterion, the model with the smallest AIC value is selected as the best fitting model. However, the model with the smallest AIC does not necessarily be a correct model. We should also check the correlation coefficient between variables. The most widely known correlation coefficients are Pearson’s , Kendall’s tau, and Spearman’s rho. Because Pearson’s only measures the linear correlation between random variables, it is not a proper measure of association in this study. We only focus on Kendall’s tau and Spearman’s rho, which does not rely on linear correlated restriction. When applying these correlation coefficients to the bivariate degradation processes, Kendall’s tau can be stated as and Spearman’s rho can be stated as where and are the CDFs of marginal random variables.

The model that fits the data well is expected to yield the maximum Kendall’s tau and Spearman’s rho. More details could be found in Nelsen [19].

5. Simulation Study

To demonstrate the quality of the EM estimators proposed in last section, a simulation study has been carried out. Assume as time scale. Without loss of generality, we fix , , , and . We choose the number of test units to be , 15, and 30 and the number of inspections times to be , 10, 15, and 20. Each unit is inspected at , for . Under each parameter setting (, ), the bias and the mean square error (MSE) of the MLE are computed based on 1000 Monte Carlo replications. The results are displayed in Table 2.

From Table 2, it is readily observed that the biases and MSEs become smaller as the number of test units and inspections times increase. The result reveals that the parameters of our model can be accurately estimated with reasonable sample sizes.

6. Illustrative Example

Lu and Meeker [14] presented a fatigue-crack-growth dataset. Degradation data from 21 samples are collected. The time unit is in million cycles; all samples are measured every 0.01 million cycles from 0 to 0.09. We choose the first 20 samples from the dataset and divide them into two groups. We assume that samples in the first group represent crack A, while samples in the second group represent crack B. The length growth of the cracks can be seen as degradation processes. Considering a system subject to fatigue cracks in two different positions, the system is defined to be failed if the length of crack A crosses 1.6 inch or the length of crack B crosses 1.3 inch. The fatigue-crack-growth data are listed in Table 3 and depicted in Figure 1. A similar idea is also adopted in other study cases to analyze their models [8]. Now, we further apply our method to the dataset.

We first apply the IG process with random drift and time scale transformation to fit the fatigue-crack-growth data. Figure 1 shows that degradation path is not linear over time. Lu and Meeker [14] suggested that Power Law is appropriate for crack size modeling, and thus is adopted as the transformed time scale. The estimates of marginal parameters are listed in Table 4. The 90% confidence intervals are estimated using the parametric percentile bootstrap method with 2000 bootstrap replication, as shown in Table 4.

Based on (6), the estimated crack length paths is given by

To demonstrate the goodness of fit, the mean crack length paths estimated from the above model and from the sample average are compared. If these two lines overlap, it suggests that the model is correct. The mean crack length paths based on (26) are plotted in conjunction with the sample average, as shown in Figure 2. In addition, the result by bivariate Birnbaum-Saunders method from Pan and Balakrishnan [8] is also depicted in the same figure for comparison. Intuitively, it can be seen that our model tally quite well with sample average, indicating that the model fits the data well. From Figure 2 we can also see that crack A has a higher degradation rate than crack B.

The marginal reliability of these two cracks obtained from the estimated IG distribution is displayed in Figure 3. Although the degradation rate of crack A is higher than crack B, the reliability curve of crack A tends to be higher than crack B. The reason is that crack B is defined to be failed if crack size exceeds 1.3 inch, which is lower than the threshold of crack A.

Given the estimates of marginal parameters in Table 4, the copula method is applied to fit the dependent crack size data. To choose the most suitable fitting copula, we test six selected copulas using the AIC as the criteria for goodness of fit. The results are summarized in Table 5.

From Table 5, we see that the model with Frank copula has the smallest AIC of −1118.13, followed by the Gumbel copula with AIC of −1114.89. Clayton copula is the worst fitting among all the copulas, with AIC of −1095.96. That is, the AIC criterion favors the Frank copula. The model with Frank copula also has both the largest Kendall’s tau and Spearman’s rho. Therefore, Frank copula may be the best choice to describe the dependence between the two degradation processes.

Given the estimated parameters, we can derive the system reliability from (12), as shown in Figure 4. By means of the bivariate Birnbaum-Saunders method proposed by Pan and Balakrishnan [8], their estimated reliability curve is also depicted in Figure 4, which is very similar to the s-dependent reliability curve from our method. By comparing these two curves we can further validate the effectiveness of our method. As we can see from this figure, s-dependent reliability curve and s-independent reliability curve are overlapped at early time. Because at the early stage both of the two degradation processes are in a good status, and unlikely to impact on each other. However, the reliability estimated by Pan and Balakrishnan [8] is a little bit high at early time. This evidence tends to suggest that the bivariate Birnbaum-Saunders method may not be the best choice for the data. After a period of time, s-independent reliability curve tends to be slightly lower than s-dependent reliability curve because degradation processes are subject to the same stress condition; small degradation amount of one degradation process tends to occur with small degradation amount of the other degradation process. That is, there exists some sort of dependence between the degradation processes.

7. Conclusions

This study has investigated reliability modeling method for systems with s-dependent degradation processes. We first introduced an inverse Gaussian (IG) process with random drift and time scale transformation to describe marginal degradation process. This is an extended family of IG process models that contains several existing models as its limiting cases. Based on copulas, we then developed a system reliability model subject to s-dependent degradation processes. The statistical inference method was developed using EM algorithm with two-stage procedure. In the crack length growth problem, we applied the model to fit the crack size data. The estimation accuracy was compared in a numerical example. The results revealed that our model is suitable when the degradation path is nonlinear. The system reliability tends to be higher than s-independent situation after early stage. This is expected because our model considers the dependence among degradation processes. If we ignore the dependent relationship, significant information loss is expected. Therefore, our model is potentially very useful for multivariate degradation data analysis, because a common way to model for multivariate data has limitation of the same marginal distribution and can only describe linear correlation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the National Natural Science Foundation of China under Agreements 61104133 and 11001005.