估计锂电池的荷电状态(SOC)主要完成了扩展卡尔曼滤波(EKF)的 实验、参数辨识和仿真完成无迹卡尔曼滤波(UKF)

电池圈聊SOC，总绕不开卡尔曼滤波这棵“老树”，但直接上标准KF肯定不行——电池这玩意儿电阻电压曲线弯得跟过山车似的，完全非线性嘛！所以今天唠唠实打实用过的扩展卡尔曼滤波（EKF）做的实验，还有顺手补的仿真里的无迹卡尔曼滤波（UKF），加参数辨识凑成一套接地气的入门SOC估计包。

---

先搭个小舞台：二阶RC等效电路模型

不管滤波多花里胡哨，先得给电池找个“替身”模型，太复杂的电化学模型算不动，太简单的零阶内阻模型又飘得像风筝，二阶RC刚好——既能算欧姆压降（R0），又能分两块（R1C1慢极化、R2C2快极化）算暂态的内阻延迟，够准够快。

先甩个Matlab里搭这个模型状态方程的片段，顺便揉点白话解释公式推导省得翻书：

% 二阶RC模型状态向量：SOC(核心要算的) + 慢极化电压U1 + 快极化电压U2
x = [SOC; U1; U2];
T = 1; % 采样周期设1秒，实验室的电池测试仪一般都能稳在10Hz-1Hz
C1 = 1e4; % 先随便蒙俩辨识前的初始值当对比
C2 = 1e3;
R0 = 0.02;
R1 = 0.01;
R2 = 0.005;
Qn = 2.5; % 标称容量，实验室是个2.5Ah的软包三元锂

% 状态方程简化离散版（EKF里不用太复杂的积分，欧拉就行，入门首选）
% 一阶RC的极化电压更新公式是 U1(k+1) = exp(-T/(R1C1))*U1(k) + R1*(1-exp(-T/(R1C1)))*I(k)
% SOC就是简单的安时积分：SOC(k+1) = SOC(k) - (I(k)*T)/(Qn*3600)
% 注意电流正负！放电I取负，充电取正（这个细节要是搞反，SOC直接反着跑）
A = [1, 0, 0;
     0, exp(-T/(R1*C1)), 0;
     0, 0, exp(-T/(R2*C2))];
B = [ -T/(Qn*3600);
      R1*(1-exp(-T/(R1*C1)));
      R2*(1-exp(-T/(R2*C2)))];
f = @(x, I) A*x + B*I; % 非线性？这里状态方程其实是线性的！非线性全在观测方程里

对，核心就在观测方程的非线性——开路电压OCV和SOC不是直线！三元锂的在10%-90%SOC区间虽然还算“缓坡”，但0%-10%和90%-100%那简直是悬崖峭壁，必须用非线性函数拟合，实验室用的是七阶多项式（别问为啥用七阶，试了五阶差点意思，九阶有点过拟合，七阶刚好卡在软包三元锂的OCV-SOC拐点附近）。

观测方程的Matlab代码和函数：

% 三元锂10%-90%拟合的七阶系数（先随便找个公开数据的系数，后面会说自己怎么用参数辨识补实测的）
OCV_coeff = [0.0002, -0.0035, 0.021, -0.062, 0.11, 0.32, 2.75, 3.0];
h = @(x, I) polyval(OCV_coeff, x(1)) + x(2) + x(3) + I*R0; % 实测电压 = OCV(SOC) + 慢极化 + 快极化 + 欧姆压降

---

接下来是EKF的实验室实测小翻车+小修复

参数辨识之前先别急着跑真实数据？不对，我是先拿实验室测试仪的循环充放电数据（HPPC+恒流恒压组合循环）试了试初始值瞎蒙的EKF，先看个初始蒙值的SOC估计图（脑补图就行，论文才放图，这里用文字说：初始值蒙的SOC0=0.9，真实值是0.87，跑了2分钟快放电（3C，三元锂一般撑得住），EKF居然直接飘到1.05了！因为观测方程的R0、C1、R1这些全不对嘛！欧姆压降算小了，慢极化电压上升算慢了，导致h函数算出来的电压比真实值高好多，EKF为了凑电压，拼命把SOC往上调——刚好悬崖峭壁的90%-100%区间OCV对SOC不敏感，往上调0.05，OCV只涨了0.02V，刚好补上瞎蒙参数的误差，结果就飘了）。

估计锂电池的荷电状态(SOC)主要完成了扩展卡尔曼滤波(EKF)的 实验、参数辨识和仿真完成无迹卡尔曼滤波(UKF)

小修复第一步：参数辨识

实验室常用的离线参数辨识方法是HPPC脉冲测试+最小二乘法（LS），脉冲测试就是给充满电的电池恒流放10s，停40s，恒流充10s，停40s，循环个10次，覆盖90%-10%SOC区间。

先把辨识R0、R1、C1、R2、C2的LS逻辑写成Matlab片段（离线辨识不用太讲究实时性，Matlab的polyfit和lsqcurvefit都能用，这里用lsqcurvefit，更灵活）：

% 假设已经把某一个SOC点（比如0.8）的HPPC停40s的电压数据存到了Y_data，电流存到了I_data，时间存到了t_data
% 要拟合的函数就是二阶RC的开路电压响应（去掉欧姆压降的实测电压，初始开路电压是停40s的起始电压V0）
% V(t) = V0 + R1*I0*(1-exp(-t/(R1C1))) + R2*I0*(1-exp(-t/(R2C2))) 这里I0是放电的负电流
V0 = Y_data(1);
I0 = I_data(1); % 整个脉冲响应期间电流都是0或者I0，这里取停40s前的放电电流
t_data = t_data - t_data(1); % 时间归一化

% 定义待拟合的参数向量：theta = [R1, C1, R2, C2]，R0可以单独算（停40s前放电瞬间的电压降除以I0的绝对值）
R0_est = abs((V0 - Y_data_before_pulse(end))/I0);

% 定义拟合函数句柄
fit_fun = @(theta, t) V0 + theta(1)*I0*(1 - exp(-t/(theta(1)*theta(2)))) + theta(3)*I0*(1 - exp(-t/(theta(3)*theta(4))));

% 初始猜测值（根据电池标称容量和经验值猜，经验值：R1是R0的1-2倍，C1是1e4-1e5F，R2是R0的0.2-0.5倍，C2是1e2-1e3F）
theta0 = [0.03, 5e4, 0.01, 5e2];

% 调用lsqcurvefit
options = optimoptions('lsqcurvefit','Display','iter');
theta_est = lsqcurvefit(fit_fun, theta0, t_data, Y_data, [], [], options);

% 提取辨识好的参数
R1_est = theta_est(1);
C1_est = theta_est(2);
R2_est = theta_est(3);
C2_est = theta_est(4);

辨识好每个SOC点的参数后，我取了中间值当全局参数（入门SOC估计不用做自适应参数，自适应的下次有空再唠），再重新拟合了自己测的OCV-SOC曲线的七阶系数，这次再跑EKF，飘的问题解决了！

---

顺便补个仿真里的UKF对比

UKF的核心就是不用算雅可比矩阵（EKF最烦人的一步，要是模型复杂点雅可比矩阵算错一个符号，估计直接崩），用“sigma点”来近似状态向量的均值和协方差，对于软包三元锂这种OCV-SOC不算特别非线性的模型，其实和EKF的精度差不了多少，但UKF的代码框架更通用——下次要是换成磷酸铁锂（OCV-SOC在30%-70%区间平得像飞机场，EKF的雅可比矩阵在这个区间可能会有奇异值问题），直接改模型函数就行，不用重新推雅可比。

甩个UKF更新sigma点和观测预测的核心片段（Matlab的unscentedTransform是个好东西，不用自己写sigma点的计算逻辑）：

% 假设已经定义好了UKF的参数：alpha（控制sigma点离均值的距离，一般取1e-3）、beta（状态向量的分布假设，高斯分布取2）、kappa（一般取0）
n = length(x); % 状态向量维度3
alpha = 1e-3;
beta = 2;
kappa = 0;
lambda = alpha^2*(n + kappa) - n;
c = n + lambda;
Wm = [lambda/c, 0.5/c + zeros(1, 2*n)]; % 均值权重
Wc = Wm;
Wc(1) = Wc(1) + (1 - alpha^2 + beta); % 协方差权重

% 第一步：生成sigma点
Xsig = unscentedTransform(x, P, lambda); % Matlab自带的生成sigma点函数，P是状态协方差矩阵

% 第二步：状态预测（用每个sigma点代入状态方程f）
Xsig_pred = zeros(n, 2*n+1);
for i = 1:2*n+1
    Xsig_pred(:,i) = f(Xsig(:,i), I_k);
end

% 第三步：计算状态预测均值和协方差
x_pred = sum(Wm.*Xsig_pred, 2);
P_pred = zeros(n,n);
for i = 1:2*n+1
    P_pred = P_pred + Wc(i)*(Xsig_pred(:,i) - x_pred)*(Xsig_pred(:,i) - x_pred)';
end
P_pred = P_pred + Q; % Q是过程噪声协方差矩阵，入门可以设成对角小矩阵，比如Q = diag([1e-6, 1e-8, 1e-9])

% 第四步：观测预测（用每个状态预测sigma点代入观测方程h）
Ysig_pred = zeros(1, 2*n+1);
for i = 1:2*n+1
    Ysig_pred(:,i) = h(Xsig_pred(:,i), I_k);
end

% 第五步：计算观测预测均值和协方差，还有互协方差
y_pred = sum(Wm.*Ysig_pred, 2);
Pyy = 0;
Pxy = 0;
for i = 1:2*n+1
    Pyy = Pyy + Wc(i)*(Ysig_pred(:,i) - y_pred)*(Ysig_pred(:,i) - y_pred)';
    Pxy = Pxy + Wc(i)*(Xsig_pred(:,i) - x_pred)*(Ysig_pred(:,i) - y_pred)';
end
Pyy = Pyy + R; % R是观测噪声协方差矩阵，实验室电池测试仪的电压误差一般是1mV，所以R = 1e-6

% 第六步：卡尔曼增益，更新状态和协方差
K = Pxy / Pyy;
x = x_pred + K*(y_k - y_pred);
P = P_pred - K*Pyy*K';

仿真结果的话（还是脑补文字：拿UDDS城市道路工况的电流数据跑，EKF的SOC最大误差是1.2%，UKF的是0.9%，确实差不了多少，但UKF的代码比我手动推二阶RC模型雅可比矩阵的EKF代码短了30行左右！而且下次要是加个温度补偿的非线性项，EKF要重新推三个雅可比矩阵的偏导，UKF只要改f和h函数就行，爽歪歪）。

---

最后唠两句小总结

• 入门SOC估计，二阶RC+全局参数辨识+离线EKF/UKF完全够用，实车测试的话再考虑加自适应参数和温度补偿。
• 参数辨识一定要用自己实验室/自己手头电池的数据！不同批次的软包三元锂参数差得还挺多的，公开数据只能当初始猜测值。
• 如果模型不算特别复杂，雅可比矩阵好推，用EKF就行，实时性比UKF稍微好一点点（因为UKF要算2n+1个sigma点）；如果模型复杂，雅可比矩阵不好推，或者是磷酸铁锂这种OCV-SOC平得像飞机场的，直接上UKF。