微生物動力學(xué)優(yōu)化算法*

陸秋琴，黃光球

西安建筑科技大學(xué) 管理學(xué)院，西安 710055

1 引言

在工程中存在很多具有大量局部最優(yōu)解的優(yōu)化問題，而且在很多情況下這些優(yōu)化問題的目標函數(shù)和約束條件的數(shù)學(xué)表達式?jīng)]有限制條件，傳統(tǒng)的優(yōu)化算法無法求解此類優(yōu)化問題。目前，求解此類優(yōu)化問題的方法是群智能優(yōu)化算法[1]。與傳統(tǒng)優(yōu)化算法不同，群智能優(yōu)化算法求解一個優(yōu)化問題時，采用若干試探解同時進行演化計算，這種“群起而攻之”的方法可求解一些非常困難的優(yōu)化問題[2]。然而，不存在一種群智能優(yōu)化算法可求解所有類型的優(yōu)化問題[3]。在群智能優(yōu)化算法中，每個試探解被比喻成具有生物特征的個體，一些特殊生物進化場景被轉(zhuǎn)換成群智能算法的邏輯結(jié)構(gòu)，活躍在該場景中的生物個體之間的相互作用關(guān)系被用來構(gòu)造生物個體的演化算子[4]。

種群動力學(xué)[5]是一門反映生物種群相互作用關(guān)系的數(shù)學(xué)理論，種群動力學(xué)包含有很多分支，每個分支均有改造成一種群智能算法的潛能。文獻[6]基于Lotka-Volterra 模型構(gòu)造出了一種特殊種群動力學(xué)優(yōu)化算法，該算法利用3個種群間的相互作用關(guān)系即互惠共存關(guān)系、競爭關(guān)系、捕食-被食關(guān)系構(gòu)建出進化算子，所有算子的演化特性可以確保算法具有全局收斂性。文獻[7]采用種群動力學(xué)理論構(gòu)造出了種群動力學(xué)優(yōu)化（population dynamics-based optimization，PDO）算法。該算法將任意兩種群間的捕食-被食、互利、融合、競爭、突變和選擇等行為用于構(gòu)造種群的進化策略。文獻[8]基于種群動力學(xué)理論提出了一種改進的PDO算法。該算法被用來快速確定井下多點集中涌出礦塵就地消除最佳引導(dǎo)路徑。文獻[9]采用種群動力學(xué)中的捕食-被食理論提出了另一種改進的PDO算法，該算法可快速確定井下礦塵收集站、井下風(fēng)機以及通風(fēng)構(gòu)筑物的最佳安裝位置。

微生物培養(yǎng)是利用培養(yǎng)室培養(yǎng)一種或多種微生物種群，微生物培養(yǎng)動力學(xué)是研究培養(yǎng)室中在人類控制下微生物種群相互作用、不斷演化的數(shù)學(xué)理論。微生物培養(yǎng)動力學(xué)是種群動力學(xué)分支之一[5]。為了利用微生物培養(yǎng)動力學(xué)的一些特殊性質(zhì)來構(gòu)造出高性能群智能優(yōu)化算法，本文采用了一種帶時滯影響的混雜食物鏈微生物脈沖培養(yǎng)動力學(xué)模型來構(gòu)造微生物動力學(xué)優(yōu)化算法（microbial dynamics optimization，MDO）。在該算法中，本文采用了與現(xiàn)有種群動力學(xué)優(yōu)化算法不同的設(shè)計思路，提出了將微生物培養(yǎng)動力學(xué)模型轉(zhuǎn)化為能求解某些優(yōu)化問題的一般方法；構(gòu)造出的算子可以充分反映培養(yǎng)系統(tǒng)中的有毒物質(zhì)和培養(yǎng)液對微生物種群的影響以及微生物種群之間的相互作用關(guān)系，從而體現(xiàn)微生物培養(yǎng)動力學(xué)理論的基本思想。MDO 算法能夠全局收斂，對求解維數(shù)較高的優(yōu)化問題具有一定優(yōu)勢。

2 算法基本原理

2.1 MDO算法的場景設(shè)計

假設(shè)在微生物培養(yǎng)系統(tǒng)E中生活有N個微生物種群，即P1,P2,…,Pi,…,PN；每個種群不但以E中的培養(yǎng)液為食，而且還以E中的其他K個種群為食。在微生物培養(yǎng)過程中，周期性地向E中注入已調(diào)制好的新鮮培養(yǎng)液；與此同時，從E中不斷移除無用的培養(yǎng)液；微生物種群吸收的培養(yǎng)液不會立即就轉(zhuǎn)化成新的微生物種群，或者說，在從培養(yǎng)液向微生物種群轉(zhuǎn)化這一過程存在一個時間段；隨著培養(yǎng)液的注入，有害物質(zhì)也隨之注入。為了調(diào)控有害物質(zhì)對種群的破壞性影響，必須不斷評估種群受到有害物質(zhì)的危害程度。培養(yǎng)系統(tǒng)的運行可通過對兩個因素進行不同控制來實現(xiàn)：一是培養(yǎng)液中營養(yǎng)物質(zhì)的濃度S；二是培養(yǎng)液的流量Q。濃度S的變化會引起微生物種群的生長產(chǎn)生變化。

微生物種群Pi的特征可用Pi={fi,1,fi,2,…,fi,n}表示，其中fi,j是Pi的特征j，n是每個種群的特征總數(shù)，i=1～N,j=1～n；E中培養(yǎng)液的營養(yǎng)物質(zhì)濃度對種群Pi的影響體現(xiàn)在對種群Pi的特征的隨機影響上。

假設(shè)所要求解的優(yōu)化模型為：

式中，H為搜索空間；f(X)為目標函數(shù)；X為n維解向量，X=(x1,x2,…,xi,…,xn)，xi為解分量。在H中任意選擇N個初始解，即H={X1,X2,…,XN}，其中Xi=(xi,1,xi,2,…,xi,n),i=1～N；培養(yǎng)系統(tǒng)E與搜索空間H相對應(yīng)，優(yōu)化模型的N個初始解與E中的N個微生物種群一一對應(yīng)，即Xi與Pi一一對應(yīng)；更精確地說，Xi的變量xi,j與Pi的特征fi,j一一對應(yīng)。

綜上，初始解與微生物種群是等價概念，以后不再區(qū)分。Pi攝入營養(yǎng)后，它的生長狀態(tài)會產(chǎn)生變化，將該變化投射到H，等價于初始解Xi從空間位置a轉(zhuǎn)移到b。一個空間位置稱作一個狀態(tài)，用它的下標描述。

假設(shè)Pi當前狀態(tài)為c，此等價于在H中的位置為Xc。若Pi攝入營養(yǎng)后，從狀態(tài)c演變到狀態(tài)d，此等價于在H中從位置Xc跳轉(zhuǎn)到位置Xd。按式（1）計算，若f(Xc)>f(Xd)，表明位置Xc比位置Xd要優(yōu)，則認為Pi強壯。反之，若f(Xc)≤f(Xd)，表明位置Xd比位置Xc差，或沒有區(qū)別（因f(Xc)=f(Xd)），則認為Pi虛弱。強壯的種群能以較高的幾率不斷生長,而虛弱的種群則有可能不再生長。

Pi的強弱可用微生物生長指數(shù)（microbial growth index，MGI）來表示。質(zhì)量好的解對應(yīng)具有較高MGI指數(shù)的種群，即強壯的種群，質(zhì)量差的解對應(yīng)具有較低MGI指數(shù)的種群，即虛弱的種群。對于式（1），Pi的MGI指數(shù)按下式計算：

在E中，微生物培養(yǎng)的生物智能特征體現(xiàn)在如下幾個方面：

（1）各種群因爭奪培養(yǎng)液而存在相互影響，此類影響會在種群的某些特征間的隨機相互作用上體現(xiàn)出來，此為“競爭的智能特征”。

（2）一個種群以另外一些種群為食，那些被食的種群的某些特征就會融入到該種群的某些特征中，此為“吃的智能特征”。

（3）任一種群的進化是通過不斷向其他強勢種群學(xué)習(xí)而取得的，此為“學(xué)習(xí)的智能特征”。

微生物培養(yǎng)的生物智能特征最終以生物演化算子的形式體現(xiàn)出來。所有微生物種群間的相互影響投射到式（1）的H中，就是解向量間存在相互作用。

MDO 算法就是采用此類策略來實現(xiàn)對式（1）全局最優(yōu)解的求解。

2.2 微生物培養(yǎng)動力學(xué)模型

假設(shè)在培養(yǎng)系統(tǒng)E中共有N種微生物種群，每個種群既以培養(yǎng)液為食，又以E中L個其他種群為食。此外，有毒物質(zhì)會毒害所有種群。時期t，種群Pi的濃度為yi(t),yi(t)≥0,i=1～N；培養(yǎng)液以流量Q流入到E中，其中營養(yǎng)物質(zhì)的濃度為S0。該培養(yǎng)系統(tǒng)具有時滯影響的混雜食物鏈微生物脈沖培養(yǎng)動力學(xué)模型[5]如下：

式中，S(t)表示時期t培養(yǎng)液的營養(yǎng)物質(zhì)濃度；μ0表示種群的相對增長系數(shù)；C0表示種群對培養(yǎng)液的消耗率；MS(t)表示時期t種群Pi捕獲到的K個種群的編號的集合；K0、α0、γ0、W0表示與種群生長相關(guān)的常數(shù)；yi(t)表示時期t種群Pi的濃度；c(t)表示有害物質(zhì)濃度，c(t)≥0；q0表示每個培養(yǎng)液投放時期輸入的有害物質(zhì)的數(shù)量；r表示種群因有害物質(zhì)而毒害死亡的減少率，r>0；T表示培養(yǎng)液的投放周期；τ表示培養(yǎng)液向微生物種群轉(zhuǎn)化的滯后時間；t+為培養(yǎng)液投放時間點；k表示正整數(shù)，k≥1。

在時期t，種群Pi的占比ri(t)按下式計算：

種群Pi的占比ri(t)越大，被其他種群捕食的幾率也會越大。記時期t參數(shù)C0、Q、K0、S0、α0、μ0、r、γ0、W0、q0的取值分別為為方便計算，將式（3）改為離散遞推形式，即：

若t≠kT，則：

若t=kT，則：

式（5）和式（6）中各參數(shù)的含義及其取值限制條件參見表1。

Table 1 Strategy of taking value for parameters of microbial dynamics model表1 具有微生物培養(yǎng)動力學(xué)模型參數(shù)的取值策略

表1中，Rand(A,B)表示在區(qū)間[A,B]生成一個服從均勻分布的隨機數(shù)，A和B為常數(shù)，A≤B。

2.3 特征微生物種群集合的生成策略

時期t，設(shè)當前種群為Pi，特征微生物種群集合的生成策略如下：

（1）產(chǎn)生被食種群集合MS(t)。時期t，種群Pi以另外K個種群為食，即從N-1 個種群中分別以{r1(t),r2(t),…,rN(t)}為概率分布隨機選擇K個種群（但Pi不能被選中），這些種群形成集合MS(t)，不妨設(shè)MS(t)={i1,i2,…,ik,…,iK}，ik為種群的編號，k=1～K。

（2）分別產(chǎn)生30%、20%、10%強壯種群的集合FMt、GMt、HMt。時期t，先從N個種群中隨機選出L個種群，這些種群的MGI指數(shù)要比當前種群Pi高30%，形成集合，其中s1,s2,…,sL是這些被選中種群的編號。然后從剩余的N-L個種群中隨機選出L個種群，這些種群的MGI指數(shù)比當前種群Pi高20%，形成集合，其中g(shù)1,g2,…,gL是這些被選中種群的編號；最后從剩余的N-2L個種群中隨機選出L個種群，這些種群的MGI指數(shù)比當前種群Pi高10%，形成集合，其中h1,h2,…,hL是這些被選中種群的編號。

（3）分別產(chǎn)生15%的強壯種群、虛弱種群、普通種群的集合ASt、BSt、CSt。時期t，從N個種群中隨機選出L個強壯種群，其MGI指數(shù)要高于當前種群Pi的15%，形成強壯種群集合其中a1,a2,…,aL是這些被選中種群的編號；然后從N個種群中隨機選出L個虛弱種群，其MGI指數(shù)要低于當前種群Pi的15%，形成虛弱種群集合，其中b1,b2,…,bL是這些被選中種群的編號；最后從剩余的N-2L個種群中隨機選L個普通種群，其MGI指數(shù)與當前種群Pi的MGI指數(shù)沒有關(guān)系，形成普通種群集合其中c1,c2,…,cL是這些被選中種群的編號。

2.4 微生物種群演化算子設(shè)計

MDO 算法利用微生物培養(yǎng)動力學(xué)模型式（3）來構(gòu)造演化算子，從而實現(xiàn)培養(yǎng)系統(tǒng)與種群之間以及種群與種群之間的信息交換，進而實現(xiàn)對式（1）的搜索空間H的搜索。

（1）吸收算子。吸收算子描述的是種群Pi與另外K個種群之間的營養(yǎng)供給關(guān)系。將集合MS(t)中被食種群的一個隨機選出的特征fs,j,s∈MS(t)傳給當前種群Pi，使種群Pi獲得生長：

式中，xs,j(t)表示時期t種群Ps的特征j的狀態(tài)值，即在時期t變量xs,j的取值；vi,j(t+1)表示時期t+1 種群Pi的特征j的狀態(tài)值，即在時期t+1變量xi,j的取值；取βs=Rand(-0.95,0.95),ηi=Rand(0.45,0.65)。

（2）攫取算子。攫取算子描述的是一個種群攫取其他種群的優(yōu)勢特征，從而使其變得強壯。讓FMt、GMt、HMt中的3L個強壯種群的特征fsl,j，sl∈{h1,h2,…,hL;g1,g2,…,gL;s1,s2,…,sL}所對應(yīng)的優(yōu)勢狀態(tài)值傳給種群Pi的特征fi,j，使Pi變得強壯。假設(shè)當前種群為Pi進行優(yōu)勢攫取的特征為fi,j，則有：

式中，d1、d2、d3，e1、e2、e3，p1、p2、p3分別是從集合FMt、GMt、HMt中隨機選出來的，且滿足d1≠d2≠d3，e1≠e2≠e3，p1≠p2≠p3。

（3）混雜算子。混雜算子描述的是種群與種群之間的相互影響關(guān)系。對于當前種群為Pi，讓ASt、BSt、CSt中的3L個種群的特征b1,b2,…,bL;c1,c2,…,cL}所對應(yīng)的狀態(tài)值經(jīng)混雜融合后傳給種群Pi的特征fi,j，則有：

式中，ρs=Rand(0.75,0.85),λi=Rand(0.65,0.85)。

（4）毒素算子。毒素算子描述的是外部系統(tǒng)中的有毒物質(zhì)對微生物的毒害作用，該毒害作用能夠在微生物種群之間傳遞。對于當前種群Pi來說，有：

式中，g1、g2、g3是從{1,2,…,N}中隨機選出來的，且滿足g1≠g2≠g3。

（5）生長算子。生長算子的定義如下：

式中，Xi(t)=(xi,1(t),xi,2(t),…,xi,n(t))；Vi(t+1)=(vi,1(t+1),vi,2(t+1),…,vi,n(t+1))；利用式（2）計算MGI(Vi(t+1))和MGI(Xi(t))。

2.5 MDO算法

（S1）初始化：①令t=0，演化時期數(shù)G=8 000～60 000，誤差要求ε=10-5～10-10，N=50～500，微生物種群被捕食的最高概率E0=1/1 000～1/100，K=2～10，T=1～10，τ=1～10，L=2～10。②隨機確定{X1(0),X2(0),…,XN(0)}，{y1(0),y2(0),…,yN(0)},S(0)。

（S2）執(zhí)行下列操作：

Fort=0 toG

按式（4）計算r1(t),r2(t),…,rN(t)。

若T不能整除t，則按式（5）計算c(t+1)、S(t+1)；若T能整除t，則按式（7）計算c(t+1)、S(t+1)。

2.6 算法特征分析

2.6.1 時間復(fù)雜度分析

MDO算法的時間復(fù)雜度計算過程如表2所示。

Table 2 Computation table of time complexity in MDO表2 MDO算法的時間復(fù)雜度計算表

2.6.2 MDO算法的收斂性分析

（1）Markov特性。從MDO算法的各算子的定義式（9）～式（12）知，時期t+1 任一新解的計算僅涉及該解在時期t時的狀態(tài)，而與該解在時期t以前是如何演變到時期t時的狀態(tài)的歷程無關(guān)，表明MDO 算法的演變過程具有Markov特性。

（2）“步步不差”特性。從生長算子的定義式（13）知，時期t+1 任一種群的MGI指數(shù)永遠不會低于其在時期t時的MGI指數(shù)，即MDO 算法的演變過程具有“步步不差”特性。

依據(jù)文獻[10]可知，滿足上述兩個特性的MDO算法具有全局收斂性，其相關(guān)證明可參見文獻[10]，本文不再贅述。

3 MDO算法與其他算法比較

CEC2013[11]是國際上通用智能優(yōu)化算法測試包，其中包括有28 個經(jīng)過精心設(shè)計的基準函數(shù)優(yōu)化問題。本文選用其中的15個基準函數(shù)優(yōu)化問題來測試MDO算法的性能，如表3所示。

表3 中，O是n維向量，其值可隨機產(chǎn)生。本文用MDO算法去求解表3中的15個函數(shù)優(yōu)化問題，其參數(shù)是N=200，n=50，G=107，ε=10-7，E0=1/200，K=3，T=4，τ=3，L=3。與MDO算法進行比較的7種智能優(yōu)化算法均選自國際著名期刊近期刊登的知名算法，這些算法如表4 所示，即RCGA（real-coded genetic algorithm）[12]、DASA（differential ant-stigmergy algorithm）[13]、NP-PSO（non-parametric particle swarm optimization）[14]、MpBBO（metropolis biogeographybased optimization）[15]、Copt-aiNet（artificial immune networks for combination optimization）[16]、SLADE（symmetric Latin-based adaptive differential evolution）[16]和GRABC（grab artificial bee colony algorithm）[17]。這7種算法的終止運行條件為：進化代數(shù)G=107或者最優(yōu)解誤差ε=10-7。

針對每個基準函數(shù)優(yōu)化問題，上述8種算法分別獨立運行51 次，以尋找全局最優(yōu)目標函數(shù)。表5～表7 給出了每種算法所求得的最優(yōu)目標函數(shù)值的平均值偏差、中值偏差、標準差、適應(yīng)度評價次數(shù)。8種算法求解每個基準函數(shù)的性能排名是這些算法基于最優(yōu)目標函數(shù)值的平均值偏差和適應(yīng)度評價次數(shù)進行的排名；最終排名是基于8 種算法求解15 個基準函數(shù)所得的排名總積分而進行的排名。

非參數(shù)弗里德曼檢驗[19-20]是基于MDO算法所得的最佳結(jié)果與7 種被比較算法所獲得的最佳結(jié)果之間進行的非參數(shù)檢驗，以確定MDO 算法產(chǎn)生的結(jié)果是否與7 種被比較算法所獲得的結(jié)果有統(tǒng)計上的不同。弗里德曼檢驗的結(jié)果顯示在表8 中，其中Significance=1 表示MDO 算法的性能與被比較算法具有99%的統(tǒng)計學(xué)差異，Significance=0 表示MDO算法的性能與被比較算法的性能沒有顯著差異。在表8中，Significance=1的案例數(shù)和Significance=0 的案例數(shù)分別表示MDO算法與7種被比較算法顯著不同和幾乎相同的求解基準函數(shù)優(yōu)化問題的數(shù)目。

Table 3 15 benchmark function optimization problems表3 15個基準函數(shù)優(yōu)化問題

Table 4 Parameters setting of compared intelligent optimization algorithms表4 參與比較的智能優(yōu)化算法的參數(shù)設(shè)置

從表5～表7 可以看出MDO、RCGA、DASA、NPPSO、MpBBO、Copt-aiNet、SLADE和GRABC按平均最優(yōu)目標函數(shù)值精度的排序如下：

MDO>SLADE>Copt-aiNet>NP-PSO>

DASA>RCGA=GRABC>MpBBO

從表8 可以知道，MDO 算法求解15 個基準函數(shù)優(yōu)化問題的顯著性案例總數(shù)為88，明顯大于不顯著性案例總數(shù)17，表明MDO 算法的性能明顯優(yōu)于7 種被比較的算法。

圖1（a）～（f）說明MDO、RCGA、DASA、NP-PSO、MpBBO、Copt-aiNet、SLADE、GRABC 算法求解表3所示的6 個典型基準函數(shù)優(yōu)化問題F3、F5、F14、F19、F26、F28 時的樣本收斂曲線，圖中水平和垂直軸采用對數(shù)刻度。從表5～表7可以看出，當MDO算法求解這6個優(yōu)化問題時，均能以很高的精度發(fā)現(xiàn)全局最優(yōu)解。綜合看來，MDO 算法的綜合性能要優(yōu)于7種被比較算法，表明其求解精度高且計算速度快。

4 MDO 算法的局部尋優(yōu)能力和全局尋優(yōu)能力分析

4.1 穿透能力和擴展能力的測量方法

穿透能力和擴展能力用于描述MDO 算法在求解過程中陷入局部最優(yōu)解陷阱后從陷阱逃逸并以一定精度收斂到全局最優(yōu)解的能力。穿透能力用于描述MDO 算法在求解優(yōu)化問題時或者快速從陷阱逃逸，或者快速收斂到全局最優(yōu)解的能力；而擴展能力用于描述MDO 算法在求解優(yōu)化問題時或者試探性跳出局部最優(yōu)解陷阱，或者試探性提升最優(yōu)解精度的能力。MDO算法的穿透和擴展能力可以通過微生物種群的行為來描述。當一個種群進行搜索時，其狀態(tài)可以分為穿透狀態(tài)和擴展狀態(tài)。一個種群每次只能停留在一個狀態(tài)：要么停留在穿透狀態(tài)，要么停留在擴展狀態(tài)。為了強調(diào)種群進入局部最優(yōu)解陷阱的行為，將擴展狀態(tài)再分為兩個狀態(tài)，一個稱為膨脹狀態(tài)，另一個稱為粘滯狀態(tài)。

Table 5 Optimal solutions of all algorithms when solving unimodal benchmark functions表5 所有算法求解單峰基準函數(shù)時所得的最優(yōu)解

當種群進行搜索時，如果種群的當前位置長期保持不變，認為種群已經(jīng)陷入粘滯狀態(tài)。兩種情況可能會導(dǎo)致粘滯狀態(tài)：一種是種群陷入局部最優(yōu)解陷阱；另一種情況是種群的當前位置不能被其他種群拋出的信息更新。

當所有種群陷入粘滯狀態(tài)時，搜索將自動停止，全局最優(yōu)解將永遠保持不變。另一方面，如果一些優(yōu)化問題具有很高的條件數(shù)，那么當種群容易陷入粘滯狀態(tài)時，提高全局最優(yōu)解的精度總是困難的。

當MDO算法求解優(yōu)化問題時，若穿透狀態(tài)總是占優(yōu)，則意味著搜索或者快速從陷阱逃逸，或者快速接近全局最優(yōu)解，因此更高的穿透能力對于MDO算法是一件好事；若膨脹狀態(tài)總是占優(yōu)，則意味著搜索或者正在試探性跳出局部最優(yōu)解陷阱，或者正在試探性提升當前最優(yōu)解的精度，但可能耗費巨大的計算時間，因此算法的膨脹能力應(yīng)該受到一定程度的控制；若粘滯狀態(tài)占優(yōu)，則意味著搜索幾乎停止，因此必須完全避免粘滯狀態(tài)。

穿透和擴展能力的評估可以基于四種收斂曲線：（1）任意種群的收斂曲線，跟蹤任意一個種群搜索最優(yōu)解的行為；（2）最佳種群的收斂曲線，其跟蹤最佳種群搜索最優(yōu)解的行為；（3）所有種群的平均收斂曲線，其跟蹤所有種群在搜索最優(yōu)解時的平均行為；（4）最佳收斂曲線，其跟蹤所有種群在搜索過程中發(fā)現(xiàn)當前最優(yōu)解的行為。

任意種群的收斂曲線可能會造成巨大的計算負擔(dān)和計算機內(nèi)存，因為在生態(tài)系統(tǒng)中有許多種群；此外，種群的行為太隨機以至于無法控制它們。最佳種群的收斂曲線難以掌握，因為最佳種群可能隨時間而變化。所有種群的平均收斂曲線與任意種群的收斂曲線相似，不僅可能造成巨大的計算負擔(dān)和計算機內(nèi)存，而且可以提高不良種群的MGI值，并壓縮優(yōu)良種群的MGI值。最佳收斂曲線可以在四條收斂曲線之間產(chǎn)生最小的計算負擔(dān)和計算機內(nèi)存，而且曲線不是關(guān)注任何單個種群，而是關(guān)注由所有種群產(chǎn)生的最佳行為。

Table 6 Optimal solutions of all algorithms when solving multimodal benchmark functions表6 所有算法求解多峰基準函數(shù)時所得的最優(yōu)解

MDO算法的穿透和擴展能力評價是基于最優(yōu)收斂曲線的，是一種簡單、直觀的方法來判斷MDO 算法的穿透和擴展能力。

假設(shè)圖2是MDO算法求解某個最優(yōu)化問題時的最佳收斂曲線，C點是演化的起始位置，對應(yīng)的最佳目標函數(shù)值為F0，CPU 時間t=0；當演化到達時間t=t1時，當前最佳位置處于A，對應(yīng)于最佳目標函數(shù)值是F1；當進化到達時間t=t2時，當前最佳位置處于B，對應(yīng)的最佳目標函數(shù)值是F2。種群的穿透和擴展能力可以通過直線AB的斜率來識別，因此有：

Table 7 Optimal solutions of all algorithms when solving composite benchmark functions表7 所有算法求解復(fù)合基準函數(shù)時所得的最優(yōu)解

顯然，若直線AB越陡，則穿透能力越高；若直線AB從陡峭向緩傾斜演變，則膨脹能力從好變壞，特別是當slope(t1)接近0時，膨脹能力變差，可能出現(xiàn)粘滯狀態(tài)。因此，必須定義兩個門限值λPs和λEs來說明這些情況：使用λPs來區(qū)分穿透狀態(tài)（penetration state）和膨脹狀態(tài)（expansion state），使用λEs來區(qū)分膨脹狀態(tài)和粘滯狀態(tài)（sticky state）。

（1）若slope(t1)>λPs，則意味著某些種群停留在穿透狀態(tài)。

（2）若λEs

（3）若slope(t1)≤λEs，則意味著大多數(shù)種群有進入粘滯狀態(tài)的趨勢，演化可能停止。

現(xiàn)在考慮種群在時間區(qū)間[tA,tB]中的穿透和擴展能力。讓CountPs(tA,tB)表示在區(qū)間[tA,tB],tA≤t≤tB滿足slope(t)>λPs的位置數(shù)；CountEs(tA,tB)表示在區(qū)間[tA,tB]滿足λEs

CountPs(tA,tB)=Count(i|slope(t)>λPs,tA≤t

CountEs(tA,tB)=Count(i|λEs≤slope(t)≤λPs,tA≤t

CountSs(tA,tB)=Count(i|slope(t)<λEs,tA≤t

其中，Count()是計算滿足給定條件的位置的函數(shù)。

定義在時間間隔[tA,tB]內(nèi)的穿透與膨脹協(xié)調(diào)度的測量方法如下：

Table 8 Comparison of Friedman test results(α=0.01)表8 Friedman檢驗結(jié)果比較(α=0.01)

其中，CPs/Es(tA,tB)稱為穿透-膨脹協(xié)調(diào)系數(shù)。

于是，有：

（1）如CPs/Es(tA,tB)>CPs，則意味著搜索在時間間隔[tA,tB]內(nèi)處于穿透占優(yōu)狀態(tài)。

（2）若CPs/Es(tA,tB)≤CPs，則意味著搜索在時間間隔[tA,tB]內(nèi)處于膨脹占優(yōu)狀態(tài)。

類似地，可以定義在時間間隔[tA,tB]內(nèi)的膨脹和粘滯之間的協(xié)調(diào)度的測量方法如下：

這里CEs/Ss(tA,tB)稱為膨脹-粘滯協(xié)調(diào)系數(shù)。

于是，有：

（1）如果CEs/Ss(tA,tB)>CEs，則意味著搜索在時間間隔[tA,tB]內(nèi)處于膨脹占優(yōu)狀態(tài)。

（2）如果CEs/Ss(tA,tB)≤CEs，則意味著搜索在時間間隔[tA,tB]內(nèi)處于粘滯占優(yōu)狀態(tài)。

4.2 穿透與擴展能力分析

選擇Michalewicz 函數(shù)來說明MDO 算法的穿透能力和擴展能力。Michalewicz函數(shù)的形式為：

該函數(shù)具有如表9 所示的特性。Michalewicz 函數(shù)有n!個局部極值點，當n=80 時，其局部極值點有7.156 9×10118個，因此該函數(shù)極難求得最優(yōu)解。

圖3展示了MDO算法在時間間隔[0,694]中求解Michalewicz函數(shù)時的穿透、膨脹和粘滯狀態(tài)的分布，計算參數(shù)為n=80,E0=0.005,K=3,T=4,τ=3,L=3,N=200,G=80 000。從圖3可以看到，當MDO算法求 解Michalewicz 函數(shù)時，CountPs(0,694)=103,CountEs(0,694)=297,CountSS(0,694)=40，需要23.41%的時間來進行穿透，67.50%進行膨脹，只有9.09%的概率落入粘滯狀態(tài)。

圖4 展示了當MDO 算法求解Michalewicz 函數(shù)時，在時間間隔[0,694]中出現(xiàn)的穿透、膨脹和粘滯狀態(tài)的數(shù)目。從圖4 可以看出，當MDO 算法求解Michalewicz函數(shù)時，穿透和膨脹狀態(tài)從0到250 s交替出現(xiàn)，但穿透狀態(tài)從0到250 s總是占優(yōu)；在250 s到350 s之間，只有膨脹狀態(tài)發(fā)生；在350 s 之后，膨脹和粘滯狀態(tài)交替出現(xiàn)，但出現(xiàn)粘滯狀態(tài)的概率增加，膨脹狀態(tài)總是在350 s之后占優(yōu)。這意味著穿透和膨脹從0到250 s依次進行；在250 s到350 s之間，搜索接近全局最優(yōu)，僅需要進行膨脹；當搜索進行了350 s 時，粘滯狀態(tài)偶爾開始出現(xiàn)。

Fig.1 Convergence curves of all algorithms when solving 6 benchmark functions圖1 所有算法求解6個基準函數(shù)時的收斂曲線

圖5 展示了當MDO 算法求解Michalewicz 函數(shù)時，在時間間隔[0,694]中的穿透和膨脹之間的協(xié)調(diào)性。從圖5可以看出，當MDO算法求解Michalewicz函數(shù)時，穿透狀態(tài)從0到250 s占優(yōu)；當t=250 s 時，穿透狀態(tài)占優(yōu)，達到最大值；250 s 后，穿透狀態(tài)占優(yōu)開始減少；相反，膨脹狀態(tài)出現(xiàn)的概率開始增加。該結(jié)論與從圖4所得到的結(jié)論一致。

圖6 展示了當MDO 算法求解Michalewicz 函數(shù)時，時間間隔[0,694]中的膨脹和粘滯之間的協(xié)調(diào)性。從圖6可以看出，當MDO算法求解Michalewicz函數(shù)時，當t=350 s 時，膨脹狀態(tài)占優(yōu)達到最大值，之后，占優(yōu)狀態(tài)出現(xiàn)次數(shù)開始減小。相反，粘滯狀態(tài)出現(xiàn)的次數(shù)開始增加，但粘滯狀態(tài)出現(xiàn)的次數(shù)總是不占優(yōu)勢。

Fig.2 Evaluation of penetration and expansion capacity of MDO algorithm圖2 MDO算法的穿透和擴展能力評價

Table 9 Characteristics of benchmark function F3 and Michalewicz function表9 基準函數(shù)F3 和Michalewicz函數(shù)的特征

Fig.3 Distribution of penetration,expansion and viscosity(ε=1.0E-11,λPs=1.0,λEs=1.0E-10)圖3 穿透、膨脹和粘滯狀態(tài)分布(ε=1.0E-11,λPs=1.0,λEs=1.0E-10)

Fig.4 Number of times of penetration,expansion and viscous state(λPs=1.0,λEs=10-10)圖4 穿透、膨脹和粘滯狀態(tài)出現(xiàn)次數(shù)(λPs=1.0,λEs=10-10)

Fig.5 Variation of coordination coefficient of penetration and expansion with time(CPs=1.0)圖5 穿透-膨脹協(xié)調(diào)系數(shù)隨時間的變化規(guī)律(CPs=1.0)

從以上分析可知，當MDO算法求解帶海量局部最優(yōu)解優(yōu)化問題和帶高條件數(shù)優(yōu)化問題時均表現(xiàn)出很好全局尋優(yōu)能力，且其局部尋優(yōu)能力和全局尋優(yōu)能力的平衡性把控得較好。

5 結(jié)束語

本文基于帶時滯影響的混雜食物鏈微生物培養(yǎng)動力學(xué)理論提出了一種具有全局收斂性的新型優(yōu)化算法，與其他典型群智能算法相比，MDO算法具有如下特點：

Fig.6 Variation of coordination coefficient of expansion and viscosity with time(CEs=1.0)圖6 膨脹-粘滯協(xié)調(diào)系數(shù)隨時間的變化規(guī)律(CEs=1.0)

（1）吸收算子、攫取算子、混雜算子和毒素算子可增加其搜索能力。

（2）采用隨機方法確定各算子中的相關(guān)參數(shù)，既大幅減少了參數(shù)輸入個數(shù)，又使模型更能表達實際情況。

（3）算法所涉及的微生物培養(yǎng)過程充分體現(xiàn)了多種微生物種群的濃度、培養(yǎng)液投放量及其營養(yǎng)物質(zhì)濃度、培養(yǎng)液定期投放量中含有的有毒物質(zhì)濃度、培養(yǎng)液投放周期、培養(yǎng)液向微生物轉(zhuǎn)化的滯后時間等參數(shù)的復(fù)雜變化情況。

（4）所有算子是通過利用具有時滯影響的混雜食物鏈微生物脈沖培養(yǎng)動力學(xué)理論以及微生物種群間的相互作用關(guān)系來進行構(gòu)造的，不需要與求解的優(yōu)化問題相關(guān)。因此，MDO算法具有很好的普適性。

（5）培養(yǎng)液及其所含有的有毒物質(zhì)的濃度的脈沖增加會導(dǎo)致試探解從一個狀態(tài)突然轉(zhuǎn)移到另外一個狀態(tài)，這種性質(zhì)有利于使搜索跳出局部最優(yōu)陷阱。

（6）算法考慮到了微生物種群培養(yǎng)過程中外界因素的不連續(xù)間斷介入的現(xiàn)象。

（7）算法所涉及的微生物種群之間的相互作用豐富多彩，體現(xiàn)了培養(yǎng)系統(tǒng)中微生物種群間的復(fù)雜捕食-被食關(guān)系和競爭關(guān)系。

（8）算法體現(xiàn)了微生物培養(yǎng)動力學(xué)模型中各參數(shù)的復(fù)雜變化情況。

（9）在進行迭代計算時，算法每次只處理每個種群特征數(shù)的1/1 000～1/100，從而使時間復(fù)雜度大幅降低。因此，MDO算法適于求解高維優(yōu)化問題。

MDO算法的下一步改進方向如下：

（1）利用微生物培養(yǎng)動力學(xué)模型優(yōu)化MDO算法的相關(guān)參數(shù)，使得這些參數(shù)設(shè)置更合理。

（2）深入研究吸收算子、攫取算子、混雜算子和毒素算子的動態(tài)特征。

（3）深入研究微生物種群的動態(tài)特征。