最近我干了一件挺磨人的事:把《基于二阶EKF的锂离子电池SOC估计的建模与仿真》这篇论文的Simulink仿真部分从头到尾复现了一遍。复现论文,尤其是电池SOC估计这种强非线性系统,难点从来不是把公式抄一遍,而是把等效电路模型、滤波算法、离散化方式、噪声参数全部放进同一套仿真环境里,让它们老老实实配合工作。
这篇论文的核心很明确:用二阶等效电路模型描述锂离子电池的动态特性,再通过二阶扩展卡尔曼滤波(EKF)算法,在电流和端电压测量存在噪声的情况下实时估计电池的荷电状态(SOC)。SOC是电池管理系统(BMS)里的核心状态,直接影响到续驶里程估算、充放电策略和均衡控制。它没法直接测量,只能靠模型和算法去"猜",所以这种仿真复现具有很强的工程参考价值。
这篇文章适合正在做电池SOC估计、BMS软件开发、论文复现或者课程设计的朋友。我尽量把从零开始的建模思路、Simulink里的模块摆放、MATLAB Function的代码封装、各种坑和调试心得都写清楚。有些细节论文不会写,但复现时必须自己补上,这些反而是最花时间的地方。
1. 项目解读与复现思路
1.1 这篇论文到底想解决什么问题
锂电池SOC估计最朴素的方法是安时积分,也就是把电流对时间积分,再用总容量折算SOC。看起来简单,但电流传感器有零漂和噪声,积分误差会随着时间不断累积,而且初始SOC必须已知。另一个方法是开路电压法,认为SOC和开路电压OCV有确定关系,测到OCV就能查表得到SOC。但OCV需要电池静置很久才能稳定,车跑在路上根本不可能频繁静置。
所以论文采用EKF的思路:把SOC作为状态变量,用一个等效电路模型去预测端电压,再把预测电压和实际测量电压的误差反馈给滤波器,修正SOC估计。这样既利用了模型的动态信息,又实时融合了电压测量,误差不会像安时积分那样无限累积,而且不需要频繁静置。仿真部分要验证的就是这个方法在典型工况下能不能把SOC误差控制在可接受的范围内。
1.2 为什么用"二阶EKF"而不是普通EKF
这里要先把概念捋清楚。论文标题里的"二阶"体现在两个地方:一是电池采用了二阶RC等效电路模型,也就是两个RC并联网络串联在欧姆内阻后面,用来分别模拟电池的快速极化响应和慢速扩散响应;二是滤波算法本身用了二阶Taylor展开的EKF,而不是只保留一阶项的经典EKF。
普通EKF在非线性不强的系统里表现不错,但锂离子电池的OCV-SOC曲线在低SOC区间变化很剧烈,比如SOC从0.2降到0.1时,OCV可能快速跌落,此时只保留一阶导数会带来明显的线性化误差。二阶EKF在算法推导里额外补偿了Hessian矩阵的影响,相当于在预测和更新时把曲线的弯曲程度也考虑进去,所以对强非线性区段的估计精度会更友好。
当然,计算代价也上去了。标准EKF只需要计算一阶雅可比矩阵,二阶EKF还要处理Hessian项,调试的时候也更容易因为矩阵运算出问题。但从仿真角度来说,MATLAB Function完全扛得住,复现这件事并没有想象中那么复杂。
1.3 整体复现框架怎么搭
我的复现方案分三块,全部放在一个Simulink模型里:
- 真实电池模型:用积分器和基本数学模块搭建二阶RC等效电路,模拟电池实际的SOC、极化电压和端电压输出,并给端电压叠加一个高斯白噪声,作为"测量值"。
- EKF估计模块:用一个MATLAB Function封装二阶EKF算法,接收电流、测量电压和时间步长,输出SOC估计值、两个极化电压估计值和协方差矩阵。
- 结果对比模块:把真实SOC和EKF估计SOC一起送到Scope,同时计算误差、RMSE等指标。
这种结构的好处是把"被估计对象"和"估计算法"完全分离。EKF模块只依赖输入的电压电流,不依赖真实模型内部的任何状态,这样验证出来的结果才有说服力:如果算法能从带噪声的电压里把SOC捞出来,说明方法本身是成立的。
2. 锂电池等效电路模型与参数处理
2.1 二阶RC等效电路模型的方程推导
二阶RC等效电路模型(也叫二阶Thevenin模型)结构很直观:正极到负极之间串一个欧姆内阻R0,然后再串两个并联的RC网络。第一个RC网络(R1、C1)描述电化学极化的快速动态,时间常数通常在几秒到几十秒;第二个RC网络(R2、C2)描述浓度极化的慢速动态,时间常数可能到几十秒甚至几分钟。
状态向量取三个变量:SOC、V1(R1C1两端的极化电压)、V2(R2C2两端的极化电压)。输入是电流I,这里约定放电为正。连续时间状态方程为:
d(SOC)/dt = -I / Q
其中Q是电池总容量,单位要换算成安秒(Ah×3600)。如果Q给的是Ah,电流单位是A,时间单位是秒,就一定要除以3600,不然SOC会以奇怪的速度飞掉,这是新手最容易忽略的点。d(V1)/dt = -V1/(R1×C1) + I/C1
d(V2)/dt = -V2/(R2×C2) + I/C2
输出方程是端电压:
Vt = OCV(SOC) - V1 - V2 - I×R0
这里OCV(SOC)是开路电压关于SOC的静态函数,通常通过HPPC实验标定,论文里一般给一条曲线或一个多项式。复现时可以直接用一个Lookup Table或者多项式函数模块。
2.2 OCV-SOC曲线与参数获取
OCV-SOC曲线是整个模型的灵魂。论文里通常会给出一个测试数据表,但复现时不一定能拿到原始数据。我的做法是用一个五阶或六阶多项式拟合典型曲线,或者直接用MATLAB的Lookup Table模块逐点查表。
这里有个重要细节:EKF算法里需要计算OCV对SOC的导数,也就是雅可比矩阵里的元素。如果直接用一维查表模块,相邻数据点之间是线性插值,导数是不连续的,在查表节点附近会产生锯齿状抖动,导致SOC估计值高频跳动。所以如果条件允许,我建议用光滑的多项式拟合,离线求好导数公式,再写进EKF代码里。比如我们假设OCV和SOC的关系满足:
OCV = a0 + a1×SOC + a2×SOC² + a3×SOC³ + a4×SOC⁴ + a5×SOC⁵
那导数就是:
dOCV/dSOC = a1 + 2a2×SOC + 3a3×SOC² + 4a4×SOC³ + 5a5×SOC⁴
二阶导就是:
d²OCV/dSOC² = 2a2 + 6a3×SOC + 12a4×SOC² + 20a5×SOC³
把这个导数函数直接内嵌进MATLAB Function里,EKF的H矩阵就是连续的,滤波也会平稳很多。
2.3 初始参数设定与模型验证
论文里未必会把所有参数贴出来,复现时我们只能先给一组合理的典型值,再根据工况微调。我用的这样一组参数:
- R0 = 0.02 Ω
- R1 = 0.005 Ω,C1 = 2000 F,时间常数约10秒
- R2 = 0.01 Ω,C2 = 10000 F,时间常数约100秒
- 额定容量 Q = 5 Ah
这些参数大致符合一个中等容量软包电池的特性。搭建完真实电池模型后,先用一个脉冲电流工况跑一遍,看端电压响应是不是"先瞬跳、再快速恢复、随后缓慢爬升"的典型形态。如果形态不对,大概率是R0、R1C1、R2C2的数值比例出了问题,而不是算法问题。
记住,模型参数和EKF里的状态转移矩阵必须一致。你的EKF里写了几组RC,真实模型就必须搭几组RC。如果真模型是二阶,EKF用一阶模型去估计,不管Q和R怎么调都会存在稳态误差,这种错误在复现时还挺容易犯的。
3. 二阶EKF算法设计与Simulink实现
3.1 从标准EKF到二阶EKF到底改了什么
经典EKF在预测和更新时,把非线性函数f(x)和h(x)只做一阶Taylor展开,也就是只使用雅可比矩阵F和H。二阶EKF则在展开式里保留二次项,并且用协方差矩阵P对方差进行修正。
具体到SOC估计这个场景,状态方程其实大部分是线性的:SOC变化是电流的线性积分,极化电压的动态也是线性一阶系统。真正非线性的是量测方程,因为OVC和SOC的关系是非线性的。所以我在实现时只把量测更新的二阶项补了进去,这样既不影响计算效率,又踩中了非线性最强的点。
量测更新的核心思想是这样:如果我们用当前预测状态x_pred算出来一个h(x_pred),直接把测量值和它做差,相当于忽略了OCV曲线弯曲带来的偏差。二阶修正项就是在h(x_pred)上加上0.5 × trace(H_xx × P_pred),其中H_xx是量测函数对状态的Hessian矩阵,P_pred是预测协方差矩阵。这样预测值更接近非线性映射后的真实统计均值,再算innovation也就更准。
3.2 状态方程离散化与矩阵计算
连续状态方程不能直接在离散滤波里用。Simulink里的MATLAB Function是离散采样执行的,每个步长dt要做一次预测和一次更新。状态转移矩阵A要按步长dt进行离散化。
放电电流下,离散状态方程为:
SOC(k) = SOC(k-1) - dt/(Q×3600) × I
V1(k) = exp(-dt/(R1C1)) × V1(k-1) + R1×(1-exp(-dt/(R1C1))) × I
V2(k) = exp(-dt/(R2C2)) × V2(k-1) + R2×(1-exp(-dt/(R2C2))) × I
写成矩阵形式就是:
A = [1, 0, 0;
0, exp(-dt/(R1C1)), 0;
0, 0, exp(-dt/(R2C2))]
B = [-dt/(Q×3600);
R1×(1-exp(-dt/(R1C1)));
R2×(1-exp(-dt/(R2C2)))]
这里的dt就是Simulink的采样周期,我一般设0.1秒。如果dt很大,欧拉离散会不稳,必须用精确指数离散,也就是上面带exp的写法,这会显著减少数值发散的概率。
3.3 在Simulink里用MATLAB Function封装EKF
Simulink里用MATLAB Function是最省事的方式,因为你在同一个模块里既能写算法,又能直接访问端口数据。我新建了一个函数块,输入是电流I、测量电压Vt、步长dt,输出是SOC_est、V1_est、V2_est。状态变量x和协方差P用persistent变量保存,保证每个采样周期能记住上一拍的值。
核心代码如下,你可以直接抄:
function [soc, v1, v2] = ekf_2nd_order(u, Vt, dt) % 二阶EKF量测更新实现(状态方程线性,量测方程非线性) persistent x P Q R if isempty(x) x = [0.6; 0; 0]; % 初始SOC假设60% P = diag([0.05, 0.001, 0.001]); % 状态协方差 Q = diag([2e-4, 1e-5, 1e-5]); % 过程噪声 R = 0.01; % 电压测量噪声方差 end % 电池参数 R0 = 0.02; R1 = 0.005; C1 = 2000; R2 = 0.01; C2 = 10000; Qnom = 5; % Ah % 状态转移矩阵 A = [1, 0, 0; 0, exp(-dt/(R1*C1)), 0; 0, 0, exp(-dt/(R2*C2))]; B = [-dt/(Qnom*3600); R1*(1-exp(-dt/(R1*C1))); R2*(1-exp(-dt/(R2*C2)))]; % 预测 x_pred = A*x + B*u; P_pred = A*P*A' + Q; % 量测预测与非线性量测函数 soc_pred = x_pred(1); h = ocv(soc_pred) - x_pred(2) - x_pred(3) - R0*u; H = [ocv_d(soc_pred), -1, -1]; % 二阶修正项:H_xx = diag(ocv_dd(soc_pred),0,0) innovation = Vt - (h + 0.5*ocv_dd(soc_pred)*P_pred(1,1)); % 卡尔曼增益 S = H*P_pred*H' + R; K = P_pred*H'/S; % 更新 x = x_pred + K*innovation; P = (eye(3) - K*H)*P_pred; soc = x(1); v1 = x(2); v2 = x(3); end function y = ocv(soc) % 这里用多项式拟合的OCV曲线示例 p = [5.2, -12.4, 9.8, -2.6, 0.6, 3.3]; % 随意一组拟合系数 y = polyval(p, soc); end function y = ocv_d(soc) p = [5.2, -12.4, 9.8, -2.6, 0.6, 3.3]; d = polyder(p); y = polyval(d, soc); end function y = ocv_dd(soc) p = [5.2, -12.4, 9.8, -2.6, 0.6, 3.3]; d = polyder(polyder(p)); y = polyval(d, soc); end注意这个代码里Q和R我给的是初始值,实际要调。function块里最好显式指定采样时间,并且把输出端口加一个Unit Delay或者在模型层级结构上避免直接形成代数环,否则Simulink会报错或者仿真极慢。
3.4 噪声协方差Q和R的调试心得
EKF结果的好坏,一半靠模型,一半靠Q和R。Q矩阵代表你对过程模型的信任程度,R代表你对电压测量传感器的信任程度。
Q太小,滤波器会过度信任模型,SOC更新得很慢,一旦真实情况和模型有偏差,估计值会滞后。Q太大,滤波器会过度信任电压测量,SOC估计会跟着噪声剧烈抖动,波形像一个毛刺很多的锯齿。
R的作用正好反过来。R太大,电压测量几乎不被信任,滤波器基本跟着模型走,SOC会平顺但可能偏离真值;R太小,电压测量的权重过大,SOC会尖峰频出。
我的调试顺序是:先固定R在0.01左右,从小到大调Q,每一步都看SOC误差曲线。如果发现SOC在纯I电流工况下还有明显的高频抖动,就把Q的第一个元素降一个数量级。如果是恒流放电阶段误差逐渐拉大,则说明Q的最后一个元素或者模型参数有偏差,要回头检查参数而不是继续增大Q。
4. 仿真流程搭建与结果分析
4.1 Simulink模型的分层搭法
整个模型我习惯分成四层:信号源层、真实电池模型层、EKF估计层、结果显示层。
信号源层用一个Signal Builder或者Step模块生成电流波形。真实电池模型层用积分器搭建三个状态方程,OCV用Lookup Table或者MATLAB Function实现。端电压输出端加一个Band-Limited White Noise模块,把噪声叠加到电压上。EKF估计层就是上一节那个MATLAB Function,直接接收电流和带噪声电压。结果显示层用几个Scope同时看SOC对比、端电压对比、极化电压估计和协方差对角线。
不建议把所有东西塞在一个个孤立的模块里来回接线,那样后续改参数太痛苦。用Simulink的Subsystem分成模块,每个Subsystem里可以再用Goto/From标签传递信号,这样主线很清楚。
4.2 真实SOC基准怎么生成
既然要验证EKF估计得准不准,就得有一个"标准答案"。在仿真环境里,真实SOC可以直接从模型里的SOC积分器拉出来,也就是用安时积分法算出来的值。虽然安时积分在真实系统里有累积误差,但仿真里所有状态都是精确已知的,我们拿它当真值没有问题。
这里提醒一下,真实模型里的SOC积分器和EKF里的SOC状态是两套独立的东西,千万别把EKF的输出直接反馈回真实模型,否则就是自己考自己,失去验证意义。真实模型的SOC只受电流和历史状态影响,EKF只能通过端电压观测间接知道SOC,这样才能体现"估计"的意义。
4.3 测试工况怎么选
论文用来验证的工况一般是恒流放电、脉冲放电或者动态工况。仿真复现时我建议从最简单的恒流放电开始,然后再上脉冲。
恒流1C放电是最容易诊断问题的工况:如果SOC从80%放20分钟到60%,EKF估计的SOC应该是一条平滑的直线,误差在几个百分点以内。如果连这种工况都稳不住,说明算法或参数一定有问题。
脉冲工况更有意思:一个电流脉冲加上去,端电压瞬间跌落,然后随着极化电压建立继续缓慢下降;电流撤掉后端电压瞬间反弹,再缓慢爬升。EKF必须在这一整套动态过程里依然保持SOC误差不发散。如果第一个脉冲刚结束时SOC估计值出现尖峰,大概率是量测方程二阶修正项没写对,或者Hessian项的符号有问题。
动态工况可以模仿城市循环里的随机加减速,电流随机变化且高频切换。这种工况对EKF是一个更全面的压力测试,误差一般会比恒流工况大一些,但均值仍然应该稳定在真值附近。
4.4 结果指标怎么算
仿真结果不能只靠眼睛看Scope,最好用To Workspace把数据存到MATLAB工作区里,再算定量指标:
- RMSE(均方根误差)反映了SOC估计的整体偏差水平;
- MAE(平均绝对误差)反映平均偏差大小;
- 最大误差则能抓到滤波器在某个工况突变时刻可能出现的瞬态大偏差。
我在复现时要求自己:恒流工况RMSE小于0.5%,脉冲工况RMSE小于1%,动态工况RMSE小于2%,最大瞬态误差不超过5%。如果这个指标要求没达成,就要回头检查模型参数、Q/R 或者二阶修正项,而不是继续改仿真时间。
5. 踩过的坑与复现经验总结
5.1 数值发散与离散化步长的坑
第一次跑模型时,我用最简单的欧拉离散,步长设成1秒,结果SOC估计第一分钟就冲到0%以下去了。后来发现问题出在状态转移矩阵:欧拉离散在步长超过最小时间常数的一半时,数值稳定性就保不住了。电池里R1C1时间常数只有几秒,步长1秒显然太大了。
解决办法有两个:要么把采样步长缩小到0.01秒,要么用精确指数离散。我后来把A矩阵写成exp(-dt/(R*C))这种形式,步长0.1秒也完全稳得住。这就是我在代码里为什么用了指数项而不是1-dt/(RC)。
5.2 初始协方差和初值带来的收敛问题
EKF的一大优点是初始SOC可以不准,但收敛时间取决于初始协方差P0。如果P0设得太小,说明你一开始对初值非常自信,滤波器就不太愿意去纠正SOC,导致收敛很慢。我一开始把P0(1,1)设为0.001,结果SOC初值设60%,真值是80%,跑了200秒都没拉回来。
把P0(1,1)改成0.05之后,滤波器大概几十秒就能收敛。如果你希望初始SOC很快跟上,P0可以甚至设成0.1,代价是最开始一段时间的SOC估计波动会比较大。工程上一般取P0等于初始不确定性的平方。
5.3 OCV查表求导不连续导致的抖动
之前提到过,不要直接用Lookup Table做EKF里的OCV导数。因为查表是线性插值,导数在数据点之间是折线段,在端点处会跳变。这个跳变传到雅可比矩阵里,卡尔曼增益就会突跳,SOC曲线出现肉眼可见的锯齿。
后来我把OCV数据用polyfit拟合成多项式,然后在代码里直接用polyder求导数,这个问题就彻底消失了。如果实在要用查表,那就用三次样条插值,并在代码里用样条的解析导数,也能缓解抖动。
5.4 给准备复现的人的建议
复现类似论文,我建议严格按照下面的顺序来:
- 先把二阶RC模型在Simulink里跑通,把OCV曲线和电流波形做好,确认端电压形态正常。
- 不要让EKF一步到位,先写一个标准一阶EKF,跑通闭环,确认能估计SOC。
- 在标准EKF基础上加入二阶修正项,对比一阶和二阶在低SOC区间的差异。
- 最后再调Q和R,做定量指标统计。
如果标准EKF都收敛不了,千万不要继续加二阶项,那只会让问题更复杂。
我个人实际复现下来最受益的一点是:把"真实模型"和"估计模型"彻底分开之后,很多算法问题都能在模型内部自检出来。如果你发现无论怎么调参SOC都偏,先怀疑自己的OCV函数写错了,再怀疑状态转移矩阵,最后才是EKF算法本身。滤波器本身只是一个递推最小二乘,真正的"主观偏见"都藏在模型里。
这个仿真做完之后,还可以继续往两个方向扩展:一是把OCV曲线换成实测电池数据,二是把二阶EKF的代码封装成S-Function,拿到其他电池管理系统硬件平台上做快速原型验证。SOC估计这行最后的落脚点永远不是漂亮的仿真曲线,而是算法在真实电池上能不能扛住温度、老化、电流突变这些不确定性。仿真是第一步,但绝对是最值得打好基础的一步。