粒子滤波算法预测电芯的寿命

今天分享一个案例，使用粒子滤波预测电芯的寿命

粒子滤波是一种滤波算法，粒子滤波器已成为求解非线性非高斯情形下最优估计问题的一类非常流行的数值方法。广泛用于机器人、车辆定位。基本原理：随机选取预测域的 N 个点，称为粒子。以此计算出预测值，并算出在测量域的概率，即权重，加权平均就是最优估计。之后按权重比例，重采样，进行下次迭代。此算法存在的问题是：粒子数量越多，计算越精确，但计算量也会增加，且随着迭代，有些粒子权重变得很小，导致粒子枯竭。

2.设置模型初始值

x=unifrnd(1.929, 1.945); %for cell three

b=unifrnd(-0.002038, -0.001947); %for cell three

X0=[x,b]';

M=500; %no of particles

p = 2; %no of parameters

Xparam=zeros(p,N); %Matrix to keep track of changing x and b values

Xparam(:,1)=X0;    %Initialization of parameter matrix

3.测试和处理噪音

4.Monte Carlo 仿真

Xm=zeros(p,M,N);

for i=1:M

   Xm(:,i,1)=X0+sqrtm(Q)*randn(p,1);

end

Xcollection(:,1) = (datasample(Xm(1,:,1),10))';

XcollectionIndex(:,1)=1;

Zm=zeros(1,M,N);

Xestimated=zeros(1,N);

W=zeros(N,M);

for k=2:N

   for i=1:M

       Xm(1,i,k)=Xm(1,i,k-1)*exp(Xm(2,i,k-1)*(k-(k-1)))+sqrt(var_x)*randn();

       Xm(2,i,k)=Xm(2,i,k-1)+sqrt(var_b)*randn();      

   end

 if(mod(k,25)==0 && k<=start )

       ind = size(XcollectionIndex,2);

       Xcollection(:,ind+1) = (datasample(Xm(1,:,k),10))';

       XcollectionIndex(:,ind+1)=k;

   end

   for i=1:M

       Zm(1,i,k)=Xm(1,i,k)+ sd_z*randn();                    

       W(k,i)=exp(-(Zmeasured(1,k)-Zm(1,i,k))^2/2/R)+1e-99; % calculate weight of each particle

   end

   W(k,:)=W(k,:)./sum(W(k,:));

   outIndex = residualR(1:M,W(k,:)');        

   Xm(:,:,k)=Xm(:,outIndex,k);

   if(mod(k,25)==0 && k<=start )

       ind = size(XcollectionIndex,2);

       Xcollection(:,ind+1) = (datasample(Xm(1,:,k),10))';

       XcollectionIndex(:,ind+1)=k+2;

   end

   Xestimated(1,k)=mean(Xm(1,:,k));

   Bestimated(1,k)=mean(Xm(2,:,k));

   Xparam(:,k)=[Xestimated(1,k);Bestimated(1,k)];    

end

5.预测容量

noOfCycles = 250;

initial_value=Zmeasured(1,1);

threshold = 0.75;

threshold_capacity = threshold*initial_value;

flag = 0;

Xextrapolated(1,1)= Xestimated(1,start);

Bmean = Bestimated(1,start);

cycle(1,1)=start;

%%

for k=start+1:noOfCycles

   Xextrapolated(1,k-start+1)=Xextrapolated(1,k-start)*exp(Bmean*dt);

   cycle(1,k-start+1)=k;

   %RUL

   if Xextrapolated(1,k-start+1)<= threshold_capacity && flag==0

   flag = 1;

   failure_index = k;

   failure_capacity = Xextrapolated(1,k-start+1);

   end      

end

noOfCyclesLeft =  failure_index-start;

6.后期图像处理

figure

hold on;box on;

plot(Xestimated,'-r*')  

plot(Zmeasured,'-b.')  

plot(cycle,Xextrapolated,'-.g.')

plot([start start]',[0 2]','g','LineWidth',2);

plot([failure_index failure_index]',[0 failure_capacity]','r','LineWidth',2);

title({'Battery Capacity vs Cycle';'Ck = Ck-1 *exp(b*dt)'});

ylabel('Battery Capacity');

xlabel('Cycle');

legend('Predicted Capacity','Measured Capacity','Extrapolated Capacity','Prediction Start Point','Failure Threshold');

RULText = ['RUL is ',num2str(noOfCyclesLeft),' Cycles'];

text(150,1.65,RULText,'FontSize',14);

text(10,2.4,'cell 3','FontSize',14);

for i=1:length(XcollectionIndex)

     for j=1:10

       if mod(i,2) == 0

           scatter(XcollectionIndex(i),Xcollection(j,i),'filled',...

               'MarkerEdgeColor',[0 .5 .5],...

               'MarkerFaceColor',[0 .7 .7],...

               'LineWidth',1.5);

       else

           scatter(XcollectionIndex(i),Xcollection(j,i),'filled',...

               'MarkerEdgeColor',[0.7 0.7 0],...

               'MarkerFaceColor',[1 1 0],...

               'LineWidth',1.5);

       end

       end

end