Extended Kalman Filter Problem [Matlab]

Prev: using QccPack
Next: Simulink:SL_TimingWarn

From: Satyanarayana on 17 Mar 2010 06:32

"Robin Roche" <robin.roche(a)removethis.gmail.com> wrote in message <gptint$lmv$1(a)fred.mathworks.com>...
> Hi,
> I am currently developping an algorithm to determine the state of charge of batteries thanks to an extended Kalman filter.
>
> My problem is that although the voltage prediction is good, the state of charge calculation is completely wrong. It varies correctely but the variations are fat to big. For example, in 10ms, according to this algorithm, the state of charge (SoC) can vary of 3%, which is impossible. It should be like a 100 times less!
>
> For this filter I use:
> x: state vector, containing the state of charge (%) and a voltage (V)
> I_meas: current measures
> U_meas: voltage measures
> h_est: voltage prediction
> f: state function, giving its values to x
>
> Does anyone have any idea of what's wrong with code?
>
> Thank you
>
> Corresponding code:
>
> % Loading of the measures
> data = csvread('mlb_test.csv'); % Opens of the measures file
> len=size(data,1); % Number of data samples
> time=data(:,2); % 2nd column of the file: time
> U_meas=data(:,5); % 5th column: voltage
> I_meas=data(:,6); % 6th column: current
>
> % Circuit and measure variables of electric battery model
> eta=1; % Coulombic efficiency: 1 for discharge, <1 for charge
> R_D=0.035; % [Ohm] Internal resistance
> R_Diff=0.0276; % [Ohm] Diffusion resistance
> C_Diff=0.5693; % [F] Diffusion capcity
> tau_Diff=R_Diff*C_Diff; % [s] Time constant of the RC ladder
> Cnom=150; % [Ah] Battery capacity
> %dt=time(2)-time(1);
> dt=0.01;
>
> % Uncertainties initialization
> n=2; % Number of states: x(1):SoC, x(2):V_Diff
> q=0.4; % Stdiv of process (voltage)
> r=0.5; % Stdiv of measurement (current)
> Q=q^2*eye(n); % Covariance of process
> R=r^2; % Covariance of measurement
>
> % Initial states
> x_k_1 = [50;0]; % Initial state: Soc=50%, U=0
> P_k_1 = eye(n); % Initial state covariance
> k = 1;
>
> % OCV(SoC) model
> OCV=@(socm)(69.4*(sqrt(socm)/7+2.75)); % Temporary model example
>
> % Function definition matrices
> Ak=[1 0;0 (1-dt/(R_Diff*C_Diff))];
> Bk=[eta*dt*100/Cnom;dt/C_Diff];
>
> % Main loop
> for k=1:len
> %for k=1:10
>
> % Observation function
> h_est=@(x)(OCV(x(1))-I_meas(k)*R_D-x(2));
>
> % State prediction function
> f=@(x)(Ak*x+Bk*I_meas(k));
>
> % *** START OF EKF ***
>
> % Calculation of the Jacobians
> [x1,A]=jaccsd(f,x_k_1); % Constant value
> [z1,H]=jaccsd(h_est,x1); % Varies with OCV(SoC)
>
> % *** PREDICTION ***
>
> % Project the state ahead
> x_k_ = Ak*x_k_1+Bk*I_meas(k);
>
> % Prediction of the plant covariance
> P_k_ = A*P_k_1*A' + Q;
>
> % *** UPDATE ***
>
> % Kalman gain
> K = P_k_*H'*inv(H*P_k_*H'+R);
>
> % Update estimate with measurement U_meas
> x_k = x_k_ + K*(U_meas(k)-h_est(x_k_));
>
> % Update the error covariance
> P_k = (eye(n)-K*H)*P_k_;
>
> % *** END OF EKF ***
>
> %[x_k, P_k] = ekf(f,x_k_1,P_k_1,h_est,U_meas(k),Q,R);
>
> % STORING OF VALUES
>
> x_store(k,:)=x_k;
> x_store_(k,:)=x_k_;
> h_store(k)=h_est(x_k_);
> %P_store(k,:)=P_k;
>
> % PREPARATION OF THE VARIABLES FOR THE NEXT K
>
> x_k_1=x_k;
> P_k_1=P_k;
>
> % if(k<(len-1))
> % dt=time(k+1)-time(k);
> % end
> % if (dt==0)
> % if (k>1)
> % dt=(time(k-1)+time(k+1))/2;
> % else
> % dt=0.01;
> % end
> % end
>
> end
>
> % Results plot
>
> % State of charge
> subplot(3,1,1)
> plot(time,data(:,14),time,x_store(:,1))
> title('State of Charge')
> h_legend = legend('Actual value','Update',2);
> set(h_legend,'Interpreter','none')
>
> % I current
> subplot(3,1,2)
> plot(time,I_meas,'-')
> title('I: battery current')
>
> % U voltage
> subplot(3,1,3)
> plot(time,U_meas,'-',time,h_store,'--')
> title('U: battery voltage')

Hi robin,
i am new to extended kalman filtering, i have started working on ekf, sorry i cannot give answer to your question. Iam implementing the soc of battery iam facing the following issue where the soc is going beyond 100% could you please let me know the issue with implementation of ekf. i will be greatly thankful to you.

Thanks and regards,
Satya

voltage = xlsread('voltage.xls');
% constants determined from experiments
K0=11.30;
K1=0.02403;
K2=0.91281;
K3=0.24648;
K4=0.335;
R = -2.5 * power(10,-2);

Rv=2;
A=1;
B=0.25; %(T=1sec sampling period; n=1 for discharge C=4 20hr rating)
D=6; %(Internal resistance appx to 6)
P=1;

Current=5; %(assume a constant input current of 5)
SOCinit=0.5; %(assume initial soc of 50%)

%--------------------------------------------------------------------%
%Extrapolate the most recent state estimate to the present state
t = 0:1:999;
for incr=1:1:1000
SOCest=(SOCinit*A - (Current/(4*3600)));
SOC = SOCest;
%--------------------------------------------------------------------%
%--------------------------------------------------------------------%
%Form the innovation vector
%Calculate the Constant C value
C = K2 - K1*(1/pow2(SOC)) + K3*(1/SOC)- K4/(1 - SOC);
%Read the voltage from the battery terminals at that corresponding moment
%VoltMeas = voltage(time);
VoltMeas = voltage(incr);
Voltest = K0 + K1/SOC + K2*SOC + K3*log(SOC) + K4*log(1 - SOC) + R*Current;
Inn= VoltMeas-Voltest;
%-------------------------------------------------------------------%

%-------------------------------------------------------------------%
%Compute the covariance of the estimation
Covar = C * P * C' + Rv;
%-------------------------------------------------------------------%

%------------------------------------------------------------------%
%Compute the kalman gain matrix
K = P*C'*(1/Covar);
%------------------------------------------------------------------%

%------------------------------------------------------------------%
%Update the state estimate
SOCest = SOC + K * Inn;
%------------------------------------------------------------------%

%------------------------------------------------------------------%
n = 1;
%Compute the covariance of the estimation error
P = (eye(n) - K * C)*P;
P = P + 1;
SOCinit = SOCest;
end
%------------------------------------------------------------------%

From: Girikrishna on 21 Mar 2010 23:40

|
Pages: 1
Prev: using QccPack
Next: Simulink:SL_TimingWarn