種群具有離散Leslie年齡結構的動力學優化算法

黃光球，陸秋琴

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

群智能優化算法[1]是啟發式算法的一種，其基本求解策略不依賴優化問題的特征（如梯度、多次可微性等），而是依賴于試探解之間的信息交換能力。群智能優化算法之所以采取此策略，是因為優化問題在很多情況下無法求得其梯度或者不具備任何可微性。為了使得試探解之間的信息交換具有合理性，群智能優化算法經常依據一些特殊的自然現象構造而得。因此，一種自然現象的質量優劣對群智能優化算法的性能影響重大。

目前，具有代表性的新一代群智能算法有烏鴉算法[2]、海洋掠食者算法[3]、完全遷徙型生物地理學算法[4]、人工藻類算法[5]、磷蝦群算法[6]、鯨魚算法[7]、雜交蜂群算法[8]、反向學習粒子群算法[9]、花授粉算法[10]等。這些算法均來自一些特殊的自然現象，但這些自然現象有一個共同的缺陷就是無法用數學模型進行描述。因此，基于這些自然現象所構造出來的算法的性能設計無法奠定在堅實的數學分析基礎之上，該缺陷也是傳統群智能算法的缺陷。此外，因這些算法沒有數學模型基礎，故相關參數選擇也相當困難。于是，若一個自然現象能夠被數學模型恰當描述，基于該自然現象的群智能算法就可以克服上述缺陷。

一個自然現象能否成為構造一個群智能優化算法的生物活動場景，必須包含有兩個條件：一是生物個體不斷進化，符合“適者生存，優勝劣汰”的自然規律；二是生物個體之間存在相互作用關系。第一個條件用來構造算法的邏輯結構并確保該算法具有全局收斂特征，第二個條件用來構造算子，個體類型越多，個體之間相互作用關系越豐富，所能構造出來的算子也就越多，算子的性能也越強。此外，若該自然現象還能被數學模型恰當描述，則據此所構造出來的群智能優化算法的性能分析和參數選擇將具有很好的科學依據[11]。

由具有不同性別、不同年齡的生物個體所組成的生物種群在一個生態環境中生存進化，是一種很常見的自然現象[3]，該現象就滿足構造一個群智能優化算法的兩個條件。不僅如此，該自然現象既能被數學模型很好描述，又能提供高度多樣化的個體間相互作用關系。本文正是依據該自然現象提出了一種新的群智能算法，即種群具有離散Leslie 年齡結構的動力學優化算法，簡稱PDO-DLAS（population dynamic optimization algorithm with discrete Leslie age structure）算法，該種群動力學模型簡稱為Leslie模型[12]。Leslie 模型具有較好的應用價值，已得到較好的研究，目前研究進展如下：

（1）種群特征分析[13-18]。在該方面，通過構建Leslie模型，研究種群的特征，如擴散特征、穩定性特征、持久和湮滅特征等。

（2）種群控制研究[19-22]。在該方面，通過構建Leslie控制模型，研究種群規模的最優控制或收獲問題。

PDO-DLAS 算法依賴于Leslie 模型，因此上述模型的研究可為PDO-DLAS 算法的構造奠定良好的理論基礎。與文獻[2-10]描述的算法相比，該算法具有如下特點：

（1）Leslie 年齡結構使得個體自動依據年齡分成多個類別，可大幅提升個體的多樣性，從而降低個體陷入局部陷阱的概率。

（2）算法參數的選擇依據Leslie 數學模型確定，因而具有堅實的理論依據。

（3）依據Leslie 模型描述的生態場景設計出來的算子具有充分的理論依據。

（4）個體進化滿足達爾文進化論規律，既可確保算法具有全局收斂性，又可使得個體進化每次只涉及極少特征的改變，可大幅提升計算速度。

由于PDO-DLAS 算法具有獨特的優勢，對其進行研究很有意義。

1 PDO-DLAS 算法構建原理

考慮一類優化問題：

式中，Rn是n維歐氏空間，n為優化問題的維數；X=(x1,x2,…,xn)是一個n維解向量；f(X)為目標函數；H為搜索空間。

1.1 決策變量與生物個體的對應關系

假設某種群在一個生態系統E中生長，該種群由具有不同年齡的生物個體組成，為簡單起見，下面將生物個體簡稱為個體。個體分為雌性（用性別變量u=1 表示）和雄性（用性別變量u=2 表示）兩種類型；種群自然地分成A個年齡組，A為該種群的最高年齡，0

在時期t，該種群中的性別為u的個體總數為Pu(t),u∈{1,2}，即：

每個個體都用唯一編號表示，即性別為u、年齡組為a的個體編號集合為{1,2},a=1～A。種群中的每個個體都具有n個特征，且特征數n不隨時間變化而改變。時期t，對于性別為u、年齡組為a的個體i來說，用其特征表示就是其中就是性別為u、年齡組為a的個體i的第j個特征，

優化問題式（1）中的決策變量與個體的對應關系如下所述。

時期t，在式（1）的搜索空間H中隨機生成的Pu(t)個試探解是令搜索空間H與生態系統E相對應，則時期t該生態系統中Pu(t)個個體就與搜索空間H中式（1）的Pu(t)個試探解一一對應，即時期t性別為u、年齡組為a的個體i與試探解一一對應；更進一步，即時期t性別為u、年齡組為a的個體i的特征與試探解的分量相對應，u∈{1,2}。

在該生態系統中，個體的演化規律如下：

（1）個體的新生是由雌性個體和雄性個體交配產生的。

（2）個體的死亡是因為個體太虛弱所導致的。

（3）在年齡相近的個體中，普通個體會不分性別地向比它優秀的個體學習。

（4）年齡較大的個體能夠不分性別地對年齡較小的施加影響。

（5）個體之間的相互作用表現在對其特征的影響上，且此相互作用及其所導致的影響是隨時間變化的。

對任何一個個體來說，按優化問題式（1）計算，其所對應的目標函數值越小，則其適應度越高，反之亦然。個體的適應度越高，其生存的概率也越高；而適應度越低的個體，其被淘汰的概率會越高。也就是說，生態系統中的個體的演化規律符合達爾文進化論規律。時期t性別為u、年齡組為a的個體i的適應度用FEI（fitness evaluation index）指數來表示，其計算方法為：

1.2 Leslie種群離散模型

假設表明第a年齡組每一個個體在一個時段內出生的后代（或新生的雌性或雄性）的平均數為出生率函數，特別地設它是一個常數，即：

由所有年齡組的后代相加而得到初始年齡組為：

假設描述由第a年齡組存活到第a+1 年齡組的比率為存活率函數它只是第a年齡組個體數的線性函數。

假設式（4）和式（6）表明不考慮參數依賴于環境條件的變化，并且忽略種群大小對出生率和死亡率的影響。

令Yu(t)=，則由式（5）、式（7）可得到方程：

這里，A×A方陣Lu形如：

Lu稱為Leslie 矩陣，其元素稱為活率或種群統計參數。稱式（8）、式（9）為Leslie模型。

向量Yu(t)以絕對值而非相對值顯示了種群大小在各年齡組中的一個分布，相對年齡分布稱為年齡普。矩陣Lu確定了A維歐氏空間的一個線性算子，把這種算子稱為Leslie算子。因為變量代表第a年齡組的個體數，它們是非負的。矩陣Lu的特征方程為：

如果矩陣Lu有模等于r的h個特征值，即λ1=r,λ2,…,λh，則所有這些特征值彼此不同，但都是λhrh=0 的根，h稱為矩陣Lu的非本原性指數，它等于那些出生率不為0 的年齡組的編號的最大公約數。若h=1，則矩陣Lu稱為素矩陣。

矩陣Lu可表示為：

式中，m為矩陣Lu的不同特征值的個數；Jj為kj×kj（kj為特征值λj的重數，k1+k2+…+km=A）Jordan塊，其形狀為：

Cu(λ)和Ru(λ)的內積滿足：

(Cu(λa),Ru(λa))=1,(Cu(λa),Ru(λb))=0,a≠b,a,b=1～A

引入極限向量函數：

式（10）、式（11）可以對Leslie 模型的軌線的漸近性質建立一些推斷。

若Lu是素矩陣，則極限向量函數與t無關，即：

定理1（關于循環的平均）[12]設h為矩陣Lu的非本原性指數，則：

定理2（關于周期性）[12]設h為矩陣L的非本原性指數，則極限向量函數L(Yu,t)關于t是周期的，周期為：

式中，lcm 表示最小公倍數；l取遍值1,2,…,h;ul/vl為不可約分數（當l=1 時，設ul=0,vl=1）。

這樣，T是非本原性指數h的因子，并依賴于向量Yu的初始分布。任何情況下，極限分布的總體大小W(Yu,t)的真正周期必須是T的因子，其中W(Yu,t)=真正周期等于T的條件，由下面的定理給出。

定理3（關于種群總大小的周期）[12]W(Yu,t)的周期一致的充要條件是：

當R的值不作為穩定性的某種度量時，是所有生物種群再生速率的廣義參數；而作為一種穩定性度量時，如果R=1，那么r=1，種群大小不會作指數性的減少和增長，而趨于極限分布狀態。R<1 和R>1 分別表明r<1 和r>1，根據式（10）和式（12），種群或滅亡或者無限制地增長。

由式（14）知，成長系數R由存活率、出生率函數和年齡組數A決定。當確定后，可確定能夠確保R=1 的A。

定理4（循環周期）[12]設Lu是一具有單位主特征值和非本原性指數h的Leslie 矩陣，則在模型式（8）的軌線之間存在循環，如果循環有周期k，則k是h的因子。

1.3 演化算子設計

時期t，對當前性別為u、年齡組為a的個體i來說，其特征個體集合的生成策略如下：

（1）優勢個體集合SIu,a(t)。從個個體中隨機選出M個個體，這些個體的FEI 指數要比當前個體要高，形成優勢個體集合SIu,a(t)，M稱為特征個體數。

（2）普通個體集合GIu,a(t)。從個個體中隨機選出M個個體，形成普通個體集合GIu,a(t)。

演化算子設計方法如下：

（1）新生算子。該算子描述種群內個體的新生。從雌性和雄性優勢個體集合SIu,a(t)中各隨機選擇一個個體進行交配，產生新個體i，其性別隨機確定，所處年齡組為第1 組。

式中，w在{1,2}隨機選擇；α=Rnd(0,1),Rnd(a,b)為在區間[a,b]內產生滿足均勻分布的隨機數；a=A0～A,A0為雌性個體生育的起始年齡；iu在中隨機選擇，u∈{1,2};J0為在1～n中隨機選擇M個特征所形成的集合。由式（15）形成的新個體為

（2）死亡算子。該算子描述種群內個體的死亡。時期t，從性別為u、年齡組為a的個體集合中，選出FEI指數最低的個體，將其刪除。

（3）學習算子。該算子描述的是處于相近年齡的普通個體會向其優勢個體學習，從而使得優勢個體的一些特征傳遞給普通個體。與年齡組a相近的年齡組取為U={a-M,a-M+1,…,a,a+1,…,a+M}，對于每個年齡組b∈U，任意選擇被學習的個體性別w，令Q(b)=w,w∈{1,2}；進一步令則對于性別為u、年齡組為a的當前個體i來說，有：

（4）影響算子。該算子描述的是年齡較大的個體的一些行為會對年齡小的個體造成影響。比年齡組a大的年齡組取為V={a+1,a+2,…,a+M}，對于每個年齡組b∈V，任意選擇施加影響的個體性別w，令Z(b)=w,w∈{1,2}；進一步令則對于性別為u、年齡組為a的當前個體i來說，有：

式中，YP=GP1?GP2。

（5）進化算子。個體的進化須滿足達爾文進化論規律。對于性別為u、年齡組為a的當前個體i，其進化算子的定義如下：

1.4 個體初始化

假設優化問題搜索空間的維數為n，每個變量的搜索區間為[li,ui],i=1～n，利用正交拉丁方生成方法產生W個初始解的正交表LM(Wn)的構造算法如算法INIT 所述。

算法1INIT(W)//產生W個初始解的正交表LM(Wn)的構造算法

步驟1計算每個變量的離散點yij：

yij=li+(j-1)(ui-li)/(W-1),i=1～n;j=1～W

步驟2根據正交拉丁方的生成方法計算初始解xij：

xij=yjk,i=1～W,j=1～n

式中，k=(i+j-1)modW；若k=0，則k=W。

算法INIT(W)所確定的W個初始解Xi=(xi1,xi2,…,xin),i=1～W，具有很好的均衡分散性和整齊可比性。

1.5 PDO-DLAS 算法構造

算法2PDO-DLAS

步驟1初始化

步驟2執行下列操作

1.6 算法特點分析

1.6.1 時間復雜度

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

1.6.2 PDO-DLAS 算法的全局收斂性證明

定理5PDO-DLAS 算法具有全局收斂性。

證明（1）從新生算子、學習算子、影響算子、生長算子的定義式（15）～式（18）可知，滿足關系表明PDO-DLAS 算法的演變過程具有Markov 特性；（2）從生長算子的定義式（18）知，任一個體的FEI 指數總是保持遞增態勢。依據文獻[23]可知，滿足上述兩個特性的PDO-DLAS算法具有全局收斂性，其相關證明可參見文獻[23]，本文不再贅述。

Table 1 Time complexity表1 時間復雜度

2 參數確定

PDO-DLAS 算法的參數由Leslie 模型參數和運行控制參數兩部分組成，前者是算法的內置參數，一旦設定，無需修改；后者需要根據所求優化問題的實際特征進行人工設置。

2.1 Leslie模型參數設置方法

依據1.2 節的結論，Leslie 模型參數設置應滿足如下要求：

（2）從定理2 和定理3 知，軌線的性質依賴于時間尺度的選擇，故應避免模型軌線的周期性。例如，可合并時間間隔或延長時段，使得有非零出生率的年齡組成為相鄰的。

（3）從定理4 知，應盡量使得Lu的主特征值不等于1，從而確保模型式（8）的軌線之間存在循環。

依據上述要求，Leslie 模型式（4）、式（6）和式（9）中的參數可按下列方式設置：A=20,A0=A/2,h=A0+1=11,n0=30,n1=50;b0=3,b1=5;s0=0.65,s1=0.846。此時，個體總數隨演化周期的變化如圖1所示。由式（14）計算得R=1.000 378 ≈1;Lu的主特征值λ=1.000 3 ≈1，式（8）的軌線之間的循環周期為11 個演化周期。從圖1 可知，各年齡組個體數趨向極限分布，且具有很好的變化特征。

Fig.1 Changing rule of total number of individuals with time圖1 個體總數隨時間的變化規律

2.2 運行控制參數設置方法

含有n個變量的Michalewicz 函數優化問題具有n!個局部極小點。例如，對于50 維優化問題，Michalewicz 函數有高達50!=3.041 41×1064個局部極小點。因此，搜索過程極易陷入局部最優解陷阱。目前，該函數的理論全局最優解還未可知，尋找該函數的全局最優解具有很大的挑戰性。Michalewicz 函數優化問題如下：

本文以求解Michalewicz 函數優化問題為例來確定PDO-DLAS 算法中的運行控制參數的適當設置，這些參數包括如下幾類：

（1）G和ε。這兩個參數是互補參數，只要滿足其中一個即可。通常取G=106～1010,ε=10-10～10-5。

（2）E0和M。這兩個參數對PDO-DLAS 算法的性能影響較大，下面重點討論這兩個關鍵參數的取值規律。

令n=50,E0=0.01,G=108,PDO-DLAS 算法運行100 次，表2 描述了M與最優目標函數值的平均值（AvgOFV）和計算時間的平均值（AvgCT）之間的關系。表2 表明，當M=3～6 時，AvgOFV的精度達到最佳，而AvgCT遞增不大。因此，建議M=3～6。

Table 2 Relationship of M with AvgOFV and AvgCT表2 M 與AvgOFV 和AvgCT 之間的關系

令n=50,M=3,G=108,PDO-DLAS 算法運行100 次，表3 描述了E0、AvgOFV和AvgCT之間的關系。結果表明，當E0=0.004～0.100 時，AvgOFV精度較高，但AvgCT增加不大；當E0>0.100 時，AvgCT增加很大，但AvgOFV精度也大大降低；尤其當E0=1.000時，無法得到最優解。因此，當E0=0.004～0.100 時，PDO-DLAS 算法性能最好。

Table 3 Relationship of E0 with AvgOFV and AvgCT表3 E0 與AvgOFV 和AvgCT 之間的關系

3 PDO-DLAS 算法與其他算法的比較

CEC2013[24]是一個國際上專門用來測試智能優化算法的測試包，該測試包中包含有28 個優化問題，每個優化問題由一些著名的傳統優化問題組合而成。這些優化問題不但克服了傳統優化問題的一些缺陷（如全局最優解關于原點對稱、條件數過低等），而且其求解難度均大幅提高。本文在該測試包中選擇6 個難度很大的優化問題來對PDO-DLAS 算法與其他算法進行比較，如表4 所示。

Table 4 6 optimization problems表4 6 個優化問題

表4 中，O是隨機產生的n維向量。求解這些優化問題時，PDO-DLAS 算法的參數設置是n=50,G=108，ε=10-8,E0=0.01,M=3。與PDO-DLAS 算法進行比較的7 種智能優化算法為RCGA（real-coded genetic algorithm）[25]、DASA（differential ant-stigmergy algorithm）[26]、NP-PSO（non-parametric particle swarm optimization）[27]、MpBBO（metropolis biogeographybased optimization）[28]、MDE-LiGO（modified differential evolution with locality induced genetic operators）[29]、SLADE（self-adaptive differential evolution）[30]、ABC（artificial bee colony）[31]，這些算法各參數設置可參見其對應文獻。

求解各個優化問題時，每個算法均獨立求解51次。表5 給出了各算法的求解結果。其中，表5 的列是所有參與比較的優化算法，表5 的行是求解每個優化問題時所獲得的關鍵對比參數，即各算法所獲得的最優解的均值、中值、標準誤、最小值、最大值、適應度評價次數。

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

PDO-DLAS>MDE-LiGO>SLADE>NP-PSO>RCGA>DASA=MpBBO>ABC

按最終排名2 進行排序所得的結果均如下：

PDO-DLAS>MDE-LiGO>SLADE>NP-PSO>RCGA=DASA>MpBBO>ABC

圖2（a）～（f）給出了各算法求解6 個優化問題時的樣本收斂曲線。各圖樣本曲線之間的差異可以用PDO-DLAS 算法的探索和求解過程來說明。

Table 5 Solution results of each algorithm表5 各算法的求解結果

表5（續）

（1）對于圖2（a），PDO-DLAS 算法在區間[0，102]、[102，103]、[103，107]、[107，108]搜索期間，搜索分別處于求精狀態、探索狀態、再求精狀態、快速探索狀態轉求精狀態，而進行“求精→探索→再求精→快速探索”的狀態轉換等價于不斷從局部最優解陷阱跳出的行為；其他算法在區間[0，102]、[102，106]、[106，108]處于初始化狀態、求精狀態、探索狀態轉再求精狀態。PDO-DLAS、MDE-LiGO、SLADE 具有最終的相同求解精度，其他算法最終的求解精度均很低。

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

（2）對于圖2（b），PDO-DLAS算法在區間[0，102]、[102，5×105]、[5×105，9×105]、[9×105，108]搜索期間，搜索分別處于求精狀態、緩慢探索狀態、快速探索狀態、再求精狀態，即進行“求精→緩慢探索→快速探索→再求精”的狀態轉換，實施相應的局部最優解陷阱逃逸操作；DASA、MpBBO、MDE-LiGO、SLADE 算法在區間[0，102]、[102，5×105]、[5×105，9×105]、[9×105，108]搜索期間，處于初始化狀態、緩慢探索狀態、快速探索狀態、求精狀態；RCGA、NP-PSO、ABC 算法在區間[0，102]、[102，9×107]搜索期間，搜索分別處于初始化狀態、緩慢求精或緩慢探索狀態，RCGA 在[9×107，108]期間進入快速探索狀態。PDO-DLAS、MDELiGO、SLADE 具有最終的相同求解精度，其他算法最終的求解精度均很低。

采用上述類似的分析方法，可以分析圖2（c）～圖2（f），本文不再重復。

4 結束語

本文基于種群具有離散Leslie 年齡結構的動力學模型提出了一種具有全局收斂性的新型優化算法，與其他典型群智能算法相比，PDO-DLAS 算法具有如下特點：

（1）個體依據其年齡組和性別被自動劃分成2A類，每類個體數依據Leslie 模型自動進行動態計算，這樣既大幅增加了個體的多樣性，又避免了人工確定個體數的困難。

（2）所有算子是通過利用Leslie 模型以及同齡組和不同齡組個體間的相互作用關系進行構造的，PDO-DLAS 算法與所求解的實際優化問題無關，故具有很好的普適性。

（3）PDO-DLAS 算法中的每個算子具有明確功能，其中學習算子可實現年齡組相近個體之間的信息交換；影響算子可實現不同齡組個體之間的信息交換；新生算子可增加強壯個體數；死亡算子可以減少虛弱個體數；進化算子可確保算法具有全局收斂性。

（4）采用正交拉丁方方法生成初始試探解可確保其具有均衡分散性和整齊可比性。

（5）采用Leslie 模型的機理確定PDO-DLAS 算法中的相關參數，大幅減少了PDO-DLAS 算法參數的人工確定個數。實際上，PDO-DLAS 算法需要人工確定的參數只有2 個。

（6）在進行迭代計算時，PDO-DLAS 算法每次只處理個體特征數的1/250～1/10，從而使時間復雜度大幅降低。因此，PDO-DLAS算法適于求解高維優化問題。

PDO-DLAS 算法的下一步改進方向如下：

（1）深入研究學習算子、影響算子、新生算子和死亡算子的動態特征。

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