基于IBOA-PF的鋰電池健康狀態(tài)預測

李 鵬，李立偉，楊玉新

(1青島大學電氣工程學院；2青島大學威海創(chuàng)新研究院；3青島大學圖書館，山東 青島 266071)

近年來，鋰電池因其高功率承受力、高能量比、高額定電壓、低自放電率及綠色環(huán)保等特點，廣泛應用于航空航天、電動汽車、智能電網(wǎng)等領域。在實際使用過程中，鋰電池會隨著充放電循環(huán)發(fā)生性能劣化現(xiàn)象[1]。電池的健康狀態(tài)(state of health，SOH)是衡量性能劣化程度的直接指標，研究鋰電池的健康狀態(tài)有助于及時更換性能劣化的單體電池，使整個電池組工作在最佳狀態(tài)。

目前提出的SOH 估計方法主要有三種：①實驗估計法，如容量測量法、電阻測量法、差分分析法；②基于模型的自適應濾波法，如卡爾曼、粒子濾波、最小二乘法；③數(shù)據(jù)驅動法，如人工神經(jīng)網(wǎng)絡、支持向量回歸。其中，粒子濾波簡單易實現(xiàn)，而且不受線性約束誤差和噪聲約束條件的影響，適用于分析非線性、非高斯的系統(tǒng)，在SOH 估計領域得到廣泛應用[2]。但粒子濾波也存在著一些缺陷，例如部分粒子的權值會隨迭代次數(shù)增大，多次迭代后只剩下極少數(shù)大權值的粒子，即粒子權值退化現(xiàn)象；而采用重采樣法復制大權值粒子來解決權值退化問題，又會使得部分有效樣本遺失，降低了采樣粒子的多樣性，導致粒子樣本貧化[3]。因此，大量學者通過不同的策略對傳統(tǒng)的粒子濾波進行改進。呂佳志等[4]提出通過無跡變換將粒子轉移，一定程度上改善了粒子權值退化問題。劉剛等[5]提出結合權值選優(yōu)的粒子濾波算法，該算法基于動態(tài)分簇結構的特點，利用簇和簇之間的傳遞關系來改善粒子分布，但沒有解決樣本貧化問題。

此外，使用元啟發(fā)式智能優(yōu)化算法來改進粒子濾波成為最近的熱門方案，朱超等[6]提出一種基于混沌螢火蟲改進粒子濾波算法，利用混沌系統(tǒng)生成初始粒子，將螢火蟲的尋優(yōu)機制引入基本粒子濾波。姬鵬等[7]提出模擬退火粒子群優(yōu)化粒子濾波，并將其用于車輛狀態(tài)估計研究。該類方案通過迭代尋優(yōu)來改善粒子分布，使其集中到更符合真實后驗分布的區(qū)域，而不是直接舍棄小權值粒子，因此可以解決粒子樣本貧化的問題。但元啟發(fā)式智能優(yōu)化算法也存在盲目搜索，導致早熟或欠收斂的問題以及未權衡好探索與開發(fā)之間的關系，導致陷入局部最優(yōu)的缺陷[8]。

劉云濤[9]提出用基本蝴蝶算法優(yōu)化粒子濾波算法，張威虎等[10]在此基礎上提出了改進蝴蝶算法優(yōu)化粒子濾波，引入吸引半徑參數(shù)降低算法的復雜度，但未能解決蝴蝶算法易陷入局部最優(yōu)和開發(fā)能力差的問題。針對這些問題，本文用混沌數(shù)替代固定的切換概率，避免算法陷入局部最優(yōu)，并引入共生生物搜索的互生階段，彌補了蝴蝶算法開發(fā)能力差的局限性，提高了算法的收斂速度。然后將改進蝴蝶優(yōu)化算法結合到粒子濾波的重要性采樣階段，有效改善了傳統(tǒng)粒子濾波(particle filter，PF)算法粒子多樣性喪失的問題，且能夠在鋰電池SOH 估計中優(yōu)化粒子的分布。實驗部分，本文基于雙指數(shù)經(jīng)驗模型和時間指標構建狀態(tài)空間模型，用單純形法改進高斯牛頓法對空間模型進行初始參數(shù)擬合，然后用提出的聯(lián)合算法對鋰電池的容量衰退曲線進行追蹤，估計每個充放電循環(huán)周期內的電池健康狀態(tài)，通過仿真實驗驗證了該方法有效且具有較高的精確度和較好的適應性。

1 主要算法

1.1 粒子濾波

粒子濾波是一種基于貝葉斯濾波和蒙特卡洛采樣的統(tǒng)計濾波算法，其主要思想是從重要性概率密度中取樣產(chǎn)生一系列帶權重的粒子集合，來近似系統(tǒng)狀態(tài)的后驗概率密度函數(shù)，并根據(jù)算法通過預測和更新兩個步驟來遞推更新粒子位置和權值，通過求取所有粒子的加權和來替代積分運算，從而得到最小方差估計[11]。

假設非線性系統(tǒng)的狀態(tài)空間模型為

式中，xk、yk分別為狀態(tài)值和量測值；f、h 分別為系統(tǒng)的狀態(tài)函數(shù)和量測函數(shù)；wk、vk分別為過程噪聲和量測噪聲。

實際應用中通常選取先驗概率密度為重要性概率密度函數(shù)，即

這樣式(4)可化簡為

1.2 基本蝴蝶優(yōu)化算法

蝴蝶優(yōu)化算法(butterfly optimization algorithm，BOA)是Arora等[12]學者提出的一種新型元啟發(fā)式全局優(yōu)化算法，與其他群智能優(yōu)化算法一樣，該算法的靈感來源也是自然界生物群體行為模式。蝴蝶優(yōu)化算法的主要思想近似于蝴蝶覓食行為，每只蝴蝶擁有不同的適應度，散發(fā)不同程度的香味。蝴蝶的香味表達式為

式中，γ為感知模態(tài)；I為刺激強度，與蝴蝶的適應度相關；α為依賴于香味的冪指數(shù)。α和γ通常在[0,1]之間取值。

每只蝴蝶都能通過自身感受器感知到其他蝴蝶的香味大小，并向著香味更大的蝴蝶移動。當蝴蝶移動到新的位置時，它擁有的適應度發(fā)生改變，散發(fā)的香味也會隨之變化，經(jīng)過數(shù)次迭代，蝴蝶會移動到最優(yōu)位置附近。

蝴蝶的移動分為兩種模式，當蝴蝶能感受到最優(yōu)蝴蝶的香味時，它會向著最優(yōu)蝴蝶移動，這被稱為全局搜索；當無法感受到最優(yōu)解時，它將會隨機移動，這被稱為局部搜索。蝴蝶的移動方式使用切換概率p 來控制，當p 大于隨機數(shù)r 時進行全局搜索，否則進行局部搜索。

全局搜索和局部搜索的公式分別為

式中，β為固定的參數(shù)；Ngen為總迭代次數(shù)。

1.3 改進蝴蝶優(yōu)化算法

任何元啟發(fā)式算法都包括探索和開發(fā)兩個過程，探索表示搜索整個解空間以找到最優(yōu)解，而開發(fā)表示在已獲得的解附近尋找更好的解[13]。蝴蝶優(yōu)化算法的全局搜索對應探索過程，局部搜索對應開發(fā)過程。當兩個過程達到平衡時，算法才能獲得良好的性能，若探索過程占主導地位，算法會浪費大量時間搜索無效區(qū)域，導致收斂速度慢，反之，若開發(fā)過程占主導地位，種群會失去多樣性，導致過早收斂。理想的算法應該兼顧早期階段的高探索能和后期階段的高開發(fā)能力[14]。所以為了提高性能，研究人員通常都會對基本算法進行一定的改進。改進方法包括融合兩種或多種局部搜索技術，或使用自適應參數(shù)來修改搜索機制，或將一種算法的最優(yōu)階段嵌入另一種算法[15]。

蝴蝶優(yōu)化算法存在易陷入局部最優(yōu)和開發(fā)能力差的局限性，因此本節(jié)提出一種改進蝴蝶優(yōu)化算法(improved butterfly optimization algorithm，IBOA)，從以下兩個方面進行改進：①使用混沌映射替換固定的切換概率p，避免陷入局部最優(yōu)；②使用共生生物搜索改進全局搜索階段，提高BOA 的開發(fā)能力。

1.3.1 混沌映射改進切換概率p

基本BOA通過固定的切換概率p 與隨機數(shù)r 比較，來決定探索還是開發(fā)，這種隨機選擇有時會導致BOA 迷失方向并偏離最佳解?；煦缬成涫翘岣咴獑l(fā)式算法性能的最佳方法之一[16]，將混沌引入到BOA中有助于算法改善探索與開發(fā)之間的權衡，從而更動態(tài)和全局地探索解空間，避免陷入局部最優(yōu)解。

混沌表示系統(tǒng)中的無序或噪聲狀態(tài)，映射表示將算法中的混沌行為與函數(shù)中的參數(shù)進行關聯(lián)，混沌圖是在非線性系統(tǒng)中顯示復雜和動態(tài)行為的圖[17]。本文從10種最廣泛使用的一維混沌圖中選擇了蟲口映射(logistic map)，用來生成一組隨迭代次數(shù)變化的混沌數(shù)pk，從而替代基本BOA 中固定的切換概率p。蟲口映射的公式如下

當參數(shù)P設為4，p的初始值設為0.8時，根據(jù)蟲口映射生成的切換概率p如圖1所示。

圖1 蟲口映射生成切換概率Fig.1 switching probability generated by logistic map

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

共生生物搜索(symbiosis organisms search，SOS)算法是Cheng 等[18]2014 年提出的一種算法，該算法包括三個階段，即共生、互生和寄生階段。共生階段一個粒子更新，另一個粒子既不更新也不降級受損；互生階段兩個粒子相互更新受益；寄生階段一個粒子通過降解其他粒子來更新。其中互生階段以其開發(fā)能力強而聞名，將其嵌入BOA 的全局搜索階段，在探索到的新值周圍開發(fā)更優(yōu)值，可以提升后者的開發(fā)能力[19]。SOS互生階段的計算式如下

1.3.3 IBOA具體實現(xiàn)步驟

改進后的蝴蝶優(yōu)化算法步驟如下：

步驟1：初始化種群及各項參數(shù)；

步驟2：由式(10)計算每只蝴蝶的香味，并計算全局最優(yōu)解；

步驟3：由式(14)計算混沌數(shù)作為切換概率p，并生成隨機數(shù)r。

當r＜p 時，根據(jù)式(11)進行全局搜索，隨機選擇蝴蝶h，由式(15)計算蝴蝶i和蝴蝶h的關系特征及互生搜索后兩只蝴蝶的新值，并與原值比較選擇更優(yōu)蝴蝶。

當r＞p時，根據(jù)式(12)進行局部搜索。

步驟5：根據(jù)式(13)更新感知模態(tài)γ的值。

步驟6：判斷是否達到最大迭代次數(shù)，如果達到則輸出最優(yōu)值和最優(yōu)解，如果不是則跳轉到步驟2。

1.3.4 改進蝴蝶優(yōu)化算法與粒子濾波結合

IBOA-PF算法的核心思想是，把粒子濾波算法重要性采樣得到的初始粒子作為初始蝴蝶種群，蝴蝶被香味吸引向香味更大的區(qū)域移動，類似于粒子改善為更符合真實后驗分布的新值，從而得到更精確的估計值。

將IBOA 引入粒子濾波存在兩個問題，一是蝴蝶的適應度值即刺激強度的選取，二是IBOA 優(yōu)化后粒子是否要進行重采樣。

針對前一個問題，文中引進最新時刻的觀測值，將刺激強度I定義為

式中，R為量測噪聲vk的方差；yi為i時刻的測量值，yi,pred為i時刻測量值的估計值。

針對重采樣問題，引入有效粒子數(shù)Neff和粒子閾值Nth，若Neff小于Nth則進行重采樣。有效粒子數(shù)Neff的表達式為

IBOA-PF的具體實現(xiàn)步驟如下：

步驟1：初始化PF 的參數(shù)，將先驗概率密度作為重要性函數(shù)，并進行重要性采樣得到采樣粒子{xi=1,2,…,N}；

步驟2：初始化IBOA 的參數(shù)，將采樣粒子作為初始蝴蝶種群；

步驟3：進行步驟2～6，對初始蝴蝶種群進行優(yōu)化，迭代后的蝴蝶種群為{,i = 1,2,...,N}，將其作為樣本粒子輸出；

步驟4：根據(jù)式(6)計算每個粒子的權值，并由式(8)進行歸一化處理；

步驟5：根據(jù)式(17)計算有效粒子數(shù)Neff，并與設定的閾值Nth進行比較，若Neff小于Nth則進行重采樣；

步驟6：由式(9)輸出t時刻的估計值。

2 實驗與分析

為了將本文提出的改進算法用于鋰電池的SOH 估計，首先要選擇電池數(shù)據(jù)集，提取充放電實驗數(shù)據(jù)，然后構建直接指標模型作為粒子濾波的系統(tǒng)狀態(tài)方程，構建間接指標模型作為系統(tǒng)觀測方程，再用IBOA-PF 算法對構建出來的系統(tǒng)狀態(tài)空間模型進行濾波計算，得到SOH的估計值。

2.1 鋰電池數(shù)據(jù)集

本文采用NASA 數(shù)據(jù)集中B0005、B0006、B0007號電池老化數(shù)據(jù)作為實驗對象。該數(shù)據(jù)集對額定容量為2 A·h 的3 種不同型號的18650 鋰電池進行連續(xù)的充電、放電和阻抗測試實驗。充電過程先用1.5 A電流恒流充電，當電壓達到4.2 V后進行恒壓充電，直到電流下降到20 mA時結束。放電過程用2 A 電流恒流放電，直到電壓下降到截至電壓，3種型號電池對應的截止電壓分別為2.7、2.5、2.2 V。充放電實驗時每隔一定的時間間隔記錄電池的電壓、電流、溫度數(shù)據(jù)。阻抗測試實驗中，用0.1～5 kHz頻率掃描測量電池阻抗。

反復進行上述測試，直到電池容量退化到失效閾值，即額定容量的70%時停止實驗。3種型號電池的容量變化曲線如圖2所示。由圖可知出，3種型號的電池都呈現(xiàn)出顯著容量衰退現(xiàn)象，且均存在局部容量再生現(xiàn)象，造成這種現(xiàn)象的原因是鋰電池內部的活性材料被激活。本文選擇退化曲線變化最明顯，容量再生最多的B0006號電池作為研究對象。

圖2 3種型號電池的容量變化曲線Fig.2 Capacity decline curves of three batteries

2.2 構建狀態(tài)空間模型

2.2.1 系統(tǒng)狀態(tài)方程

電池衰減最直觀的表現(xiàn)是最大可用容量的減小，因此SOH 可定義為每次充放電循環(huán)，電池從充滿電狀態(tài)以恒流放電到截至電壓，釋放出的最大容量與新電池的額定容量的比值[20]，表達式為

故電池的當前容量可作為SOH 估計的直接指標，雙指數(shù)經(jīng)驗退化模型能夠較好地擬合非線性系統(tǒng)，可以用來表示電池容量[21]，公式如下

式中，k 為充放電次數(shù)；ak、bk、ck、dk為第k次模型參數(shù)。

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

根據(jù)式(18)～式(20)可構建系統(tǒng)的狀態(tài)方程為

2.2.2 系統(tǒng)觀測方程

根據(jù)粒子濾波的相關定義，在得到當前時刻的預測值后，必須用該時刻的觀測值來進行修正，因此需要一項可以表征SOH 變化的間接指標。Liu等[22]研究發(fā)現(xiàn)各充放電周期內，電池當前的SOH與放電電壓曲線存在著一定關系，因此可以用提取到的電壓和時間數(shù)據(jù)來構建與SOH 相關的時間指標TI。TI 的定義為不同充放電周期的恒流放電期間，電池變化相同的電壓所對應的時間間隔，即等放電電壓差時間間隔，表達式為

式中，TkV,max和TkV,min分別是電壓下降到采樣最大值和最小值對應的時間。時間指標TI旨在用片段放電時間去推斷整個放電時間，電壓區(qū)間越小該指標的實用性越高，但同時精確度會下降。本文設定電壓采樣區(qū)間的最大值為4 V，最小值為3 V，從B0006號電池提取出的時間指標TI如圖3所示。TI與SOH的關系可表示為

圖3 B0006號電池的時間指標TIFig.3 Time index of B0006 battery

或

其中，q1、q2、q3和p1、p2是各自的模型參數(shù)。

式(23)可用來衡量構建的時間指標TI 與SOH的相關性，式(24)作為系統(tǒng)的觀測方程，即

綜上，構建的系統(tǒng)狀態(tài)空間模型為

式中，SOHk、TIk分別為狀態(tài)估計值和觀測值；f、g分別為系統(tǒng)的狀態(tài)函數(shù)和量測函數(shù)；wk、vk分別為過程噪聲和量測噪聲。

2.3 參數(shù)擬合

SOH 估計中的模型參數(shù)可用非線性最小二乘法得出，如二階梯度法、高斯牛頓法、阻尼牛頓法等，其中高斯牛頓法原理簡單且精度較高，經(jīng)常被用來擬合參數(shù)。但高斯牛頓法存在受初值影響較大的缺點，因此本文提出一種改進高斯牛頓法，將單純形法與高斯牛頓法結合。單純形法在擬合時可以大范圍搜索，雖然精度較低但受初值影響很小[23]，將單純形法得到的結果作為高斯牛頓法的初始值，可以有效改善后者因初始值選擇不當導致不能收斂的缺點，同時大大減少迭代次數(shù)，收斂效果好。

改進高斯牛頓法的步驟如下：

步驟1：初始化，設定初始值，設有n個參數(shù)，則定義n+1個點；

步驟2：計算各點對應的目標函數(shù)值進行比較；

步驟3：按照以下規(guī)則，去除壞點，替換新點；

規(guī)則1：去除最壞點，替換為對稱反射點；規(guī)則2：若新點仍為最壞點，則將次壞點替換為對稱反射點；規(guī)則3：若某一點始終保留，則縮短步長，以該點為新起點。替換公式為

步驟4：達到迭代次數(shù)，或得到滿意結果時停止，否則返回步驟2；

步驟5：輸出當前最優(yōu)點作為高斯牛頓法的初始值，x0=(l1l2…ln)；

步驟6：求誤差ek和雅可比矩陣Jk

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

步驟8：令xi+1=xi+Δxi，如果Δxi足夠小，則停止迭代輸出xi+1，否則返回步驟6。

用改進高斯牛頓法對B0006 號電池的容量、SOH、TI 模型進行參數(shù)擬合，得到的各參數(shù)見表1。

表1 模型參數(shù)Table 1 Model parameters

2.4 相關性分析

B0006 號電池雙指數(shù)模型擬合函數(shù)如圖4 所示。擬合曲線與容量退化曲線之間的皮爾遜相關系數(shù)R=0.9886，說明雙指數(shù)模型能夠表征B0006 號電池容量衰退過程。

圖4 B0006號電池雙指數(shù)模型擬合Fig.4 Double exponential model of B0006 battery

根據(jù)式(23)及擬合出的q1、q2、q3參數(shù)可得TI 到SOH 的映射圖像如圖5、圖6 所示給出了每次循環(huán)兩者的誤差絕對值。兩者的最大誤差為0.0058，皮爾遜相關系數(shù)R=0.9998，因此提取的時間指標TI 可以作為SOH 的間接指標，用來修正預測值。

2.5 實驗設計

本文使用的仿真實驗平臺為MATLAB R2019a，仿真程序中算法使用的各參數(shù)如下：粒子總數(shù)N為100，有效粒子數(shù)閾值Nth為2N/3，狀態(tài)噪聲方差Qw和觀測噪聲方差Qv均為0.001，改進蝴蝶算法的迭代次數(shù)Ngen為300，感知模態(tài)γ0=0.01，冪指數(shù)α為0.1，常數(shù)β為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數(shù)據(jù)集中提取B0006號電池每個充放電循環(huán)周期的容量、時間指標數(shù)據(jù)，根據(jù)式(26)構建系統(tǒng)狀態(tài)空間模型，并用改進高斯牛頓法進行參數(shù)擬合，得到模型中的各個初始參數(shù)，最后用IBOA-PF 算法進行估計，得到SOH的預測值。

2.6 結果分析

分別用PF、BOA-PF、IBOA-PF 對B0006 號電池的前100 次充放電循環(huán)周期進行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 項評價指標均優(yōu)于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在不同型號的電池中仍取得了較好的估計效果，證明本文提出的方法具有較好的適應性和穩(wěn)定性。

3 結 論

為了提高預測的準確率，本工作提出了一種基于IBOA-PF 的鋰電池SOH 估計方法。IBOA-PF 算法的主要思想是用改進的蝴蝶算法優(yōu)化粒子濾波的重要性采樣過程，能有效改善傳統(tǒng)PF 算法粒子多樣性喪失的問題。在實驗階段，該方法基于NASA電池數(shù)據(jù)集，首先構建B0006號電池的雙指數(shù)退化模型，提取時間指標TI即等放電電壓差時間間隔作為觀測值，然后用改進高斯牛頓法進行參數(shù)擬合，并根據(jù)所建立的狀態(tài)空間模型，使用IBOA-PF 對電池的容量衰退進行追蹤，估計每次充放電循環(huán)周期內的SOH。通過MATLAB 仿真實驗表明，該方法在預測鋰電池SOH 上可行，精確度優(yōu)于傳統(tǒng)的PF 算法，且適用于不同型號的電池，具有較好的適用性和穩(wěn)定性，對電池健康管理系統(tǒng)具有重要意義。

需要說明的是，本文用到的雙指數(shù)經(jīng)驗模型追蹤電池容量退化的精度不高，可能會影響SOH 的估計結果，另外，研究中用到的時間指標TI是在恒流放電條件下提取的，但在實際運行中，鋰電池不會始終以恒流狀態(tài)工作，這也可能會帶來一定的估計誤差。本文接下來將從這兩方面展開研究，進一步提高SOH預測的精確度和穩(wěn)定性。