扩展卡尔曼滤波(EKF)估计电池SOC 含MATALAB代码

x小歌

已于 2024-05-14 14:53:38 修改

阅读量3k

点赞数 25

文章标签： matlab

于 2024-05-14 14:48:23 首次发布

本文链接：https://blog.youkuaiyun.com/xsnhj/article/details/138851384

版权

一、卡尔曼滤波

卡尔曼滤波理论（Kalman Filter, KF）是在1960年由鲁道夫·卡尔曼提出，该滤波是基于最小均方差准则，是一种利用上一时刻的估计值和当前时刻的观测值来得出当前时刻的最优状态值估计值的线性离线系统滤波器，其状态方程和观测方程如式(1-1)和式(1-2)所示。

$x_{k+1}=Ax_{k}+Bu_{k}+w_{k}$ (1-1)

$y_{k}=Cx_{k}+Du_{k}+v_{k}$ (1-2)

式中，xk为系统的状态变量；yk为系统的观测向量；uk为系统的输入量；k与k+1指的是时刻；A为B分别为系统的传递矩阵和输入矩阵；C和D为系统的输出矩阵和前馈矩阵；wk和vk分别为系统的过程噪声和测量噪声，均为高斯白噪声且互不相关。

卡尔曼滤波算法的计算流程如下所示

（1）初始化

定义初始值：

$\hat{x}=E(x)$ (1-3)

$P_{0}=E\left \{ \left ( x_{0}-\hat{x_{0}} \right )\left ( x_{0}-\hat{x_{0}} \right ) \right \}$ (1-4)

（2）时间更新

先验状态估计 :

$\hat{x}_{k+1}=A\hat{x}_{k}+Bu_{k}$ (1-5)

计算先验状态协方差矩阵：

$P_{k+1}=AP_{k}A^{T}+Q_{k}$ (1-6)

（3）测量更新

计算卡尔曼增益 :

$K_{k+1}=P_{k+1}C_{k+1}^{T}\left (C_{k+1} P_{k+1}C_{k+1}^{T}+R_{k}\right )^{-1}$ (1-7)

计算k+1时刻系统的最优估计值

$\hat{x}_{k+1}=\hat{x}_{k+1}+K_{k+1}(y_{k+1}-C_{k+1}\hat{x}_{k+1}-Du_{k})$ (1-8)

计算后验状态协方差矩阵

$P_{k+1}=(I-K_{k+1}C_{k+1})P_{k+1}$ (1-9)

中，I为单位矩阵。

二、扩展卡尔曼滤波

卡尔曼滤波是基于线性离散系统的滤波器，为了可以将卡尔曼滤波同样适用于非线性系统，R.S.Bucy与K.D.Senne于20世纪70年代初提出了扩展卡尔曼滤波（Extended Kalman Filter，EKF），该方法的思想是将非线性系统线性化逼近，转为线性化问题。EKF是将非线性系统通过泰勒展开，保留一阶项，去除高阶项来达到近似线性化的目的，然后利用标准的卡尔曼滤波体系对非线性系统的状态进行估计。

非线性系统的表达式如式(2-1)、(2-2)所示

$x_{k+1}=f(x_{k},u_{k})+w_{k}$ (2-1)

$y_{k+1}=g(x_{k},u_{k})+v_{k}$ (2-2)

其中，f为非线性系统的传递函数，g为非线性测量函数，wk和vk为高斯白噪声。

将非线性函数进行一阶泰勒展开：

$f(x_{k},u_{k})\approx f(\hat{x}_{k},u_{k})+\frac{\partial f(x_{k},u_{k})}{\partial x_{k}}|_{x_{k}=\hat{x}_{k}}(x_{k}-\hat{x}_{k})$ (2-3)

$g(x_{k},u_{k})\approx g(\hat{x}_{k},u_{k})+\frac{\partial g(x_{k},u_{k})}{\partial x_{k}}|_{x_{k}=\hat{x}_{k}}(x_{k}-\hat{x}_{k})$ (2-4)

令

$A_{k}=\frac{\partial f(x_{k},u_{k})}{\partial x_{k}}|_{x_{k}=\hat{x}_{k}}$ (2-5)

$C_{k}=\frac{\partial g(x_{k},u_{k})}{\partial x_{k}}|_{x_{k}=\hat{x}_{k}}$ (2-6)

式（2-1），（2-2）可改写为：

$x_{k+1}\approx A_{k}x_{k}+[ f(\hat{x}_{k},u_{k})-A_{k}\hat{x}_{k}]+w_{k}$ (2-7)

$y_{k+1}\approx C_{k}x_{k}+[ g(\hat{x}_{k},u_{k})-C_{k}\hat{x}_{k}]+v_{k}$ (2-8)

扩展卡尔曼滤波(EKF)算法的计算流程如下所示:

（1）定义初始值：

$\hat{x}=E(x)$ (2-9)

$P_{0}=E\left \{ \left ( x_{0}-\hat{x_{0}} \right )\left ( x_{0}-\hat{x_{0}} \right ) \right \}$ (2-10)

（2）状态变量估计

$\hat{x}_{k+1}=f(\hat{x}_{k},u_{k})$ (2-11)

（3）计算协方差矩阵的先验值

$P_{k+1}=AP_{k}A^{T}+Q_{k}$ (2-12)

（4）卡尔曼增益

$K_{k+1}=P_{k+1}C_{k+1}^{T}\left (C_{k+1} P_{k+1}C_{k+1}^{T}+R_{k}\right )^{-1}$ (2-13)

（5）修正状态变量

$x_{k+1}=\hat{x}_{k+1}+K_{k+1}(y_{k+1}-g(\hat{x}_{k},u_{k}))$ (2-14)

（6）更新协方差矩阵

$P_{k+1}=(I-K_{k+1}C_{k+1})P_{k+1}$ (2-15)

中，I为单位矩阵。

三、MATLAB代码

clear;
clc;
close all;

load("HPPC.mat")
Cn=98.15*3600;%电池容量，单位Ah
delta_t=1;%采样时间
ro=0.0026;%电池内阻
rp1=9.227E-4;
cp1=4.943E4; 
rp2=4.715E-4;
cp2=7.290E6;%RC环节的内阻和电容值，单位欧姆和法

nita=1;%放电系数
i=-DATA.I(4466:55718);
soc=DATA.SOCT;
ocv=DATA.U(4466:55718);
SOC=DATA.SOC/100;
OCV=DATA.OCV;
p=polyfit(SOC,OCV,7);  %%开路电压七阶曲线拟合系数
fn=poly2sym(p);  %%函数
g=diff(fn);   %%函数求导
p1=sym2poly(g);   %%求导系数

A=[1 0 0;0 exp(-delta_t/(rp1*cp1)) 0;0 0 exp(-delta_t/(rp2*cp2))];
B=[-nita*delta_t/Cn;(1-exp(-delta_t/(rp1*cp1)))*rp1;(1-exp(-delta_t/(rp2*cp2)))*rp2];%A，B矩阵
t=DATA.T(4466:55718)-4466;
N=length(t);

size=[3,N];%数组大小
Q=4e-10*eye(3);%状态方程中的误差
R=1e-2; %观测方程中的误差
X(:,1)=[1;0;0];%状态变量初值

for k=2:N
    X(:,k)=A*X(:,k-1)+B*i(k)+Q*randn(3,1);
end
for k=1:N
    Soc(k)=X(1,k);
    Up1(k)=X(2,k);
    Up2(k)=X(3,k);
    E(k)=polyval(p,Soc(k));
    Uoc(k)=E(k)-i(k)*ro-Up1(k)-Up2(k)+R*randn(1);
end   %求取仿真所用的真值

P0=eye(3);%协方差矩阵初值
soc_input=1;   %%%soc初值
Xekf(:,1)=[soc_input;0;0];%EKF状态变量初值，可更改

for k=2:N
    Xn=A*Xekf(:,k-1)+B*i(k-1);%计算先验估计
    P1=A*P0*A'+Q;%更新协方差矩阵
    Zm=polyval(p,Xn(1))-Xn(2)-Xn(3)-i(k)*ro;%计算观测值
    H=[polyval(p1,Xn(1)) -1 -1];%雅克比矩阵计算
    K=(P1*H')/(H*P1*H'+R);%计算卡尔曼增益
    Xekf(:,k)=Xn+K*(Uoc(k)-Zm);%更新状态
    P0=(eye(3)-K*H)*P1;%更新协方差方程
end

for k=1:N
    Socekf(k)=Xekf(1,k);
    Eekf(k)=polyval(p,Socekf(k));
    Uocekf(k)=Eekf(k)-i(k)*ro-Xekf(2,k)-Xekf(3,k);
end%计算三个值

%依次绘制先验数据，后验数据和实验数据做比较
%这里绘制的是soc-ocv曲线
figure;
plot(X(1,:),Uoc,'b');
hold on;
plot(Xekf(1,:),Uocekf,'g');
hold on;
plot(soc,ocv,'r');
grid on;
legend('先验数据','后验数据','实验数据');
xlabel('soc');ylabel('Uocv/V');
title('拓展卡尔曼滤波仿真soc-ocv曲线对比');

%依次绘制先验数据，后验数据和实验数据做比较
%这里绘制的是t-ocv曲线
figure;
plot(t,Uoc,'r');
hold on;
plot(t,ocv,'b');
hold on;
plot(t,Uocekf,'g');
grid on;
xlabel('t/s');ylabel('Uocv/V');
legend('先验数据','实验数据','后验数据');
title('拓展卡尔曼滤波仿真t-Uov曲线对比');

figure;
plot(t,X(1,:),'r');
hold on
plot(t,Xekf(1,:),'b');
hold on
plot(t,soc,'g');
hold on
legend('先验数据','后验数据','实验数据');
title('拓展卡尔曼滤波仿真t-soc曲线对比');

a=Xekf(1,:)-soc';
figure;
plot(t,a,'g');
legend('soc估计误差');
title('soc估计误差');

四、运行结果图