电池SOC估算中的EKF与UKF：MATLAB双滤波工程案例拆解

简介针对电池荷电状态估算需求这份MATLAB实现资源聚焦扩展卡尔曼滤波EKF与无迹卡尔曼滤波UKF两种算法的落地对比适合电池管理系统研究人员、电力电子方向学生及从事SOC估算开发的工程师。压缩包内共3个文件包含2个m脚本与1个slx仿真模型分别对应EKF/UKF核心算法实现、主程序入口及BBDST工况模型可支撑从算法原理验证到工况仿真的完整流程。已有247人学习下载。资料提供可直接运行的 MATLAB 代码并内置了电池Thevenin等效模型参数与典型驾驶工况帮助使用者在同一框架下比较EKF与UKF的估算精度、收敛性和鲁棒性通过调整噪声协方差与模型参数还可快速完成算法调优为实际BMS中的SOC估计提供实验参考。1. 电池SOC估算为什么绕不开EKF和UKF拆解一个MATLAB双滤波工程做电池SOC估算的人迟早会撞上同一个问题电压、电流都能直接测唯独SOC测不了只能靠模型和滤波去推。EKF和UKF是应用最广的两条非线性滤波路径一个每步做泰勒展开线性化一个用sigma点采样绕开线性化。这个MATLAB工程把两种算法放在同一个Thevenin模型下实现包含两套滤波递推函数、BBDST动态工况Simulink模型和可直接运行的主脚本mainc.m跑一遍就能比较SOC估算的收敛速度与稳态误差。它面向正在做BMS算法、写毕业设计或者想对比EKF/UKF实际表现的工程师拿到就能复现不需要额外装工具箱。2. 先把模型立住Thevenin等效电路、OCV查表与BBDST工况2.1 为什么资源里选的是Thevenin一阶RC模型电池等效电路模型从简单到复杂排最基础的是Rint模型一个理想电压源串一个欧姆内阻R0只能描述电流突变瞬间的电压跌落完全抓不住电流变化之后的电压弛豫回弹。Thevenin模型在Rint基础上加了一组RC并联网络Rp描述极化内阻Cp描述极化电容这一组RC正好对应电池在电流阶跃之后端电压的指数恢复过程。这个动态特性对SOC估算不是锦上添花而是必要前提。原因在于滤波器的测量方程写成Vt OCV(SOC) - R0·I - Vp极化电压Vp是叠加在端电压上的干扰项。如果模型里没有Vp这个状态滤波器会把端电压里的动态分量全部当成OCV变化去修正SOC电流一变SOC估算就跟着抖。把Vp显式建出来之后滤波器才分得清端电压里哪部分是极化电压、哪部分是真实的OCV变化这是整个估算精度能否成立的基础。一阶RC和二阶RC的取舍也值得说明。二阶模型用两组RC分别描述电化学极化和浓差极化精度上限更高但要辨识Rp1、Cp1、Rp2、Cp2四个参数加上R0和OCV曲线标定工作量翻倍。对SOC估算这个目标一阶Thevenin的极化建模已经能覆盖大多数动态工况。资源里采用一阶模型状态量只有SOC和Vp两个协方差矩阵是2×2EKF的雅可比简单UKF每拍只要传播5个sigma点计算开销和调试难度都在可控范围内。参数表如下参数含义典型示例值Cn额定容量(Ah)20R0欧姆内阻(mΩ)15Rp极化内阻(mΩ)8Cp极化电容(F)800tau极化时间常数(s)6.42.2 状态方程离散化与OCV查表实现一阶Thevenin模型的连续状态方程是d(SOC)/dt -I / (Cn·3600)dVp/dt -Vp / (Rp·Cp) I / Cp其中Cn的单位是Ah所以要乘3600换算成安秒电流I约定放电为正。按零阶保持法离散化即假设一个采样周期内电流恒定得到两个递推式SOC(k1) SOC(k) - I(k)·dt / (Cn·3600)Vp(k1) exp(-dt/tau)·Vp(k) (1 - exp(-dt/tau))·Rp·I(k)第一个式子是安时积分的离散形态SOC对时间的依赖很弱这是EKF里状态转移矩阵A第一行为[1, 0]的原因。第二个式子里exp(-dt/tau)是RC回路的自然响应系数tau Rp·Cp是极化时间常数锂电池的tau一般在几秒到几十秒量级。这里有个容易被忽略的设计约束采样周期dt应该远小于tau工程经验是dt不超过tau的十分之一否则离散化误差会吃掉极化动态。资源里dt取0.1s对tau约6.4s的电池是合适的如果换电池后tau只有0.5sdt就要跟着缩到0.05s甚至更小。OCV与SOC的关系是测量方程里唯一的非线性来源代码里用查表加插值实现function OCV OCV_SOC_Lookup(soc) % soc: 0~1的小数ocv_table: 第一列SOC第二列OCV单位V global ocv_table; OCV interp1(ocv_table(:,1), ocv_table(:,2), soc, linear, extrap); endinterp1默认线性插值。OCV曲线在SOC两端斜率变化很快线性插值会带一点误差但对SOC估算影响不大真正要注意的是查表方向SOC输入超过表范围时用extrap外推避免返回NaN。如果滤波过程中SOC冲到0或1附近的边界区外推能保证计算继续但边界外的OCV值已经失真此时估算只能依赖模型积分误差会明显增大属于算法本身的边界限制。2.3 参数来源与BBDST工况数据入口R0、Rp、Cp这三个参数资源里是直接给好的但换到别的电池上必须重新辨识。常见做法是HPPC脉冲测试电池静置后放一个10s左右的电流脉冲再静置40s以上。电流阶跃瞬间的电压突降对应R0断电瞬间的电压突升也对应R0弛豫阶段的指数恢复曲线拟合出tau脉冲幅值结合电流推Rp再由tau Rp·Cp算Cp。拟合用lsqcurvefit或者cftool都行关键是每段脉冲前的静置要足够久让极化电压完全归零。BBDST工况在资源里是BBDST_workingcondition.slx一个Simulink模型。BBDST是北京市公交车动态测试工况一段电流序列里包含启动加速、匀速、制动回馈和怠速停靠电流从负值充电到2C以上的放电峰值都有。它的价值在于持续激励极化回路让Vp反复充放滤波器的动态响应在低速平台区、回充切换点这些场景下差异会暴露得最明显。模型输出到工作空间的变量一般是时间向量和电流向量mainc.m里load进来之后按索引逐拍读取。手头有实测的BBDST电流数据时替换模型里的数据源就能复用整套算法不需要改滤波代码。3. EKF实现从雅可比矩阵到SOC估算曲线3.1 EKF五步递推与代码骨架EKF的思路是把非线性测量方程在每步估计点做一阶泰勒展开用雅可比矩阵代替标准卡尔曼里的常数矩阵。核心函数是EKF_Step每调用一次完成一个采样周期的递推。下面是从资源里的EKF_UKF_SoCEstimation_Theveninc.m中抽出的骨架文件名里Thevenin被截断成Theveninc不影响运行function [soc_ekf, x_est, P] EKF_Step(x_est, P, I, Vt, Q, R, param) % 状态: x [SOC; Vp]输入: I量测: Vt Cn param.Cn; R0 param.R0; Rp param.Rp; Cp param.Cp; dt param.dt; tau Rp * Cp; % 极化时间常数 % 1. 状态预测 x_pred zeros(2,1); x_pred(1) x_est(1) - I * dt / (Cn * 3600); x_pred(2) exp(-dt/tau) * x_est(2) (1 - exp(-dt/tau)) * Rp * I; % 2. 协方差预测 A [1, 0; 0, exp(-dt/tau)]; P_pred A * P * A Q; % 3. 测量预测 OCV OCV_SOC_Lookup(x_pred(1)); V_pred OCV - R0 * I - x_pred(2); % 4. 卡尔曼增益 dOCV dOCV_dSOC_Lookup(x_pred(1)); % OCV曲线斜率 C [dOCV, -1]; K P_pred * C / (C * P_pred * C R); % 5. 测量更新 x_est x_pred K * (Vt - V_pred); P (eye(2) - K * C) * P_pred; soc_ekf x_est(1); end这段代码里最值得推敲的是C矩阵。C [dOCV, -1]表明端电压对SOC的偏导来自OCV曲线斜率对Vp的偏导恒为-1。dOCV_dSOC_Lookup是OCV斜率查表函数OCV曲线在中间平台区斜率只有0.1到0.3 mV/%两端会陡增。这个非线性正是EKF线性化误差的主要来源SOC在0.5附近时dOCV很小C矩阵第一列就小卡尔曼增益偏向Vp通道SOC的修正幅度变小表现出来就是平台区SOC收敛变慢。这也是后面UKF对比时最明显的差异点。协方差更新那行P (I - KC)P_pred是教科书写法数值稳定性一般实际调试时如果发现协方差慢慢变得不对称建议换成Joseph形式这一条在第5章的踩坑记录里会展开说。3.2 Q/R与P0的设定先物理后微调Q和R设不好EKF要么发散要么迟钝。我习惯先从物理含义给初值R对应电压传感器噪声方差车规级电压采集精度通常在±5mV以内把R设为(0.005)² 2.5e-5即可如果测量通道有明显噪声用一段静置数据算标准差再平方比拍脑袋给值可靠得多。Q矩阵里q_soc对应SOC状态方程的不确定性。SOC递推本质是电流积分电流测量有偏置时误差会线性累积这个不确定性就体现在q_soc上。经验量级是q_soc取1e-8到1e-6q_vp取1e-5到1e-3两者差两到三个数量级是正常的因为SOC是慢变量、模型可信度高Vp是快变量、受参数偏差影响大需要给它更大的过程噪声自由度。P0是协方差初始值SOC初值没把握时直接把P0(1,1)设到0.05以上让滤波器初期有足够自由度去修正SOC。参数建议初值调整方向R2.5e-5 ~ 1e-4电压噪声大就调大q_soc1e-8 ~ 1e-6收敛太慢就调大q_vp1e-5 ~ 1e-3动态响应弱就调大P0diag([0.05, 0.001])初值不确定就调大3.3 mainc.m主循环怎么跑通主脚本mainc.m的结构很固定加载参数和工况数据初始化EKF与UKF各自的状态和协方差然后逐拍调用两个滤波函数最后画SOC真值、两个估计值和误差曲线。核心循环如下% mainc.m 主循环 param load_battery_params(thevenin_params.mat); data load(bbdst_data.mat); % 包含 I, Vt, soc_true N length(data.I); soc_ekf zeros(N,1); soc_ukf zeros(N,1); x_ekf [0.7; 0]; P_ekf diag([0.05, 0.001]); % SOC初值故意给偏 x_ukf x_ekf; P_ukf P_ekf; Q diag([1e-7, 1e-4]); R 2.5e-5; for k 1:N I data.I(k); Vt data.Vt(k); [soc_ekf(k), x_ekf, P_ekf] EKF_Step(x_ekf, P_ekf, I, Vt, Q, R, param); [soc_ukf(k), x_ukf, P_ukf] UKF_Step(x_ukf, P_ukf, I, Vt, Q, R, param); end注意SOC初值被设成0.7而不是真值附近的0.9这是刻意为之目的是观察两种滤波器从偏差初值收敛的速度。跑完把soc_true、soc_ekf、soc_ukf画在同一张图再单独画两条误差曲线。如果两条估算曲线完全重合先检查是不是初值没给偏或者滤波器压根没收到有效量测更新这两种情况我调试时都遇到过。4. UKF实现sigma点传播与双算法结果对比4.1 sigma点生成参数含义与Cholesky分解UKF的核心是用一组确定性采样点近似状态分布。二维状态下生成5个sigma点每个点带均值权重Wm和协方差权重Wc。alpha、beta、kappa三个参数控制采样散布alpha决定sigma点离均值的距离通常取1e-3beta对高斯分布取2最优kappa一般取0。lambda由alpha、kappa和状态维度n共同算出采样代码如下function [X, Wm, Wc] sigma_points(x, P, alpha, beta, kappa) n length(x); lambda alpha^2 * (n kappa) - n; Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2 : 2*n1 Wm(i) 1 / (2 * (n lambda)); Wc(i) 1 / (2 * (n lambda)); end S chol((n lambda) * P, lower); % Cholesky分解P必须正定 X zeros(n, 2*n1); X(:,1) x; for i 1:n X(:,i1) x S(:,i); X(:,in1) x - S(:,i); end endCholesky分解要求P正定这是UKF最容易翻车的点。一旦数值误差让P变成半正定或非正定chol直接报错整个仿真中断。alpha取得太小会让sigma点离均值过近非线性传播后估计的方差偏小滤波器过度自信实际数据里alpha取1e-3是稳妥起点往下调要谨慎。4.2 UKF递推流程没有雅可比的非线性传播UKF不做线性化而是把每个sigma点直接送进非线性的状态方程和测量方程再对传播后的点做加权统计。状态传播是对5个点分别做SOC积分和Vp更新测量传播是把预测点送进OCV查表和端电压方程。完整的UKF_Step实现如下function [soc_ukf, x_est, P] UKF_Step(x_est, P, I, Vt, Q, R, param) n length(x_est); [X, Wm, Wc] sigma_points(x_est, P, 1e-3, 2, 0); dt param.dt; tau param.Rp * param.Cp; % 状态传播每个sigma点过非线性方程 Y zeros(n, 2*n1); for i 1:2*n1 Y(1,i) X(1,i) - I * dt / (param.Cn * 3600); Y(2,i) exp(-dt/tau) * X(2,i) (1 - exp(-dt/tau)) * param.Rp * I; end x_pred sum(Y .* Wm, 2); P_pred (Y - x_pred) * diag(Wc) * (Y - x_pred) Q; % 测量传播sigma点过OCV查表和端电压方程 Z zeros(1, 2*n1); for i 1:2*n1 Z(i) OCV_SOC_Lookup(Y(1,i)) - param.R0 * I - Y(2,i); end z_pred Z * Wm; Pzz (Z - z_pred) * diag(Wc) * (Z - z_pred) R; Pxz (Y - x_pred) * diag(Wc) * (Z - z_pred); % 增益与更新 K Pxz / Pzz; x_est x_pred K * (Vt - z_pred); P P_pred - K * Pzz * K; soc_ukf x_est(1); end对比第3章的EKF代码UKF里没有A矩阵也没有C矩阵所有非线性都在传播里直接完成。EKF在dOCV小的平台区会低估SOC的可观测性导致SOC修正慢UKF的sigma点直接携带了OCV曲线在估计点附近的变化信息天然规避了这个问题。代价是计算量每拍要做2n1次状态传播和测量传播。二维状态5次不算多但换成二阶RC模型变成四维状态就是9次算力受限的嵌入式平台上这个差异要提前评估。4.3 同一组数据下EKF与UKF的差异用资源自带数据和相同Q/R跑完对比能稳定复现三个现象。第一SOC初值给偏时UKF收敛比EKF快直观表现是误差曲线更快收窄原因就是平台区dOCV小导致EKF增益保守。第二稳态误差上EKF在OCV平台区会出现周期性小偏差误差曲线呈锯齿状UKF的误差更平缓。第三计算量上EKF单拍耗时大约只有UKF的一半数据量大的长工况仿真里差距肉眼可见。对比项EKFUKF非线性处理一阶泰勒展开sigma点直接传播雅可比计算需要dOCV/dSOC不需要平台区收敛速度偏慢快单拍计算量低约2n1次传播实现难度低中实际选型时如果控制器算力充足、对收敛速度有要求UKF更省心算力紧张且模型线性度尚可就选EKF。这个结论在这套资源里可以直接用数据验证比空对空讨论有说服力。5. 避坑与排查发散、失配与工况适配的踩坑记录5.1 协方差非正定导致UKF崩溃现象UKF跑着跑着在chol那一行报错提示矩阵必须为正定矩阵有时EKF也会出现SOC突然跳到离谱值然后回不来的情况。 原因协方差矩阵在数值上失去正定性。常见诱因是Q给得过大、长时间没有量测更新导致P_pred被截断误差污染以及协方差更新用了P (I-KC)P_pred这种简化形式它在数值上不保证对称正定。 解决把协方差更新改成Joseph形式P (I-KC)P_pred(I-KC) KRK数值稳定性好很多在每次chol之前加一行P (PP)/2 1e-9*eye(n)强制对称且对角占优。改完之后同样的数据和参数不会再现。5.2 SOC初值给偏导致收敛慢现象数据真值从SOC 0.9开始主脚本里初值设0.5EKF跑了很长时间才追上平台区曲线几乎是一条平线。 原因OCV曲线在中间平台很平dOCV接近0测量方程对SOC的敏感度极低卡尔曼增益把更新量几乎全部分给Vp通道SOC修正非常微弱。 解决用静置端电压反查OCV曲线来初始化SOC。静置足够久的端电压等于OCV查表得到的SOC初值误差能压到5%以内。如果就是要验证大初值误差下的收敛行为把P0(1,1)调到0.1以上给足修正自由度但要做好收敛慢的心理准备。5.3 Q/R调参的玄学边界现象R调得很小估算SOC曲线剧烈抖动看起来跟得很紧但全是噪声R调大曲线平滑但真实SOC变化被抹平动态工况下误差变大。 原因R表征传感器噪声真实电压噪声没那么小R过小等于告诉滤波器测量绝对可信于是每个噪声尖峰都被当成真值修正状态SOC自然抖。R过大的话滤波器干脆不听测量只信模型积分积分误差完全无法校正。 解决先用一段静置数据算电压噪声标准差平方后作为R的基准值。Q按第3章表格给初值然后观察端电压残差残差高频成分多说明R偏小或Q偏大残差低频慢漂说明R偏大或Q偏小。调参要有方向不要同时乱动两个矩阵。5.4 Simulink采样步长与脚本时间不同步现象BBDST_workingcondition.slx导出的电流序列在主脚本里递推SOC积分偏差越来越大同一份数据在不同电脑上跑出不同结果。 原因Simulink模型默认变步长求解器导出到工作空间的电流数据时间间隔不固定后处理脚本却用固定dt做离散化递推两者对不上另外模型输出首尾时间戳和脚本索引对齐也容易出问题。 解决把Simulink求解器改成固定步长步长设置与mainc.m里的dt一致比如0.1s输出端用零阶保持。数据载入后先检查时间向量的差分是否恒等于dt不等就先重采样。第一拍之前补一个静置点避免从第一个采样点就直接开始递推。5.5 电流符号约定不统一现象估算的SOC变化方向和真值相反放电时SOC反而上升或者误差曲线整体反向。 原因电流符号约定在不同模块间没对齐。模型里是放电为正实测数据处理或读取Simulink输出时有人乘了-1导致滤波器收到的电流方向和模型假设相反。 解决在mainc.m加载数据后加一行断言检查前几拍SOC真值变化方向与电流积分方向是否一致不一致就翻转电流符号。这种低级错误靠肉眼盯曲线很难发现加个自动检查一劳永逸。6. 进阶验证残差分析、极限初值与抗噪测试代码跑通只是第一步真正判断一套滤波参数和模型是否可靠我习惯补上三个验证动作。第一个是极限初值测试。把SOC初值分别设成0.1和0.9其他条件不变跑完整段工况。合格的滤波器应该在一段工况内收敛到真值附近收敛时间用于量化初值敏感性。这个测试下EKF和UKF的差距最直观也是写报告时最值得放出的对比图。第二个是端电压残差分析。把Vt_predicted - Vt_measured画出来观察残差均值是否接近0、是否白噪声。残差有偏说明OCV或R0有偏差残差有慢漂说明Rp、Cp参数失配残差高频抖动说明R设小了。残差不干净SOC曲线再平滑都不值得信。我调参时都是先看残差再看SOC误差顺序反了容易被平滑曲线误导。第三个是叠加测量噪声的鲁棒性测试。在端电压上加标准差2mV的高斯白噪声再跑一遍把R相应调整到(0.002)²量级观察SOC误差是否还能稳定在可接受范围。这个测试能暴露滤波器对噪声的过敏程度EKF和UKF的表现差异在噪声叠加后会更明显。从那以后我每次拿到SOC估算代码第一件事就是先跑残差分析残差干净了才开始动Q/R这个习惯帮我挡掉了不少看似准确实则偶然的估算结果。调参和验证的顺序本身就是工程经验的一部分按这三个验证动作走一遍这套资源里EKF和UKF的实际差异会看得非常清楚。希望帮到你。本文还有配套的精品资源点击获取