基于IBOA-PF的锂电池健康状态预测

李 鹏，李立伟，杨玉新

(1青岛大学电气工程学院；2青岛大学威海创新研究院；3青岛大学图书馆，山东 青岛 266071)

近年来，锂电池因其高功率承受力、高能量比、高额定电压、低自放电率及绿色环保等特点，广泛应用于航空航天、电动汽车、智能电网等领域。在实际使用过程中，锂电池会随着充放电循环发生性能劣化现象[1]。电池的健康状态(state of health，SOH)是衡量性能劣化程度的直接指标，研究锂电池的健康状态有助于及时更换性能劣化的单体电池，使整个电池组工作在最佳状态。

目前提出的SOH 估计方法主要有三种：①实验估计法，如容量测量法、电阻测量法、差分分析法；②基于模型的自适应滤波法，如卡尔曼、粒子滤波、最小二乘法；③数据驱动法，如人工神经网络、支持向量回归。其中，粒子滤波简单易实现，而且不受线性约束误差和噪声约束条件的影响，适用于分析非线性、非高斯的系统，在SOH 估计领域得到广泛应用[2]。但粒子滤波也存在着一些缺陷，例如部分粒子的权值会随迭代次数增大，多次迭代后只剩下极少数大权值的粒子，即粒子权值退化现象；而采用重采样法复制大权值粒子来解决权值退化问题，又会使得部分有效样本遗失，降低了采样粒子的多样性，导致粒子样本贫化[3]。因此，大量学者通过不同的策略对传统的粒子滤波进行改进。吕佳志等[4]提出通过无迹变换将粒子转移，一定程度上改善了粒子权值退化问题。刘刚等[5]提出结合权值选优的粒子滤波算法，该算法基于动态分簇结构的特点，利用簇和簇之间的传递关系来改善粒子分布，但没有解决样本贫化问题。

此外，使用元启发式智能优化算法来改进粒子滤波成为最近的热门方案，朱超等[6]提出一种基于混沌萤火虫改进粒子滤波算法，利用混沌系统生成初始粒子，将萤火虫的寻优机制引入基本粒子滤波。姬鹏等[7]提出模拟退火粒子群优化粒子滤波，并将其用于车辆状态估计研究。该类方案通过迭代寻优来改善粒子分布，使其集中到更符合真实后验分布的区域，而不是直接舍弃小权值粒子，因此可以解决粒子样本贫化的问题。但元启发式智能优化算法也存在盲目搜索，导致早熟或欠收敛的问题以及未权衡好探索与开发之间的关系，导致陷入局部最优的缺陷[8]。

刘云涛[9]提出用基本蝴蝶算法优化粒子滤波算法，张威虎等[10]在此基础上提出了改进蝴蝶算法优化粒子滤波，引入吸引半径参数降低算法的复杂度，但未能解决蝴蝶算法易陷入局部最优和开发能力差的问题。针对这些问题，本文用混沌数替代固定的切换概率，避免算法陷入局部最优，并引入共生生物搜索的互生阶段，弥补了蝴蝶算法开发能力差的局限性，提高了算法的收敛速度。然后将改进蝴蝶优化算法结合到粒子滤波的重要性采样阶段，有效改善了传统粒子滤波(particle filter，PF)算法粒子多样性丧失的问题，且能够在锂电池SOH 估计中优化粒子的分布。实验部分，本文基于双指数经验模型和时间指标构建状态空间模型，用单纯形法改进高斯牛顿法对空间模型进行初始参数拟合，然后用提出的联合算法对锂电池的容量衰退曲线进行追踪，估计每个充放电循环周期内的电池健康状态，通过仿真实验验证了该方法有效且具有较高的精确度和较好的适应性。

1 主要算法

1.1 粒子滤波

粒子滤波是一种基于贝叶斯滤波和蒙特卡洛采样的统计滤波算法，其主要思想是从重要性概率密度中取样产生一系列带权重的粒子集合，来近似系统状态的后验概率密度函数，并根据算法通过预测和更新两个步骤来递推更新粒子位置和权值，通过求取所有粒子的加权和来替代积分运算，从而得到最小方差估计[11]。

假设非线性系统的状态空间模型为

式中，xk、yk分别为状态值和量测值；f、h 分别为系统的状态函数和量测函数；wk、vk分别为过程噪声和量测噪声。

实际应用中通常选取先验概率密度为重要性概率密度函数，即

这样式(4)可化简为

1.2 基本蝴蝶优化算法

蝴蝶优化算法(butterfly optimization algorithm，BOA)是Arora等[12]学者提出的一种新型元启发式全局优化算法，与其他群智能优化算法一样，该算法的灵感来源也是自然界生物群体行为模式。蝴蝶优化算法的主要思想近似于蝴蝶觅食行为，每只蝴蝶拥有不同的适应度，散发不同程度的香味。蝴蝶的香味表达式为

式中，γ为感知模态；I为刺激强度，与蝴蝶的适应度相关；α为依赖于香味的幂指数。α和γ通常在[0,1]之间取值。

每只蝴蝶都能通过自身感受器感知到其他蝴蝶的香味大小，并向着香味更大的蝴蝶移动。当蝴蝶移动到新的位置时，它拥有的适应度发生改变，散发的香味也会随之变化，经过数次迭代，蝴蝶会移动到最优位置附近。

蝴蝶的移动分为两种模式，当蝴蝶能感受到最优蝴蝶的香味时，它会向着最优蝴蝶移动，这被称为全局搜索；当无法感受到最优解时，它将会随机移动，这被称为局部搜索。蝴蝶的移动方式使用切换概率p 来控制，当p 大于随机数r 时进行全局搜索，否则进行局部搜索。

全局搜索和局部搜索的公式分别为

式中，β为固定的参数；Ngen为总迭代次数。

1.3 改进蝴蝶优化算法

任何元启发式算法都包括探索和开发两个过程，探索表示搜索整个解空间以找到最优解，而开发表示在已获得的解附近寻找更好的解[13]。蝴蝶优化算法的全局搜索对应探索过程，局部搜索对应开发过程。当两个过程达到平衡时，算法才能获得良好的性能，若探索过程占主导地位，算法会浪费大量时间搜索无效区域，导致收敛速度慢，反之，若开发过程占主导地位，种群会失去多样性，导致过早收敛。理想的算法应该兼顾早期阶段的高探索能和后期阶段的高开发能力[14]。所以为了提高性能，研究人员通常都会对基本算法进行一定的改进。改进方法包括融合两种或多种局部搜索技术，或使用自适应参数来修改搜索机制，或将一种算法的最优阶段嵌入另一种算法[15]。

蝴蝶优化算法存在易陷入局部最优和开发能力差的局限性，因此本节提出一种改进蝴蝶优化算法(improved butterfly optimization algorithm，IBOA)，从以下两个方面进行改进：①使用混沌映射替换固定的切换概率p，避免陷入局部最优；②使用共生生物搜索改进全局搜索阶段，提高BOA 的开发能力。

1.3.1 混沌映射改进切换概率p

基本BOA通过固定的切换概率p 与随机数r 比较，来决定探索还是开发，这种随机选择有时会导致BOA 迷失方向并偏离最佳解。混沌映射是提高元启发式算法性能的最佳方法之一[16]，将混沌引入到BOA中有助于算法改善探索与开发之间的权衡，从而更动态和全局地探索解空间，避免陷入局部最优解。

混沌表示系统中的无序或噪声状态，映射表示将算法中的混沌行为与函数中的参数进行关联，混沌图是在非线性系统中显示复杂和动态行为的图[17]。本文从10种最广泛使用的一维混沌图中选择了虫口映射(logistic map)，用来生成一组随迭代次数变化的混沌数pk，从而替代基本BOA 中固定的切换概率p。虫口映射的公式如下

当参数P设为4，p的初始值设为0.8时，根据虫口映射生成的切换概率p如图1所示。

图1 虫口映射生成切换概率Fig.1 switching probability generated by logistic map

1.3.2 共生生物搜索改进全局搜索

共生生物搜索(symbiosis organisms search，SOS)算法是Cheng 等[18]2014 年提出的一种算法，该算法包括三个阶段，即共生、互生和寄生阶段。共生阶段一个粒子更新，另一个粒子既不更新也不降级受损；互生阶段两个粒子相互更新受益；寄生阶段一个粒子通过降解其他粒子来更新。其中互生阶段以其开发能力强而闻名，将其嵌入BOA 的全局搜索阶段，在探索到的新值周围开发更优值，可以提升后者的开发能力[19]。SOS互生阶段的计算式如下

1.3.3 IBOA具体实现步骤

改进后的蝴蝶优化算法步骤如下：

步骤1：初始化种群及各项参数；

步骤2：由式(10)计算每只蝴蝶的香味，并计算全局最优解；

步骤3：由式(14)计算混沌数作为切换概率p，并生成随机数r。

当r＜p 时，根据式(11)进行全局搜索，随机选择蝴蝶h，由式(15)计算蝴蝶i和蝴蝶h的关系特征及互生搜索后两只蝴蝶的新值，并与原值比较选择更优蝴蝶。

当r＞p时，根据式(12)进行局部搜索。

步骤5：根据式(13)更新感知模态γ的值。

步骤6：判断是否达到最大迭代次数，如果达到则输出最优值和最优解，如果不是则跳转到步骤2。

1.3.4 改进蝴蝶优化算法与粒子滤波结合

IBOA-PF算法的核心思想是，把粒子滤波算法重要性采样得到的初始粒子作为初始蝴蝶种群，蝴蝶被香味吸引向香味更大的区域移动，类似于粒子改善为更符合真实后验分布的新值，从而得到更精确的估计值。

将IBOA 引入粒子滤波存在两个问题，一是蝴蝶的适应度值即刺激强度的选取，二是IBOA 优化后粒子是否要进行重采样。

针对前一个问题，文中引进最新时刻的观测值，将刺激强度I定义为

式中，R为量测噪声vk的方差；yi为i时刻的测量值，yi,pred为i时刻测量值的估计值。

针对重采样问题，引入有效粒子数Neff和粒子阈值Nth，若Neff小于Nth则进行重采样。有效粒子数Neff的表达式为

IBOA-PF的具体实现步骤如下：

步骤1：初始化PF 的参数，将先验概率密度作为重要性函数，并进行重要性采样得到采样粒子{xi=1,2,…,N}；

步骤2：初始化IBOA 的参数，将采样粒子作为初始蝴蝶种群；

步骤3：进行步骤2～6，对初始蝴蝶种群进行优化，迭代后的蝴蝶种群为{,i = 1,2,...,N}，将其作为样本粒子输出；

步骤4：根据式(6)计算每个粒子的权值，并由式(8)进行归一化处理；

步骤5：根据式(17)计算有效粒子数Neff，并与设定的阈值Nth进行比较，若Neff小于Nth则进行重采样；

步骤6：由式(9)输出t时刻的估计值。

2 实验与分析

为了将本文提出的改进算法用于锂电池的SOH 估计，首先要选择电池数据集，提取充放电实验数据，然后构建直接指标模型作为粒子滤波的系统状态方程，构建间接指标模型作为系统观测方程，再用IBOA-PF 算法对构建出来的系统状态空间模型进行滤波计算，得到SOH的估计值。

2.1 锂电池数据集

本文采用NASA 数据集中B0005、B0006、B0007号电池老化数据作为实验对象。该数据集对额定容量为2 A·h 的3 种不同型号的18650 锂电池进行连续的充电、放电和阻抗测试实验。充电过程先用1.5 A电流恒流充电，当电压达到4.2 V后进行恒压充电，直到电流下降到20 mA时结束。放电过程用2 A 电流恒流放电，直到电压下降到截至电压，3种型号电池对应的截止电压分别为2.7、2.5、2.2 V。充放电实验时每隔一定的时间间隔记录电池的电压、电流、温度数据。阻抗测试实验中，用0.1～5 kHz频率扫描测量电池阻抗。

反复进行上述测试，直到电池容量退化到失效阈值，即额定容量的70%时停止实验。3种型号电池的容量变化曲线如图2所示。由图可知出，3种型号的电池都呈现出显著容量衰退现象，且均存在局部容量再生现象，造成这种现象的原因是锂电池内部的活性材料被激活。本文选择退化曲线变化最明显，容量再生最多的B0006号电池作为研究对象。

图2 3种型号电池的容量变化曲线Fig.2 Capacity decline curves of three batteries

2.2 构建状态空间模型

2.2.1 系统状态方程

电池衰减最直观的表现是最大可用容量的减小，因此SOH 可定义为每次充放电循环，电池从充满电状态以恒流放电到截至电压，释放出的最大容量与新电池的额定容量的比值[20]，表达式为

故电池的当前容量可作为SOH 估计的直接指标，双指数经验退化模型能够较好地拟合非线性系统，可以用来表示电池容量[21]，公式如下

式中，k 为充放电次数；ak、bk、ck、dk为第k次模型参数。

令Xk=[ak,bk,ck,dk]，且

根据式(18)～式(20)可构建系统的状态方程为

2.2.2 系统观测方程

根据粒子滤波的相关定义，在得到当前时刻的预测值后，必须用该时刻的观测值来进行修正，因此需要一项可以表征SOH 变化的间接指标。Liu等[22]研究发现各充放电周期内，电池当前的SOH与放电电压曲线存在着一定关系，因此可以用提取到的电压和时间数据来构建与SOH 相关的时间指标TI。TI 的定义为不同充放电周期的恒流放电期间，电池变化相同的电压所对应的时间间隔，即等放电电压差时间间隔，表达式为

式中，TkV,max和TkV,min分别是电压下降到采样最大值和最小值对应的时间。时间指标TI旨在用片段放电时间去推断整个放电时间，电压区间越小该指标的实用性越高，但同时精确度会下降。本文设定电压采样区间的最大值为4 V，最小值为3 V，从B0006号电池提取出的时间指标TI如图3所示。TI与SOH的关系可表示为

图3 B0006号电池的时间指标TIFig.3 Time index of B0006 battery

或

其中，q1、q2、q3和p1、p2是各自的模型参数。

式(23)可用来衡量构建的时间指标TI 与SOH的相关性，式(24)作为系统的观测方程，即

综上，构建的系统状态空间模型为

式中，SOHk、TIk分别为状态估计值和观测值；f、g分别为系统的状态函数和量测函数；wk、vk分别为过程噪声和量测噪声。

2.3 参数拟合

SOH 估计中的模型参数可用非线性最小二乘法得出，如二阶梯度法、高斯牛顿法、阻尼牛顿法等，其中高斯牛顿法原理简单且精度较高，经常被用来拟合参数。但高斯牛顿法存在受初值影响较大的缺点，因此本文提出一种改进高斯牛顿法，将单纯形法与高斯牛顿法结合。单纯形法在拟合时可以大范围搜索，虽然精度较低但受初值影响很小[23]，将单纯形法得到的结果作为高斯牛顿法的初始值，可以有效改善后者因初始值选择不当导致不能收敛的缺点，同时大大减少迭代次数，收敛效果好。

改进高斯牛顿法的步骤如下：

步骤1：初始化，设定初始值，设有n个参数，则定义n+1个点；

步骤2：计算各点对应的目标函数值进行比较；

步骤3：按照以下规则，去除坏点，替换新点；

规则1：去除最坏点，替换为对称反射点；规则2：若新点仍为最坏点，则将次坏点替换为对称反射点；规则3：若某一点始终保留，则缩短步长，以该点为新起点。替换公式为

步骤4：达到迭代次数，或得到满意结果时停止，否则返回步骤2；

步骤5：输出当前最优点作为高斯牛顿法的初始值，x0=(l1l2…ln)；

步骤6：求误差ek和雅可比矩阵Jk

步骤7：求解增量方程，得到增量Δxi

步骤8：令xi+1=xi+Δxi，如果Δxi足够小，则停止迭代输出xi+1，否则返回步骤6。

用改进高斯牛顿法对B0006 号电池的容量、SOH、TI 模型进行参数拟合，得到的各参数见表1。

表1 模型参数Table 1 Model parameters

2.4 相关性分析

B0006 号电池双指数模型拟合函数如图4 所示。拟合曲线与容量退化曲线之间的皮尔逊相关系数R=0.9886，说明双指数模型能够表征B0006 号电池容量衰退过程。

图4 B0006号电池双指数模型拟合Fig.4 Double exponential model of B0006 battery

根据式(23)及拟合出的q1、q2、q3参数可得TI 到SOH 的映射图像如图5、图6 所示给出了每次循环两者的误差绝对值。两者的最大误差为0.0058，皮尔逊相关系数R=0.9998，因此提取的时间指标TI 可以作为SOH 的间接指标，用来修正预测值。

2.5 实验设计

本文使用的仿真实验平台为MATLAB R2019a，仿真程序中算法使用的各参数如下：粒子总数N为100，有效粒子数阈值Nth为2N/3，状态噪声方差Qw和观测噪声方差Qv均为0.001，改进蝴蝶算法的迭代次数Ngen为300，感知模态γ0=0.01，幂指数α为0.1，常数β为0.025。

图5 SOH与TI映射比较Fig.5 Comparison of SOH with TI mapping

图6 SOH与TI映射误差Fig.6 Error of SOH with TI mapping

仿真实验步骤为：首先从NASA数据集中提取B0006号电池每个充放电循环周期的容量、时间指标数据，根据式(26)构建系统状态空间模型，并用改进高斯牛顿法进行参数拟合，得到模型中的各个初始参数，最后用IBOA-PF 算法进行估计，得到SOH的预测值。

2.6 结果分析

分别用PF、BOA-PF、IBOA-PF 对B0006 号电池的前100 次充放电循环周期进行SOH 估计，得到的结果如图7 所示。图8 给出了3 种算法的估计误差，可以看出IBOA-PF 算法具有较小的估计误差。

本文采用最大相对误差(maximum relative error，MRE)、均方根误差(root mean square error，RMSE)和平均绝对百分比误差(mean absolute percentage error，MAPE)3 个 指 标 来 评估预测结果。MRE、RMSE 和MAPE 的计算式如下

图7 3种算法对SOH的估计Fig.7 SOH estimation with three algorithms

图8 SOH估计误差Fig.8 Error of SOH estimation

各算法估计结果定量分析见表2。

从表2 中可以看出，IBOA-PF 对SOH 的估计最大相对误差在5%以内，且3 项评价指标均优于PF和BOA-PF。

表2 B0006号电池SOH估计比较Table 2 Comparison of SOH estimation of B0006 battery

为验证所提出方法的适应性，使用同样的方法对B0007 号电池进行实验，得到的估计结果分析见表3。

表3 B0007号电池SOH估计比较Table 3 Comparison of SOH estimation of B0007 battery

由表3可知，IBOA-PF在不同型号的电池中仍取得了较好的估计效果，证明本文提出的方法具有较好的适应性和稳定性。

3 结 论

为了提高预测的准确率，本工作提出了一种基于IBOA-PF 的锂电池SOH 估计方法。IBOA-PF 算法的主要思想是用改进的蝴蝶算法优化粒子滤波的重要性采样过程，能有效改善传统PF 算法粒子多样性丧失的问题。在实验阶段，该方法基于NASA电池数据集，首先构建B0006号电池的双指数退化模型，提取时间指标TI即等放电电压差时间间隔作为观测值，然后用改进高斯牛顿法进行参数拟合，并根据所建立的状态空间模型，使用IBOA-PF 对电池的容量衰退进行追踪，估计每次充放电循环周期内的SOH。通过MATLAB 仿真实验表明，该方法在预测锂电池SOH 上可行，精确度优于传统的PF 算法，且适用于不同型号的电池，具有较好的适用性和稳定性，对电池健康管理系统具有重要意义。

需要说明的是，本文用到的双指数经验模型追踪电池容量退化的精度不高，可能会影响SOH 的估计结果，另外，研究中用到的时间指标TI是在恒流放电条件下提取的，但在实际运行中，锂电池不会始终以恒流状态工作，这也可能会带来一定的估计误差。本文接下来将从这两方面展开研究，进一步提高SOH预测的精确度和稳定性。