A General Accelerated Degradation Model Based on the Wiener Process.

Accelerated degradation testing (ADT) is an efficient tool to conduct material service reliability and safety evaluations by analyzing performance degradation data. Traditional stochastic process models are mainly for linear or linearization degradation paths. However, those methods are not applicable for the situations where the degradation processes cannot be linearized. Hence, in this paper, a general ADT model based on the Wiener process is proposed to solve the problem for accelerated degradation data analysis. The general model can consider the unit-to-unit variation and temporal variation of the degradation process, and is suitable for both linear and nonlinear ADT analyses with single or multiple acceleration variables. The statistical inference is given to estimate the unknown parameters in both constant stress and step stress ADT. The simulation example and two real applications demonstrate that the proposed method can yield reliable lifetime evaluation results compared with the existing linear and time-scale transformation Wiener processes in both linear and nonlinear ADT analyses.

1. Introduction

Due to the highly competitive market, nowadays many products are requested to have long lifespans and high reliability. In order to quantify their lifetime and reliability characteristics, accelerated degradation testing (ADT) is conducted under harsh environmental conditions to accelerate the performance degradation and obtain sufficient data for reliability analysis in a short time with a limited budget [1,2,3]. Thus, it has been widely used in many engineering applications, e.g., batteries [4,5], light emitting diodes (LEDs) [6,7], etc.

The aim of statistical analysis of ADT data is to extrapolate product characteristics that are of interest at normal stress levels which, in general, comes from two directions. In the time direction, the deterministic degradation trend varying with time should be modeled at all accelerated levels for linear or nonlinear scenarios, i.e., the degradation model. In the stress direction, the relationship between acceleration variables and degradation-related parameters need to be established based on physical mechanisms or empirical experience, i.e., an acceleration model which will be used for the lifetime and reliability extrapolation at the normal stress level. In addition, the uncertainties from the temporal variation of the degradation process and unit-to-unit variation from the inherent heterogeneity of the tested products should be properly considered when analyzing ADT data. In the literature, two classes of models are widely used for degradation modeling, i.e., degradation path models [8,9] and stochastic process models [10,11,12], with examples for acceleration models including the Arrhenius model for temperature stress and the Eyring model for voltage stress [3,8,13]. With different combinations of degradation and acceleration models from these two directions, extensive work has been given to the statistical analysis of ADT data.

Since the stochastic process models can describe the temporal variation of the degradation process in a finite time interval, more attention has been given to them. In linear ADT analysis with a single acceleration variable, Park and Padgett [14] proposed two accelerated degradation models based on geometric Brownian motion and the Gamma process, and analyzed the constant stress ADT (CSADT) data of carbon-film resistors under temperature stress. Pan and Balakrishnan [15] proposed two linear modeling methods based on the Wiener and Gamma processes, and simulated one step stress ADT (SSADT) data with the Arrhenius model for model verification. Wang, et al. [16] used a linear Wiener process to model the degradation process for reliability evaluation of products with the integrated ADT and field information. In the cases with multiple acceleration variables, Liao and Elsayed [17] selected Brownian motion with a linear drift as the degradation model to infer the field reliability and applied it to LED CSADT data with electric current and ambient temperature.

In reality, the degradation processes of some products may experience nonlinear due to the inner deterioration mechanism of the product material, e.g., crack growth. In order to analyze this kind of ADT data, a time-scale transformation is generally used based on linear stochastic process models with a single acceleration variable. To our knowledge, Whitmore and Schenkelberg [11] were the first to accomplish this work in the area of ADT and proposed a time-scale transformation Wiener process model for the cable CSADT data analysis under temperature stress. A similar transformation is introduced in the Gamma process [10] and inverse Gaussian process models [12] for ADT modeling. In consideration of the unit-to-unit variation, Tang, et al. [18] incorporated a random variable into the acceleration model to capture the random effect and also used the time-scale transformation Wiener process model for nonlinear LED CSADT data analysis under electric current stress. However, the problem of using time-scale transformation is the implied assumption that the degradation processes can be linearized. It may not be suitable for nonlinear ADT analysis where the degradation paths cannot be linearized. In traditional degradation analysis, Wang, et al. [19,20] proposed a general Wiener model which utilizes two time-scale parameters to extend the above transformation and can solve the nonlinear degradation modeling to a greater extent. In our previous work [21], we introduced this model into ADT analysis, but without the consideration of unit-to-unit variation.

Based on the above research, it can be concluded that the existing methods partially solve the modeling problems of some ADT data with their applications in both the time and stress directions. However, more efforts are still needed, especially for nonlinear ADT modeling with multiple acceleration variables. Therefore, in this paper, a general ADT model is proposed to fill this gap based on the Wiener process. The general ADT model will consider the unit-to-unit variation and temporal variation of the degradation processes, and is suitable for both linear and nonlinear ADT analyses with single or multiple acceleration variables, simultaneously.

The rest of this paper is organized as follows: Section 2 introduces the general ADT model and derives its failure time distribution under the given stress level; Section 3 deals with the statistical inference of unknown parameters in the proposed model for both CSADT and SSADT; Section 4 uses one simulation example and two real applications to illustrate the superiority of the proposed method over other existing Wiener processes for both linear and nonlinear ADT analyses; and some concluding remarks are given in Section 5.

2. The General ADT Model

2.1. Models

The general ADT model based on the Wiener process is given by: where X(t) is the performance degradation value of product at time t, σ denotes the diffusion coefficient, which is assumed to be constant, B(.) is the standard Brownian motion to describe the temporal variation of the degradation process, Λ(.) and τ(.) are the monotonous functions with time, θ and γ are the time-scale parameters to present the linear or nonlinear modeling. Without loss of generality, the initial degradation value is set to be zero. If not, X(t) = X(t) − X(0) will be used.

For the drift coefficient μ(;), it is assumed to be dependent with the acceleration variable as the acceleration model. The acceleration model for both single and multiple acceleration variables is denoted as:

If there are p acceleration variables, then = [s1, s2, …, s] and = [η0, η1, …, η], where s is the vth acceleration stress type and η is the vth constant coefficient. While φ(.) is the continuous function of acceleration variables s, i = 1, 2, …, p. For instance, if φ(s) = exp(1/s), Equation (2) presents the Arrhenius relationship (exponential type); while if φ(s) = s, Equation (2) presents the Eyring relationship (power rule type).

Considering the unit-to-unit variation due to the inherent heterogeneity among the tested products, the drift coefficients μ(;) are varied from product to product. Hence, we consider it as a random variable to present this kind of variation. Similar methods can be found in Peng and Tseng [22], Wang [23], and Si, et al. [24]. Here, for simplicity, the coefficient η0 in Equation (2) is assumed to follow a normal distribution with mean value a and variance b, i.e., η0 ~ N(a, b). The parameter values should be such that Pr(μ(;) < 0) is nearly zero to avoid negative values existing in the drift coefficient of Equation (1). Thus, if such a random effect is not considered (b = 0), Equation (2) will become the traditional acceleration model in Park and Padgett [25].

The attraction of Equation (1) is that it can cover two commonly used Wiener processes in the ADT field as its limiting cases, which are:

Case 1. If Λ(t;θ) = t and τ(t;γ) = t, Equation (1) is simplified to the traditional linear Wiener process, which is widely used for ADT analysis, see [16,25,26,27] and so on:

Case 2. If Λ(t;θ) = τ(t;γ), Equation (1) reduces to a time-scale transformation Wiener process as in [11,18]:

Herein, if the time-scale transformation z = τ(t;γ) is used in Equation (4), then Case 2 will become a similar form of Case 1 as Y(z) = μ(;)z + σB(z).

Furthermore, if the acceleration variables in Equation (2) are set to be at normal stress levels where μ(;) ~ N(μ0, ), Equation (1) and its limiting cases (i.e., Equations (3) and (4)) are the degradation models used in traditional degradation analysis [19,20,28].

As described above, the general ADT model in Equation (1) can present the uncertainties from the unit-to-unit variation (b ≠ 0) and temporal variation of the degradation process, and can be used in the situations with single (p = 1) and multiple (p > 1) acceleration variables for linear and nonlinear degradation processes. For clarity, Equation (1) is named M0, while Equations (3) and (4) are M1 and M2. Since model M1 and M2 are widely used in the ADT field, in this paper, we concentrate on the comparison of M1 and M2 with M0 to verify the effectiveness of the proposed model for both linear and nonlinear ADT analyses.

For the evaluation of lifetime and reliability, the probability distribution function (PDF) of the failure time should be given, which will be derived in the following section.

2.2. Derivation of the Failure Time Distribution under the Given Stress Level

Although ADT is implemented at accelerated stress levels, the lifetime and reliability evaluation are conducted at the given normal stress level 0 = [, ,⋯, ], where is the vth acceleration stress type under normal level. Thus, the drift coefficient μ(0;) follows a normal distribution. Therefore, we simplify the notation of μ(0;) to μ in the derivation of the failure time distribution. From Equation (2), it is known that μ ~ N(μ0, ), where and .

In general, the failure time is defined as the time when degradation path X(t) first exceeds the failure threshold ω, i.e., the first passage time (FPT):

Let a time transformation z = τ(t;γ); thus t = τ−1(z;γ). We define ρ(z;θ) = Λ(τ−1(z;γ);θ). So, Equation (1) becomes [20]: and its drift coefficient is:

For simplicity, let h(z;θ) = dρ(z;θ)/dz. Under some mild assumptions, the PDF of the FPT for the new degradation process Y(z) is (see Theorem 2 in Si, et al. [24]):

In consideration of the random effect due to unit-to-unit variation, the uncertainty of μ should be included in the PDF of the FPT. We compute the result by the law of total probability, i.e.,:

In order to explicitly obtain the formula of Equation (9), Theorem 3 of Si, et al. [24] is introduced and modified accordingly; that is:

If μ ~ N(μ∈ :

Therefore, we substitute Equations (8) into (9) according to Equation (10). The result is: where Q = ρ2(z;θ) + σ2z.

Hence, the PDF of the FPT for the general model can be obtained by the reverse process of the time transformation z = τ(t;γ) through Equation (11); that is: where Q = Λ2(t;θ)+ σ2τ(t;γ); G = Λ(t;θ) − h(τ(t;γ);θ)τ(t;γ).

Supposing that Λ(t;θ) = t and τ(t;γ) = t, Equation (12) becomes:

Herein, the relationship ∫p(t)dt = 1 should be satisfied. Therefore, the PDF and cumulative distribution function (CDF) of the FPT for model M0 are modified as:

For M1, i.e., Λ(t;θ) = t and τ(t;γ) = t, the PDF of the FPT for M1 is known as an inverse Gaussian distribution [29]. Considering the random effect, the PDF and CDF of the FPT are [22,28]:

Obviously, Equation (15) is a limiting case of Equation (14) if we substitute θ = γ = 1 into Equation (13), where ∫p(t)dt = 1.

For M2, i.e., Λ(t;θ) = τ(t;γ), the PDF of the FPT for M2 is in accordance with Equation (15) through a time-scale transformation by replacing t into Λ(t;θ) or τ(t;γ) [18,30], which is also a limiting case of Equation (14).

In the following section, the problems of parameter estimation will be addressed for CSADT and SSADT. After that, the PDF and CDF of the FPT at normal stress levels will be computed through Equation (14) for the lifetime and reliability evaluation of long lifespan products.

3. Statistical Inference

With different loading profiles of accelerated stresses, CSADT and SSADT are widely used in engineering applications. Figure 1a,b show the schematic of these two types of ADT under three stress levels. In the left, samples are separated into different groups and tested under different constant stress levels, i.e., 1, 2 and 3, while, on the right, all samples are in the same group and tested from the lowest stress level to the highest in a step-by-step manner, i.e., 1 → 2 → 3. Thus, the advantage of SSADT is that it can save the number of samples with a shorter test time [31].

Figure 1

The schematic of (a) CSADT and (b) SSADT under three stress levels.

Herein, we provide the method for estimating the unknown parameters in CSADT. Given in Section 2, the unknown parameter vector is Ω = [θ, γ, σ, a, b, η], v = 1, 2, …, p. The analytic expressions of those parameters are hard to obtain directly. Hence, a two-stage maximum likelihood estimation (MLE) method is proposed to address this issue. In the first stage, the parameters related to the degradation process model are estimated, i.e., Ω1 = [θ, γ, σ] in Equation (1). In the second stage, the rest of the parameters related to the acceleration model, i.e., Ω2 = [a, b, η] in Equation (2), are given accordingly.

3.1. Estimation of Ω1 for CSADT

The observed CSADT data x is the kth degradation value of unit j under the ith stress level and t is the corresponding measurement time, i = 1, 2, …, K, j = 1, 2, …, n, k = 1, 2, …, m, where K is the number of stress levels, n is the number of test samples under the ith stress level, and m is the number of measurements for unit j under the ith stress level. Let = (x1, x2, …, x)′ and = (Λ(t1;θ), Λ(t2;θ), …, Λ(t;θ))′. According to the properties of the Wiener process, the degradation value follows a multivariate normal distribution: where μ is the drift coefficient of unit j under the ith stress level:

Let . The likelihood function of the CSADT data can be easily obtained by Equation (16) and the logarithm function is:

The first partial derivative of Equation (17) to and are:

The maximum value of the log-likelihood function in Equation (17) is obtained when Equations (18) and (19) are equal to zero. Thus, the estimates of and relying on θ and γ are:

Substituting Equations (20) and (21) into (17), the log-likelihood function is only a function of θ and γ, i.e.,:

Thus, and can be obtained by the two-dimensional search for the maximum value of Equation (22) [32]. Then, the estimation of and can be easily computed by substituting and into Equations (20) and (21).

3.2. Estimation of Ω2 for CSADT

The estimate of Ω2 is related to the acceleration model in Equation (2). From Section 3.1, the estimates of the drift coefficients for unit j under the ith stress level are given and the corresponding stresses are , i = 1, 2, …, K, j = 1, 2, …, n. With the consideration of unit-to-unit variation in Equation (2), the relationship among them is denoted as: where will be replaced by its estimate from Equation (20).

For simplicity, let . Thus, the log-likelihood function for Ω2 is:

Similar to the estimation procedure of Ω1, we compute the estimation of and relying upon , which are:

Then, substituting Equations (25) into (24), the can be obtained by a p dimensional search for the maximum value of Equation (24). After that, the estimates of and can be given by Equation (25), accordingly.

3.3. Estimation of Ω1 and Ω2 for SSADT

The estimates of unknown parameters in SSADT are slightly different from that in CSADT. If we divide SSADT data by the number of the stress levels K and set the initial values at each stress level to be zero, the SSADT data can be transformed into CSADT data. For instance, in Figure 1b, the degradation values from 0 to x1 during the time interval is obviously a subset of CSADT data. For the time intervals and , a transformation will be given since the initial values of both time and degradation value are not zero, as they are in CSADT.

Keeping this in mind, we assume that and its corresponding degradation value is . According to the properties of the Wiener process, the degradation value x satisfies that x − x1 ~ N(μ2(Λ(t;θ) − Λ(t1;θ)), σ2(τ(t;γ) − τ(t1;γ))) where μ2 is the drift coefficient when . Given that x′ = x − x1, Λ′(t;θ) = Λ(t;θ) − Λ(t1;θ) and τ′(t;γ) = τ(t;γ) − τ(t1;γ). The SSADT data in the time interval can be interpreted as a subset of CSADT data where the transferred degradation value is from 0 to x′ during the time interval , where is Λ′(t;θ) or τ′(t;γ) accordingly. This procedure is the same for the data in the time interval (t2, t3].

Following this idea, let the observed SSADT data be x, which is the kth degradation value of unit j at the stress level i and t is the corresponding measurement time, i = 1, 2, …, K, j = 1, 2, …, n, k = 1, 2, …, m. The transformations are given in Equations (26)–(28): where I× is the m × n all ones matrix; t0 = 0 and x0 = 0.

Then, the unknown parameters in SSADT can be given by the proposed two-stage MLE method for CSADT, as shown in Section 3.1 and Section 3.2.

4. Case Study

In this section, one simulation example and two real cases are used to demonstrate the validity and superiority of the proposed models and methods over other existing forms of the Wiener processes on linear and nonlinear ADT analyses.

The literature has shown that the expectation of the degradation path can be formulated by an exponential function [33]. Therefore, the forms Λ(t;θ) = t and τ(t;γ) = t will be used in this paper. In order to compare the s-fits of model M0 with M1 and M2, the Akaike information criterion (AIC) is introduced: where l is the maximum value of the log-likelihood function in Equation (17), and n is the number of unknown parameters. The lower the AIC value is, the better the model fits.

According to Burnham and Anderson [34], it is imperative to rescale AIC to: where AIC is the minimum of AIC values. Then, Δ = 0 for the best model, Δ ≤ 2 for models having substantial support, 4 ≤ Δ ≤ 7 for models having considerably less support and Δ > 10 for models having essentially no support compared to the best models.

In addition, for illustration purpose of the fitting results, a quantile-quantile (Q-Q) plot is used to graphically present the goodness of fit of each model on ADT data with the following standard normal distribution: where the plot is approximately linear if the fitting is satisfactory.

4.1. Simulation Example

SSADT data is simulated under three temperature stress levels, T = 60 °C, 100 °C, and 120 °C, and the normal temperature of 25 °C. The observation time interval is one hundred hours, while the number of observations under three stress levels are 15, 10 and 5, respectively. Hence, the total test time is 3000 h. In order to evaluate the influence of the sample size on parameter estimation and the reliability evaluation, we conduct three simulation analyses under the sample size n = 5, 10 and 30, respectively. The failure threshold ω is 100. Without a loss of generality, the parameters are set to be θ = 1.5, γ = 0.4, σ2 = 0.01, a = 20, b = 5, and η1 = −1500. The Arrhenius model is selected as the acceleration model, which is:

First, we check the normal assumption of μ0 at the normal stress level. Given by Equation (32) and the parameter settings, Pr(μ0 < 0) = 1.8720 × 10−19, which is approximately near zero. Thus, the non-negative assumption can be satisfied to compute the PDF of the FPT. In the following, both the model comparison and sensitivity analysis are conducted with the abovementioned parameter values.

4.1.1. Model Comparison

Given the ADT data with different sample sizes, the relative errors (RE) of the parameter estimation in percentage are computed by:

To investigate the variance and the bias of the estimator, the relative square error (RSE) of parameters in percentage is also given by Equation (34):

From the above settings, M0, M1 and M2 have 6, 4 and 5 parameters, respectively. Meanwhile, the absolute error (AE) of the candidate model M (i = 0, 1, 2) to the real model M is given by Equation (35) to quantitatively analyze the reliability evaluation results: where F(t) is the CDF of the FPT at the normal stress level given by Equation (13), t = 0.1, 1.1, …, 699.1, for hundreds of hours in this study, N = 700. Herein, if AE > 0, it means that model M overestimates the reliability evaluation results comparing with the true values, otherwise it underestimates.

Table 1 gives the estimates of unknown parameters with their REs and RSEs, and AEs of reliability estimation for different models at different sample sizes. Figure 2 illustrates the fitting results of each model on the simulation data when n = 10. Clearly, from l and Δ values at each sample size, M0 is the most suitable model, then is model M2, and the worst is model M1. The reason is that the time-scale transformation model M2 can capture the nonlinear property of the degradation process to some extent, but still perform worse than M0. As to M1, it tries to linearize the nonlinear degradation process, which leads to dreadful fitting results and poor parameter estimation compared with the true results.

Table 1

Simulation example: parameter estimates with REs and RSEs in percentage (in parentheses), and AEs of reliability estimation for three candidate models under three sample sizes.

Mi	n	θ	γ	σ2	η0		η1	lmax	np	AIC	Δ	AE
					a	b
M0	5	1.491	0.273	0.0177	21.62	4.33	–1502	317	6	–621	0	–2.9×10−3
		(−0.57, 3.3 × 10−3)	(−32, 10)	(77, 59)	(8.1, 0.65)	(−13, 1.8)	(0.13, 1.7 × 10−4)
	10	1.502	0.381	0.0097	18.65	4.42	–1477	635		–1259	0	–1.2 × 10−3
		(0.14, 2.0 × 10−4)	(−4.7, 0.22)	(−3.1, 0.096)	(−6.7, 0.45)	(−12, 1.3)	(−1.5, 0.023)
	30	1.501	0.352	0.0114	20.90	5.18	–1515	1898		–3784	0	2.7 × 10−4
		(0.055, 3.0 × 10−5)	(−12, 1.4)	(14, 2.1)	(4.5, 0.20)	(3.6, 0.13)	(1.0, 0.011)
M1	5	1 (fixed)	1 (fixed)	0.0607	8.5 × 103	7.9 × 105	–3045	−2.7	4	13.4	635	0.343
		(−33, 11)	(150, 225)	(507, 2.6 × 103)	(4.2 × 104, 1.8 × 107)	(1.6 × 107, 2.5 × 1012)	(103, 106)
	10	1 (fixed)	1 (fixed)	0.0579	8.2 × 103	8.3 × 105	–3049	1.7		4.6	1263	0.365
		(−33, 11)	(150, 225)	(479, 2.3 × 103)	(4.1 × 104, 1.7 × 107)	(1.7 × 107, 2.8 × 1012)	(103, 107)
	30	1 (fixed)	1(fixed)	0.0578	9.2 × 103	1.1 × 106	–3087	5.5		–3.0	3781	0.374
		(−33, 11)	(150, 225)	(478, 2.3 × 103)	(4.6 × 104, 2.1 × 107)	(2.2 × 107, 4.6 × 1012)	(106, 112)
M2	5	1.461	=θ	7.29 × 10−4	29.55	8.27	–1577	217	5	–425	197	3.2 × 10−3
		(−2.6, 0.068)	(265, 704)	(−93, 86)	(48, 23)	(65, 43)	(5.1, 0.26)
	10	1.476	=θ	4.79 × 10−4	24.54	7.64	–1544	490		–971	288	4.2 × 10−3
		(−1.6, 0.026)	(269, 724)	(−95, 91)	(23, 5.2)	(53, 28)	(2.9, 0.087)
	30	1.477	=θ	5.87 × 10−4	26.96	8.58	–1577	1380		–2750	1034	5.3 × 10−3
		(−1.6, 0.025)	(269, 724)	(−94, 89)	(35, 12)	(72, 51)	(5.2, 0.27)

Figure 2

The Q-Q plots for three candidate models for the simulated SSADT data when n = 10.

From the AE values in Table 1, it is known that both model M1 and M2 overestimate the reliability evaluation results at normal stress levels. With the sample size increasing, the AEs for M1 and M2 become larger since that they are not the right models for ADT analysis and the level of error will be amplified when more data are available for model validation. Meanwhile, model M0 slightly underestimates the reliability evaluation results from AE indexes, and the accuracies are improved by several orders of magnitude when the sample size goes from five to 30. The results demonstrate that model M0 is the most applicable model for nonlinear ADT analysis and can provide accurate reliability and lifetime evaluation results.

Figure 3 presents the PDFs and CDFs of the FPT for models M0 and M2 with the real values when the sample size n = 10. It is clear that model M0 is closer to the real values than M2. The mean time to failure (MTTF) for M0 and M2 are 8340 and 8720 h, while the true value is 8430 h. The results verify the effectiveness of model M0 than M2 with respect to nonlinear ADT analysis.

Figure 3

The (a) PDFs and (b) CDFs of the FPT for models M0 and M2 with the real values when n = 10.

4.1.2. Sensitivity Analysis

In this section, sensitivity analysis is conducted to analyze the robustness of the general model M0 with different values of model parameters [θ, γ, σ2, η0 (i.e., a, b), η1] for the simulation example.

Here, we set the parameters’ real values to 90%, 95%, 100%, 105% and 110% as the five factor levels. If model M0 is robust, its relative error of reliability evaluation results at normal stress levels should be as small as possible when compared with the real model M. Herein, we repeated the simulation procedure of SSADT data for N = 100 times and n = 10 samples will be generated at each time point. Then, the relative error for model M0 (RE of M0) is given through Equation (36): where is the CDF of the FPT at the normal stress level for the kth simulation.

Due to the extensive pairs of parameter combinations, say 56 = 15,625 and N = 100 repeated simulations, it will lead to heavy computational effort to compute the results. Thus, the orthogonal design of experiments is introduced to reduce the number of combinations, but still be able to find the sensitive parameters from the response (RE of M0) at each factor level [35]. The orthogonal array L25(56) is selected with 25 overall tests rather than 56.

The sensitivity analysis results are listed in Table 2, where the numbers from 1 to 5 refer to the first to fifth factor levels of the real values. For instance, in Test No. 1, the values are all 1, which means that 90% of all of the real values will be used to compute the corresponding response, which is 0.026812. The means of the responses for the five factor levels are listed at the bottom of Table 2, i.e., MRj, j = 1, 2, …, 5. The absolute biases are then calculated through δ = max(MRj) − min(MRj) which presents the influence of each parameter to the RE of M0. From Table 2, it is known that the sensitivity of the parameters is ranked as θ > b > η1 > γ > σ2 > a for the simulation study. Hence, special attention should be given to those parameters when using them for reliability and lifetime evaluation at normal stress levels. Furthermore, the absolute bias δ ranged from 0.01806 to 0.05734, indicating that model M0 is quite robust for lifetime and reliability evaluation.

Table 2

Sensitivity analysis of M0 with five levels of parameters through the orthogonal array L25(56) and Taguchi analysis.

Test No.	θ	γ	σ2	η0		η1	RE of M0
				a	b
1	1	1	1	1	1	1	0.026812
2	1	2	2	2	2	2	0.049430
3	1	3	3	3	3	3	0.076352
4	1	4	4	4	4	4	0.107904
5	1	5	5	5	5	5	0.149784
6	2	1	2	3	4	5	0.101002
7	2	2	3	4	5	1	0.017453
8	2	3	4	5	1	2	0.005818
9	2	4	5	1	2	3	0.042759
10	2	5	1	2	3	4	0.068059
11	3	1	3	5	2	4	0.014726
12	3	2	4	1	3	5	0.060681
13	3	3	5	2	4	1	0.032238
14	3	4	1	3	5	2	0.018415
15	3	5	2	4	1	3	0.006743
16	4	1	4	2	5	3	0.019738
17	4	2	5	3	1	4	0.008050
18	4	3	1	4	2	5	0.010650
19	4	4	2	5	3	1	0.053845
20	4	5	3	1	4	2	0.031300
21	5	1	5	4	3	2	0.053440
22	5	2	1	5	4	3	0.044617
23	5	3	2	1	5	4	0.021416
24	5	4	3	2	1	5	0.009040
25	5	5	4	3	2	1	0.061280
MR1	0.08206	0.04314	0.03371	0.03659	0.01129	0.03833	T = 1.09155
MR2	0.04702	0.03605	0.04649	0.03570	0.03577	0.06138
MR3	0.02656	0.02929	0.02977	0.05302	0.06248	0.03804
MR4	0.02472	0.04639	0.05108	0.03924	0.06341	0.04403
MR5	0.03796	0.06343	0.05725	0.05376	0.04536	0.06623
δ	0.05734	0.03414	0.02748	0.01806	0.05212	0.03455
Rank	1	4	5	6	2	3

Since the parameters a and b contribute to the normal drift coefficient μ0, we may be interested in the performances of M0 on a wide range of the coefficients of variations (CVs). Five new factor levels are given to a and b with other parameters fixed, i.e., 20%, 60%, 100%, 140% and 180%. Twenty five tests were conducted with N = 100 repeated simulations according to the values at the first three columns in Table 2. The results are shown in Figure 4. As we can see, the relative errors of M0 remain as significantly lower values than 0.05 with the CVs ranging 0~0.25. While the relative errors rise to the values around 0.3 with the CVs ranging 0.25~0.75, the robustness of M0 is still shown.

Figure 4

The correlation between CVs and RE of M0.

4.2. Real Applications

In the next two sections, two real ADT applications are used to further illustrate the advantages of M0 over M1 and M2 for both linear and nonlinear ADT analyses.

4.2.1. LED Application

Light emitting diodes (LEDs) have the merit of longer lifetime, lower power consumption, and higher brightness than traditional light sources and, thus, are widely used in the area of lighting systems. The CSADT for LEDs is conducted under two electric current levels: 35 mA and 40 mA. The normal stress level is 25 mA. Twelve LEDs are tested at each stress level and the degradation data are recorded at 50, 100, 150, 200 and 250 h. The LED will be considered to have failed when the relative degradation value of the light intensity exceeds ω = 50. For the original data, refer to Chaluvadi [36] (Table 6.3). Figure 5 shows the degradation paths for the twenty-four tested samples under two stress levels.

Figure 5

The degradation paths for twenty four LEDs under two electric current levels: (a) 35 mA and (b) 40 mA.

In [36], a linear model is used to extrapolate the pseudo-failure time of each LED at two stress levels. Then, the inverse power law is used to determine the relation between the failure time data and the electric current stresses. This procedure is similar to model M1 in our approach. As we can see from Figure 5, the tested LEDs experience nonlinear degradation paths. Therefore, a linear model like M1 may not be appropriate for ADT analysis in this case. Thus, we use model M0 to fit the data and compare the results with model M1 and M2 to verify its effectiveness on nonlinear ADT analysis. The inverse power model is used as the acceleration model, i.e.,:

Table 3 gives the estimates of unknown parameters according to the procedure in Section 3.1 and Section 3.2. It can be calculated that Pr(μ0 < 0) = 6.1870 × 10−11, which verifies that there is no danger of there being any negative drift coefficient. Figure 6 presents the fitting results. It is obvious that model M0 is the most applicable model with the maximum log-likelihood and lowest AIC value, which is the same as the fitting results. As with the linear model M1, the fitting results are worse than the other two models. Thus, the linear model in [36] is not appropriate for the LED application, while the performance of model M2 lies between model M0 and M1.

Table 3

LED application: parameter estimates for three candidate models with random drift coefficients (b ≠ 0).

Model	θ	γ	σ2	η0		η1	lmax	np	AIC	Δ
				a	b
M0	0.442	0.117	73.784	0.273	0.0018	0.677	–310	6	633	0
M1	1 (fixed)	1 (fixed)	0.761	1.69 × 10−6	5.62 × 10−14	3.112	–389	4	785	152
M2	0.450	=θ	5.840	3.52 × 10−5	2.43 × 10−11	3.112	–317	5	644	11

Figure 6

The Q-Q plots for three candidate models for the LED data.

Additionally, as shown in Table 3, the b values are close to zero, which suggests a deterministic drift coefficient without the consideration of unit-to-unit variation. Hence, we consider this situation (b = 0) for all three models. The parameters are estimated through a one-stage MLE (see Table 4).

Table 4

LED application: parameter estimates for three candidate models with deterministic drift coefficients (b = 0).

Model	θ	γ	σ2	η0 = a	η1	lmax	np	AIC	Δ
M0	0.448	0.171	45.429	0.012	1.179	–314.6634	5	639	0
M1	1 (fixed)	1 (fixed)	0.776	8.67 × 10−7	3.297	–389.7358	3	785	146
M2	0.4477	=θ	6.238	1.83 × 10−5	3.297	–319.9119	4	648	9

Compared with the AIC values in Table 3 and Table 4, it can be concluded that the models considering unit-to-unit variation (b ≠ 0) fit the LED data better than without consideration (b = 0), which indicates the presence of such variation. Additionally, among the three models, M0 should be selected for the LED ADT data analysis with the lowest AIC value.

The PDFs of the FPT for different models (b ≠ 0) at normal stress levels are shown in Figure 7. The MTTF for M0, M1, and M2 are 1167.2, 1345.0, and 13,169.2 h, while the 95% confidence intervals are [202.1, 3459.1], [364.1, 3744.1], and [1660.1, 52,787.1] h, respectively. However, the estimated lifetime using the degradation-path model in Chaluvadi [36] (p. 103) is 1346 h. It is interesting to see that the intervals of M0 and M1 can capture this value to show the consistency of evaluation results under different models, while M2 computes significantly larger values, meaning that model M2 is unreliable. A reasonable explanation may be that the time-scale transformation is inapplicable for the nonlinear ADT analysis for LEDs. With respect to model M1, its lifetime evaluation results are closer to model M0, although it fits worse on the ADT data than model M2. In Tang, et al. [18], the time-scale transformation model M2 is used to analyze the same set of data and their 95% confidence interval is [1672, 53,466] h, which, however, is also not valid in the LED case, as discussed. For the proposed model M0, its evaluation results are reliable with the best nonlinear ADT data fitting.

Figure 7

The PDFs of the FPT for different models in the LED case: (a) for M0 and M1, (b) for M2.

4.2.2. Resistor Application

Carbon-film resistors are a fixed-form type of resistor and have superior characteristics to carbon composite resistors in terms of much closer tolerances, higher maximum ohmic values, and being used in high voltage and high temperature applications. Thus, the resistance of such a resistor is affected by the temperature (s1) and applied voltage (s2). Hence, in order to evaluate their lifetimes, the CSADT is conducted under nine constant stress levels with two acceleration variables, i.e., s1 = 3.5, 4.0 and 5.0 in hundred Kelvin, s2 = 10, 15 and 20 volts. The normal stress level is s1 = 3.2315 and s2 = 5. Ten resistors are tested at each stress level. For more details about the description of resistor data, readers are referred to Park and Padgett [25]. The original data is modified to ensure that the initial degradation values are equal to zero and the threshold ω is, therefore, set to be 0.2. The CSADT data at the first stress level is presented in Figure 8. It is clear that the degradation processes follow linear paths which are the same as the data in other stress levels.

Figure 8

The degradation paths for ten resistors in CSADT (s1 = 3.5 and s2 = 10).

Regarding the acceleration model with the unit-to-unit variation, the following model with two acceleration variables is used [25]:

Table 5 presents the estimates of the unknown parameters for different models on resistor CSADT data. We also calculated that Pr(μ0 < 0) = 6.0094 × 10−6, which verifies that there is no danger of there being any negative drift coefficient. Intuitively, compared with M1 and M2, model M0 displays the best fit with the maximum log-likelihood value and minimum AIC value. However, according to the Δ values, there is substantial support for model M2 and considerably less support for model M1. The fitting results for the different models are shown in Figure 9, which imply that the performances of those models are similar in the resistor case since that the degradation paths are approximately linear.

Table 5

Resistor application: parameter estimates for three candidate models with random drift coefficients (b ≠ 0).

Model	θ	γ	σ2	η0		η1	η2	lmax	np	AIC	Δ
				a	b
M0	1.076	0.918	3.02 × 10−4	8.96 × 10−4	4.19 × 10−8	0.462	0.108	1698	7	–3383	0
M1	1 (fixed)	1 (fixed)	2.59 × 10−4	0.0012	7.02 × 10−8	0.440	0.102	1694	5	–3378	4.8
M2	1.046	=θ	2.35 × 10−4	0.0010	5.12 × 10−8	0.451	0.106	1697	6	–3381	1.6

Figure 9

The Q-Q plots for three candidate models for the resistor data.

Herein, we also checked the models without the consideration of unit-to-unit variation. The results are listed in Table 6. Compared with the results in Table 5, it can be seen that models with the normal drift coefficients (b ≠ 0) have lower AIC values, which means that the unit-to-unit variation should be considered. It is also interesting to see that, from Table 6, M1 is the best model with a deterministic drift coefficient, while supports are given to the other two models.

Table 6

Resistor application: parameter estimates for three candidate models with deterministic drift coefficients (b = 0).

Model	θ	γ	σ2	η0 = a	η1	η2	lmax	np	AIC	Δ
M0	1.020	0.9498	3.33 × 10−4	0.0011	0.445	0.104	1647	6	–3282	2.7
M1	1 (fixed)	1 (fixed)	3.02 × 10−4	0.0012	0.439	0.103	1646	4	–3284	0
M2	1.008	=θ	2.97 × 10−4	0.0011	0.442	0.103	1646	5	–3283	1.8

We may also be interested in whether M0 can be simplified into M2 with the assumption (θ = γ) since M2 has substantial support. Hence, the likelihood ratio (LR) test is implemented with the log-likelihood values in Table 5. The resulting LR statistic is 3.57 (< = 3.84). Thus, we accept the assumption and choose M2.

The PDFs and CDFs of the FPT for different models are given in Figure 10a,b, which are almost identical with the MTTF at around 2400 h, and the PDF of the FPT for model M0 is slightly sharper than that of M1 and M2.

Figure 10

The (a) PDFs and (b) CDFs of the FPT for three candidate models in the resistor case.

Comparing Section 4.1 with Section 4.2, it can be concluded that the proposed method can be effectively used for ADT analysis for both linear and nonlinear scenarios with single and multiple acceleration variables.

5. Conclusions

This paper has proposed a general ADT model based on the Wiener process and provided statistical analysis methods for unknown parameter estimation. The general ADT model is suitable for both linear and nonlinear ADT analyses, with single and multiple acceleration variables, and considers the unit-to-unit variation among the tested samples and temporal variation of the degradation processes simultaneously. The simulation example demonstrates that the general model is robust with respect to ADT analysis and its reliability evaluation results become more accurate when the sample size increases. Furthermore, the LED and resistor cases have verified its effectiveness and superiority to real engineering applications over the commonly used linear and time-scale transformation Wiener process models, which not only fit well with the degradation data at all accelerated stress levels, but also computes reliable lifetime and reliability evaluation results.

This study focused on the modeling procedure for ADT data based on the Wiener process. However, other methods, e.g., the Gamma or inverse Gaussian processes, are worth further research when the degradation paths are monotonic.