Abstract
In this work, a state-space battery model is derived mathematically to estimate the state-of-charge (SoC) of a battery system. Subsequently, Kalman filter (KF) is applied to predict the dynamical behavior of the battery model. Results show an accurate prediction as the accumulated error, in terms of root-mean-square (RMS), is a very small value. From this work, it is found that different sets of and values (KF’s parameters) can be applied for better performance and hence lower RMS error. This is the motivation for the application of a metaheuristic algorithm. Hence, the result is further improved by applying a genetic algorithm (GA) to tune and parameters of the KF. In an online application, a GA can be applied to obtain the optimal parameters of the KF before its application to a real plant (system). This simply means that the instantaneous response of the KF is not affected by the time consuming GA as this approach is applied only once to obtain the optimal parameters. The relevant workable MATLAB source codes are given in the appendix to ease future work and analysis in this area.
1. Introduction
Battery Management System (BMS) [1–3] comprises of mechanism that monitors and controls the normal operation of a battery system so as to ensure its safety while maintaining its State-of-Health (SoH). The BMS, in essence, measures the voltage, current, and temperature of each cell in a battery pack. These data are then analyzed by a management system that guarantees safe and reliable operations. A common example of an independent battery pack is portrayed in Figure 1. The battery is an essential component and should be accurately modeled in order to design an efficient management system [4]. Hence, a generic tool to describe the battery performance under a wide variety of conditions and applications is highly desirable [4]. As such, electrical modeling is necessary to provide such a tool that enables visualization of the processes occurring inside rechargeable batteries. These generic models are necessary for the development of battery management algorithms. These algorithms control the operation and maintain the performance of battery packs. In short, the ultimate aim of BMS is to prolong battery life, while ensuring reliable operation alongside many applications, especially in photovoltaic systems [5–7].
Battery modeling is performed in many ways depending on the types of battery. In general, the resulting battery model is a mathematical model comprising numerous mathematical descriptions [8]. Ultimately, battery models aim to determine the state-of-charge (SoC) of the battery system. However, the complexity of the nonlinear electrochemical processes has been a great barrier to modeling this dynamic process accurately. The accurate determination of the SoC will enable utilization of the battery for optimal performance and long-life and prevent irreversible physical damage to the battery [9]. The solution to the SoC via neural networks [10] and fuzzy logic [11] has been difficult and costly for online implementation due to the large computation required, causing the battery pack controller to be heavily loaded. This may however be a good alternative in the near future as the computational power of processing chips increase alongside their declining cost.
The state-estimation process usually leads to some state variables in a dynamical system. The SoC is a measure of a battery’s available power and thus it is important to calculate this value accurately from BMS by the cell voltage, temperature, and polarization effect caused by the electrolyte concentration gradient during high rate charging/discharging cycle [12]. Recently, the battery state-of-charge (SoC) is one of the significant elements attracting significant attention [13, 14]. By definition, SoC is the ratio of remaining capacity to the nominal capacity of the battery. Here, the remaining capacity is the number of ampere-hours (Ah) that can be extracted at normal operating temperature. The mathematical expression for the SoC is given in [13, 14], which is where is time, is battery SoC, in amphere-hours (Ah), is current flowing through the battery (passing through ), illustrated in Figure 2, and is nominal battery capacity. For charging, and for discharging, .
From this mathematical expression, it is noted that the SoC cannot be explicitly measured. In the literature, there is a myriad of methods dealing with predicting and estimating of SoC. The most popular of these methods are described in the following paragraphs.
At present, Coulomb-counting [15], also known as charge counting, or current integration is the most commonly used technique; it requires dynamic measurement of battery current. This is an open-loop method; however, it suffers from a reliance on the mathematical integration, and errors (noise, resolution, and rounding) are cumulative, which can lead to large SoC errors at the end of the integration process in (1). On the positive side, if an accurate current sensor is incorporated, the implementation will be much easier.
Another prominent SoC estimator is the well-known Kalman filter (KF), invented by Kalman in 1960. Although his popular work was published almost 54 years ago in [16], it remains as an important citation source in the literature. Readers who are new to this method can refer to an excellent KF tutorial by Faragher in [17]. The KF method is a well-known technique for dynamic state estimation of systems such as target tracking, navigation, and also battery systems [18, 19]. The state-of-the-art method provides recursive solution to linear filtering for state observation and prediction problems. The key advantage of the KF is that it accurately estimates states affected by external disturbances such as noises governed by Gaussian distribution. On the contrary, the disadvantage of KF is that it requires highly complex mathematical calculations. This can be realized by a state-space model, as shown in previous work by the author in [20, 21]. The modeling is a heavy duty task and is also presented in this work to ease the understanding of readers. As such, there exist some possibilities of divergence due to an inaccurate model and complex calculation. In the case of a slow processor, the calculation results may be delayed and exceed the sampling interval, thereby result in an inaccurate tracking.
Various artificial intelligence (AI) methods, mainly the neural networks and fuzzy logic, are being applied in the estimation of battery’s SoC [10, 22]. Neural networks are computationally expensive, which can overload the BMS and thus this approach, though theoretically feasible, may not be suitable for online implementation that requires instantaneous feedback and action. Also, neural networks require huge datasets in the time-consuming training process. Other techniques for SoC, include the sliding mode observer, are reported in [12].
In this work, a mathematical derivation leading to a state-space model is presented. The basic schematic model is adopted from [18, 20]. A thorough analysis in the form of state variables with the application of the Kalman filter is presented. The rest of the paper is organized as follows. The mathematical model is derived in Section 2, resulting in a state-space model. Further, in Section 3, the KF is applied to estimate the SoC of a battery system. This is followed by the tuning of KF’s parameter by adopting a metaheuristic approach, namely, a genetic algorithm in Section 4. Relevant results are presented in Section 5, and finally the conclusions are derived in Section 6.
2. Battery Modeling
Many model-based state-estimations have been proposed in [18, 20, 21, 23]. In [18, 21], the well-known Kalman filter [16] had been applied successfully for both state observation and prediction of the SoC. Work in [24] utilized manufacturers’ data in modeling the dynamic behavior of the battery. Several battery models exist from many works over the past few years. Each of these models varies in terms of its complexity and applications. In this work, a dynamical battery model is adopted, consisting of state variable equations, from [20, 21]. The schematic representation of this model is shown in Figure 2. In this model, there exists a bulk capacitor that acts as an energy storage component in the form of the charge, a capacitor that models the surface capacitance and diffusion effects within the cell , a terminal resistance , a surface resistance , and an end resistance . The voltages across both capacitors are denoted as and , respectively.
2.1. Mathematical Derivations of Battery Model
In this derivation, we aim to form a state-space model consisting of the state variables , , and . State variables are mathematical descriptions of the “state” of a dynamic system. In practice, the state of a system is used to determine its future behavior. Models that consist of a paired first-order differential equations are in the state-variable form.
Following the voltages and currents illustrated in Figure 2, the terminal voltage can be expressed as which is similar to By (2) and (3), and following straightforward algebraic manipulation, can be written as From Kirchoff’s current law, , Thus, substituting (5) into (4) yields By assuming a slow varying , that is, (from basic formula of ), and then substituting it into (6), and after rearranging results in By applying a similar derivation, the rate of change of the surface capacitor voltage , derived also from (2) and (3), gives Further, by assuming and , (7) and (8) can be written as respectively. Further, (9) and (10) can be combined to form a state variable relating voltages and and current flow , which is Next, the output voltage is derived from (2) and (3). By adding both equations, we obtain Then by substituting and into (11), it is further simplified as Taking the time derivative of the output voltage and assuming (this simply means that the change rate of terminal current can be ignored when implemented digitally), we obtain Substituting the formulae obtained in (9) and (10) into (13) results in Then, solving for from (12) we obtain and after substituting it into (14), it yields Finally, the complete state variable network is obtained by substituting (16) into (10), represented in matrix form as whereby constants , , and were derived previously but are hereby restated as This completes the initial derivation of the battery model.
2.2. Numerical Example
By substituting all capacitor and resistor values from Table 1 into (18), the following values are obtained:
By defining matrix , Again by substituting the values from Table 1 to calculate , , and , we obtain the value of as
and as
As such (17) can be rewritten as or numerically as
2.3. State-Space Modeling
Based on control theories, a lumped linear network can be written in the form where in this work, the state variable is Obviously, with Therefore, by restating the previous calculation values in (21) and (23), we should note that the values of , , , and are as follows: respectively. As such, with (32), the output is in fact This means that the output of the system is the open terminal voltage , as expected. Note that this is an important observation from this state-space modeling.
Further, the above state-space variables are transformed to a transfer function, . This is done by using function in MATLAB, which after factorization yields The plot of the unit step response for the gain in (35) is given in Figure 3. Basically, it shows that the open circuit terminal voltage in Figure 2 increases linearly during the charging operation in a very slow manner after transient behaviour for few seconds.
2.4. Observability of the RC Battery Model
In control theory, observability is the degree to which the internal states of a system can be predicted via its external outputs. As such, for an observable system, the behavior of the entire system can be predicted via the system’s outputs. On the other hand, if a system is not observable, the current values of some of its states cannot be estimated through the output signal. This means the controller does not know the states’ values. In theory, the observability of a system can be determined by constructing an observability matrix . and a system is observable if the row rank of is equal to (this is also known as a full rank matrix). The ultimate rationale of such a test is that if rows are linearly independent, then each of the states is viewable through linear combinations of the output .
Further, substituting and values from (30) and (32) into (36) yields Clearly, in this case is a full rank matrix, which concludes that this system is observable.
3. Kalman Filter for SoC Estimation
A continuous time-invariant linear system can be described in the state variable form as where is the input vector, is the state vector, is the output vector, is the time invariant dynamic matrix, is the time invariant input matrix, and is the time invariant measurement matrix.
If we assume that the applied input is constant during each sampling interval, a discrete-time equivalent model of the system will now be where whereby is the identity matrix and is the sampling period. As for this system, two independent process noises are present which are additive Gaussian noise, vector representing system disturbances and model inaccuracies and vector representing the effects of measurement noise.
Both and have a mean value of zero and the following covariance matrices: where denotes the expectation (or mean) operator and superscript means the transpose of the respective vectors. In usual case, and are normally set to a constant before simulation; in our case both are set to one (see Section 5). Further, a genetic algorithm (GA) is applied in order to obtain a better set of and values resulting in lower RMS error from the KF’s output.
By inclusion of these noises, the resulting system can now be described by which is illustrated in Figure 4,
3.1. Property of Kalman Filter
An important property of the KF is that it minimizes the sum-of-squared errors between the actual value and estimated states , given as To understand the operations of the KF, the meaning of the notation is crucial. Simply stated, it means that the estimate of at event takes into account all discrete events up to event . As such, (43) can include such information, now expanded as In recursive implementation of the KF, the current estimate , together with the input and measurement signals , is used for further estimating . This means that in this discrete system, the input for each sample step will be , , with respect to the output of , , . The recursive KF algorithm is obtained with the predictor and corrector stages.
3.2. KF Online Implementation
For the case of a battery, it is well understood that only the terminal quantities can be measured (terminal voltage and current ). Assuming that battery parameters are time-invariant quantities, the recursive KF algorithm is applied. By applying (40) into (30)–(32), we obtain the following values of updated matrices, with : Note that and remain similar to their previous values, as given in (31) and (32). By utilizing MATLAB’s control toolbox, the KF is placed in parallel to the state-space model, hereby represented by plant in Figure 5. The complete source code is given in the Appendix.
4. Genetic Algorithm
The genetic algorithm (GA), introduced by John Holland, is an approach based on biological evolution [25]. The algorithm is developed based on Charles Darwin’s theory of survival of the fittest. The GA has a very powerful encoding mechanism that enables the representation of a solution vector as either a real coded or binary string. Both encodings serve a different purpose in the context of different problem space. GAs are regarded as the global optimizer that often spot or locate the potential area or even accurately obtain the best solution, known as the global minimum [26–28]. Underneath this popular algorithm are the three operators that contribute to its success in performing optimization task.(i)Selection. In selection, offsprings with higher fitness have better chance for survival to the following generation in the evolutionary process. This basically is based on the theory of “survival of the fittest.”(ii)Crossover. The crossover increases and maintains the diversity of the entire population over the entire run. This is due to the fact that a population with higher diversity has the ability to explore a wider range of search space.(iii)Mutation. This enables chromosomes (potential solutions) to jump to a wider range than crossover. Again, mutation also increases the diversity of the entire population.
The pseudocode of the GA is presented in Algorithm 1, and relevant figure depicting the algorithm flow is illustrated in Figure 6. To avoid extended discussions on GA, we include here a complete workable source code in the appendix. All parameter settings for the GA are available in the source code.
| ||||||||||||||||||
In the context of Kalman filter, GA is applied to tune the and parameters, which was explained in detail in Section 3.
5. Results
The program, implemented in MATLAB, is given in the appendix to clarify the results obtained in this work. Take note that the and mentioned in (41) are both set to a numerical value of one ( = = 1) in the first simulation. The results obtained are tabulated in Table 2. From these results, RMS of the estimated error, which is the error from KF, is far smaller in comparison to the measured error, with values of V and V, respectively. This RMS error is further minimized by utilizing and , found using the GA approach. A graph depicting the convergence characteristic is shown in Figure 7.
The time plot of the estimated error from 0 s to 60000 s is shown in Figure 8, depicting a very small amplitude, in blue color ( V) along the timeline. This is observed through the zoomed display of the MATLAB graph. On the contrary, the measurement error, portrayed by green color noise in Figure 8, has an absolute magnitude of V.
5.1. Charging Behaviour
The charging characteristic is illustrated in Figure 9 whereby the initial terminal voltage starts from 0 V and rises up to approximately 0.5 V within 60000 seconds (which is 100 minutes). Charging impulses of 1.53 A are applied in this case as shown in Figure 9. As expected, this is a time consuming process as in general case it may take hours for a car battery (lead acid) to be completely charged.
(a)
(b)
5.2. Discharging
For the discharging process, the initial value of terminal voltage, , is set to V in the MATLAB program. Again, in this case, impulses of 1.53 A are applied, but now in negative form. The dynamic behavior of the discharging characteristic is shown in Figure 10. From this figure, it is observed that the discharging process is similar to the charging process, but now with a linearly decreasing terminal voltage slope. The open terminal voltage drops from 2.2 V to 1.7 V in 60000 seconds (100 minutes); this is similar to the charging process as it takes 100 minutes to reach V from zero potential.
(a)
(b)
6. Conclusion
In this work, we successfully obtained the state variables of the RC model representing a battery in terms of mathematical derivations. The derivations lead to the conclusion that there exist three state variables relevant to a battery’s state-space model. With this state-estimation model, a prominent technique known as the Kalman filter is applied in the aim of estimating state-of-charge for a Battery Management System. From numerical results, the KF is shown to be accurate in predicting the dynamic behavior of the system. This is shown by a very small RMS error of the estimation in comparison to its measurement. The estimated error is further reduced after incorporating the optimized values of and through the offline GA approach. As such, the efficacy of such an approach is, thus, validated.
Appendix
A. MATLAB Source Code, in a Script File
format long;
randn(‘seed’,0) %Same sequence
%Value for resistors and capacitors
Csurface=82.11;
Cbk=88372.83;
Re=0.00375;
Rs=0.00375;
Rt=0.002745;
a=1/(Cbk(Re+Rs));
b=1/(Csurface(Re+Rs));
d=(ReRs)/(Re+Rs);
%State variable matrices
A=−a a 0; b −b 0; (−a+b) 0 (a−b) ;
B=[aRe; bRe;
a (0.5Rs−Rt−d)+ b(0.5Re+Rt+d);
C=0 0 1 ;
D=0;
%Transfer function
figure1; %Figure 1
num, den=ss2tf(A,B,C,D,1);
G=tf(num,den)
step(G),grid;
% For Kalman filter:
% Identity matrix + diagonal element
A=[1−a a 0; b 1−b 0; (−a+b) 0 1+(a−b) ]
B=[aRe; bRe;
a (0.5Rs−Rt−d)+ b(0.5Re+Rt+d)
C= 0 0 1
Tc=1;
A=ATc;
B=BTc;
C=C;
x=2,2
options = gaoptimset(…
‘PopulationType’, ‘doubleVector’,…
‘PopInitRange’, 0;1,…
‘PopulationSize’, 5,…
‘EliteCount’, 1,…
‘CrossoverFraction’, 0.8,…
‘ParetoFraction’, 0.35, …
‘MigrationDirection’, ‘both’, …
‘MigrationInterval’, 20,…
‘MigrationFraction’, 0.2,…
‘Generations’, 50,…
‘TimeLimit’, Inf,…
‘FitnessLimit’, −Inf,…
‘StallGenLimit’, 50,…
‘StallTimeLimit’, Inf,…
‘TolFun’, 0,…
‘TolCon’, 0,…%1e−6,…
‘InitialPopulation’, x,…
‘InitialScores’, ,…
‘InitialPenalty’, 10,…
‘PenaltyFactor’, 100,…
‘CreationFcn’, @gacreationuniform,…
‘FitnessScalingFcn’, @fitscalingrank,…
‘SelectionFcn’, @selectionroulette,…
‘CrossoverFcn’, @crossoverarithmetic,…
‘MutationFcn’, @mutationadaptfeasible,…
‘DistanceMeasureFcn’, @distancecrowding,…
‘HybridFcn’, ,…
‘Display’, ‘final’,…
‘OutputFcns’, ,…
‘PlotFcns’, @gaplotbestf, …
‘PlotInterval’, 2, …
‘Vectorized’, ‘off’, …
‘UseParallel’, ‘always’);
Fmin =@(x)fobj(x, A, B, C);
lb=0.001,0.001; %Set lower bounds
vX,fval,exitflag =…
ga(Fmin, 2, ,,,…
, lb, ,, options );
% After simulation,
% repeat simulation all plot all graphs
if fval<0.0001
% Sample time=−1 for discrete model
Plant = ss(A,B B,C,0,−1,…
‘inputname’,{‘u’ ‘w’},…
‘outputname’, ‘y’);
kalmf,L,P,M = kalman(Plant,Q,R);
kalmf = kalmf(1,:);
Kalmf
a = A;
b = B B 0B;
c = C;C;
d = 0 0 0;0 0 1;
P = ss(a,b,c,d,−1,…
‘inputname’,{‘u’ ‘w’ ‘v’},…
‘outputname’,{‘y’ ‘yv’});
sys = parallel(P,kalmf,1,1,,)
% Close loop for input #4 and output #2
SimModel = feedback(sys,1,4,2,1)
% Delete yv from I/O list
SimModel = SimModel(1 3,1 2 3)
SimModel.inputname
%This part generates rectangular waveform
A=1.53; %Amplitude of negative pulses
t = 0:1:60000;
T=30000/5.5;
% I here represents current
I = (A/2)square(2pi(1/T)t);
I = I+(A/2);
figure1; %Figure 2
subplot(211),plot(t,I),
axis(0 60000 −0.5 2.0);
grid on; %charging
xlabel(‘Time (sec)’);
ylabel(‘Cell terminal current (A)’);
u = I’;
n = length(t);
%randn(‘seed’,0) %Same sequence
w = sqrt(Q)randn(n,1);
v = sqrt(R)randn(n,1);
[out,x] = lsim(SimModel,[w,v,u]);
y0=0; %This is initial terminal voltage
y = out(:,1)+y0; % true response
ye = out(:,2)+y0; % filtered response
yv = y + v; % measured response
figure2;
plot(t,y, ‘g−−’,t,ye, ‘b−’), grid on;
figure1;
subplot(212),plot(ye, ‘b−’),
axis(0 60000 0 0.7); grid on; %charging
xlabel(‘Time (s)’),
ylabel(‘Cell terminal voltage (V)’)
%Kalman filter response
figure3; %Figure 3
plot(t,y−yv, ‘g−’, t, y−ye, ‘b−’), grid on;
xlabel(‘Time (s)’), ylabel(‘Error’)
%Calculate Errors
MeasErr = y−yv; %Measurement error
MeasErrCov=...
sum(MeasErr.MeasErr)/length(MeasErr);
EstErr = y−ye; %Estimated error
EstErrCov =…
sum(EstErr.EstErr)/length(EstErr);
%Display onto screen
MeasErrCov %Measurement error
EstErrCov %Estimated error
End
B. MATLAB Source Code, in a Function File
function f = fobj(x, A, B, C)
Q=x(1);
R=x(2);
% Sample time=−1 for discrete model
Plant = ss(A,B B,C,0,−1,…
‘inputname’,{‘u’ ’w’},…
‘outputname’, ‘y’);
kalmf,L,P,M = kalman(Plant,Q,R);
kalmf = kalmf(1,:);
Kalmf
a = A;
b = B B 0B;
c = C;C;
d = 0 0 0;0 0 1;
P = ss(a,b,c,d,−1,…
‘inputname’,{‘u’ ‘w’ ‘v’},…
‘outputname’,{‘y’ ‘yv’});
sys = parallel(P,kalmf,1,1,,)
% Close loop around input #4 and output #2
SimModel = feedback(sys,1,4,2,1)
% Delete yv from I/O list
SimModel = SimModel(1 3,1 2 3)
SimModel.inputname
%This part generates rectangular waveform
A=1.53; %Amplitude of negative pulses
t = 0:1:60000;
T=30000/5.5;
% I here represents current
I = (A/2)square(2pi(1/T)t);
I = I+(A/2);
u = I’;
n = length(t);
%randn(‘seed’,0) %Same sequence
w = sqrt(Q)randn(n,1);
v = sqrt(R)randn(n,1);
out,x = lsim(SimModel,w,v,u);
y0=0; %This is initial terminal voltage
y = out(:,1)+y0; % true response,
ye = out(:,2)+y0; % filtered response
yv = y + v; % measured response
%Calculate Errors
MeasErr = y−yv; %Measurement error
MeasErrCov= …
sum(MeasErr.MeasErr)/length(MeasErr);
EstErr = y−ye; %Estimated error
EstErrCov = …
sum(EstErr.EstErr)/length(EstErr);
f=EstErrCov;
end
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the support from Xi’an Jiaotong-Liverpool University under RDF-13-01-13. The authors would also like to thank their colleague, Sanghyuk Lee, for the sharing of his expertise, contributing to the success of this project.