上一篇文章讲了如何使用扩展卡尔曼滤波完成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再画了（手动狗头）。