基于平滑性假设的锂电池开路电压-荷电状态映射提取

摘要

开路电压（OCV）是准确估计电动汽车锂离子电池荷电状态（SoC）的重要参数。由于OCV-SoC的映射关系随着电池老化而持续变化，因此在某一特定阶段确定的OCV-SoC函数无法适用于电池全生命周期内的SoC估计，由此需要对OCV进行定期测试及老化校准。受OCV-SoC曲线迟滞现象的影响，传统的OCV测试通常需要数天时间才能获得一个或多个完全充放电周期的数据，因此从电动车实际运行的维度上缺乏OCV实时测试和校准的可实现性。本文提出了一种快速灵活的OCV-SoC提取方法，主要基于锂电池放电过程的OCV-t曲线平滑性假设，利用非支配排序遗传算法（NSGA-II）实现基于任意电流-电压测量数据的OCV-SoC关系提取；随后在UDDS工况下结合拓展卡尔曼滤波（EKF）进行了SoC验证。结果表明，基于平滑性假设可以有效地构建OCV-SoC的映射关联，其中SoC的最大估计误差为2%，且不受滤波器SoC初值的影响。

关键词

锂离子电池; 开路电压; 荷电状态; 开路电压-荷电状态关系; 遗传算法; 平滑性假设

与其它蓄电设备（EES）相比，锂离子电池具有更高的能量密度和功率密度、长日历（calendar）寿命和循环（cycling）寿命、高库伦效率以及低自放电率等优点^［1-3］，因此在混合动力汽车（HEV）、插电式混合动力汽车（PHEV）和纯电动汽车（EV）上得到了广泛应用^［4-6］。实际应用场景下，电动车主要通过电池管理系统（BMS）进行单体电池或电池模组的充放电控制来维持车辆系统安全可靠的运行。而电池的荷电状态（SoC）是BMS系统中一个十分重要的参数，准确的SoC估计有利于BMS实现合理的充放电控制，从而延长电池的使用寿命，同时可防止因电池过充过放而引发的热失控^［7］，以及辅助驾驶员合理地规划行程。

由于SoC难以直接测量得到，因此学者们提出了一系列SoC的估计方法^［8， 9］。其中，基于模型的方法通常具有较高的精度和较强的自校正能力^［10-11］。与复杂的电化学模型相比，等效电路模型（ECM）可以用更少的参数表示电池的动态响应而被广泛应用^［12-14］。典型的ECM通常包括一个串联电压源，用于表示开路电压（OCV），以及一个电阻和一个或多个电阻-电容器（RC）支路^［15-18］。其中，OCV是SoC估计的关键特征变量，通常被视为SoC的函数^［19-20］。电阻和RC参数可以通过最小二乘方法^［17，21］、遗传算法^［22，23］离线识别。此外，递归最小二乘方法^［13，15］常用于ECM参数的实时更新，由此提高SoC的估计精度。值得注意的是，在这些参数识别算法中，OCV- SoC关系通常被视为先验知识，需要通过额外的OCV测试确定。

由于OCV需要在电池平衡状态下测量得到，且电池存在迟滞效应，因此采用传统的方法来确定OCV和SoC之间的映射关系需要较长的时间，例如采用递增OCV（IO）测试^［24］和小电流OCV（LO）测试^［25］通常需要数天才能获得所需的数据^［26］。此外，由于OCV-SoC映射还会随着电池老化而变化^［27-29］，因此需要反复地测试才能保证OCV-SoC关系的准确性，由此将极大地影响车辆和电池的正常使用。

在以往的一些研究中，OCV有时被视为模型参数的一部分，可在放电过程中通过递归最小二乘法^［30-32］、自适应扩展卡尔曼滤波器^［33］或自适应无迹卡尔曼滤波器^［34］进行识别。然而，在SoC的估计过程中，上述研究仍然需要通过增量OCV测试来确定OCV-SoC的关联。同时，由于OCV和模型参数在识别过程中存在相互耦合的关系，因此其识别结果将在很大程度上受到模型初值和建模误差的影响。此外，即使在同一充放电周期内，电池的欧姆阻抗和低频阻抗也会随环境温度和SoC产生明显的变化^［35］，由此导致参数辨识的鲁棒性因车辆实际运行工况的差异而显著降低。

本文提出了一种快速灵活的OCV-SoC映射提取方法，主要基于锂电池放电过程的OCV-t曲线平滑性假设，利用非支配排序遗传算法（NSGA-II）实现了基于任意电流-电压测量数据的OCV-SoC关系提取，随后在UDDS工况下结合拓展卡尔曼滤波（EKF）进行了SoC验证。研究结果表明，基于平滑性假设可以有效地构建OCV-SoC的映射关联，且不受滤波器SoC初值的影响。本研究旨在为电动车复杂运行工况下的在线参数（OCV-SoC）评估和校准提供一种新的解决思路和处理方法。

1 NSGA-II优化提取OCV-SoC映射

图1所示为基于平滑性假设的OCV-SoC映射提取方法流程图。由图可见：首先从测得的电池电压-电流曲线中任取一段电压、电流序列作为输入；然后通过ECM计算得到该段序列对应的电池OCV曲线并评估其平滑性；接着利用NSGA-II根据平滑性指标更新ECM的内部参数，并重新计算OCV；重复进行以上过程实现迭代优化。

图 1  基于平滑性假设的OCV-SoC映射提取流程

Fig. 1  OCV-SoC mapping extraction process based on smoothness assumption

1.1　等效电路模型

由于本方法涉及采用ECM计算电池的OCV曲线，因此ECM的选择将直接影响OCV曲线提取的准确性和效率。一阶等效电路（见图2）是一种在兼顾准确性和模型复杂度的有效分析方法。它通过将复杂电路简化以此来提高模型的分析和设计效率^［36］。研究表明^［36］，尽管二阶等效电路在提高准确性方面有一定效果，但将显著增加计算成本。基于此，本文采用一阶等效电路模型对单体锂电池进行建模。由图2一阶等效电路可得如下电压关系：

图 2  一阶等效电路

Fig. 2  First order equivalent circuit

$${\mathit{U}}_{\mathrm{O}\mathrm{C}}={\mathit{U}}_{\mathit{t}}+{\mathit{U}}_{\mathrm{0}}+{\mathit{U}}_{\mathit{p}}$$

（1）

式中：$U_{\mathrm{O}\mathrm{C}},U_t$， $U_{\mathrm{0}}$和$U_p$分别表示开路电压、端电压、内阻分压和电极极化分压，$R_{\mathrm{0}},C_p,R_p$分别表示等效电路中的电阻、电极极化电容和电极极化电阻。

定义锂电池的系统参数$\theta =[U_{\mathrm{0}},U_p]$。在给定一组等效电路参数$R_{\mathrm{0}},C_p,R_p$并且忽略其在同一充放电循环内随SoC和其他因素变化的前提下，有以下状态转移方程：^［26］

$$\left\{\begin{array}{c}{\dot{U}}_p=\frac{i_{\mathrm{0}}}{C_p}-\frac{U_p}{\tau }\\ U_t=U_{\mathrm{O}\mathrm{C}}-U_p-i_{\mathrm{0}}{R}_{\mathrm{0}}\end{array}\right.$$

（2）

式中：$i_{\mathrm{0}}$为电流负载，$U_t$为测量得到的端电压，$\tau$为RC单元的时间常数，$\tau =R_p{C}_p$。取时间步长为$\mathrm{\Delta }t$，将式（2）离散化后可得：

$$\left\{\begin{array}{c}{U}_{p,k+\mathrm{1}}=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{\Delta }t}{{\tau }_k}\right)U_{p,k}+\left(\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{\Delta }t}{{\tau }_k}\right)\right)R_p{i}_{\mathrm{0},k}\\ U_{t,k}=U_{\mathrm{o}\mathrm{c},k}-U_{p,k}-i_{\mathrm{0},k}{R}_{\mathrm{0},k}\end{array}\right.$$

（3）

1.2　平滑性优化理论推导与分析

放电开始前将电池静置，一段时间后可以认为电池系统参数初值${\theta }_{\mathrm{0}}=[\mathrm{0},\mathrm{0}]$。。根据式（3）可知，只要给定一段连续的充电电压-电流序列$[U_{t,\mathrm{0}},\dots ,U_{t,k};i_{\mathrm{0,0}},\dots ,i_{\mathrm{0},k}]$，就可以通过递推计算获得该时段内任意时刻的OCV。然而，当等效电路参数$R_{\mathrm{0}},C_p,R_p$未知时，使用错误的参数计算OCV将导致其曲线波动剧烈，如图3所示。这是因为基于错误参数的等效电路模型无法正确计算各个电路部分的分压变化。而对于同一段电流电压信号和相同的等效电路模型，使用正确的参数却可以得到相对平滑的曲线。因此，OCV-t曲线的平滑性能一定程度上验证等效电路参数的正确性。基于这一观察，本文提出了基于平滑性假设的等效电路参数求解方法，并在下文中给出具体论证过程。

图 3  一段放电过程内的OCV-t曲线

Fig. 3  OCV-t curve in a discharge process

注：   浅色线表示错误的ECM参数计算结果，深色线表示正确的ECM参数计算结果。

由式（3）可以得到：

$$\begin{array}{l}∆U_{\mathrm{o}\mathrm{c},k}=U_{\mathrm{o}\mathrm{c},k+\mathrm{1}}-U_{\mathrm{o}\mathrm{c},k}=\\ ∆U_{t,k}+∆i_{\mathrm{0},k}{R}_{\mathrm{0}}+\left(\frac{i_{\mathrm{0},k}}{C_p}-\frac{U_{p,k}}{\tau }\right)∆t\end{array}$$

（4）

式中：$∆U_{\mathrm{o}\mathrm{c},k}$表示$U_{\mathrm{o}\mathrm{c}}$在第$k$和$k+\mathrm{1}$个时间步之间的差值，如图4所示。OCV通常可以看作SoC的单调函数^［37］，在电池的整个放电过程中，OCV呈递减趋势（见图3）。通过累加SoC放电从100%至0%整个过程的$∆U_{oc}$，可以得到：

$${\sum }_{\mathrm{s}\mathrm{o}\mathrm{c}=\mathrm{1}}^{\mathrm{s}\mathrm{o}\mathrm{c}=\mathrm{0}}\left|∆U_{\mathrm{o}\mathrm{c}}\right|=C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)$$

（5）

其中：$U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}$和$U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}$分别表示电池的最高截止开路电压和最低截止开路电压；系数$C_{\mathrm{B}\mathrm{R}}$（$\ge \mathrm{1}$）用于补偿放电过程中制动回收（brake recovery）导致的电池电动势的反向变化（升高），如$C_{\mathrm{B}\mathrm{R}}=\mathrm{1}$表示放电过程无制动回收。

图 4  ΔU_oc 示意图

Fig.4  ΔU_oc schematic diagram

由式（4）、式（5）可以得到：

$$\sum \left|∆U_t+∆i_{\mathrm{0}}{R}_{\mathrm{0}}+\left(\frac{i_{\mathrm{0}}}{C_p}-\frac{U_p}{\tau }\right)∆t\right|=C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)$$

（6）

根据测量值可以直接计算得到$∆U_t$和$∆i_{\mathrm{0}}$，因此等式的左边可以视为以等效电路参数$R_{\mathrm{0}}\mathrm{、}{R}_p\mathrm{、}{C}_p$为自变量的函数。为了简化表示，定义$X={\left[R_{\mathrm{0}},R_p,C_p\right]}^{\mathrm{T}}$，则式（6）可表示为：

$$f\left(X\right)=C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)$$

（7）

此时，引入新的方程（8）：

$$f\left(X\right)-C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)=\mathrm{0}$$

（8）

正如前文所述，由于等效电路参数$X$是未知的，因此在实际求解OCV的过程中只能使用$X$的估计值，即$\widehat{X}={\left[\widehat{R_{\mathrm{0}}},\widehat{R_p},\widehat{C_p}\right]}^{\mathrm{T}}$，因此可将方程（8）重新表达为：

$$f\left(\widehat{X}\right)-C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)\approx \mathrm{0}$$

（9）

为了便于表述，此处定义：

$$F\left(\widehat{X}\right)=f\left(\widehat{X}\right)-C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)$$

（10）

$F\left(\widehat{X}\right)$表示估计的OCV曲线和理想曲线的偏离程度，可被视为OCV曲线的平滑性指标。其取值主要受以下因素的影响： ① 一阶等效电路仅使用一个RC单元来简化迟滞效应和极化现象，因此可能引入一定的模型误差； ② 等效电路参数估计的不准确引起的OCV曲线计算误差； ③ 传感器采样频率低导致的“频率失真”（frequency distortion）。基于此，可通过等效电路参数的优化以降低上述第二个因素对$F\left(\widehat{X}\right)$的贡献，从而使$F\left(\widehat{X}\right)$更接近0。注意到$F\left(\widehat{X}\right)\ge \mathrm{0}$，因此本文的目标是将$F\left(\widehat{X}\right)$最小化，即等效电路参数$X$的求解可以表示为：

$$X=\mathrm{A}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\left(f\left(\widehat{X}\right)-C_{\mathrm{B}\mathrm{R}}\left(U_{\mathrm{o}\mathrm{c},\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}-U_{\mathrm{o}\mathrm{c},\mathrm{l}\mathrm{o}\mathrm{w}}\right)\right)=\mathrm{A}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}(f\left(\widehat{X}\right))$$

（11）

由此，将式（11）描述的优化问题代入NSGA-II^［38］（Nondominated Sorting Genetic Algorithm II）算法，便可以得到最优的等效电路参数，从而实现OCV曲线的最优提取。NSGA-II是一种遗传算法，旨在解决单/多目标优化问题，是NSGA的改进版。上述OCV-SoC的最优提取过程可总结如图1所示。

2 优化过程

2.1　实验设置和数据来源

本实验使用斯坦福大学能源控制实验室提供的电池老化数据集^［39］。该数据集使用了石墨/硅阳极和镍-锰-钴阴极的INR21700-M50T电池，采用了UDDS（urban dynamometer driving schedule）放电工况和恒流（CC）-恒压（CV）充电协议，同时以不同的充电速率（从C/4到3C）进行了为期23个月的电池老化实验。实验在恒温环境中（23℃）进行，分别对10个电池进行了测试。实验中采用参考性能测试（RPTs）定期评估电池的老化表现，包括容量、混合脉冲功率表征（HPPC）和电化学阻抗谱（EIS）测试。

需要注意的是，该数据集并没有提供SoC和真实值。因此，本文采用安时积分（AHC）来计算SoC的参考值，定义放电截止时刻的SoC为0。式（12）所示为AHC的计算公式，其中$\eta$为库伦效率，$C_{\mathrm{m}\mathrm{a}\mathrm{x}}$表示从放电截止时刻到电池充满电时对应的电荷量。

$$\mathrm{S}\mathrm{o}{\mathrm{C}}_k=\mathrm{S}\mathrm{o}{\mathrm{C}}_{k-\mathrm{1}}+\frac{\eta }{C_{\mathrm{m}\mathrm{a}\mathrm{x}}}{I}_{k-\mathrm{1}}∆T$$

（12）

2.2　等效电路参数的优化求解

初值和取值范围通常对优化问题的求解至关重要。根据先验知识水平可将此类问题分为以下两种情况：几乎没有先验知识，需要广泛地搜索参数空间；根据已知参数的近似值以此加速优化的收敛。实验设置一和实验设置二（见表1）分别针对上述两种情况进行设计，并使用NSGA-II优化算法求解等效电路参数。

表 1  等效电路参数初值和范围设定

Tab. 1  Initial value and range setting of equivalent circuit parameters

实验设置	变量名	初值	取值范围
实验设置一	$R_0$	1×10-5 Ω	（0， 1］
	$R_p$	1×10-5 Ω	（0， 1］
	$C_p$	1 F	（0， 1×105］
实验设置二	$R_0$	0.028 5 Ω	（0， 0.1］
	$R_p$	0.012 7 Ω	（0， 0.1］
	$C_p$	1 270 F	（0， 1×104］

图5和表2展示了目标函数随迭代次数的变化，从中可以发现不同初值条件下仍能得到相近的最优解，然而合适的初值可显著加速NSGA-II的收敛，从而减少所需的迭代次数（见图6）。在此优化过程中算法所探索的所有参数组合在图6中给出。电池常规检测时可用上一次检测时标定的参数初始化，以获得快速和稳定的计算结果。

图 5  目标函数随迭代次数的变化曲线

Fig.5  The variation of optimization objective with the number of iterations

表2  目标函数随迭代次数的变化

Tab.2  The change of the objective function with the number of iterations

实验设置	迭代次数	$R_0$	$R_p$	$C_p$	目标函数值
实验设置一	1	0.001 2	0.160 4	971	186.85
	25	0.029 2	0.015 6	1 427	6.59
	50	0.028 6	0.011 7	1 429	5.13
	100	0.028 6	0.013 1	1 297	5.01
实验设置二	1	0.027 2	0.012 0	7 217	9.51
	25	0.028 6	0.014 5	1 332	5.02
	50	0.028 5	0.013 4	1 308	5.01
	100	0.028 5	0.012 7	1 273	5.00

图 6  NSGA-II的求解过程

Fig. 6  The solving process of NSGA-II

(a) 经任意初始化（实验一）； (b) 经恰当初始化（实验二）

图7展示了优化目标在变量约束空间内的取值情况。可以观察到，优化目标对$R_{\mathrm{0}}$的取值非常敏感。这是因为电阻是线性元件，它的值决定了ECM的高频响应特性。由于$R_{\mathrm{0}}$可根据电池电压的瞬时变化值较容易地进行计算，因此可以将其固设为已知，由此ECM变量数将减少为两个，即$R_p$和$C_p$。此时函数值与两个变量的关系如图7（b）所示，在可行域内，目标函数有且仅有一个局部极值。由此说明在已知$R_{\mathrm{0}}$的情况下，此优化问题具有良好的可解性和全局最优解的存在性。

图 7  优化目标在变量约束空间内的取值

Fig. 7  Optimize values of the target in the variable constraint space

(a) 优化目标在变量约束空间内的所有取值 (b) 优化目标F与变量R₁、C₁的关系

3 优化结果验证

3.1　等效电路参数验证

首先，通过傅里叶变换可以获取放电工况（UDDS）内电流信号的频域特性，如图8所示。由图可以观察到，由于放电起始阶段为持续约250 0 s的恒流工况，因此在0 Hz处存在一个较大的幅值。除此之外，频率成分主要分布在0~0.2 Hz的范围内。其次，在EIS测试数据中筛选频率在0~0.2 Hz的部分，并使用优化得到的等效电路参数计算对应频率的阻抗，将其绘制在奈奎斯特图上如图 9所示。由图可以看到，优化求解的等效电路参数能够较好地拟合主频范围内的阻抗特性。

图 8  一个放电周期内的电流信号的时频特性

Fig. 8  Time/frequency characteristics of the current signal in a discharge cycle

(a) 时域特性 (b) 频域特性

图 9  奈奎斯特图：样本点来自EIS测试，曲线是优化得到的等效电路的计算结果

Fig. 9  Nyquist diagram: The sample points are derived from EIS tests, and the curves are the calculated results of the optimized equivalent circuits

3.2　OCV-SoC曲线预处理

图10（d）展示了优化计算得到的OCV-t曲线，其与SoC-t曲线呈现近似的下降趋势。相比于直接测得的端电压U_t-t曲线（图10（b）），OCV-t曲线在平滑性方面得到了显著的改善，这表明等效电路模型和优化得到的等效电路参数能够很好地拟合电池的动态响应特性。然而，OCV-t曲线的平滑性还受其他诸多因素的影响（见1.2节），由此导致优化提取的OCV-t曲线仍然存在一些小幅度振荡，从而影响基于OCV-t曲线的SoC预测的稳定性。基于此，本文进一步使用八次多项式差值对OCV-SoC曲线进行平滑性处理，如图 11所示。

图10  优化计算得到的OCV-t曲线

Fig.10  Optimized OCV-t curve obtained from calculation

(a) 电流-时间曲线 (b) SoC参考值⁃时间曲线 (c) 电压⁃时间曲线 (d) OCV估计值⁃时间曲线

图11  提取的OCV-SoC曲线

Fig.11  Extracted OCV - SoC curve

3.3　结合EKF的SoC估计

基于提取所得的OCV-SoC曲线，本研究进一步结合EKF^［26］对下一个循环的放电过程进行SoC估计如图12所示，可以看到SoC的最大估计误差基本控制在2%以内。其中，将干扰的EKF工况定义为在其中某一时刻（图12（b）中t=10 000 s附近），人为地将滤波器内部存储的SoC重置为0.8（偏移量约0.4），以评估优化提取的OCV曲线对EKF估计值的修正性能，如图12（b）所示。结果表明，基于平滑性优化提取的OCV曲线可以快速有效地修正滤波器的SoC估计，即在约1000 s后SoC重新回到了真值附近，由此说明本研究提出的OCV-SoC提取方法在抗传感器噪音或故障方面具有较强的鲁棒性。

图 12  EKF放电过程的SoC估计

Fig. 12  SoC estimation results of EKF discharge process

（a) 无干扰 (b) 有干扰

4 结语

本研究基于OCV曲线的平滑性假设理论推导，结合优化求解器NSGA-II将等效电路参数的求解和OCV曲线的提取转化为多维空间的全局最优问题，提出了一种可基于任意电压电流测量序列准确地提取OCV-SoC曲线的方法，由此避免了耗时的传统OCV测试。经复杂的UDDS工况验证，本方法的OCV提取结果有较高的精度及鲁棒性，即SoC的最大估计误差为2%，且不受滤波器SoC初值的影响，同时在信号干扰下可快速有效地对滤波器（EKF）的SoC估计进行修正。此外，本方法无须专用仪器设备获取数据，易于集成在车载系统之内，可望为电动车的在线参数校准提供实际意义上的参考。