求解加權(quán)MTSP問題的CUDA并行群智能方法

蘇守寶， 趙 威， 李 智

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

0 引言

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

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

1 相關(guān)研究

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

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

Cij=dij·vij。

(1)

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

(2)

(3)

1.2 粒子群K-means聚類算法

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

(4)

聚類中心更新公式為

(5)

粒子群優(yōu)化(particle swarm optimization，PSO)算法中每個(gè)粒子根據(jù)自身經(jīng)驗(yàn)數(shù)據(jù)和全局經(jīng)驗(yàn)數(shù)據(jù)利用迭代來逼近全局最優(yōu)解。每個(gè)粒子用如下計(jì)算式更新自己的聚類中心坐標(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個(gè)粒子所經(jīng)歷過的最優(yōu)位置；Gbest為目前全局最優(yōu)位置；粒子的運(yùn)動(dòng)用位置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問題時(shí)，初始化各個(gè)城市間的信息素濃度，在時(shí)刻t螞蟻k從當(dāng)前城市i轉(zhuǎn)移到下一個(gè)城市j的概率為

(8)

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

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

(9)

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

2 本文方法

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

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

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

(10)

(11)

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

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

(12)

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

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

(13)

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

GPSO-AC算法如下：

初始化粒子群

do

parallel_for 粒子 in 粒子群;

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

parallel { // 重分類

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

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

sort 樣本 by 距離;

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

end

}

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

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

初始化蟻群

do

parallel_for 螞蟻 in 蟻群;

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

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

end

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

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

end

}

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

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

end

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

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

2.2.1 算法實(shí)現(xiàn)

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

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

GPU的1個(gè)SM上有32個(gè)CUDA核心，每32個(gè)線程組成1個(gè)線程束(warp)，CUDA核心只能以線程束作為基本單元執(zhí)行指令，所以1個(gè)block內(nèi)的線程束必須為32的倍數(shù)。在設(shè)計(jì)kernel函數(shù)時(shí)，線程塊數(shù)通過計(jì)算得到，例如假定每個(gè)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)致指令異化，那么這時(shí)一個(gè)線程束中的線程串行運(yùn)行在CUDA核心上，一個(gè)線程采用多個(gè)kernel函數(shù)分解。

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

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

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

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

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

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

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

Step5計(jì)算已遍歷的路徑長度，與全局最優(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)時(shí)需要做內(nèi)存優(yōu)化，可用2種編碼方式定義n個(gè)粒子，編碼方式如下。

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

4 結(jié)論

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