具脈沖出生和季節性捕殺的種群系統優化算法

黃光球，陸秋琴

西安建筑科技大學 管理學院，西安 710055

非線性規劃問題的全局最優解的求解是一個非常困難的問題，因為只有凸規劃且其數學表達式可導時，采用傳統優化理論才能獲得其全局最優解。若非線性規劃是非凸非凹或者是其數學表達式不可導，則傳統優化理論失效。為了解決該問題，人們發明了群智能優化算法來求取非線性優化問題的全局最優解[1-3]。此類算法對非線性優化問題的結構沒有特殊限制條件，因而具有廣泛的適應性。

傳統的群智能優化算法是依據一些特殊的自然現象構造而得，自然現象本身的特征對群智能優化算法的性能影響很大。因此，尋找具有優良特性的自然現象來構造群智能優化算法是提升群智能優化算法的有效手段。最近幾年，人們已提出了很多新型的群智能優化算法，如蝙蝠算法[4]、蜂群算法[5]、鯨魚算法[6]、蛙跳算法[7]、灰狼算法[8]、布谷鳥算法[9]等，這些算法的共同特征是算法所依據的自然現象非常簡單，無法用數學模型進行描述。由于存在此缺陷，關于算法的參數確定和性能分析變得異常困難。

為了解決上述問題，人們開始尋求能被數學模型很好描述的自然現象來開發新一代群智能優化算法，其中能夠被種群動力學數學模型[10-11]描述的自然現象就是其中的一種。文獻[12]采用兩個種群Lotka-Volterra 種群動力學模型構建出了第一個種群動力學算法，種群間的作用關系有競爭、互利、捕食-被食、融合、突變和選擇；以文獻[12]為基礎，文獻[13]構建出了基于3 個種群的Lotka-Volterra 種群動力學優化算法。

在自然界中，很多種群的增長是通過脈沖出生新生個體來完成的，一段時間后，從幼年種群成長而來的成年種群的捕獲也是季節性的。該自然現象非常普遍，且能被具脈沖出生和季節性捕殺的種群動力學數學模型[14-15]很好描述。本文正是依據該自然現象構造出了一種新的群智能優化算法——具脈沖出生和季節性捕殺的種群系統優化算法（population system optimization with impulsive birth and seasonal killing，PSO-IBSK）。與現有的種群動力學優化算法相比，本文算法具有如下特點：

（1）種群個體自然地劃分成幼年個體(簡稱幼體)和成年個體（簡稱成體）兩類，每類個體數是依據具脈沖出生和季節性捕殺的種群動力學模型自動進行動態計算而得，避免了人工確定個體數的困難。

（2）所有算子是通過幼體的脈沖出生、幼體演變為成體、成體的季節性捕殺、虛弱個體的死亡而構造出來的，符合種群動力學規律，且與所求解的實際優化問題無關，故具有很好的普適性。

（3）每個算子具有明確功能，其中出生算子和成長算子可分別實現成體向幼體瞬時和延遲傳遞信息，有助于搜索跳出局部陷阱；捕殺算子可周期性地將不良成體清除，死亡算子可將虛弱個體隨機清除，該兩個算子有利于提升算法的求精能力；強勢算子可實現強壯個體向虛弱個體擴散強壯信息，競爭算子可實現幼年和成體之間的有效信息交換，該兩個算子有利于提升算法的探索能力；進化算子可確保算法具有全局收斂性。

（4）采用具脈沖出生和季節性捕殺的種群動力學模型確定PSO-IBSK 算法中的相關參數，使算法的參數確定具有科學性。

（5）計算過程中，PSO-IBSK 算法每次只處理個體特征數的6‰～8%，從而使時間復雜度大幅降低。

1 PSO-IBSK 算法設計方法

假設要求解的優化問題如下：

式中，X=(x1,x2,…,xn)為變量，在求解過程中X的不同取值稱為試探解，Rn是n維歐氏空間；f(X)為目標函數；H為搜索空間。

1.1 優化問題試探解的生物學解釋

在一個生態系統G中生活著一個生物種群，該種群由具有兩種階段狀態的若干生物個體組成，即階段狀態s=1 的幼體和階段狀態s=2 的成體，其中幼體經過一段時間的生長后會長大成為成體，而幼體是由成體出生的，幼體的出生持續時間很短，故稱之為脈沖出生。為了節省資源，降低種群中的個體密度，需要周期性（季節性）地對該生態系統中一些生長狀況不良的成體進行捕殺，以便提升種群的整體質量。

在時期t，假設該種群的幼體數和成體數分別為N1(t)和N2(t)，個體總數為N(t)=N1(t)+N2(t)；每個幼體和每個成體都用唯一編號表示，于是所有幼體和所有成體的編號集合分別為種群中的每個幼體和每個成體都具有n個特征，于是對于階段狀態為s的個體i來說，用其特征表示就是，其中就是階段狀態為s的個體i的第j個特征，i∈Cs(t)，j=1～n，s=1～2。

搜索空間H中優化問題式（1）的試探解與生態系統G中種群的幼體和成體的對應關系如下所述。

時期t，在式（1）的搜索空間H中隨機生成N(t)個試探解，即∈Cs(t),s=1～2}，其中。將搜索空間H視為生態系統G，則時期t在生態系統G中階段狀態為s的個體i就與搜索空間H中的試探解一一對應，也就是階段狀態為s的個體i的特征與試探解的分量相對應。在生態系統G中，生物個體的動力學演化規律總結如下：

（1）幼體是由成體以脈沖方法突然產生的；

（2）幼體經過時間段T后成長為成體；

（3）在某個季節內，某些生長狀況不良的成體會被捕殺掉；

（4）某些虛弱的幼體和成體會死亡；

（5）幼體和成體的生長狀況越好，繼續生存下去的概率會越高，即個體的演化規律符合達爾文進化論規律。

個體的生長狀況由優化問題式（1）的目標函數值描述，個體的生長狀況越好，其對應的目標函數值就越小。時期t階段狀態為s的個體i的生長狀況用IGI（individual growth index）指數來表示，其計算方法為：

1.2 種群系統動力學模型

在生態系統G中，假設幼體的出生具有規則脈沖性，則具有階段結構和脈沖出生的單種群離散模型為[16]：

其中，t∈Z+，Z+為非負整數集；令從t到t+1 的幼體的出生數=0，b>0，當t=kω（k∈Z+，ω為正數）時，幼體數N1(t)增加了N2(t)；δ為幼年種群的成長率，0<δ<1；幼年和成年的死亡率分別為d1和d2，0

若考慮對成體的捕殺，E為對成體進行捕殺的捕殺率，0

考慮脈沖出生，并考慮對式（4）的成體進行季節性捕殺，不失一般性，設對成體的捕殺發生在(T1/ω,T2/ω]（0 ≤T1

迭代N1(t)和N2(t)，可得式（5）在脈沖區間內的解析解，即：

式（6）在脈沖區間成立且α≠μ，α≠β。為簡單起見，下面總是設α≠μ，α≠β。當t=(m+1)ω時，由式（6）得：

定義內稟再生數R0為R0=bp/(1-r)(1-q)，令b0=(1-r)(1-q)/p。當R0<1 時，如果在平均數量上，個體在死亡前沒有得到替補和補充，則種群會走向滅絕。式（8）滿足：

當R0>1 時，存在一個正平衡點E*(u*，v*)，其中：

定理1[17]令γ=若b>b0，則式（8）存在一個正不動點E*(u*,v*) ；若b0bc，則E*(u*,v*) 是不穩定的。若b作為一個分支參數，則b=bc是一個Flip 分支。

1.3 算子設計

PSO-IBSK 算法的算子是通過各生物個體及其相互之間的作用關系構造而成。令：

（1）出生算子。該算子描述的是由成體產生幼體。首先生成階段狀態為s的強壯個體集合SIs(t)，該集合中的個體的IGI 指數高于同一階段狀態的IGI 指數的平均值，SIs(t)集合中的個體數為L，L稱為特征個體數；然后，從SI2(t)中隨機選出兩個成體i1和i2，由其產生的幼體數為B(t+1)個，B(t+1)=，即：

式中，i=1～B(t+1) ；MF(t) 為從{1,2,…,n} 中以概率Z0隨機選擇所形成的特征編號集合；Z0稱為個體特征受影響的最大概率；λ=Rand(-1,1)，Rand()為均勻分布隨機取值函數；ρ=Rand(0,1)。

由式（11）產生的B(t+1) 個幼體的編號集合為BI(t+1)={i0,i0+1,…,i0+B(t+1)-1}，i0=|C2(t)|+1。

一方面，因幼體是通過兩個強壯個體雜交產生的，故幼體的質量普通較高；另一方面，因幼體是瞬時產生的，故一定數量的優質個體的瞬時投放有助于搜索跳出局部最優解陷阱。

（2）成長算子。該算子描述的是幼體經過一段時間后長大為成體。假設當前幼體i以概率δ要轉變為一個成體，為了使該幼體具有成體的一些特征，首先，在階段狀態為s的個體中隨機選擇L個個體，由其編號形成的集合為GIs(t)，然后，將GI2(t)中的成體的部分特征傳給該幼體，使其成為成體，即：

然后，將幼體i從幼體集合中刪除，即C1(t+1)=C1(t)-{i}，將該個體的階段狀態改為成年狀態，其編號為|C2(t)|+1，即C2(t+1)=C2(t)+{|C2(t)|+1}。

成長算子描述了質量較高的幼體演變為成體的過程。成體是由高質量的幼體演變而來，若成體的質量不斷得到提升，則意味著搜索向全局最優解不斷靠近；若成體質量沒得到提升，則會被死亡算子和捕殺算子清除掉。因此，在這兩個算子輔佐下，成長算子有利于提升本算法的全局尋優能力。

（3）捕殺算子。該算子描述的是將生長狀態不良的成體進行人為刪除，刪除數目為E個。從成體集合選擇IGI 指數最低的E個個體，其編號形成的集合為WI2(t)，將WI2(t)中的所有個體進行人為清除：

（4）死亡算子。該算子描述的是虛弱個體的自然死亡。對當前階段狀態為s的個體i，若該個體的IGI 指數低于處于該階段狀態的個體平均IGI 指數，則以概率ds令其死亡，即：

（5）強勢算子。該算子描述的是強壯個體的特征向虛弱個體擴散的現象，即強壯的個體對虛弱的個體產生影響。對階段狀態為s，個體編號為i的當前個體來說，有：

（6）競爭算子。該算子描述的是幼年和成體之間的相互競爭現象。對階段狀態為s，個體編號為i的當前個體來說，有：

（7）進化算子。該算子描述的是個體的進化需滿足達爾文進化論規律。對于階段狀態為s的當前個體i，其進化算子的定義如下：

式中，函數IGI()按式（2）計算。

1.4 PSO-IBSK 算法構造

算法1PSO-IBSK

步驟1初始化：

（1）令演化時期數G=105～108，誤差要求ε=10-5～10-10，L，Z0；N1(0)=N2(0)=200，m=0。

（2）依據2.1 節確定參數b、δ、d1、d2、ω、T1、T2。

（3）在搜索空間H中隨機生成試探解集：

（4）確定個體編號集合C1(0)={1,2,…,N1(0)}，C2(0)={1,2,…,N2(0)}。

（5）以X(0)為基礎，找出初始全局最優解Y*0。

步驟2執行下列操作：

步驟3結束。

1.5 算法特點分析

1.5.1 時間復雜度

PSO-IBSK 算法的時間復雜度計算如表1 所示，其中

1.5.2 PSO-IBSK 算法的全局收斂性分析

Table 1 Time complexity表1 時間復雜度

不失一般性，令f1即為所求的全局最優解。由式（18）的下標形成的集合為U={1,2,…,Z(t)}。

?X∈H，有f1≤f(X)≤fZ(t)，將H劃分為如下非空子集={X|X∈H且f(X)=fi}，i=1～Z(t)，顯然有：

令Xi,j(i=1～Z(t),j=1～表示中的第j個狀態；個體從狀態(i,j) 轉移到狀態(k,l) 表示為Xi,j→Xk,l；設pij,kl、pij,k為從Xi,j分別到Xk,l、中任一狀態的轉移概率，pi,k為從中任一狀態到中任一狀態的轉移概率，則：

引理1在PSO-IBSK 算法中，,i=1～Z(t),j=，滿足：

（1）引理式（20）的證明：設狀態i為時期t個體i的狀態，其空間位置為，由式（17）知，該個體具有適應度遞增特性，故在時期t+1 有：

（2）引理式（21）的證明：設狀態i為時期t個體i的狀態，在時期t+1，個體i隨機選各算子進行演化以便轉移到更好的狀態k上。此時，存在有如下兩種情況：

①若狀態i=1，即全局最優狀態，因下一步不會轉移到較差的狀態上去，故必以概率p1,1=1 留在原狀態i上。因p1,1=1>0。命題得證。

②若狀態i≠1，則在狀態1 和當前狀態i之間必至少存在一個中間狀態k，使得f1≤fk0。命題得證。

綜上所述，可得?k0。證畢。

定理2[18]設P′是一n階可歸約隨機矩陣，即通過相同的行和列變換后可得到，其中C是m階本原隨機矩陣，且T≠0，R≠0，則有：

上述矩陣是一個穩定隨機矩陣，P′∞=1′P′∞，P′∞=P′0P′∞唯一確定且與初始分布無關，P′∞滿足條件：

定理3PSO-IBSK 算法具有全局收斂性。

證明從各算子的定義式（11）～式（16）可知，與滿足關系，表明PSO-IBSK算法的演變過程具有Markov特性。每個,i=1～Z(t)是有限Markov 鏈上的一個狀態，根據式（20）可得狀態轉移矩陣為：

且P′的每行均滿足式（19）。根據式（21）可得：

于是，P′是一個Z(t)階可歸約隨機矩陣，滿足定理2 的條件，故有：

因C∞=C=(1)，T∞=0，故必有R∞=(1,1,…,1)T，于是：

上式表明，當k→∞時，pi,1=1，i=1～Z(t)，于是：

因此，PSO-IBSK 算法具有全局收斂性，證畢。

目前，進行群智能算法收斂性分析的方法有圖搜索法[19]、代數方法、解析分析法、狀態空間模型法和馬爾科夫分析法[20]等，這些方法的證明過程復雜，僅適合于特定算法的收斂性分析，沒有通用性。本文提出的群智能算法收斂性分析法是基于有限馬爾科夫鏈的可歸約隨機矩陣穩定性定理，證明過程非常簡單，可適用于所有群智能算法的收斂性分析，具有通用性。同時，本證明過程表明，對一個群智能算法來說，只要其種群個體演化符合達爾文進化論規律，該群智能算法就是全局收斂的。

2 算法參數確定

2.1 種群系統模型的參數設置方法

種群系統模型的參數是算法的內置參數，由1.2節介紹的理論進行設置，用戶無需修改。下面討論這些參數的設置方法。

首先討論捕殺時間T1和T2選取對式（9）描述的幼年和成體數量的影響，以及種群最大可持續捕殺量E。由定理1 可知，式（9）的正平衡點依賴于收獲的時間，盡管捕殺量E相同，太遲的捕殺將導致種群的滅絕。在每個繁殖季節后，收獲成體越早，系統可承受的捕殺量越大，對生存下來的個體來說，其存活的概率也越大。從生物學角度來看，捕殺成體可以減少種內競爭，使得其他個體得到充足的食物和空間，進而增加存活概率。因此，在繁殖季節結束后，收獲成體對具有脈沖出生的種群是比較有利的。同時也可得到成年種群的平衡態v*是關于T1的減函數。事實上，如果τ為一常數，則：

由于01 和β<μ，則則v*是關于T1的減函數。

由定理1 知，若b0

因此，只需考慮一個周期內的可持續捕殺量。不失一般性，取m=0，則可持續捕殺量為：

Fig.1 Changes of the number of juveniles and adults with time圖1 幼年和成體數隨時間的演變規律

2.2 人工控制參數設置方法

人工控制參數包括Z0和L，需要用戶根據所求優化問題的實際特征進行人工設置。下面通過Awad等人最新發布的智能優化算法測試包中所介紹的基準函數F6[21]，該函數由Scaffer 函數經旋轉和平移擴展而得，其表達式如下：式中，M為n×n維旋轉矩陣；O為n維向量，其中每一維的值在[-80,80]中隨機選取。

令n=50，Z0=0.01，G=108，PSO-IBSK 算法運行100 次，表2 描述了L與最優目標函數值的平均值(OFV)和計算時間的平均值(CT)之間的關系。表2 表明，當L=3～7 時，OFV的精度達到最佳，而CT遞增不大，建議L=3～7。

Table 2 Relationship of L with OFV and CT表2 L 與OFV 和CT 之間的關系

令n=50，L=3，G=108，PSO-IBSK 算法運行100 次，表3 描述了Z0、OFV和CT之間的關系。結果表明，當Z0=0.006～0.080 時，OFV精度較高，但CT增加不大，建議Z0=0.006～0.080。

Table 3 Relationship of Z0 with OFV and CT表3 Z0 與OFV 和CT 之間的關系

3 PSO-IBSK 算法與其他算法的比較

本文選用CEC2013[22]智能優化算法測試包中12個難度很大優化問題來對PSO-IBSK 算法與其他算法進行比較，如表4 所示。

Table 4 12 benchmark functions in CEC2013表4 12 個CEC2013 基準函數優化問題

求解這些優化問題時，PSO-IBSK 算法的參數設置是n=50，G=108，ε=10-8，Z0=0.01，M=3。與PSOIBSK 算法進行比較的7 種智能優化算法為BRKGA（biased random-key genetic algorithm）[23]、ACO-IM（ant colony optimization-influence maximization）[24]、SRPSO（self-regulating particle swarm optimization）[25]、SLIWBBO（self learned invasive weed-mixed biogeography based optimization）[26]、DE-DMSC（differential evolution with dual mutation strategies collaboration）[27]，DARSRWO（double adaptive random spare reinforced whale optimization）[28]、GABC（genetic artificial bee colony）[29]，這些算法各參數設置可參見其對應文獻。

求解各個優化問題時，每個算法均獨立求解51次。表5 給出了各算法的求解結果，表中計算結果是各算法求解各優化問題時的計算值與理論值之間的偏差。

Table 5 Results obtained by each algorithm表5 各算法的求解結果

Fig.2 Sample convergence curves圖2 樣本收斂曲線

從表5 可以看出，這8 個算法按最終排名1 和最終排名2 排序所得的結果均如下：

PSO-IBSK>DARSR-WO>DE-DMSC>BRKGA>ACO-IM=GABC>SRPSO>SLIW-BBO

圖2（a）～（f）給出了各算法求解優化問題F3、F7、F14、F20、F22、F23 時的樣本收斂曲線。

4 結束語

本文基于具脈沖出生和季節性捕殺的種群動力學模型提出的PSO-IBSK 算法能夠動態自動計算幼年和成體數，從而使得種群個體數的確定具有科學性；PSO-IBSK 算法中的算子是依據幼體的脈沖出生、成體的長成和季節性捕殺、個體間的競爭關系和虛弱個體的死亡而構造出來的，符合種群動力學規律。各個算子分工明確，個體間信息交換充分，出生算子和成長算子的綜合作用有助于搜索跳出局部最優解陷阱，強勢算子和競爭算子的綜合作用有利于提升算法的求精能力和探索能力，而進化算子能確保算法具有全局收斂性。PSO-IBSK 算法中的大量參數采用具脈沖出生和季節性捕殺的種群動力學模型確定，提升了算法參數確定的合理性和科學性。PSO-IBSK 算法每次處理的個體特征數很少，適應于求解維數較高的優化問題。

PSO-IBSK 算法的下一步改進方向如下：

（1）深入研究各算子的動態特征；

（2）深入研究個體的動態特征。