Particle Filter for SoC estimation

6 次查看（过去 30 天）
Amjad Abdulla
Amjad Abdulla 2020-5-31
回答： Lorenzo 2024-12-17
Hi I'm currently trying to estimate SoC for lithium-ion battery using Particle Filter.
the code is as follows:
%% clear everything
clear all
close all
clc
%% ******* Battery Parameter ****
dt = .5;
C_nom = 0.850; %% nominal Ah
R0 = 0.07446;
R1 = 0.04669;
C1 = 703.6;
R2 = 0.04984;
C2 = 4475;
SOC0 = 1;
V1=1;
V2=1;
%% load data of nonlinear model (4th order model in the paper)
load('Vb_Data_Measurment.mat');
load('Voc_Figure1.mat');
load('t_time.mat');
load('Current_input.mat');
load('Vsoc_figure.mat');
load('V_B_Figure1.mat');
load('V_soc_Figure1.mat');
%% ******* Generating process noise ,sensor noise & noisy plant response *******
x = 2; % initial actual state
x_N = 10^-6; % Noise covariance in the system (i.e. process noise in the state update, here, we'll use a gaussian.)
x_R = 10^-6; % Noise covariance in the measurement
T = 75; % ( number of iterations).
N = 10; % The number of particles the system generates.
%initilize our initial, prior particle distribution as a gaussian around
%the true initial value
V = [2]; %define the variance of the initial esimate
x_P = []; % define the vector of particles
% make the randomly generated particles from the initial prior gaussian distribution
for i = 1:N
x_P(i) = x + sqrt(V) * randn;
end
%{
%show the distribution the particles around this initial value of x.
%generate the observations from the randomly selected particles, based upon
%the given function
F=[1 -1 -1];
z_out = [ F*x - R0*I + sqrt(x_R) * randn]; %the actual output vector for measurement values.
x_out = [x]; %the actual output vector for measurement values.
x_est = [x]; % time by time output of the particle filters estimate
x_est_out = [x_est]; % the vector of particle filter estimates.
for t = 1:T
SOC(t) = x_P(1);
V1(t) = x_P(2);
V2(t) = x_P(3);
Vf(t)=(-1/(R1*C1)*V1(t)) + (1/C1)*T(t);
Vc(t)= (-1/(R2*C2)*V1(t)) + ((1/C2)*T(t));
OCV(t) = -1.031*exp(-35*SOC(t))+0.2156*SOC(t)-0.1178*SOC(t)^2+0.3201*SOC(t)^3 + 3.685;
Vb_h(t) = OCV(t)-V1(t)-V2(t)-R0*I(t);
A = [1,0,0; 0,exp(-dt/(R1*C1)),0;0,0,exp(-dt/(R2*C2))];
B = [-dt/(3600*C_nom) ; R1*(1-exp(-dt/(R1*C1))) ; R2*(1-exp(-dt/(R2*C2)))];
x = A*x + B*I(t)+ sqrt(x_N)*randn ;
z = F*x - R0*I(t) + sqrt(x_R) * randn;
%Here, we do the particle filter
for i = 1:N
SOC(i) = x_P(1);
V1(i) = x_P(2);
V2(i) = x_P(3);
Vf(i)=(-1/(R1*C1)*V1(i)) + (1/C1)*I(i);
Vc(i)= (-1/(R2*C2)*V1(i)) + ((1/C2)*I(i));
OCV(i) = -1.031*exp(-35*SOC(i))+0.2156*SOC(i)-0.1178*SOC(i)^2+0.3201*SOC(i)^3 + 3.685;
Vb_h(i) = OCV(i)-V1(i)-V2(i)-R0*I(i);
x_P_update(i) = A*x_P(i) + B*I(i) ;
% update the observations
%for each of these particles.
z_update(i) = F*x_P_update(i) - R0*I(i);
P_w(i) = (1/sqrt(2*pi*x_R)) * exp(-(z - z_update(i))^2/(1*x_R));
end
% Normalize to form a probability distribution (i.e. sum to 1).
P_w = P_w./sum(P_w);
%% Resampling: From this new distribution, now we randomly sample from it to generate our new estimate particles
%next iteration
for i = 1 : N
x_P(i) = x_P_update(find(rand <= cumsum(P_w),1));
end
%The final estimate is some metric of these final resampling, such as
%the mean value or variance
x_est = mean(x_P);
% Save data in arrays for later plotting
x_out = [x_out x];
z_out = [z_out z];
x_est_out = [x_est_out x_est];
end
I'm havving error says :
Unable to perform assignment because the indices on the left
side are not compatible with the size of the right side.
Error in Particle_Filter (line 110)
x_P_update(i) = A*x_P(i) + B*I(i)
I can't solve this issue for so long

请先登录，再回答此问题。

回答（1 个）

Lorenzo
Lorenzo 2024-12-17
Hello Amjad,
Is difficult to say which variable is exactly causing the problem since I am not able to simulate your code. However the error is telling you that you are trying to assign a vector to x_P_update that is not matching the dimensions of x_P_update. Your current implementation would work if x_P_update is a 1xn or a nx1 vector (and I guess this is not the case).
In case your x_P_update is something like a 3xn vector, the assignement should be at each n something like:
x_P_update(:,i) = A*x_P(i) + B*I(i)
And you also have to make sure that P(i) and B(i) are correct in dimensionality. Best way to see what happens is to debug the code by setting a debug point on the line causing the error.
Hope this helps,
Best,
Lorenzo

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

尚未输入任何标签。

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by