基于種群個體數自適應的多尺度量子諧振子優化算法

焦育威 王 鵬

伴隨人工智能的發展,大量的啟發式算法相繼提出,如遺傳算法、粒子群優化算法等,這些算法在優化領域中有著出色的表現.隨著優化問題復雜程度的加深,這些算法的弊端開始顯現,如算法“早熟”或種群多樣性喪失導致精度下降、計算成本上升等問題,為解決這些問題近年來大量的改進策略被提出,這些改進策略可以劃分為兩大類.

第1 類改進方案屬于經典力學范疇,如對算法的參數進行調節,文獻[1-3] 對粒子群優化算法(Particle swarm optimization,PSO) 的慣性系數進行研究,分析慣性系數對全局搜索和局部搜索的影響,以此提高 PSO 的性能.或是融合其他算法的某種機制以提高自身性能,如文獻[4]將粒子群優化算法的局部快速收斂性與模擬退火算法的全局收斂性結合,提出一種協同進化的算法,新算法在全局收斂能力和收斂速度上都有所提高;或是以算法迭代過程中的種群為入手點,改變種群的更新方式或將單一種群進行分組或劃分,使其多種群化,如基于生物學上物種劃分的小生境遺傳算法[5].2010 年Luo 等[6]基于人工智能體的原始結構,將人工智能體劃分為幾個獨立的子群體,不同的子群體中進行通信,這一方法被證實有效.同時也有學者基于拓撲結構或層次結構對 PSO 中粒子進行種群劃分[7-9],以此平衡算法的勘探能力和開采能力[10].

第2 類改進方案是在原算法基礎上引入量子機制,從量子物理與優化問題在概率意義上的相似性出發,通過引入量子機制提高算法的性能,在優化算法中應用較為廣泛的量子機制可分為如下幾種:

1) 量子波動的引入.1994 年 Finnila 等[11]在模擬退火算法(Simulated annealing,SA)中引入量子波動,利用系統中逐級遞減的量子波動尋找經典系統的基態,提出量子退火算法(Quantum annealing,QA).量子波動的引入將函數優化問題轉化為求解系統基態的問題,目標函數f(x) 則對應量子系統中的勢阱,伴隨量子系統的演化,粒子有一定的概率發生隧道效應,在量子系統中粒子能通過勢壘貫穿出現在經典禁區(可映射為最優解區間),這稱為隧道效應.隧道效應增強了算法跳出局部最優的能力,但量子退火算法并未對隧道效應的發生概率進行主觀控制.

2) 勢阱等效.通過薛定諤方程建立采樣函數與目標函數f(x) 之間的關系,文獻[12-13]以 Delta勢阱下的波函數模方引導粒子位置更新,提出量子粒子群優化算法(Quantum PSO,QPSO),將粒子尋優空間等效為勢阱,粒子的搜索行為受勢阱約束,全局最優位置映射為勢能最低點位置.王鵬等[14]通過模擬量子諧振子從高能態向基態收斂過程,將目標函數以諧振子勢進行等效,提出用于解決高維函數優化問題的多尺度量子諧振子優化算法(Multiscale quantum harmonic oscillator algorithm,MQHOA),并于 2019 年通過模擬量子系統下自由粒子的運動過程解決函數優化問題,提出了基于自由粒子模型的優化算法框架[15].通過勢阱等效,采用某一勢阱下的波函數對粒子搜索過程進行概率描述,較經典力學的種群確定性演化過程而言該體系下算法的隨機性更強.

3) 疊加態原理的應用.首先是 Narayanan 等[16]在 1996 年通過多種群的方式將量子多宇宙的概念引入遺傳算法;隨后,Han 等[17]將量子編碼和量子旋轉門應用在種群的編碼和更新過程中,提出量子遺傳算法,疊加態原理在遺傳算法中的應用提升了算法的尋優能力,采用量子編碼的染色體可表達多個態的疊加,通過量子手段增加了種群的多樣性,類似的改進工作有量子蟻群算法[18];同時文獻[19]采用量子編碼的方式對差分進化算法進行改進,并在電力系統的優化問題中取得了不錯的效果,將量子疊加態的概率特性以量子編碼的形式應用在已有的算法改進工作中已被證明是有效的.

無論是量子退火算法中以遞減的量子波動尋找量子系統基態,還是通過勢阱等效的方法使粒子在量子系統下尋找勢能最低點,在算法的迭代過程中,粒子都有一定的概率發生量子隧道效應,發生隧道效應的概率在函數優化問題中可衡量采樣點從局部最優跳到全局最優的能力.隧道效應的發生是量子系統特有的現象,目前應用量子機制改進算法的方案中雖然涉及隧道效應的發生,但未對隧道效應發生概率有效控制以及隧道效應概率大小對算法性能影響進行研究.勢壘與粒子簇對量子隧道效應的發生概率具有較大的影響,因此在算法迭代過程中對這兩個變量信息的獲取與利用變得尤為重要.本文以經典力學中多種群搜索方式為基礎,將子種群的個體數映射為某一粒子簇中粒子的數量,將子種群間分布的相對位置用于評估當前粒子簇(采樣點)到函數最優解間的勢壘寬度,將經典力學下的多種群搜索策略賦予量子意義,通過動態調節迭代過程中子種群規模而增加種群整體量子隧道效應的發生概率,達到增強算法性能的目的,從控制隧道效應發生概率的角度為算法改進提供一種新思路.

本文結構安排如下.第1 節闡述優化算法的隧道效應,并建立多種群搜索策略與隧道效應發生概率的關系;第2 節和第3 節以 MQHOA 為例,說明其優化過程背后的物理背景,同時從發生隧道效應的角度分析該算法的不足之處,證明通過多種群手段增加隧道效應發生概率的必要性,然后基于蒙特卡洛評估的思想提出子種群規模自適應的策略;第4 節給出改進后具體的算法MQHOA-d (Multiscale quantum harmonic oscillator algorithm based on dynamic subpopulation),并以雙阱函數為例設計實驗跟蹤算法在尋優過程中發生量子隧道效應點的數量變化情況,同時說明隧道效應對算法尋優的影響;第5 節設計實驗并對改進后的算法(MQHOA-d) 從尋優誤差、計算資源消耗等方面進行分析.

1 多種群策略與隧道效應

在經典力學體系下,將優化算法的采樣方式由單一種群采樣改進為多種群采樣,從某種程度上增加了算法迭代過程中解的多樣性,由此提高原算法尋優能力.但該方法仍存在弊端,對于將原算法單一種群采樣變為多種群采樣這一類改進方案,普遍存在一個共性,即當算法對子種群的劃分完成后,其子種群的數目是固定的,隨著迭代進行,子種群的個體數并不會動態變化.

對于類似多種群化的改進策略,除種群間的協作或學習等方式對算法性能的提高外,采樣過程中子種群規模的變化對算法性能也具有影響,大量的算法改進專注于種群間的行為方式,而忽略子種群規模自身對性能的影響.如果從量子力學的角度看待多種群策略,可將以多種群的形式代表的整體采樣點中子種群的數量映射為粒子簇的數量,某一子種群規模(個體數) 映射為粒子簇中粒子的數量.當前子種群到最優解區間的距離用于勢壘的評估,目標函數f(x) 用于勢能評估,詳細對應關系見表1.由此可將函數優化過程看作種群尋找勢能最低點位置的過程,在量子體系下粒子向勢能最低點的搜索過程中,粒子存在穿越勢壘到達勢能最低點的可能性,即采樣點從局部最優跳到全局最優區間,這種現象稱為隧道效應,也是量子系統特有的現象.而隧道效應的發生概率受子種群個體數與當前子種群到最優解區間的距離(勢壘) 所影響,通過勢壘信息自適應地調節子種群規模,可改變總體隧道效應的發生概率,達到增強算法尋優能力的目的,同時彌補了經典力學下多種群改進方案忽略子種群規模自身影響算法性能的弊端.

表1 量子理論與優化算法的對應關系Table 1 The relationship between quantum theory and optimization algorithm

通過定義優化問題的波函數,可將最優區域[x0-ε,x0+ε] (x0為最優解,ε為求解誤差) 發生隧道效應的概率通過波函數進行計算,此處記為Pi.則通過量子隧道效應在下一次迭代中出現在最優區域的個體數為×Pi,其中k為種群數,dyni為第i個子種群的規模(個體數).在尋優過程中,發生隧道效應的粒子將出現在最優解區域,由于多種群采樣方式具有一定的信息交流,發生量子隧道效應的個體可引導所屬的子種群向最優解方向遷移,以此達到由各子種群組成的種群整體向最優區域收斂.

這里不妨假設一種情況,某一子種群i的規模(個體數) 為dyni,且當前區域內發生隧道效應的概率為Pi,則該子種群在下一次迭代過程中通過隧道效應達到全局最優區間的個體數期望為Pi×dyni.無限制的增大子種群規模dyni可以增強該種群向最優解區間的遷移能力,但該操作會改變種群采樣率.文獻[20]指出增加采樣點(種群規模)的數量會提高采樣率,而采樣率越小算法效率越高,出于上述考慮,基于蒙特卡洛估計的思想[21],提出可以根據種群適應度動態調節子種群個體數的策略.該策略在不改變算法采樣率的基礎上改變了迭代過程中種群發生隧道效應的概率,由此提高了算法性能,同時這一策略可以廣義地應用在各種優化算法的改進工作中.

第2 節將以 MQHOA 為例,從控制隧道效應的角度對該算法進行改進.選擇 MQHOA 作為改進案例的原因在于該算法是基于量子理論提出的優化算法,若選擇經典力學體系下的優化算法作為目標算法,量子機制的引入容易從多個維度增強算法隨機性(如種群在量子體系下的疊加性),雖然能增強算法性能達到改進的目的,但無法準確地驗證通過子種群規模自適應改變隧道效應的發生對提高算法尋優能力的有效性.而 MQHOA 以薛定諤方程為基礎,建立了優化問題與量子系統之間的關系,將優化問題轉化為求解粒子在約束下的基態波函數.MQHOA在迭代過程中會自發地產生隧道效應,以 MQHOA為目標算法,通過動態調節子種群個體數的方式改變其迭代過程中發生隧道效應的粒子數,可以更好地驗證通過隧道效應改進優化算法的有效性和可行性.

2 標準 MQHOA

2.1 MQHOA 物理模型

量子力學中波函數可以描述粒子在某位置出現的概率,通過文獻[22]可知,薛定諤方程可以描述量子系統,不隨時間變化的薛定諤方程表示為

通過薛定諤方程的變換,可將求解目標函數最小值問題轉化為求解諧振子勢阱下的基態波函數的問題.經過近似操作后,波函數在基態的概率密度函數為

由式(3) 可知,目標函數的最優解的概率分布為高斯分布.MQHOA 算法以高斯采樣為基本操作,k個高斯采樣點作為粒子,利用k個高斯采樣疊加作為波函數描述當前最優解的概率分布.在 MQHOA 的收斂過程中,諧振子波函數從高能態逐步變向基態,從高能態多個正態分布的概率疊加,逐步收斂到基態單一正態分布的穩定狀態,此時算法獲得最優解.

2.2 MQHOA 基本流程

MQHOA 算法流程如算法1 所示.

算法 1.MQHOA 算法

算法1 中,參數k為采樣中心數,S為第1 次采樣區間,σ0和σmin分別為初始尺度和最小尺度,其中σmin對算法解的精度具有重要影響.f(x) 為目標函數.偽代碼中從內到外的3 層循環分別對應MQHOA 的3 種物理過程,分別是能級穩定、能級下降、尺度下降.第1 層循環即對應能級穩定過程,在此過程中以k個采樣點為中心分別按正態分布N(xi,)在每一個維度進行采樣,以此生成新的采樣點.計算當前采樣點的目標函數值f() 是否小于當前采樣中心的目標函數值f(xi),滿足則進行替換,不滿足則暫不替換當前采樣中心.此外,如采樣過程中受優化空間限制,則使用映射操作將越界點映射到優化空間內,本文不討論上述映射機制.

當k個位置的采樣中心全部完成采樣后,計算x0～xk的方差σk,然后與迭代前的方差做差,計算出迭代前后方程的變化量 Δσk.如果Δσk ＜σs則判定算法能級穩定過程結束,進行能級下降操作,以k個采樣點的均值替換當前目標函數值最差的采樣點,即xworstxmean.當算法基態判據滿足σk ＜σs時,可認為算法達到了尺度σs下的能量基態,將進行尺度下降過程,這個過程通過尺度減半的操作,即σs/2 降低當前尺度,這一過程決定了搜索精度,隨著尺度下降這一步驟的進行,算法從全局搜索逐步過渡到局部搜索.

3 標準 MQHOA 存在的問題及解決方案

3.1 標準 MQHOA 存在的問題

原算法(MQHOA)忽視了不同采樣中心在整個迭代過程中的適應度值是變化的,不同采樣中心發生隧道效應概率不同,因此在每個中心位置采取相同的采樣行為是不合理的.在解空間中,當部分采樣點陷入局部最優區域時應增強該區域種群通過隧道效應跳出該區域的能力,由于原算法存在上述缺陷,由此對部分函數尋優時容易發生算法“早熟”現象.

MQHOA 搜索過程為k個正態概率分布不斷向最優解位置迭代的過程,即算法的波函數為k個正態分布的疊加,可表示為

以一維函數為例,xi(i1,2,3,···,k) 為當前各正態分布期望,ε為誤差上限,存在k個正態分布Ni(xi,).對k個正態分布約束下的任一采樣點xi,在當前采樣中能找到最優解的概率為

此時以該正態分布進行采樣的種群發生量子隧道效應采樣點的數量計算式為

其中,x0為最優解位置,dyni為以當前正態分布采樣的子種群個體數.MQHOA 為單一種群采樣,可看作k個規模為1 (采樣點) 的子種群集合所構成的種群整體(單一種群),因此對 MQHOA 而言,dyni1,i1,2,3,···,k.由于采樣中心xi的不同,進行 Ni(xi,) 采樣時,采樣點落在區間[x0-ε,x0+ε]的概率具有差異性(在高維函數解空間中任意兩采樣中心x1,x2相等的概率極低,因此不考慮這種情況).圖1 所示的2 條高斯曲線分別代表采樣中心x12,x28 的概率密度分布情況,X軸上的黑點代表最優解x0的位置,陰影部分為采樣點落在區間 [x0-ε,x0+ε] 的概率P(xi).以x2 為中心的高斯曲線下陰影部分面積大于以x8 為中心的高斯曲線下的陰影部分面積,即在這次采樣過程中以x2 為中心的采樣區域能獲得最優解位置的概率大于以x8 為采樣中心的區域.因此,采樣成功的概率會隨采樣中心的變化而產生差異.

圖1 不同采樣中心概率密度曲線Fig.1 Probability density curves of different sampling centers

本文將采樣中心在當前迭代中能找到最優解的概率稱作適應度值,適應度的大小以采樣點與最優解間的勢壘寬度評估.在采樣點所構成的集合中,采樣點xi通過隧道效應跳到最優解區域的概率最大,本文僅需比較k個采樣中心適應度的排名,不需要進行精確計算,且當前適應度值最大的采樣中心xbest應滿

式中,xr為除具有最大適應度采樣點xbest外的任一采樣點.原算法忽略了對適應度不同的采樣中心進行區別采樣,導致種群在適應度較小的區域采樣時,不能準確獲取最優解位置信息,即能級穩定過程緩慢且易陷入局部最優解,同時影響收斂速度.實驗表明,當目標函數的維度較高時,在規定的迭代次數中,MQHOA 的求解成功率會大幅度降低.當算法隨著迭代次數的增加,從全局搜索過渡到局部搜索后,這種采樣策略在搜索能力上的不足開始體現,尋優誤差較大.

3.2 基于蒙特卡洛估計的種群數分配

函數的優化問題從某方面而言是具有黑盒特性的,在優化問題中目標函數的表達式f(x) 是不確定的,在迭代過程中,每一個子種群發生mi次測量,mi為當前子種群的個體數(規模),實際上可以將采樣映射為對函數f(x) 積分的評估.通過第2 節可知,增加以當前中心xi的采樣次數mi可以增大找到最優解概率,本節內容具體分析對不同mi之間是否應該具有差異.基于蒙特卡洛思想[23],可知采樣越多,越近似于最優解.對于函數優化,在最優解xo所在誤差范圍 [x0-ε,x0+ε] 內,應滿足

由于MQHOA 采樣過程中尺度σs是不斷遞減的,對于k個子種群搜索過程而言,實際上是評估k個子種群在當前尺度下積分最小的問題,表達為

由于波函數概率密度函數為正態分布,函數采樣區間有限,這里做近似處理

式(9)是以xi(i1,···,k) 為中心的子種群對當前區域的評估蒙特卡洛算子,mi(i1,···,k)為該子種群的規模(個體數),P(xi)(i1,···,k)為當前第i個子種群的概率密度函數,其中,N為采樣次數,它與mi相對應,f(x)為目標函數,P(x)為x的概率密度函數,記f(xi)/P(xi)g(xi),據大數定律可知

由此可以推出區間內積分在N次重復采樣下的表現形式為

對于不同中心Xi而言,定義域內的采樣概率密度函數具有差異,而固定區間內對函數的積分為定值,此時應分析概率密度函數對積分評估的影響.我們知道所有采樣中心點服從正態分布的方差是相同的,不同的是期望,而期望值為采樣中心的坐標信息.此時計算方差為

記g(xi)方差為σ2,此時g(xi) 相互獨立,因此蒙特卡洛方差與σ2的關系為

4 基于多種群的MQHOA-d

4.1 MQHOA-d 流程

k個種群以當前種群中心的波函數為基準進行采樣,每個子種群的個體數mi(i1,···,k) 的更新是根據當前中心到適應度最高種群中心的歐氏距離決定的,更新式為

其中,n為維度大小,l為當前維度,Xi為當前種群中心,Xbest為適應度最高的采樣中心,依據歐氏距離的種群規模更新式為

式 (14)通過各子種群中心Xi到當前最優子種群中心Xbest的歐氏距離評估各子種群所處解空間的適應度,距當前最優子種群中心距離越遠,則判斷該子種群適應度越差,由此通過式 (15) 更新該子種群規模(個體數).該操作在量子力學上的意義在于通過評估粒子簇與最優解區間的勢壘寬度(式(14)對子種群中心歐氏距離的計算) 調節粒子數量(式(15)對子種群規模的調節),從而改變該區域發生量子隧道效應的概率.

能級穩定判據對算法在當前尺度下種群進行充分搜索有重要影響,當k個種群遷移速度過快,達到能級穩定判據標準時,會進行能級下降,種群分布聚集,由此容易導致算法早熟,全局搜索不充分.因此,需要在每輪迭代過程中對當前尺度下的種群空間進行限定,即實現種群遷移步長控制.MQHOAd 在迭代過程中,首先初始化k個種群中心,以當前尺度σs2生成n維協方差矩陣,具體計算為

其中,cov(g)表示第g代的協方差矩陣,在當前種群中心Xi處按照正態分布 N(Xi,cov(g))生成個采樣點,完成當前種群個體初始化.當種群方差滿足能級穩定判據時 (Δσk ≤σs),MQHOA-d 進行能級下降操作,需要將適應度最差的種群所有個體淘汰.取剩下k個種群中心位置的均值坐標,初始化一個新的種群中心m(g),具體為

以Xi為中心通過式 (14)和式(15) 生成種群個體,用新種群替換當前最差種群,完成替換過程,實現能級下降.當所有的個體滿足能級下降判據,此時在當前尺度的搜索過程完成,下一步在更小的尺度進行搜索,進行尺度下降操作.MQHOAd 的尺度下降步驟與MQHOA 相同,仍為尺度減半操作.當滿足基態判據時,輸出當前最優解,尋優結束.

圖2 表示k為30、λ為1 000 時,MQHOA-d所生成的種群在 Ellipsoidal 函數上種群分布示意圖.函數最優解位置為(1,2,3,···,n),粗點代表當前種群中心Xi,細點為采樣點(當前子種群個體).本文中子種群是指以當前位置Xi為均值,根據正態分布 N(Xi,cov(g))生成的個采樣點,即屬于同一個子種群中的個體 (采樣點) 變更滿足N(Xi,cov(g)).由圖2 可見,依據式 (14)和式(15)生成的種群特點如下: 適應度差的區域分布子種群規模大,即個體數多;適應度優的區域子種群分布規模小.這一操作的目的是降低第2.1 節中對適應度差的子種群空間的評估誤差,以此保證整體評估的穩定性.

圖2 MQHOA-d 所生成種群在Ellipsoidal 函數二維空間分布示意圖Fig.2 Schematic diagram of spatial distribution of subpopulations generated by MQHOA-d in Ellipsoidal function of 2D

MQHOA-d 偽代碼見算法2,k為種群的數量,S為第1 次采樣區間,σs和σmin分別為當前尺度和最小尺度,其中,σs初始化為σs=UB-LB,初始化的k個采樣中心.(g1,2,···,k) 代表每一個種群的規模 (個體數).

算法 2.MQHOA-d 算法

4.2 量子隧道效應跟蹤實驗

多種群化可以增加種群發生量子隧穿的概率,為了更好地分析該策略與粒子發生量子隧道效應的關系,用 MQHOA/MQHOA-d 做對照實驗.實驗組1 中 MQHOA 采樣方式為單一種群采樣,種群數為1,個體數通常為30,這里為控制參數一致,將個體數設置為200;實驗組2 中 MQHOA-d 采樣方式為多種群采樣,設置子種群數為30,個體總數為200.在迭代過程中可以根據當前波函數的概率密度計算發生量子隧道效應點的數量,即

圖3(a) 中實線為實驗所用函數的一維示意圖,其中加粗部分為最優區間,圖3(b)為子種群中部分采樣個體發生量子隧道效應,貫穿勢壘跳出局部最優示意圖.由式(18)可知在下一次迭代過程中,該子種群發生隧道效應的個體數量與當前子種群規模正相關,式(15) 通過適應度的差異性調節各子種群的規模,由此增加適應度差的子種群個體數,以此提高當前能級下該子種群量子隧穿到最優解區域的概率,增強該子種群向最優解區域的遷移能力.隧道效應的發生是粒子數量與勢壘(本文映射為適應度) 共同影響的結果,固定子種群規模的搜索策略忽略了子種群個體數量對隧道效應的影響,因此,基于子種群自適應策略改進的 MQHOA-d 在相同的環境下較固定種群規模搜索策略的算法發生隧道效應的概率更大,由此可以提升算法搜索能力.

圖3 一維雙阱函數圖與隧道效應示意圖Fig.3 One dimensional double well function diagram and tunnel effect diagram

實驗的目的之一是計算種群中通過量子隧道效應出現在最優區間的點的數量,虛線為隨迭代次數增加采樣概率密度函數尺度遞減示意圖,當迭代次數增加到一定數量時,概率密度函數中σs收斂于0,此時發生量子隧道效應的概率也趨近于0,因此會發生圖4 中隨迭代次數增加,量子隧道效應的點的數量先增加后遞減甚至趨于0.圖4 中三角形折線和圓形折線分別代表 MQHOA-d和MQHOA種群中發生量子隧道效應的點數量與迭代次數的關系圖,MQHOA-d 與 MQHOA 的迭代過程中,在相同的迭代次數下,MQHOA-d 發生量子隧道效應點的數量要大于 MQHOA,實驗設置的輸出條件為迭代次數達到10 000.通過圖4 可知,MQHOA-d 在迭代次數達到約2 500 時尺度σs收斂于0,即種群已完成搜索,而 MQHOA 需要迭代約9 500 次才完成搜索.從發生量子隧道效應點數量與迭代次數的關系圖可知,快速地達到隧道效應點數量最大化有利于縮減算法搜索時間,提高種群收斂速度,這也是對第3.1 節中通過歐氏距離評估適應度而動態調節子種群規模使整體發生量子隧道效應最大化的有效性驗證.

圖4 MQHOA-d 與MQHOA 發生隧道效應點數量與迭代次數關系圖Fig.4 The relationship between the number of tunneling points and the number of iterations between MQHOA-d and MQHOA

發生量子隧道效應的點數量越多可映射為當前種群分布尋找到最優解區域的概率越大.MQHOA-d在迭代次數較小時即可達到整體種群發生隧道效應的最大狀態,此時尺度的縮減程度低,可以說明MQHOA-d 在全局搜索時即可快速定位全局最優區域,面對復雜函數或多峰函數時更具有優勢;而MQHOA 在整體種群發生隧道效應達到最大狀態時需要的迭代次數較大,此時尺度σs縮減程度高,在大尺度下對全局的評估不如 MQHOA-d,由此可以推測當測試函數為結構復雜的函數時,MQHOA容易陷入局部最優區域,且 MQHOA 的尋優能力不如 MQHOA-d.從蒙特卡洛估計積分的角度上看,合理降低適應度差的子種群采樣方差可以提高對當前種群整體分布的評估準確度,有利于解空間的全局評估.

5 實驗及分析

為了全面分析 MQHOA-d 的性能,選取煙花算法 (Fireworks algorithm,FWA)[24-25]的衍生EFWA (Enhanced FWA)[26]、AFWA (Adaptive FWA)[27]以及 MQHOA、SPSO2011[28]、QPSO[29-30]作為對照組.

所有的算法在 CEC2013 標準函數測試集[31]上進行,每組實驗重復51 次,單次運算的最大迭代次數為10 000 ×d,其中d為維度,實驗維度設置為30.測試集一共包括28 個測試函數,根據函數結構不同可劃分為單峰函數 (f1～f5)、多峰函數(f6～f20)和復合函數 (f21～f28).所有函數的定義域為[-100,100],后續對實驗結果的誤差進行分析.實驗環境為Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz,2.21 GHz,16.0 GB,程序采用 MATLAB2018a 編寫.其中,MQHOA、MQHOA-d 的種群數設置為30,MQHOA-d 的個體總數λ設置為200;QPSO和SPSO2011 粒子總數設置為30.具體實驗數據見表1(EFWA、AFWA 實驗數據來自于文獻[32]).表1中 MeanErr 表示誤差均值,Std 代表解的標準差,最優的實驗結果在表2 中用粗體標出.

5.1 實驗數據分析

為更直觀地分析數據,將6 組實驗數據的誤差均值進行等級劃分(1 級～ 6 級),均值誤差最小的為1 級,最大的為6 級,同時計算對照組的平均等級并繪制成圖5.

圖5 CEC2013 測試集上各算法誤差等級圖(維度為30,重復51 次)Fig.5 Rank sum test results of all algorithms on CEC2013 test set (repeated 51 times in the 30-dimensional case)

與原 MQHOA 算法相比較,通過表2和圖6可以看出,MQHOA-d 在單峰函數上 (除f3) 表現較差,這是由于 MQHOA-d 每一個高斯采樣空間中產生的個采樣點并不能實現向最優解空間的快速遷移,一次采樣所需要的計算次數為而這一操作對應算法的能級穩定過程: MQHOA 僅在每一個高斯采樣空間中生成一個采樣點,即一次采樣僅需要計算k次.而結構簡單的單峰函數并無過多局部最優解,多種群搜索的策略雖可以提高跳出局部最優解的概率,但針對這類型的函數,需要的是整個種群快速向最優解區域收斂,即搜索空間從大尺度向小尺度的過渡,這一過程對應 MQHOA/MQHOA-d 的能級下降過程,算法需要盡快從全局過渡到局部.而 MQHOA-d 的多種群生成決定了它將大量的計算資源消耗在能級穩定的過程中,由于實驗計算次數存在上限,因此在有限的計算中,MQHOA-d 進行能級下降的次數比 MQHOA 少,即尺度縮減程度不如 MQHOA.

從多峰函數和復合函數上看,MQHOA-d 的求解誤差要遠小于原 MQHOA,這也驗證了第4.2 節中根據量子隧道效應實驗數據對算法性能的推測是可靠的,MQHOA-d 在大尺度的狀態下盡可能多地進行搜索,全局搜索更加全面.此策略在適應度低的區域生成更多的采樣個體,對整個種群的計算能力進行合理分配,提高了種群多樣性,增加跳出局部最優解概率.因此,MQHOA-d 對結構復雜的多峰函數和復合函數更具有針對性.相比于 SPSO2011、QPSO、EFWA、AFWA 等算法,通過圖6 可知,MQHOA-d 在單峰函數下處于劣勢,但多峰函數和復合函數的平均等級最小,且總體平均誤差等級最小.綜合28 個測試函數分析,MQHOA-d 的平均誤差等級值最小,僅為1.93.從表2 中可知,MQHOAd 的標準差在多峰函數和復合函數中較小,具有較強的魯棒性.通過每組實驗重復51 次的數據繪制箱體圖(具體見圖6).這些箱體圖清晰地顯示出MQHOA-d 在多峰函數和復合函數中的準確性和穩定性占據主導地位.

MQHOA-d 性能的增強得益于搜索隨機性的增強,在量子模型下該算法的隨機性可由量子隧道效應衡量.種群個體數自適應過程擴大適應度差的子種群規模,以此增強搜索過程的隨機性,同時降低種群對目標函數的評估誤差.在量子模型下解釋為當勢壘寬度較寬時,增加粒子數量可以加大粒子簇發生隧道效應的概率,非量子模型下的優化算法通常以增加擾動項或設置變異操作的方式增強算法的隨機性.如文獻[30]中通過評估種群所處的階段實現變異策略的反饋調節,由此提出基于狀態估計反饋的策略自適應差分進化算法 SEFDE,且 SEFDE 在與 CEC2013 中性能最好的 DE 改進算法SHADE 的對比實驗中表現出較強競爭力.將表2 中MQHOA-d 的 CEC2013 測試數據與文獻[33]中SEFDE 在同等條件下的 CEC2013 測試數據進行對比,結果顯示,MQHOA-d 在16 個函數優化上優于SEFDE,該16 個函數除f2外,都為多峰函數或復合函數,與此同時,MQHOA-d 在函數f8的優化上與SEFDE 具有相似的優化效果.總體而言,MQHOAd 在多峰函數和復合函數的優化上優于 SEFDE.這也驗證了基于蒙特卡洛評估提出的動態種群策略的合理性,增大適應度差的種群規模可以提高算法搜索隨機性,有利于降低整體評估誤差.

統計51 次重復實驗的平均迭代時間,繪制成圖7.在圖7 中,方框中數字代表算法在對應測試函數中所消耗的平均時間,單位為s.將時間歸一化后進行涂色,顏色越深的區域代表消耗時間越多.

本次實驗中,算法MQHOA 在測試函數f1上的平均迭代次數為 3.49×104,小于上限D×10 000,其他實驗組的平均迭代次數均達到上限D×10 000.針對測試函數f1而言,MQHOA 所消耗時間較其他3 種算法少一個數量級,且圖6 中箱體圖顯示MQHOA 的51 次求解數據分布不穩定,可見 MQHOA 在f1測試函數上存在算法“早熟” 現象,過早地收斂在某一局部最優無法跳出,也驗證了增大種群初期發生量子隧道效應的概率更不易陷入局部最優.在測試函數f9上,MQHOA-d 所消耗時間較其他3 種算法具有極大的優勢,平均誤差僅次于 QPSO,且誤差范圍與 QPSO 在同一個數量級,當實際工業生產面對結構類似于f9函數時,可在節約計算成本的情況下使用 MQHOA-d.總體而言,MQHOA 算法消耗時間最少,但尋優精度遠不如其他3種算法;MQHOA-d 消耗時間僅次于 MQHOA,但尋優精度和魯棒性強于 MQHOA、QPSO和SPSO2011.

6 結束語

本文基于蒙特卡洛估計中不同子種群對總體評估誤差的差異性,提出動態調節子種群規模的尋優策略,將此策略應用于 MQHOA 的改進,改進后的算法(MQHOA-d)在未增加隨機化特征,或者添加擾動等多余步驟的前提下解決了原算法的“早熟”現象[34].核心在于動態調節種群規模方式增大了種群迭代初期發生量子隧道效應的概率,使其對解空間全局評估更加精準.在迭代次數相同的情況下MQHOA-d 對多峰函數和復合函數尋優誤差更小,根據適應度評估的多種群采樣方式動態地調節了種群規模,在解空間中更加合理地進行勘探和開發,同時種群利用率得到了提高.改進后的算法雖然結構上較之前相比變得復雜,但時間復雜度上仍優于幾大主流算法(SPSO2011、QPSO).目前,MQHOAd 仍有局限性,如: 1) MQHOA 系列算法的隨機性和隱含并行性尚不明確[35-36];2) MQHOA-d 種群個體行為所抽象的隨機游走模型[37-38]尚未研究.