Reliability Modeling for Systems with Multiple Degradation Processes Using Inverse Gaussian Process and Copulas

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 twostage 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.


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 sindependence.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.

Mathematical Problems in Engineering
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.

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 {();  ≥ 0}.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.
Note that, conditional on  −1  , the unconditional PDF of x  can be obtained by integrating  −1  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.

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  1 ,  2 , . . .,   , then there exists a unique copula  such that, for all  1 ,  2 , . . .,   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 V 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 Table 1: Summary of some one-parameter Archimedean copulas.

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.
Using MLE, the full log-likelihood function of the copula with marginal function can be expressed as ln  (, ) 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.
Based on (9), the log-likelihood function of each degradation process can be expressed as ln   (  ;   ) ,  = 1, 2, . . ., , (15) 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   = [ 1 ,  2 , . . .,   ] as the unobserved data and x  as the observed data.Given the complete data including   and x  , the log-likelihood of unknown parameters   = (  ,   ,   ,   ), up to a constant, can be expressed as For the E-Step, we need to estimate ((  ,   ,   ,   )) given the observed data and current estimated parameters.Ye and Chen [6] where It is noted that  1 (  ,   ) only depends on   ,   , and  2 (  ,   ) only depends on   ,   .So, in the M-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 M-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 loglikelihood function.where Θ is the parameter space.( 1 (⋅;  1 ),  2 (⋅;  2 ), . ..,   (⋅;   ); ) is the density function derived from copula (⋅; ).
The copula parameter  can be estimated by maximizing the full log-likelihood function in (22).

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 ln() is the maximized value of the loglikelihood 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 V 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].

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  = 10,  = 0.5,  = 1, and  = 1.5.We choose the number of test units to be  = 10, 15, and 30 and the number of inspections times to be  = 5, 10, 15, and 20.Each unit is inspected at   = , for  = 1, 2, . . ., .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.

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 sdependent 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.

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.

Figure 2 :Figure 3 :
Figure 2: Estimated mean paths based on sample average and on our model.

Figure 4 :
Figure 4: Comparison of the system reliability under both sdependent and s-independent assumptions.
By using the property of truncated normal distribution, we can easily obtain   = [ −1  | x  ] and   = [ −2  | x  ] for each unit.After substituting them into (16), we have  (ln  (  ,   ,   ,

Table 2 :
Biases and standard errors of the estimates based on the proposed EM approach.

Table 4 :
Estimation and the 90% confidence intervals.
Figure 1: The development of crack length over time: (a) crack A, (b) crack B.

Table 5 :
Goodness-of-Fit for the four copulas.