莫仕勛， 楊 皓， 蔣坤坪， 梁振燊

(廣西大學電氣工程學院， 廣西 南寧 530004)

1 引言

變壓器作為電力系統(tǒng)中的關(guān)鍵設(shè)備，其建模方法一直是學者研究的熱點，而目前電力系統(tǒng)各類設(shè)備的建模趨于精細化，對變壓器的建模要求也越來越高[1，2]。

變壓器建模的一個關(guān)鍵點在于磁滯現(xiàn)象的表示，由于磁滯特性在實際運行中會影響到波形形狀，經(jīng)過傅里葉變換后得出的各次諧波含量都會有很大差異，所以會對變壓器保護的二次諧波判據(jù)等閉鎖條件產(chǎn)生影響，可見，磁滯能否正確且準確地在變壓器建模中得以表征，無論對于電力系統(tǒng)穩(wěn)態(tài)暫態(tài)的精確計算，還是繼電保護的準確整定都有較大影響。

在目前的研究中，磁滯模型也分為多種，而相比于SW(Stoner-Wolhfarth)模型，Preisach模型等[3]，J-A(Jlies-Atherton)磁滯模型具有參數(shù)較少、物理意義清晰等優(yōu)點，在電磁材料磁特性模擬研究中被廣泛應(yīng)用[4]。

J-A模型由代數(shù)-微分方程組描述，其中有5個未知參數(shù)，未知參數(shù)的大小直接決定了磁滯回線的形狀，從而影響到變壓器模型的精確度[5]。無論是需要利用J-A建模還是在仿真軟件中利用帶磁滯的變壓器，都需要用到這5個參數(shù)，而要得到這5個參數(shù)，就要通過已有的曲線數(shù)據(jù)進行參數(shù)辨識。

目前，智能算法在優(yōu)化領(lǐng)域是研究的熱點之一，使用智能優(yōu)化算法對微分方程的參數(shù)進行辨識也具有優(yōu)越性，在磁滯參數(shù)辨識領(lǐng)域，模擬退火算法(Simulated Annealing, SA)[6]及遺傳算法(Genetic Algorithm,GA)等是較為常用的智能算法。利用智能優(yōu)化算法進行參數(shù)辨識，或者將智能算法與傳統(tǒng)的方法結(jié)合[7]具有更高的精度，所以得到了更廣泛的應(yīng)用。

雖然這些算法能基本辨識出所需的參數(shù)，但是精度依然不高，并且往往需要迭代多次，消耗很長時間才能尋找到最優(yōu)解，也有可能陷入局部最優(yōu)，精確度和實用性均達不到要求，在此基礎(chǔ)上結(jié)合傳統(tǒng)方法進行改進的文獻[7]中算法復雜，步驟較為繁瑣。針對這些問題，本文提出了一種改進禿鷹搜索(Bald Eagle Search， BES)算法，以實現(xiàn)對 J-A 磁滯模型參數(shù)的快速、精確辨識。該方法將Tent混沌映射、精英反向?qū)W習策略與自適應(yīng)權(quán)重融入BES算法，將待求值快速收斂于全局最優(yōu)值。結(jié)合實驗數(shù)據(jù)，使用該算法對 J-A 磁滯模型參數(shù)進行辨識，然后基于所得到的參數(shù)對磁性材料的磁滯特性進行模擬，并將模擬結(jié)果與其他算法的模擬結(jié)果及實驗結(jié)果進行了對比分析。結(jié)果表明，使用該算法辨識模型參數(shù)時收斂速度更快、辨識精度更高，基于該算法辨識參數(shù)的模擬磁滯回線與實測曲線吻合較好，驗證了該算法的實用性和有效性。

2 J-A磁滯模型

磁滯模型有Preisach模型、J-A模型及SW模型等，其中由量子理論描述的J-A磁滯模型以其物理意義明確、參數(shù)易于獲取等優(yōu)點得到廣泛應(yīng)用[8]，J-A模型中，無磁滯磁化曲線由改進的Langevin函數(shù)表示：

(1)

式中，Man為無磁滯磁化強度；He為有效磁場強度；Ms為飽和磁化強度；a為無磁滯磁化曲線的形狀參數(shù)；He和磁場強度H又有如下關(guān)系：

He=H+αM

(2)

式中，α為平均場參數(shù)，反應(yīng)磁疇間的耦合程度；磁化強度M又可分為可逆磁化強度Mrev及不可逆磁化強度Mirr兩部分：

M=Mrev+Mirr

(3)

可逆分量Mrev、不可逆分量Mirr與無磁滯磁化強度Man之間存在如下關(guān)系：

Mrev=c(Man-Mirr)

(4)

式中，c為磁化因數(shù)。

最終根據(jù)能量平衡方程可得：

(5)

式中，δ為方向系數(shù)，當dH/dt>0時取1，dH/dt<0時取-1；k為牽制系數(shù)；μ0為真空磁導率。

式(1)～式(5)即為J-A理論描述的磁滯現(xiàn)象，求出dM/dH后，由式(6)便可求得鐵心磁導率μFe：

(6)

由以上分析可知，J-A磁滯模型共包含5個未知參數(shù)，分別為Ms，α，a，c，k。對于用于制造變壓器鐵心的硅鋼片，已知的數(shù)據(jù)往往為此硅鋼片的B-H或M-H磁滯回線，若需建模，需使用的是5個未知參數(shù)。如此問題便轉(zhuǎn)化為微分方程組的參數(shù)辨識問題，以實測值與仿真值的均方根誤差Fitness作為適應(yīng)度函數(shù)，表達式如下：

(7)

式中，n為B-H磁滯回線所包含點的個數(shù)；g為其中的第g個點；Bmea為硅鋼片已知的B-H曲線的磁感應(yīng)強度；Bcal為通過已知的H數(shù)據(jù)由J-A模型計算得到的磁感應(yīng)強度值。式(7)的值越小，說明計算出的點與實際的點越接近，磁滯曲線越吻合。此問題便轉(zhuǎn)化為最小值優(yōu)化問題，利用優(yōu)化算法不斷優(yōu)化便能達到最優(yōu)參數(shù)。

3 禿鷹搜索算法

禿鷹主要生活在美洲地區(qū)，有出色的視力和觀察技能，主要選擇魚作為食物。在捕食獵物時，禿鷹會從各種方面評估狩獵的成本及成功概率，再進行移動與捕捉。其首先根據(jù)個體和種群對鮭魚的集中程度選擇搜索空間，向特定區(qū)域飛行;其次，在選定的搜索空間內(nèi)對水面進行搜索，找到合適的獵物后，禿鷹逐漸改變高度，迅速俯沖，成功地從水中捕獲鮭魚等獵物[9]。

馬來西亞學者Alsattar 受到禿鷹捕食過程的啟發(fā)，于2020年提出了一種新型智能算法——禿鷹搜索優(yōu)化算法, 此算法以上述禿鷹的捕食習慣為基礎(chǔ), 將禿鷹的捕食過程分為三個階段[10], 數(shù)學模型為：

(1)確定搜索區(qū)域：禿鷹在棲息地通過判斷獵物在各個區(qū)域的密度來確定最佳搜尋位置，便于搜索獵物,并通過觀察結(jié)果來確定位置變更，更新由隨機搜索的先驗信息Pi,new乘以β來確定。該階段的數(shù)學模型為:

Pi,new=Pbest+β·rand·(Pmean-Pi)

(8)

式中，β為控制位置變化參數(shù)，變化區(qū)間為(1.5,2)；rand為 (0,1) 間隨機數(shù)；Pbest為當前禿鷹搜索確定的最佳搜索位置；Pmean為先前搜索結(jié)束后禿鷹的平均分布位置；Pi為第i只禿鷹位置。

(2)搜索空間獵物：在上一步選定的搜索區(qū)域中，禿鷹呈螺旋形狀盤旋飛行，尋找最佳俯沖捕獲位置，螺旋飛行數(shù)學模型采用極坐標方程表示,位置更新公式如式(9)～式(12)所示:

θ(i)=Pbest+bπ·rand

(9)

r(i)=θ(i)+R·rand

(10)

(11)

(12)

式中，θ(i)與r(i)分別為螺旋運動方程的極角與極徑；b、R控制螺旋運動軌跡，變化區(qū)間分別為(0,5)、(0.5,2)；x(i)與y(i)為極坐標中禿鷹位置。禿鷹位置按式(13)更新：

Pi,new=Pi+x(i)(Pi-Pmean)+y(i)(Pi-Pi+1)

(13)

(3)向獵物發(fā)起猛撲：確定獵物位置后，禿鷹從搜索空間的最佳位置快速俯沖飛向目標獵物, 此時種群其他個體也跟著向最佳位置移動并攻擊獵物, 這一階段的運動方程為式(14)～式(16)、式(12)。

θ(i)=bπ·rand

(14)

r(i)=θ(i)

(15)

(16)

俯沖中禿鷹位置更新公式如式(17)、式(18)所示：

(17)

Pi,new=rand·Pbest+δx+δy

(18)

式中，c1與c2為禿鷹向最佳與中心位置的運動強度，取值范圍均為(1,2)。

根據(jù)多個測試函數(shù)表明，BES算法具有較強的全局搜索能力, 能夠有效地解決各類復雜數(shù)值優(yōu)化問題，因此可以式(7)作為適應(yīng)度函數(shù)，不斷尋優(yōu)，理論上可以得到J-A磁滯模型參數(shù)。但同時，標準的BES算法和其他元啟發(fā)式算法一樣，存在收斂速度慢、解的搜索精度相對較低且容易陷入局部最優(yōu)值的缺點[11]，需要針對這些缺點對算法進行改進。

4 BES算法的改進

4.1 Tent混沌映射

混沌是自然界中的一種現(xiàn)象，智能算法初始化種群時使用rand函數(shù)隨機生成種群個體，無法保證種群個體分布在各個位置，即無法保證種群的多樣性，降低了算法的效率，且更容易陷入局部最優(yōu)。而采用混沌映射方法可以解決這一問題。在多種混沌映射方法中，Tent混沌映射各方面性能優(yōu)于其他一些常用混沌映射方法，如Logistic 混沌映射[12]，可以用于產(chǎn)生優(yōu)化算法的混沌序列[13]。Tent映射結(jié)構(gòu)簡單，且映射呈現(xiàn)的結(jié)果分布密度比較均勻，具有很好的遍歷性，因此本文采用Tent混沌映射增加BES算法的初始種群多樣性。

Tent 混沌映射的表達式如式(19)所示：

(19)

式中，sn為第n次運算求出的混沌序列值；are用于調(diào)整范圍，對are與s的初始值進行調(diào)整，就可以得到在不同范圍內(nèi)的按Tent混沌分布的隨機值，從而增加BES 算法的初始種群多樣性，求出Tent混沌序列后，利用混沌序列對初始種群進行映射，即將式(19)產(chǎn)生序列的集合載波到待求解空間，最終的種群生成公式如式(20)所示：

odim,new=omin+(omax-omin)odim

(20)

式中，dim為混沌序列的維度；odim,new為混沌序列載波到待求解空間后的新值；odim為式(19)產(chǎn)生的混沌序列；omax和omin分別為序列中變量的最大值、最小值。

使用 Tent混沌映射得到的新值xdim,new作為 BES算法初始種群的分布，增加了種群的多樣性，提高了BES算法的全局搜索能力。

4.2 精英反向?qū)W習

即便BES算法能有不錯的效果，也存在一定的局限性[10]，禿鷹選擇搜索空間利用上個階段的可用信息來確定下個搜索區(qū)域, 在隨機選擇另一個搜索區(qū)域時, 會依據(jù)之前的搜索域來確定相應(yīng)的鄰域進而完成本次禿鷹第一階段的搜索過程, 因此針對禿鷹選擇搜索空間而言, 如果禿鷹群體在上一次的迭代過程陷入了局部最優(yōu)狀態(tài), 那么便無法跳出此局部最優(yōu)，整個尋優(yōu)過程便無法更為準確, 即無法在某類特定優(yōu)化問題中求得最優(yōu)解, 降低了啟發(fā)式智能優(yōu)化算法的最優(yōu)化效果。

反向?qū)W習機制(Opposition-Based Learning，OBL)可以有效增加種群的多樣性和質(zhì)量，目前已被較好地應(yīng)用于多種算法的改進中[14，15]，反向?qū)W習機制首先計算當前解的反向解，然后從當前解和其反向解的種群中選取最優(yōu)解更新個體。反向?qū)W習策略的定義如下：

反向?qū)W習的原理為，對于一個解來說，其自身和其反向解總有一個更接近最優(yōu)解，反向解的引入，可以擴大算法的搜索區(qū)域，但是反向?qū)W習存在其生成的反向解可能比當前搜索空間更難搜索到最優(yōu)值的問題，并且對那些原解適應(yīng)度值大于反向解適應(yīng)度值的個體，對其進行反向區(qū)域的搜索，則相當于在做無用功，而對原解適應(yīng)度值小于反向解適應(yīng)度值的個體，對其進行反向區(qū)域的搜索價值要高于其領(lǐng)域的開發(fā)價值。

針對反向?qū)W習的缺點，精英反向?qū)W習(Elite Opposition-Based Learning，EOBL)被提出，由前文分析可知，精英個體比普通個體攜帶更多有效信息，EOBL利用這一特點，通過種群中精英個體形成反向種群，然后從反向種群和當前種群中選取優(yōu)秀個體構(gòu)成新的種群[16]，避免了對低價值信息的反復計算及利用，所以相對于反向?qū)W習來說，精英反向?qū)W習提高了效率。精英反向?qū)W習在每次迭代中增加了種群多樣性與局部搜索能力，可以有效地避免算法落入局部最優(yōu)解。

精英反向解的定義如下：

(21)

4.3 自適應(yīng)權(quán)重

慣性權(quán)重是智能算法改進的重要方法，迭代初期，為了快速達到最優(yōu)解附近，對全局搜索能力的要求高些；迭代后期，需要對區(qū)域進行精細的搜索，所以對局部搜索能力的要求高些。如果按照基礎(chǔ)的智能算法，每一步的搜索方式都是一樣的，無法達到以上要求，降低了算法的效率。借鑒慣性權(quán)重的思想，在禿鷹位置更新方式中繼續(xù)引入動態(tài)權(quán)重因子，使其在迭代初期具有較大的值，能夠更好地進行全局探索，在迭代后期自適應(yīng)地減小，從而更好地進行局部搜索，同時提高收斂速度[17]。

自適應(yīng)權(quán)重的公式如式(22)所示：

(22)

式中，t為算法本次迭代的迭代次數(shù)；max_iter為算法的最大迭代次數(shù)。

改進后的位置更新由式(18)變?yōu)槭?23)所示：

Pi,new=rand·Pbest+ω(δx+δy)

(23)

經(jīng)過以上三點改進得到新算法，在此命名為改進禿鷹搜索(Improved Bald Eagle Search，IBES)算法。

5 J-A模型參數(shù)辨識及比較

5.1 J-A模型參數(shù)辨識流程

根據(jù)上文所提算法改進，可得以下具體步驟：

(1)初始化參數(shù)，如種群數(shù)量、最大迭代次數(shù)等，并利用Tent混沌映射初始化禿鷹種群。

(2)結(jié)合自適應(yīng)權(quán)重因子計算每只禿鷹的適應(yīng)度值，找出當前最優(yōu)適應(yīng)度值和最差適應(yīng)度值，以及相對應(yīng)的位置。

(3)從適應(yīng)度值較優(yōu)的個體中，選取部分作為精英個體。

(4)對精英個體求其反向解，計算精英反向解所對應(yīng)的適應(yīng)度值及相應(yīng)的位置。

(6)將反向解對應(yīng)的適應(yīng)度與原來的解進行比較與排序，更新出較優(yōu)的解。

(7)判斷是否達到結(jié)束條件，若是，則進行下一步，否則跳轉(zhuǎn)步驟(2)。

(8)程序結(jié)束，輸出最優(yōu)結(jié)果。

J-A模型參數(shù)辨識流程圖如圖1所示。

圖1 IBES算法流程圖Fig.1 IBES algorithm flow chart

5.2 參數(shù)辨識及算法的比較

為驗證IBES算法在 J-A 磁滯模型參數(shù)辨識問題上的正確性與有效性。本文使用文獻[5]的參數(shù)生成一條仿真曲線作為基準磁滯曲線，包含500個數(shù)據(jù)點。

SA是J-A模型參數(shù)辨識的常用方法，本文采用SA、BES、IBES分別對磁滯回線進行參數(shù)辨識并進行對比，IBES和BES種群數(shù)量為30，迭代次數(shù)為100次，SA算法最大溫度為100，最大容忍度為10-8。

本仿真測試環(huán)境為：操作系統(tǒng) Windows 10，CPU為AMD A10-9630P RADEON R5, 2.60 GHz，內(nèi)存為4 GB，仿真軟件為Matlab2018b。

比較智能算法J-A參數(shù)辨識主要參考幾方面：磁滯曲線吻合程度、迭代次數(shù)、迭代時間、最小適應(yīng)值以及具體辨識出參數(shù)的相對誤差。

三種方法辨識得到參數(shù)計算出的磁滯曲線與待辨識磁滯曲線的對比如圖2所示。

圖2 各算法辨識參數(shù)生成的磁滯回線與基準磁滯回線的對比Fig.2 Hysteresis loops generated by parameters identified by each algorithm as compared with reference hysteresis loops

由圖2可見，幾種算法所辨識的結(jié)果計算出的曲線基本與實際磁滯曲線相吻合。證明了利用IBES算法對J-A模型參數(shù)進行辨識的有效性。

三種算法的均方根誤差(適應(yīng)度函數(shù))隨迭代次數(shù)的變化如圖3所示。

圖3 各算法均方根誤差隨迭代次數(shù)的變化Fig.3 Change of root mean square error of each algorithm with number of iterations

由圖3可見，SA算法經(jīng)過多次迭代后依然未收斂至全局最優(yōu)值，收斂速度最慢。BES算法快速收斂于全局最優(yōu)值附近，但是從第6代開始收斂速度減慢，在第20代收斂于局部最優(yōu)值；而IBES算法在10代前收斂速度均高于BES算法，收斂至相同結(jié)果附近所用的迭代次數(shù)遠少于BES和SA，最終得到的適應(yīng)度值也優(yōu)于其余兩種算法，三種算法的最終均方根誤差見表1。

表1 IBES與BES、SA算法辨識結(jié)果比較Tab.1 Comparison of identification results between IBES and BES and SA algorithms

由表1可見，根據(jù)仿真結(jié)果，IBES算法的均方根誤差優(yōu)于SA與BES算法，適應(yīng)度函數(shù)越小，說明曲線對應(yīng)的點距離越近，即辨識出參數(shù)所計算的曲線與實際曲線越吻合，而由于使用IBES算法所辨識參數(shù)與理論值較為貼近，所生成的仿真B-H曲線和基準曲線的吻合度較高。

各個算法辨識出的5個參數(shù)見表2，參數(shù)辨識所花時間見表3。

表2 IBES與BES、SA算法辨識結(jié)果比較Tab.2 Comparison of identification results between IBES and BES and SA algorithms

表3 IBES與BES、SA算法辨識耗時比較Tab.3 Comparison of identification time consuming between IBES and BES and SA algorithms

由表2展示的辨識結(jié)果可見，SA算法在每個參數(shù)的辨識上誤差都最大，其中參數(shù)c的相對誤差達到了9%，其余參數(shù)的誤差均在5%以上，表明J-A參數(shù)辨識所使用的傳統(tǒng)智能算法的精度不高；未改進的BES算法優(yōu)于SA，同樣以參數(shù)c為例，相對誤差為2.5%，其余參數(shù)的相對誤差基本為1%～2%左右；IBES算法辨識參數(shù)c的相對誤差為0.19%，參數(shù)a的相對誤差為0.25%，其余參數(shù)的相對誤差也均為0.2%以下，最為精準地辨識出了各種參數(shù)。

由表3可見，SA算法耗時最長，達到了2 912 s，BES算法1 714 s，IBES算法由于需要進行精英反向?qū)W習過程，耗時稍長于BES算法，達到了1 833 s，但是達到相同的適應(yīng)度值所需要的時間更短，并且最終能收斂于更優(yōu)的值，實際應(yīng)用中，通常不需設(shè)置這么高的迭代次數(shù)來進行參數(shù)辨識，所以時間相差會更短，如迭代20次左右，BES和IBES已到達最優(yōu)解附近，時間僅相差16 s，而IBES達到了更佳的目標函數(shù)值，能避免算法陷入局部最優(yōu)，綜合來看，IBES算法具有更高的效率。

在實測中，測量結(jié)果常常受到儀器、人為觀測誤差等方面的影響，與理論曲線有一定的偏差，為了驗證本文在誤差影響下是否仍具有有效性，本文對噪聲曲線進行參數(shù)辨識。

在曲線中加上信噪比為50 dB與70 dB的噪聲，局部圖如圖4所示。再對曲線進行參數(shù)辨識，可得辨識結(jié)果見表4。

圖4 帶噪聲曲線Fig.4 Curves with noise

表4 利用IBES算法辨識含噪聲曲線結(jié)果Tab.4 Using IBES algorithm to identify curve containing noise

由表4可見，加噪聲后辨識出的5個參數(shù)依然與原參數(shù)接近，各參數(shù)誤差基本在1%以內(nèi)，幾乎不受噪聲影響，本文方法可以正確辨識帶噪聲數(shù)據(jù)，這也更加證明了方法的實用性。

6 結(jié)論

智能優(yōu)化算法求取J-A模型參數(shù)是常用的方法，但是傳統(tǒng)的方法如SA、GA、粒子群優(yōu)化算法存在求解速度慢、算法后期種群多樣性和搜索能力均有所下降等缺陷，使得其應(yīng)用起來效率不高。本文使用最新的智能算法——禿鷹搜索算法，并且在此基礎(chǔ)上采用Tent混沌映射，自適應(yīng)權(quán)重與精英反向?qū)W習，形成IBES算法，提高了BES算法的各方面性能，解決了現(xiàn)有智能算法求J-A模型參數(shù)易陷入局部最優(yōu)、求解精度與計算效率較低的問題。為J-A磁滯模型的參數(shù)辨識提供了新方法。得到以下結(jié)論:

(1) IBES算法具有求解精度高、收斂速度快的特點，對求解 J-A 磁滯模型參數(shù)這一工程問題有較好的效果，且不受噪聲影響，具有很強的實用性。

(2) 無論與傳統(tǒng)的J-A模型參數(shù)辨識方法，還是與未改進的BES算法相比，IBES均具有更高的效率與精確度。