單成鑫，李立偉，楊玉新

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

近年來，能源危機和環(huán)境污染的影響，使得電動汽車的研究和應用實現(xiàn)了迅速發(fā)展[1]。鋰離子電池憑借能量比高、壽命長、自放電率低的特點，在電動汽車領域得到了廣泛應用。由于鋰離子電池的失效可能會導致嚴重的后果，為了保障電池的安全應用，提出了一種對電池狀態(tài)進行實時監(jiān)控的電池管理系統(tǒng)[2](battery management system，BMS)。電池荷電狀態(tài)[3](state of charge，SOC)是BMS 的重要參數(shù)，由于SOC 是電池的內部狀態(tài)，不能直接測量。因此，SOC的估計成為了許多研究者的關注的問題，目前已經提出了許多估計方法，可以分為三大類：開環(huán)法、數(shù)據(jù)驅動法和基于模型的方法[4]。開環(huán)法包括安時積分法和開路電壓法，數(shù)據(jù)驅動的方法主要包括神經網(wǎng)絡[5]模糊算法，卡爾曼濾波[6]和粒子濾波[7](particle filter，PF)是基于模型方法中常見的算法。上述估算方法的優(yōu)缺點見表1。

表1 電池SOC估算方法的優(yōu)缺點Table 1 Advantages and disadvantages of battery SOC estimation method

PF 是一種高精度的基于概率的估計算法，PF利用蒙特卡洛方法采樣的一組加權粒子來近似系統(tǒng)的后驗分布，對分布的形式沒有任何明確的假設，該方法適用于具有較強非線性[8]的復雜模型的狀態(tài)估計。對于非高斯分布系統(tǒng)噪聲條件下具有復雜形式的非線性鋰離子電池模型，利用PF 進行荷電狀態(tài)估計，具有更好的優(yōu)勢。傳統(tǒng)PF 算法在迭代求解的過程中，往往會出現(xiàn)粒子多樣化降低，產生粒子貧化現(xiàn)象，造成估算結果產生較大誤差降低算法的估計精度[9]。近年來智能優(yōu)化算法在處理優(yōu)化問題中具有獨特優(yōu)勢，可以用于解決粒子濾波粒子貧化的問題[10]。蟻群算法(ant colony optimization，ACO)是一種類似于遺傳算法的迭代優(yōu)化方法[11]。該算法在許多領域都有良好的優(yōu)化性能，在PF 算法過程中，采用ACO 算法在更新步驟之前對粒子分布進行優(yōu)化，從而使粒子貧化問題最小化。本文采用改進ACO算法，通過改進傳統(tǒng)ACO算法中的信息素濃度和啟發(fā)式的相對影響參數(shù)以及信息素揮發(fā)速率，解決該算法的局部最優(yōu)問題，再對PF 進行優(yōu)化，改善PF中粒子貧化問題。

1 電池模型及參數(shù)辨識

1.1 鋰離子電池模型

等效電路模型的建立，不僅考慮其精確性，也要兼顧復雜性。電池的充放電動態(tài)性能由RC 環(huán)節(jié)描述，并且考慮了電池的極化現(xiàn)象。隨著更多RC環(huán)節(jié)的累積，模型更準確，但更復雜。考慮到模型的精確性和簡單性，選擇二階Thevenin 電池等效模型來描述鋰離子電池。模型電路如圖1 所示，UOC為電池開路電壓；R0表示電池歐姆內阻；UL為負載電壓；R1、R2、C1、C2為極化內阻、極化電容；I為電路中流過的電流，兩個RC環(huán)節(jié)的環(huán)路電壓分別為U1、U2。

圖1 二階Thevenin等效模型Fig.1 Second-order Thevenin equivalent model

電池SOC 選取安時積分法進行計算，并根據(jù)基礎電路定律，對等效電路模型進行分析，可列出以下關系式

式中，SOC0表示為SOC 的初始值；I 表示瞬時電池端電流，在放電過程中被認為是正的，充電過程則為負；t 為時間；η 為庫侖效率；QN為電池額定容量；UL、I由試驗測量得出。

對式(1)～(3)進行離散化處理，選取SOC作為狀態(tài)變量，電池端電壓UL作為觀測變量，得到電池的非線性狀態(tài)空間模型。

狀態(tài)方程

觀測方程

式(5)中

式中，?t 為采樣間隔；ω、v 為過程和觀測噪聲；UOC(SOCt)所表示的UOC與SOC 的函數(shù)擬合關系以及兩個RC環(huán)節(jié)中R1、C1、R2、C2參數(shù)，將由下一節(jié)的參數(shù)辨識給出。

1.2 參數(shù)辨識

鋰電池的測試平臺由直流電子負載、dSPACE MicroLabBox試驗箱以及計算機組成。本次試驗采用額定電壓3.7 V，額定容量為1.5 A·h的鋰電池組在室溫25°C下進行測試。對選取的二階Thevenin電池等效模型通過試驗測得的數(shù)據(jù)采用指數(shù)擬合方法進行參數(shù)辨識[12]。

首先將電池充滿電后靜置1 h，測得開路電壓。然后采用1 C的間歇脈沖電流對電池進行放電180 s，然后靜置180 s，間歇周期為360 s。每次放電過程，電池的SOC減少5%。并將靜置結束時的電壓值記錄下來，近似為開路電壓。放電過程中電壓的部分曲線變化如圖2所示。

圖2 脈沖放電電壓部分波形Fig.2 Partial waveform of pulse discharge voltage

通過圖2 可以發(fā)現(xiàn)在放電過程中，在AB 段與CD段電池電壓出現(xiàn)了跳躍式突變，由等效電路圖可知，C1、C2端電壓不會突變，這是因為電池本身歐姆內阻造成的。在BC段與DE段電池電壓則出現(xiàn)了指數(shù)式緩慢變化，這是因為電池內部的極化現(xiàn)象造成的，對應等效電路圖中兩個RC環(huán)節(jié)影響的結果。

結合脈沖放電試驗數(shù)據(jù)，對電池歐姆內阻與RC環(huán)節(jié)進行參數(shù)辨識。在這里對于歐姆內阻R0的值，選取如圖CD 段電池放電結束時突變的壓降與放電電流之比，即

在脈沖放電后的靜置時間段中(DE 段)，RC 環(huán)節(jié)可看作零輸入響應，電池的輸出方程為

U1(0)為R1C1在放電結束瞬間電壓初始值，U2(0)為放電的初始值。在MATLAB的Curve Filting工具箱中利用式(8)跟DE段電壓數(shù)據(jù)曲線擬合，求出時間常數(shù)τ1、τ2。在脈沖放電過程中BC 段，兩個RC環(huán)節(jié)可看作為零狀態(tài)響應，電池的輸出方程為

將上一步通過擬合求得的時間常數(shù)τ1、τ2代入，再次使用MATLAB Curve Filting 工具箱利用式(9)跟BC段電壓數(shù)據(jù)曲線擬合，得出R1、R2，再根據(jù)τ1=R1C1，τ2=R2C2，分別求出C1、C2。

參數(shù)辨識結果見表2。

表2 二階模型辨識參數(shù)Table 2 Second-order model identification parameter

1.3 OCV-SOC擬合

根據(jù)1.2 節(jié)記錄的放電試驗過程中的近似開路電壓，以及對應SOC 值的試驗數(shù)據(jù)。在MATLAB中的cftool 工具箱進行7 次多項式的曲線擬合，得到開路電壓UOC與SOC的函數(shù)關系式

7次擬合曲線方程的判定系數(shù)為0.9997，效果較好。既能保證擬合精度，又能減輕后續(xù)運行負擔。

2 PF算法

粒子濾波是基于蒙特卡洛方法遞推的貝葉斯濾波[13]，通過一組隨機樣本表示所需的狀態(tài)概率分布，并給出解決貝葉斯估計問題的方法。它通過一組加權粒子來逼近后驗，其形式沒有任何明確的假設，因此可以用于任何非線性非高斯系統(tǒng)。

首先，需要求得后驗分布的最優(yōu)貝葉斯解。假設用于貝葉斯估計的狀態(tài)值(xk)和觀測值(yk)由以下方程式控制。

式中，fk為先前狀態(tài)xk和狀態(tài)噪聲ωk的非線性函數(shù)；hk為狀態(tài)xk和觀測噪聲vk的非線性函數(shù)；ωk、vk為相互獨立的系統(tǒng)狀態(tài)和觀測噪聲。假設初始狀態(tài)p(x0)概率密度函數(shù)已知，貝葉斯估計目的是遞歸計算每個狀態(tài)xk的后驗概率密度p(xk|y1:k)。

假設系統(tǒng)狀態(tài)服從一階馬爾可夫系統(tǒng)模型，得到先驗概率分布的一般表達式，再由此求出后驗概率分布

上式中，分母是一個常數(shù)值

式(13)和式(14)的遞歸關系構成了最佳貝葉斯解的基礎。但是，后驗概率密度的這種遞歸傳播在理論上只是一個最佳解決方案，而且由于式(13)中的積分通常很難處理，因此無法通過求解獲得。所以在粒子濾波算法中使用蒙特卡洛模擬和重要性采樣執(zhí)行遞歸貝葉斯估計，假設已知重要性概率密度函數(shù)q(xk|yk)，從中采樣抽取粒子集{xi0:k}Ni=0，x0:k={xj,j=0,1…,k}，是0～k時刻的狀態(tài)集，其對應權值{wik}Ni=0更新式如下

將權值歸一化得到下式

其中后驗密度通過狀態(tài)空間附近的粒子的相對密度來近似

其中δ(?)為狄拉克函數(shù)。當N足夠大時，近似值(18)就越接近真實后驗密度p(xk|y1:k)。

粒子濾波的具體實現(xiàn)步驟如下。

（1）初始化狀態(tài)：從先驗概率分布p(x0)中隨機選取N 個初始狀態(tài)粒子xi0，i=(1,2,...N)，并將特定粒子的初始粒子權重指定為1/N。粒子數(shù)量取決于系統(tǒng)類型和計算成本。重采樣的閾值可以初始化為Nth=2N/3。

（2）重要性采樣：使用狀態(tài)方程更新狀態(tài)值，并根據(jù)式(16)更新粒子權重。

（3）權重歸一化：使用式(17)將所有粒子的總權重歸一化。

（4）重采樣：重采樣的閾值定義為有效樣本Neff，可以將其定義為

如果有效粒子低于給定的Nth，則對粒子重新采樣以得到權重為1/N的新粒子。

（5）狀態(tài)估計：將狀態(tài)計算為

3 IACO-PF算法

3.1 ACO算法

蟻群算法本質上是一個基于生物的系統(tǒng)，該系統(tǒng)基于模擬螞蟻自然行為的媒介[14]。ACO源自螞蟻展示的自然優(yōu)化能力。當尋找食物時，它們傾向于沿著巢穴i 和一些外部地標j 之間的最佳路線(或路徑)前進，因為如果更多的螞蟻選擇相同的路徑，它們的特定信息素濃度τij就會變得更高。地標離巢穴越近，螞蟻可進行的往返次數(shù)就越高。信息素濃度τij越高，更多的螞蟻會選擇這種途徑而不是其他可能的途徑。該迭代過程實現(xiàn)了端點之間的最佳路徑。因此，ACO 算法使用數(shù)學公式來模擬此自然優(yōu)化過程。螞蟻k(k=1,2,...,m)的移動方向基于式(21)中定義的概率函數(shù)。當一只螞蟻尚未完成移動時，式(21)用于確定要采取的下一條路徑，直到它們收斂到最佳解。

式中，allowedk為可供選擇地標的集合；ηij是可作為先驗啟發(fā)式值，定義為i、j 兩點間距離dij的倒數(shù)，ηij(t)=1/dij。該概率分布受參數(shù)α 和β 的偏倚影響，這些參數(shù)分別確定信息素濃度和啟發(fā)式信息的相對影響。

完成一次循環(huán)后，信息素濃度τij(t)受螞蟻的運動影響，根據(jù)下式更新

式中，ρ是信息素揮發(fā)速率；?τij(t)為此期間信息素濃度的增強值。

一般來說，ACO算法按照以下規(guī)則工作[15]。

（1）對于給定的問題，其候選解與對應于螞蟻正在訪問的路徑或頂點相關聯(lián)。例如，在PF 框架中，如果一個粒子被認為是一只螞蟻，那么，這個粒子(螞蟻)就會在狀態(tài)空間中沿著一條路徑運動，最終到達它的新狀態(tài)。具有新狀態(tài)的粒子被視為頂點。于是新粒子集形成候選解。

（2）當螞蟻沿著一條路徑移動時，沿著該路徑的信息素濃度以一定的水平增加。

（3）當螞蟻需要在兩條或多條路徑中選擇移動時，往往信息素濃度較高、距離較短的路徑被螞蟻選擇的概率較大。

（4）信息素濃度會隨著時間的推移而蒸發(fā)，這一特性會使螞蟻最多的路徑始終保持信息素的最高水平。

3.2 IACO算法

3.2.1 調整信息素濃度與啟發(fā)式相對影響參數(shù)

在計算轉移概率Pkij時，信息素濃度τij隨著迭代次數(shù)的增加而累積，使得τij的影響逐漸增大，期望啟發(fā)式的影響逐漸減小，傳統(tǒng)蟻群算法往往會陷入局部最優(yōu)。參數(shù)的選擇直接影響模型求解的速度和效率。傳統(tǒng)蟻群算法中信息素濃度與啟發(fā)式相對影響參數(shù)α、β 一般取1、2.5，為提高ACO 的效率，定義α、β為

式中，Nc為當前迭代次數(shù)。α 隨著迭代次數(shù)的增加逐漸減小，使得τij的影響逐漸降低，同時β 隨著迭代進行而逐漸增加，從而增加期望啟發(fā)式的影響，避免算法陷入局部最優(yōu)。

3.2.2 改進信息素揮發(fā)速率

蒸發(fā)速率ρ 的選擇對算法的性能有很大的影響。ρ太小會導致信息素的快速積累。因此，也會很容易使該算法陷入局部最優(yōu)。ρ過大會降低算法的效率。為了提高解的質量，提出了非線性增量蒸發(fā)率。ρ的新的更新公式定義為

式中，ρ0為初始定義信息素揮發(fā)速率；Nc,max為最大迭代次數(shù)。

4 IACO-PF算法設計

4.1 算法思想

為了解決傳統(tǒng)粒子濾波中粒子多樣性貧化現(xiàn)象，優(yōu)化粒子濾波的重采樣步驟，在更新步驟之前，采用改進的蟻群優(yōu)化算法[16]。在此之中，螞蟻將替代隨機生成的粒子，它們將根據(jù)下面的蟻群算法重新定位，以便在重采樣步驟中接近最優(yōu)解。

在迭代過程中，τij受每次粒子運動的影響，其值的增加或減少取決于粒子j是否位于運動路徑中，如下式所示

信息素路徑定義意味著路徑中信息素濃度會隨著更多螞蟻通過同一路徑時而增加，而較少螞蟻經過的路徑的信息素濃度將隨時間成比例下降。螞蟻沿著路徑行進并再次返回所花費的時間越長，信息素的蒸發(fā)量就越大。相比之下，由于大多數(shù)螞蟻最終選擇了該路徑，因此信息素的濃度在最短路徑中保持較高。

優(yōu)化步驟基于從式(28)獲得的概率函數(shù)迭代運行。它表示粒子i在N?1個粒子中選擇粒子j作為移動方向的概率。

啟發(fā)式函數(shù)ηij現(xiàn)在定義為兩個粒子間距離的倒數(shù)

參數(shù)τij的初始值為目標粒子j的權重，即

因此，目標粒子j的權值越大，則粒子i向目標粒子j移動的概率就越大，反之越小。α和β是與原始ACO 中定義的參數(shù)相同的參數(shù)，在粒子濾波中定義為確定權值信息和位置信息的相對影響。

當Pij接近1 時，IACO 算法會收斂，這意味著粒子i會重新定位到粒子j的更近距離處。當在每次迭代過程中都保持這種收斂性時，大多數(shù)粒子會收斂到該粒子j，它表示為基于Pij的可能性更高(權重更高)的鄰域。如果直接通過概率式(28)將粒子向最優(yōu)權值的粒子轉移，經過多次迭代后會導致粒子多樣性降低。因此，本文引入距離閾值λL和概率閾值λP來控制較優(yōu)權值粒子間的轉移。IACO 算法直到所有移動的粒子的位置收斂到在式(31)中定義的距離閾值內的高似然區(qū)域(一般或局部最優(yōu)解)時終止。概率閾值的具體取值通過試驗獲得。

式中，wj是目標粒子j的權重；randn是正態(tài)分布的隨機數(shù)。

4.2 算法步驟

PF 算法將根據(jù)以下步驟重新定位，以便在重采樣步驟時接近最優(yōu)解。

（1）初始化粒子狀態(tài)。從先驗概率分布p(x0)中隨機選取N 個SOC 初始狀態(tài)粒子{SOCi0,i=1,2,...,N},粒子的初始權重大小{wik}N

i=0=1/N。

（2）重要性采樣。根據(jù)系統(tǒng)狀態(tài)方程更新粒子，得到下一時刻先驗粒子{SOCik,i=1,2,...,N}，并根據(jù)下列方程更新粒子權重

（3）采用改進蟻群算法優(yōu)化粒子狀態(tài)，當粒子與其目標粒子之間的距離不在閾值λL內，并且迭代時間不超過最大值，則選擇距離在閾值λL以內的粒子i根據(jù)概率Pij選擇移動目標j，設置概率閾值λP，當Pij<λP時，粒子位置保持不變，當Pij>λP時，粒子i轉移至目標粒子j位置。之后更新ACO的參數(shù)(如η、τ等)。

（4）權重歸一化。將采樣到的N個粒子的權重進行歸一化，得到

（5）重采樣。根據(jù)有效粒子數(shù)Neff比較重采樣的閾值，判斷是否進行重采樣。

（6）狀態(tài)輸出。

5 實驗設計與結果分析

5.1 脈沖放電試驗

為了驗證本文提出的IACO-PF 算法的準確性和有效性，首先根據(jù)電池在1C 脈沖放電狀態(tài)的試驗數(shù)據(jù)以及辨識參數(shù)，并代入電池模型的狀態(tài)方程(4)與觀測方程(5)中，分別采用IACO-PF算法和PF算法對鋰離子電池進行SOC 估算，并將這兩種算法的預測SOC值與真實SOC值進行比較。

各算法均采用過程噪聲方差Q=1，觀測噪聲R=0.01，并取SOC 的初始值SOC0=100，采用粒子數(shù)N=100，有效粒子閾值Nth=2N/3。IACO-PF算法中取最大迭代次數(shù)Nc,max=100，ρ0=0.5，概率閾值取λP=0.0025。試驗得到SOC估算結果以及誤差如圖3和圖4所示。

圖3 SOC試驗仿真曲線Fig.3 Experiment simulation curve of SOC

圖4 SOC試驗仿真曲線誤差Fig.4 Experiment simulation curve error of SOC

為了更直觀地對比兩個算法的SOC估算性能，分別計算了兩種算法的均方根誤差和最大絕對誤差來評估預測結果

從圖3、圖4可以看出，由IACO-PF算法與PF算法得到的SOC預測值都與SOC的真實值非常接近，說明了這兩種算法都能對電池的SOC 值做出很好的預測。并且從表3中的數(shù)據(jù)得出IACO-PF算法和PF算法的SOC預測值與真實值的最大誤差都能控制在2.5%以內，說明了對于非高斯分布系統(tǒng)噪聲條件下的非線性鋰離子電池SOC 的估算，粒子濾波算法具有較強的適用性。從圖4可以直觀看出，PF 算法由于產生了粒子貧化問題，使得SOC的估計誤差開始增大，反觀IACO-PF 算法在整個預測過程中誤差值都能保持穩(wěn)定，并且最大誤差減小到了1%以內，由此可以說明IACO-PF算法具有更高的有效性和準確性。

表3 兩種算法的數(shù)據(jù)統(tǒng)計Table 3 Data statistics of two algorithms

5.2 DST工況試驗

模擬電動汽車實際行駛過程中電流變化情況，并進一步驗證IACO-PF 算法的準確性，進行DST方法對電池進行充放電仿真試驗[17]，DST工況電流如圖5所示，試驗前先將電池充分靜置，試驗參數(shù)不變，仿真時間15000 s，試驗得到SOC估算結果以及誤差如圖6和圖7所示。

圖5 DST工況電流分布圖Fig.5 DST operating current profile

圖6 SOC試驗仿真曲線Fig.6 Experiment simulation curve of SOC

圖7 SOC試驗仿真曲線誤差Fig.7 Experiment simulation curve error of SOC

從圖6、圖7可知，在DST工況電流下，IACOPF算法與PF算法得到的SOC預測值同樣也能很好地趨于SOC真實值，做出良好的預測。由表4中的數(shù)據(jù)可知，在復雜的工況電流下，兩種算法的精度會相應降低，但相比于PF 算法，IACO-PF 算法的SOC估計值仍保持良好的有效性和準確性。

表4 兩種算法的數(shù)據(jù)統(tǒng)計Table 4 Data statistics of two algorithms

6 結 論

結合二階Thevenin 電池等效模型搭建鋰電池狀態(tài)方程，對傳統(tǒng)蟻群優(yōu)化算法進行調整信息素濃度與啟發(fā)式相對影響參數(shù)以及改進信息素揮發(fā)速率來解決傳統(tǒng)粒子濾波的粒子貧化問題。試驗結果表明，基于IACO-PF 算法的SOC 估計結果相比傳統(tǒng)PF 算法更具準確和穩(wěn)定，對鋰電池健康管理具有重要意義。擴展了智能算法SOC估計的研究領域，具有廣泛的應用價值。