基于扩展卡尔曼滤波的SOC估计代码

上篇文章讲解卡尔曼滤波如何估计SOC的原理、计算过程，这篇文章将上述过程变为代码展示出来。看完这篇文章您将能快速掌握基于matlab软件实现基于EKF估计SOC的具体代码。

导入数据

将电压、电流、初始SOC值导入

T=1;    %采样时间为1s
N=20042/T;  %总步长
x=zeros(3,N);   %定义状态向量x,定义一个3行N列的零矩阵。
x(:,1)=[0;0;1]; %状态向量x初值设定，初始soc=1
I = xlsread('i.xlsx');     %读取电流数据
u = xlsread('v.xlsx');     %读取电压数据
soc = xlsread('soc.xlsx');  %读取soc数据
I = -I';       %将电流数据变成行向量,放电为正
u = u';        %将电压数据变成行向量
soc = soc';    %将soc数据变成行向量
time = 1:1:N;
Cap = cumtrapz(time, I);
Cn = 33.1;

展示电压、电流图

OCV、二阶等效电路模型参数与SOC的关系

这些参数将在后续使用，由于这部分是可以提前单独完成的，因此单独说一下，如何获取这些参数与SOC的关系，在我之前的文章中已经详细介绍过了。链接：https://blog.csdn.net/a2267897749/article/details/134719450

SOC-OCV曲线

如何获取SOC-OCV曲线在我之前的文章已经介绍过，这里就不重复了，直接列出拟合后的多项式。

 Usoc = 11.08*soc^6-25.58*soc^5+17.54*soc^4-1.159*soc^3-2.386*soc^2 + 1.263*soc+3.422;   %开路电压Usoc是关于SOC的函数

RC参数与SOC的关系

% 阻容参数与SOC的关系
R1 = 0.03023*x(3,j-1)^6-0.1141*x(3,j-1)^5+0.1689*x(3,j-1)^4-0.1243*x(3,j-1)^3+0.04728*x(3,j-1)^2-0.008527*x(3,j-1)+0.0006967;
R2 = 0.003971*x(3,j-1)^6-0.005341*x(3,j-1)^5-0.006872*x(3,j-1)^4+0.01435*x(3,j-1)^3-0.007282*x(3,j-1)^2+0.001231*x(3,j-1)+0.0002082;
C1 = 510100*x(3,j-1)^6-1276000*x(3,j-1)^5+1031000*x(3,j-1)^3-176300*x(3,j-1)^2+69310*x(3,j-1)+909;
C2 = -12940000*x(3,j-1)^6+39640000*x(3,j-1)^5-45050000*x(3,j-1)^4+23320000*x(3,j-1)^3-5467000*x(3,j-1)^2+523200*x(3,j-1)+52340;
if I(1,:) >= 0   %R0区分充放电来拟合，也是关于SOC的函数
    R0=(-0.1313*x(3,j-1)^3+0.4812*x(3,j-1)^2-0.5452*x(3,j-1)+2.96)/1000;%充电时欧姆内阻三阶拟合,除1000表示毫Ω化Ω
else
    R0=(9.033*x(3,j-1)^6-39.36*x(3,j-1)^5+65.72*x(3,j-1)^4-54.46*x(3,j-1)^3+24.21*x(3,j-1)^2 -5.774*x(3,j-1)+2.58)/1000;%放电时欧姆内阻6阶拟合 除1000表示毫Ω化Ω
end

这其中的x就是上篇文章中的状态变量，x的第一行是U1、第二行是U2、第三行是SOC。 j 时刻的阻容参数由 j - 1 时刻的SOC决定，因此 x(3,j-1) 代表了上一个时刻的SOC，这个在之后的代码中会有更加深刻的体现。

基于EKF的SOC估计

卡尔曼滤波的超参数

卡尔曼估计中要先定义P、Q、R矩阵，这三个超参数可以很大程度上决定SOC估计结果的好坏，很多论文中的对比结果往往都暗改这个参数来修改结果。
P：决定卡尔曼滤波的收敛速度，越大则收敛速度越快
Q：过程噪声，过程噪声越小越信任安时积分过程
R：测量噪声，测量噪声越小越信任利用端电压获取的SOC结果

卡尔曼滤波过程

P=eye(3);
Q = 0.00005;
R = 1;

for j = 2:N
    % 导入R、C与SOC的关系
    R1 = 0.03023*x(3,j-1)^6-0.1141*x(3,j-1)^5+0.1689*x(3,j-1)^4-0.1243*x(3,j-1)^3+0.04728*x(3,j-1)^2-0.008527*x(3,j-1)+0.0006967;
    R2 = 0.003971*x(3,j-1)^6-0.005341*x(3,j-1)^5-0.006872*x(3,j-1)^4+0.01435*x(3,j-1)^3-0.007282*x(3,j-1)^2+0.001231*x(3,j-1)+0.0002082;
    C1 = 510100*x(3,j-1)^6-1276000*x(3,j-1)^5+1031000*x(3,j-1)^3-176300*x(3,j-1)^2+69310*x(3,j-1)+909;
    C2 = -12940000*x(3,j-1)^6+39640000*x(3,j-1)^5-45050000*x(3,j-1)^4+23320000*x(3,j-1)^3-5467000*x(3,j-1)^2+523200*x(3,j-1)+52340;
    if I(1,:) >= 0   %R0区分充放电来拟合，也是关于SOC的函数
        R0=(-0.1313*x(3,j-1)^3+0.4812*x(3,j-1)^2-0.5452*x(3,j-1)+2.96)/1000;%充电时欧姆内阻三阶拟合,除1000表示毫Ω化Ω
    else
        R0=(9.033*x(3,j-1)^6-39.36*x(3,j-1)^5+65.72*x(3,j-1)^4-54.46*x(3,j-1)^3+24.21*x(3,j-1)^2 -5.774*x(3,j-1)+2.58)/1000;%放电时欧姆内阻6阶拟合 除1000表示毫Ω化Ω
    end
    H = 66.48*x(3,j-1)^5 - 127.9*x(3,j-1)^4 + 70.16*x(3,j-1)^3 - 3.48*x(3,j-1)^2....
           - 4.58*x(3,j-1) + 1.26;
    A=[exp(-1/(R1*C1)) 0 0;0 exp(-1/(R2*C2)) 0;0 0 1];%系数矩阵A  3行3列
    B=[R1*(1-exp(-1/(R1*C1))) R2*(1-exp(-1/(R2*C2))) 1/(3600*Cn)]';%系数矩阵B 3行1列
    C=[-1 -1 H];%系数矩阵C 1行3列

    %KF过程，理论过程见上一篇文章
    xpre(:,j-1) = A*x(:,j-1)+B*I(1,j-1);  %状态x的预测更新
    Ppre = A*P*A'+eye(3)*Q;   %协方差预测更新  
    K(:,j-1) = Ppre*C'*inv(C*Ppre*C'+R);  %卡尔曼增益更新
    Usoc = 11.08*xpre(3,j-1)^6-25.58*xpre(3,j-1)^5+17.54*xpre(3,j-1)^4-1.159*xpre(3,j-1)^3-2.386*xpre(3,j-1)^2....
        +1.263*xpre(3,j-1)+3.422;   %开路电压Usoc是关于SOC的函数
    
    Um(1,j-1) = Usoc - xpre(1,j-1) - xpre(2,j-1) - I(1,j-1)*R0;%预测的电池端电压
    e(1,j-1) = u(1,j-1)-Um(1,j-1);%新息
    x(:,j) = xpre(:,j-1) + K(:,j-1)*e(1,j-1);
    P=(eye(3)-K(:,j-1)*C)*Ppre;%协方差更新
end

对A、B矩阵进行展开展开

上述过程已经完成了基于KF的SOC估计。为了减少指数运算带来的麻烦，对A、B矩阵进行线性化，提高运算速度。

A=[1-(1/(R1*C1)) 0 0;0 1-(1/(R2*C2)) 0;0 0 1];  %系数矩阵A  3行3列
B=[1/C1 1/C2 1/(3600*Cn)]'; %系数矩阵B 3行1列

SOC估计结果展示

这样基于EKF的SOC估计的整个过程展示完毕。在实际BMS应用时往往会综合各种情况考虑，不会仅使用单一的算法。
本文主要介绍了基于卡尔曼滤波的SOC估计代码。
如果有什么问题交流可以留言或QQ：2267897749一起讨论，如果有您希望直接购买到可以直接使用的代码和数据，可以在闲鱼中搜索用户 BMS小旭 进行购买，价格十分便宜。