基于等效电路与扩展卡尔曼滤波的SOC估算方法

等效电路模型

简介

参考：锂电池的二阶RC网络等效电路模型 - 車馬炮的文章 - 知乎 https://zhuanlan.zhihu.com/p/407572989

二阶RC网络等效电路模型

• 锂离子动力电池正负极之间存在电位差，外在表现为电池的电动势能，用电容器模拟；电池电动势能与电池SOC之间存在关系，在模型中用EMF=f(SOC)表达式描述；

• 锂离子动力电池通过锂离子在正负极之间的嵌入、脱嵌和移动来实现电流的输出，在这个电化学反应和离子运动的过程中，必然会产生一定的阻力，因此在建模时，我们使用电池的内阻来模拟。研究表明，在电池工作时，由于形成的机理不同，电池的内阻又可以进一步划分为欧姆内阻和极化内阻：

1. 欧姆内阻是由于离子在运动过程中受到的阻力而产生的；
2. 极化内阻是由于材料本身的一些性质而产生的，可细分为电化学极化内阻和浓差极化内阻 。

• 电化学极化也称活化极化，是由于正负极活性物质发生的电化学反应速率小于电子运动速率引起的极化，响应时间微秒级；
• 浓差极化是由于反应物消耗引起电极表面得不到及时补充（或是某种产物在电极表面积累，不能及时疏散），例如锂离子在电池正极的积累，导致电极电势偏离通电前按总体浓度计算的平均值，响应时间秒级。

在模型中，通过电阻来描述欧姆内阻，通过两个阻容网络（RC 网络）来描述极化内阻。通过和来描述电池工作时的负载电流和电池的端电压。从图中可以看出改进模型通过二极管的单向导电性对欧姆内阻区别充电与放电设置。放电时对应放电欧姆内阻，充电时对应充电欧姆内阻。考虑到电池的开路电压应该包括滞回电压和电动势两部分，所以增加了等效电压源部分：其中 EMF 是电池平衡电势，是滞回电压（历史电流方向会影响的正负号）。

基于等效电路与扩展卡尔曼滤波的SOC估算方法

1. 数据需求

HPPC测试数据：三元锂与磷酸铁锂

基于HPPC测试数据对等效电路模型中的SOC-OCV，内阻，极化电阻和极化电容进行识别，在测试前先对电芯进行恒流恒压满充，然后根据如下步骤进行测试:

步骤	电池状态	时间（s）	电流倍率（C）	SOC（%）
1	静置	7200	0	100
2	放电	360	0.5	95
3	静置	1200	0	95
4	充电	10	1	95
5	静置	40	0	95
6	放电	10	0.75	95
7	静置	180	0	95
8	放电	360	0.5	90
9	静置	1200	0	90
…	…	…	…	…

2. 方法介绍

2.1 等效电路模型

二阶等效电路的动态方程如下：
$$\begin{array}{rl}{U}_L& =U_{OCV}-I{R}_0-U_{p1}-U_{p2}\\ I& =\frac{U_{p1}}{R_{p1}}+C_{p1}\frac{d{U}_{p1}}{dt}=\frac{U_{p2}}{R_{p2}}+C_{p2}\frac{d{U}_{p2}}{dt}\end{array}$$
求解微分方程$I=\frac{U_p}{R_p}+C_p\frac{d{U}_p}{dt}$有：
$$\begin{array}{rl}{U}_p& =I{R}_p\left(1-e^{-\frac{t}{R_p{C}_p}}\right)+U_p(0)e^{-\frac{t}{R_p{C}_p}}\\ U_L& =U_{OCV}-I{R}_0-\left(I{R}_{p_1}\left(1-e^{-\frac{t}{R_{p_1}{C}_{p_1}}}\right)+U_{p_1}(0)e^{-\frac{t}{R_{p_1}{C}_{p_1}}}\right)-\left(I{R}_{p_2}\left(1-e^{-\frac{t}{R_{p_2}{C}_{p_2}}}\right)+U_{p_2}(0)e^{-\frac{t}{R_{p_2}{C}_{p_2}}}\right)\end{array}$$

2.2 参数识别

2.2.1 SOC-OCV

选择100%，95%，…，0%的测试节点中长时间静置的末端电压作为对应的OCV.

2.2.2 内阻

$$R_0=\frac{\Delta U_{AB}+\Delta U_{DC}}{2I}$$

2.2.3 极化电阻&极化电容

当BC段放电开始时，为零状态响应，有$t=0,U_p(0)=0$：
$$U_L=U_{OCV}-I{R}_0-I{R}_{p_1}\left(1-e^{-\frac{t}{R_{p_1}{C}_{p_1}}}\right)-I{R}_{p_2}\left(1-e^{-\frac{t}{R_{p_2}{C}_{p_2}}}\right)$$
当DE段搁置时，为零输入响应，有$I=0$：
$$U(t)=U_{OCV}-U_{p1}(0)e^{-\frac{t}{R_{p_1}{C}_{p_1}}}-U_{p2}(0)e^{-\frac{t}{R_{p_2}{C}_{p_2}}}$$
用指数拟合的方法可以求出$U_{p1}(0),U_{p2}(0),-\frac{1}{R_{p_1}{C}_{p_1}},-\frac{1}{R_{p_2}{C}_{p_2}}$。因为DE段开始为BC段结尾，因此：
$$\begin{array}{rl}{U}_{p1}(0)& =I{R}_{p_1}\left(1-e^{-\frac{T_C}{R_{p_1}{C}_{p_1}}}\right)\\ U_{p2}(0)& =I{R}_{p_2}\left(1-e^{-\frac{T_C}{R_{p_2}{C}_{p_2}}}\right)\end{array}$$
从而求出具体$R_p，C_p$。对与更高阶的等效电路模型，求得的方法一致，只是参数更多。

2.3 基于扩展卡尔曼滤波的SOC估算

根据等效电路的参数识别有下表

SOC(%)	$U_{OCV}$	$R_{p1}$	$R_{p2}$	$C_{p1}$	$C_{p2}$	$R_0$
100	$u_0$	$r_{11}$	$r_{21}$	$c_{11}$	$c_{21}$	$r_{01}$
95	$u_1$	$r_{12}$	$r_{22}$	$c_{12}$	$c_{22}$	$r_{02}$
…	…	…	…	…	…	…

状态变量设置为$X_k={\left(SOC(k),U_{p1}(k),U_{p2}(k)\right)}^T$，电芯最大可用容量C，时间间隔为t，则状态方程为
$$X_k=A{X}_{k-1}+B{I}_k+{\omega }_k$$
其中，${\omega }_k\sim N(0,Q_k)$，$A=\left(\begin{array}{ccc}1& 0& 0\\ 0& e^{-\frac{t}{C_{p1}{R}_{p1}}}& 0\\ 0& 0& e^{-\frac{t}{C_{p1}{R}_{p1}}}\end{array}\right)$，$B=\left(\begin{array}{c}-\frac{t\times 100}{C}\\ R_{p1}(1-e^{-\frac{t}{C_{p1}{R}_{p1}}})\\ R_{p2}(1-e^{-\frac{t}{C_{p2}{R}_{p2}}})\end{array}\right)$，矩阵中参数根据$SOC(k)$进行匹配.

因为观测方程为
$$\begin{array}{rl}{U}_k& =U_{OCV}\left(SOC(k)\right)-U_{p1}(k)-U_{p2}(k)-I_k{R}_0+v_k\end{array}$$
其中，$v_k\sim N(0,R_k)$.

把观测方程进行泰勒展开：
$$\begin{array}{rl}{U}_k& =g(X_k)=U_{OCV}\left(SOC(k)\right)-U_{p1}(k)-U_{p2}(k)-I_k{R}_0+v_k\\ U_k& \approx g(X_{k-1})+g^{\prime }(X_{k-1})(X(k)-X(k-1))\\ & \\ H& =g^{\prime }(X_{k-1})=\left[\frac{\alpha g}{\alpha SOC(k-1)}\frac{\alpha g}{\alpha U_{p1}(k-1)}\frac{\alpha g}{\alpha U_{p1}(k-1)}\right]\\ & =\left[U_{OCV}^{\prime }(SOC(k-1))-1-1\right]\\ & \\ g(X_{k-1})& =U_{OCV}\left(SOC(k-1)\right)-U_{p1}(k-1)-U_{p2}(k-1)-I{R}_0\\ & \\ & \\ & \\ U_k& \approx U_{OCV}\left(SOC(k-1)\right)-U_{p1}(k-1)-U_{p2}(k-1)-I{R}_0\\ & +\left[U_{OCV}^{\prime }(SOC(k-1))-1-1\right]\left[\begin{array}{c}SOC(k)\\ U_{p1}(k)\\ U_{p2}(k)\end{array}\right]\\ & -\left[U_{OCV}^{\prime }(SOC(k-1))-1-1\right]\left[\begin{array}{c}SOC(k-1)\\ U_{p1}(k-1)\\ U_{p2}(k-1)\end{array}\right]\\ & =H{X}_k-I{R}_0+D_{k-1}\end{array}$$

上述过程可简化为以下流程，由于
$$U_{OCV}\left(SOC(k)\right)\approx U_{OCV}\left(SOC(k-1)\right)+U_{OCV}^{\prime }(SOC(k-1))(SOC(k)-SOC(k-1))$$
因此观测方程可写成
$$\begin{array}{rl}{U}_k& \approx H{X}_k-I_k{R}_0+D_{k-1}+v_k\end{array}$$
其中
$$H=\left(\begin{array}{c}{U}_{OCV}^{\prime }(SOC(k-1)),-1,-1\end{array}\right)$$

$$D_{k-1}=U_{OCV}\left(SOC(k-1)\right)-U_{OCV}^{\prime }(SOC(k-1))SOC(k-1)$$

根据卡尔曼滤波算法，算法步骤如下:

1.初始化：

初始状态为$X_0={\left(SOC(0),U_{p1}(0),U_{p2}(0)\right)}^T$，其协方差矩阵为$P_0$，$Q_k$设置为固定值$Q$，$R_k$设置为固定值$R$.

2.预测阶段：
$$X_k=A{X}_{k-1}+B{I}_k$$

$$P_k=A{P}_{k-1}{A}^T+Q_k$$

3.修正阶段：

卡尔曼增益为
$$K_k=\frac{P_k{H}^T}{H{P}_k{H}^T+R_k}$$
修正状态为
$$X_k^{\prime }=X_k+K_k(z_k-U_k)$$
修正状态的协方差矩阵为
$$P_k^{\prime }=\left(I-K_{k}H\right)P_k$$

2.4 基于无迹卡尔曼滤波的SOC估算

算法步骤：

1. 初始化

初始状态为$X_0={\left(SOC(0),U_{p1}(0),U_{p2}(0)\right)}^T$，其协方差矩阵为$P_0$，$Q_k$设置为固定值$Q$，$R_k$设置为固定值$R$.

2. 预测阶段

状态变量设置为$X_k={\left(SOC(k),U_{p1}(k),U_{p2}(k)\right)}^T$，电芯最大可用容量C，时间间隔为t，则状态方程为
$$X_k=A{X}_{k-1}+B{I}_k+{\omega }_k$$
其中，${\omega }_k\sim N(0,Q_k)$，$A=\left(\begin{array}{ccc}1& 0& 0\\ 0& e^{-\frac{t}{C_{p1}{R}_{p1}}}& 0\\ 0& 0& e^{-\frac{t}{C_{p1}{R}_{p1}}}\end{array}\right)$，$B=\left(\begin{array}{c}-\frac{t\times 100}{C}\\ R_{p1}(1-e^{-\frac{t}{C_{p1}{R}_{p1}}})\\ R_{p2}(1-e^{-\frac{t}{C_{p2}{R}_{p2}}})\end{array}\right)$，矩阵中参数根据$SOC(k)$进行匹配.

• 构造sigma点集，X有3个变量，n=3。

$$\begin{array}{rl}{X}_k^{(0)}& =X_{k-1}\\ & \\ X_k^{(1)}& =X_{k-1}+{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第1列}\\ X_k^{(2)}& =X_{k-1}+{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第2列}\\ X_k^{(3)}& =X_{k-1}+{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第3列}\\ & \end{array}$$

$$\begin{array}{rl}{X}_k^{(4)}& =X_{k-1}-{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第1列}\\ X_k^{(5)}& =X_{k-1}-{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第2列}\\ X_k^{(6)}& =X_{k-1}-{\left(\sqrt{(n+\lambda )P_{k-1}}\right)}_{第3列}\\ & \end{array}$$

• 各个sigma点代入状态方程计算

$$\begin{array}{rl}{X_k^{(i)}}^{-}& =A{X}_k^{(i)}+B{I}_k\\ X_k^{-}& =\sum _{i=0}^{2n}{w}_m^i{X_k^{(i)}}^{-}\\ P_k^{-}& =\sum _{i=0}^{2n}{w}_c^i[{X_k^{(i)}}^{-}-X_k^{-}][{X_k^{(i)}}^{-}-X_k^{-}{]}^T+Q\end{array}$$

其中，
$$\{\begin{array}{c}{w}_m^{(0)}=\frac{\lambda }{n+\lambda }\\ w_m^{(i)}=\frac{\lambda }{2(n+\lambda )},i=1\sim 2n\\ w_c^{(0)}=\frac{\lambda }{n+\lambda }+\left(1-a^2+\beta \right)\\ w_c^{(i)}=\frac{\lambda }{2(n+\lambda )},i=1\sim 2n\end{array}$$
无迹变换引入四个参数：$\alpha ,\beta ,\kappa ,\lambda$。其中$\lambda$满足：
$$\lambda ={\alpha }^2(n+\kappa )-n$$
应为是线性变换，上述的计算结果与下面先验计算一样
$$X_k^{-}=A{X}_{k-1}+B{I}_k$$

$$P_k^{-}=A{P}_{k-1}{A}^T+Q_k$$

3. 修正阶段

观测方程为
$$\begin{array}{rl}{U}_k& =U_{OCV}\left(SOC(k)\right)-U_{p1}(k)-U_{p2}(k)-I_k{R}_0+v_k\end{array}$$
其中，$v_k\sim N(0,R_k)$.

• 构造sigma点集

$$\begin{array}{rl}{X}^{(0)}& =X_k^{-}\\ & \\ X^{(1)}& =X_k^{-}+{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第1列}\\ X^{(2)}& =X_k^{-}+{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第2列}\\ X^{(3)}& =X_k^{-}+{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第3列}\\ & \end{array}$$

$$\begin{array}{rl}{X}^{(4)}& =X_k^{-}-{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第1列}\\ X^{(5)}& =X_k^{-}-{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第2列}\\ X^{(6)}& =X_k^{-}-{\left(\sqrt{(n+\lambda )P_k^{-}}\right)}_{第3列}\\ & \end{array}$$

• 计算各个sigma点的非线性变换，sigma点集代入观测方程，得到预测的观测量：

$$\begin{array}{rl}{X}^{(i)}& ={\left(SO{C}^{(i)},U_{p1}^{(i)},U_{p2}^{(i)}\right)}^T\\ Z(k{)}^{(i)}& =U_{OCV}(SO{C}^{(i)})-U_{p1}^{(i)}-U_{p2}^{(i)}-I_k{R}_0\end{array}$$

• 通过加权求得观测量新的均值及协方差：

$$\begin{array}{rl}\bar{Z}(k)& =\sum _{i=0}^{2n}{w}_m^{i}Z(k{)}^{(i)}\\ P_{z_k{z}_k}& =\sum _{i=0}^{2n}{w}_c^i[Z(k{)}^{(i)}-\bar{Z}(k)][Z(k{)}^i-\bar{Z}(k){]}^T+R\\ P_{x_k{z}_k}& =\sum _{i=0}^{2n}{w}_c^i[{X^{(i)}}_k^{-}-X_k^{-}][Z(k{)}^i-\bar{Z}(k){]}^T\end{array}$$

其中，
$$\{\begin{array}{c}{w}_m^{(0)}=\frac{\lambda }{n+\lambda }\\ w_m^{(i)}=\frac{1}{2(n+\lambda )},i=1\sim 2n\\ w_c^{(0)}=\frac{\lambda }{n+\lambda }+\left(1-a^2+\beta \right)\\ w_c^{(i)}=\frac{1}{2(n+\lambda )},i=1\sim 2n\end{array}$$

• 更新K值

$$K_k=P_{x_k{z}_k}{P}_{z_k{z}_k}^{-1}$$

• 系统状态更新和协方差更新

$$\begin{array}{rl}\hat{X}(k)& =X_k^{-}+K_k[Z(k)-\bar{Z}(k)]\\ P(k)& =P_k^{-}-K_k{P}_{z_k{z}_k}{K}_k^T\end{array}$$

2.3 基于自适应卡尔曼滤波

参考论文：《Adaptive Adjustment of Noise Covariance in Kalman Filter for Dynamic State Estimation》

一般情况下，kalman滤波中的Q、R都是根据经验、实验或数据手册得到的，但是有些参数是无法获得的，尤其是过程噪声，就需要通过不断试凑确定参数，显然是不靠谱的。

自适应卡尔曼滤波提供一种自适应调整Q、R的方法，不需要精确的初值，就可以获得较高的滤波精度。

• 只要Q、R的比值和Qtrue、Rtrue相等，就能获得最优的性能，并不需要保证Q、R和真实值接近。