Reliability estimation from two types of accelerated testing data considering measurement error

Reliability Engineering and System Safety 193 (2020) 106610

Contents lists available at ScienceDirect

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

Reliability estimation from two types of accelerated testing data considering measurement error

T

Ma Zhonghaia, Wang Shaopinga, Cesar Ruizb, Zhang Chaoa, , Haitao Liaob, Edward Pohlb ⁎

a b

School of Automation Science and Electrical Engineering, Beihang University, Beijing Haidian district, Xueyuan Road No. 37, Beijing, 100191, China Department of Industrial Engineering, University of Arkansas, Fayetteville, AR, USA

ARTICLE INFO

ABSTRACT

Keywords: Accelerated life testing (ALT) Accelerated degradation testing (ADT) Measurement error Inverse Gaussian (IG) process Reliability estimation Expectation-maximization (EM)

Reliability testing is an indispensable tool for evaluating the lifetime of a product. However, for a highly reliable product, it is quite common that a large proportion of test units will be censored in a regular life test or even in accelerated life testing (ALT) when the total testing time is too short. As an alternative, accelerated degradation testing (ADT) can be conducted to collect degradation data of a highly reliable product under accelerated conditions. For a reliability practitioner, it will be very valuable to use both ALT and ADT data for reliability estimation. In practice, degradation data are often contaminated by measurement error, which may affect the accuracy of reliability estimation. Therefore, a statistical procedure is needed when using both ALT data and ADT data with measurement error for evaluating the reliability of a highly reliable product. In this paper, an Inverse Gaussian (IG) process is used to model the degradation process of a product considering measurement error. To incorporate the two types of accelerated testing data, a new expectation-maximization (EM) algorithm is developed to estimate the model parameters by taking advantage of the parameter structure. A simulation study and a case study on a hydraulic piston pump are presented to illustrate the practical value of the proposed method in improving the accuracy of reliability estimation for a highly reliable product.

1. Introduction Accelerated testing is usually adopted to speed up failures of highly reliable products for reliability analysis. In such tests, products are exposed to harsher-than-normal conditions (e.g., load, strain, temperature, voltage, vibration, pressure) in an effort to collect failure time data in a short amount of time without changing the failure mechanisms [1,2]. For example, the lifetimes of some key components in an aircraft hydraulic system, such as the hydraulic piston pump, hydraulic actuator and landing gear, are expected to be in the order of thousands of hours. It is often difficult, if not impossible, to observe failures in a short period of time under normal operating conditions. In order to reduce the product development time, it is necessary to collect failure time data in a quicker way using some form of accelerated testing. Accelerated testing can be classified into accelerated life testing (ALT) and accelerated degradation testing (ADT) [3–6]. Elsayed et al. [7] provided an overview of different reliability testing methods along with prediction models and briefly discussed ADT and ALT plans. The most important factors that affect the accuracy of reliability prediction are the underlying failure time distribution, the ways stress factors are applied and the understanding of underlying life-stress relationship. In

⁎

practice, four types of stress loadings have been widely used: constantstress, step-stress, progressive-stress and cyclic-stress [7,8], where a constant-stress test is the most commonly applied. The life-stress relationship also plays an important role in modeling accelerated testing data. Such a relationship can be modeled using a statistical approach or a physics-based model such as the Arrhenius for temperature, Eyring for temperature and humidity, Blattau for temperature cycling and the exponential-power-law (EPL) model for a mechanical stress [9,10]. The lifetime of a highly reliable product is usually in the range of thousands of hours. An example of such a product is a hydraulic piston pump which is a key component of hydraulic systems. Several researchers have conducted studies on the failure mechanisms, prognostics health management (PHM) and ALT of such pumps. Ma et al. studied the typical failure modes and methods to conduct ALT of hydraulic piston pumps [11]. Wang et al. proposed a statistical methodology for analyzing ALT data of hydraulic pumps [12]. However, it is difficult to obtain failure time data in a short testing time even using ALT which makes it difficult to accurately estimate the reliability of interest. An alternative to overcome this difficulty is to conduct degradation tests. As an example, Peng et al. proposed a method to calibrate the initial model parameter estimates based on test data once field

Corresponding author. E-mail addresses: [email protected] (Z. Ma), [email protected] (C. Zhang).

https://doi.org/10.1016/j.ress.2019.106610 Received 4 April 2019; Received in revised form 16 July 2019; Accepted 11 August 2019 Available online 17 August 2019 0951-8320/ © 2019 Elsevier Ltd. All rights reserved.

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

Notation

NL , NL* L(•), E[•] k M0 M1

number of failed units tested under ALT number of censored units tested under ALT (•) likelihood function and log-likelihood function expected value number of cycles method that only considers ALT data method that only considers ADT data without measurement error M2 method that only considers ADT data with measurement error M3 method that considers ALT data and ADT data without considering measurement error M4 method that considers ALT data and ADT data with measurement error B10 time at which 10 percent of units have failed (i.e., 10 percent quantile) Pa, Va accelerated discharge pressure and rotating speed of hydraulic piston pump P0, V0 discharge pressure and rotating speed of pump under normal operation ν shape parameter of hydraulic piston pump time transformation ω,α scale parameters of hydraulic piston pump time transformation AF acceleration factor

N number of units tested Y(t) true value of the degradation signal at time t G(•|Λ(tq), μ, λ) cumulative distribution function (CDF) of the Inverse Gaussian (IG) process μ, λ mean and diffusion parameters of the IG process S stress level Λ(t) a monotonically increasing function of time with (0) = 0 ɛ measurement error term σ2 Variance of Z(t) performance degradation measurement at time t Δzi, j observed degradation increment of unit i in time interval [tj – 1; tj] ΔZ vector of all degradation increments for all units. Σ covariance matrix of the degradation measurements. φ(•) probability density function of the truncated multivariate normal distribution Φ(•) CDF of the standard normal distribution D failure threshold of the degradation process TD failure time, first passage time of the level of degradation D tq failure time of unit q tr* censoring time ND number of units tested under ADT

data with measurement error. Section 4 provides the simulation study to illustrate the accuracy of the proposed method compared to the alternatives of using only ALT or ADT data. Then, the case study on a hydraulic piston pump is provided to illustrate the use of the proposed method in practice. Finally, a summary and conclusions are given in Section 5.

degradation data becomes available [13]. In addition, Wang et al. assess the reliability based on the Wiener degradation process while considering multiple sources of information (prior degradation information, life data information, and current degradation data) [14]. Compared to ALT, ADT is more suitable for highly reliable products as it can further speed up reliability testing without waiting for actual failures to occur. Moreover, degradation data may be abundant for each test unit, which may contain more reliability information than failure time data [3]. Product degradation is often modeled by a stochastic process {Y(t); t ≥ 0} to account for the cumulative damage along with inherent randomness of the process [15–17]. The most commonly used stochastic process models for degradation data analysis are the Wiener process, Gamma process and Inverse Gaussian (IG) process [18–20]. For a degradation process that is irreversible in nature (e.g., mechanical wear, leakage), the IG process is a natural choice for modeling the corresponding degradation data. Ye et al. [21] and Méndez-González et al. [22] adopt the Inverse Gaussian (IG) process to analyze constantstress ADT, and a random-effects IG process model was considered in designing the optimal ADT plans. However, it is inevitable that some measurement errors may be introduced during degradation data collection. Such measurement errors are usually caused by sensor inaccuracy, random environments, or human error. Such errors may significantly influence the prediction accuracy of a degradation model. Li et al. [23], Pan et al. [24] and Ye et al. [25] proposed a Wiener process model for ADT analysis with measurement error. In this paper, both ADT data and ALT data are considered in order to improve our capability in reliability estimation. A statistical method is proposed to integrate these two types of data. In addition, measurement error for the ADT data is considered to improve the accuracy of parameter estimation. The IG process is used to describe a monotonically increasing and irreversible degradation process that is affected by complex working conditions. An expectation-maximization (EM) algorithm is developed for estimating the model parameters. The proposed method is illustrated via a simulation study as well a case study. The remainder of this paper is organized as follows. Section 2 presents the IG process and the models for ALT data and ADT data with measurement error. The EM algorithm is presented in Section 3. It is used to estimate the model parameters based on both ALT data and ADT

2. Inverse Gaussian process model for accelerated testing data Several components in hydraulic systems experience irreversible damage as time passes, such as fatigue and wear of hydraulic piston pumps [26]. Ultimately, such damage will lead to a failure when the system performance reaches a failure threshold. In contrast to the Wiener process, whose degradation path is not necessarily monotone, the IG process is more suitable for modelling the degradation of a mechanical component because it is a monotonically increasing function. 2.1. Inverse Gaussian process model for ADT data with measurement error The degradation process Y(t) of a product follows an IG process if it has the following properties:

• Y(0) ≡ 0; • Y(t) has •

independent increments, that is, Y (t2) Y (t1) and Y (s2 ) Y (s1) are independent for any two non-overlapping time intervals, ∀t2 > t1 ≥ s2 > s1; Y (t ) Y (s ) follows an IG distribution IG(μΔΛ, λΔΛ2) with = (t ) (s ), ∀t > s ≥ 0, where μ and λ are the mean and diffusion parameters, respectively, and Λ(t) is a monotonically time transformation function with (0) = 0 .

In the ADT analysis, degradation-stress relationships are usually utilized to describe the correlation between the accelerated data and accelerated stress(s). In this paper, a constant-stress ADT is considered. In the constant-stress ADT, let S0 be the normal operation stress level, and S1 < S2 < ,…, < SL be the stress levels used in ADT where L is the highest stress level used in the test. Consider N test units whose performance measures degrade over time under stress Sl. The performance 2

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

degradation of unit i at each point at time 0 = ti,0 < ti,1 < ti, m is denoted by Y(ti,j|Sl), and the degradation increment y (ti, j |Sl ) = y (ti, j |Sl ) y (ti, j 1 |Sl ) in time interval [t j 1, t j] follows an IG distribution [20]. The increment time transformation is (ti, j |Sl ) (ti, j 1 |Sl ) which is related to the stress level. In this i, j = paper, define Y(ti,j|Sl) = Y(ti,j) and (ti, j |Sl ) = (ti, j ) for simplicity. By defining yi . j = y (ti, j ) , the probability density function (PDF) of degradation increment can be expressed as:

fIG ( yi, j µ , ,

(

i, j) =

i, j

)2

( yi, j

exp

2 ( yi, j )3

when tq is large), Y(tq) is approximately normally distributed with mean μΛ(tq) and variance μ3Λ(tq)/λ [27, 28]. Then, the PDF and CDF of failure time can be approximately expressed as:

2µ2 yi, j

(1)

f ( R IG

Zi

z i, m

=

i

µ, ) (

zi,1

i) d

j =1

i) d

i,1

d

m

(2 ) 2

1 2 exp

P(

(

i

1 ( 2

i)

1(

i)

)

for

Z, t , t *) = LA DT (µ , ,

Zi)

(8)

Z)*LA LT (µ ,

f ( Zi µ , , ) i=1

NL*

R (tr* µ , )

f (tq µ, ) q=1

r=1

(9)

Z, t , t *) = ln L (µ , ,

(µ , ,

=

ADT (µ ,

Z, t , t *) ,

Z) +

A LT (µ ,

t) +

Note that, the log-likelihood function for ADT data without considering measurement error is: ADT (µ ,

ND i=1

y) =

ln fIG ( yi µ , ) =

ND i=1

m j=1

1 2

ln + ln(

ND i=1

m j=1

(ti, j ))

ln fIG ( yi, j µ , ) (D

µ

(11) where yi = { yi,1 , ..., yi, m } . On the other hand, the log-likelihood function considering measurement error can be described as:

2

=

(4)

ADT (µ ,

,

ND

Z) =

ln f ( Zi µ , , )

(12)

i=1

2.2. Inverse Gaussian process model for ALT data

where f(ΔZi|μ, λ, σ) is given by Eq. (2). Moreover, for the ALT data, the corresponding log-likelihood function is given by:

Due to the monotonicity of the IG process, once the degradation process of unit q reaches failure threshold D at tq, the unit fails. According to Eq. (1), the cumulative distribution function (CDF) of failure time TD can be expressed as:

P(TD < tq) = P(Y (tq) > D) = 1 =

D

(tq)

(ti, j ))2

2µ2 yi, j

2

j, j

t *)

ALT * (µ ,

(10)

where P(Δɛi ≤ ΔZi) is the fraction after truncation, and the covariance matrix Σ is a positive definite tridiagonal matrix with:

j=j =1 2 j=j >1 2 j j =1 0, otherwise

t *)

t )*LALT * (µ ,

NL

=

(3)

0, otherwise

.

(tr*)/

ND

i, m

Zi

i

µ (tr*)

µ3

Suppose ADT data and ALT data with right censoring are available. In ADT, ND samples are tested under an accelerated stress level. In ALT, NL failure times tq (q = 1,2,…, NL) are obtained and NL* units are right censored at time tr* (r = 1,2,…, NL*). The point estimates of the model parameters can be obtained through direct maximization of the loglikelihood function for the data from different sources. The likelihood function L(μ, λ, σ|ΔZ, t, t*) and log-likelihood function (µ, , | Z , t ,t *) of the data can be expressed as:

where Zi = { z i,1, ..., z i, m}T is the vector of degradation increments of unit i, Z = { Z1, ..., ZN }T contains the degradation data of all N test T is the vector of increments in the units, and i = { i,1, ..., i, m} measurement error for the ith unit. Note that Δɛi follows a truncated multivariate normal distribution, Δɛi ∼ N(0, Σ), with a joint PDF: i) =

(7)

D

F (tr* µ , ) =

(2)

(

,

(tq)/

3.1. Maximum likelihood estimation

L (µ , , µ, ) (

(6)

3. Model parameter estimation

i

i, j

,

2µ2D

µ (tq)

µ3

R (tr* µ , ) = 1

m

fIG ( z i, j

D

µ (tq))2

(D

exp

respectively. Then, the reliability function evaluated at a right censoring time tr* can be expressed as:

It is quite common that measurement error is unavoidable when measuring the degradation process with a sensor. To incorporate the impact of measurement error in degradation-based reliability analysis, a stochastic process model explaining both the evolution of actual performance degradation and measurement error must be developed. Let Z be the performance degradation measurement, which is composed of the true degradation measure Y and an error term ɛ that follows the Normal distribution with mean 0 and unknown variance σ2, that is Z (t ) = Y (t ) + . Because the real degradation measure satisfies Y(t) ≥ 0, it should be noted that a Truncated Normal distribution with an upper bound of z must be considered for the error term, that is < z . Under this model, let Zi,j be the jth degradation measurement of unit i. Then, the observed degradation increment of unit i in time interval [tj – 1; tj] can be expressed as: z i, j = yi, j + i, j , where i , j = i, j i, j 1 (j = 2,3,…,m). Note that Δɛ is not ini,1 = i,1 and dependent since there exists a dependence constraint among Δɛi, j. For the ith unit, the PDF of observed degradation increment is: f ( Zi µ , , ) =

2 D3

F (tq µ , ) = 1

2 i, j)

µ

(tq ) 2

fIG (tq µ , ) =

D µ

G (D exp

2

NL ALT (µ ,

q=1

ln r=1

(tq ) D

1 ln + ln (tq) 2

NL*

+

(tq), µ , ) µ

t ,t *) =

D (tq ) + µ

D µ3

(D

µ (tq)) 2 2µ2 D

µ (tr*) (tr*)/

(13)

It is worth pointing out that estimating the model parameters by considering measurement error is computationally expensive. To overcome the challenge, an Expectation-Maximization (EM) algorithm is developed in this paper. The EM algorithm iteratively searches for the maximum likelihood estimates by a two-step process:

(5) where G(•|Λ(tq), μ, λ) is the CDF of IG(μΛ(tq), λΛ(tq)2), and Φ(•) is the CDF of the standard normal distribution. When λΛ(tq) is large (e.g., 3

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

(1) The expectation (E) step - creates a function for the expectation of the log-likelihood evaluated using the current estimates of the model parameters. (2) The maximization (M) step - updates the parameters estimated by maximizing the expected log-likelihood function from the E-step.

parameters of the proposed model. 4. Illustrative examples In this section, a simulation study is conducted to validate the proposed parameter estimation methodology. We consider four alternative models (M0–M3) by using only one source of data or ignoring measurement errors. Table 1 shows the assumptions and log-likelihood function of each method. The data sources and log-likelihood functions of these methods are shown in Table 1. Appendix A details the statistical procedures used for estimating the parameters of all the methods.

The two-step process is iteratively carried out until the log-likelihood value is not significantly improved compared to that of the previous iteration (i.e., the parameter estimates converge). Note that for the E-step, the initial guesses µ = µ 0 , = 0 and = 0 for the model parameters can be selected based on the data, then y^i, j = E [ yi, j ] can be found by using σ0, where E [ yi, j ] = z i, j E [ i, j] [29]. For a truncated multivariate normal distribution with upper bound ΔZ, the expected value of Δɛi for unit i can be obtained from the moment generating function as [30]:

4.1. Simulation study By assuming the underlying degradation process follows an IG process, both ADT data with measurement error and ALT data are generated. In generating ADT data, the degradation with time transformation (t ) = t 0.7 is applied, so that the IG model has a mean parameter µ = 4 and a diffusion parameter = 10 . In addition, it is assumed that the measurement error has mean 0 and variance 2 = 0.01 . A test unit is considered failed if its degradation process reaches a failure threshold of D = 70. Forty sets of ADT and ALT data are generated. In each data set, 5 units are tested under ADT and inspected every hour up to 50 h, 10 units are tested under ALT with censoring time of t*. Note that the sample size of ADT data is ND = 5, and ALT data includes NL failure times and NL* right censoring times where the values of NL and NL* may be different for those 40 data sets. As an example, Fig. 2 shows the simulated ADT samples with measurement error for one such replication, Fig. 3 shows the measurement errors of ADT data, and the simulated ALT data are given in Table 2. The average values for parameter estimates from the 40 data sets along with their errors compared to the true values are shown in Table 3.

m

E[

i

i

< zi ] =

(Fj ( zj ))

(14)

j=1

where

Fj (x ) z1

=

zj 1

zj + 1

zm

...

(x1,. ..,xj 1, x , x j + 1, ...x m) d xm. ..dxj + 1 (15)

dxj 1. ..dx1,

in which, φ(•) is the PDF of a truncated multivariate normal distribution. For the M-step, after obtaining the expected value of Δɛi, the parameters μ and λ can be estimated using Eq. (11) instead of Eq. (12) with known parameter σ0. In this paper, Eq. (10) was numerically maximized by using ADT (µ , | y) from Eq. (11) with y^ij , and coded using R software. Then, the estimated parameters µ^ and ^ are used for the next iteration.

However, after obtaining the values of µ^ and ^ , parameter σ needs to be calculated with the updated values. To remedy this, we separate the model parameters involved in Eq. (12) to find ^ . Because the loglikelihood function in Eq. (12) involves a multidimensional integration of the form given by Eq. (2), a Monte Carlo method is utilized, whose convergence rate (i.e., O (N 1/2 ) ) is independent of the order of in^ and Δɛi as: tegration [31]. Define the following function of y i

g ( zi ,

i)

= fIG ( zi

i

^, µ^ )

(16)

where parameters µ^ and ^ are updated in the first step. Eq. (2) can be described as an expectation and included in the log-likelihood function as follows: ADT (

µ^ , ^) =

ND

ln i=1 ND

=

ln i=1 ND

=

f ( R IG

R

Zi

i

g ( zi ,

i)

ln(E [g ( zi ,

^, µ^ ) (

(

i) d

i) d

i

i

i )])

(17)

i=1

where the integration region R is defined by Δɛi ≤ ΔZi, i = 1,2,…ND, and E[g(Δzi, Δɛi)] is estimated using Ns sets of sample points drawn from a truncated multivariate normal distribution with Δɛi ∼ N(0, Σ). Given the parameter estimates µ^ and ^ obtained in the first part, the ( |µ^ , ^) can be expressed as follows: log-likelihood function ADT

ADT (

=

µ^ , ^)

ND i=1

ln

(

ND i=1

ln

=

ND i=1

ln

1 Ns

Ns s=1

(

(

1 Ns

Ns g( s=1

1 Ns

Ns f ( s = 1 IG

m f ( j = 1 IG

z i, j

zi ,

i

(s ) )

zi

)

i i, j

(s ) ;

(s );

µ^ , ^)

µ^ , ^)

)

) (18) Fig. 1. The flowchart of EM algorithem for parameter estimation.

Fig. 1 shows the flowchart of the EM algorithm for estimating 4

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

Table 1 Alternative methods considering different sources of data. Method

ADT data

Measurement error (ɛ)

ALT data

M0

×

×

✓

M1

✓

M2

×

✓

M3

✓

✓

M4

× ×

×

✓

✓

✓

✓

Fig. 3. Simulated measurement errors in ADT data. Table 2 Simulated ALT data (10 units). Time

No

Time

1 2 3 4 5

46 64 54 69 51

6 7 8 9 10

66 61 67 47 70*(censoring)

⁎

M0

=

ALT (µ ,

|t ,t *)

M1

=

ADT (µ ,

| y)

M2

=

ADT (µ ,

, | Z)

M3

=

ADT (µ ,

| y) +

M4

=

ADT (µ ,

, | Z) +

ALT (µ ,

|t ,t *)

ALT (µ ,

|t ,t *)

error of λ is reduced, which benefits from the fair amount of degradation measurements. It is important to indicate that method M2 improves the estimation of λ compared to method M1 due to separating the variability of the observed degradation increments into the variabilities of underlying degradation process and measurement error. Method M3 significantly improves the estimation of μ without significantly affecting the estimation of λ. The best method (i.e., M4) takes full advantages of both types of data and obtains the most accurate parameter estimates. This is one of the main contributions of this paper. As shown in Fig. 4, the B10 life estimated by methods M0 and M1 are both larger than the true value, and the reliability function estimated using method M0 is larger than the true function in the product's early age and much less in the product's later age. The reliability function estimated by method M1 is always large than the true function. In practice, a maintenance decision made based on such reliability estimation results may have a significantly negative economic influence. It is obvious that only method M4 adequately captures the real values throughout the system lifetime. The main reason is that method M4 uses ALT data as well as ADT data and considers the measurement error to increase the estimation accuracy of the model parameters, and especially the diffusion parameter. Moreover, due to the use of both ALT and ADT data, method M3, which ignores the measurement error, gives a more accurate estimate of μ compared to method M2 that only considers ADT data with measurement error. One can see from Table 3 that if the measurement error is not considered, the diffusion λ will have a greater error regardless of whether ADT data or both data types are used. In addition, by considering ALT data, the estimate of μ can be corrected to reduce the estimation error. As shown in Fig. 5, by combining two types of data and introducing the measurement error, the reliability function can be estimated more accurately. Figs. 6 and 7 show the boxplots of relative errors of parameter estimations using methods M1–M4. Fig. 6 clearly shows that methods M3 and M4 yield small errors due to using both types of data. In other words, when using both ALT and ADT data the estimation accuracy of the parameter μ will be improved. Similarly, in Fig. 7, methods M2 and M4, which take into account the measurement error are more accurate in estimating λ. In summary, the proposed method M4 outperforms the other alternatives in estimating model parameters.

Fig. 2. Simulated degradation paths of 5 units with measurement error (5 units).

No

Log-likelihood function

4.2. Case study on hydraulic piston pump In this section, a set of accelerated testing data from a hydraulic piston pump is used to demonstrate the applicability and validity of the proposed method. A typical structure for a hydraulic piston pump is shown in Fig. 8. The main function of the hydraulic piston pump is to generate hydraulic energy for the hydraulic system. The drive shaft is driven by an engine, and the cylinder block rotates with the pistons reciprocally within the cylinder block, while the valve plate and the swash plate are fixed. For each revolution, the pistons complete a backand-forth reciprocating motion in the piston chambers and convert mechanical energy into hydraulic energy. Previous studies have shown that the three main failure mechanisms of the hydraulic piston pump are fatigue, wear, and ageing. Among these failure mechanisms, the most important one is the wear of three

censoring time in ALT.

One can see from Table 3 that method M4 results in the lowest estimation errors in the B10 life (the time at which 10% of units have failed) and parameters among the five methods. Method M0, which only uses ALT data, has a better estimate of μ than both methods M1 and M2, while the estimation error of λ is greater than 150%, which is too high to be credible. In other words, because of the insufficient number of ALT data, this method is inadequate in estimating model parameters. However, if only ADT data are used (i.e., method M1), the estimation 5

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

Table 3 Average values of parameter estimates with different methods. Type

μ

Errors (%)

λ

Errors (%)

σ

B10

Errors (%)

True value M0 (ALT) M1 (ADT) M2 (ADT (ɛ)) M3 (ADT + ALT) M4 (ADT(ɛ) + ALT)

4 3.9996 3.7269 3.7662 3.9021 3.9143

– 0.0096 6.8281 5.8443 2.4479 2.1423

10 25.1023 9.5316 9.7078 9.6412 9.8201

– 151.0227 4.6838 2.9222 3.5884 1.7992

0.1 – – 0.041 – 0.0413

42.94 49.23 48.35 47.61 44.6 44.5

– 14.65 12.60 10.88 3.87 3.63

Fig. 6. Boxplots of estimation errors of μ for four methods. Fig. 4. True and estimated reliability functions using methods M0 M1 and M4.

Fig. 7. Boxplots of estimation errors of λ for four methods.

= 1.663, = 0.6527 and = 0.00122 . In the basic wear process, it is possible to obtain the relationship between the wear and life of the component. As shown in Fig. 10, the basic wear process can be generally divided into three stages. Stage I is called running stage with a usually high wear rate. Adhesion, abrasion, occlusion or severe wear appear during this stage. Unstable working conditions or a larger impact of the assembly gap will make the wear in stage I turn into progressive wear that accelerates the process. Stage II is called stable wear with a fairly constant wear rate under normal operating conditions. Therefore, it can be considered that the wear process linearly evolves over time during this stage. Stage III is called sharp wear and presents a severe increase of the wear rate. In this stage, as the lubrication of the system deteriorates and the accumulation of abrasive particles increases, the surface temperature raises which accelerates the wear process. When the performance of the product reaches stage III, the system quickly becomes unstable and eventually reaches to the failure threshold D (at the failure time TD) when it stops working. From a practical point of view, define the soft failure time as the time to reach stage III, TDS. In particular, define the DS as the failure threshold of hydraulic piston pump, which is 2.15 L/min. This is different from what was used in the

Fig. 5. True and estimated reliability functions with methods M2 M3 and M4.

friction pairs (i.e., swash plate/slipper, valve plate/cylinder block, and piston/cylinder bore) [11]. As the wear and tear of friction pairs accumulates, the pump will fail eventually due to fluid leakage. The return oil flow is often used to indicate the performance degradation of the hydraulic pump. For an ADT, a flow sensor was mounted on the hydraulic pump to record the return oil flow during testing. Three pumps were tested under an accelerated discharge pressure Pa = 27.6 MPa (4000 psi) with a rotating speed Va = 6000 r/min, and the observations data were collected every 10 h for 220 h [32]. The return oil flow rates of the 3 pumps over time are shown in Fig. 9. In addition to the ADT data, ALT was conducted on another 8 pumps, and the test results are shown in Table 4. The testing condition was the same as the one applied in ADT. In this study, method M4 is applied to ADT and ALT data presented in Fig. 9 and Table 4. For this kind of pump, the time transformation t t (t ) = exp( ( ) ( ) 1) is utilized according to a previous study [33], where ν is the shape parameter, ω (0 < ω < 1) and α are the scale parameters. In [33], the parameters in time transformation are 6

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

Fig. 8. Typical structure of a hydraulic piston pump.

Fig. 10. Typical mechanical wear process.

Table 5 Estimated parameters and the 95% confidence intervals. Parameters

Point estimate

Lower bound

Upper bound

µ^ ^

2.236

1.9054

2.4726

47.152

19.9116

67.29

^

0.0103

0.0101

0.0196

Fig. 9. ADT data of three hydraulic piston pumps.

Table 4 The ALT data of hydraulic piston pumps (8 units). No.

Time (hours)

No.

Time (hours)

1 2 3 4

77.833 140.5 200 125.833

5 6 7 8

172.583 203.333 157.667 300*(censoring)

⁎

Fig. 11. The CDF estimate using method M4 with a 95% bootstrap confidence interval.

censoring time in ALT.

reference. Using the EM algorithm proposed in this paper, the parameter estimates were obtained as µ^ =2.236, ^ = 47.152, ^ = 0.0103 . To quantify the estimation uncertainty, a parametric bootstrap method was used to find the confidence intervals (CIs) of parameters and the reliability of this type of pump. In this paper, 1000 bootstrap replications were generated, and the 95% CIs of the model parameters are shown in Table 5. The CDF of failure time can be calculated from Eq. (5) using the estimated parameters, and the CI of the CDF was obtained and is shown in Fig. 11. Wang et al. [12] investigated the life-stress relationship for the same hydraulic piston pump. When effects of oil temperature are neglected, the life-stress relationship for the hydraulic piston pump can be described by a power law model, and the acceleration factor is given by:

AF =

Pab Vac P0b V 0c

P0 = 20.7 MPa (3000 psi) and V0 = 4000 r/min, the acceleration factor of the accelerated tests can be calculated as AF=6.342 . In addition, the mean-time-to-failure (MTTF) can be calculated as follows:

R (t|µ^ , ^,

MTTF =

(t )) dt

0

^

=

D

0

^ D

D µ^ (t ) +

(t ) D µ^

+ exp

2 ^ (t ) µ^

dt (20)

Using the estimated model parameters, the MTTF of this type of pump under the accelerated condition was obtained as 1037.5 h, the supplier of this type of pump recommends a first maintenance after 1200 h of operation under this accelerated condition. The prediction error is 13.5%, the prediction accuracy will further improved with the addition of more samples. As a result, the MTTF under normal operation will be 6580 h. It is worth pointing out that it is extremely important to use both ALT and ADT data for improving the accuracy of reliability estimation for this type of pump. Moreover, the estimation

(19)

where P0 and V0 are the discharge pressure and rotating speed of pump under normal operation. The value of parameters in Eq. (19) were obtained from an experiment where b = 6.4011 and c = 0.0141 [8]. Since the normal operation condition of this kind of pump is 7

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

results can be used to predict the remaining useful life of an operating pump based on its current degradation and a predetermined failure threshold. During operation, such results are quite valuable for guiding maintenance activities to ensure operational safety and for suggesting further reliability improvement.

information about the pump's performance. The results show the viability of the proposed model in improving reliability estimation by fully utilizing all the available data. Such results can be used for planning preventive maintenance actions during operation and to inform design changes during the product design phase. There are different ways to extend the research reported in this paper. For example, future work may incorporate random effects (e.g., variation between the degradation parameters across units) into the proposed IG process model by using the two types of accelerated testing data. Another line of research is to consider ALT and ADT with multiple stresses and non-constant stress profiles. These research topics will be studied in our future work.

5. Conclusion In this paper, a new reliability estimation method is proposed to fully utilize both ALT data and ADT data with measurement error. To this end, a novel model based on an IG process was introduced. An EM algorithm was developed to effectively estimate the model parameters. Numerical studies show that the proposed method improves the accuracy of parameter estimation using different sources of accelerated testing data. In addition, it is worth pointing out that considering ADT data with measurement error significantly reduces the estimation error of the diffusion parameter. To demonstrate the proposed approach, the ALT and ADT data of a hydraulic piston pump were analyzed. In practice, the ALT data set usually has a limited sample size. To overcome the shortcomings of using a small sample size, the ADT data was utilized to obtain more

Acknowledgments This work is supported by the National Natural Science Foundation of China (grant nos. 51620105010, 51575019), Natural Science Foundation of Beijing Municipality (L171003), Program 111 of Chinaand the program of China Scholarship Council (no. 201806020035).

Appendix A The five different methods which consider different factors are tabulated in Table 1. The log-likelihood functions of the methods can be derived as follows: M0: one of the simplest method M0 considering only ALT data is applied. According to Eq. (13), the log-likelihood function M0 can be simplified as: M0

=

t ,t *)

ALT (µ, NL q= 1

µ (tq))2

(D

1 ln + ln (tq) 2

=

2µ2D

NL*

+

ln r=1

D µ3

µ (tr*) (tr*)/

(A1)

Note that only ALT data are used for estimating the model parameters μ and λ. The first derivatives of calculated as follows. M0

µ

NL

µ (tq))2

(D

=

µ3

q= 1

NL*

+

(z ) (z )

r=1

1 µ 2

3D 2 (tr ) µ3

µ (tr ) ,

M1

with respect to the parameters are

(A2)

and M0

NL

= q= 1

µ (tq))2

(D

1 2

2µ2D

NL*

+

r=1

(z ) D µ (tr ) (z ) 2 µ3 (tr )

(A3)

M1: considers ADT date while ignoring measurement error, the log-likelihood function of this model is: M1

=

ADT (µ,

y) 2

ND

m

= i=1 j=1

1 ln 2

yi, j + ln

i, j

µ

i, j

2µ2 yi , j (A4)

Note that method M1 only considers ADT data for estimating the model parameters without taking into account measurement error. The first derivatives of M1 with respect to the parameters are calculated as follows. M1

µ

ND

m

=

( yi, j

µ

i, j)

2

µ3

i=1 j=1

,

(A5)

and M1

ND

m

= i=1 j=1

1 2

( yi, j 2µ2

µ

i, j)

2

yi, j

(A6)

N m N m ^ = ND m (N *m)/( y µ i, j ) . Then, the MLE of model parameters are µ^ = i =D1 j = 1 yi, j / i =D1 j = 1 i, j and D i, j i=1 j =1 M2: this method considers ADT data with measurement error. Similarly to method M0, the log-likelihood function M2 can be simplified as:

8

Reliability Engineering and System Safety 193 (2020) 106610

Z. Ma, et al.

M2

=

ADT (µ, ND

=

,

Z)

ln f ( Zi µ, , )

(A7)

i=1

as.

M3: Method M3 uses ALT data and ADT data without considering measurement error. The corresponding log-likelihood function can be described

M3

=

ADT (µ ,

y) +

t ,t *)

ALT (µ ,

2

ND i=1

=

m j =1

+

NL q=1

+

NL* ln r=1

1 2

1 2

ln

ln

yi, j

+ ln

i, j

i, j

2µ2 yi, j µ (tq))2

(D

+ ln (tq)

µ

2µ2D

µ (tr*)

D

(A8)

µ3 (tr*)/

M4: The proposed method M4 considers both ADT data with measurement error and ALT data. M4

=

ADT (µ,

=

ND i=1

+

,

Z) +

ALT (µ ,

ln f ( Zi µ , , ) + NL* ln r=1

D

NL q= 1

t ,t *) 1 2

ln + ln (tq)

(D

µ (tq))2 2µ2D

µ (tr*)

(A9)

µ3 (tr*)/

References [1] Nelson W. Accelerated testing: statistical models, test plans, and data analysis 344. John Wiley & Sons; 2009. [2] LuValle M. A theoretical framework for accelerated testing. Recent advances in reliability theory. Boston, MA: Birkhäuser; 2000. p. 419–33. https://doi.org/10. 1007/978-1-4612-1384-0_27. [3] Limon S, Yadav OP, Liao H. A literature review on planning and analysis of accelerated testing for reliability assessment. Qual Reliab Eng Int 2017;33(June (8)):2361–83https://doi.org/10.1002/qre.2195. [4] Luo W, Zhang C, Chen X, Tan Y. Accelerated reliability demonstration under competing failure modes. Reliab Eng Syst Saf 2015;136(April):75–84https://doi. org/10.1016/j.ress.2014.11.014. [5] Nelson WB. A bibliography of accelerated test plans. IEEE Trans Reliab 2005;54(2):194–7https://ieeexplore.ieee.org/iel5/24/30936/01435710. [6] Escobar LA, Meeker WQ. A review of accelerated test models. Stat Sci 2006;21(Sept. (4)):552–77https://doi.org/10.1214/088342306000000321. [7] Elsayed EA. Overview of reliability testing. IEEE Trans Reliab 2012;61(June (2)):282–91https://doi.org/10.1109/TR.2012.2194190. [8] Ma Z, Wang S, Zhang C, Tomovic MM, li T. A load sequence design method for hydraulic piston pump based on time-related markov matrix. IEEE Trans Reliab 2018;67(May (3)):1237–48https://doi.org/10.1109/TR.2018.2830330. [9] Meeker WQ, Escobar LA, Lu CJ. Accelerated degradation tests: modeling and analysis. Technometrics 2012;40(2):89–99https://doi.org/10.1080/00401706.1998. 10485191. [10] Han L, Wang S, Zhang C. A partial lubrication model between valve plate and cylinder block in axial piston pumps. ARCHIVE Proc Inst Mech Eng Part C 2015;229(17):1989–96https://doi.org/10.1177/0954406214568824. [11] Ma J, Ruan L, Fu Y, Chen J, Luo J. Research on current situation and methods of accelerated life test of aircraft hydraulic pump (serial 2)—typical failure modes and accelerated life test methods for aircraft hydraulic pump. Chin Hydraul Pneum 2015;7:1–6https://doi.org/10.11832/j.issn.1000-4858.2015.07.001. [12] Wang S, Shi J. Accelerated life testing under variable synthetic stresses for hydraulic pump. Mater Sci Forum Jan. 2006;505:1165–70https://doi.org/10.4028/ www.scientific.net/MSF.505-507.1165. [13] Peng C, Li X, Yuan Q. Life evaluation methods based on laboratory and field degradation data. 2017 annual reliability and maintainability symposium (RAMS). 2017. p. 1–6https://doi.org/10.1109/RAM.2017.7889741. [14] Wang X, Guo B, Cheng Z. Reliability assessment of products with Wiener process degradation by fusing multiple information. Acta Electron Sinica 2012;40(5):977–82. [15] Peng W, Li Y, Yang Y, Mi J, Huang H. Bayesian degradation analysis with inverse Gaussian process models under time-varying degradation rates. IEEE Trans Reliab 2017;66(March (1)):84–96https://doi.org/10.1109/TR.2016.2635149. [16] Pan Z, Balakrishnan N. Reliability modeling of degradation of products with multiple performance characteristics based on gamma processes. Reliab Eng Syst Saf 2011;96(8):949–57https://doi.org/10.1016/j.ress.2011.03.014. [17] Rodríguez-Picón LA, Rodríguez-Picón AP, Méndez-González LC, Rodríguez-Borbón

[18] [19] [20] [21] [22]

[23] [24] [25] [26] [27] [28] [29]

[30] [31] [32] [33]

9

MI, Alvarado-Iniesta A. Degradation modeling based on gamma process models with random effects. Commun Stat-Simul Comput 2018;47(6):1796–810https:// doi.org/10.1080/03610918.2017.1324981. Lim H, Yum BJ. Optimal design of accelerated degradation tests based on Wiener process models. J Appl Stat 2011;38(2):309–25https://doi.org/10.1080/ 02664760903406488. Tsai CC, Tseng ST, Balakrishnan N. Optimal design for degradation tests based on gamma processes with random effects. IEEE Trans Reliab 2012;61(2):604–13. 2012 https://doi.org/10.1109/TR.2012.2194351. Peng W, Li Y, Yang Y, Huang H, Zuo MJ. Inverse Gaussian process models for degradation analysis: a Bayesian perspective. Reliab Eng Syst Saf 2014;130:175–89https://doi.org/10.1016/j.ress.2014.06.005. Ye Z, Chen L, Tang L, Xie M. “Accelerated degradation test planning using the inverse Gaussian process. IEEE Trans Reliab 2014;63(3):750–63https://doi.org/10. 1109/TR.2014.2315773. Méndez-González LC, Rodríguez-Picón LA, Valles-Rosales DJ, Romero-López R, Quezada-Carreón AE. Reliability analysis for electronic devices using beta-Weibull distribution. Qual Reliab Eng Int 2017;33(8):2521–30https://doi.org/10.1002/qre. 2214. Li J, Wang Z, Liu X, Zhang Y, Fu H, Liu C. Wiener process model for accelerated degradation analysis considering measurement. Microelectron Reliab 2016;65(October):8–15https://doi.org/10.1016/j.microrel.2016.08.004. Pan D, Wei Y, Fang H, Yang W. A reliability estimation approach via Wiener degradation model with measurement errors. Appl Math Comput 2018;320(March (1)):131–41https://doi.org/10.1016/j.amc.2017.09.020. Ye Z, Wang Y, Tsui KL, Pecht M. Degradation data analysis using Wiener process with measurement errors. IEEE Trans Reliab 2013;62(Dec. (4)):772–80https://doi. org/10.1109/TR.2013.2284733. Ma J, Song Y, Lu Y, Qi X, Fu Y. Affecting factors on life time of aircraft hydraulic pump. Chin Hydraul Pneum 2016;0(07):6–12https://doi.org/10.11832/j.issn. 1000-4858.2016.01.002. Chhikara R, Folks L. The inverse Gaussian distribution: theory, methodology, and applications. New York: Marcel Dekker; 1989. Ye Z, Chen N. The inverse Gaussian process as a degradation model. Technometrics 2014;56(3):302–11https://doi.org/10.1080/00401706.2013.830074. Ruiz C, Liao H, Pohl EA. Reliability demonstration tests considering performance degradation with measurement error. 10th IWA international conference on modelling in industrial maintenance and reliability. 2018https://doi.org/10.19124/ ima.2018.001.23. Manjunath BG, Wilhelm S. Moments calculation for the double truncated multivariate normal density. 2009 Available at SSRN 1472153. Caflisch RE. Monte carlo and quasi-monte carlo methods. Acta Numer 1998;7:1–49https://doi.org/10.1017/S0962492900002804. Aerospace—Military type variable delivery, pressure compensated hydraulic pump, AS19692B, Jun. 22, 2016. https://doi.org/10.4271/AS19692B. Ma Z, Wang S, Liao H, Zhang C. Engineering driven performance degradation analysis of hydraulic piston pump based on the inverse Gaussian process. Qual Reliab Int 2019:1–19https://doi.org/10.1002/qre.2502.