微擾動和混合變異的差分生物地理學優化算法

葛 強，李玉晶，喬保軍，左憲禹，王更科

(1.河南大學 數據與知識工程研究所，河南 開封 475004；2.河南大學 河南省大數據分析與處理重點實驗室，河南 開封 475004；3.中國科學院 空天信息創新研究院，北京 100101；4.中國科學院 中國科學院大學，北京 100101)

0 引 言

生物地理學優化算法是由Simon教授等[1]基于生物地理學數學模型提出的一種群智能優化算法。BBO不同于其它群智能優化算法之處在于其主要通過遷移算子來實現各個棲息地之間的資源通信，并通過變異算子來提高種群的多樣性。BBO算法在提出之后便得到了國內外智能優化算法研究學者的密切關注，眾多學者將其應用到實際工程優化問題的求解中，如電力系統[2]、圖像處理[3]、模式識別[4]、排列流水間調度[5]等。BBO算法的優點是收斂速度快、結構簡單，但尋優后期易過早收斂，陷入局部最優。針對這一缺點，眾多學者從不同的角度改進BBO算法來提高優化性能。王雅萍等[6]發現基于雙曲正切函數遷移率模型的BBO尋優能力更強。鄭夏等[7]引入Fibonacci迭代思想來消除種群重復項，差分進化來改進變異策略。張文輝等[8]提出基于自適應遷移算子和差分變異策略的BBO。唐繼勇等[9]提出動態選擇遷入地與自適應遷入策略。這些對BBO算法的改進均在某種程度上提升了算法的尋優能力，但在收斂精度、收斂速度和穩定性上還有很大的提升空間。本文提出的MDEBBO算法采用差分變異策略和自適應的局部擾動因子作為新的遷移算子，并采用一種基于高斯變異和柯西變異的根據迭代次數自調整的混合變異算子。仿真測試實驗數據顯示，MDEBBO算法在收斂速度和尋優精度上的優勢明顯。通過將MDEBBO算法應用到對Richards模型的參數估計上，實驗結果表明其在參數估計方面也有很好的表現力。

1 生物地理學優化算法

生物地理學是根據生物有機體在空間和時間上的地理分布對其進行研究的學科，探索了不同地理位置的生態系統中不同物種的遷移和突變。在BBO算法中，優化問題的可能解被稱為棲息地(Habitat)，每個棲息地都根據棲息地適應性指數(habitat suitability index，HSI)進行評估。其中，影響棲息地適應性指數的因素稱為適應性指數變量(suitability index variables，SIVs)。BBO算法被用來處理實際工程問題時，首先隨機產生n個棲息地作為最初解，然后使用遷移操作來實現棲息地之間的信息共享，并通過變異操作來增加解的多樣性，最后經過多次迭代操作獲取最優解。

1.1 初始化

首先隨機產生NP個可以表示為N維向量Hk(k=1,2,…,NP) 的候選個體，Hk(SIVj)(j=1,2,…,N) 表示為棲息地Hk的第j個SIV，Hk(SIVj) 初始化如式(1)所示

Hk(SIVj)=H(SIVjmin)+rand(H(SIVjmax)-H(SIVjmin))

(1)

其中，H(SIVjmax) 和H(SIVjmin) 分別為第j維變量的最大、最小限度。

1.2 遷移操作

在BBO算法中，棲息地通過遷移操作與其它棲息地共享信息。一個HSI值高的棲息地比HSI值低的棲息地容納更多的物種。將棲息地Hk按HSI值以優至劣的順序排序，棲息地Hk的物種數量按式(2)映射獲取

Sk=Smax-k

(2)

其中，Smax為棲息地最多可生存的物種數量。在BBO中，每個棲息地都有各自的遷入率λ和遷出率μ，這兩個參數控制著物種從一個棲息地向另一個棲息地的遷移。如式(3)和式(4)[7]所示，遷入率和遷出率取決于棲息地的物種數量

(3)

(4)

其中，I、E分別表示最大遷入率、最大遷出率，Sk為棲息地Hk的物種數量。首先根據遷入率λk挑選進行遷移操作的棲息地，然后對棲息地Hk的每個SIV根據遷出率μk按輪盤賭選擇法獲取待遷出的棲息地。最后按式(5)[7]完成遷移操作

Hk(SIVj)←Hm(SIVj)

(5)

其中，Hk(SIVj) 和Hm(SIVj) 分別為遷移棲息地和待遷出棲息地的j維分量。在BBO算法中，遷移算子的每次調用都會導致單個SIV從一個棲息地遷移到另一個棲息地中[10]。

1.3 變異操作

除了遷移，BBO還具有大多數進化算法所具有的概率特性，在算法的迭代過程中可能會發生突變[11]。變異算子根據棲息地存在的先驗概率隨機修改棲息地的SIV，其中棲息地Hk的物種概率Pk可由式(6)計算得出

(6)

通過對生物地理學分析得知，棲息地的突變概率mk和物種概率Pk成反比，棲息地Hk的突變概率mk如式(7)[12]所示

(7)

其中，mmax是棲息地Hk的最大突變概率，Pmax是所有物種概率的最大值。對于每個棲息地Hk，若隨機數rand不大于該棲息地的突變率mk， 則隨機選取一個規定界限內的值來替換棲息地Hk的某維解變量Hk(SIVj)。 BBO算法流程如圖1所示。

圖1 BBO算法流程

2 基于微擾動和混合變異的差分生物地理學優化算法

群智能優化算法的過程一般由勘探和開發兩部分組成，勘探就是在全局找到最優解可能出現的大致位置的過程，而開發則是在局部最優解附近進行精確搜索的過程[13]。為了增強生物地理學優化算法的優化能力并克服該算法不能很好平衡開發能力與局部最優解之間的矛盾，本文從遷移率模型、遷移算子、變異算子3個方面改進，下面進行詳細介紹。

2.1 非線性遷移率模型

BBO通過遷移操作實現棲息地之間的資源共享，遷入率、遷出率是決定棲息地進行遷移操作的主要因素，故不同的遷移率模型會直接影響該算法的尋優能力。文獻[6]介紹了3種新的非線性遷移率模型，并將其應用到BBO算法中。實驗結果表明越傾向于自然法則的模型越能提高算法的性能且基于雙曲正切變體遷移率模型的BBO算法比余弦遷移模型尋優能力強，其解更接近函數的全局最小值，故本文選取如圖2所示的雙曲正切變體遷移率模型。

圖2 雙曲正切變型遷移率模型[6]

從圖2可以看出，在雙曲正切變型遷移率模型下，棲息地中物種數量相對較少或較多時，遷移率變化較余弦遷移率模型平緩些，當該棲息地中的物種數量增加到一定程度后，遷移率變化速率隨之變快。雙曲正切遷移率模型的遷入、遷出率表達式分別如式(8)、式(9)所示，其中a為遷移率模型的影響因子，n為棲息地Hk的物種數量，文獻[6]將其取值為1.1，I和E分別為分別表示最大遷入率、最大遷出率

(8)

(9)

2.2 局部擾動和差分遷移算子

低開發能力和陷入局部最小化可以被認為是標準BBO的弱點，遷移算子的改進是增強算法優化能力的最佳方法之一。差分進化算法(DE)是一個優秀而強大的進化算法，它的一些差異策略提供了出色的全局搜索能力，引起了眾多學者的廣泛關注。通過分析式(5)得知，BBO的遷移算子只是棲息地之間SIV的簡單替換，這種遷移方式單一，導致進化過程中影響算法的探索性能。在DE的啟發下，MDEBBO將DE的變異策略引入到遷移算子中，新的遷移算子如式(10)所示。其中，Hbest是當代種群中HSI值最優的棲息地，Hindex1、Hindex2、Hindex3、Hindex4是從種群中隨機挑選的4個棲息地，且滿足index1≠index2≠index3≠index4≠i，F為縮放因子，通常取值為0.5。Hk(SIVj) 的某些信息來自Hbest， 因此Hk的位置能夠快速移動到搜索空間中的最佳位置，而且可從4個不同的棲息地獲取多樣化信息，增加了種群的多樣性。在迭代后期，BBO算法的種群多樣性會隨著尋優過程中逐步收斂而減小，且陷入局部極值的可能性增大

Hk(SIVj)=Hbest(SIVj)+F*(Hindex1(SIVj)-Hindex2(SIVj))+F*(Hindex3(SIVj)-Hindex4(SIVj))

(10)

為了進一步提高全局收斂精度，避免算法陷入局部極值，本文引入了一個微擾動算子來優化算法的性能。微擾動算子如式(11)所示

(11)

其中，σ為微擾動算子，Hi為輪盤賭方式選擇的待遷出棲息地，Gindex為當前迭代次數，Gmax為算法的最大迭代次數，rand(0,1) 表示從0到1的隨機數。由σ可以得到一個介于-0.5到0.5之間的隨機值，將這個隨機值作為權重與Hi和Hk之間的差分值相乘，得到擾動值。該策略通過引入正弦函數使σ的取值更具靈活性，可以實現以Hi(SIVj) 為中心的空間的精確搜索，提高算法的查找精度。在滿足BBO遷移規則的基礎下，添加一個約束條件(rand<0.2)，如果達到了條件則采用全局搜索能力更強的差分遷移算子(式(10))進行尋優，若不滿足，則采用微擾動算子來進行遷移操作(式(11))。改進的遷移算子結合DE的全局搜索和微擾動算子的局部搜索，在整體上提高了算法的開發能力，又由于通過引入正弦函數得出的微擾動因子增加了種群多樣性從而避免算法陷入局部最優，故能夠很好平衡開發能力與局部最優解之間的矛盾。局部擾動和差分遷移算子如算法1所示，其中輸入參數Population為種群，NP為種群大小，N為種群維數，λ和μ分別為棲息地的遷入率和遷出率，Pmod為棲息地修改概率；輸出參數Population為經過遷移操作后的種群。

算法1：遷移算子算法

輸入： (Population,NP,N,λ,μ,Pmod)

輸出： (Population)

(1) fork=1 toNPdo

(2) ifrand>Pmod

(3) continue;

(4) end if

(5) forj=1 toNdo //根據λ判斷Hk是否發生遷移操作

(6) ifrand<λthen

(7) ifrand<0.2 then

(8) 通過式(10)更新Hk(SIVj)

(9) else

(10) //根據μ選擇待遷出的棲息地Hi

(11) ifrand<μthen

(12) 通過式(11)更新Hk(SIVj)

(13) end if

(14) end if

(15) end if

(16) end for

(17) end for

2.3 混合變異算子

BBO算法的突變算子是通過隨機選取一個規定界限內的值來代替原棲息地的某維解變量。在這種策略下，雖然該算子可以生成不同的解，但現有的良好解很有可能會被隨機產生的值破壞而影響算法的收斂速度。為了提高解決方案的質量，本文引用進化規劃算法的高斯、柯西變異，將兩種變異算子結合迭代次數設置一種新的自調整的混合變異算子。高斯分布、柯西分布的概率密度函數分別如式(12)、式(13)所示

(12)

(13)

其中，μ和σ2分別表示高斯分布的均值和方差，x(x∈R)和t(t>0)都屬于尺度參數。當σ=1，μ=0時N(0,1)表示標準高斯分布， C(0,1) 表示標準柯西分布。在BBO算法中，高斯、柯西變異就是將一個服從其分布的隨機擾動項加到棲息地的某一分量上。根據高斯分布和柯西分布的特點，本文引入自調整混合變異算子，以此來提升算法的尋優速度和精度。混合變異算子如式(14)所示，其中Gindex表示當前迭代次數，Gmax表示算法的最大迭代次數。該策略將迭代次數作為權重系數來控制高斯變異和柯西變異之間的比例。該變異算子可以隨著進化過程自啟發的動態更新，在迭代初期，Gindex較小，通過柯西變異會獲取大幅度的變異步長，增加了候選解的多樣性，引領候選解快速向全局最優解移動，避免算法陷入局部最優解。在迭代后期，Gindex較大，高斯變異杰出的局部探索能力使候選解在局部范圍進行精確搜索，提高算法了的尋優精度。本文采取的變異策略只對整個種群中HSI值較差的后一半棲息地改變其適應性指數變量，以免變異算子破壞迭代周期內的較優解

(14)

混合變異算子如算法2所示，其中輸入參數Population為種群，NP為種群大小，N為種群維數，mk為棲息地的突變概率；輸出參數Population為經過變異操作后的種群。

算法2： 混合變異算子

輸入： (Population,NP,N,mk)

輸出： (Population)

(1) //只對整個種群中HSI值較差的后一半棲息進行突變

(2) fork=NP/2 toNPdo

(3) forj=1 toNdo

(4) ifrand

(5) 通過式 (14) 更新Hk(SIVj)

(6) end if

(7) end for

(8) end for

2.4 MDEBBO總流程

綜上所述，MDEBBO算法的流程如下：

MDEBBO算法步驟:

步驟1 初始化BBO參數：最大物種數量Smax， 最大遷出率E，最大遷入率I，最大變異概率mmax， 棲息地修改概率Pmod， 最大迭代次數Gmax， 種群大小NP，候選個體向量維度N等。隨機生成NP個候選個體，并按式(1)計算其適應性指數變量SIVs。

步驟2 計算每個候選個體的適應性指數HSI值，并將種群按優至劣順序進行排序。按式(2)計算每個候選個體可容納的物種數量S，分別按式(8)、式(9)計算其遷入率λ和遷出率μ。保留該迭代周期內最優的兩個候選個體作為精英個體。

步驟3 按算法1的遷移策略進行遷移操作。

步驟4 按式(6)計算候選個體的物種概率Pk， 按式(7)計算候選個體的變異率mk。 按算法2中的變異策略進行變異操作。

步驟5 重新計算候選個體的HSI值，將其重新進行排序。將保存的精英個體替換本迭代周期內最差的兩個候選個體。

步驟6 判斷是否達到終止條件，若是，結束并輸出最優候選個體；若否，則跳到步驟2。

3 仿真實驗與結果分析

為了驗證MDEBBO算法的尋優能力，通過對12個不同復雜度的基準函數進行大量的實驗。實驗分為兩個部分，第一組實驗是MDEBBO算法與經典、新型智能優化算法的對比，第二部分選取部分學者先進的改進算法與MDEBBO進行比較。

3.1 測試函數與參數設置

在這12個基準函數中，函數Sphere、Sumsquare、Rosenbrock、Schwefel1.2、Schwefel2.22屬于單峰函數，Step是具有一個最小值且不連續階梯函數，Quartic是一個噪聲四次函數，多用來測試收斂速度和尋優精度；Ackley、Griewank、Penalty1、Penalty2、Levy是多模函數，其局部最小值隨著維數的增加呈指數增長，多被用來測試算法避免局部最優解的能力以及全局探索性能。這12個基準函數已經被廣泛用于不同的優化算法中以評估算法的性能。

在實驗部分，本文所采用的硬件配置：CPU為2.0 GHz，內存24 GB的筆記本，實驗所用操作系統為Windows10，編程語言采用MATLAB R2018a。本文中對算法的參數設置為：種群大小NP為50，維數N為30，棲息地最大變異概率Pmax為0.005，棲息地修改概率Pmod為1，最大遷入率I和最大遷出率E都為1，最大迭代次數Gmax為1000，精英個體保留數為2。ACO中信息啟發式因子α為1，期望啟發式因子β為5，局部信息素蒸發系數ρ為0.5，全局信息素蒸發系數q為0.9，CS中隨機調整發現概率pa為0.25。BBO/DEs選取的是DE/best/2策略，其中F取值為0.5。因為算法運行存在一些不可控的隨機因素，所以本文將每個函數獨立運行30次來避免誤差。本文選取最優值、平均值以及標準差3項來評估算法的性能。其中，最優值是算法找到的最優解；平均值是算法所有解的平均值，能夠顯示算法的優化精度；標準差能夠體現算法的魯棒性。

3.2 與不同智能優化算法的比較

第一部分選取經典算法ACO(蟻群算法)、新型智能優化算法CS(布谷鳥搜索算法)、基本BBO和MDEBBO進行比較，結果見表1。其中，字體加粗的是算法結果中對應的最優項。根據機器對雙精度實數的取值，精度大于 1e-15 的結果即被認為零[14]。由表1可以看出，MDEBBO算法在單峰、多峰函數上與其它智能優化算法相比，無論是最優值、平均值還是標準差都有很明顯的數量級優勢，且大部分函數最優解達到了函數的理論全局最優解0。對于函數Step、Schwefel2.22和Ackley，MDEBBO算法達到了最優值、平均值和標準差同時都為0，這說明遷移算子中加入微擾動因子，變異算子引入自適應混合變異思想可以提高算法的尋優精度。算法初始化具有隨機性且標準差近乎為0，說明多次尋優結果的誤差很小，即MDEBBO算法的探索能力優勢明顯，不易陷入局部最優。同時，4種算法在不同基準函數上的收斂曲線如圖3所示，因為篇幅問題，本文只展示了其中6個不同類型的測試函數收斂圖。收斂圖中的橫坐標為迭代周期數，縱坐標為最優適應度函數值。由圖3中可以看出，其它智能優化算法在迭代后期中最優值下降速度明顯緩慢，且迭代次數為200的時候就有陷入局部最優的趨勢。而MDEBBO算法在迭代前期曲線近乎垂直，顯示了其更強的勘探能力；且在迭代次數為100左右的時候就幾乎達到了全局最優，展現了其優異的開發能力。綜上所述，MDEBBO算法在收斂精度和收斂速度上都較不同智能優化算法有很明顯的優勢，在滿足卓越開發能力的同時也避免了算法陷入局部最優。

3.3 與不同改進BBO算法的比較

為了進一步驗證MDEBBO算法的尋優能力，第二組實驗是MDEBBO算法與部分學者先進改進算法的對比，對比算法分別為基于雙曲正切變型遷移率模型的BBO(MBBO)[6]、基于Fibonacci迭代的差分變異BBO(IDEBBO)[7]、基于自適應遷移和差分變異的BBO(SDBBO)[8]，實驗結見表2。由表2可看出，MDEBBO算法在最優值、平均值與標準差3個評價標準上都較其它3個改進算法有明顯的優勢。對于全域最小值附近平緩難以尋求最優值的病態二次函數Rosenbrock，MDEBBO算法的最優值結果比其它算法的求解精度有了一定程度的提升。對于具有一個最小值且不連續階梯函數Step來說，改進的IDEBBO算法的最優值為0，但MDEBBO的平均值和標準差均為0，這說明MDEBBO算法通過微擾動因子、變異算子使算法的收斂精度和穩定性有大幅提升。多模函數常被用來測試算法避免局部最優解的能力，由表2可看出MDEBBO算法在多模函數上的尋優能力遠遠優于其它3個改進算法。同時，4種算法在不同基準函數上的收斂曲線如圖4所示，MDEBBO算法的尋優能力在同類競爭力很強的先進改進算法面前也絲毫不遜色。

表1 不同優化算法對12個函數的測試結果

圖3 不同優化算法對6個函數的收斂曲線

表2 不同改進BBO算法對12個函數的測試結果

圖4 不同改進BBO算法對6個函數的收斂曲線

3.4 算法復雜度分析

本節主要針對MDEBBO算法與傳統BBO算法的計算復雜度進行理論分析。在傳統BBO的遷移算子中，滿足遷入率條件后需通過輪盤賭方式選擇遷出棲息地，此時需要平均計算 (1+N)/2次。而在MDEBBO算法中，遷移算子添加了一個約束條件，若滿足約束條件則不需要通過輪盤賭選擇棲息地，而是通過DE的變異策略來進行遷移操作。若不滿足約束條件則進行輪盤賭，此時平均計算次數小于等于 (1+N)/2次。MDEBBO算法的差分遷移算子中雖然多了一些判斷步驟，但在數量級上并沒有增加額外的復雜度。MDEBBO算法與傳統BBO算法的變異算子的循環步驟一樣，則計算復雜度相同。綜上分析可得，MDEBBO算法不僅在收斂速度、精度提高的同時，計算復雜度并沒有增加。

4 Richards模型的參數估計

Richards模型是一個包含4個參數的非線性回歸方程，該模型通過時間變化來描述生物的生長過程[15]。Richards模型的生長方程如式(15)所示，其用來描述生長過程中的增長量

(15)

其中，yt表示谷氨酸菌體在t時刻的生長濃度值，α表示該菌體生長濃度的飽和值，β表示初始生長濃度的值，γ表示生長速率參數，t表示時間間隔，δ描述谷氨酸菌體的生長曲線。通過MDEBBO算法對Richards模型進行參數估計，其中每個棲息地代表一個解，棲息地維度為4，即α、β、γ、δ這4個參數值。非線性參數估計的解決過程是在模型給定后，通過實際觀測數據來估計各個參數的值。MDEBBO算法中的適應度函數用偏方平方和表示，如式(16)所示。實際觀測值yi與預測值的差值越小說明擬合出來的函數越接近實際的值，即預測的參數越準確[16]

(16)

以預測谷氨酸菌體生長濃度為例，使用MEDBBO算法通過仿真實驗，計算出Richards模型中的參數α、β、γ、δ分別為0.8985、5.2206、0.6309、3.4197。由MDEBBO算法估算不同時刻谷氨酸菌體生長濃度與實際觀測值之間的差值的平方和為0.009 091 1，這個值越小代表算法擬合的生長曲線越好，參數估計得越有效。圖5表示的為由MDEBBO算法對Richards模型擬合的谷氨酸菌體生長曲線。從圖中可以看出隨著時間的增長，擬合的曲線與實際觀測值之間的誤差越來越小，表明了MDEBBO算法對參數估計的有效性和適用性。

圖5 擬合的谷氨酸菌體生長曲線

為驗證MDEBBO算法對Richards模型參數估計的適用性，本文選取多個算法和MDEBBO算法進行對比。將各個算法分別獲取的谷氨酸菌體生長濃度預測值與實際觀測值比較分析，實驗結果見表3。其中對比算法包括粒子群優化算法(PSO)、遺傳算法(GA)和變步長果蠅優化算法(VS-FOA)算法。表3中實際觀測值和對比算法的預測數據均來自文獻[15]。

表3 谷氨酸菌體生長濃度的實際觀測值與不同算法的預測值

同時，本文將均方根誤差(RMSE)、決定系數(R2)以及平均誤差(MAE)作為評價指標來驗證MDEBBO算法的有效性，其分別如式(17)～式(19)[17]所示

(17)

(18)

(19)

表4 不同算法基于不同評價指標的對比

由表4可看出MEDBBO在3個評價指標下的值更具優勢，其預測精度、擬合優度、預測質量都較對比算法更高。這說明采用MDEBBO對Richards模型參數估計來預測谷氨酸菌體生長濃度更具有效性。

5 結束語

生物地理學優化算法作為一種新型的群智能優化算法，其全局搜索能力與避免早熟的性能方面仍有很大的進步空間。本文提出的MDEBBO算法通過在遷移階段將差分變異算子和自適應的局部擾動操作結合形成新的遷移算子，增強了算法的全局探索能力；且通過高斯變異和柯西變異的混合變異算子來提高算法的開發性能和避免算法陷入局部最優解的能力。實驗結果表明，MDEBBO較對比算法在收斂速度、尋優精度和穩定性3個方面都有很大程度上的提升。同時，通過仿真驗證MDEBBO算法對Richards模型參數估計的適用性，對非線性模型的參數估計提供了一種新的方法。在接下來的工作中，重點是將BBO算法應用到對地觀測衛星任務調度實際工程領域中并通過改進該算法來解決其全局優化問題。