Automobile Battery Remaining Useful Life Prediction
Jianguo Wu
May 12, 2014
1
INTRODUCTION
Nowadays, the hybrid and full electric automobiles are becoming more and more popular with the
development of battery technologies and increasing demand of sustainable energy. For these
types of automobiles, the reliability and life of the battery is crucial to guarantee the overall
functional capabilities and performance. Failure of such a unit can be catastrophic. Therefore, the
battery must be well-maintained for the overall system reliability and to improve end user
satisfaction. In traditional maintenance strategies, the battery is either repaired after its failure
(run-to-failure maintenance) or scheduled for time-based preventive maintenance (planned
maintenance).
The rapid development of sensing technology makes it possible to monitor the status of a battery
in the field for diagnostics and prognostics, enabling a more cost-effective strategy called
condition-based-maintenance. This study is focused on the remaining useful life prediction of the
automobile battery using the battery time-to-failure data and resistance signal, which is strongly
associated with failure of the battery and contains rich information about its health status.
Specifically, the goal of this study is to address the following three questions:
1) How does the resistance of the battery change with time?
2) What is the relationship between battery health status and resistance?
3) How to predict the remaining useful life?
The rest of this report is organized as follows. Section 2 introduces the experiment design and
data collection. Section 3 presents the statistical analysis. The conclusion is given in Section 4.
1
2
DATA DESCRIPTION
Since the real aging process of automotive batteries is typically very slow and it takes several
years for a battery to fail, the data is obtained from an accelerated aging test performed by
General Motors based on the aging cycle defined in SAE J2801. In this test, a battery is
considered as dead when it fails to crank the engine. The dataset contains total 14 lead acid
batteries of two different types: 8 batteries of type JBI and 6 of type JCI. The failure time of each
battery and the resistance in each week are recorded. Figure 1 shows the resistance of each
battery as a function of time and the failure time.
5
JBI
10
JCI
Resistance (mΩ)
8
7
6
5
4
3
5
10
Time (week)
Figure 1. Resistance of 14 lead acid batteries. All batteries are failed and there are no right
censored data.
3
3.1
STATISTICAL ANALYSIS
Change of Resistance
From Figure 1 we can see that the resistance increases with time and the increasing speed
becomes higher and higher. The second order polynomial regression is applied to study the
change of the battery resistance in the aging process.
resistance = ∗ + ∗ + 2
Since there are only 14 batteries, the influence of battery type is not included in the model. In the
remaining useful life prediction the battery type will be added to differentiate these two types of
batteries. Figure 2 shows the 95% confidence intervals on the estimated coefficients of the
polynomial model for each battery. As we can see, these coefficients are different among
different batteries. Therefore the linear mixed effects model are used where the coefficients , and follow multivariate normal distribution: , , ~ , , , Σ.
time
battid
JCI6
JCI5
JCI4
JCI3
JCI2
JCI1
JBI7
JBI2
JBI16
JBI15
JBI14
JBI13
JBI10
JBI1
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
-0.4
|
|
|
|
|
-0.6
|
|
|
|
-0.2
0.0
0.2
0.4
(Intercept)
JCI6
JCI5
JCI4
JCI3
JCI2
JCI1
JBI7
JBI2
JBI16
JBI15
JBI14
JBI13
JBI10
JBI1
|
|
|
|
I(time^2)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
3.0
|
|
|
|
|
4.0
4.5
|
|
5.0
|
|
|
|
|
|
|
|
3.5
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
-0.02 0.00 0.02 0.04 0.06 0.08
Figure 2. Comparison of 95% confidence intervals on the coefficients of second order polynomial
model
Figure 3 shows the comparison of the fitted resistances of the linear mixed effects model with the
observed resistances. Figure 4 shows the residuals with the fitted resistance. From these two
figures we can see that for most batteries, the second order polynomial model is very accurate to
describe the change of the resistance. The estimated coefficients and confidence intervals are
given in Table 1.
3
JBI1
JBI10
JBI13
JBI14
JBI15
JBI16
JBI2
JBI7
JCI1
JCI2
JCI3
JCI4
JCI5
JCI6
8
7
6
5
4
3
8
r
7
6
5
4
3
8
7
6
5
4
3
5
10
5
10
5
10
5
10
obstime
Figure 3. Fitted resistance (line) of the linear mixed effects model VS. the observations (dot)
resid(., type = "pearson")
3 4 5 6 7
0.6
0.2
-0.2
-0.6
0.6
0.2
-0.2
-0.6
JCI5
JCI6
JCI1
JCI2
JBI15
JBI1
0.6
0.2
-0.2
-0.6
JCI3
JCI4
JBI16
JBI2
JBI7
JBI10
JBI13
JBI14
3 4 5 6 7
3 4 5
0.6
0.2
-0.2
-0.6
6 7
fitted(.)
Figure 4. fitting residuals VS. fitted resistance
Table 1. Estimated coefficients and 95% confidence intervals
Coefficients
Estimated mean 0.011
0.001
3.781
4
95% CI
[0.0038, 0.019]
[-0.0908, 0.093]
[3.44,
4.11]
3.2
The Relationship between Resistance and Health Condition
The Cox survival model is used to study the relationship between battery resistance and health
condition and to predict the remaining life. The hazard rate is expressed as
ℎ = ℎ exp[# + $ % 0 + $' % − % 0]
where ℎ is the hazard rate for battery * at time , ℎ is the nonparametric baseline hazard
function, # is the intercept for the battery type of battery *, % 0 is the fitted initial resistance for
battery * and % is the fitted resistance at time for battery *. Note here we use the fitted
resistance obtained in Section 3.1 to fit the Cox model. The reason is that in the linear mixed
effects model in Section 3.1 we assume that the observed resistance have measurement errors or
noises.
Table 2. Estimated coefficients and 95% confidence intervals for the Cox model
Estimated
-0.055
-1.20
1.61
P-value
0.94
0.27
0.053
95% C.I.
[-1.5, 1.39]
[-3.32, 0.92]
[-0.02, 3.24]
30
20
10
0
Baseline Hazard Rate
40
Coefficient
# (type JBI)
$
$'
2
4
6
8
10
12
14
Time
Figure 5. Fitted baseline hazard rate function
The fitted coefficients for the Cox model are given in Table 2. Figure 5 shows the fitted baseline
hazard function ℎ . As we can see, the battery type has almost no effect on the hazard rate.
The estimated coefficient $ for the initial resistance is negative. From Figure 1 we can see that
there are three batteries of type JBI with high initial resistance and long battery life, which may
5
cause the estimated coefficient $ to be negative. The estimated $' is positive and significant,
indicating that the increase of resistance is positively correlated with hazard rate of the battery
failure.
3.3
Remaining Useful Life Prediction
In section 3.1, we estimated the mean and covariance matrix for , , + = , , ~, Σ
0.00015 −0.0016 0.007
where , = 0.011, 0.001, 3.781 and 2 = 3−0.0016
0.021
−0.0628.
0.007
−0.062
0.34
To predict the remaining useful life of a battery (indexed as 9) at time ∗ , we can calculate the
survival function as follows:
S (t | T > t * , w p , rph* )
t
= ∫ exp{− ∫ hˆ0 (u ) exp[ w p + βˆ0 x p (0) + βˆ1 ( x p (u ) − x p (0))]du}p (θ p | rph* )dθ p
t*
where :; = <
; , ; , ; = is the resistance coefficients for battery 9. From the formula above we
can see that the future resistance needs to be predicted to calculate the survival function. To
predict the future resistance, the estimated and Σ are used as the population prior and given the
resistance of battery 9 up to the current time ∗ , the posterior distribution of :; can be updated
sequentially (Figure 6).The survival function can be estimated using Monte Carlo simulation.
Figure 6. Sequential prediction of the resistance
6
Figure 7. Estimated survival functions at ∗ = 5 and ∗ = 10 for battery JBI15 as shown in
Figure 3 (true life is 12 weeks)
Figure 7 shows the estimated survival function for battery JBI15 at time ∗ = 5 and ∗ = 10
respectively. As we can see, as we observe more resistance data, estimation interval gets narrower.
The expected remaining useful life can be calculated by
>?@A
∗
= >B − ∗ |B
>
∗
I
= E FG| ∗ HG
J∗
Table 3 shows the comparison between the expected RUL and the true RUL for part of the
batteries. As we can see, the expected value is reasonable, though not very close to the true value.
Table 3. Comparison between the expected RUL and true RUL
Battery ID
JCI1
JCI4
JCI5
JCI6
JBI1
JBI7
∗ = 4
ERUL
True value
6.51
6
6.5
7
6.52
7
6.50
9
6.53
8
6.53
9
7
∗ = 7
ERUL
True value
6.04
3
4.97
4
6.32
4
6.04
6
6.26
5
6.32
6
4
CONCLUSION
In this report, the linear mixed effects model is used to study the change of battery resistance
along time in the aging process. The exploratory analysis shows that the resistance increases
quadratically in the aging process. Therefore the second order polynomial with random
coefficients is used to model the resistance change.
The fitted true resistances are used to fit the survival Cox model for the time-to-failure events.
The fitted results indicate that the battery type has almost no influence on the battery life. The
increase of resistance could increase the failure hazard rate.
To predict the remaining useful life, the estimated mean and variance of the coefficients for the
resistance are used as prior and the posterior distribution of the resistance of the target battery is
sequentially updated. The survival function and the expected remaining useful life could be
calculated using the fitted Cox model. The estimated expected remaining useful life is reasonable.
8
Appendix
R code
library(survival)
library(nlme)
library(lattice)
library(ggplot2)
library(lme4)
library(mgcv)
library(glmmML) # Gauss-Hermite quadrature
library(AlgDesign) # factorial design
library(statmod) #Gauss-Legendre quadrature
rm(list=ls())
setwd("D:/Jianguo/study/2014Spring/Stat998-Consulting/project3");
source("R_functions_new_new.R")
Data = read.table("DATA.txt",header=TRUE)
xyplot(r~time|type,groups=battid,type=c("l","p"),
lty=c(1,3),lwd=c(2,2),
data=Data,layout=c(2,1),
xlab="Time (week)",ylab="Resistance (mΩ)",cex.lab=1.5)
listFit=nlme::lmList(r~(I(time^2)+time+1)|battid,Data)
plot(intervals(listFit))
lmer.fit=lmer(r~I(time^2)+time+1+(I(time^2)+time+1|battid),data=Data,REML=T)
summary(lmer.fit)
confint(lmer.fit,method="Wald")
coef(lmer.fit)
plot(lmer.fit,resid(., type = "pearson") ~ fitted(.)|as.factor(battid))
plot(lmer.fit,r~ fitted(.)|as.factor(type))
confint(lmer.fit,method="Wald")
################################################
TEMP = read.table("DATA.txt",header=TRUE)
tr_temp = length(TEMP[,1])
TYPE = c("JCI1","JCI2","JCI3","JCI4","JCI5","JCI6","JBI1","JBI2","JBI7","JBI10","JBI13","JBI14","JBI15","JBI16")
b1 = TEMP[TEMP$battid==TYPE[1],]
b2 = TEMP[TEMP$battid==TYPE[2],]
b3 = TEMP[TEMP$battid==TYPE[3],]
b4 = TEMP[TEMP$battid==TYPE[4],]
b5 = TEMP[TEMP$battid==TYPE[5],]
b6 = TEMP[TEMP$battid==TYPE[6],]
b7 = TEMP[TEMP$battid==TYPE[7],]
b8 = TEMP[TEMP$battid==TYPE[8],]
b9 = TEMP[TEMP$battid==TYPE[9],]
b10 = TEMP[TEMP$battid==TYPE[10],]
b11 = TEMP[TEMP$battid==TYPE[11],]
9
b12 = TEMP[TEMP$battid==TYPE[12],]
b13 = TEMP[TEMP$battid==TYPE[13],]
b14 = TEMP[TEMP$battid==TYPE[14],]
l
=
c(length(b1[,1]),length(b2[,1]),length(b3[,1]),length(b4[,1]),length(b5[,1]),length(b6[,1]),length(b7[,1]),length(b8[,1]),l
ength(b9[,1]),length(b10[,1]),length(b11[,1]),length(b12[,1]),length(b13[,1]),length(b14[,1]))
ol=l
odeath
=
c(rep(l[1]+1,l[1]),rep(l[2]+1,l[2]),rep(l[3]+1,l[3]),rep(l[4]+1,l[4]),rep(l[5]+1,l[5]),rep(l[6]+1,l[6]),rep(l[7]+1,l[7]),rep(l[
8]+1,l[8]),rep(l[9]+1,l[9]),rep(l[10]+1,l[10]),rep(l[11]+1,l[11]),rep(l[12]+1,l[12]),rep(l[13]+1,l[13]),rep(l[14]+1,l[14]))
odata = rbind(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14)
D = c(11,9,9,12,12,14,13,11,14,10,15,14,13,7)
data = odata
V = odeath
tr = length(data[,1])
rNA = rep(NA,tr)
bdata = data.frame(battid=data$battid, x=data$x, r=data$r, status=rep(1,tr), obstime=data$time, eventtime=V[1:tr],
start=data$time, stop=data$time+1, died=rNA, rini = rNA, rfit = rNA, rinifit = rNA)
for (a in 1:14)
{
ts = sum(l[0:(a-1)])
bdata[(ts+1):(ts+l[a]),9] = c(rep(0,l[a]-1),1)
bdata[(ts+1):(ts+l[a]),10] = rep(bdata$r[ts+1],l[a])
}
bdata_surv=Surv(bdata$start,bdata$stop,bdata$died)
fitlme=
lme(r~obstime+I(obstime^2),random=~obstime+I(obstime^2)|battid,data=bdata,control=lmeControl(returnObject=TR
UE))
summary(fitlme)
sigma20 = (fitlme$sigma)^2
MUB0 = fitlme$coefficients$fixed
SIGMAB0 = var(fitlme$coefficients$random$battid)
##### estimation of survival intial parameters #####
# add fitted resistance values into bdata (step 2 of the two-step approach)
bdata$rfit = fitted(fitlme)
ggplot(bdata, aes(x=obstime,y=r)) + facet_wrap(~battid, nrow =3, ncol = 5)+geom_point() +
geom_line(aes(x=obstime,y=rfit),colour="red") # second layer
for (i in 1:14)
{
ir = ceiling(D-1)[i] # number of rows for the ith battery
ts = sum(ceiling(D-1)[0:(i-1)])
bdata$rinifit[(ts+1):(ts+l[i])] = rep(bdata$rfit[(ts+1)],l[i]) # fitted initial resistance
}
fitsurv = coxph(bdata_surv~as.factor(x)+rinifit+I(rfit-rinifit),data=bdata)
#summary(fitsurv)
gam0 = coef(fitsurv)[1]
10
beta00 = coef(fitsurv)[2]
beta10 = coef(fitsurv)[3]
# spline smoothed cumulative baseline hazard function (take derivative to get h0t)
H0t_est = basehaz(fitsurv,centered=FALSE)
plot(H0t_est$time,H0t_est$hazard,xlab="Time",ylab="Baseline Hazard Rate",type='b')
H0t_sm = smooth.spline(H0t_est$time,H0t_est$hazard,nknots=5)
summary(fitsurv)
confint(fitsurv)
test=cox.zph(fitsurv)
par(mfrow=c(1,2))
plot(test[2])
plot(test[3])
# function definitions for Weibull distribution
# default par: lam=0.001, alpha=1.05, gam=2.5, beta=c(0.15,0.50)
lambda = 0.001
alpha = 1.05
gam = 0.20 #coefficient to a fixed variable x (e.g., manufacturer A/B)
beta0 = 0.15 # coefficient to initial resistance
beta1 = 0.50 # coefficient to resistance increase
Zt = function(t,B) B[1]+t*B[2]+t*t*B[3] # resistance degradation function,t/month,mohm
h0t = function(t) lambda*alpha*(t^(alpha-1)) #baseline hazard
H0t = function(t) lambda*(t^alpha)
ht = function(t,x,B) h0t(t)*exp(gam*x+beta0*B[1]+beta1*(Zt(t,B)-B[1])) #hazard function
S0t = function(t) exp(-lambda*(t^alpha)) # baseline survival
St_tp = function(t,x,B) exp(-integrate(f=ht, lower=1E-3, upper=t, x,B, subdivisions=1000)$value)
function
St = Vectorize(St_tp, "t") # (vectorized!) survival function
ft = function(t,x,B) St(t,x,B)*ht(t,x,B) # density function
#temp survival
### define the function for estimated conditional survival function
#PAR = list(H0t_sm=H0t_sm,BiMAT=BiMAT,beta00=beta00,
#
beta10=beta10, gam0=gam0, SIGMABp=SIGMABp, MUBp=MUBp, wtMAT=wtMAT, Nff=Nff)
St_cond_EST_tp = function(t, tstar, x, PAR)
{
Itp = 0
for (iNff in 1: PAR$Nff)
{
Bi = PAR$BiMAT[,iNff]
Itp = Itp + G(PAR$H0t_sm, tstar, t, PAR$Bi, PAR$beta00, PAR$beta10,PAR$gam0, x)*PAR$wtMAT[iNff]
}
return(Itp*(sqrt(2)^3)/((2*pi)^(3/2)))
# sqrt(2) is for change of variable: exp(-0.5*x^2)->exp(-x^2)
# ((2*pi)^(3/2)) is for the constant part in the normal distribution
}
St_cond_EST = Vectorize(St_cond_EST_tp, "t") # vectorised function
### define the function for the ALTERNATIVE estimated conditional survival function
Eht_tp = function (t,x,PAR)
# expected value of hazard rate function
{
h0t = predict(PAR$H0t_sm,t,deriv=1)$y
betaz = t(t(c(PAR$beta00,PAR$beta10*t,PAR$beta10*t*t))) # a column vector
return(h0t*exp(PAR$gam0*x + t(betaz)%*%PAR$MUBp + t(betaz)%*%PAR$SIGMABp%*%betaz/2))
11
}
Eht = Vectorize(Eht_tp, "t") # vectorised function
St_cond_EST2 = function(t,tstar,x,PAR)
{
exp(-sum(Eht(seq(tstar,t,length.out=1000),x,PAR)*(t-tstar)/999))
}
tstar = 10 # time of prediction
tp = 1:tstar
#select battery b12
Rp =b13$r[1:tstar]
Zp = matrix(0,length(tp),3)
Zp[,1] = rep(1,length(tp))
Zp[,2] = tp
Zp[,3] = tp^2
# posterior dist of B
SIGMABp = solve(solve(SIGMAB0)+t(Zp)%*%Zp/sigma20)
MUBp = SIGMABp%*%(t(Zp)%*%Rp/sigma20+solve(SIGMAB0)%*%MUB0)
# # Gauss_Hermit integration
nnode = 10
ghxw = ghq(nnode, FALSE) # a list of Gauss-Hermite nodes and weights
ZO = ghxw$zeros
Nff = nnode^3 # 3 is the number of random parameters in B
Ap=t(chol(SIGMABp))
FF = gen.factorial(c(nnode,nnode,nnode),center=FALSE) # a ff design for index
BiMAT = t(t(MUBp))%*%matrix(1,1,Nff) + Ap%*%rbind(ZO[FF[,1]],ZO[FF[,2]],ZO[FF[,3]])*sqrt(2)
# Bi calculated based on Z
# sqrt(2) is for change of variable: exp(-0.5*x^2)->exp(-x^2)
wtMAT = ghxw$weights[FF[,1]]*ghxw$weights[FF[,2]]*ghxw$weights[FF[,3]]
PAR = list(H0t_sm=H0t_sm,BiMAT=BiMAT,beta00=beta00, beta10=beta10, gam0=gam0, SIGMABp=SIGMABp,
MUBp=MUBp, wtMAT=wtMAT, Nff=Nff)
glxw = gauss.quad(n=10,kind="legendre")
mrl_EST = 0 #mean residual life (estimation)
mrl_EST2 = 0 #mean residual life (alternative estimation)
xp=1
G = function (H0t_sm, tstar, t, B, beta0, beta1,gam0,x)
{
tp = ht_sm(seq(tstar,t,length.out=1000),H0t_sm,B,beta0,beta1,gam0,x)
return(exp(-sum(tp*(t-tstar)/999)))
}
for (igl in 1:10)
{
mrl_EST
=
mrl_EST
+
tstar)/2+(200+tstar)/2,tstar,xp,PAR)*glxw$weights[igl]*(200-tstar)/2
mrl_EST2
=
mrl_EST2
+
tstar)/2+(200+tstar)/2,tstar,xp,PAR)*glxw$weights[igl]*(200-tstar)/2
}
12
St_cond_EST(glxw$nodes[igl]*(200St_cond_EST2(glxw$nodes[igl]*(200-
© Copyright 2024 ExpyDoc
