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

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

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

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

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

1 基本知識

1.1 多旅行商問題

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

(1)

1.2 K-means聚類方法

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

算法的基本步驟如下。

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

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

步驟3聚類前后樣本點不再發生變化，聚類達到穩態，則算法終止，否則進入步驟4。

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

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

1.3 因子圖

因子圖[7]是一個特殊的二分圖，代表由許多變量組成的全局函數分解成局部函數的圖示。信息傳播算法以因子圖為工具，可以直觀地表現一個多種約束條件的復雜問題。如函數p(x)可以分解成fA、fB、fC、fD、fE5個局部函數的乘積，構造對應關系為式(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年人工智能領域著名學者Pearl提出了置信傳播(belief propagation, BP)算法，BP算法把全局的概率推理過程轉變為局部變量間的消息傳遞，從而降低了推理的復雜度。因為BP算法的優點，自提出起，就受到了國際上眾多學者的關注，掀起了研究的熱潮。

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

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

(4)

(5)

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

(6)

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

2.1 旅行商問題的線性規劃方程

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

minimizewTx

xe∈[0,1]|E|

(7)

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

2.2 旅行商問題的因子圖

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

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

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

(8)

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

(9)

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

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

2.3 旅行商問題的算法優化

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

(10)

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

(11)

2.4 模擬退火算子

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

其中模擬退火算子是以一定的概率接受新的解，當在溫度T時，從當前解curr改變為新解new；若curr

2.5 局部搜索操作

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

2.5.1 交換操作及示例

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

2.5.2 逆序操作及示例

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

圖4 交換操作Fig.4 Exchange operation

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

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

算法1:求解旅行商問題的信息傳播算法輸入: 圖G=(V,E),加權鄰接矩陣A,最大迭代Tmax,閾值εmax輸出:巡視中路徑Tour?E1 構造初始化因子圖2 初始化勢函數信息^μ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 數值實驗及分析

3.1 K-means聚類分組分析

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

聚類方法的偽代碼如下。

算法2:K-means質心選擇輸入:數據集D,聚類個數K,鏈長M輸出:K個質心的矩陣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)) #根據概率隨機選擇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++的聚類效果，發現改進后避免了隨機化中心點對算法造成的影響。在聚類效果上進行對比，對比效果如圖6、圖7所示。

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

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

3.2 算法性能測試

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

3.2.1 與經典式啟發算法性能對比

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

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

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

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

表1 算法參數設置Table 1 Algorithm parameter settings

表2 經典算法效果對比Table 2 Comparison of effects of classical algorithms

從圖8中可以看出，與經典的三種算法比較可知，人工蜂群算法和粒子群算法波動過大，算法性能不穩定，本文算法在區間[2,4]波動，算法平均解和最優解之差優于其他三種算法，且幅度較小算法穩定性高。結合表2，本文算法解質量好，且針對不同數據集結果穩定。

圖8 四種算法穩定性分析Fig.8 Stability analysis of four algorithms

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

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

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

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

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

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

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

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

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

圖11 不同數據集的最優路徑效果Fig.11 The optimal path renderings of different sets

4 總結與展望

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