基于卡尔曼滤波器 (Kalman Filter) 的 剩余电池电量 (SOC) 估计算法的仿真设计

背景

在电源管理系统 (BMS) 中通常使用电流积分法来计算电池当前电量,但是电流采集的累计误差、以电流积分为基础的系统未建模噪声会导致 SOC 计算产生明显偏差，比如电量计算到达 0% 或者 100% 时，电池电压却未到达预定值，或超过了预定值，虽然通过添加一些软件策略可以应对以上问题，但是效果始终不那么完美。

设计核心思想：电流、电压测量有累计噪声，电流积分模型有过程噪声，电压同电池电量在充电电压平台期变化微小，几乎被噪声淹没。单靠其中一个很难准确估计电池电量，所以需要将两套有噪声的估计通过卡尔曼滤波器进行数据融合，减小总的估计误差的同时，使电压和电池电量的对应变化更加平滑。

设计过程

建立系统预测方程和观测方程

以电流积分方法建立电池充电过程的系统状态差分方程为：

\[x_{[k]} = A x_{[k-1]} + w_{[k-1]} \]

状态变量\(x_{[k]} = \begin{bmatrix} Q_{[k]} \\ i_{[k]} \end{bmatrix}\), \(Q_{[k]}\)表示电量，\(i_{[k]}\)表示电流，\(k\) 表示采样序列。状态转移矩阵 \(A = \begin{bmatrix} 1 & T_s \\ 0 & 1 \end{bmatrix}\)，\(T_s\) 表示一个采样周期。\(w_{[k-1]}\) 表示过程噪声，过程噪声服从零均值的正态分布 \(N(0, C_w)\)，过程噪声协方差 \(C_w = \begin{bmatrix}\sigma^2_Q & 0 \\ 0 & \sigma^2_i \end{bmatrix}\) 。这也是电池电量的预测方程。需要说明的是，预测方程中需要使用真实电流，这在实际中无法直接得到，因此这里使用电流的后验估计值。

观测方程需要电压和电流的测量值，而电池电压和电池剩余电量可以建立一个函数关系：

\[v_{[k]} = OCV(Q_{[k]}) \]

在工程上这通常需要实测电池电压和电池剩余电量的对应值表，通过查表法和插值法可以得到电池电量的估计值。

实际的电压电流测量值要考虑测量噪声以及电压降：

\[v_{mea,[k]} = v_{[k]} + q_{v,[k]} + i_{mea,[k]}r \\[10pt] i_{mea,[k]} = i_{[k]} + q_{i,[k]} \]

\(q_{v,[k]}\) 表示电压噪声，\(q_{i,[k]}\) 表示电流噪声，\(i_{[k]} r\) 表示电流流进电池时，由电池内阻、线路板产生的电压降，\(q_{v,[k]}\) 服从零均值的正态分布 \(N(0, \sigma^2_v)\), \(q_{i,[k]}\) 服从零均值的正态分布 \(N(0, \sigma^2_i)\)

假设包含电池电量和电流的观测向量为：

\[z_{[k]} = h(x_{[k]}) + q_{[k]}\\[10pt] z_{[k]} = \begin{bmatrix} V_{mea,[k]} \\ i_{mea,[k]} \end{bmatrix} \]

其中，\(h(x)\) 是观测方程的函数，它是一个非线性函数，\(h(x_k) = \begin{bmatrix} OCV(Q_{k})-i_{k}r \\ i_{k} \end{bmatrix}\)。观测噪声 \(q_{[k]}\) 服从零均值的正态分布 \(N(0, C_q)\)，\(C_q\) 表示观测噪声协方差矩阵，\(C_q = \begin{bmatrix} \sigma^2_v & 0 \\ 0 & \sigma^2_i \end{bmatrix}\)。

通过对 \(h(x)\) 求导得到观测矩阵 \(H_{[k]}\)：

\[H_{[k]} = \dfrac{\partial h}{\partial x} \bigg|\hat{x}^-_{[k]} = \begin{bmatrix} \dfrac{\partial OCV}{\partial Q} \bigg|\hat{Q}^-_{[k]} & -r \\ 0 & 1 \end{bmatrix} \]

\(H_{[k]}\) 是观测矩阵，值得注意的是，由于电压和电池电量的关系是非线性的，因此这个矩阵是时变的，在每一轮估计时都应该根据静态工作点 \(\hat{Q}^-_{[k]}\) 去更新。

计算先验估计值以及先验估计误差方差值

先验估计是用上一拍的后验估计，通过系统方程，在不考虑过程噪声的情况下得到的估计：

\[\hat{x}^-_{[k]} = A \hat{x}_{[k-1]} \]

由先验估计得到先验估计误差，它是与真实值之间的偏差：

\[e^-_{[k]} = x_{[k]} - \hat{x}^-_{[k]} \\[10pt]= (Ax_{[k-1]} + w_{[k-1]}) - (A \hat{x}_{[k-1]}) \\[10pt]= A e_{[k-1]} + w_{[k-1]} \]

先验估计误差方差：

\[P^-_{[k]} = E\{(e^-_{[k]})(e^-_{[k]})^T\} \\[10pt]= E\{(A e_{[k-1]} + w_{[k-1]})(A e_{[k-1]} + w_{[k-1]})^T\} \\[10pt]= A P_{[k-1]} A^T + C_w \]

这意味着可以通过上一拍的后验估计误差方差 \(P_{[k-1]}\) 来计算当前拍的先验估计误差方差 \(P^-_{[k]}\)。

计算最优卡尔曼增益值

先列写出后验估计误差方差的表达式：

\[P_{[k]} = E\{(e_{[k]})(e_{[k]})^T\}\\[10pt] = K_{[k]} H_{[k]} P^-_{[k]} H_{[k]}^T K_{[k]}^T+ K_{[k]} C_q K_{[k]}^T \]

它是一个关于卡尔曼增益 \(K_{[k]}\) 的函数，目标是使它最小化。

由\(\dfrac{\partial Tr(P_{[k]})}{\partial K_{[k]}} = 0\)，得到最优卡尔曼增益：

\[K_{[k]} = P^-_{[k]}H_{[k]}^T({H_{[k]} P^-_{[k]}H_{[k]}^T + C_q})^{-1} \]

计算后验估计值以及后验估计误差方差值

电量值的后验估计，同时也是卡尔曼滤波器的输出值：

\[\hat{x}_{[k]} = \hat{x}^-_{[k]} + K_{[k]}\Big(z_{[k]} - h(\hat{x}^-_{[k]})\Big) \]

接着计算后验估计误差：

\[e_{[k]} = x_{[k]} - \hat{x}_{[k]}\\[10pt] = x_{[k]} - \bigg(\hat{x}^-_{[k]} + K_{[k]}\Big(z_{[k]} - h(\hat{x}^-_{[k]})\Big)\bigg)\\[10pt] = e^-_{[k]} - K_{[k]}\Big(h(x_{[k]}) + q_{[k]} - h(\hat{x}^-_{[k]})\Big)\\[10pt] \]

对 \(h(x)\) 在 \(\hat{x}^-_{[k]}\) 处进行泰勒展开，忽略二阶项，取一阶近似：

\[h(x_{[k]}) \approx h(\hat{x}^-_{[k]}) + H_{[k]} e^-_{[k]} \]

带入上式得：

\[e_{[k]} = (I_2 - K_{[k]}H_{[k]}) e^-_{[k]} - K_{[k]} q_{[k]} \]

计算后验估计误差方差 \(P_{[k]}\) ：

\[P_{[k]} = E\{(e_{[k]})(e_{[k]})^T\}\\[10pt] = E\{\Big((I_2 - K_{[k]}H_{[k]}) e^-_{[k]} - K_{[k]} q_{[k]}\Big)\Big((I_2 - K_{[k]}H_{[k]}) e^-_{[k]} - K_{[k]} q_{[k]}\Big)^T\} \\[10pt] = (I_2 - K_{[k]}H_{[k]})P^-_{[k]}(I_2 - K_{[k]}H_{[k]})^T + K_{[k]} C_q K_{[k]}^T \\[10pt] = (I_2 - K_{[k]}H_{[k]})P^-_{[k]} - (I_2 - K_{[k]}H_{[k]})P^-_{[k]}(K_{[k]}H_{[k]})^T + K_{[k]} C_q K_{[k]}^T \\[10pt] \]

因为最优卡尔曼增益满足 \(K_{[k]}{H_{[k]} P^-_{[k]}H_{[k]}^T + K_{[k]}C_q} = P^-_{[k]}H_{[k]}^T\)，所以上式右边两项为零，可以简化为：

\[P_{[k]} = (I_2 - K_{[k]}H_{[k]})P^-_{[k]} \]

以上构成了卡尔曼滤波器单个估计周期的计算过程，通过重复递归计算，可以得到电池电量的实时估计。
下面就需要对以上设计进行仿真验证了。

仿真验证

仿真采用以下参数：

参数	取值	含义
\(T_s\)	1 s	采样周期
\(C_N\)	10 Ah	电池容量
\(R_0\)	0.005 Ω	内阻
I	0.5A	电流
\(\sigma_Q\)	\(10^{-6}\)	SOC 过程噪声标准差
\(\sigma_i\)	\(10^{-6}\)	电流过程噪声标准差
\(\sigma_V\)	3 mV	电压测量噪声标准差
\(\sigma_i^{meas}\)	5 mA	电流测量噪声标准差
\(P_0\)	diag(0.01, 0.001)	初始协方差
N	10000	迭代次数

OCV-SOC 模型采用三次多项式代替查表：

\[OCV(SOC) = 2.5 + 1.2\,SOC - 0.5\,SOC^2 + 0.3\,SOC^3 \]

仿真结果：

可以看到 SOC 估计对比曲线栏，黑色实线是用真实电流作出的真实 SOC，红色虚线是用 EKF 估计的 SOC，蓝色点是用电流测量值积分估计的 SOC。初始时刻，两种估计同时给予 5% 的误差，电流积分的估计值全程保持误差没变，EKF 估计值则会随着时间的增加而减少误差，逐渐逼近真实值。