求解多起點多旅行商問題的K-means聚類信息傳播算法

程亞南， 王曉峰,2*， 劉凇佐， 莫淳惠

( 1. 北方民族大學計算機科學與工程學院, 銀川 750021； 2. 北方民族大學,圖像圖形智能處理國家民委重點實驗室, 銀川 750021 )

旅行商問題(traveling salesman problem,TSP)作為NP難問題，近年來研究的熱度只增不減，多旅行商問題(multiple traveling salesman problem, MTSP)是其擴展問題。相對于傳統(tǒng)的TSP問題來說，MTSP更符合現(xiàn)實社會的復(fù)雜條件，在現(xiàn)實生活中應(yīng)用更加廣泛，在無人機路徑規(guī)劃、不同社區(qū)快遞配送、調(diào)度等實際問題中比較貼合。相較于TSP問題，同樣作為NP難問題的MTSP的求解難度在于求解的過程中需要進行分組優(yōu)化，組內(nèi)進行TSP的遍歷。由于MTSP貼合實際情況的特性，近年來在學術(shù)界關(guān)注度較高，所以有效地求解MTSP的方法有助于解決現(xiàn)實生活中的問題。但是目前對于求解MTSP的算法研究還不夠深入。在以往的研究中大致可以分為兩類：精確算法和啟發(fā)式算法。分支限界是現(xiàn)階段為止最為精確的算法，可以用來求解TSP和小規(guī)模的對稱型MTSP問題，精確算法的發(fā)展雖然解決了一些問題，但是其本身嚴格的數(shù)學定義，在求解時受到很大的限制，隨著城市規(guī)模的增大，MTSP問題的維度會迅速增加，精確算法求解的難度逐漸提升，難以在合理的時間找到最優(yōu)解?；诖耍瑔l(fā)式算法越來越受到研究者的青睞。文獻[1]提出一種改進分組遺傳算法，設(shè)計了新的分組編碼構(gòu)造了一種快速交叉算子，同時，設(shè)計了一種新的局部交叉算子用來提高算法的求解精度。文獻[2]提出一種新的求解MTSP的斷點算子-初等斷點算子，主要是通過初等矩陣運算來生成MTSP的解空間。文獻[3]利用模糊C均值聚類按照城市隸屬度劃分類別，用單親遺傳算法求解分類的旅行商問題。近年來，MTSP問題引起眾多研究人員的關(guān)注，一些新的啟發(fā)式算法被改進解決MTSP問題。文獻[4]利用改進的人工蜂群算法(artificial bee colony，ABC)與入侵雜草算法(invasive weed algorithm, IWA)相結(jié)合求解MTSP問題。文獻[4]針對MMTSP問題，提出了一種改進的兩部分狼群搜索算法(modified two-part wolf pack search, MTWPS)，算法的全局性得到了提高。

上述啟發(fā)式算法在局部搜索與全局搜索的平衡和求解質(zhì)量上存在不足?，F(xiàn)主要針對多旅行商問題中的多起點進行研究，提出了一種基于K-means聚類方法的信息傳播算法來求解。信息傳播算法的思想起源于統(tǒng)計物理學，信息傳播算法可以解決多種圖模型上的概率計算問題，在消息傳播的過程是并行實現(xiàn)的，時間復(fù)雜度得到了降低，所以選擇信息傳播算法進行分組后的旅行商問題進行求解。

1 基本知識

1.1 多旅行商問題

MTSP[5]是傳統(tǒng)TSP問題的一種擴展，與傳統(tǒng)的TSP不同之處在于，MTSP是m個推銷員訪問n個城市(n>m)，保證每個城市只訪問一次。目標依舊是訪問距離最短也就是m個旅行商旅行的距離之和最小。多旅行商問題會根據(jù)不同的起點城市和每個旅行商訪問的城市個數(shù)分為不同的類型?，F(xiàn)主要針對多起點閉合路線的MTSP問題。m個旅行商從不同的城市出發(fā)，訪問一定數(shù)量的城市后返回出發(fā)城市，要求每個城市只能被一個旅行商訪問一次。

(1)

1.2 K-means聚類方法

K-means[6]屬于無監(jiān)督學習的聚類算法，算法本身比較簡單，對于給定的數(shù)據(jù)集，算法按照樣本點之間的距離大小進行劃分，劃分結(jié)果是K個簇也可稱為k類，同簇之間相似度較高，簇與簇之間相似度較低。算法迭代終止的條件是類簇中心點變化較少或者達到預(yù)先設(shè)定的迭代次數(shù)?；驹硎墙o定的數(shù)據(jù)集樣本D={x1,x2,…,xm}，初始中心點，從D中隨機選擇k個樣本作為初始均值向量，計算樣本與均值之間的距離，根據(jù)距離最近均值向量確定的簇標記，重新計算更新均值向量，均值向量未更新則輸出簇劃分C={C1,C2,…,CK}。

算法的基本步驟如下。

步驟1選定聚類的個數(shù)k，隨機初始化k個中心點。

步驟2針對數(shù)據(jù)集中的每個樣本點，找到距離其最近的中心點，按照樣本點距離每個中心點的距離選擇聚類。

步驟3聚類前后樣本點不再發(fā)生變化，聚類達到穩(wěn)態(tài)，則算法終止，否則進入步驟4。

步驟4對步驟3聚類后的類別進行中心點的更新，轉(zhuǎn)至步驟3。

將聚類算法應(yīng)用到旅行商問題上，假設(shè)城市個數(shù)n=9,旅行商個數(shù)m=3，9個城市進行編號聚類設(shè)置k=3，這里每一個旅行商訪問的城市數(shù)目并不均等，就可得到{{1,3,8,1}{2,4,7,2}{5,6,9,5}} 3組數(shù)目不等的城市組合，在第二階段對每一組按照傳統(tǒng)的旅行商問題進行遍歷排序就可以。

1.3 因子圖

因子圖[7]是一個特殊的二分圖，代表由許多變量組成的全局函數(shù)分解成局部函數(shù)的圖示。信息傳播算法以因子圖為工具，可以直觀地表現(xiàn)一個多種約束條件的復(fù)雜問題。如函數(shù)p(x)可以分解成fA、fB、fC、fD、fE5個局部函數(shù)的乘積，構(gòu)造對應(yīng)關(guān)系為式(2)的因子圖如圖1所示。

p(x1,x2,x3,x4,x5)=fA(x1,x2,x3)fB(x2)·

fC(x3,x4,x5)fD(x4)fE(x5)

(2)

例如，對于一個CNF公式：

(3)

其因子圖如圖2所示。

圖1 因子圖實例Fig.1 Example of factor graph

α1=(x1∨x2∨x3)，α2=(x1∨x2∨x4)， α3=(x2∨x3∨x5) 圖2 因子圖實例Fig.2 Example of factor graph

1.4 最小和信息傳播算法

1988年人工智能領(lǐng)域著名學者Pearl提出了置信傳播(belief propagation, BP)算法，BP算法把全局的概率推理過程轉(zhuǎn)變?yōu)榫植孔兞块g的消息傳遞，從而降低了推理的復(fù)雜度。因為BP算法的優(yōu)點，自提出起，就受到了國際上眾多學者的關(guān)注，掀起了研究的熱潮。

置信傳播算法[8]的特征是基于因子圖上的邊傳遞消息，當這種消息傳遞達到一種穩(wěn)態(tài)時，可計算因子圖上節(jié)點取值的邊際概率，從而以高概率地確定變元的某種取值。研究結(jié)果表明，信息傳播算法求解組合優(yōu)化問題性能較好。當前，置信傳播算法除了以“和積”的形式傳遞信息，近年來出現(xiàn)由和積算法簡化的最小和算法，目前國內(nèi)外主要是針對和積算法進行研究。

最小和信息傳播算法[9]是一個信息逐步迭代的過程，在因子圖的邊 (i,α)上信息隨著過程的迭代進行更新。μi→α(di)表示由變量節(jié)點發(fā)出給因子節(jié)點的信息；代表變量i取值為di的信息；μα→i(di)表示由因子節(jié)點發(fā)給變量節(jié)點的信息。信息迭代方程為

(4)

(5)

式中：V(i)為與i相連接的因子節(jié)點集合；V(i)a為不包含a的因子節(jié)點集合；V(a) 為與a相連接的變量節(jié)點集合；V(a)i為不包含i的變量節(jié)點集合；fa(dj)是對應(yīng)問題的描述函數(shù)，也稱為勢函數(shù)。當BP收斂時，邊際信念表示為

(6)

2 求解MMTSP的K-means聚類信息傳播算法

2.1 旅行商問題的線性規(guī)劃方程

旅行商問題可以用非負賦權(quán)無向圖G(V,E)表示，旅行商問題的約束條件為每一個城市只能一條邊進，一條邊出，即旅行商問題的線性規(guī)劃方程為

minimizewTx

xe∈[0,1]|E|

(7)

式(7)中：圖中各邊的權(quán)重集合為w=(w1,w2,…,wm)，邊的標簽屬性集合為x=(x1,x2,…,xm)T，標簽x1=1，該邊加入回路，反之，不加入回路；δ(v)是頂點v所連接的邊；xe是頂點v加入回路邊數(shù)量；|E|是邊數(shù)目的模。

2.2 旅行商問題的因子圖

圖結(jié)構(gòu)用G(V,E)表示，i和j兩個節(jié)點之間的距離表示為wi:j，V={v1,v2,…,vn}為節(jié)點集合，節(jié)點之間存在邊相連則(i:j)=1。令d={d1,d2,…,dM}∈{0,1}M是一組M維二進制變量，其中M=|E|。如果節(jié)點i和節(jié)點j所在的邊在TSP問題的可行解中則di:j賦值為1。

根據(jù)TSP問題的特性定義兩種條件約束。

(1)成本約束：代表因子圖上邊的成本，當i和j兩個節(jié)點存在邊，信息為非負距離，反之為0。

(8)

(2)勢函數(shù)：確保每個結(jié)點剛好都與其他兩個結(jié)點連接。

(9)

這里用小規(guī)模的旅行商問題給出無向圖轉(zhuǎn)換因子圖的示例，如圖3所示，圖3(a)為四個城市的簡單無向圖，城市分別用a、b、c、d表示，城市之間的權(quán)重用無向邊上的數(shù)字表示。邊上的數(shù)字代表城市之間的權(quán)重，圖3(b)中變量節(jié)點為城市之間的邊構(gòu)成，因子節(jié)點為城市和約束條件(勢函數(shù))，其中白框和黑框分別表示城市和約束條件。

圖3 四個城市的無向圖示例Fig.3 Example factor diagrams for four cities

2.3 旅行商問題的算法優(yōu)化

文獻[10]針對消息傳遞的復(fù)雜性,提出簡化信息的方法。根據(jù)TSP的特性， 當di:1=1時，表示節(jié)點i和節(jié)點j的邊選入路徑，反之則不選入。令min[k]A表示集合A中第k個最小值，將式(5)進行重寫得到消息的更新公式為

(10)

在BP算法迭代過程中，消息基于因子圖進行傳遞，會產(chǎn)生消息震蕩，進而影響算法的收斂速度，動態(tài)阻尼[11]是將最近兩次的迭代信息加和求取平均值，經(jīng)處理后提高BP算法信息更新的收斂速度。動態(tài)阻尼的公式為

(11)

2.4 模擬退火算子

模擬退火的出發(fā)點是在固體退火的過程中和組合優(yōu)化算法的求解是類似的，作為一種隨機迭代思想是在搜索過程中進行隨機搜索，具有突變性質(zhì)，解決容易陷入局部最優(yōu)的缺陷，模擬退火算子具有跳出局部最優(yōu)的能力，全局搜索能力較強，能夠有效求解旅行商問題。

其中模擬退火算子是以一定的概率接受新的解，當在溫度T時，從當前解curr改變?yōu)樾陆鈔ew；若curr

2.5 局部搜索操作

局部搜索(local search, LS)的過程中采用一定概率選擇交換、逆序兩種操作，在信息傳播算法全局迭代之后采用交換、逆序進行局部搜索，便于找到最好的解。

2.5.1 交換操作及示例

交換操作是在原始解中任意選擇兩個城市進行交換，如圖4所示，左側(cè)是原始路徑經(jīng)過選擇城市3和城市7，變換后城市序列如右側(cè)所示。

2.5.2 逆序操作及示例

逆序操作和交換操作一樣同樣任意選取兩個城市，將以兩個城市為首尾的城市逐一逆序輸出，如圖5所示，左側(cè)為原始路徑，右側(cè)為經(jīng)過逆序操作的路徑圖。

圖4 交換操作Fig.4 Exchange operation

圖5 逆序操作Fig.5 Reverse order operation

時間復(fù)雜度為O(n3)，算法具體步驟如下。

算法1:求解旅行商問題的信息傳播算法輸入: 圖G=(V,E),加權(quán)鄰接矩陣A,最大迭代Tmax,閾值εmax輸出:巡視中路徑Tour?E1 構(gòu)造初始化因子圖2 初始化勢函數(shù)信息^μi:j→?vi←0?i∈V,j∈?vi 3 初始化^φi:j←wi:j?(i:j)∈E 4 While True:4.1 ε←0,T←04.2 Whileε<εmax and T

3 數(shù)值實驗及分析

3.1 K-means聚類分組分析

K-means在聚類算法思想簡單，聚類時間快，但是也存在一些缺點，初始點的隨機產(chǎn)生對算法會造成一定影響，其中改進的K-means算法較多，基于K-means++[12]的基礎(chǔ)上使用長度為m的獨立馬爾可夫鏈在每一次迭代中進行中心采樣。隨機選擇一個中心點，在迭代的過程中通過概率計算進行中心點的更新。

聚類方法的偽代碼如下。

算法2:K-means質(zhì)心選擇輸入:數(shù)據(jù)集D,聚類個數(shù)K,鏈長M輸出:K個質(zhì)心的矩陣Cc1←random.choice(D),C←c1For x in D:q(x)←d(x,c1)2/2[∑x'∈Dd(x',c1)2]+1/2DFor i in range(2,K+1):x←random.choice(D,p=q(x)) #根據(jù)概率隨機選擇dx←d(x,ci-1)2For j in range(2,M+1): y←random.choice[D,p=q(y)] dy←d(y,ci-1)2 If d(y)×q(y)>Unif(0,1)×d(x)×q(x): x←y,dx←dy ci←ci-1∪{x}Return C

比較了K-means和改進的K-means++的聚類效果，發(fā)現(xiàn)改進后避免了隨機化中心點對算法造成的影響。在聚類效果上進行對比，對比效果如圖6、圖7所示。

在對數(shù)據(jù)集的處理上，基本的聚類方法，聚類點不均，存在離散點分簇不明，經(jīng)過多次實驗對比，測試改進算法針對旅行商問題分組明朗，不存在異常點，求解距離更短，所以選擇改進的K-means++進行MMTSP的分組。

圖6 Kroa200分組效果圖Fig.6 Kroa200 group renderings

3.2 算法性能測試

根據(jù)不同實驗組的解的精度和算法穩(wěn)定性，實驗分析為兩個部分。使用TSPLIB中標準數(shù)據(jù)集作為實驗數(shù)據(jù)進行測試。實驗環(huán)境為：Python3.7,Inter(R)i7-9750H 2.60 GHz,內(nèi)存為8 GB的個人計算機(personal computer，PC)進行實驗分析。經(jīng)多次試驗總結(jié)，交換全局搜索能力強，逆序局部局部搜索能力強，本文算法需要局部搜索優(yōu)化解，故設(shè)置交換概率設(shè)為0.4，逆序概率設(shè)為0.6。

3.2.1 與經(jīng)典式啟發(fā)算法性能對比

在眾多的經(jīng)典啟發(fā)式算法中選擇了人工蜂群算法[13-14](artificial bee colnony，ABC)。粒子群算法[15](particle swarm optimization，PSO)。蟻群算法[16](ant colony optimization，ACO)三種經(jīng)典算法和本文算法進行實驗對比。選擇這三種經(jīng)典算法進行比較，蟻群算法[17]和粒子群算法[18]在處理多旅行商問題效果較好，而人工蜂群算法又是最新應(yīng)用于多旅行商問題的新算法，表1為對比算法的參數(shù)設(shè)置，其中參數(shù)設(shè)置根據(jù)文獻中對應(yīng)算法進行設(shè)置。

在實驗過程中每個旅行商必須訪問城市，即單個旅行商訪問的城市總和大于1，四種算法每個實例獨立運行15次的實驗結(jié)果如表2所示。

從表2中可知，本文算法和經(jīng)典的三種算法的比較來看，在所測試數(shù)據(jù)集上，旅行商的個數(shù)分別為3、5、8三種情況下，本文算法得到的結(jié)果都比對比算法要好。同時引入對比參數(shù)PAB(bercentage average best)，PAB是平均解與最優(yōu)解的差值與最優(yōu)解的比值，PAB小表示平均解和最優(yōu)解之差小，波動小表示算法穩(wěn)定。可以看到本文算法的PAB集中在2左右變化，說明算法在不同測試集的上的性能較為穩(wěn)定。四種算法PAB變化曲線如圖11所示。

圖7 Lin318分組效果圖Fig.7 Lin318 group renderings

表1 算法參數(shù)設(shè)置Table 1 Algorithm parameter settings

表2 經(jīng)典算法效果對比Table 2 Comparison of effects of classical algorithms

從圖8中可以看出，與經(jīng)典的三種算法比較可知，人工蜂群算法和粒子群算法波動過大，算法性能不穩(wěn)定，本文算法在區(qū)間[2,4]波動，算法平均解和最優(yōu)解之差優(yōu)于其他三種算法，且幅度較小算法穩(wěn)定性高。結(jié)合表2，本文算法解質(zhì)量好，且針對不同數(shù)據(jù)集結(jié)果穩(wěn)定。

圖8 四種算法穩(wěn)定性分析Fig.8 Stability analysis of four algorithms

3.2.2 與近年文獻改進的算法實驗對比

在近年文獻中選取改進效果較好的灰狼優(yōu)化算法[19]和雜草入侵算法[20-21]與本文算法進行實驗結(jié)果對比如表3所示。

在和近年文獻對比的過程中，旅行商個數(shù)m進行了增加，分別對m為5、8、10進行測試，所得到的效果本文算法依舊為最優(yōu)解，改進的灰狼優(yōu)化算法與本文求得的結(jié)果相差近2倍，雜草入侵算法更是相差更多，關(guān)于算法求解穩(wěn)定性方面，只有n=150、m=10和n=51、m=10產(chǎn)生波動，為了更形象地展示對比算法的穩(wěn)定性，將得到的結(jié)果進行了可視化，圖9為PAB的變化曲線。

圖9 三種對比算法的穩(wěn)定性分析Fig.9 Stability analysis of three contrast algorithms

從圖9可以明顯地看出，本文算法整體的變化曲線較為平緩，PAB小于其他兩種算法，平均解和最優(yōu)解之差小。如表3所示本文算法在三種對比算法中求得最優(yōu)解的同時，對于每個數(shù)據(jù)集不同的旅行商個數(shù)也較為穩(wěn)定。

針對算法的性能,選取較有公共數(shù)據(jù)集a280和lin318進行測試，在n=280、m=10和n=318、m=10時四種對比算法運行10次的結(jié)果變化曲線，由于經(jīng)典的算法在求解過程中與本文算法相差甚大，根據(jù)參考文獻[19]的數(shù)據(jù)對比，此處選擇近年來最新改進效果較好的算法進行比較，結(jié)果圖10所示。

從圖10中可以清晰地看到，在n=280、m=10的情況下，求解精度方面本文算法是最好的，運行10次折線圖波動的幅度較小，算法求解較為穩(wěn)定。對比算法AC-PGA(ant colony partheno genetic algorithms)在四種算法中求解質(zhì)量和穩(wěn)定性都比較差，STASA-2opt求解效果有個別點比較好，但是整體波動的幅度較大，結(jié)果具有隨機性。加入2opt的STASA算法運行10次的效果較為穩(wěn)定，但是求解能力不如本文算法。在n=318、m=10情況下，STASA(state transition simulated annealing algorithm)求解的結(jié)果相差太大，n=318時求解效果相對于在n=280時大幅降低，旅行商個數(shù)不變，城市增加40個，算法性能驟減，該算法不適合求解大規(guī)模問題。由此可知，本文算法在求解精度和穩(wěn)定性方面要優(yōu)于對比算法。圖11分別是不同數(shù)據(jù)集求得最優(yōu)路徑，分類清晰且無交叉，求解的結(jié)果是對比算法中的最優(yōu)解。

表3 改進算法效果對比Table 3 Comparison of improved algorithm effects

圖10 四種算法的效果對比Fig.10 Comparison of the four algorithms

圖11 不同數(shù)據(jù)集的最優(yōu)路徑效果Fig.11 The optimal path renderings of different sets

4 總結(jié)與展望

多旅行商問題具有強大的實用背景和理論價值，提出一種基于改進K-means++聚類的信息傳播算法求解多起點的多旅行商問題(MMTSP)。采用兩段式方法解決多旅行商問題，采用聚類算法改進的K-means++進行分組，其算法本身簡單，聚類時間較快，對大規(guī)模的數(shù)據(jù)集運行效率較高且具有伸縮性，改進的K-means++算法中k值根據(jù)旅行商的個數(shù)進行確定較為簡單。分組之后采用改進的信息傳播算法進行優(yōu)化，信息傳播算法在圖模型上進行迭代，在解決的過程中利用因子圖的特性和BP算法的傳播性，隨機選擇開始城市，在因子圖上進行迭代運算。實驗表明在城市數(shù)目較少的情況下，信息傳播算法可以在較少的迭代次數(shù)中直接求得最優(yōu)解。但是隨著城市數(shù)目的增加，因子圖會變的復(fù)雜從而影響解決問題的效率，加入局部搜索算法可以很好的解決這一問題。經(jīng)多數(shù)據(jù)集測試和多種算法比較，本文算法求得最優(yōu)解和平均解效果優(yōu)于其他算法，PAB更小，并且波動區(qū)間小算法穩(wěn)定。但在旅行商問題規(guī)模較大時，城市無向圖較為復(fù)雜，因子圖更加龐大，計算量進而增加，在個別旅行商問題上，精度會有略微下降，這是算法的不足之處。如何在問題規(guī)模較大時，縮減因子圖的規(guī)模這將是下一步工作首要解決的問題。