Particle Filter for SoC estimation

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

Answer 1

[Lorenzo]

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