鄭慶徽，韓海山

（內(nèi)蒙古民族大學(xué)數(shù)理學(xué)院，內(nèi)蒙古通遼 028043）

特征值互補(bǔ)問(wèn)題（簡(jiǎn)稱EICP）[1］是指：尋找一個(gè)實(shí)數(shù)λ和一個(gè)非零向量x∈Rn{0} ，使得

其中，矩陣A,B∈Rn×n，B為正定矩陣。EICP是由COSTA等[2］于1999年對(duì)機(jī)械系統(tǒng)進(jìn)行靜態(tài)平衡研究而提出的，并在動(dòng)力學(xué)系統(tǒng)等科學(xué)工程領(lǐng)域應(yīng)用廣泛。

COSTA等[3］證明了EICP可以用錐約束特征值問(wèn)題來(lái)表示：對(duì)于給定的矩陣A∈Rn×n和正定矩陣B∈Rn×n，求實(shí)數(shù)λ和向量x∈Rn{0} ，使

其中，Rn+={(x1,x2,…,xn)T∈Rn|xi≥0,i=1,…,n} 稱為Pareto錐。若存在x∈Rn{0 } 使得公式（3）成立

則實(shí)數(shù)λ稱為A的Pareto特征值，非零向量x稱為A的Pareto特征向量，A的所有Pareto特征值構(gòu)成的集合稱為A的Pareto譜[4］。

COSTA等[3］定義了Mn(R)的Pareto勢(shì)

COSTA等[3］和SEEGER等[4］給出了Pareto勢(shì)較為精確的下界。對(duì)于n≥2 時(shí)，有

JUDICE等[5］證明了EICP等價(jià)于單純形Ω={x∈Rn|eTx=1,x≥0}（e是Rn中的單位向量）上的有限維變分不等式VI(F,Ω)。令F:Rn{0 } →Rn是由下式所定義的映射

若矩陣A是嚴(yán)格共正（SC）的，結(jié)合式（1）可將求解EICP 轉(zhuǎn)化為求解變分不等式VI(F,Ω)，即求∈Ω，使得F()(x-)≥0,?x∈Ω。

近年來(lái)，EICP 的求解算法有：ADLY 等[6］的半光滑牛頓法，ABDI 等[7］的逆互補(bǔ)冪迭代法（ICPIM），韓海山等[8］和黃迪帥等[9］的非線性ABS 算法，張美玲等[10］的經(jīng)典遺傳算法（GA）和自適應(yīng)遺傳算法（AGA），趙銳等[11］的基本粒子群優(yōu)化方法（BPSO）和二階粒子群優(yōu)化方法（SECPSO），趙寒等[12］的交替方向乘子法（ADMM）等。

人工蜂群算法（簡(jiǎn)稱ABC）是2005年由KARABOGA[13］提出的一種群智能算法，該算法在全局優(yōu)化方面表現(xiàn)出了較好的性能。由于ABC不需要求出目標(biāo)函數(shù)的導(dǎo)數(shù)，且對(duì)初始點(diǎn)的依賴度小，算法的魯棒性好，所以應(yīng)用范圍廣。因此，ABC是一種具有優(yōu)異性能的群智能優(yōu)化算法，同時(shí)，ABC為求解非線性方程及非線性方程組提供了新思路。

用ABC求解EICP，首先借助于FB函數(shù)把EICP轉(zhuǎn)化為一個(gè)非線性方程組，用適應(yīng)度函數(shù)將非線性方程組轉(zhuǎn)化為一個(gè)約束優(yōu)化問(wèn)題，然后，用ABC對(duì)優(yōu)化問(wèn)題進(jìn)行求解，并對(duì)其收斂性進(jìn)行證明。最后，用數(shù)值實(shí)驗(yàn)對(duì)所建立的EICP-ABC進(jìn)行驗(yàn)證，說(shuō)明其有效性。

1 預(yù)備知識(shí)

記z=(x,y,λ)∈Rn×Rn×R。通過(guò)Fischer-Burmeister（FB）函數(shù)[6］

把（1）轉(zhuǎn)化成非線性方程組：

ABC可以求解非線性方程組，其基本思想是將非線性方程組轉(zhuǎn)換成最優(yōu)化問(wèn)題，通過(guò)解最優(yōu)化問(wèn)題進(jìn)而得到非線性方程組的解。即將非線性方程組（6）轉(zhuǎn)換成最小化問(wèn)題

因此，非線性方程組（6）的解是滿足F(z)=0 的全局最優(yōu)解z*，其中，F(xiàn)為適應(yīng)度函數(shù)[14-15］。同樣，也可將非線性方程組（6）轉(zhuǎn)化成最大化問(wèn)題進(jìn)行求解[10-11］。

在ABC中，蜜源位置代表著最優(yōu)化問(wèn)題的解。蜂群中有被雇用的蜜蜂和沒(méi)有被雇用的蜜蜂，這些被雇用的蜜蜂也被稱作引領(lǐng)蜂，而非雇用的蜜蜂則有跟隨蜂和偵查蜂。引領(lǐng)蜂具有記憶功能，當(dāng)其搜索到一處蜜源時(shí)可將蜜源信息進(jìn)行存儲(chǔ)，并以一定的概率傳遞給跟隨蜂，其中，蜜源數(shù)量等于偵察蜂的數(shù)量。根據(jù)引領(lǐng)蜂所提供的蜜源信息，跟隨蜂采用輪盤賭選擇跟隨那只引領(lǐng)蜂。若經(jīng)過(guò)有限次迭代后還沒(méi)有產(chǎn)生更好的蜜源，即放棄當(dāng)前的蜜源，跟隨蜂轉(zhuǎn)變成偵察蜂，并且會(huì)隨機(jī)地尋找新蜜源。

定義1[16］對(duì)于一個(gè)蜂群，其中，全部蜜蜂的狀態(tài)組成了此蜂群的狀態(tài)，記為s=(Y1,Y2,…,YN)。蜂群中所有可能狀態(tài)集構(gòu)成蜂群的狀態(tài)空間，用S={s=(Y1,Y2,…,YN) |Y1∈Y,1 ≤i≤SN} 表示。

定義2[16］在ABC的迭代中，對(duì)于?si,sj∈S，蜂群狀態(tài)由si經(jīng)一次迭代轉(zhuǎn)移到sj，記Hs(si)=sj。由si經(jīng)一次迭代將蜂群狀態(tài)移到sj的轉(zhuǎn)移概率為：

對(duì)于優(yōu)化問(wèn)題（7），有隨機(jī)優(yōu)化算法T，第k次迭代的結(jié)果為zk，則下一次迭代的結(jié)果為zk+1=T(zk,ξ)。其中，D為可行解的集合，F(xiàn)為適應(yīng)度函數(shù)，ξ為算法T迭代中曾經(jīng)搜索到的解[17］。

引理1[17］設(shè)F是可測(cè)的，可測(cè)空間D是Rn上的可測(cè)度的子集是算法T產(chǎn)生的解序列，若算法T滿足下述2個(gè)條件：

條件1F(T(z,ξ))≤F(z)，若ξ∈D，則有F(T(z,ξ))≤F(ξ)。

條件2對(duì)于D的任意Borel子集B，使得v[B] >0，則有其中，uk[B] 是集合B中算法T第k次迭代搜索解的概率測(cè)度，則有

其中，p(zk∈Rε,M)是在Rε,M中算法T第K步的解zk的概率測(cè)度。

定義3[16］設(shè)優(yōu)化問(wèn)題（10）的最優(yōu)解是g?，定義人工蜂群的最優(yōu)狀態(tài)集

若G=S，在求解空間內(nèi)，各個(gè)解既是可行的又是最優(yōu)的，這時(shí)再做最優(yōu)的就沒(méi)有意義了。因此，以下是關(guān)于G?S的情況的討論。

引理2[16］假定馬爾科夫鏈有一個(gè)非空閉集E，且不存在另一個(gè)非空閉集O，使得E?O=?，則當(dāng)j∈E時(shí)，有l(wèi)ni→m∞p(Yn=j)=πj，且當(dāng)j?E時(shí)，有l(wèi)ni→m∞p(Yn=j)=0。

引理3[16］人工蜂群狀態(tài)空間S，不存在非空閉集M，使得M?G=?。

2 求解EICP的人工蜂群算法（EICP-ABC）及收斂性

考慮非線性方程組

設(shè)隨機(jī)產(chǎn)生向量xk∈Rn,k=1,2,…,N，令

并記zk=(xk,yk,λk)。根據(jù)ABC求解非線性方程組的方法，定義適合度函數(shù)為[14］：

即將式（7）的非線性方程組轉(zhuǎn)換成帶有約束的最優(yōu)化問(wèn)題：

其中，D={(x1,x2,…xn)T|ai≤xi≤bi,i=1,2,…,n} ，且式（9）與式（10）是等價(jià)的[15］。

2.1 算法描述

EICP-ABC算法的具體步驟如下：

Step0：輸入控制參數(shù)。蜜源個(gè)數(shù)、引領(lǐng)蜂個(gè)數(shù)與跟隨蜂個(gè)數(shù)是相等的，并用SN表示，迭代次數(shù)為M，蜜源試驗(yàn)次數(shù)上限為L(zhǎng)，蜜源變換的加速度系數(shù)為a，收斂精度為ε，蜜源的位置范圍為[zmin,zmax]。

Step1：初始化。初始化蜜源的位置隨機(jī)產(chǎn)生向量xij∈Rn,i=1,2,…SN，j=1,2,…,D，令

記zij=(xij,yij,λij)，并計(jì)算其適應(yīng)度值。

Step2：引領(lǐng)蜂依式（11）對(duì)鄰近區(qū)域進(jìn)行搜索，得到新的蜜源Vi，并求出其適應(yīng)度值

如果Vi的適應(yīng)度值優(yōu)于zi，以新蜜源Vi取代舊蜜源zi，將新蜜源Vi作為當(dāng)前最好解，反之保留原zi不變。

Step3：計(jì)算所有zi的適應(yīng)度值，并使用輪盤賭計(jì)算與zi相關(guān)的概率值Pi

Step4：跟隨蜂依概率值Pi選擇跟隨哪只引領(lǐng)蜂出去采蜜并依式（11）對(duì)鄰近區(qū)域進(jìn)行搜索產(chǎn)生新蜜源Vi，計(jì)算其適應(yīng)度值。如果Vi的適應(yīng)度值優(yōu)于zi，以新蜜源Vi取代舊蜜源zi，將新的蜜源Vi作為當(dāng)前最好解，反之保持zi不變。

Step5：判斷是否要放棄蜜源，也就是說(shuō)，在經(jīng)過(guò)L次循環(huán)后蜜源沒(méi)有得到改善，則偵察蜂將基于下式產(chǎn)生一個(gè)新蜜源zi來(lái)替換舊蜜源，即

Step6：記錄到目前為止最好的解。

Step7：檢驗(yàn)是否滿足終止條件。若|F(zi) |<ε，算法終止輸出最優(yōu)解zi，否則轉(zhuǎn)至Step3。

2.2 EICP-ABC的收斂性證明

由于在EICP-ABC 中蜜蜂的運(yùn)動(dòng)過(guò)程本質(zhì)上是一個(gè)隨機(jī)過(guò)程，因此，可以利用馬爾科夫理論及蜂群的一步轉(zhuǎn)移概率，同時(shí)，結(jié)合隨機(jī)優(yōu)化算法的收斂性準(zhǔn)則，對(duì)EICP-ABC的收斂性進(jìn)行證明。

引理4在EICP-ABC中蜂群狀態(tài)序列{s(t):t>0} 為有限齊次的馬爾科夫鏈（MC）。

證明證明與文獻(xiàn)[17］中的定理2相類似。

引理5對(duì)于EICP-ABC，在其狀態(tài)序列{s(t);t≥0} 中，最優(yōu)人工蜂群的最優(yōu)狀態(tài)集G是狀態(tài)空間S上的一個(gè)閉集。

證明設(shè)si∈G,?sj?G，對(duì)于任意的轉(zhuǎn)移步長(zhǎng)l,l≥1，由查普曼-科莫高洛夫方程（CK方程）可得

其中，Plsi,sj表示經(jīng)過(guò)l步蜂群狀態(tài)由si轉(zhuǎn)移到狀態(tài)sj的概率，在式（12）中有P(Hs(src-1)=src)，使得src-1∈G且src?G，其中，1 ≤c≤l，由定義2 得蜂群的轉(zhuǎn)移概率為由src-1∈G且src?G，有F(Yc)>F(Yc-1)=F(g*)=inf(F(b))，b∈B，則至少存在P(Hs(src-1)=src)=0，使得=0，所以G是狀態(tài)空間S上的一個(gè)閉集。

引理6當(dāng)EICP-ABC中的蜂群內(nèi)部迭代趨于無(wú)窮時(shí)，其狀態(tài)序列必然為最優(yōu)狀態(tài)集G。

證明由引理2、引理3、引理5即可得出引理6成立。

引理7EICP-ABC同時(shí)滿足條件1與條件2。

證明EICP-ABC在每次迭代中都使用貪婪算法，即

因此，EICP-ABC在每次迭代中保留最好的蜂群位置，所以條件1成立。若滿足條件2，則要求算法T經(jīng)過(guò)無(wú)限次的連續(xù)搜尋未找到B中的元素的概率為0。由引理6可知，在無(wú)限次的連續(xù)搜索中，算法T不能找到全局最優(yōu)解的概率為0，有0

定理1EICP-ABC是依概率1收斂的。

證明由于EICP-ABC同時(shí)滿足條件1與條件2，通過(guò)引理1可得出EICP-ABC收斂到全局最優(yōu)。因此，EICP-ABC是依概率1全局收斂的。

3 數(shù)值實(shí)驗(yàn)

數(shù)值實(shí)驗(yàn)是在處理器Intel（R）Core（TM）i5-7200U CPU@2.50 GHz，內(nèi)存8.0 GB等Windows 10環(huán)境下用MATLAB軟件運(yùn)行的，數(shù)值實(shí)驗(yàn)中，取B=In。

例1求系數(shù)矩陣為

的EICP。此矩陣已被ADLY等[6］證明存在3個(gè)Pareto特征值。

表1中的所有結(jié)果都是在控制參數(shù)為SN=20，M=200，L=36，a=1時(shí)求得的，文中所構(gòu)造的人工蜂群算法求得了此EICP所有的Pareto特征值。

表1 用EICP-ABC求解關(guān)于系數(shù)矩陣的EICPTab. 1 Solving EICP of coefficient matrix by EICP-ABC

例2求系數(shù)矩陣為

的EICP。此矩陣已被ADLY等[6］證明存在9個(gè)Pareto特征值。

表2中的所有結(jié)果都是在控制參數(shù)為SN=20，M=400，L=108，a=1 時(shí)求得的，文中所構(gòu)造的人工蜂群算法求得了此EICP所有的Pareto特征值。