简介这份资源面向电池管理系统BMS方向的学习者与工程人员提供基于扩展卡尔曼滤波器EKF与二阶RC等效电路模型的电池SOC预测仿真方案适合已掌握MATLAB基础、希望深入理解SOC估算算法的中高级读者。压缩包共3个文件包含2个mat实验数据文件与1个m脚本文件整体约23KB其中mat文件存放不同SOC状态下的电压与SOC对应数据m脚本则完整实现二阶RC模型的动态方程、观测方程以及EKF的初值、噪声设定与状态更新流程。目前已有4218人学习下载说明该仿真在同类资源中具有较高参考热度。读者可借助这套代码与数据直观理解EKF如何对非线性电池模型进行线性化处理掌握二阶RC模型参数拟合与SOC递推估计的完整链路并在此基础上尝试模型校准、算法移植或实时性优化为BMS算法研究与工程实践提供可复用的仿真起点。1. 从一条“电压对不上”的曲线说起EKF 二阶 RC 模型到底在算什么如果你在电池 SOC 估算这个方向待过一阵大概率见过这样的场景同一组 HPPC 工况数据有人用安时积分跑出来 SOC 一路飘有人套了个一阶 RC 加扩展卡尔曼滤波末端误差还是能到 5% 以上。问题往往不在滤波算法本身而在电池模型——一阶 RC 把极化过程压成了一根指数衰减而真实电芯在放电中后期存在明显的双时间尺度极化电化学极化快、浓差极化慢。二阶 RC 模型就是给这两个过程各配一组 RC 并联环节再用扩展卡尔曼滤波EKF在非线性观测方程下做状态递推。这套方案解决的核心诉求很具体在嵌入式算力有限、电流采样有噪声、初始 SOC 不确定的条件下把 SOC 估算误差压到 3% 以内并且能在线修正。它适合做 BMS 算法验证的工程师、做储能 EMS 策略的研究生以及需要给仿真平台搭一套可信电池模型的开发者。往下我会把建模、参数辨识、EKF 离散化、仿真验证和踩坑点一层层拆开代码全部用 Python 写能直接跑。2. 二阶 RC 等效电路建模从电芯到状态空间方程2.1 为什么不是一阶 RC也不是电化学模型选模型本质是在精度和算力之间找平衡点。电化学 P2D 模型精度高但偏微分方程组在单片机上根本跑不动一阶 RC 只有一个极化时间常数在脉冲工况下能凑合但遇到连续变电流或者低温大倍率放电端电压拟合残差会明显变大。二阶 RC 多了一组 RC能同时刻画秒级和分钟级的极化恢复过程状态量从 2 个变成 4 个SOC、两个极化电压EKF 的雅可比矩阵规模仍然很小算力开销对 Cortex-M4 级别的主控是可接受的。我一般会这样判断如果你的工况以恒流为主、温度稳定一阶 RC 够用如果要做动态工况、要在线辨识、要覆盖低温直接上二阶 RC省得后面返工。2.2 电路方程与状态空间形式二阶 RC 等效电路的结构是开路电压源 Uoc(SOC) 串联欧姆内阻 R0再串联两组并联 RCR1//C1、R2//C2。端电压方程U Uoc(SOC) - I·R0 - U1 - U2两个极化电压的动态dU1/dt -U1/(R1·C1) I/C1 dU2/dt -U2/(R2·C2) I/C2SOC 的动态就是安时积分dSOC/dt -I/(Qn·3600)把这三个微分方程写成状态空间状态向量 x [SOC, U1, U2]^T输入 u I输出 y U。这就是后面 EKF 要离散化的连续模型。2.3 用 Python 搭出模型骨架先把模型写成可复用的类参数后面辨识出来直接塞进去。import numpy as np class SecondOrderRC: 二阶 RC 等效电路模型输入电流约定放电为正 def __init__(self, R0, R1, C1, R2, C2, Qn, ocv_func): self.R0 R0 self.R1, self.C1 R1, C1 self.R2, self.C2 R2, C2 self.Qn Qn # 额定容量单位 Ah self.ocv_func ocv_func # OCV-SOC 查表函数 def step(self, soc, u1, u2, I, dt): # 欧姆压降与极化压降 tau1 self.R1 * self.C1 tau2 self.R2 * self.C2 # 一阶后向差分离散化保证数值稳定 u1_new u1 * np.exp(-dt / tau1) self.R1 * I * (1 - np.exp(-dt / tau1)) u2_new u2 * np.exp(-dt / tau2) self.R2 * I * (1 - np.exp(-dt / tau2)) soc_new soc - I * dt / (self.Qn * 3600.0) u_ocv self.ocv_func(soc_new) u_term u_ocv - I * self.R0 - u1_new - u2_new return soc_new, u1_new, u2_new, u_term这段代码里有两个关键点。第一极化电压的离散化用的是精确解形式指数项比前向欧拉稳定得多dt 稍大也不会发散。第二SOC 更新用放电为正的约定如果你拿到的数据充电为正记得在入口统一符号否则 SOC 会反向跑。参数 dt 建议和实际采样周期一致常见是 0.1s 或 1s不要为了“看起来平滑”去插值改 dt会引入额外误差。3. 参数辨识OCV-SOC 曲线和 RC 参数怎么拿3.1 OCV-SOC 曲线的获取与拟合OCV-SOC 是模型的基准它错了后面全错。常见做法是小电流0.02C~0.05C间歇放电每次放 5% SOC 静置 2 小时以上记录静置末端电压作为该 SOC 点的 OCV。如果条件有限也可以用低倍率充放电的平均值近似但精度会差一些。拿到数据点后用单调三次样条或者分段线性插值建表。我一般用 numpy 的 interp1d 做线性插值简单可靠不会出现样条过冲。from scipy.interpolate import interp1d # soc_pts 从 0 到 1 递增ocv_pts 是对应静置电压 soc_pts np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) ocv_pts np.array([3.00, 3.45, 3.55, 3.60, 3.65, 3.70, 3.78, 3.88, 4.00, 4.12, 4.20]) ocv_func interp1d(soc_pts, ocv_pts, kindlinear, fill_valueextrapolate)注意 OCV 曲线必须单调递增如果辨识出来有回环说明静置时间不够或者温度变化太大要重新测。另外 SOC 点不要只取 10 个中段电压平台区建议加密到 5% 一个点否则 EKF 在平台区会“感觉不到”SOC 变化。3.2 用脉冲工况数据辨识 R0、R1、C1、R2、C2R0 从脉冲瞬间的电压跳变直接算R0 ΔU / ΔI。两个 RC 参数用脉冲结束后的电压弛豫曲线拟合常见做法是先用最小二乘拟合双指数函数再反解出 R1、C1、R2、C2。from scipy.optimize import curve_fit def dual_exp(t, u_inf, a1, tau1, a2, tau2): return u_inf - a1 * np.exp(-t / tau1) - a2 * np.exp(-t / tau2) # t_relax 是弛豫时间序列u_relax 是对应端电压 popt, _ curve_fit(dual_exp, t_relax, u_relax, p0[u_relax[-1], 0.05, 10, 0.05, 200], bounds([0, 0, 0.1, 0, 1], [5, 1, 100, 1, 2000])) u_inf, a1, tau1, a2, tau2 popt R1, C1 a1 / I_pulse, tau1 / (a1 / I_pulse) R2, C2 a2 / I_pulse, tau2 / (a2 / I_pulse)这里 p0 的初值很关键tau1 给 10 秒左右、tau2 给 200 秒左右符合大多数锂电芯的双极化时间尺度。如果拟合不收敛先检查弛豫数据有没有去掉前几个毫秒的采样毛刺。辨识出来的参数建议按 SOC 分段存表因为 R0 和极化参数随 SOC 变化明显尤其是两端。3.3 参数随 SOC 和温度变化的处理策略工程上不会只辨识一组参数。常见做法是按 SOC 每 10% 一个区间、温度每 10℃ 一个区间做二维查表EKF 运行时根据当前 SOC 和温度插值取参。如果算力紧张至少把 R0 做成 SOC 的函数因为欧姆压降对端电压影响最直接。提示参数表不要用高阶多项式去拟合外推到 SOC 两端容易出负值查表加线性插值最稳。4. EKF 离散化与仿真闭环从方程到能跑的代码4.1 状态转移矩阵和雅可比矩阵的推导把第 2 章的连续模型离散化状态转移用精确解形式SOC(k1) SOC(k) - I(k)·dt/(Qn·3600) U1(k1) U1(k)·exp(-dt/τ1) R1·I(k)·(1-exp(-dt/τ1)) U2(k1) U2(k)·exp(-dt/τ2) R2·I(k)·(1-exp(-dt/τ2))状态转移矩阵 F 对状态求偏导SOC 对 SOC 的偏导是 1极化电压对自身的偏导是指数项其余交叉项为 0。观测方程 y Uoc(SOC) - I·R0 - U1 - U2观测矩阵 H 对 SOC 的偏导是 dUoc/dSOC对 U1、U2 的偏导是 -1。dUoc/dSOC 用 OCV 表的数值差分算不要解析求导查表函数没有解析式。4.2 EKF 的五个方程和代码实现EKF 的标准流程是预测状态、预测协方差、算增益、更新状态、更新协方差。下面是一个可直接调用的实现。class EKF: def __init__(self, model, soc01.0, P0None, QNone, RNone): self.m model self.x np.array([soc0, 0.0, 0.0]) # [SOC, U1, U2] self.P np.diag([1e-3, 1e-4, 1e-4]) if P0 is None else P0 # 过程噪声SOC 通道小极化通道稍大 self.Q np.diag([1e-8, 1e-6, 1e-6]) if Q is None else Q self.R np.array([[1e-4]]) if R is None else R # 电压测量噪声 def predict(self, I, dt): tau1 self.m.R1 * self.m.C1 tau2 self.m.R2 * self.m.C2 a1 np.exp(-dt / tau1) a2 np.exp(-dt / tau2) # 状态预测 soc, u1, u2 self.x soc_p soc - I * dt / (self.m.Qn * 3600.0) u1_p u1 * a1 self.m.R1 * I * (1 - a1) u2_p u2 * a2 self.m.R2 * I * (1 - a2) self.x np.array([soc_p, u1_p, u2_p]) # 状态转移矩阵 F np.array([[1, 0, 0], [0, a1, 0], [0, 0, a2]]) self.P F self.P F.T self.Q return F def update(self, I, u_meas): soc, u1, u2 self.x # 数值差分求 dUoc/dSOC eps 1e-4 dUoc (self.m.ocv_func(soc eps) - self.m.ocv_func(soc - eps)) / (2 * eps) H np.array([[dUoc, -1.0, -1.0]]) u_pred self.m.ocv_func(soc) - I * self.m.R0 - u1 - u2 y u_meas - u_pred S H self.P H.T self.R K self.P H.T np.linalg.inv(S) self.x self.x (K y).flatten() # SOC 限幅防止越界 self.x[0] np.clip(self.x[0], 0.0, 1.0) self.P (np.eye(3) - K H) self.P return self.x[0]Q 和 R 的整定是玄学也是血泪经验。R 取电压传感器的实际噪声方差一般 1e-4 到 1e-3 之间。Q 的 SOC 通道取小一点1e-8 量级表示对安时积分比较信任极化通道可以大一些因为极化参数辨识误差大。如果发现 SOC 收敛慢先调大 Q 的 SOC 通道如果 SOC 抖动厉害调大 R。4.3 用动态工况数据跑闭环验证拿一段 DST 或 US06 工况的电流电压数据把 EKF 跑起来对比估算 SOC 和参考 SOC。def run_simulation(current, voltage, dt, model, ekf): soc_est [] for k in range(len(current)): ekf.predict(current[k], dt) soc ekf.update(current[k], voltage[k]) soc_est.append(soc) return np.array(soc_est) # 假设 current、voltage 已按放电为正对齐dt 为采样周期 soc_est run_simulation(current, voltage, 1.0, model, ekf) rmse np.sqrt(np.mean((soc_est - soc_ref) ** 2)) print(fSOC RMSE: {rmse*100:.2f}%)验证时重点看三个指标RMSE、最大误差、末端误差。RMSE 反映整体最大误差反映极端工况末端误差反映有没有累积漂移。如果末端误差比最大误差还大说明模型有系统性偏差回去查 OCV 曲线和 R0 的 SOC 分段。注意参考 SOC 如果是安时积分算的只能用来对比趋势不能当绝对真值。有条件的话用充放电机的 SOC 记录做参考。5. 避坑与排查仿真跑不通时先看这几条5.1 SOC 一路发散到负值现象是仿真跑几十秒后 SOC 估算值冲到 0 以下或者超过 1。原因通常是电流符号约定和模型不一致放电为正的模型喂了充电为正的数据安时积分方向反了。解决方法是统一在数据入口做符号翻转并在 EKF 里加 clip 限幅作为兜底。另外检查 Qn 单位如果容量给的是 mAh 而公式里按 Ah 算SOC 变化速度会差 1000 倍。5.2 端电压拟合很好但 SOC 误差很大这是最迷惑人的情况。端电压拟合好说明模型输出跟得上测量但 SOC 可能被 OCV 曲线的平台区“吞掉”了。在电压平台区dUoc/dSOC 很小H 矩阵对 SOC 的敏感度低EKF 几乎不修正 SOC。解决办法是在平台区适当调大 Q 的 SOC 通道或者引入安时积分的强约束别让滤波把 SOC 拉偏。另一个原因是 OCV 曲线在平台区取点太稀加密到 2% 一个点会有改善。5.3 协方差矩阵不收敛或报奇异现象是 np.linalg.inv(S) 报奇异矩阵或者 P 矩阵出现负特征值。原因通常是 R 给得太小或者 H 矩阵在某一步全零。检查 dUoc 是否算出来是 0如果 OCV 表在某个 SOC 区间是平的数值差分就是 0。解决方法是给 dUoc 设一个下限比如 1e-3保证 H 矩阵始终有非零元素。另外 P 更新用 Joseph 形式更稳但标准形式在大多数情况下够用。5.4 参数辨识出来的 C1 或 C2 是负数双指数拟合没有加边界约束时tau 可能拟合出负值反解 C 就成负数。解决方法是给 curve_fit 加 boundstau 下限设 0.1上限设几千。如果还是负检查弛豫数据是不是包含了脉冲期间的电压那段数据不符合弛豫模型要去掉。5.5 仿真步长和实际采样周期不一致有人拿 0.1s 采的数据仿真时为了“平滑”重采样到 1s结果 SOC 误差变大。原因是重采样改变了电流的积分量安时积分对电流的积分是线性的但极化电压的指数响应不是。正确做法是仿真步长严格等于采样周期如果数据有缺失用插值补点而不是降采样。6. 进阶技巧让 EKF 在嵌入式上跑得更稳的两个习惯第一个习惯是给 EKF 加一个“SOC 开路电压回正”的慢修正。车辆静置时电流为 0端电压约等于 OCV这时候可以用 OCV 反查 SOC 去缓慢修正 EKF 的 SOC 估计时间常数给几百秒避免停车后 SOC 慢慢漂。实现上就是在电流绝对值小于阈值时把 OCV 反查的 SOC 作为一个虚拟观测用很小的增益更新。第二个习惯是参数在线辨识和 EKF 解耦。不要试图用 EKF 同时估 SOC 和所有 RC 参数状态维数一高协方差矩阵条件数就爆炸。常见做法是用递推最小二乘RLS单独在线辨识 R0 和极化参数辨识结果每隔一段时间喂给 EKF。RLS 的遗忘因子取 0.98 到 0.995 之间太小会导致参数抖动太大则跟踪不上。class RLS: def __init__(self, n, lam0.99, delta1e4): self.theta np.zeros(n) self.P np.eye(n) * delta self.lam lam def update(self, phi, y): # phi 是回归向量y 是观测值 K self.P phi / (self.lam phi self.P phi) self.theta self.theta K * (y - phi self.theta) self.P (self.P - np.outer(K, phi self.P)) / self.lam return self.theta验证方法上我习惯做两组对比一组用固定参数跑 EKF一组用 RLS 在线辨识参数跑 EKF看 SOC 的 RMSE 和最大误差有没有改善。如果改善不明显说明参数时变不是主要误差源问题在 OCV 曲线或者电流采样。另外做蒙特卡洛测试给初始 SOC 加 ±20% 的随机偏差跑 100 次看收敛时间分布收敛时间超过 200 秒的工况要单独分析。最后说个我自己的教训早期做这套仿真时我花了两周调 EKF 的 Q 和 R结果发现真正的问题是 OCV 曲线测的时候静置时间不够曲线整体偏低 20mV导致 SOC 在低端一直偏大。后来养成习惯任何 SOC 估算的异常先回去查 OCV 曲线和电流符号这两个对了EKF 的参数整定才有意义。希望帮到你。本文还有配套的精品资源点击获取