简介压缩包内是用于电池荷电状态SOC估算的MATLAB实现聚焦扩展卡尔曼滤波EKF与无迹卡尔曼滤波UKF两种算法面向电池管理系统研究人员、电力电子方向学生及从事SOC估算的工程师尤其适用于电动汽车、储能系统等动态工况下的状态估计对比研究。包体共3个文件包含2个m源文件与1个slx仿真模型整体仅105KB轻量易用m文件承载EKF/UKF及Thevenin模型估算主程序slx提供BBDST动态工况仿真环境可直接运行并观察滤波效果。已有247人学习代码结构清晰注释完整下载后能快速在MATLAB中打开通过切换滤波器对比SOC估算轨迹并依据误差曲线调整噪声参数以提升精度。资源小巧却覆盖建模、滤波、仿真全流程是开展电池SOC估计研究的实用参考。1. 第一次在 MATLAB 里同时跑通 EKF 和 UKF 做电池 SOC 估算时我的直觉是 EKF 精度不如 UKF实际数据出来两个滤波器在 NEDC 工况下的 SOC 误差都在 1.5% 以内反而是一次把 OCV 曲线分段拟合的失误让整条 SOC 曲线偏移了近 4%。这个反差说明EKF 和 UKF 在 SOC 估算里的差距远小于模型误差和噪声参数引入的偏差。对做 BMS、储能系统或电池仿真的人来说卡尔曼滤波家族是 SOC 估计绕不开的主流方案。标题里的 EKF扩展卡尔曼滤波和 UKF无迹卡尔曼滤波本质都是同一套“预测 更新”框架差别只在非线性误差怎么处理。这篇不抄公式按在 MATLAB 里从零写 EKF、再升级到 UKF 的顺序把状态方程怎么建、初值怎么给、参数怎么调、坑在哪里说一遍保证你能直接对照自己的数据上手。2. 电池模型先于滤波器SOC 估算精度上限来自状态方程很多新手把注意力全放在滤波器的五个公式上结果跑出来的曲线又抖又偏。我的经验是先标定 OCV-SOC再辨识 RC 参数最后用同一组数据复核模型没问题才轮到滤波器。EKF 和 UKF 再怎么高级也是在模型给的物理约束里做修正。2.1 一阶 RC 和二阶 RC 模型怎么选从离线辨识实验说起电池戴维南等效电路模型里最常用的是一阶 RC 和二阶 RC。一阶 RC 是恒压源 OCV(SOC) 串联内阻 R0再接一个 R1-C1 并联网络二阶 RC 多一组 R2-C2用来描述扩散极化较慢的过程。模型参数数量辨识难度动态精度适用场景一阶 RCR0、R1、C1低脉冲实验一次拟合中动态切换时误差偏大磷酸铁锂、稳态充放电、简单 BMS二阶 RCR0、R1、C1、R2、C2中需要双时间常数拟合高电流突变后收敛更快三元锂电池、车载动态工况、储能调频我的选择习惯先试一阶 RC如果 SOC 误差在动态工况下超过 2%再升到二阶。一阶模型对电压响应描述简单计算量小在嵌入式上跑得动二阶模型精度高但参数辨识容易陷入局部最优需要用不同倍率的脉冲数据一起拟合。参数来源是 HPPC 混合脉冲特性实验。操作不复杂容量标定后用 1C 放电 10s、静置 40s再用 1C 充电 10s然后在不同 SOC 点重复。放电和静置段的电压回弹能同时拟合出 R0、R1、C1。注意OCV-SOC 曲线必须用静置后的电压而不是放电过程中的电压。静置时间不够测出来的 OCV 会带极化电压整个滤波器都会跟着偏。2.2 把模型写成离散状态空间这是 MATLAB 代码的第一块砖以一阶 RC 为例状态向量选 SOC 和极化电压 V_RC系统状态方程连续时间形式d(SOC)/dt -η·I / C_capd(V_RC)/dt -V_RC/(R1·C1) I/C1观测方程V_t OCV(SOC) - V_RC - R0·I离散化后写成递推形式MATLAB 里最好封装成两个函数一个算状态预测一个算观测预测后面 EKF 和 UKF 都共用这一对函数。function x_next stateFcn(x, I, params) % 状态 x [SOC; V_RC] % 输入 I 为电流放电为正 dt params.dt; cap params.capacity * 3600; % Ah 转 As eta params.coulomb_efficiency; tau params.R1 * params.C1; SOC_next x(1) - eta * dt / cap * I; VRC_next exp(-dt / tau) * x(2) params.R1 * (1 - exp(-dt / tau)) * I; x_next [SOC_next; VRC_next]; end这段代码里SOC 更新用的是安时积分的形式但后面滤波器会对它做修正所以不会像纯安时积分那样漂移。V_RC 更新用了一阶 RC 电路的零输入响应加零状态响应等效于反向欧拉离散稳定性比纯前向欧拉好适合大电流突变。观测函数如下function Vt measurementFcn(x, I, params) SOC x(1); ocv polyval(params.ocv_coeff, SOC); Vt ocv - x(2) - params.R0 * I; endOCV 多项式建议用 polyfit 拟合成 5 到 8 次多项式。次数太低拟合不出磷酸铁锂中段的平台次数太高容易在端点出现龙格现象。常见的做法是把 SOC 范围限制在 0 到 1 之间并在两端做夹持否则 polyval 在超出范围时会给出完全离谱的 OCV。2.3 为什么是 EKF 和 UKF安时积分和开路电压法各自卡在哪安时积分的问题是开环初始 SOC 猜错就永远错电流传感器只要有微小偏置误差就随时间累积。开路电压法虽然不积累误差但必须等电池静置足够久动态工况下基本没法用。EKF 和 UKF 做 SOC 估算时把安时积分放在状态预测里再用端电压观测做闭环修正扬长避短。EKF 对非线性函数做一阶泰勒展开计算快但线性化误差大UKF 用一组确定性采样点直接经过非线性函数传播精度更高代价是多了 Sigma 点的计算。我的判断标准很简单如果模型强非线性比如低温大倍率脉冲或者 OCV-SOC 曲线平台很平直接用 UKF如果只做常温下的常规充放电EKF 足够应付而且调参压力小。3. 在 MATLAB 手写 EKF 做 SOC 估算一套能直接落地的初始化与递推循环EKF 的实现分三块初始化、预测、更新。很多流程图把这三步画得花里胡哨实际代码就几十行。最怕的不是公式而是初值和噪声矩阵给得不符合物理意义。3.1 状态初值、协方差 P、噪声矩阵 Q 和 R先用一组能跑的数值%% 电池与采样参数 params.R0 0.05; % 欧姆内阻单位 Ohm params.R1 0.02; % 极化电阻 params.C1 2000; % 极化电容单位 F params.capacity 40; % 电池容量单位 Ah params.dt 1; % 采样周期单位 s params.coulomb_efficiency 1; % 库仑效率 %% 状态向量初值 x [SOC; V_RC] x0 [0.5; 0]; % 先猜 50% SOC极化电压为 0 %% 协方差矩阵与噪声矩阵 P0 diag([0.1^2, 0.5^2]); % SOC 初值不确定度 10%极化电压不确定度 0.5V Q diag([1e-6, 1e-6]); % 过程噪声 R 1e-4; % 端电压测量噪声方差单位 V^2P0 反映的是对初值的信任程度。如果初始 SOC 完全未知把 P0 的第一个对角元设到 0.2^2让滤波器一开始就大胆修正。如果 P0 设得太小滤波器会盲目信任初值收敛会非常慢。Q 和 R 的比例才是 EKF 调参的核心。Q 大代表状态方程不可信卡尔曼增益会偏大估计值容易跟着电压噪声抖R 大代表测量不可靠滤波器会更依赖模型预测响应变慢。工程上 R 可以从传感器数据手册给的标准差平方起步比如 5mV 的噪声对应 R (0.005)^2 2.5e-5。3.2 EKF 预测和更新两步走线性化这一步要特别小心EKF 需要计算两个雅可比矩阵状态转移矩阵 A 和观测矩阵 H。A 是状态函数对状态变量的偏导H 是观测函数对状态变量的偏导。对模型稍加变形A 可以直接写成常数矩阵。function [soc_est, v_est] runEKF(current_data, voltage_data, params, x0, P0, Q, R) dt params.dt; tau params.R1 * params.C1; A [1, 0; 0, exp(-dt / tau)]; n length(x0); x x0; P P0; N length(current_data); soc_est zeros(N, 1); v_est zeros(N, 1); for k 1:N I current_data(k); %% 预测 x_pred stateFcn(x, I, params); P_pred A * P * A Q; %% 更新 Vt_meas voltage_data(k); Vt_pred measurementFcn(x_pred, I, params); % 观测雅可比 H [dOCV/dSOC, -1] dOCVdSOC polyval(polyder(params.ocv_coeff), x_pred(1)); H [dOCVdSOC, -1]; S H * P_pred * H R; K P_pred * H / S; x x_pred K * (Vt_meas - Vt_pred); P (eye(n) - K * H) * P_pred; soc_est(k) x(1); v_est(k) x(2); end endA 矩阵的第一行 [1, 0] 对应 SOC 的预测不依赖极化电压第二行第二列是 RC 网络的自衰减系数。这里直接用了离散化后的精确表达比用expm(Ac*dt)更直观。H 矩阵里的 dOCVdSOC 来自多项式导数这个值在 SOC 平台期接近 0卡尔曼增益会变小滤波器会更多依赖安时积分这是正常现象。平台期出现轻微误差是物理限制不是 bug。3.3 跑一条实际 SOC 曲线代码跑通后先看这几张图soc_ref load(soc_ref.mat).soc_ref; time (0:length(soc_ref)-1) * params.dt; [soc_est, ~] runEKF(current_data, voltage_data, params, x0, P0, Q, R); figure(Color, w); subplot(2,1,1); plot(time, soc_ref, k--, LineWidth, 1.5); hold on; plot(time, soc_est, r-, LineWidth, 1.2); legend(参考SOC, EKF估计); ylabel(SOC); grid on; ylim([0 1]); subplot(2,1,2); plot(time, soc_est - soc_ref, b-, LineWidth, 1); ylabel(误差); xlabel(时间(s)); grid on;跑通后先看四个点SOC 是否单调误差是否在头几个点内快速收敛收敛后是否围绕零线小幅波动放电快结束时有没有明显发散。如果 SOC 曲线总是贴着急停电压掉下去多半是 R0 偏大或者 OCV 曲线末端拟合过头。4. 从 EKF 升级到 UKF无迹变换绕开雅可比矩阵也把非线性误差压小一点EKF 的短板是线性化。OCV-SOC 曲线在磷酸铁锂平台期非常平泰勒展开在这些点上的偏差会反馈到卡尔曼增益里导致滤波性能下降。UKF 不展开非线性函数而是用多个 Sigma 点去“碰”它最后加权统计。4.1 无迹变换到底在做什么用一组 Sigma 点代替一个求导先解释最常见的比例无迹变换。设状态维数为 n均值为 x协方差为 P尺度参数 lambda 的计算方式是lambda alpha^2 * (n kappa) - n其中 alpha 决定 Sigma 点离均值多远通常取 1e-3 到 1。kappa 是次级缩放参数n2 时取 1 即可。然后生成 2n1 个 Sigma 点sigma_0 xsigma_i x sqrt((n lambda) * P)_isigma_{in} x - sqrt((n lambda) * P)_i这组点经过非线性函数后重新加权计算均值和协方差精度至少达到二阶比 EKF 的一阶线性化更稳。4.2 MATLAB 写 UKF 的完整递推采样、权重、更新一步不少function [soc_est, v_est] runUKF(current_data, voltage_data, params, x0, P0, Q, R) n length(x0); alpha 1e-3; kappa 1; beta 2; 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) 0.5 / (n lambda); Wc(i) 0.5 / (n lambda); end x x0; P P0; N length(current_data); soc_est zeros(N, 1); v_est zeros(N, 1); for k 1:N I current_data(k); %% 生成 Sigma 点 sqrtP chol((n lambda) * P, lower); X zeros(n, 2*n1); X(:, 1) x; for i 1:n X(:, i1) x sqrtP(:, i); X(:, ni1) x - sqrtP(:, i); end %% 状态预测 X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:, i) stateFcn(X(:, i), I, params); end x_pred zeros(n, 1); P_pred zeros(n, n); for i 1:2*n1 x_pred x_pred Wm(i) * X_pred(:, i); end for i 1:2*n1 dx X_pred(:, i) - x_pred; P_pred P_pred Wc(i) * (dx * dx); end P_pred P_pred Q; %% 观测预测 Y_pred zeros(2*n1, 1); for i 1:2*n1 Y_pred(i) measurementFcn(X_pred(:, i), I, params); end y_pred 0; for i 1:2*n1 y_pred y_pred Wm(i) * Y_pred(i); end Pxy zeros(n, 1); Pyy 0; for i 1:2*n1 dx X_pred(:, i) - x_pred; dy Y_pred(i) - y_pred; Pxy Pxy Wc(i) * dx * dy; Pyy Pyy Wc(i) * dy * dy; end Pyy Pyy R; %% 更新 K Pxy / Pyy; x x_pred K * (voltage_data(k) - y_pred); P P_pred - K * Pyy * K; soc_est(k) x(1); v_est(k) x(2); end end这段代码里最需要注意的就是 chol 函数。它要求 (n lambda) * P 是正定矩阵一旦 P 因为数值问题变得不正定程序会直接报错。工程上更稳的写法是先做一次对称化P (P P) / 2;4.3 EKF 和 UKF 的对比表精度、耗时、调参难度对比项EKFUKF雅可比矩阵需要且要保证导数连续不需要单步计算量小适合嵌入式实时大约 2 倍于 EKF非线性适应能力弱OCV 平台期误差偏大强能够捕捉二阶矩SOC 初值收敛速度慢一点但平稳通常快 2-3 个采样周期调参难度低重点调 Q/R多出 alpha、beta、kappa 三个参数实际测试中两条曲线的 RMSE 差距通常在 0.2% 到 0.8% 之间并没有网上吹得那么悬殊。UKF 的真正价值在初值错得离谱或者电流突变剧烈的场景它们的优势更容易测出来。rmse_ekf sqrt(mean((soc_est_ekf - soc_ref).^2)); rmse_ukf sqrt(mean((soc_est_ukf - soc_ref).^2)); mae_ukf mean(abs(soc_est_ukf - soc_ref));如果两条曲线差异小于 0.3%不如把时间花在标定 OCV 和辨识 RC 参数上那才是决定 SOC 估算精度的主要矛盾。5. SOC 估算里 EKF 和 UKF 的避坑指南能把结果画成鬼图的几个细节这一章全是血泪经验。很多问题看起来像滤波器写错了但真正原因藏在数据标定和噪声处理里。按现象、原因、解决方案的顺序一条条排查比重新抄一遍公式有用。5.1 P 矩阵发散成 NaN现象程序跑几秒后SOC 突然变成 NaN或者画出来的曲线直接断掉。原因最常见的是状态量单位不一致。比如电流用了 A容量却用了 Ah没有乘 3600导致 SOC 步长超过实际物理范围。另一类是 OCV 多项式在边界外给出极大值H 矩阵出现天文数字协方差很快失去正定性。解决方案在递推循环里加一道防御if any(~isfinite(P_pred(:))) || any(~isfinite(x_pred(:))) error(数值发散检查单位换算和 H 矩阵); end另外每次更新后把 P 强制对称化P (P P) / 2能挡掉大部分数值问题。5.2 OCV-SOC 标定不准再好的滤波器也救不回来现象SOC 估计值在某个区间反复横跳尤其是磷酸铁锂的 20% 到 80% 平台期。原因OCV-SOC 曲线用分段多项式拟合段与段之间导数不连续H 矩阵这里的 dOCVdSOC 会跳变导致卡尔曼增益突变。解决方案不要用两段独立多项式改用一条完整的高次多项式或者用平滑样条。polyfit 前先把 SOC 归一化到 [0, 1]拟合后检查各点残差。对于平台期宁可让多项式在中间平一点也不要为了拟合噪声点增加局部抖动。提示OCV 标定时每个 SOC 点静置时间至少 1 到 2 小时否则残余极化会让平台期被抬高。5.3 Q 和 R 的物理意义没对上现象SOC 曲线要么迟钝得不像话要么抖得跟白噪声一样。原因Q 和 R 的比例失当。Q 设得太大滤波器疯狂相信电压测量电流突变时 SOC 会跟着电压一起剧烈抖动R 设得太大滤波器几乎忽略测量退化成纯安时积分。先说个落地习惯R 首先由传感器噪声决定电压采集板的静音噪声是多少就用多少。Q 则要看模型误差量级先用 Q diag([1e-6, 1e-6])让滤波器正常收敛再缓慢调整。判断标准很简单把误差曲线画出来如果误差波动幅度接近传感器噪声说明 R 偏小。如果误差收敛到某个偏置后长期不消失说明 Q 偏小。调参没有玄学就是看误差曲线上噪声和偏置哪个更碍眼。5.4 采样时间不同步数据预处理里埋的雷现象SOC 估计在电流突变的瞬间出现一个明显尖峰然后慢慢拉回来。原因电压和电流两列数据时间基准不一致或者用了重采样但没对齐。滤波器会用当前时刻的电压去更新上一时刻的状态相当于把未来信息混进来表现就是超前或尖峰。解决方案在递推前先统一时间轴。常见做法是time (0:length(voltage_data)-1) * dt; current_data interp1(time_current, current_raw, time, linear, extrap);除非万不得已别用二次插值。线性插值在动态工况下已经够用高阶插值反而会在电流跳变处引入过冲。5.5 只做了常温标定低温工况直接翻车现象常温下 SOC 误差 1%拿到低温或者大倍率脉冲下误差直接到 5% 以上。原因R0、R1、C1 和 OCV 都是温度的函数尤其低温下欧姆内阻和极化阻抗会显著增大常温参数在低温下完全失真。解决方案至少做 0 度、25 度、45 度三组 HPPC 标定把参数表写成温度插值。如果硬件资源有限可以先在 MATLAB 里用 lookup table 离线做验证完再移植到嵌入式。注意温度对 EKF 和 UKF 的影响是等价的不要指望换 UKF 就能解决模型参数漂移。先把参数表补全再谈滤波算法。6. 把“能跑”变成“可信”我的 SOC 估计验证顺序和现在坚持的一条习惯代码跑出一条光滑的 SOC 曲线只是起点。我见过最多的翻车是换一组数据或者换一个初值后估计结果就面目全非。所以我现在不会拿单条曲线当验收标准而是固定的三步验证。第一步是初值扰动测试。分别用 0.8、0.5、0.2 三个 SOC 初值跑同一段工况看滤波器能否在 200 秒内收敛到同一参考值。init_list [0.8, 0.5, 0.2]; for i 1:length(init_list) x0 [init_list(i); 0]; [soc_est, ~] runUKF(current_data, voltage_data, params, x0, P0, Q, R); rmse(i) sqrt(mean((soc_est - soc_ref).^2)); end如果三条曲线最终都收敛说明初值依赖在可控范围。如果有一条收敛到错误平台问题多半出在 OCV 曲线拟合而不是滤波器本身。第二步是工况覆盖测试。恒流放电、动态工况、间歇静置各跑一组记录 RMSE 和最大绝对误差顺便检查 SOC 末端是否出现提前跳零或冲过 100% 的现象。第三步是传感器偏置注入。把电流信号故意加一个 0.01C 的偏置看 SOC 估计能否依旧保持在真实曲线附近。这一步直接检验滤波器修正能力纯安时积分在这里会明显漂移。我现在坚持的习惯是每次改完模型参数先把 EKF 和 UKF 在同一个数据集上同时跑一遍把两条误差曲线叠在一张图里保存成基线。以后再改任何东西只对比基线的 RMSE 和最大误差。这个习惯让我少走了很多弯路也避免在“看起来更高级”的算法上浪费无谓的时间。希望帮到你。本文还有配套的精品资源点击获取