无迹卡尔曼滤波估计SOC的simulink模型详解

      上一篇文章讲了如何使用扩展卡尔曼滤波完成SOC的估计，今天就来讲一个更高阶的方法-------无迹卡尔曼滤波（UKF），也有叫sigma点卡尔曼滤波的，实际上都是一样的，下面我都用UKF来表示。simulink模型如图1所示。

 图1 UKF估计SOC的simulink模型

这个模型主要分为两大模块：都是用f函数来完成，

1..首先介绍第一个function函数：function [vt,socdot,vpedot,vpcdot,y,voc]  = fcn(s,vpe,vpc,I)

function [vt,socdot,vpedot,vpcdot,y,voc]  = fcn(s,vpe,vpc,I)
%% Battery Parameters
Rpe=0.011097142857143;
Rpc=0.002654285714286;
Cpe=860.680852654388;
Cpc=381904.908951141;
Rin=0.03807;
Cn=14*3600;
k=1.193;
%%
a1 = 1/(Rpe*Cpe)+1/(Rpc*Cpc);
a2 = 1/(Rin*Cn);
a3 = 1/(Rpc*Cpc);
a4 = 1/(Rpe*Cpe);
b1 = k/Cn + 1/Cpe + 1/Cpc + Rin/(Rpe*Cpe)+ Rin/(Rpc*Cpc);
b2 = 1/Cpe;
b3 = 1/Cpc;
% uncer=0.2;
%%
     %voc = -201.8*s^6 + 739.3*s^5 -1063*s^4 + 763.2*s^3 + -285.9*s^2 + 52.94*s + 9.514;
%voc=k*s+12.67;
% voc=3.608-1.209*s^4 + 3.055*s^3 - 2.215*s^2 +0.9301*s;
voc=12.21*s^3 -21.4*s^2 + 11.74*s + 11.45;
vt=voc-vpe-vpc-Rin*I;
socdot=a2*vt - a2*(voc) + a2*vpc + a2*vpe;%-0.0326*uncer;
vpedot=-a4*vpe + b2*I;%+0.0122*uncer;
vpcdot=-a3*vpc+b3*I;%+0.0204*uncer;
y=vt;

输入：

S：SOC

vpe：R1两端极化电压

vpc：R2两端极化电压

I：电流

输出：

Vt：R0两端的电压

Socdot：下一时刻的SOC

Vpedot：下一时刻R1两端极化电压

Vpcdot：下一时刻R1两端极化电压

y：y=vt,代码中令y=vt，即R0两端的电压

voc：求得的开路电压

需要修改参数的地方如下：

1.Rpe应对应我们实验的R1，

2.Rpc应对应我们实验的R2，

3.Cpe应对应我们实验的C1，

4.Cpc应对应我们实验的C2，

5.Rin应对应我们实验的R0，

6.Cn应对应我们实验的电池容量，

7.k是SOC-OCV的系数，代码中令voc=k*s+3.625;k是拟合得到的。

2.其次是右边的function函数：function [P_o,State_e_o] = fcn(V_t,I,P_i,state_e_i)

function [P_o,State_e_o] = fcn(V_t,I,P_i,state_e_i)
%%
%{
sigma_xx(:,k)=[
    sigma_x(1,k)=soc(k-1);
    sigma_x(2,k)=Vpe(k-1);
    sigma_x(3,k)=Vpc(k-1);
]

sigma_xxx(:,i)-->
sigma_xxx(1,k)= soc(k-1);
sigma_xxx(2,k)= Vpe(k-1);
sigma_xxx(3,k)= Vpc(k-1);
sigma_z(1,k)=vt;
%}
%% Battery Parameters
R_pe=0.011097142857143;
R_pc=0.002654285714286;
C_pe=860.680852654388;
C_pc=381904.908951141;
R_in=0.03807;
C_n=14*3600;
%% 
n=3;
n_m=2*n+1;
kappa=0;
alpha=0.6;
beta=300;
%% 7.34
lambda=alpha^2*(n+kappa)-n;
wm = ones(n_m,1)*1/(2*(n+lambda));
wc = wm;
wm(1)= lambda/(lambda+n); 
wc(1)= lambda/(lambda+n)+1-alpha^2+beta;
sigma_x=zeros(3,n_m);
sigma_xx=zeros(3,n_m);
sigma_xxx=zeros(3,n_m);
sigma_z=zeros(1,n_m);
P_sqrt=chol(P_i,'lower');
%%
% Q = diag([0.02 0.1 0.2]);R=0.2;   
%% 7.39 Calculate sigma point
for i=1:1:n_m
    if i==1
        sigma_x(:,i)= state_e_i;
    elseif i>=2 && i<=n+1
        sigma_x(:,i)= state_e_i+sqrt(n+lambda)*(P_sqrt(:,i-1));
    elseif i>n+1
        sigma_x(:,i)= state_e_i-sqrt(n+lambda)*(P_sqrt(:,i-1-n));
    end
end
%% Time Update 7.41
for k=1:1:n_m 
    sigma_xx(:,k)=[sigma_x(1,k)-0.01*(I/C_n);
                   sigma_x(2,k)+0.01*((-1/(R_pe*C_pe))*sigma_x(2,k)+I/C_pe);
                   sigma_x(3,k)+0.01*((-1/(R_pc*C_pc))*sigma_x(3,k)+I/C_pe)
                   ];
end
%% 7.41
x_n=sigma_xx*wm;
P_n=zeros(3,3);
%% 7.42
for k=1:n_m
    P_n=P_n+wc(k)*(sigma_xx(:,k)-x_n)*(sigma_xx(:,k)-x_n)'  ;
end

P_sqrt=chol(P_n,'lower');
%% 7.39 Calculate sigma point
for i=1:1:n_m
    if i==1
        sigma_xxx(:,i)= x_n;
    elseif i>=2 && i<=n+1
        sigma_xxx(:,i)= x_n+sqrt(n+lambda)*(P_sqrt(:,i-1));
    elseif i>n+1
        sigma_xxx(:,i)= x_n-sqrt(n+lambda)*(P_sqrt(:,i-1-n));
    end
end

%% Time Update 
for k=1:1:n_m
% V_oc=3.608-1.209*sigma_xxx(1,k)^4 + 3.055*sigma_xxx(1,k)^3 - 2.215*sigma_xxx(1,k)^2 +0.9301*sigma_xxx(1,k);
     V_oc = 12.21*sigma_xxx(1,k)^3 -21.4*sigma_xxx(1,k)^2 + 11.74*sigma_xxx(1,k) + 11.45;
%V_oc=1.193*sigma_xxx(1,k)+12.67;
sigma_z(:,k)=V_oc-sigma_xxx(2,k)-sigma_xxx(3,k)-I*R_in;
end
z_n=sigma_z*wm;
%% Measurement update equations
Pyy=10^-5;
Pxy=zeros(n,1);
%% 7.45 , 7.46
for k=1:1:n_m
    Pyy=Pyy+wc(k)*(sigma_z(:,k)-z_n)*(sigma_z(:,k)-z_n)'  ;
    Pxy=Pxy+wc(k)*(sigma_xxx(:,k)-x_n)*(sigma_z(:,k)-z_n)';
end
%% 7.47-49
K=Pxy/Pyy;
State_e_o=x_n+K*(V_t-z_n);
P_o=P_n-K*Pyy*K';

输入：

V_t：R0电压

I：电流

P_i：协方差矩阵，为对角矩阵，初值为主元素为5的3阶对角函数

state_e_i：状态矩阵，1行3列，分别为SOC、R1两端的端电压U1、R2两端的端电压U2.

输出：P_o：下一时刻的协方差矩阵

State_e_o：下一时刻的状态

需要修改的参数有：

Rpe对应我们实验的R1，

Rpc对应我们实验的R2，

Cpe对应我们实验的C1，

Cpc对应我们实验的C2，

Rin对应我们实验的R0，

Cn对应我们实验的电池容量，

k是SOC-OCV的一个系数，代码中令voc=k*s+3.625;k是拟合得到的。

3.仿真结果验证

         SOC对比如下所示：

        给定SOC初始值为0.2，可见UKF算法在初始值不精确的情况下也能够迅速收敛。相比于扩展卡尔曼滤波具有更高的精度。图片有点不太清晰懒得用visio再画了（手动狗头）。扩展卡尔曼和无迹卡尔曼都讲了，下期讲点什么呢？也许是神经网络吧，最简单的BP神经网络。欢迎大家在评论区提出问题，谢谢，今天就写到这里该吃饭啦啦啦啦啦啦。

硕博期间所有的程序代码，一共2个多g，可以给你指导，赠送半个小时的语音电话答疑。电池数据+辨识程序+各种卡尔曼滤波算法都在里面了，后续还会有新模型的更新。快速入门BMS软件。某鹅：2629471989