Accelerated degradation tests with inspection effects

Delft University of Technology

Accelerated degradation tests with inspection effects

Zhao, Xiujie; Chen, Piao; Gaudoin, Olivier; Doyen, Laurent

DOI

10.1016/j.ejor.2020.11.041

Publication date

2021

Document Version

Final published version

Published in

European Journal of Operational Research

Citation (APA)

Zhao, X., Chen, P., Gaudoin, O., & Doyen, L. (2021). Accelerated degradation tests with inspection effects.

European Journal of Operational Research, 292(3), 1099-1114. https://doi.org/10.1016/j.ejor.2020.11.041

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

European

Journal

of

Operational

Research

journalhomepage:www.elsevier.com/locate/ejor

Innovative

Applications

of

O.R.

Accelerated

degradation

tests

with

inspection

effects

Xiujie

Zhao

a

_,

_Piao

_Chen

b ,∗

_,

_Olivier

_Gaudoin

c

_,

_Laurent

_Doyen

c

a College of Management and Economics, Tianjin University, Tianjin, China

b Department of Applied Mathematics, Technische Universiteit Delft, Delft, the Netherlands c Laboratoire Jean Kuntzmann, Université Grenoble Alpes, Grenoble, France

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 6 January 2020 Accepted 27 November 2020 Available online 15 December 2020

Keywords:

Reliability

Accelerated degradation tests Confidence density Degradation modeling Wiener process

a

b

s

t

r

a

c

t

Thisstudyproposesaframeworktoanalyzeaccelerateddegradationtesting(ADT)datainthepresence ofinspectioneffects.Motivatedbyarealdatasetfromtheelectricindustry,westudytwotypesofeffects inducedbyinspections.Aftereachinspection,thesystemdegradationlevelinstantaneouslyreducesbya randomvalue.Meanwhile,thedegradingrate iselevatedafterwards.Consideringtheabsenceof obser-vationsduetopracticalreasons,weemploytheexpectation–maximization(EM)algorithmtoanalytically estimatetheunknownparametersinastepwiseWienerdegradationprocesswithcovariates.Moreover, tomaintainthelevelofgeneralityfortheadaptionofthemethodinvariousscenarios,aconfidence den-sityapproachisutilizedtohierarchicallyestimatetheparametersintheaccelerationlinkfunction.The proposedmethodscanprovideefficientparameterestimationundercomplexlinkfunctionswithmultiple stressfactors.Further,confidenceintervalsarederivedbasedonthelarge-sampleapproximation. Simu-lationstudiesandacasestudyfromSchneiderElectricareusedtoillustratetheproposedmethods.The resultsshowthattheproposedmodelyieldsaremarkablybetterfittotheSchneiderdataincomparison totheconventionalWienerADTmodel.

© 2020TheAuthor(s).PublishedbyElsevierB.V. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introductionandmotivation

Reliability tests arewidely used to predictproduct lifetime in various industries.A successfullyplanned andconducted reliabil-itytest can provideimportantinformationsupporting managerial decisions under a reasonable test budget, thereby reducing both prospectivecosts andrisks.In orderto shortenthe testduration, conventional life tests are commonly conducted under elevated stresses toaccelerate the failuresoftest units.Inrecent decades, withtheadvancesinsensors andmonitoringtechnologies, degra-dation testsbecomepreferabletolife testsinthesensethat they can predict reliability characteristics over time without the con-cernofcensoring(Meeker, Escobar, & Lu, 1998 ).Inatypical degra-dation test,discrete degradationmeasurementsare takenandthe observeddegradationpathsarethenemployedtomakeinferences aboutthe product reliability. Degradationpaths areusually mod-eled basedonaqualitycharacteristic(QC),suchasthebrightness ofdisplaysandbatterylifeofelectronicdevices(Wang, Tang, Bae, & Xu, 2018 ).

∗_{Corresponding author.}

E-mail address: p.chen-6@tudelft.nl (P. Chen).

When utilized for lifetime prediction, the structure and pat-tern of degradation data are desired to be as simplistic as pos-sible to relievethe modeling complexity andcomputational bur-den.However,underpracticalusageorevencontrolled experimen-tal conditions, degradation paths may inevitably behave atypical patternsfrom time to time.In degradation tests, since thestress levels are strictly controlled, the degradation paths of test sys-temsareusuallydeemedtobestableduringmostofthetesttime. Nevertheless,some interventions tothe test systemsmay be un-avoidableduetovariouspracticalconcerns.Acommonexerciseis thatengineershavetoalterthetestconditionstemporarilyto ob-taindegradationmeasurements.Forexample,manyreliabilitytests areconductedundercertaincombinationsoftemperatureand hu-midity,where test chambersthat provide such environments are employed. Typicaldegradation testsof thiskindcan be found in Meeker and Escobar (1998 , Chapter 21) and references therein. To take effective degradation measurements, the involvement of manual inspection or/and preciseinstruments are mandatory,yet the exercise isdifficult,if not impossibleunderthe test environ-ment.In such cases,the test unitshave to be removedfromthe chambertemporarily,whichresultsinthechangeoftest environ-ments.Although thedurationof measurementisusually shortor even negligible compared to the whole test duration, the drastic change of test environments may still cause substantial changes

https://doi.org/10.1016/j.ejor.2020.11.041

onsystemdegradationlevels.Anotherreasonofatypical degrada-tion paths lies in the interveningnature of inspections. In other words, certain types of inspectionsmay inherently influence the system degradation (Zhao, Gaudoin, Doyen, & Xie, 2019 ). One of the examplesis thedestructive test,whereinspectionscan cause directly destructive effects to the test systems (Shi, Escobar, & Meeker, 2009 ).Consideringtheaforementionedissues,wepropose to modelthe effects broughtby inspections in degradation tests andafterwardsinvestigateparameterestimationuponsuchtesting data.

Theresearch tobe proposedismotivatedbyarealexperiment carried out by SchneiderElectricwiththe objectiveto revealthe degradation characteristicofanelectrical distributiondevice.Asa key partof the device, a mechanical linkage corrodesover time, whichisadominantcauseofperformancedegradation.Toquantify thedegradation levelofthedevice, engineersmeasurethetorque that is neededtoseparate the linkage.Ahigher torqueimpliesa moresevereconditionofcorrosion.Ontheonehand,sincethe in-spection separates the linkage, the grown corrosion is physically disassembled,leadingtoareductioninthedegradationlevel dur-ingtheinspection.Ontheotherhand,theinspectioncauses dam-age to the integrity of surface treatments in the linkage, which leads to a higher rate of corrosion. Obviously, the two types of effects are opposite with respect to the system health. Another concern ofthe problemisthe difficulty inrevealingthe accurate degradationreductionduringtheinspection.Tofollowabasicrule toinspectsystems,engineerstendtominimizetheinfluenceof in-spectionandthereforeonlymeasurethetorquethatseparatesthe linkage.Onceameasurementisobtained,theinspectionis termi-nated immediately. Consequently, the measurement process may only capture the degradation level before inspections yet fail to observedegradationreduction.Althoughtheengineerscangivean approximationofthedegradationreductionfromtheprior knowl-edgeorotherexperiments,theaccuratevaluecanneverbeknown. Thus, thereduction effectis ahiddenvariablethat cannot be di-rectlyutilizedforstatisticalinference.

ByemployingtheexperimentfromSchneiderElectricasan il-lustrativeexample,thestudyaimsatestablishinga frameworkto analyzeaccelerateddegradationtestswithcomplexinspection ef-fects. In general, the proposed method can be applied to model degradation datainthepresenceofenvironmentalcovariates and interventionsthatexertbothpositiveandnegativeeffects.

2. Literaturereview

With the fast emergence of system monitoring technologies, modelingandinferenceofdegradationdatanowplayavitalrolein the researcharea ofreliabilityengineeringandits interfaceswith other areas,such asmechanicalengineering(Wang & Tsui, 2017 ), energy(Lin et al., 2017 ),electricalengineering(Si, 2015 )etc. Degra-dation analysisnot onlysubsumes approaches to modelrelevant data, but also creates alternative planning methods of reliability testsforlifeprediction.

Interest inaccelerateddegradationtest(ADT)hasgrown in re-centdecadesowing toitssuccessfulapplicationstovarious prod-ucts and systems, such as LED lamps, lithium-ion batteries and rail tracks (Ye & Xie, 2015 ). Initiated by Meeker et al. (1998) , regression-basedgeneralpathmodelsarewidelyusedfor degrada-tion modelingin ADT(Hong, Duan, Meeker, Stanley, & Gu, 2015 ). Further, dueto clear physicalexplanations and appealing mathe-maticaltractability,stochasticprocessessuch asWienerprocesses (Hu, Lee, & Tang, 2015 ),gamma processes(Tsai, Sung, Lio, Chang, & Lu, 2016 )andinverseGaussianprocesses(Ye, Chen, Tang, & Xie, 2014 ) startto playan influentialrole inADT modelingand plan-ning. Foramoredetailedoverview,onecanbe referredtoLimon, Yadav, and Liao (2017) .Todate,severalnewconsiderationsonADT

planning and analysis emerged in the literature to tackle more practical issues.Tonameafew,Tseng and Lee (2016) proposed a generalexponential-dispersion model to characterizedegradation andgavetheoptimalADTplansinanalyticalforms.InLi, Wu, Ma, Li, and Kang (2018) ,random fuzzytheory was adoptedto model the uncertainty in ADT data. Wang and Tsui (2017) considered multiplestressesinADTforrubbersealedO-rings.Withregardto modeluncertainty,Liu, Li, Zio, Kang, and Jiang (2017) appliedthe BayesianmodelingaveragingapproachtomodelADTdata.

In real problems, it is common that a degradation path can-not be characterized by models in regular forms such as linear or typical nonlinear ones (e.g., polynomial, logarithm and expo-nentialmodels).The reasonsbehindan atypicaldegradation path canbe rathercomplex. Forexample,Hong et al. (2015) proposed a degradation modeling approach by utilizing dynamic weather-ingcovariatestocharacterizeirregulardegradationpaths.Insome other studies (Bae, Yuan, Ning, & Kuo, 2015; Wang et al., 2018 ), change-pointdetectionandmodelingwerediscussedfor degrada-tion data.Adaptive anddynamicmethods foronline degradation modelinghavealsoprevailedintheliterature(Si, 2015; Zhai & Ye, 2018 ).

Apart from these, human interventions to industrial systems can also be a pivotal cause of atypical degradation path, and a common example is imperfect maintenance (Mercier & Cas- tro, 2019 ). Surprisingly,despiteplentiful extantworksonatypical degradationpaths,wecanonlyfindveryscantresearchinthe lit-erature that addressed similar issues in ADT problems. Xiao and Ye (2016) discussed the ADT planningproblemwith random ini-tial degradation levels. In Ye, Hu, and Yu (2019) , the initial per-formanceof testunits wereconsidered to allocateunitstostress levels.Nevertheless,theseworksdidnotincorporatetheeffectof inspectionasdescribedinSection 1 .

Therestofthepaperisorganizedasfollows.InSection 3 ,a sys-tematicapproachtomodelconstructionandparameterestimation isestablished forADTdatawithinspection effects.Section 4 dis-cussestheuncertaintyquantificationofparameterestimators. Sim-ulationstudiesarecarriedoutinSection 5 .Section 6 presentsthe casestudyfromSchneiderElectric.Finally,conclusionsare drawn inSection 7 .

3. Degradationmodelswithinspectioneffects

3.1. Preliminaries

Consider a degradation test with M stress levels and Ni test

unitsareallocatedtostressleveli_,wherei₌1_,...,M.Forthe jth

testunitunderstressi,aswhichwecallunit

(

i,j

)

forsimplicity,a totalofOi jinspectionsarecarriedoutattimeepochs

τ

i j1,...,

τ

i jO_{i j}.

Since an instant degradation reduction occursupon each inspec-tion, we denotethe degradation level before the reductiveeffect byy−_{i jk} forthekthinspectionforunit

(

i,j

)

,withk=0representing theinitialinspectionpriortothetest.Meanwhile,thedegradation levelafterthe reductiveeffectisdenoted byy+ _{i jk}.Weintroduce a variablezi jktomodeltheproportionofdegradationreductionwith

respecttoy−_{i jk} forthekthinspectionofunit

(

i,j

)

asfollows:

zi jk=

y−_{i jk}− y+ i jk

y−_{i jk} , 0≤ zi jk≤ 1. (1)

Thenotationaldetailsareillustrated inFig. 1 .Asmentionedin Section 1 ,y−_{i jk} isusuallyobservable whereasy+ _{i jk} isrelatively diffi-culttoreveal.Thus,fornotationalsimplicity,wesubsetthedatato

y− and zthat are given by

y−=

{

y−_{i jk},i=1,...,M,j=1,...,Ni,k=0,...,Oi j

}

,

Time

Degradation path

Fig. 1. Illustration for notations by a possible sample path of degradation for unit j tested under environment i .

Tofurtherfacilitatedegradationmodeling,we let



y=

{



y_{i jk},i= 1,...,M,j=1,...,Ni,k=1,...,Oi j

}

, where



yi jk=y−_{i jk}− y+ _{i j(k−1)}.

Consideringthenotationin(1) ,wecanrewrite



yi jk as



yi jk=y−_{i jk}− y−_{i j(k−1)}

(

1− zi j(k−1)

)

.

3.2. Wienerdegradationmodelandinspectioneffect

We employtheWienerprocessasthebaselinemodelto char-acterizetheinherentgenerationofdegradationforsystemof inter-est.Moreconcretely,ifweassume thatthesystemoperates with-out anyintervention, the degradation path of thesystem can be modeledbyadriftedWienerprocess.TheWienerprocessfeatures independent Gaussian increments over non-overlappingtime pe-riods, which enables its wide application in degradation model-ing.AdriftedWienerprocess

{

W

(

t

)

; t≥ 0

}

canbecharacterizedby drift anddiffusionparameters, denotedby

μ

and

σ

_,respectively. In thismanner, itgivesW

(

t

)

=

μ

t+

σ

B(t

)

,whereB(·

)

isa stan-dard Brownian motion. The increment



W

(

t− s

)

=W

(

t

)

− W

(

s

)

foranyt>sfollowsanormaldistributionwithmean

μ

(

t− s

)

and variance

σ

2

(

t− s

)

, i.e.,



W

(

t− s

)

∼ N

(

μ

(

t− s

)

,

σ

2

₍

_{t− s}

₎₎

_{.In the} presence ofinspection effects, i.e., nonzero zi jk’s, thedegradation

pathnolongerfollowstheconventionalWienerprocess. Neverthe-less,accordingtotheaforementionedproperties,givenz_{i jk}’s,



y_{i jk}

are independentincrementsofa Wienerprocess andthey follow:



yi jk∼ N

(

μ

i jk



τ

i jk,

σ

2



τ

i jk

)

. (2)

where

μ

_{i jk} isthedegradation ratebetweenthe

(

k− 1

)

thandkth inspectionfortestunit

(

i,j

)

and



τ

i jk=

τ

i jk−

τ

i j(k−1).Notethatif

theinspectionhasnoeffectonthesystemdegradation,thenzi jk’s

are0foralli_,j andk_,andthedegradationpathcanbemodeledby a conventional driftedWienerprocess.Note thatinthe proposed model, we use a flexible notational convention in the sense that

τ

i jk andOi j can bedifferentfordifferenttestunits,implyingthat

theproposedmethodscanbewellappliedtounbalanceddataset intermsoftimeandnumberofobservations.

Next,experimentalfactorsascovariatesareintroducedintothe model.Themajorityofresearchondegradationtestshassuggested touseparametricmodelsto link

μ

andcovariatesxi (Jakob, Kim-

melmann, & Bertsche, 2017 ).Weusefacc

(

xi

)

todenotethebaseline

degradation rate under stress xi and the acceleration model can

beformulatedineitherparametricornonparametricmanners.The

parametricaccelerationmodelcantakevariousformsbasedonthe physicalmechanismofdegradationandfactorsinvolvedinthetest. The selection of facc

(

xi

)

is not the main focus of the study,and

somebriefdiscussionsregardingtheSchneiderexamplearegiven inAppendix A . Moreover,the complexityin facc

(

xi

)

mayimpede

theinferentialefficiencyundervariousspecificmodels.Under dif-ferentformsof f_acc

(

x_i

)

_,we wishtomaximize thegeneralizability oftheproposedmodel.Towardsthisend,wetreat facc

(

xi

)

as

sep-arateparameters first andthen employ a hierarchicalanalysis in Section 3.6 forfurtherinferences.

Tocapturetheeffectofdegradationrateincreaseafterthekth inspection,afunctiong

(

k;

ω

)

isintroduced,where

ω

isavectorof unknown parameters.Tobenchmarktheeffect,we useg

(

k;

ω

)

as an addedtermtothe baselinedegradation rateandtherefore we assume g

(

0;

ω

)

≡ 0.Then, the degradationratebetween the

(

k− 1

)

thandkthinspectionfortestunit

(

i,j

)

isgivenby

μ

i jk=facc

(

xi

)

+g

(

k;

ω

)

. (3)

Oneimplicitassumptionfrom(3) isthatthe increaseeffectof degradationratebroughtbyinspectionsandstressvariablesxiare

independent.Inotherwords,weassumethattheinspectioneffect ondegradationrateonlydependsonkanddoesnotinteractwith environmentalstresses. Sincethe degradation rateincreaseswith the number of inspections, g

(

k;

ω

)

is a non-decreasing function. Thesimplestformsarepolynomial,e.g.,thefirstandsecondorder polynomialmodelsaregivenasfollows:

g

(

k;

ω

)

=

ω

k, g

(

k;

ω

)

=

ω

1 k2 +

ω

2 k.

Without lossofgenerality, we assume g

(

k_;

ω

)

₌

ω

k foranalytical simplicityinthefollowinganalysis.

Recallthat wehaveintroducedzi jk earliertodescribethe

pro-portionofdegradation reduction.Abeta distributionisemployed tomodeltheunconditionedz_{i jk},foranyrealizationofz_{i jk},denoted byz,thedensityfunctionisgivenby

fBeta

(

z; u,

v

)

=

_(

(

_uu

_)(

+

v

_v

)

₎

zu−1

(

1− z

)

v−1 , 0≤ z≤ 1, (4) whereuand

v

aretheshape parameters ofthebetadistribution. Itisworthmentioningthatbetadistributionhasbeenwidelyused tomodeltheeffectoftheimperfectrepair(Zhang, Gaudoin, & Xie, 2015 ).

Remark. Ify+ _{i jk}isunobservable,uand

v

cannotbeestimatedfrom the modeldueto the absenceof zunder thefrequentist setting. However,engineersmaymanagetoapproximatelyquantifythe re-ductive effectfromdomainexpertiseor preliminaryexperiments. Forthe examplefromthe SchneiderElectric,engineers cancarry outadifferenttypeofexperimenttomeasurethedegradation lev-elsbeforeandafterthefirstinspectionatt=0andmodelthe re-ductiveeffects,thoughitisinapplicableduringtheADTasittakes much longer time to obtain the reductive measurements. In the presence ofthedescribedavailable data att=0,it is easyto fit abetadistributiontotheobservedreductivemeasurements.Other typesof models to characterizez can also be usedafter uncom-plicatedmodificationstothelikelihoods inthefollowingcontents in the section and density functions in Appendix B . It is note-worthythat theBayesianmethodcan beanappealingalternative, where u and

v

can be characterized by some prior distributions at first. The computational burden can be an important issue in Bayesian inference, especially when z is unobservable (Bernardo et al., 2003 ). Further,inthispaper,wepresume thatz_{i jk}’sare in-dependentrandomvariables.Theassumptionisvalidifthe inspec-tioneffectsofdegradationreductionareinstantaneousanddonot interactover time. In realistic applications, the effects can inter-actovertimeunderstressenvironments.Insuchcases,ajoint dis-tributioncanbe employedto characterizetheinterdependenceof

z_{i jk}’s.However,theintroductionofinterdependenceleadstomore complicated likelihoods,which hindersthe tractability of estima-tors.MonteCarloEMalgorithmcanbeusefulinimplementingthe aforementionedextensions(Levine & Casella, 2001 ).

InSections 3.3 and3.4 ,we willdiscussthemodelingof degra-dation data under two different scenarios. In the first one, the degradation levels before andafter the inspection are observable sothatalltheparametersinthemodelcanbeestimatedthrougha completelikelihood.Inthesecondone,onlythedegradationlevels before theinspection can be obtainedforinference,underwhich caseparametersuand

v

aregivenaprioritofacilitatethe estima-tionofotherunknownparameters.

3.3. Inspectioneffectmodelingwithcompleteobservations

Aspreliminariesforthefollowinganalyses,thispartaimsat es-tablishing data modeling framework via maximum likelihood es-timation (MLE).The complete log-likelihoodfunction of

θ

c under

data

(

y−,z

)

canberepresentedby logL

(

θ

c

|

y−,z

)

=logp



y−

|

z,

θ

c



+logp(z

|

θ

c

)

= M  i=1 N_i  j=1 Oi j  k=1



−1₂log2

π

−1₂log

τ

i jk− log

σ



− M  i=1 Ni  j=1 Oi j  k=1





y−_{i jk}− y− i j(k−1)



1− zi j(k−1)



−

μ

i jk

τ

i jk

2 2

σ

2

τ

i jk

+ M  i=1 Ni  j=1 O_i j+1 k=1



log

(

u+

v

)

− log

(

u

)

− log

(v

)

+

(

u− 1

)

logz_{i j(k}−1)+

(v

− 1

)

log



1− zi j(k−1)



, (5) where

μ

_{i jk=} facc

(

xi

)

+

ω

k. Note againthat we assume that zi jk’s

aremutuallyindependentandtheyarealsoindependentofyand the environmental factors. Under the current model setting, the unknown parameters in the model can be summarized by

θ

c=

(

facc

(

xi

)

,i=1,...,M,

ω

,

σ

2 ,u,

v

)

T.Ifthe effect ofdegradation

re-duction can be observed,the problemissimplified toa standard MLE problemwithallobservations available in(5) .Further,since parameters u and

v

are independent of other parameters in the model,thelikelihoodinvolvinguand

v

canbeindependently max-imized viathe observationsz_{i jk}.The remaining partofthe likeli-hoodcanalsobemaximizedbynumericalmethods.

3.4. Inspectioneffectmodelingwithhiddeneffectobservations

When y+ _{i jk} isunobservable,theeffectof degradationreduction cannot be captured, which makes u and

v

in the model ines-timable.Thiskindofinformationcanalsobequantifiedbya beta distribution asdescribed in (4) . As discussed before, we assume that the pilot distribution ofz is knownand characterized by u0 and

v

0 , respectively.By holdingthe property ofindependence of

zi jk’s.the followinglog-likelihoodisto bemaximizedto estimate

θ

=

(

f_acc

(

x_i

)

_,i₌1_,...,M_,

ω

_,

σ

2

)

: logL(

θ|

y−,z)=logp



y−

|

z,

θ



= M  i=1 Ni  j=1 Oi j  k=1



−1₂log2

π

−1₂log

τ

i jk− log

σ



− M  i=1 N_i  j=1 O_{i j}  k=1





y−_{i jk}− y− i j(k−1)



1− zi j(k−1)



−

μ

i jk

τ

i jk

2 2

σ

2

τ

i jk . (6) Duetotheabsenceofzi jk,thelog-likelihoodcannotbemaximized

in its current form. Alternatively, we resort to the expectation– maximization (EM) algorithm to obtain the parameter estimates.

TheEM algorithmisan iterativemethod tofindtheMLEfor sta-tistical models with latent variables. Following its first introduc-tion byDempster, Laird, and Rubin (1977) ,the EMalgorithm has beenstudied extensively fromboth theoretical andpractical per-spectives (McLachlan & Krishnan, 2007 ).Specifically, in reliability engineering,theEMalgorithmiscommonlyusedtocapturelatent random effect in life anddegradation models (Chen & Ye, 2017; Duan & Wang, 2018 ). Apart from the EM algorithm, the hidden semi-Markovmodel wasused forhealth diagnosisandprognosis withlatenteffects inDong and He (2007) .A two-stageapproach wasproposed inLee, Hu, and Tang (2017) to estimate themodel fromtime-censoredADTdata.TheKalmanfilteringtechniquehas also beenprevailingly employed inthe remaining life estimation based on degradation models (Si, Wang, Hu, & Zhou, 2014 ). The EM algorithm consistsof two steps: (1) the E-stepin which the conditionalexpectationofthecompletelog-likelihoodwithrespect toincompletedataiscompleted;(2)theM-stepinwhichthe ex-pectedlog-likelihoodismaximizedtogenerateparameter estima-tionatthecurrentiteration.Denotetheconditionalexpectationof thecompletelog-likelihoodatiterationnbyQ



θ|

θ

(n)



_,wehave Q



θ|

θ

(n)



= M  i=1 Ni  j=1 Oi j



−1 2log2

π

− 1 2log



τ

i jk− log

σ



− 1 2

σ

2 M  i=1 Ni  j=1 Oi j  k=1 1



τ

i jk



_

y−_{i jk}− y− i j(k−1) − facc

(

xi

)

τ

i jk−

ω

k



τ

i jk



2 +2y−_{i j(k−1)}



y−_{i jk}− y− i j(k−1) − facc

(

xi

)

τ

i jk−

ω

k



τ

i jk



E



z_{i j(k−1)}

|

y−,

θ

(n)



+



y−_{i j(k−1)}



2 E



z2 _{i j(k−1)}

|

y−,

θ

(n)



. (7) In the expectation step,we need to compute E



z_{i j(k−1)}

|

y−,

θ

(n)



andE



z2 _{i j(k−1)}

|

y−,

θ

(n)



,Unfortunately,theconditionaldistribution ofzi j(k−1)cannotbeidentifiedasaknownrandomdistribution.We

havetoresorttonumericalmethodstoevaluatetheexpected val-ueswithrespecttoz_{i j(k−1)}.RelateddetailsaregiveninAppendix B . Next,the firstorder derivativesofQ



θ|

θ

(n)



are providedfor its maximization. Following previous analyses, we also treat f_acc

(

x_i

)

asunknownparameters.Thus,

∂

Q



θ|

θ

(n)



∂

f_acc

(

x_i

)

=− 1

σ

2

facc

(

xi

)

N_i  j=1 O_{i j}  k=1



τ

i jk − N_i  j=1 O_{i j}  k=1



y−_{i jk}− y− i j(k−1)+y−i j(k−1)E



zi j(k−1)

|

y−,

θ

(n)



+

ω

Ni  j=1 O_{i j}  k=1 k



τ

i jk



. (8)

By solving

∂

Q



θ|

θ

(n)



/

∂

facc

(

xi

)

=0, we obtain the following

equation: f_acc (n+1)

(

xi

)

=A_i(n)− Bi

ω

(n), (9) where A(_in)= Ni j=1 Oi j k=1



y−_{i jk}− y− i j(k−1)+y−i j(k−1)E



zi j(k−1)

|

y−,

θ

(n)



Ni j=1 Oi j k=1



τ

i jk ,

Bi= Ni j=1 Oi j k=1 k



τ

i jk Ni j=1 Oi j k=1



τ

i jk .

Furtherfor

ω

,wehave

∂

Q



θ|

θ

(n)



∂ω

=− 1

σ

2

− M  i=1 Ni  j=1 Oi j  k=1 k



y−_{i jk}− y− i j(k−1) +y−_{i j(k−1)}E



z_{i j(k−1)}

|

y−,

θ

(n)



+ M  i=1 Ni  j=1 Oi j  k=1 k



τ

i jkfacc

(

xi

)

+

ω

M  i=1 Ni  j=1 Oi j  k=1 k2



τ

i jk



, (10) ofwhichthesolutionisfollowedby

C

ω

(n)=D(n)− M  i=1 facc

(

xi

)

Fi, (11) where C= M  i=1 Ni  j=1 Oi j  k=1 k2



τ

i jk, D(n)₌M i=1 Ni  j=1 Oi j  k=1 k



y−_{i jk}− y− i j(k−1)+y−i j(k−1)E



z_{i j(k−1)}

|

y−,

θ

(n)



_, Fi= Ni  j=1 Oi j  k=1 k



τ

i jk.

ToplugEq. (9) toEq.(11) ,thefollowingequationisobtained

C

ω

(n)=D(n)− M  i=1



A_i(n)− Bi

ω



Fi, (12)

whichyieldstheestimatesof

ω

and facc

(

xi;

δ

)

givenby

ω

(n+1)₌ D(n)− M i=1 A(in)Fi C−M i=1 BiFi , facc (n+1)

(

xi

)

=Ai(n)− Bi

ω

(n+1). (13)

Thenfor

σ

2_{,likewisewehave}

∂

Q



θ|

θ

(n )



∂σ

2 =− 1 2

σ

2 M  i=1 Ni  j=1 Oi j+ 1 2

σ

4 M  i=1 Ni  j=1 Oi j  k=1 1

τ

i jk



_

y−_{i jk}− y− i j(k−1) − facc

(

xi

)τ

i jk−

ω

k

τ

i jk

2

+2y−_{i j(k−1)}



y−_{i jk}− y− i j(k−1) − facc

(

xi

)τ

i jk−

ω

k

τ

i jk



E



zi j(k−1)

|

y−,

θ

(n )



+



y−_{i j(k−1)}

2

E



z2 i j(k−1)

|

y−,

θ

(n )



. (14)

Therootoftheequationcanbeeasilyobtainedbysolving



σ

2



(n+1)M i=1 Ni  j=1 Oi j= M  i=1 Ni  j=1 Oi j  k=1 1



τ

i jk



_

y−_{i jk}− y− i j(k−1) − facc

(

xi

)

τ

i jk−

ω

k



τ

i jk



2 +2y−_{i j(k−1)} ×



y−_{i jk}− y− i j(k−1)− facc

(

xi

)

τ

i jk −

ω

k



τ

i jk



E



z_{i j(k−1)}

|

y−,

θ

(n)



+



y−_{i j(k−1)}



2 E



z2 _{i j(k−1)}

|

y−,

θ

(n)



. (15) Through(9) –(15) ,thecurrentiterationoftheEMalgorithmis real-ized.Theiterationsarecontinueduntiltheconvergenceof param-eterestimators.

3.5. Guessofinitialestimatesandendingofiterations

To start the aforementioned EM algorithm, starting estimates

θ

(0) _{are needed. The convergence speed ofthe algorithm hinges}

on theselection of

θ

(0). Here, we utilizethe meanof the zi jk to

approximatelyobtain

θ

(0).ItisobviousthatE

(

z_{i jk}

)

=u₀ /

(

u₀ +

v

0

)

. Therefore,touseE

(

zi jk

)

ratherthanunobservablezi jk,wehave



y−_{i jk}−

v

0 y − i j(k−1) u0 +

v

0



∼ N



facc

(

xi

)

τ

i jk+

ω

k



τ

i jk,

σ

2



τ

i jk



, (16) where the left-hand-side term is completely observable and af-ter the manipulation of a typical MLE,

θ

(0) can be given as in Appendix C .

Another issue of the EM algorithm is when to terminate the iterations. The question poses a tradeoff between the estimat-ingprecision andcomputationalefficiency.Generally, itisa com-mon criterion to terminate the iterations when the proportions ofchangesinabsolutevaluesofparameterestimators aresmaller than criticalvalues

ε

. It isa vector because different parameters mayhavedifferentcriticalvalues.Fortheproblemdescribedinthe paper,parametersplay differentrolesdependingonhowdecision makers wouldutilizethe estimates.Forexample,in termsof life prediction, facc

(

xi

)

’sare importantfortheextrapolationto

under-stand the degradation rateunder normalusage conditions, espe-cially in the presence of conditionfluctuations, while

ω

is more usefulifthe deviceisfrequentlyinspected.Regarding these relia-bilityissues,

ε

canbeproperlydeterminedto satisfytherequired estimatingaccuracy.

3.6. Ahierarchicalanalysistoestimate facc

(

xi

)

Inaforementionedanalyses, facc

(

xi

)

,i=1,...,M are treatedas

unknownparametersforestimation.Asoneofthemainobjectives of the study, the estimation of degradation rate under normal usage condition isrealized by a hierarchicalmethod. Due to the possible complexity in f_acc

(

x_i

)

_, it could be onerous to derive analytical iterative solutions to the parameters herein to the EM algorithm. In view of this, the hierarchical method can pro-videreasonableestimationandmeanwhilekeepthemathematical derivationsinthepaperdirectlyadaptableinvariousscenarios.Liu, Liu, and Xie (2015) reported a method to conduct meta-analysis ofindependentstudiesviaaconfidencedensity(CD)approach.In thepaper,we employarevisedversion oftheconfidencedensity to estimate the parameters

δ

in facc

(

xi

)

, and we denote facc

(

xi

)

by facc

(

xi;

δ

)

inthefollowingcontext.AsdiscussedinAppendix A ,

the formof f_acc

(

x_i;

δ

)

dependsoncertain knownphysical mech-anisms. To ensure facc

(

xi;

δ

)

to be greater than zero, we take

naturallogarithmonit.DuetotheinvariancepropertyofMLE,the MLEoflogf_acc

(

x_i;

δ

)

isreadilygivenbylogfˆ_acc

(

x_i;

δ

)

.Threetimes differentiable mapping functionsM=

(M

1 ,...,MM

)

 is used to

linklogfacc

(

xi;

δ

)

andtheunknownparametervector

δ

:

logfacc

(

xi;

δ

)

=Mi

(

δ

)

. (17)

Further, following Xie and Singh (2013) , we can construct a CD for logfacc

(

xi;

δ

)

,i=1,...,M, which is a multivariate normal

(MN) distribution with mean



logfˆacc

(

x1 ;

δ

)

,...,logfˆacc

(

xM;

δ

)



and covariance matrix

{

[I

(

θ

ˆ

)

]−1  ˆ

θ

_θ

ˆT

_}

_{1: M, where}

_{

_·

_}

_{1: M denotes} the M× M square matrix partitioned from the upper left in the original matrix and _{ denotes} an element-wise division. Note that

δ

is only involved in facc

(

xi;

δ

)

, thus we only use the first

M rows and columns in the original covariance matrix in the

model described in Section 3.4 , and the parameters

ω

and

σ

2 are assumed known in the subsection. The covariance matrix is approximated by the delta methods. For notational convenience,

we let Vˆ =

{

[I

(

θ

ˆ

)

]−1  ˆ

θ

_θ

ˆT

_}

1: M. Additionally, we let V be the

covariance matrix under true parameters. Then the CD for

δ

is giveninaformofMNdistributionby

h

(

δ

)

=

(

2

π

)

−M/ 2 det

(

Vˆ

)

−1 /2 exp



−1 2a

T_Vˆ−1 _a



_{, (18)} where a is a M-dimensional column vector with each element givenby

ai=logfacc

(

xi;

δ

)

− logfˆacc

(

xi;

δ

)

, i=1,...,M.

Bymaximizing(18) ,weobtainthepointestimatorof

δ

underCD, whichwedenoteby

δ

ˆ_CD :

ˆ

δ

CD=argmax

δ h

(

δ

)

. (19)

The reason for using the CD estimation is twofold. First, as mentioned previously, it facilitates the derivation of closed-form estimators undertheEMframework.Second, theefficiencyof es-timation isnotcompromisedbythehierarchicaloperationsby CD estimation.Thefollowinganalysesarepresentedtojustifythe lat-terstatement.UnderaconventionalMLEframework,if

δ

isdirectly usedtomaximizethelikelihoodfunctionin(6) ,wecanobtainthe MLEof

δ

givenby

ˆ

δ

DIR =argmax δ L

(

δ

; y

−_,z

₎

_{. (20)}

For notational convenience, we let

ζ

=



logf_acc

(

x_{1 ;}

δ

)

_,...,logf_acc

(

xM;

δ

)



 _{and use L}

₍

_δ

₎

_{and L}

₍

_ζ

₎

_to

respectively represent the likelihood functions L

(

δ

; y−_,z

₎

_and

L

(

ζ

; y−_,z

₎

_,whereL

₍

_ζ

_{; y}−_,z

₎

_{denotesthelikelihoodunder}

_ζ

_.Thus, we have

ζ

=M

(

δ

)

. Moreover, let n=M

i=1 Nj=1 i Oi j be the total

numberofobservationsinthetest.

Lemma1. Asn_→∞,thedirectestimator

δ

ˆ_DIR obtainedfrom(20)is consistentandnormallydistributed.Specifically,

n1 /2



δ

ˆDIR −

δ



_d − →MN



0,



J

(

δ

)

TI˜

(

ζ

)

J

(

δ

)



−1



, (21)

where ˜I

(

ζ

)

=V−1 and J

(

δ

)

=

∂

M

(

δ

)

/

∂

δ

istheJacobianofM with respectto

δ

.

Theproofsofthelemmaandthefollowingresultsareprovided inAppendix D .Lemma 1 impliestheasymptoticpropertiesofthe directestimators

δ

ˆ_DIR viathedeltamethod.Next,wefocusonthe CDestimators

δ

ˆ_CD withthefollowinglemma.

Lemma 2. The first-order derivative of the log-confidence density function logh

(

ζ

)

≡ logh

(

δ

)

with respect to

ζ

, is asymptotically equivalenttothescorefunctions

(

δ

)

₌

∂

logL

(

δ

)

_/

∂

δ

.

According to Lemma 2 , the CD estimator ˆ

δ

_CD and direct esti-mator

δ

ˆ_DIR shareexactlyidenticalasymptoticproperties.Therefore, we canintroduceTheorem 1 inanalogytoLemma 1 immediately followingLemma 2 .

Theorem1. Asn→∞,theCDestimator

δ

ˆ_CD isconsistentand nor-mallydistributed.Specifically,

n1 /2



δ

ˆ_CD −

δ



−→d MN



0,



J

(

δ

)

T˜I

(

ζ

)

J

(

δ

)



−1



, (22) where˜I

(

ζ

)

=V−1 andJ

(

δ

)

=

∂

M

(

δ

)

/

∂

δ

.

The statements inthe theoremshow that theCD approachis asymptotically as efficient as the direct estimation approach. To quantify the uncertainty in the CD estimators, we put forward Corollary 1 basedonLemma 2 andTheorem 1 .

Corollary1. Thecovariance matrixofn1 /2



δ

ˆ_CD −

δ



canbe consis-tentlyestimatedbyn

ˆ_CD ,where

ˆ

CD =

−

∂

2

∂

δ∂

δ

T logh



ˆ

δ

CD



−1 , (23)

Remark. Ontheonehand,theconfidencedensitybasedmethods canaddresstheaforementionedproblemtohierarchicallyestimate parameterswithout imposingdifficultiesinthe EMalgorithm.On the other hand,aswasadvised in Liu et al. (2015) , the informa-tionfromindependentstudiescan bewellcombinedviathe con-fidence density.Inthe problemwe havebeenfocusing on inthe paper,theproposedhierarchicalanalysiscan beappliedto degra-dationtestsunderdifferentaccelerationfunctionsandfinallyyield integratedresultsfortheparametersofinterest,whichcould bea subset or transformationof parameters that are already involved inthosetests.

Inlightofthepreviousanalyses onthehierarchicalestimation, we have justified the efficiencyof the CDapproach. The estima-tion of

δ

can be readily obtained via (19) . Due to the nonlinear andnon-additivepropertiesofthePeckmodel,itisnotintuitiveto comparetheeffectsofenvironmentalfactorsundertheusage envi-ronments,whichcouldbeofgreatinteresttoreliabilityengineers anddecisionmakers.Thefollowingpropositionisgivenunderthe Peckmodeltocomparetheeffectsbroughtbyasingle-unitchange intemperatureandrelativehumidity.

Proposition1. (IntuitivecomparisonofeffectsunderthePeckmodel)

The baseline degradation rate at the reference environment

f_acc

(

x_{ref ;}

δ

)

under parameters

δ

₌

(

δ

_{0 ,}Ea,

δ

1

)

is more sensitive to

temperatureif Ea 11605· 1 TK_ref 2 · RH_ref

δ

1 >1,

andismoresensitivetorelativehumidityotherwise.

Thepropositioncanbestraightforwardlyjustifiedbytakingthe first-orderderivativeon facc

(

x;

δ

)

withrespecttoTK andRH,thus theproofisomitted.Thepropositionwillbeusedforillustrationin Section 6.2 .Ifaccelerationmodels other thanthePeck modelare used,similarstatementscanalsobeentertainedforeffect compar-ison.

4. Uncertaintyquantificationoftheestimatedparameters

Toquantify theuncertainties inthe parameterestimators is a vital taskto enableandjustify theadoption of theestimators in knowledgecreationanddecisionmaking.Comparedtopoint esti-mation,interval estimationis usuallypreferable inrealproblems. Inthissection,wediscussthelarge-samplebasedmethodto con-structconfidenceintervalsforestimatedparameters.

Theassumptionsoflarge-sampleapproximationarecommonly utilizedtoprovideasymptoticcovariancematrixofparameter esti-matorsfromwhichconfidenceintervalsareobtained.Forcomplete datasets, a routine practice is to compute the Fisher information fromthe log-likelihooddirectly, theelementsin theFisher infor-mationmatrixI

(

θ

)

aregivenby



I

(

θ

)

_{i j}=E

−

∂

2

∂

θ

i

∂

θ

j logL

(

θ|

y−,z

)

|

θ



. (24)

Asymptotically,

θ

ˆfollowsaMNdistribution,i.e.,

θ

ˆ∼ N

(

θ

,



I

(

θ

)

−1

)

. Since

θ

isalsounknown,wecanalternativelyemploytheobserved informationI

(

θ

ˆ

)

tocomputetheasymptoticcovariancematrixof

θ

. Forincompletedatasetswithhiddenobservationsasdiscussed inSection 3.4 , we adopttheapproach fromOakes (1999) ,where the Fisher information can be directly calculated via function Q.

Table 1

Experimental design πfor the test under treatment 1–9 (numbered in parenthesis). RH level Temperature level Total

303.15 Kelvin 318.15 Kelvin 333.15 Kelvin (30 degree Celsius) (45 degree Celsius) (60 degree Celsius) 60 (1) 16/49 (4) 8/49 (7) 4/49 4/7 75 (2) 8/49 (5) 4/49 (8) 2/49 2/7 90 (3) 4/49 (6) 2/49 (9) 1/49 1/7

Total 4/7 2/7 1/7 1

Table 2

Parameters as input to simulation study.

Model Parameter Value

Peck model δ0 1

Ea 2 × 10 7

δ1 3

Increased degradation rate ω 1 Diffusion parameter σ2 _0.3 Degradation reduction (known) u0 2

v0 3

Accordingly, theobservedFisher informationcan becomputedby

I



θ

ˆ



=−

⎡

⎣

∂

2 _Q



_θ|

_θ

ˆ



∂

θ∂

θ

T +

∂

2 _Q



_θ|

_θ

ˆ



∂

θ∂

_θ

ˆT

⎤

⎦

θ= _θˆ . (25)

Note that thesecond terminthe r.h.s. oftheequation is viewed asthemissinginformationduetothe absenceofz.Likewise,the asymptotic confidence intervals for unknown parameters can be constructed by theFisher information.The derivation of (25) re-quires manipulations based on Appendix B , and the analytical details of (25) are given in Appendix E . As a side note, for the parameter

σ

2 , thenormal approximated confidence intervalsare usually inappropriate. Alternatively,we build theconfidence inter-valsbasedonlog

σ

2 viathedeltamethod.Further,theconfidence intervalsoflog

σ

2 aretransformedbyanexponentialoperationto quantify the uncertainty in

σ

2 . More approaches foruncertainty quantification based on EM algorithm can be found in Louis (1982) andMeng and Rubin (1991) .

5. Simulationstudy

Tofacilitatethe simulation,we needtospecifyan experimen-taldesignofthedegradationtest.BylettingN=iNibethetotal

testunitsand

π

_i=N_{i/N be}theproportionoftestunitsallocatedto stressi,wesupposethat

π

=

(

π

1 ,...,

π

M

)

isapre-specified

exper-imentaldesignforthesimulationstudy.Withoutlossofgenerality, thetestisassumedtoallowforthechangeintemperatureand hu-midityunderthePeckmodelintroducedinAppendix A . Addition-ally,threelevelsforeachfactorarespecifiedandwefollowa4:2:1 allocationruleforeachfactor(Meeker & Escobar, 1998; Meeker & Hahn, 1977 ).Tobespecific,thetestplanallocatesmoretestunits tolowerstresslevelstoavoidoverwhelmingextrapolation.The de-tailedplanisshowninTable 1 .

Parameters are setas showninTable 2 forthe purposeof il-lustration. As a side note, the referencetemperature andrelative humidityarefixedatTK_ref =293.15Kelvin(20degreeCelsius)and

RH_{ref =}50%_,respectively.

Toexploretheeffectofsamplesizes,wewillshowresults un-derN=49, 98and147inthesimulationstudies.Testunitsare as-sumedtobeinspectedfor3times.Thefollowingfoursub-studies constitutethissectiontoexploretheeffectivenessoftheproposed model.

5.1. EstimatesfromtheEMalgorithmandconfidenceintervals

First,point estimatesareobtainedfromthe EMalgorithm un-der1000simulationreplicatesforeachsamplesizeofinterest.The convergencecriterion isset as

ε

=0.001(1‰) forall parameters. In Table 3 , the mean bias and root mean squared error (RMSE) are shown under N₌49_, 98 and 147. By observing the results, we can implythat the EMalgorithm can accurately estimate the unknown parameters,withratherlow meanbias andRMSEeven under moderate sample sizes. Moreover, the accuracy of the es-timation enhanceswith the increase in sample size. Specifically, compared to other parameters, the estimation accuracy of

σ

2 is relatively low but drastically improves over the sample size. A possible reasonbehindthisis that the estimationof

σ

2 involves bothE

(

zi j(k−1)

|

y−,

θ

(n)

)

andE

(

z2 i j(k−1)

|

y−,

θ

(n)

)

asindicatedin(15) ,

wheremore uncertaintyofhidden variablesare broughtinto the estimators.

Further,withthe extantpoint estimates,the confidence inter-valsareconstructedviathemethodproposedinSection 4 . Specif-ically,the

(

1−

α

)

× 100%confidenceintervalisgivenby

ˆ

θ

i± zα/2



[I

(

θ

ˆ

)

]ii



1 /2 ,

where[·]ii denotes theithdiagonalelement ofa matrixandzα/2 is the 1−

α

/2 quantile of the standard normal distribution. In Table 4 , the coverage probabilities and the average lengths are listedunderthe simulateddatasetsunderthreesamplesizes. We computethe95%confidenceintervalsunderlarge-sample approx-imation.Asthesamplesizeincreases,thecoverageprobability be-comescloserto0.95andtheaveragelengthisshortened.Asseen, evenunderasampleof49,theconfidenceintervalsperformwell andformostparametersover90%ofthemcancoverthetrue val-ues. Again, influenced by the relatively large bias, the coverage probability of

σ

2 is moderately lower under small samplesizes. Recallthat weusethelog

σ

2 toconstructconfidenceintervalsfor

σ

2 _{. The trick isproven to benefitthe performance. A supporting} exampleisthatunderN=49,weobtainacoverageprobabilityof 0.841comparingto0.750where

σ

2 isdirectlyusedtoquantifythe uncertainty.

5.2. Hierarchicalanalysis

We now consider the estimation of

δ

in the Peck model by meansofthe proposed hierarchicalanalysis.Likewise, the perfor-manceofpointestimationanduncertaintyquantificationarelisted inTable 5 andTable 6 ,respectively.

Theresults implygoodperformanceswithlow meanbiasand RMSE for the point estimators as well as coverage probabilities that are closeenough to 0.95. Itis worth notingthat the hierar-chicalanalysisconsumeslimitedcomputational efforts.For exam-ple,pointestimation togetherwithuncertaintyquantification un-der N=147 only takes less than 1 seconds on a single Intel i5 core. If hierarchical analysis is not employed, the complexity of

facc

(

·

)

willhinderthederivationofanalyticalresultsintheEM al-gorithm,which couldintroduceenormouscomputational burdens totheproblem.

5.3. Sensitivityanalysiswithrespecttothedegradationreduction effect

AsdiscussedinSection 3 ,thechoiceofu₀ and

v

0 leansonthe experienceofengineers.Misspecificationofthedistributionforz₀

mayoccurandleadto higherbiasofparameterestimation.Here, by assuming thetrue valuesofu0 and

v

0 to be 2and3, respec-tively,we changethe assumedvalues toexplore the influenceof misspecificationunder samplesize N=49. Note that E

(

z₀

)

=0.4

Table 3

Mean bias and RMSE of unknown parameters under N = 49 , 98 and 147. Parameter θ True

value

N = 49 N = 98 N = 147 Bias RMSE Bias RMSE Bias RMSE

facc( x 1) 2.104 0.011 0.264 −0.001 0.190 0.008 0.155 facc( x 2) 4.110 0.010 0.283 0.005 0.202 0.004 0.162 facc( x 3) 7.101 −0.005 0.344 0.004 0.250 0.006 0.200 facc( x 4) 2.751 −0.001 0.288 0.003 0.206 0.001 0.173 facc( x 5) 5.373 0.003 0.344 0.006 0.248 0.009 0.196 facc( x 6) 9.284 −0.002 0.439 0.007 0.313 0.006 0.266 facc( x 7) 3.511 0.013 0.343 0.003 0.247 0.011 0.204 facc( x 8) 6.857 0.006 0.444 −0.006 0.316 0.012 0.256 facc( x 9) 11.849 −0.001 0.647 −0.012 0.425 −0.004 0.358 ω 1.000 −0.009 0.178 −0.003 0.130 −0.008 0.106 σ2 _{0.300 −0.044 0.075 −0.020 0.049 −0.009 0.037} Table 4

Coverage probability and average length of the 95% confidence intervals of θunder N = 49 , 98 and 147. Parameter θ N = 49 N = 98 N = 147

Cov. prob. Avg. len. Cov. prob. Avg. len. Cov. prob. Avg. len.

facc( x 1) 0.916 0.944 0.922 0.693 0.929 0.572 facc( x 2) 0.915 1.032 0.930 0.758 0.944 0.627 facc( x 3) 0.919 1.267 0.938 0.924 0.954 0.766 facc( x 4) 0.918 1.044 0.933 0.766 0.948 0.633 facc( x 5) 0.926 1.254 0.933 0.922 0.948 0.762 facc( x 6) 0.927 1.665 0.951 1.215 0.949 1.002 facc( x 7) 0.933 1.249 0.933 0.917 0.948 0.757 facc( x 8) 0.927 1.647 0.946 1.201 0.941 0.993 facc( x 9) 0.933 2.380 0.949 1.676 0.946 1.367 ω 0.913 0.638 0.925 0.465 0.937 0.383 σ2 _{0.841 0.220 0.901 0.167 0.936 0.140} Table 5

Mean bias and RMSE of estimated δunder N = 49 , 98 and 147.

Parameter δ True value N = 49 N = 98 N = 147 Bias RMSE Bias RMSE Bias RMSE

δ0 1 0.050 0.176 0.025 0.112 0.017 0.092 Ea ( ×10 7 ) 2 −0.039 0.2481 −0.020 0.160 −0.012 0.138 δ1 3 −0.081 0.315 −0.031 0.187 −0.027 0.171 Table 6

Coverage probability and average length of the 95% confidence intervals of δunder N = 49 , 98 and 147. Parameter δ N = 49 N = 98 N = 147

Cov. prob. Avg. len. Cov. prob. Avg. len. Cov. prob. Avg. len.

δ0 0.933 0.542 0.944 0.389 0.957 0.319 Ea ( ×10 7 ) 0.944 0.692 0.941 0.508 0.958 0.419 δ1 0.905 0.686 0.938 0.506 0.948 0.417

andvar

(

z₀

)

=0.04 holdundertruevalues.In Table 7 ,we listthe biasandRMSEunderfoursettingsofmisspecfiedvaluesofu₀ and

v

0 asfollows:

• HighMean:u₀ =3,

v

0 =2,i.e.,E

(

z0

)

=0.6,var

(

z0

)

=0.04; • LowMean:u₀ =0.6,

v

0 =2.4,i.e.,E

(

z0

)

=0.2,var

(

z0

)

=0.04; • High Variance: u0=1.2,

v

0=1.8, i.e., E

(

z0

)

=0.4,var

(

z0

)

=

0.06;

• Low Variance: u₀ =9.2,

v

0 =13.8, i.e., E

(

z0

)

=0.4,var

(

z0

)

= 0.01.

As seen from the result, the misspecification of mean of z₀

brings considerablebiastotheestimatorsof facc

(

xi

)

and

ω

,while

theestimationof

σ

2 suffersmorewhenthevarianceis misspeci-fied.ForextrapolatinganalysisbasedontheADTdata, facc

(

xi

)

and

ω

play more important roles.For thispurpose, engineersshould focusonevaluatingthemeaneffectofdegradationreductionfora

betterestimationaccuracy.Toexploretheinfluence of misspecifi-cationontheestimationofparametersinthePeckmodel,wecarry outthehierarchicalmethodstoestimate

δ

andshowthemeanbias andRMSEinTable 8 .The estimated

δ

suffersa considerablebias whenthe meanofz₀ is misspecified,whilethe influenceof mis-specifiedvarianceexertsrelativelysmallerinfluenceonthe estima-tionaccuracy.

6. Casestudy

6.1. DatafromSchneiderElectric

Schneider Electric conducted a degradation test for a type of electrical distribution device.A total of104 test unitsunderwent the test underfourdifferent settings oftemperatureand humid-ity.Specifically,twolevelsoftemperature(313.15Kelvinand333.15

Table 7

Mean bias and RMSE of unknown parameters θunder misspecified u 0 and v0 .

Parameter θ High mean Low mean High variance Low variance Bias RMSE Bias RMSE Bias RMSE Bias RMSE

facc( x 1) −0.662 0.712 0.644 0.704 0.065 0.273 −0.008 0.287 facc( x 2) −0.566 0.642 0.539 0.625 0.052 0.305 −0.008 0.310 facc( x 3) −0.531 0.640 0.465 0.595 0.032 0.352 0.009 0.397 facc( x 4) −0.614 0.676 0.600 0.674 0.063 0.294 −0.010 0.312 facc( x 5) −0.535 0.639 0.487 0.608 0.039 0.346 −0.008 0.366 facc( x 6) −0.498 0.700 0.458 0.684 0.041 0.486 0.007 0.556 facc( x 7) −0.574 0.671 0.550 0.662 0.054 0.352 −0.006 0.377 facc( x 8) −0.529 0.699 0.441 0.654 0.025 0.450 −0.005 0.528 facc( x 9) −0.471 0.835 0.394 0.811 −0.006 0.647 0.001 0.762 ω 0.755 0.776 −0.733 0.761 −0.058 0.192 0.008 0.191 σ2 _{−0.025 0.071 0.005 0.071 −0.078 0.096 0.213 0.239} Table 8

Mean bias and RMSE of unknown parameters δunder misspecified u 0 and v0 .

Parameter δ High mean Low mean High variance High variance Bias RMSE Bias RMSE Bias RMSE Bias RMSE

δ0 −0.181 0.221 0.178 0.202 0.096 0.176 0.035 0.162 Ea ( 10 7 ) 0.152 0.262 −0.146 0.233 −0.052 0.187 −0.009 0.221 δ1 0.252 0.328 −0.234 0.321 −0.105 0.210 −0.023 0.208

0

0.5

1

1.5

Time

0.2

0.4

0.6

0.8

1

1.2

Degradation level

TK=313.15, RH=60

0

0.5

1

1.5

Time

0.2

0.4

0.6

0.8

1

1.2

TK=313.15, RH=95

0

0.5

1

1.5

Time

0.2

0.4

0.6

0.8

1

1.2

Degradation level

TK=333.15, RH=60

0

0.5

1

1.5

Time

0.2

0.4

0.6

0.8

1

1.2

TK=333.15, RH=95

Fig. 2. Degradation test data under 4 stress levels.

Kelvin)andtwolevelsofrelativehumidity(60%and90%)are con-sidered. Notethat 313.15 Kelvin and 333.15 Kelvin are equivalent to 40degree Celsiusand60degree Celsius, respectively.The test chambersprovide generallyhigherstresses thantheusage condi-tion of thedevice, thus the test canbe regarded asan ADT. The test dataundereachstress conditionsareplottedinFig. 2 .Itcan

beobservedthattheobservationepochs(

τ

i jk)andnumberof

ob-servations(Oi jk)varyacrossthetest unitsinthetest.The

specifi-cationsofthetestareshowninTable 9 .Basedonthediscussions withengineersfromSchneider,weproposetomodelthe degrada-tionreductioneffectbyabetadistributionwithparametersu₀ =1 and

v

0 =3.