一種求解團(tuán)簇優(yōu)化問題的混合進(jìn)化算法*

向 垚 鄭斯瑞 賴向京

（南京郵電大學(xué)先進(jìn)技術(shù)研究院 南京 210023）

1 引言

全局優(yōu)化問題廣泛存在于各個(gè)領(lǐng)域，其中通過搜索原子團(tuán)簇勢(shì)能函數(shù)的全局最小值，是計(jì)算化學(xué)領(lǐng)域中一個(gè)著名的全局優(yōu)化問題。在計(jì)算化學(xué)領(lǐng)域，原子團(tuán)簇的許多物理和化學(xué)性質(zhì)是由其幾何結(jié)構(gòu)決定的，因此確定它們的基態(tài)結(jié)構(gòu)對(duì)于理解它們的各種性質(zhì)具有重要意義［1］。由于團(tuán)簇的局部最優(yōu)結(jié)構(gòu)的數(shù)量隨團(tuán)簇規(guī)模呈指數(shù)型增長(zhǎng)，確定原子團(tuán)簇基態(tài)結(jié)構(gòu)（即團(tuán)簇的幾何優(yōu)化）是一項(xiàng)非常具有挑戰(zhàn)性的任務(wù)。事實(shí)上，原子團(tuán)簇優(yōu)化問題也已經(jīng)被證明為NP-hard問題［2］。

原子團(tuán)簇的結(jié)構(gòu)優(yōu)化問題由于其NP-hard 特性和應(yīng)用價(jià)值，受到了許多研究者的廣泛研究。在文獻(xiàn)中已有大量相關(guān)算法被提出，這些方法可以分為兩種，包括基于單解的算法和基于種群的算法。對(duì)于基于種群算法，具有代表性的算法包括隨機(jī)隧道優(yōu)化算法（Random Tunneling Algorithm，RTA）［3］，自適應(yīng)免疫優(yōu)化算法（Adaptive Immune Optimization Algorithm，AIOA）［4］，構(gòu)型空間退火算法（Conformational Space Annealing，CSA）［5］，基于種群的盆地跳躍算法（Population Basin Hopping，PBH）［6］等。對(duì)于基于單解的算法，具有代表性的算法包括basin-hopping 算法（BH）［2］，帶有表面算子和內(nèi)部算子的啟發(fā)式算法（Heuristic Algorithm with the Surface and Interior Operators，HA-SIO）［7～8］，動(dòng)態(tài)格子搜索算法（Dynamic Lattice Searching，DLS）及其變種［9～11］，和模糊全局優(yōu)化算法（Fuzzy Global Optimization，F(xiàn)GO）［12］等。

PBH算法是一種混合進(jìn)化算法，在許多全局優(yōu)化問題上表現(xiàn)出了很高的性能，例如圓形packing問題［13］。本文通過整合PBH 算法的框架和基于連通表（Connectivity Table，CT）［4］的距離函數(shù)，為Gupta 勢(shì)能函數(shù)建模的原子團(tuán)簇的結(jié)構(gòu)優(yōu)化提出了PBH算法的一個(gè)變種，并將該算法成功應(yīng)用于銀團(tuán)簇的結(jié)構(gòu)優(yōu)化。計(jì)算結(jié)果很好地表明了該算法的有效性。

2 Gupta勢(shì)能函數(shù)模型

Gupta 勢(shì)能［14］常用于給金屬團(tuán)簇建模，其勢(shì)能函數(shù)具有如下形式：

其中N為原子數(shù)目，rij為原子i與原子j之間的歐式距離，r0為團(tuán)簇兩原子之間的平衡距離，r0=1.0，A，p，q，UN為勢(shì)能函數(shù)的參數(shù)，這些參數(shù)的不同組合能表示不同類型的原子團(tuán)簇。對(duì)于銀團(tuán)簇，其相應(yīng)的參數(shù)取值如下：

3 基于連通表的種群盆地跳躍算法（PBH-CT）

3.1 PBH-CT算法的主框架

本文提出的PBH-CT 算法是一種混合進(jìn)化算法，它由一個(gè)種群初始化方法、一個(gè)局部?jī)?yōu)化方法，一種擾動(dòng)算子和一個(gè)種群更新策略組成，其偽代碼如算法1所述。

從一個(gè)隨機(jī)生成的m個(gè)個(gè)體構(gòu)成的初始種群開始，算法執(zhí)行若干次迭代，直到連續(xù)MaxNoImp次迭代無法改進(jìn)當(dāng)前最好解為止，其中MaxNoImp 是一個(gè)參數(shù)。在每個(gè)迭代，首先找到種群X中與當(dāng)前解Yk距離最近的構(gòu)型Xr，然后通過將兩個(gè)構(gòu)型之間的距離與接收罰值dcut進(jìn)行比較，來判斷當(dāng)前構(gòu)型Yk是否與種群X的解Xr相似。

若d(Yk，Xr)＞dcut，則表示兩個(gè)構(gòu)型Yk與Xr是不相似的。此時(shí)，我們比較Yk與種群X中的最差構(gòu)型Xw的目標(biāo)函數(shù)值；若E(Yk)＜E(Xw)，則用Yk替換種群X中的Xw；若d(Yk，Xr)≤dcut，則認(rèn)為兩個(gè)構(gòu)型是相似的，此時(shí)，我們比較這兩個(gè)構(gòu)型的目標(biāo)函數(shù)值；若E(Yk)＜E(Xr)，則用Yk替換種群X中的Xr。

3.2 初始種群

在基于種群的算法中，種群的多樣性對(duì)算法性能起著關(guān)鍵作用。因此，為了保證初始種群X的多樣性，本文采用了如下的初始種群策略：隨機(jī)生成m個(gè)構(gòu)型來構(gòu)成初始種群，每個(gè)個(gè)體都是在一個(gè)半徑為R=（3N/（4Π））1/3·r0的球體中隨機(jī)撒入N個(gè)原子中心點(diǎn)，從而得到一個(gè)隨機(jī)構(gòu)型，并通過局部?jī)?yōu)化方法對(duì)該構(gòu)型進(jìn)行局部?jī)?yōu)化，從而得到初始種群X。在本文算法中，局部?jī)?yōu)化方法采用的是L-BFGS［16］算法，其中L-BFGS 算法是一種高效的局部?jī)?yōu)化算法，已被廣泛應(yīng)用于可導(dǎo)函數(shù)的局部?jī)?yōu)化。

3.3 擾動(dòng)

對(duì)于團(tuán)簇優(yōu)化問題，其勢(shì)能函數(shù)具有大量的局部最優(yōu)解，且局部最優(yōu)值的個(gè)數(shù)隨著原子個(gè)數(shù)N呈指數(shù)關(guān)系增長(zhǎng)。因此，對(duì)于規(guī)模較大的算例，如何跳出這些局部陷阱變成了解決團(tuán)簇優(yōu)化問題的主要難點(diǎn)。為了避免陷入局部最優(yōu)陷阱，本文采用經(jīng)典的隨機(jī)擾動(dòng)策略。具體而言，給定一個(gè)候選解x=(x1，x2，…，x3N)，通過隨機(jī)小擾動(dòng)來生成后代解，即xi←xi+rand(-Δ，Δ) (1 ≤i≤3N)，其中rand(-Δ，Δ)表示[-Δ，Δ]中的一個(gè)隨機(jī)數(shù)，Δ 是一個(gè)參數(shù)，在本文中設(shè)置成0.7。

3.4 構(gòu)型間的距離函數(shù)

在該算法中，采用團(tuán)簇的最近鄰接數(shù)（記為nnn）和連通表（CT）［4］來定義兩個(gè)團(tuán)簇構(gòu)型的距離函數(shù)d（A，B），具體而言，對(duì)于一個(gè)團(tuán)簇來說，nnn的定義如下：

其中rij表示原子i和原子j之間的歐式距離，r0為原子間的平衡距離。

另外，CT 是一個(gè)基于團(tuán)簇構(gòu)型的拓?fù)湫畔⒌南蛄?，由團(tuán)簇中原子的最近鄰居個(gè)數(shù)生成。CT 為14 維向量，其中CT［i］（1 ≤i≤14）表示在團(tuán)簇中有i個(gè)鄰居原子的原子個(gè)數(shù)。值得注意的是，在所研究的原子團(tuán)簇中，原子的最近鄰居數(shù)均不超過14。

借助連通表（CT）和最近鄰接數(shù)（nnn），可以如下定義團(tuán)簇構(gòu)型A和團(tuán)簇構(gòu)型B之間的距離函數(shù)d（A，B）：

其中，（或）是構(gòu)型A（或B）的最近鄰接數(shù)。CT（Ai）（或CT（Bi））是構(gòu)型A（或B）中有i個(gè)鄰居原子的原子個(gè)數(shù)。

4 實(shí)驗(yàn)結(jié)果與比較

4.1 實(shí)驗(yàn)條件

本文算法通過C 語言實(shí)現(xiàn)，并在Windows 10操作系統(tǒng)、主頻為2.5GHz、內(nèi)存為8GB 的PC 機(jī)上進(jìn)行實(shí)驗(yàn)。算法參數(shù)設(shè)置如下：m=20，MaxNoImp=400。另外，由于算法的隨機(jī)性，對(duì)于n∈［10，61］范圍內(nèi)的每個(gè)銀團(tuán)簇，本文算法獨(dú)立地運(yùn)行了10次。

4.2 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)結(jié)果總結(jié)在表1 中，其中第1 列給出了團(tuán)簇中的原子個(gè)數(shù)（N），第2 列給出了文獻(xiàn)中已知的最優(yōu)解，本文所提出的算法的計(jì)算結(jié)果在隨后3 列中給出，包括算法的最優(yōu)解，達(dá)到最優(yōu)解所需要的平均計(jì)算時(shí)間和達(dá)到最好解的成功率。

表1 計(jì)算結(jié)果與當(dāng)前文獻(xiàn)中最好的結(jié)果比較，其中最好的結(jié)果來自參考文獻(xiàn)［3］和［9］

從表1 中可以看出，所提出的算法在大多數(shù)測(cè)試的算例中都可以找到當(dāng)前文獻(xiàn)中最好的結(jié)果，這表明所提出的算法在10 ≤N≤61 的范圍內(nèi)具有強(qiáng)的搜索能力。然而，該算法在4 個(gè)難例（即Ag53、Ag54、Ag58、Ag59）上未找到當(dāng)前已知的最優(yōu)結(jié)果，其計(jì)算結(jié)果在表1中用斜體表示。

對(duì)于10 ≤N≤29 的小團(tuán)簇來說，所提出的算法在平均計(jì)算時(shí)間和算法成功率上都具有較高的性能。但是隨著原子個(gè)數(shù)N的增加，算法的計(jì)算時(shí)間會(huì)顯著增加。

另外，為了直觀地表示計(jì)算結(jié)果，我們?cè)趫D1中分別給出4 個(gè)代表性團(tuán)簇構(gòu)型的最優(yōu)構(gòu)型的圖形。如圖1所示，Ag33的構(gòu)型是基于十面體的結(jié)構(gòu)，Ag39的構(gòu)型是基于面心立方體（FCC）的結(jié)構(gòu)，Ag46的構(gòu)型是無序形態(tài)的，而Ag49的構(gòu)型是基于二十面體的結(jié)構(gòu)。

圖1 4個(gè)代表性算例的最優(yōu)構(gòu)型

5 結(jié)語

本文針對(duì)Gupta 勢(shì)能描述的銀團(tuán)簇的結(jié)構(gòu)優(yōu)化問題，提出了一個(gè)PBH 算法的新變種，并對(duì)原子數(shù)在10 ≤N≤61范圍內(nèi)的銀團(tuán)簇進(jìn)行了優(yōu)化。與文獻(xiàn)中已有的PBH 算法不同的是，本文提出的PBH 算法的變種利用了連通表來定義團(tuán)簇構(gòu)型間的距離函數(shù)。計(jì)算結(jié)果表明，該算法對(duì)于Gupta 勢(shì)能描述的銀團(tuán)簇的結(jié)構(gòu)優(yōu)化是有效的。具體而言，在10 ≤N≤61范圍內(nèi)，提出的算法能夠獲得除4個(gè)難例以外的已知全局最優(yōu)構(gòu)型。實(shí)驗(yàn)結(jié)果再次表明，連通表是定義團(tuán)簇構(gòu)型間距離函數(shù)的一種很好的工具。在將來，我們打算在其他原子團(tuán)簇上驗(yàn)證它的有效性。