蘇守寶， 趙 威， 李 智

(1.江蘇科技大學(xué) 計算機(jī)學(xué)院,江蘇 鎮(zhèn)江 212003；2.金陵科技學(xué)院 數(shù)據(jù)科學(xué)與智慧軟件江蘇省重點(diǎn)實驗室,江蘇 南京 211169)

0 引言

多旅行商(MTSP)問題是旅行商(TSP)問題的拓展，屬于NP-hard組合優(yōu)化問題，即m個旅行商前往n個城市銷售商品，每個城市有且只能有一個旅行商經(jīng)過，要求所有旅行商經(jīng)過的總距離之和最短。壽濤等[1]通過Delaunay三角剖分及樹分解算法，將MTSP問題轉(zhuǎn)化為多個TSP問題。劉楠[2]采用哈密頓圈分割覆蓋方法，設(shè)計了求解MTSP問題的近似精確解法。這些算法性能穩(wěn)定、時間復(fù)雜度小，但只能解決小規(guī)模MTSP問題。于是，啟發(fā)式算法得到了廣泛應(yīng)用，主要包括混合迭代法、分治法、協(xié)同進(jìn)化策略等。Trigui等[3]將多機(jī)器人任務(wù)分配(MRTA)看作MTSP的實例，同時優(yōu)化最大行駛距離和總行駛距離，問題成本是這2個指標(biāo)的組合，以簡化為單一目標(biāo)優(yōu)化問題；張美燕等[4]將自主式水下潛器(AUV)的能量消耗和能量均衡作為MTSP問題路徑上的代價，運(yùn)用MTSP-GA算法模型在水下三維空間內(nèi)進(jìn)行AUV路徑規(guī)劃。在確保工作負(fù)載均衡的情況下，有研究者提出了MTSP問題的兩階段啟發(fā)式算法TPHA，通過改進(jìn)的K-means分組確保城市數(shù)量均衡，再用改進(jìn)的遺傳算法(GA)求解各聚類城市點(diǎn)TSP問題[5-7]。也有學(xué)者使用分治策略(如單親GA)將MTSP問題分組，再用ACO算法逐個求解TSP問題[8-10]，如在GA迭代過程中加入ACO算法作為算法收斂約束，則將GA改為單親GA以提高算法速度[11-12]。分治策略算法之間無交互，因而簡化了問題的復(fù)雜度，但算法魯棒性不強(qiáng)，尤其是大規(guī)模數(shù)據(jù)仍面臨搜索空間指數(shù)增長和維數(shù)災(zāi)難問題[13-14]。

針對混合迭代法算法精度高但執(zhí)行時間長的問題，本文基于并行軟硬件體系CUDA運(yùn)算平臺，提出一種基于GPU的混合粒子群聚類蟻群的并行算法，在保證算法精度的同時，提高了混合迭代法的執(zhí)行速度。由于實際場景中旅行商在每個路段上開銷不同，將其抽象為每段路程區(qū)間上都有一個與之對應(yīng)的代價，將路程代價考慮到MTSP問題中。

1 相關(guān)研究

1.1 加權(quán)MTSP問題建模

單點(diǎn)出發(fā)的加權(quán)MTSP問題：m個旅行商(記為b1,b2，…，bm)前往n個城市(記為T1,T2，…，Tn)售貨，所有旅行商從同一個城市T1出發(fā)并最終在城市T1會合，要求每個城市有且只能有1個旅行商經(jīng)過。考慮到每條路段的開銷不同，記連接i、j兩地的路段上的單位代價為vij，則i、j兩地旅行代價Cij為

Cij=dij·vij。

(1)

求所有旅行商總代價和最短路線方案，目標(biāo)函數(shù)為

(2)

(3)

1.2 粒子群K-means聚類算法

粒子群K-means聚類是將聚類中心設(shè)為粒子的位置X，維度是聚類個數(shù)，則粒子群K-means聚類粒子最佳位置即為聚類最優(yōu)解。將N個樣本聚類到K個類中，需要滿足：

(4)

聚類中心更新公式為

(5)

粒子群優(yōu)化(particle swarm optimization，PSO)算法中每個粒子根據(jù)自身經(jīng)驗數(shù)據(jù)和全局經(jīng)驗數(shù)據(jù)利用迭代來逼近全局最優(yōu)解。每個粒子用如下計算式更新自己的聚類中心坐標(biāo)和速度[6-7]：

Vi(t+1)=wVi(t)+C1·rand()·(Pbesti(t)-

Xi(t))+C2·rand()·(Gbest(t)-Xi(t))；

(6)

Xi(t+1)=Xi(t)+Vi(t+1)。

(7)

式中：Pbesti為第i個粒子所經(jīng)歷過的最優(yōu)位置；Gbest為目前全局最優(yōu)位置；粒子的運(yùn)動用位置X和速度V表示；w為慣性因子，其值為(0,1)，影響著粒子的全局搜索能力，w越大粒子全局搜索能力越強(qiáng)，反之越弱；rand()為(0,1)的隨機(jī)值；C1和C2通常取值在2附近。粒子群K-means聚類算法流程如圖1所示。

圖1 粒子群K-means聚類算法流程圖Figure 1 Flow chart of PSO-based K-means method

1.3 蟻群優(yōu)化算法

蟻群優(yōu)化(ant colony optimization，ACO)算法是一種模擬螞蟻覓食行為的群智能算法，具有可并行執(zhí)行、易于與其他算法混合、較強(qiáng)魯棒性等特點(diǎn)[8-9]。蟻群優(yōu)化算法求解TSP問題時，初始化各個城市間的信息素濃度，在時刻t螞蟻k從當(dāng)前城市i轉(zhuǎn)移到下一個城市j的概率為

(8)

式中：α為信息素因子，表示各個路徑被選擇的比重大小；β為啟發(fā)式因子，表示啟發(fā)函數(shù)ηij的重要程度。螞蟻k周游一周后形成路徑，計算路徑長度更新各邊信息素。當(dāng)所有螞蟻都走完后，由最短路徑對信息素和禁忌表進(jìn)行更新：

τij(t+n)=(1-ρ)τij(t)+Δτij，

(9)

式中：ρ∈(0,1)為信息素?fù)]發(fā)系數(shù)；Q為信息素總量；Lk為螞蟻k周游一次的路徑長度。

2 本文方法

為了利用CUDA實現(xiàn)并行計算，先用粒子群K-means算法將各個城市點(diǎn)按旅行商個數(shù)進(jìn)行分類，再使用蟻群算法進(jìn)行分類評估，n個粒子搜索最優(yōu)解過程中有n個蟻群參與迭代，二者混合迭代過程中會出現(xiàn)大量獨(dú)立運(yùn)算過程，即粒子、蟻群間和蟻群內(nèi)各個螞蟻間的獨(dú)立運(yùn)算。這些獨(dú)立運(yùn)算過程中使用的數(shù)據(jù)重新排列為向量化數(shù)據(jù)，GPU的流處理器單個時鐘周期內(nèi)同時調(diào)度多個數(shù)據(jù)，如果向量的元素彼此獨(dú)立則可以并行計算向量，這也與GPU利用SIMT指令架構(gòu)加速原理是一致的[15]，因此，提出一種基于CUDA的并行群智能方法，記為GPSO-AC。

2.1 GPSO-AC算法設(shè)計

GPSO-AC算法前半部分采用基于粒子群的K-means算法，聚類中心為粒子位置，維度為聚類個數(shù)，粒子速度為聚類中心的偏移量。由于K-means算法初始聚類中心對結(jié)果影響很大，選擇合適的初始聚類中心可以加快算法收斂速度和提升聚類質(zhì)量。本文采用粗分類方法對初始城市分類，則從各個城市前往其他城市的所有代價和全局城市路段互通的代價總和分別為

(10)

(11)

式中：n為城市數(shù)量；vij為城市i和j之間的單位代價。

將城市到聚類中心距離排序，再依次遍歷，若遍歷條件滿足式(12)，則將已經(jīng)遍歷的城市歸為一類，按照式(5)更新聚類中心。

(12)

式中：m為旅行商個數(shù)即聚類個數(shù)；f為布爾標(biāo)記，如果旅行商k經(jīng)過城市i則為1，反之為0。

在粒子群迭代流程中，所有重新分類的過程均使用式(12)進(jìn)行約束，目的是確保每類城市之間遍歷代價總和大致均衡。接著細(xì)分類，計算適應(yīng)度，使用ACO算法計算遍歷某一類城市的最小代價，使用式(9)更新信息素和禁忌表，得出的最小代價為Ck。當(dāng)前粒子的適應(yīng)度為

(13)

式中：t為粒子個數(shù)；S(·)為計算方差的函數(shù)；S(Ck)表示聚類結(jié)果之間的離散程度，S(Ck)越小則各類之間遍歷代價越均衡；q∈[0,1]為調(diào)諧系數(shù)，當(dāng)q=0時，表示不考慮代價均衡，只計算最低代價，可根據(jù)實際問題調(diào)整取值。

GPSO-AC算法如下：

初始化粒子群

do

parallel_for 粒子 in 粒子群;

parallel{更新位置和速度,使用式(5)～(7)};

parallel { // 重分類

for 聚類中心 in 所有聚類結(jié)果;

計算樣本到當(dāng)前聚類中心距離;

sort 樣本 by 距離;

聚類中心選擇樣本,使用式(4);

end

}

parallel { // 計算適應(yīng)度

for 聚類中心 in 所有聚類結(jié)果;

初始化蟻群

do

parallel_for 螞蟻 in 蟻群;

parallel{搜索路徑,使用式(8)};

parallel{更新信息,使用式(9)};

end

while 達(dá)到最大迭代次數(shù)

計算最優(yōu)路徑適應(yīng)度,使用式(10)～(13);

end

}

parallel{更新個體最優(yōu)};

parallel{更新群體最優(yōu)};

end

while 達(dá)到最大迭代次數(shù)。

2.2 GPU并行算法實現(xiàn)與優(yōu)化

2.2.1 算法實現(xiàn)

在CUDA架構(gòu)中，將算法拆解成多個并行算子，單個并行算子稱為線程(thread)，多個線程組成1個線程塊(block)。根據(jù)粒子群算法固有的并行特點(diǎn)，可將圖1中的算法流程進(jìn)行拓展，如圖2所示。

圖2 基于GPU的粒子群聚類算法流程圖Figure 2 Flowchart of PSO clustering based on GPU

GPU的1個SM上有32個CUDA核心，每32個線程組成1個線程束(warp)，CUDA核心只能以線程束作為基本單元執(zhí)行指令，所以1個block內(nèi)的線程束必須為32的倍數(shù)。在設(shè)計kernel函數(shù)時，線程塊數(shù)通過計算得到，例如假定每個block內(nèi)線程數(shù)為32，設(shè)粒子數(shù)即并行線程數(shù)op-size為O，則線程塊數(shù)B=(O+32-1)/32。考慮到GPU使用SIMT指令架構(gòu)，線程內(nèi)指令相似度越高其執(zhí)行速度越快，如果線程內(nèi)算法流程過于復(fù)雜導(dǎo)致指令異化，那么這時一個線程束中的線程串行運(yùn)行在CUDA核心上，一個線程采用多個kernel函數(shù)分解。

在更新Gbest前，粒子間互不干擾，線程獨(dú)立運(yùn)行；更新Gbest時，線程間有數(shù)據(jù)交換。CUDA上無鎖機(jī)制，但是有同步原語syncthreads()函數(shù)，同步實際上是將線程串行運(yùn)行，時間復(fù)雜度為O(N)。使用并行規(guī)約算法求解Gbest，規(guī)約算法對傳入的N個數(shù)據(jù)，使用一個二元的符合結(jié)合律的操作符?，生成1個結(jié)果。求Gbest的規(guī)約可表示為Gbest=Pbest(1)?Pbest(2)?Pbest(3)?…?Pbest(N)。首先把粒子所有的Pbest數(shù)據(jù)放至共享內(nèi)存中并編號，之后用1個線程計算前2個粒子的Pbest，再用1個線程計算中間2次結(jié)果，以此迭代最終求得Gbest，時間復(fù)雜度為O(lgN)。圖3為并行規(guī)約算法的示意圖。

圖3 并行規(guī)約算法Figure 3 Sketch of parallel protocol

ACO算法中螞蟻在搜索時互不依賴，各螞蟻個體搜索邏輯一致僅數(shù)據(jù)不同，可將螞蟻個體獨(dú)立看待，每個個體對應(yīng)CUDA中的一個線程，并行中的線程指令一致但處理的數(shù)據(jù)不一致。多個粒子同時進(jìn)行適應(yīng)度計算時，實際上是并行了多個蟻群，每個蟻群內(nèi)部再進(jìn)行二次并行。ACO算法并行過程如下。

Step1計算需要的內(nèi)存并分配給device端；

Step2數(shù)據(jù)預(yù)處理，初始化隨機(jī)種子、計算距離矩陣、初始化禁忌表等；

Step3單個線程初始化，重復(fù)執(zhí)行n-1次，計算轉(zhuǎn)移概率寫入概率表，選擇下一城市；

Step4線程同步，等待所有線程完成路徑選擇；

Step5計算已遍歷的路徑長度，與全局最優(yōu)路徑比較，如果優(yōu)于全局最優(yōu)則對其更新，根據(jù)當(dāng)前路徑更新禁忌表；

Step6更新信息素矩陣，清空禁忌表；

Step7如果滿足程序結(jié)束條件則停止，否則跳轉(zhuǎn)Step 3。

2.2.2 針對GPU的編程優(yōu)化

CUDA程序75%的性能瓶頸在內(nèi)存交互上[15]，定義數(shù)據(jù)結(jié)構(gòu)時需要做內(nèi)存優(yōu)化，可用2種編碼方式定義n個粒子，編碼方式如下。

方式1：

struct Particle_t{

doublex;

doubley;

double fitness;

};

Particle_t*particle

=new Particle_t[n]。

方式2：

struct Particle_t{

double*x;

double*y;

double*fitness;

} particle;

particle.x=new double[n];

particle.y=new double[n];

particle.fitness=new double[n]。

首先要減少CUDA線程讀取數(shù)據(jù)時產(chǎn)生的cache miss, cache miss是指warp從L1 Cache尋址失敗繼而請求從全局內(nèi)存(global memory)中讀取數(shù)據(jù)的過程。warp是CUDA中最小指令執(zhí)行單元，如圖4(a)所示，線程從global memory中讀取數(shù)據(jù)不是單個線程依次讀取，而是讀取整個warp所需要的數(shù)據(jù)，經(jīng)由L2 Cache到達(dá)L1 Cache，然后按照數(shù)據(jù)地址線性訪問L1 Cache中的數(shù)據(jù)塊。

圖4 CUDA線程讀取內(nèi)存過程Figure 4 CUDA thread read memory process

使用第1種編碼方式，內(nèi)存訪問情況如圖4(b)所示。假設(shè)此刻n個粒子線程同時操作x變量，線程2對應(yīng)的x變量地址位于線程1的x地址+size of(Particle_t)處；若粒子線程n從數(shù)據(jù)塊中找不到對應(yīng)的x變量，則會從全局內(nèi)存重新讀取一塊數(shù)據(jù)，這就產(chǎn)生了一次cache miss。最壞的情況下，n個線程能產(chǎn)生n次cache miss，這會對程序性能產(chǎn)生嚴(yán)重影響。使用第2種編碼方式則不會產(chǎn)生頻繁的內(nèi)存加載、內(nèi)存訪問，如圖4(c)所示，能極大提高程序性能，同理，蟻群的數(shù)據(jù)結(jié)構(gòu)定義也如此。另外，優(yōu)先使用CUDA中的共享內(nèi)存(shared memory)，一個block內(nèi)所有線程是共享內(nèi)存的，相比于全局內(nèi)存速度更快，ACO算法的禁忌表和轉(zhuǎn)移概率表需要放在全局內(nèi)存中。

3 實驗結(jié)果與分析

實驗所用CPU為Intel? CoreTMi5-8400，RAM8 GB，GPU為NVIDIA GeForce GTX 1060， Ubuntu 18.04.5 LTS操作系統(tǒng)， GCC7.5 CUDA10.1編譯器。

使用TSPLIB中eil51數(shù)據(jù)集，旅行商數(shù)為3，不考慮代價均衡的情況，各個城市點(diǎn)間的權(quán)值設(shè)為1，要解決單點(diǎn)出發(fā)的MTSP問題，分別用GPSO-AC算法(使用GPU加速)和PSO-AC算法(未使用GPU加速)作對比。由表1可知，PSO-AC算法運(yùn)行時間遠(yuǎn)遠(yuǎn)大于GPSO-AC。在實時性要求較高的系統(tǒng)中使用，PSO-AC算法無實用性，并且隨著群規(guī)模增大，其運(yùn)行時間呈線性增長。而GPSO-AC隨群規(guī)模增大，算法運(yùn)行時間增幅較小，加速比增大，這是因為算法聚類過程趨于穩(wěn)定后，質(zhì)心和聚類結(jié)果在算法中后期不再變化，算法程序在GPU中指令異化情況大大減少，進(jìn)一步提高了并行的執(zhí)行效率。

表1 算法運(yùn)行時間對比Table 1 Algorithm running time comparison

使用TSPLIB中6個數(shù)據(jù)集，將GPSO-AC和PSO-AC[7]、TPHA[5]及K-means-AC[11]算法進(jìn)行對比，旅行商數(shù)為3，群規(guī)模為64，最大迭代次數(shù)為500，根據(jù)經(jīng)驗，學(xué)習(xí)因子C1和C2均為1.97，α=1，β=3，實驗數(shù)據(jù)如表2所示，GPSO-AC算法部分實驗最優(yōu)路線如圖5、6所示。

表2 4種算法實驗結(jié)果對比Table 2 Experimental results comparisons of four algorithms

圖5 GPSO-AC在eil101上的運(yùn)行結(jié)果Figure 5 Result of GPSO-AC on eil101

圖6 GPSO-AC在ch150上的運(yùn)行結(jié)果Figure 6 Result of GPSO-AC on ch150

分析表2數(shù)據(jù)可知，未使用GPU加速的PSO-AC算法收斂結(jié)果優(yōu)于TPHA和K-means-AC這2種兩步式算法，但算法運(yùn)行時間較長。兩步式算法遵循先分類再計算的思路，運(yùn)行速度較快，但收斂精度較差。使用GPU加速的GPSO-AC算法雖然運(yùn)行時間比兩步式算法長，但仍處于合理可接受的范圍內(nèi)且收斂精度優(yōu)于兩步式算法。2種算法均采用聚類對城市點(diǎn)進(jìn)行有效分組，TPHA只考慮分組而算法后期未考慮對分組進(jìn)行優(yōu)化，導(dǎo)致其收斂結(jié)果不佳，這也是兩步式算法共同的弊端。

使用chn31數(shù)據(jù)集，在不考慮代價均衡、代價均衡約束、加權(quán)代價均衡的情況下，當(dāng)旅行商數(shù)分別取1、2、3、4時GPSO-AC運(yùn)行結(jié)果(遍歷代價)如表3所示。在旅行商為4時，3種情況下的運(yùn)行結(jié)果如圖7所示。由表3可知，算法考慮代價均衡約束后總代價會變大，雖然各個旅行商的最優(yōu)解都接近，但標(biāo)準(zhǔn)差遠(yuǎn)小于不考慮代價均衡的情況?？紤]代價均衡的協(xié)調(diào)系數(shù)q=0.05，如圖7(b)所示。加權(quán)代價均衡時，v06和v46設(shè)為30，即C0～C6路段、C4～C6路段上的單位代價為30，協(xié)調(diào)系數(shù)q=0.1，此時各個旅行商最優(yōu)解偏差更小，同時繞開了C0～C6和C4～C6路段，選擇了相鄰代價小的路線,如圖7(c)所示。各旅行商間的周游代價離散程度越小，總代價則越大，總代價最低和旅行商間的代價均衡不可能同時滿足，需要根據(jù)實際問題進(jìn)行取舍。

表3 GPSO-AC對各旅行商的運(yùn)行代價結(jié)果Table 3 Results of GPSO-AC for traveling salesmen

圖7 旅行商為4時GPSO-AC在chn31上的運(yùn)行結(jié)果Figure 7 Results of GPSO-AC on chn31 with 4 salesmen

4 結(jié)論

針對混合迭代法算法運(yùn)行時間長的問題，本文根據(jù)粒子群和蟻群算法良好的并行性，提出一種基于CUDA的混合算法GPSO-AC。該算法充分利用GPU多流處理器的設(shè)計和單指令多線程的指令架構(gòu)特點(diǎn)，將算法中大量獨(dú)立運(yùn)算過程同時執(zhí)行。GPSO-AC算法在求解一般MTSP問題及其衍生加權(quán)、代價均衡MSTP問題上，不僅加快了混合迭代法的執(zhí)行速度同時確保了收斂精度。最后，討論了加權(quán)MTSP問題中代價均衡和總代價最優(yōu)之間的關(guān)系。在實際場景中旅行商在各個路段上開銷不同，可抽象為代價權(quán)重加權(quán)MTSP問題。對于各個旅行商如何均衡開銷即代價均衡的加權(quán)MTSP問題，在保證所有旅行商總路程最短的前提下，代價均衡和總代價最優(yōu)難以同時滿足，仍需進(jìn)一步研究。