# Reliability Analysis with Multiple Dependent Features from a Vibration-Based Accelerated Degradation Test

Authors: Fuqiang Sun; Jingcheng Liu; Xiaoyang Li; Haitao Liao
Journal: Shock and Vibration
 (2016)

---

## Abstract

Accelerated degradation testing (ADT) has been widely used for reliability prediction of highly reliable products. In many applications, ADT data consists of multiple degradation-related features, and these features are usually dependent. When dealing with such ADT data, it is important to fully utilize the multiple degradation features and take into account their inherent dependency. This paper proposes a novel reliability-assessment method that combines Brownian motion and copulas to model ADT data obtained from vibration signals. In particular, degradation feature extraction is first carried out using the raw vibration signals, and a feature selection method quantifying feature properties, such as trendability, monotonicity, and robustness, is adopted to determine the most suitable degradation features. Then, a multivariate s-dependent ADT model is developed, where a Brownian motion is used to depict the degradation path of each degradation feature and a copula function is employed to describe the dependence among these degradation features. Finally, the proposed ADT model is demonstrated using the vibration-based ADT data for an electric motor.


---

## Body

## 1. Introduction

Nowadays, products are made more reliable due to the advances of design and manufacturing and to the improvement of materials technology. As a result, traditional reliability analysis methods based on lifetime data are becoming less capable of evaluating the reliability of such highly reliable products. In practice, the failure mechanism of a product often can be traced to the degradation of the product’s key performance parameters, such as wear, fatigue crack, and material corrosion. The product is considered to be failed if any of the performance parameters exceeds its critical value. Therefore, an alternative way for evaluating the product’s reliability is to analyze the degradation data obtained from a degradation test. To collect degradation data in a short time, accelerated degradation testing (ADT) technique has been proposed, where some units of the product are exposed to harsher-than-normal conditions in order to speed up product degradation. The resulting ADT data will be used to develop an ADT model that enables the extrapolation of the product’s reliability under the normal operating conditions.In order to find the degradation trend of a product and predict its lifetime using ADT, the product’s performance indicators must be defined first. Such indicators may include not only the product’s functional or performance parameters, but also those special features extracted from raw sensor data [1], such as vibration, force and acoustic signals, temperature, and voltage. In particular, vibration monitoring has been widely used in fault diagnostics and remaining useful life (RUL) prediction for mechanical and electromechanical systems [2]. Nectoux et al. [3] developed a bearing experimental platform, called PRONOSTIA, to perform bearing ADT, where both vibration and temperature signals were monitored. Liu et al. [4] introduced the confidence value (CV) as the degradation index for a product and used a fuzzy regression model to analyze the bearing vibration data given in [3]. However, in the field of prognostic and health management (PHM), little effort has been focused on the mathematical modeling of multiple degradation features obtained from vibration-based ADT.Generally, raw vibration data may not be directly used for degradation modeling. As a solution, signal processing techniques are usually utilized to extract some features from the vibration data [5]. It is not uncommon that quite a few features can be extracted in the time, frequency, and time-frequency domains as well as via information entropy analysis [3]. However, since not all features are required in describing the product’s degradation process, it is necessary to perform dimension reduction for the purpose of degradation modeling. There are two popular ways to do that. The first one is to perform principal component analysis, singular value decomposition, self-organizing map, or clustering [2]. The other one is feature selection that selects a subset of features by assessing the fitness of candidate features [6–10]. The purpose of performing feature selection is to identify those features that can best describe the fault progression process [8]. Coble and Hines [6] proposed three metrics for feature selection, that is, trendability, monotonicity, and prognosability. Javed et al. [9] performed feature selection by evaluating feature fitness using the monotonicity metric defined in [6] and a new trendability metric, and Liao [8] only used the monotonicity criterion. Camci et al. [10] divided a degradation feature into segments and the average separability of all segments was defined for feature selection. Zhang et al. [7] defined three other metrics, that is, correlation, monotonicity, and robustness, for selecting more relevant degradation features.It is quite common that there is an s-dependent relationship among the stochastic processes describing different degradation features. Clearly, ignoring such dependency in ADT analysis will lead to inaccurate reliability estimates. For modeling multiple degradation processes, most work considers an independence assumption or multivariate normal distributions [11, 12]. However, these assumptions may not be appropriate in many engineering applications. Copulas provide a convenient way to model the dependence of random variables characterized by their marginal distributions [13]. In addition, no constraints on the univariate marginal distributions are required in implementing copula. Because of its advantages, copula has become a popular approach in multivariate statistical analysis. Sari et al. [14] used a copula function to describe the correlation between two performance characteristics of light emitting diode (LED) and combined it with a generalized linear regression model to model bivariate constant-stress ADT. Pan et al. [15] discussed the bivariate degradation modeling approaches based on Wiener processes and copulas under constant-stress accelerating scenarios as well. Similarly, Peng et al. [16], Liu et al. [17], Wang and Pham [18], and Li and Xue [19] also employed the copula method to model the dependency between bivariate degradation features. However, previous research is limited on modeling the dependence between bivariate degradation processes, and little effort has been focused on multiple dependent degradation processes.Another key issue in ADT modeling is how to choose appropriate mathematic models to depict the degradation paths. In the literature, two kinds of degradation models that have been widely used are stochastic process model (e.g., Brownian motion, gamma process, and inverse Gaussian process) and general path model [20]. Ye and Xie [20] conducted a comprehensive literature review on Brownian motion based degradation models. Liao and Elsayed [21] developed an ADT model based on Brownian motion to predict the field reliability of LED considering stress variations. Guan et al. [22] adopted Brownian motion to model the constant-stress ADT (CSADT) and estimated the model parameters using an objective Bayesian method.Based on the related literature, one can see that some technical problems have not been resolved in the study of vibration-based ADT with multiple features, including fitness analysis of degradation indicators, development of multivariate s-dependent ADT model, and parameter estimation. To overcome these challenges, a novel reliability modeling method for vibration-based CSADT with multiple dependent features is proposed in this paper. To the best of our knowledge, this study is the first attempt in the related areas. The remainder of this paper is organized as follows. The framework of the proposed method is given in Section2. Section 3 introduces the feature exaction and selection methods. Section 4 provides a multivariate ADT model based on copulas. In Section 5, a case study is used to illustrate the validity of the proposed model. Finally, conclusions are drawn in Section 6.


## 2. Framework of the Proposed Model

Figure1 provides the proposed framework for reliability and life evaluation using vibration-based ADT data. The framework consists of two parts: fitness analysis of degradation indicators and modeling of multivariate ADT data.Figure 1
The framework of reliability and life evaluation for vibration-based ADT.The first part is to extract and select the degradation indicators from the raw vibration data. It is often necessary to first denoise the vibration data to reduce the negative influence of environmental noise on the subsequent analysis. Then, features are extracted in four domains, that is, time, frequency, time-frequency, and information domains. Since the dimension of the features is high, three goodness metrics are used for the selection of the most suitable degradation indicators from these candidate features.The other part is to model the multiple accelerated degradation processes describing the selected features. At first, a Brownian motion model is adopted to depict the degradation path of each feature and to extrapolate the reliability related to each feature under the normal stress with an acceleration model and a specified failure threshold. Afterwards, a copula function is used to describe the dependence of different degradation features resulting in a multivariate s-dependent ADT model. A two-stage statistical inference procedure, called Inference Functions for Margins (IFM) method, is employed to estimate the unknown parameters of the proposed model. Finally, Akaike information criterion (AIC) is used to determine the best fitting copula function, and the system joint reliability under the normal operating condition can be obtained.The following sections will elaborate on the analysis of vibration-based ADT with multiple features.


## 3. The Degradation Indicator Fitness Analysis

### 3.1. Signal Processing and Feature Extraction

Denoising is an important step in effective fault detection, diagnosis, and prognostics [23]. In this paper, wavelet automatic denoising technique [24] is used to remove the noise from vibration signals prior to feature extraction. This technique is capable of deconstructing complex signals into basis signals of finite bandwidth and then reconstructing them again with very little loss of information, which is best for nonstationary signal analysis.Then, different signal processing methods are utilized to find the most valuable features relevant to the product’s performance state. Indeed, feature extraction via signal processing is considered one of the major elements and the most challenging aspect in analyzing vibration-based ADT data [25]. The existing signal processing and feature extraction techniques dealing with vibration signals mainly include four categories [2].(  1) Time Domain. Time domain feature extraction is performed directly on the waveforms of vibration signals, which are suitable for fault detection. Usually, the extracted features include mean, variance, waveform factor, crest factor, clearance factor, and kurtosis factor.(  2) Frequency Domain. Frequency domain analysis can identify and isolate those important frequency components of the signals. The most popular method is the fast Fourier transform, with which features such as center frequency, mean square frequency, and frequency variance can be obtained.(  3) Time-Frequency Domain. Time-frequency analysis aims at investigating the signals in both time and frequency domains [9]. In particular, Hilbert-Huang Transform (HHT), based on empirical mode decomposition, is the most powerful technique for the analysis of nonstationary signals [26]. HHT marginal energy is calculated as the feature in this domain.(  4) Information Entropy. Information entropy is a measure of uncertainty of signals [4]. A higher value indicates a higher uncertainty of the signal. Following this idea, the Hilbert marginal energy spectrum entropy is computed in feature extraction.Letx1,x2,…,xN be the recorded vibration signals, f be the frequency components of the signals, and S(f) be the power spectrum. Table 1 gives various features obtained by processing the original vibration signals.Table 1
Candidate degradation features.
Feature Formula Feature Formula Mean x - = 1 N  ∑ i = 1 N x i Clearance factor L = max ⁡ x i 1 / N ∑ i = 1 N x i 2 Variance σ 2 = 1 N  ∑ i = 1 N x i - x - 2 Frequency center F C = ∫ 0 + ∞ f S f d f ∫ 0 + ∞ S f d f Crest factor C = max ⁡ x i 1 / N ∑ i = 1 N x i 2 Frequency variance V F = ∫ 0 + ∞ f - F C 2 S f d f ∫ 0 + ∞ S f d f Kurtosis factor α = 1 N  ∑ i = 1 N x i - x - 4 σ 4 HHT marginal energy (HME) H M E ω =  ∫ 0 T H 2 ω , t d t H ω , t = R e  ∑ j = 1 m c j t e i ∫ 2 π ω j t d t Waveform factor S = 1 / N ∑ i = 1 N x i 2 1 / N ∑ i = 1 N x i HHT marginal energy spectrum entropy (HMESE) H M E E = - ∑ j = 1 m p j log 2 ⁡ p j log 2 ⁡ m, p j = H M E j ∑ j = 1 m H M E j, ⁢  ∑ j = 1 m p j = 1 Note:m is the number of frequency components in Hilbert-Huang Transform (HHT).


### 3.2. Feature Selection

The selection of a subset of appropriate features is essential for prognostic modeling. It is often performed on the basis of feature characteristics by comparing candidate features with respect to their fitness. In other words, the features with high attribute values, such as monotonicity, correlation, trendability, and robustness, should lead to more accurate prognostic results and thus be included for further degradation modeling.Technically, after removing unwanted noise, a feature is decomposed into its trend and random part by applying a robust locally weighted scatter plot smooth (RLOESS) with a span value of 0.5 (i.e., 50%), as shown in (1). Unlike locally weighted scatter plot smooth (LOESS) that is a popular smoothing method based on a locally weighted regression function and a quadratic polynomial, RLOESS uses a robust regression weight function to make the LOESS process insensitive to outliers(1)x=xT+xR,where x is the actual feature value, xT is its trend value, and xR is the residual.Then, a set of mathematical criteria in favor of the features with high attribute values are considered in feature selection.(1)
Monotonicity characterizes the underlying increasing or decreasing trend of feature. It can be measured by the average absolute difference of derivatives for each feature [6] as follows:(2)M=meanNo.ofd/dx>0n-1-No.ofd/dx<0n-1,wheren is the number of observations, d/dx is the derivative, and mean represents averaging. The value of M is from 0 to 1, with M=1 indicating highly monotonic features and M=0 indicating nonmonotonic features.(2)
Trendability is related to the linear relationship between the feature and the time. Clearly, it shows how the degrading state varies with time. The metric of trendability is given by [9](3)T=n∑xT·t-∑xT∑tn∑xT2-∑xT2n∑t2-∑t2,whereT∈[0,1] is the absolute correlation coefficient between feature trend xT and time t. Obviously, a constant function of feature has zero correlation with time, while a linear function indicates a strong correlation between the feature and the time.(3)
Robustness reflects the tolerance of the feature to outliers [7], which is given by(4)R=1n∑exp⁡-xRx.One can see that all the three metrics of features are confined in the range of[0,1]. Due to the fact that one metric only partially measures the suitability of candidate feature, a weighted sum of the three metrics is used in this paper as the fitness function for feature selection(5)max⁡fitness=w1·M+w2·T+w3·Rs.t.ωi>0∑iωi=1,i=1,2,3,where constants ω1, ω2, and ω3 control the importance of each metric in the optimization problem.


## 3.1. Signal Processing and Feature Extraction

Denoising is an important step in effective fault detection, diagnosis, and prognostics [23]. In this paper, wavelet automatic denoising technique [24] is used to remove the noise from vibration signals prior to feature extraction. This technique is capable of deconstructing complex signals into basis signals of finite bandwidth and then reconstructing them again with very little loss of information, which is best for nonstationary signal analysis.Then, different signal processing methods are utilized to find the most valuable features relevant to the product’s performance state. Indeed, feature extraction via signal processing is considered one of the major elements and the most challenging aspect in analyzing vibration-based ADT data [25]. The existing signal processing and feature extraction techniques dealing with vibration signals mainly include four categories [2].(  1) Time Domain. Time domain feature extraction is performed directly on the waveforms of vibration signals, which are suitable for fault detection. Usually, the extracted features include mean, variance, waveform factor, crest factor, clearance factor, and kurtosis factor.(  2) Frequency Domain. Frequency domain analysis can identify and isolate those important frequency components of the signals. The most popular method is the fast Fourier transform, with which features such as center frequency, mean square frequency, and frequency variance can be obtained.(  3) Time-Frequency Domain. Time-frequency analysis aims at investigating the signals in both time and frequency domains [9]. In particular, Hilbert-Huang Transform (HHT), based on empirical mode decomposition, is the most powerful technique for the analysis of nonstationary signals [26]. HHT marginal energy is calculated as the feature in this domain.(  4) Information Entropy. Information entropy is a measure of uncertainty of signals [4]. A higher value indicates a higher uncertainty of the signal. Following this idea, the Hilbert marginal energy spectrum entropy is computed in feature extraction.Letx1,x2,…,xN be the recorded vibration signals, f be the frequency components of the signals, and S(f) be the power spectrum. Table 1 gives various features obtained by processing the original vibration signals.Table 1
Candidate degradation features.
Feature Formula Feature Formula Mean x - = 1 N  ∑ i = 1 N x i Clearance factor L = max ⁡ x i 1 / N ∑ i = 1 N x i 2 Variance σ 2 = 1 N  ∑ i = 1 N x i - x - 2 Frequency center F C = ∫ 0 + ∞ f S f d f ∫ 0 + ∞ S f d f Crest factor C = max ⁡ x i 1 / N ∑ i = 1 N x i 2 Frequency variance V F = ∫ 0 + ∞ f - F C 2 S f d f ∫ 0 + ∞ S f d f Kurtosis factor α = 1 N  ∑ i = 1 N x i - x - 4 σ 4 HHT marginal energy (HME) H M E ω =  ∫ 0 T H 2 ω , t d t H ω , t = R e  ∑ j = 1 m c j t e i ∫ 2 π ω j t d t Waveform factor S = 1 / N ∑ i = 1 N x i 2 1 / N ∑ i = 1 N x i HHT marginal energy spectrum entropy (HMESE) H M E E = - ∑ j = 1 m p j log 2 ⁡ p j log 2 ⁡ m, p j = H M E j ∑ j = 1 m H M E j, ⁢  ∑ j = 1 m p j = 1 Note:m is the number of frequency components in Hilbert-Huang Transform (HHT).


## 3.2. Feature Selection

The selection of a subset of appropriate features is essential for prognostic modeling. It is often performed on the basis of feature characteristics by comparing candidate features with respect to their fitness. In other words, the features with high attribute values, such as monotonicity, correlation, trendability, and robustness, should lead to more accurate prognostic results and thus be included for further degradation modeling.Technically, after removing unwanted noise, a feature is decomposed into its trend and random part by applying a robust locally weighted scatter plot smooth (RLOESS) with a span value of 0.5 (i.e., 50%), as shown in (1). Unlike locally weighted scatter plot smooth (LOESS) that is a popular smoothing method based on a locally weighted regression function and a quadratic polynomial, RLOESS uses a robust regression weight function to make the LOESS process insensitive to outliers(1)x=xT+xR,where x is the actual feature value, xT is its trend value, and xR is the residual.Then, a set of mathematical criteria in favor of the features with high attribute values are considered in feature selection.(1)
Monotonicity characterizes the underlying increasing or decreasing trend of feature. It can be measured by the average absolute difference of derivatives for each feature [6] as follows:(2)M=meanNo.ofd/dx>0n-1-No.ofd/dx<0n-1,wheren is the number of observations, d/dx is the derivative, and mean represents averaging. The value of M is from 0 to 1, with M=1 indicating highly monotonic features and M=0 indicating nonmonotonic features.(2)
Trendability is related to the linear relationship between the feature and the time. Clearly, it shows how the degrading state varies with time. The metric of trendability is given by [9](3)T=n∑xT·t-∑xT∑tn∑xT2-∑xT2n∑t2-∑t2,whereT∈[0,1] is the absolute correlation coefficient between feature trend xT and time t. Obviously, a constant function of feature has zero correlation with time, while a linear function indicates a strong correlation between the feature and the time.(3)
Robustness reflects the tolerance of the feature to outliers [7], which is given by(4)R=1n∑exp⁡-xRx.One can see that all the three metrics of features are confined in the range of[0,1]. Due to the fact that one metric only partially measures the suitability of candidate feature, a weighted sum of the three metrics is used in this paper as the fitness function for feature selection(5)max⁡fitness=w1·M+w2·T+w3·Rs.t.ωi>0∑iωi=1,i=1,2,3,where constants ω1, ω2, and ω3 control the importance of each metric in the optimization problem.


## 4. Modeling of CSADT with Multiple Features

### 4.1. The Univariate Accelerated Degradation Model

The premise of ADT-based reliability evaluation is to choose an appropriate model based on the degradation physics or the degradation data [20]. We first focus on each degradation feature.


#### 4.1.1. Basic Assumptions

Among many stochastic process models, Brownian motion is the most widely used in degradation modeling and analysis. This paper employs this model to describe each degradation process. To evaluate the reliability characteristics of the product under the normal operating condition based on ADT data, the relationship between lifetime characteristics and accelerated stress levels should be modeled. To this end, the following assumptions are made.(A1)
The test specimens are independent, and no catastrophic failures occur during the test.(A2)
Under either the normal conditionS0 or accelerated conditions S1<S2<⋯<SK, degradation process Y(t) over time can be described by(6)Yt=σBt+μ·t+y0,wheret is time, B(t) is a standard Brownian motion, y0 is the initial degradation level at time zero, μ>0 is the drift parameter, and σ>0 is the diffusion parameter describing the variability due to unit-to-unit variation and the variation of operating and environmental conditions.The drift parameterμ, commonly known as the degradation rate, is regarded as a function of stress conditions. Usually, such acceleration relationships, such as Arrhenius model, Eyring model, or inverse power model [27], are obtained from physical mechanism analysis or empirical experiences. In particular, the following log-linear relationship is considered: (7)ln⁡μ=ln⁡dS=a+b·φS,wherea and b are unknown parameters, d(S) is the degradation rate under stress level S, and φ(S) is a function of stress level depending on the stress type. For example, if S is absolute temperature, φ(S)=1/S; if S is voltage, φ(S)=ln⁡S.(A3)
The diffusion parameterσ does not change with stress, that is, σ0=σ1=⋯=σK=σ for the normal condition S0 and accelerated conditions S1<S2<⋯<SK.


#### 4.1.2. Reliability Function for Each Degradation Process

Because of the independent increment property of Brownian motion, the degradation incrementΔy over nonoverlapped interval Δt is normally distributed with mean μΔt and variance σ2Δt. Its probability density function (PDF) is(8)fΔy=12πσ2Δtexp⁡-Δy-μ·Δt22σ2Δt.Generally, a product is supposed to be failed if its performance degradation exceeds a critical threshold for the first time. In other words, the first passage time (FPT) distribution is used as the product’s lifetime distribution. For a given critical thresholdD, the lifetime TD of the product is the instant when degradation process Y(t) exceeds D for the first time; that is,(9)TD=inf⁡t>0∣Yt≥D.It is well-known that the first passage time follows the inverse Gaussian distribution with the following PDF [28]:(10)ft;y0,D=D-y0σ2πt3exp-D-y0-μ·t22σ2t.Combined with (7), the associated reliability function under normal operational condition S0 can be expressed as(11)Rt=ΦD-y0-dS0tσt-exp⁡2dS0D-y0σ2Φ-D-y0+dS0tσt,where Φ(·) is the cumulative distribution function of standard normal distribution.


### 4.2. Use of Copulas for Multiple Features

#### 4.2.1. Definition and Basic Properties

Sklar’s theorem lays the theoretical foundation for copulas and their applications [13].Theorem 1 (Sklar’s theorem [29]).
LetX=(X1,…,XN) be a random vector with marginal distributions F1(x1),…,FN(xN) and F be their joint distribution function. Then, there exists a copula function C such that(12)Fx1,…,xn,…,xN=CF1x1,…,Fnxn,…,FNxN.If F1(x1),…,FN(xN) are continuous, the copula function C is unique. Conversely, if F1(x1),…,FN(xN) are univariate distributions, the function F defined by (12) is the joint distribution function associated with margins F1(x1),…,FN(xN).Based on the theorem, a multivariate copula function can be defined as follows.Definition 2 (n-dimensional copula or n-copula [13]).
AnN-dimensional copula is a function C mapping from IN=0,1N to I with the following properties:(1)
Grounded: for everyu=(u1,…,uN) in IN, C(u)=0 if at least one coordinate of u is 0.(2)
Uniform marginals: if all coordinates ofu are ones except uk, then(13)Cu=C1,…,1,uk,1,…,1=uk.(3)
N-increasing: for each hyperrectangleB=∏i=1Nui,vi⊆[0,1]N, the C-volume of B is nonnegative:(14)∫BdCu,v=∑z∈×i=1Nui,vi-1NzCz≥0,whereN(z)=#k:zk=uk.
Furthermore, the PDF corresponding to joint distribution functionF(x1,…,xn,…,xN) can be calculated as(15)fx1,…,xn,…,xN=cF1x1,…,Fnxn,…,FNxN∏n=1Nfnxn,where fn(·) is the PDF of a marginal distribution Fn(·), and c(·) is the density of copula function C:(16)cu1,…,un,…,uN=∂NCu1,…,un,…,uN∂u1,…,∂un,…,∂uN∀u=u1,…,un,…,uN∈IN.


#### 4.2.2. Selection of Copula Family

It is important to select an appropriate copula family for the given data, as the selection significantly influences the prediction performance. Table2 presents five popular multivariate copula families, where u1,…,un are n random variables.Table 2
Typical multivariate copulas.
Copulas C u 1 , u 2 , … , u n Parameters Gaussian copula Φ ρ Φ - 1 u 1 , Φ - 1 u 2 , … , Φ - 1 u n ρ ∈ - 1,1 Student’st-copula t ρ , v t v - 1 u 1 , t v - 1 u 2 , … , t v - 1 u n ρ ∈ - 1,1 k > 2 Clayton copula ∑ i = 1 n u i - δ - n + 1 - 1 / δ δ > 0 Frank copula - δ - 1 log ⁡ 1 +  ∏ i = 1 n exp ⁡ - δ u i - 1 exp ⁡ - δ - 1 n - 1 δ ∈ - ∞ , ∞ ∖ 0 δ > 0 for n ≥ 3 Joe copula 1 - 1 -  ∏ i = 1 n 1 - 1 - u i δ 1 / δ δ ∈ 1 , ∞To balance the statistical fitting and model complexity, the values of classic Akaike information criterion (AIC) are computed for the candidate copula families:(17)AIC=2k-2ln⁡L,where k is the number of estimated parameters in the model, and L is the value of the corresponding maximum likelihood. By comparing the AIC values, the family with the minimum value will be selected.


### 4.3. Multivariate ADT Model

To analyze the multivariate CSADT data, additional assumptions are considered.(A4)
For a product with multiple degradation indicators, the dependency among the features can be characterized by a copula function.(A5)
The degradation measurements of all specimens are collected at the same time.(A6)
The specimen is considered to be failed if one of the features reaches its corresponding failure threshold for the first time.Suppose that a product hasP degradation features and each of the accelerated degradation processes is modeled by a drift Brownian motion. Yp(t), p=1,2,…,P, denotes the degradation process of the pth characteristic index at operating condition S0, and its corresponding failure threshold is Dp. Let the failure time of the pth characteristic index be Tp, and the lifetime of the product be T, so T=min(T1,…,Tp). Clearly, the system reliability under normal condition S0 can be expressed as(18)Rt=PT>t=PT1>t,…,TP>t=PY1t<D1,…,YPt<DP.Apparently, if theP degradation features are independent, (18) can be rewritten as(19)Rt=PY1t<D1×⋯×PYPt<DP=R1t×⋯×RPt,where Rp(t), p=1,2,…,P is the pth characteristic index’s reliability function under normal operating condition S0, as defined in (11). However, this is often not true in many engineering applications. To address the dependence of P degradation features, a copula function can be utilized. In particular, the product’s reliability at time t under normal operating condition S0 can be performed as(20)Rt=CR1t,R2t,…,RPt;θ,where θ is the parameter set of the copula.


### 4.4. Parameter Estimation

Consider a copula-based multivariate distribution for the random vector[Y1,…,YP], with joint PDF(21)fy;α1,…,αP,θ=cF1y1;α1,…,FPyP;αP;θ∏p=1Pfpyp;αp,where Fp(·;αp) and fp(·;αp), p=1,…,P, are the pth marginal CDF and PDF with parameter set αp, and c(·;θ) is the density of copula function C(·;θ) with parameter θ.For a sample of sizen, the full log-likelihood function can be expressed as(22)ln⁡Lα1,…,αP,θ=∑i=1nln⁡cF1y1i;α1,…,FPyPi;αP;θ︸dependence  structureLC+∑p=1P∑i=1nfpypi;αp︸marginals∑p=1PLp.The maximum likelihood estimates (MLE) of the parameters(α1,…,αp,θ) can be obtained by maximizing (22). However, as the dimension of parameters increases, it is difficult to get the optimal solution. Joe [30, 31] proposed the IFM method as a computationally attractive alternative to estimate the parameters in multivariate copula models. Technically, (22) is first decomposed into two parts. The first part, denoted by LC, is the contribution from the dependence structure in data, and the other part is the contributions from each margin, Lp, p=1,…,P. Then, The IFM method separately estimates the parameters of marginal distributions and the parameters of the copula in two stages. In the first stage, the estimates of the marginal distribution parameters are obtained from the log-likelihood Lp of each margin. In the second stage, the copula parameter estimates are computed by maximizing LC based on the marginal distribution parameters estimated in the first stage. In this paper, we use this method to estimate the parameters of the proposed multivariate ADT model.


#### 4.4.1. Marginal Distribution Parameters Estimation

Suppose that in a CSADT withK stress levels the total sample size is n. There are nl specimens under stress level Sl. During the CSADT, all specimens are measured once in every Δt time interval and there are Ml inspections under Sl. Then, the observation of the pth characteristic index at time tlij is yp(tlij), p=1,…,P, l=1,…,K, i=1,…,nl, j=1,…,Ml, where tlij is the time of the jth measurement of the ith unit under the lth stress level. According to (8), the likelihood function of the pth characteristic indices is given by(23)Lp=∏l=1K∏i=1nl∏j=1Ml-112πσp2Δtexp⁡-yptlij+1-yptlij-μpl·Δt22σp2Δt,p=1,…,P.Considering the acceleration model given by (7), μpl can be expressed as μpl=dp(Sl)=exp⁡(ap+bpφ(Sl)), where ap and bp are unknown parameters of the pth characteristic index. So the model parameter set of the characteristic index is αp=(ap,bp,σp), p=1,…,P. Let Δyplij=yp(tlij+1)-yp(tlij); then (23) can be expressed as(24)ln⁡L=-12∑l=1K∑i=1nl∑j=1Ml-1ln⁡2πσp2Δt+Δyplij-exp⁡ap+bpφSl·Δt2σp2Δt,p=1,…,Pand the parameter set αp=(ap,bp,σp) of each marginal distribution can be estimated by maximizing the log-likelihood function.


#### 4.4.2. Copula Parameter Estimation

By replacing the marginal parametersαp, p=1,…,P, with their first-stage estimateα^p, the copula parameters can be estimated by maximizing the copula likelihood contribution. Then, the copula parameter θ in (20) is given by(25)θ^=argmax∑p=1P∑i=1nln⁡cR1y1i;α^1,…,RPypi;α^P;θ.In addition, the AIC method mentioned in Section 4.2.2 is used to quantitatively select the best fitting model from candidate copula functions.


## 4.1. The Univariate Accelerated Degradation Model

The premise of ADT-based reliability evaluation is to choose an appropriate model based on the degradation physics or the degradation data [20]. We first focus on each degradation feature.


### 4.1.1. Basic Assumptions

Among many stochastic process models, Brownian motion is the most widely used in degradation modeling and analysis. This paper employs this model to describe each degradation process. To evaluate the reliability characteristics of the product under the normal operating condition based on ADT data, the relationship between lifetime characteristics and accelerated stress levels should be modeled. To this end, the following assumptions are made.(A1)
The test specimens are independent, and no catastrophic failures occur during the test.(A2)
Under either the normal conditionS0 or accelerated conditions S1<S2<⋯<SK, degradation process Y(t) over time can be described by(6)Yt=σBt+μ·t+y0,wheret is time, B(t) is a standard Brownian motion, y0 is the initial degradation level at time zero, μ>0 is the drift parameter, and σ>0 is the diffusion parameter describing the variability due to unit-to-unit variation and the variation of operating and environmental conditions.The drift parameterμ, commonly known as the degradation rate, is regarded as a function of stress conditions. Usually, such acceleration relationships, such as Arrhenius model, Eyring model, or inverse power model [27], are obtained from physical mechanism analysis or empirical experiences. In particular, the following log-linear relationship is considered: (7)ln⁡μ=ln⁡dS=a+b·φS,wherea and b are unknown parameters, d(S) is the degradation rate under stress level S, and φ(S) is a function of stress level depending on the stress type. For example, if S is absolute temperature, φ(S)=1/S; if S is voltage, φ(S)=ln⁡S.(A3)
The diffusion parameterσ does not change with stress, that is, σ0=σ1=⋯=σK=σ for the normal condition S0 and accelerated conditions S1<S2<⋯<SK.


### 4.1.2. Reliability Function for Each Degradation Process

Because of the independent increment property of Brownian motion, the degradation incrementΔy over nonoverlapped interval Δt is normally distributed with mean μΔt and variance σ2Δt. Its probability density function (PDF) is(8)fΔy=12πσ2Δtexp⁡-Δy-μ·Δt22σ2Δt.Generally, a product is supposed to be failed if its performance degradation exceeds a critical threshold for the first time. In other words, the first passage time (FPT) distribution is used as the product’s lifetime distribution. For a given critical thresholdD, the lifetime TD of the product is the instant when degradation process Y(t) exceeds D for the first time; that is,(9)TD=inf⁡t>0∣Yt≥D.It is well-known that the first passage time follows the inverse Gaussian distribution with the following PDF [28]:(10)ft;y0,D=D-y0σ2πt3exp-D-y0-μ·t22σ2t.Combined with (7), the associated reliability function under normal operational condition S0 can be expressed as(11)Rt=ΦD-y0-dS0tσt-exp⁡2dS0D-y0σ2Φ-D-y0+dS0tσt,where Φ(·) is the cumulative distribution function of standard normal distribution.


## 4.1.1. Basic Assumptions

Among many stochastic process models, Brownian motion is the most widely used in degradation modeling and analysis. This paper employs this model to describe each degradation process. To evaluate the reliability characteristics of the product under the normal operating condition based on ADT data, the relationship between lifetime characteristics and accelerated stress levels should be modeled. To this end, the following assumptions are made.(A1)
The test specimens are independent, and no catastrophic failures occur during the test.(A2)
Under either the normal conditionS0 or accelerated conditions S1<S2<⋯<SK, degradation process Y(t) over time can be described by(6)Yt=σBt+μ·t+y0,wheret is time, B(t) is a standard Brownian motion, y0 is the initial degradation level at time zero, μ>0 is the drift parameter, and σ>0 is the diffusion parameter describing the variability due to unit-to-unit variation and the variation of operating and environmental conditions.The drift parameterμ, commonly known as the degradation rate, is regarded as a function of stress conditions. Usually, such acceleration relationships, such as Arrhenius model, Eyring model, or inverse power model [27], are obtained from physical mechanism analysis or empirical experiences. In particular, the following log-linear relationship is considered: (7)ln⁡μ=ln⁡dS=a+b·φS,wherea and b are unknown parameters, d(S) is the degradation rate under stress level S, and φ(S) is a function of stress level depending on the stress type. For example, if S is absolute temperature, φ(S)=1/S; if S is voltage, φ(S)=ln⁡S.(A3)
The diffusion parameterσ does not change with stress, that is, σ0=σ1=⋯=σK=σ for the normal condition S0 and accelerated conditions S1<S2<⋯<SK.


## 4.1.2. Reliability Function for Each Degradation Process

Because of the independent increment property of Brownian motion, the degradation incrementΔy over nonoverlapped interval Δt is normally distributed with mean μΔt and variance σ2Δt. Its probability density function (PDF) is(8)fΔy=12πσ2Δtexp⁡-Δy-μ·Δt22σ2Δt.Generally, a product is supposed to be failed if its performance degradation exceeds a critical threshold for the first time. In other words, the first passage time (FPT) distribution is used as the product’s lifetime distribution. For a given critical thresholdD, the lifetime TD of the product is the instant when degradation process Y(t) exceeds D for the first time; that is,(9)TD=inf⁡t>0∣Yt≥D.It is well-known that the first passage time follows the inverse Gaussian distribution with the following PDF [28]:(10)ft;y0,D=D-y0σ2πt3exp-D-y0-μ·t22σ2t.Combined with (7), the associated reliability function under normal operational condition S0 can be expressed as(11)Rt=ΦD-y0-dS0tσt-exp⁡2dS0D-y0σ2Φ-D-y0+dS0tσt,where Φ(·) is the cumulative distribution function of standard normal distribution.


## 4.2. Use of Copulas for Multiple Features

### 4.2.1. Definition and Basic Properties

Sklar’s theorem lays the theoretical foundation for copulas and their applications [13].Theorem 1 (Sklar’s theorem [29]).
LetX=(X1,…,XN) be a random vector with marginal distributions F1(x1),…,FN(xN) and F be their joint distribution function. Then, there exists a copula function C such that(12)Fx1,…,xn,…,xN=CF1x1,…,Fnxn,…,FNxN.If F1(x1),…,FN(xN) are continuous, the copula function C is unique. Conversely, if F1(x1),…,FN(xN) are univariate distributions, the function F defined by (12) is the joint distribution function associated with margins F1(x1),…,FN(xN).Based on the theorem, a multivariate copula function can be defined as follows.Definition 2 (n-dimensional copula or n-copula [13]).
AnN-dimensional copula is a function C mapping from IN=0,1N to I with the following properties:(1)
Grounded: for everyu=(u1,…,uN) in IN, C(u)=0 if at least one coordinate of u is 0.(2)
Uniform marginals: if all coordinates ofu are ones except uk, then(13)Cu=C1,…,1,uk,1,…,1=uk.(3)
N-increasing: for each hyperrectangleB=∏i=1Nui,vi⊆[0,1]N, the C-volume of B is nonnegative:(14)∫BdCu,v=∑z∈×i=1Nui,vi-1NzCz≥0,whereN(z)=#k:zk=uk.
Furthermore, the PDF corresponding to joint distribution functionF(x1,…,xn,…,xN) can be calculated as(15)fx1,…,xn,…,xN=cF1x1,…,Fnxn,…,FNxN∏n=1Nfnxn,where fn(·) is the PDF of a marginal distribution Fn(·), and c(·) is the density of copula function C:(16)cu1,…,un,…,uN=∂NCu1,…,un,…,uN∂u1,…,∂un,…,∂uN∀u=u1,…,un,…,uN∈IN.


### 4.2.2. Selection of Copula Family

It is important to select an appropriate copula family for the given data, as the selection significantly influences the prediction performance. Table2 presents five popular multivariate copula families, where u1,…,un are n random variables.Table 2
Typical multivariate copulas.
Copulas C u 1 , u 2 , … , u n Parameters Gaussian copula Φ ρ Φ - 1 u 1 , Φ - 1 u 2 , … , Φ - 1 u n ρ ∈ - 1,1 Student’st-copula t ρ , v t v - 1 u 1 , t v - 1 u 2 , … , t v - 1 u n ρ ∈ - 1,1 k > 2 Clayton copula ∑ i = 1 n u i - δ - n + 1 - 1 / δ δ > 0 Frank copula - δ - 1 log ⁡ 1 +  ∏ i = 1 n exp ⁡ - δ u i - 1 exp ⁡ - δ - 1 n - 1 δ ∈ - ∞ , ∞ ∖ 0 δ > 0 for n ≥ 3 Joe copula 1 - 1 -  ∏ i = 1 n 1 - 1 - u i δ 1 / δ δ ∈ 1 , ∞To balance the statistical fitting and model complexity, the values of classic Akaike information criterion (AIC) are computed for the candidate copula families:(17)AIC=2k-2ln⁡L,where k is the number of estimated parameters in the model, and L is the value of the corresponding maximum likelihood. By comparing the AIC values, the family with the minimum value will be selected.


## 4.2.1. Definition and Basic Properties

Sklar’s theorem lays the theoretical foundation for copulas and their applications [13].Theorem 1 (Sklar’s theorem [29]).
LetX=(X1,…,XN) be a random vector with marginal distributions F1(x1),…,FN(xN) and F be their joint distribution function. Then, there exists a copula function C such that(12)Fx1,…,xn,…,xN=CF1x1,…,Fnxn,…,FNxN.If F1(x1),…,FN(xN) are continuous, the copula function C is unique. Conversely, if F1(x1),…,FN(xN) are univariate distributions, the function F defined by (12) is the joint distribution function associated with margins F1(x1),…,FN(xN).Based on the theorem, a multivariate copula function can be defined as follows.Definition 2 (n-dimensional copula or n-copula [13]).
AnN-dimensional copula is a function C mapping from IN=0,1N to I with the following properties:(1)
Grounded: for everyu=(u1,…,uN) in IN, C(u)=0 if at least one coordinate of u is 0.(2)
Uniform marginals: if all coordinates ofu are ones except uk, then(13)Cu=C1,…,1,uk,1,…,1=uk.(3)
N-increasing: for each hyperrectangleB=∏i=1Nui,vi⊆[0,1]N, the C-volume of B is nonnegative:(14)∫BdCu,v=∑z∈×i=1Nui,vi-1NzCz≥0,whereN(z)=#k:zk=uk.
Furthermore, the PDF corresponding to joint distribution functionF(x1,…,xn,…,xN) can be calculated as(15)fx1,…,xn,…,xN=cF1x1,…,Fnxn,…,FNxN∏n=1Nfnxn,where fn(·) is the PDF of a marginal distribution Fn(·), and c(·) is the density of copula function C:(16)cu1,…,un,…,uN=∂NCu1,…,un,…,uN∂u1,…,∂un,…,∂uN∀u=u1,…,un,…,uN∈IN.


## 4.2.2. Selection of Copula Family

It is important to select an appropriate copula family for the given data, as the selection significantly influences the prediction performance. Table2 presents five popular multivariate copula families, where u1,…,un are n random variables.Table 2
Typical multivariate copulas.
Copulas C u 1 , u 2 , … , u n Parameters Gaussian copula Φ ρ Φ - 1 u 1 , Φ - 1 u 2 , … , Φ - 1 u n ρ ∈ - 1,1 Student’st-copula t ρ , v t v - 1 u 1 , t v - 1 u 2 , … , t v - 1 u n ρ ∈ - 1,1 k > 2 Clayton copula ∑ i = 1 n u i - δ - n + 1 - 1 / δ δ > 0 Frank copula - δ - 1 log ⁡ 1 +  ∏ i = 1 n exp ⁡ - δ u i - 1 exp ⁡ - δ - 1 n - 1 δ ∈ - ∞ , ∞ ∖ 0 δ > 0 for n ≥ 3 Joe copula 1 - 1 -  ∏ i = 1 n 1 - 1 - u i δ 1 / δ δ ∈ 1 , ∞To balance the statistical fitting and model complexity, the values of classic Akaike information criterion (AIC) are computed for the candidate copula families:(17)AIC=2k-2ln⁡L,where k is the number of estimated parameters in the model, and L is the value of the corresponding maximum likelihood. By comparing the AIC values, the family with the minimum value will be selected.


## 4.3. Multivariate ADT Model

To analyze the multivariate CSADT data, additional assumptions are considered.(A4)
For a product with multiple degradation indicators, the dependency among the features can be characterized by a copula function.(A5)
The degradation measurements of all specimens are collected at the same time.(A6)
The specimen is considered to be failed if one of the features reaches its corresponding failure threshold for the first time.Suppose that a product hasP degradation features and each of the accelerated degradation processes is modeled by a drift Brownian motion. Yp(t), p=1,2,…,P, denotes the degradation process of the pth characteristic index at operating condition S0, and its corresponding failure threshold is Dp. Let the failure time of the pth characteristic index be Tp, and the lifetime of the product be T, so T=min(T1,…,Tp). Clearly, the system reliability under normal condition S0 can be expressed as(18)Rt=PT>t=PT1>t,…,TP>t=PY1t<D1,…,YPt<DP.Apparently, if theP degradation features are independent, (18) can be rewritten as(19)Rt=PY1t<D1×⋯×PYPt<DP=R1t×⋯×RPt,where Rp(t), p=1,2,…,P is the pth characteristic index’s reliability function under normal operating condition S0, as defined in (11). However, this is often not true in many engineering applications. To address the dependence of P degradation features, a copula function can be utilized. In particular, the product’s reliability at time t under normal operating condition S0 can be performed as(20)Rt=CR1t,R2t,…,RPt;θ,where θ is the parameter set of the copula.


## 4.4. Parameter Estimation

Consider a copula-based multivariate distribution for the random vector[Y1,…,YP], with joint PDF(21)fy;α1,…,αP,θ=cF1y1;α1,…,FPyP;αP;θ∏p=1Pfpyp;αp,where Fp(·;αp) and fp(·;αp), p=1,…,P, are the pth marginal CDF and PDF with parameter set αp, and c(·;θ) is the density of copula function C(·;θ) with parameter θ.For a sample of sizen, the full log-likelihood function can be expressed as(22)ln⁡Lα1,…,αP,θ=∑i=1nln⁡cF1y1i;α1,…,FPyPi;αP;θ︸dependence  structureLC+∑p=1P∑i=1nfpypi;αp︸marginals∑p=1PLp.The maximum likelihood estimates (MLE) of the parameters(α1,…,αp,θ) can be obtained by maximizing (22). However, as the dimension of parameters increases, it is difficult to get the optimal solution. Joe [30, 31] proposed the IFM method as a computationally attractive alternative to estimate the parameters in multivariate copula models. Technically, (22) is first decomposed into two parts. The first part, denoted by LC, is the contribution from the dependence structure in data, and the other part is the contributions from each margin, Lp, p=1,…,P. Then, The IFM method separately estimates the parameters of marginal distributions and the parameters of the copula in two stages. In the first stage, the estimates of the marginal distribution parameters are obtained from the log-likelihood Lp of each margin. In the second stage, the copula parameter estimates are computed by maximizing LC based on the marginal distribution parameters estimated in the first stage. In this paper, we use this method to estimate the parameters of the proposed multivariate ADT model.


### 4.4.1. Marginal Distribution Parameters Estimation

Suppose that in a CSADT withK stress levels the total sample size is n. There are nl specimens under stress level Sl. During the CSADT, all specimens are measured once in every Δt time interval and there are Ml inspections under Sl. Then, the observation of the pth characteristic index at time tlij is yp(tlij), p=1,…,P, l=1,…,K, i=1,…,nl, j=1,…,Ml, where tlij is the time of the jth measurement of the ith unit under the lth stress level. According to (8), the likelihood function of the pth characteristic indices is given by(23)Lp=∏l=1K∏i=1nl∏j=1Ml-112πσp2Δtexp⁡-yptlij+1-yptlij-μpl·Δt22σp2Δt,p=1,…,P.Considering the acceleration model given by (7), μpl can be expressed as μpl=dp(Sl)=exp⁡(ap+bpφ(Sl)), where ap and bp are unknown parameters of the pth characteristic index. So the model parameter set of the characteristic index is αp=(ap,bp,σp), p=1,…,P. Let Δyplij=yp(tlij+1)-yp(tlij); then (23) can be expressed as(24)ln⁡L=-12∑l=1K∑i=1nl∑j=1Ml-1ln⁡2πσp2Δt+Δyplij-exp⁡ap+bpφSl·Δt2σp2Δt,p=1,…,Pand the parameter set αp=(ap,bp,σp) of each marginal distribution can be estimated by maximizing the log-likelihood function.


### 4.4.2. Copula Parameter Estimation

By replacing the marginal parametersαp, p=1,…,P, with their first-stage estimateα^p, the copula parameters can be estimated by maximizing the copula likelihood contribution. Then, the copula parameter θ in (20) is given by(25)θ^=argmax∑p=1P∑i=1nln⁡cR1y1i;α^1,…,RPypi;α^P;θ.In addition, the AIC method mentioned in Section 4.2.2 is used to quantitatively select the best fitting model from candidate copula functions.


## 4.4.1. Marginal Distribution Parameters Estimation

Suppose that in a CSADT withK stress levels the total sample size is n. There are nl specimens under stress level Sl. During the CSADT, all specimens are measured once in every Δt time interval and there are Ml inspections under Sl. Then, the observation of the pth characteristic index at time tlij is yp(tlij), p=1,…,P, l=1,…,K, i=1,…,nl, j=1,…,Ml, where tlij is the time of the jth measurement of the ith unit under the lth stress level. According to (8), the likelihood function of the pth characteristic indices is given by(23)Lp=∏l=1K∏i=1nl∏j=1Ml-112πσp2Δtexp⁡-yptlij+1-yptlij-μpl·Δt22σp2Δt,p=1,…,P.Considering the acceleration model given by (7), μpl can be expressed as μpl=dp(Sl)=exp⁡(ap+bpφ(Sl)), where ap and bp are unknown parameters of the pth characteristic index. So the model parameter set of the characteristic index is αp=(ap,bp,σp), p=1,…,P. Let Δyplij=yp(tlij+1)-yp(tlij); then (23) can be expressed as(24)ln⁡L=-12∑l=1K∑i=1nl∑j=1Ml-1ln⁡2πσp2Δt+Δyplij-exp⁡ap+bpφSl·Δt2σp2Δt,p=1,…,Pand the parameter set αp=(ap,bp,σp) of each marginal distribution can be estimated by maximizing the log-likelihood function.


## 4.4.2. Copula Parameter Estimation

By replacing the marginal parametersαp, p=1,…,P, with their first-stage estimateα^p, the copula parameters can be estimated by maximizing the copula likelihood contribution. Then, the copula parameter θ in (20) is given by(25)θ^=argmax∑p=1P∑i=1nln⁡cR1y1i;α^1,…,RPypi;α^P;θ.In addition, the AIC method mentioned in Section 4.2.2 is used to quantitatively select the best fitting model from candidate copula functions.


## 5. Case Study

A case study is provided in this section to demonstrate the use of the proposed framework in real-world industrial applications.


### 5.1. Description of the ADT Experiment

DC motor converts DC power into mechanical energy. It is a key part of robots, electric vehicle, numerical control machines and household appliances, and so forth. The failure of motors will significantly decrease the availability and safety of such machines and even lead to big economic losses.To assess the life and reliability of PZB131 brushed DC motors, a dedicated CSADT platform was designed as shown in Figure2. The platform is composed of DC power, electric motors, vibration sensors, charge amplifier, data acquisition card, and computer. The key failure mechanism of this type of DC motors is the wear-out between brush and commutator, which is mainly related to the rotating speed and load. To speed up the degradation process, the input voltage was employed as the accelerated stress that controls the rotating speed of motors under a nonload condition.Figure 2
Overview of the motor vibration-based ADT platform.The rated voltage of the motor is 2.0 V, and the accelerated stress levels are 3.5 V, 4.5 V, and 6.0 V. The sample size under each stress level was 3. The vibration signals on the surface of motor were collected every hour. The sampling time is 2 seconds and the total number of sampling points is 3000.


### 5.2. Feature Extraction and Selection

First, denoising of raw vibration signals was performed using the Daubechies (db8) wavelet method based on the principle of Stein’s unbiased risk with soft-thresholding [24]. Ten features of each motor were exacted using the signal processing methods addressed in Section 3.1. Figure 3 shows all the candidate features collected from one of the motors.Figure 3
Feature set of a motor.Afterwards, feature selection was implemented using the proposed feature selecting procedure given in Section3.2, and the most suitable features were selected from the ten candidate features. Because prognostic features showing significant monotonic trends are more desired [7], different weights, 0.3, 0.6, and 0.1, were given to trendability, monotonicity, and robustness, respectively. The feature selection results are presented in Table 3.Table 3
Feature evaluation results.
No. Feature Stress T M R Fitness Mean (f) Rank 1 Mean 3.5 V 0.2491 0.5556 0.9500 0.5425 0.5102 8 4.5 V 0.1616 0.3333 0.9839 0.3469 6.0 V 0.7032 0.5556 0.9682 0.6411 2 Variance 3.5 V 0.6042 0.5556 0.8549 0.6001 0.5350 7 4.5 V 0.4839 0.3333 0.7852 0.4237 6.0 V 0.5634 0.5556 0.7892 0.5813 3 Crest factor 3.5 V 0.8614 0.5556 0.9808 0.6898 0.7509 3 4.5 V 0.7359 1 0.9263 0.9134 6.0 V 0.7430 0.5556 0.9319 0.6494 4 Kurtosis factor 3.5 V 0.3349 0.5556 0.9438 0.5282 0.4236 10 4.5 V 0.3380 0.3333 0.7947 0.3809 6.0 V 0.3179 0.3333 0.6620 0.3616 5 Waveform factor 3.5 V 0.5856 0.3333 0.9987 0.4756 0.6896 4 4.5 V 0.7208 0.7778 0.9962 0.7825 6.0 V 0.8182 0.7778 0.9858 0.8107 6 Clearance factor 3.5 V 0.8603 1 0.9436 0.9525 0.8961 1 4.5 V 0.8570 1 0.9460 0.9517 6.0 V 0.7649 0.7778 0.8804 0.7842 7 Frequency center 3.5 V 0.6258 0.5556 0.8631 0.6074 0.6260 5 4.5 V 0.6474 0.5556 0.7387 0.6014 6.0 V 0.8213 0.5556 0.8940 0.6691 8 Frequency variance 3.5 V 0.6124 0.5556 0.8679 0.6039 0.6194 6 4.5 V 0.6907 0.5556 0.7590 0.6164 6.0 V 0.7157 0.5556 0.8978 0.6378 9 HHT marginal energy (HME) 3.5 V 0.5675 0.5556 0.7430 0.5779 0.5057 9 4.5 V 0.4637 0.5556 0.8415 0.5566 6.0 V 0.3360 0.3333 0.8192 0.3827 10 HHT marginal energy spectrum entropy (HMESE) 3.5 V 0.7513 0.7778 0.9545 0.7875 0.8569 2 4.5 V 0.8404 1 0.9643 0.9486 6.0 V 0.9094 0.7778 0.9518 0.8347Because the clearance factor, crest factor, and HHT marginal energy spectrum entropy (HMESE) fitness were all greater than 0.7, they were selected as the degradation indicators for the motors. The degradation processes of three selected optimal features are illustrated in Figures4–6.Figure 4
The CSADT data of Hilbert marginal energy spectrum entropy.Figure 5
The CSADT data of clearance factor.Figure 6
The CSADT data of crest factor.


### 5.3. Reliability Modeling

Defining the failure threshold for each feature is the first step for degradation-based reliability modeling. To this end, two motors were tested until failure, and the failure thresholds were determined by analyzing the clearance factor, Hilbert marginal energy spectrum entropy, and crest factor of the motors upon failure. As can be seen in Figure7, the relative failure threshold D-y0 of Hilbert marginal energy spectrum entropy, clearance factor, and crest factor is 0.5, 3, and 3, respectively.Figure 7
Determination of the failure threshold.Then, the univariate accelerated degradation models for the three features were established, respectively, using the Brownian motion models addressed in Section4.1. The model parameters were estimated using the method given in Section 4.4.1. Figure 8 shows the reliability functions of the three features calculated by (11). One can see that the reliability curve of crest factor tends to be higher than those of the other two features.Figure 8
Reliability of three features under normal operating condition.Next, the multivariate ADT model was developed using the copula function that describes the dependence among three features, and the IFM method was used to estimate the copula parameters. Figure9 shows the system reliability curves of the motor using different copulas. To determine the best copula function, the AIC values of the five candidate copula families were obtained and compared. From Table 4, one can see that the AIC value of Frank copula function is the minimum, so it was adopted as the best choice among the candidates.Table 4
Goodness-of-fit for the candidate copulas.
Copula family AIC value Ranking Gaussian copula 42961 5 t-copula −2278 2 Clayton copula −1179 3 Frank copula −17136 1 Joe copula −899 4Figure 9
Comparison of different copula functions.Figure10 shows the marginal reliability and system reliability functions under both s-dependent and s-independent assumptions. One can see that the reliability function by assuming s-independent degradation features tends to be much lower than that based on the s-dependent assumption. Clearly, the dependence among the multiple features cannot be ignored.Figure 10
Comparison of the system reliability.


## 5.1. Description of the ADT Experiment

DC motor converts DC power into mechanical energy. It is a key part of robots, electric vehicle, numerical control machines and household appliances, and so forth. The failure of motors will significantly decrease the availability and safety of such machines and even lead to big economic losses.To assess the life and reliability of PZB131 brushed DC motors, a dedicated CSADT platform was designed as shown in Figure2. The platform is composed of DC power, electric motors, vibration sensors, charge amplifier, data acquisition card, and computer. The key failure mechanism of this type of DC motors is the wear-out between brush and commutator, which is mainly related to the rotating speed and load. To speed up the degradation process, the input voltage was employed as the accelerated stress that controls the rotating speed of motors under a nonload condition.Figure 2
Overview of the motor vibration-based ADT platform.The rated voltage of the motor is 2.0 V, and the accelerated stress levels are 3.5 V, 4.5 V, and 6.0 V. The sample size under each stress level was 3. The vibration signals on the surface of motor were collected every hour. The sampling time is 2 seconds and the total number of sampling points is 3000.


## 5.2. Feature Extraction and Selection

First, denoising of raw vibration signals was performed using the Daubechies (db8) wavelet method based on the principle of Stein’s unbiased risk with soft-thresholding [24]. Ten features of each motor were exacted using the signal processing methods addressed in Section 3.1. Figure 3 shows all the candidate features collected from one of the motors.Figure 3
Feature set of a motor.Afterwards, feature selection was implemented using the proposed feature selecting procedure given in Section3.2, and the most suitable features were selected from the ten candidate features. Because prognostic features showing significant monotonic trends are more desired [7], different weights, 0.3, 0.6, and 0.1, were given to trendability, monotonicity, and robustness, respectively. The feature selection results are presented in Table 3.Table 3
Feature evaluation results.
No. Feature Stress T M R Fitness Mean (f) Rank 1 Mean 3.5 V 0.2491 0.5556 0.9500 0.5425 0.5102 8 4.5 V 0.1616 0.3333 0.9839 0.3469 6.0 V 0.7032 0.5556 0.9682 0.6411 2 Variance 3.5 V 0.6042 0.5556 0.8549 0.6001 0.5350 7 4.5 V 0.4839 0.3333 0.7852 0.4237 6.0 V 0.5634 0.5556 0.7892 0.5813 3 Crest factor 3.5 V 0.8614 0.5556 0.9808 0.6898 0.7509 3 4.5 V 0.7359 1 0.9263 0.9134 6.0 V 0.7430 0.5556 0.9319 0.6494 4 Kurtosis factor 3.5 V 0.3349 0.5556 0.9438 0.5282 0.4236 10 4.5 V 0.3380 0.3333 0.7947 0.3809 6.0 V 0.3179 0.3333 0.6620 0.3616 5 Waveform factor 3.5 V 0.5856 0.3333 0.9987 0.4756 0.6896 4 4.5 V 0.7208 0.7778 0.9962 0.7825 6.0 V 0.8182 0.7778 0.9858 0.8107 6 Clearance factor 3.5 V 0.8603 1 0.9436 0.9525 0.8961 1 4.5 V 0.8570 1 0.9460 0.9517 6.0 V 0.7649 0.7778 0.8804 0.7842 7 Frequency center 3.5 V 0.6258 0.5556 0.8631 0.6074 0.6260 5 4.5 V 0.6474 0.5556 0.7387 0.6014 6.0 V 0.8213 0.5556 0.8940 0.6691 8 Frequency variance 3.5 V 0.6124 0.5556 0.8679 0.6039 0.6194 6 4.5 V 0.6907 0.5556 0.7590 0.6164 6.0 V 0.7157 0.5556 0.8978 0.6378 9 HHT marginal energy (HME) 3.5 V 0.5675 0.5556 0.7430 0.5779 0.5057 9 4.5 V 0.4637 0.5556 0.8415 0.5566 6.0 V 0.3360 0.3333 0.8192 0.3827 10 HHT marginal energy spectrum entropy (HMESE) 3.5 V 0.7513 0.7778 0.9545 0.7875 0.8569 2 4.5 V 0.8404 1 0.9643 0.9486 6.0 V 0.9094 0.7778 0.9518 0.8347Because the clearance factor, crest factor, and HHT marginal energy spectrum entropy (HMESE) fitness were all greater than 0.7, they were selected as the degradation indicators for the motors. The degradation processes of three selected optimal features are illustrated in Figures4–6.Figure 4
The CSADT data of Hilbert marginal energy spectrum entropy.Figure 5
The CSADT data of clearance factor.Figure 6
The CSADT data of crest factor.


## 5.3. Reliability Modeling

Defining the failure threshold for each feature is the first step for degradation-based reliability modeling. To this end, two motors were tested until failure, and the failure thresholds were determined by analyzing the clearance factor, Hilbert marginal energy spectrum entropy, and crest factor of the motors upon failure. As can be seen in Figure7, the relative failure threshold D-y0 of Hilbert marginal energy spectrum entropy, clearance factor, and crest factor is 0.5, 3, and 3, respectively.Figure 7
Determination of the failure threshold.Then, the univariate accelerated degradation models for the three features were established, respectively, using the Brownian motion models addressed in Section4.1. The model parameters were estimated using the method given in Section 4.4.1. Figure 8 shows the reliability functions of the three features calculated by (11). One can see that the reliability curve of crest factor tends to be higher than those of the other two features.Figure 8
Reliability of three features under normal operating condition.Next, the multivariate ADT model was developed using the copula function that describes the dependence among three features, and the IFM method was used to estimate the copula parameters. Figure9 shows the system reliability curves of the motor using different copulas. To determine the best copula function, the AIC values of the five candidate copula families were obtained and compared. From Table 4, one can see that the AIC value of Frank copula function is the minimum, so it was adopted as the best choice among the candidates.Table 4
Goodness-of-fit for the candidate copulas.
Copula family AIC value Ranking Gaussian copula 42961 5 t-copula −2278 2 Clayton copula −1179 3 Frank copula −17136 1 Joe copula −899 4Figure 9
Comparison of different copula functions.Figure10 shows the marginal reliability and system reliability functions under both s-dependent and s-independent assumptions. One can see that the reliability function by assuming s-independent degradation features tends to be much lower than that based on the s-dependent assumption. Clearly, the dependence among the multiple features cannot be ignored.Figure 10
Comparison of the system reliability.


## 6. Conclusions

This study is focused on the reliability modeling of vibration-based ADT data. The proposed framework for data analysis focuses on feature selection and multivariate accelerated degradation modeling. A set of feature selection criteria were used to select the most suitable features after extracting features from raw vibration signals. Then, a multivariate s-dependent ADT model was presented, where Brownian motion was applied to establish accelerated degradation model for each feature and a copula function was employed to characterize the dependence among multiple degradation processes. A case study on motor vibration-based CSADT data was conducted to demonstrate the usefulness and validity of the proposed model and statistical method. The results show that the proposed framework can provide a powerful tool for modeling the reliability of products exhibiting multiple dependent degradation features.



---

*Source: 2315916-2016-12-07.xml*