混合算法求解著色瓶頸旅行商問題

董學士 董文永 蔡永樂

(武漢大學計算機學院 武漢 430072) (dxs_cs@163.com)

著色旅行商問題(colored traveling salesman problem, CTSP)是旅行商問題(traveling salesman problem, TSP)與多旅行商問題(multiple traveling salesman problem, MTSP)的一種擴展模型，CTSP是一種來源于但不局限于多機器工程系統(multi-machine engineering system, MES)的模型，可應用在具有部分重合工作區域的規劃[1].為建模有部分重合區域的人員與車輛調配的路線優化問題，本文給出了一種新的模型——著色瓶頸旅行商問題(colored bottleneck traveling salesman problem, CBTSP)，可應用于有合作與單獨任務的人員與車輛的路線規劃等問題，其目標是最小化所有旅行路線的最大邊，具體的應用場景參見1.3節的論述.在智能交通、多任務協作等領域，一些實際問題可用CBTSP來建模，所構建模型的尺度(對應著CBTSP的城市數量)往往趨于大尺度(城市數量1 000或以上)，因此，研究大尺度的CBTSP及其求解算法有一定的意義.

由于CBTSP是本文提出的模型，尚沒有相應的求解算法，算法方面的研究主要集中在CTSP模型.東南大學Li等人[2]提出CTSP模型，并用遺傳算法求解該問題；之后，Li等人[1]將貪心遺傳算法(genetic algorithm with greedy initialization, GAG)、爬山法遺傳算法(hill-climbing genetic algorithm, HCGA)與模擬退火遺傳算法(simulated annealing genetic algorithm, SAGA)應用在求解小規模的CTSP，即CTSP的城市數量小于等于101.瓶頸旅行商(bottleneck traveling salesman problem, BTSP)[3-10]相關工作有：Vairaktarakis[3]給出了 BTSP是NP完全問題，可應用在工作流的規劃領域；BTSP也可應用在重構相鄰信息的排序序列[4]；Garfinkel等人[5]將BTSP應用在集成線路的排序；BTSP另一種應用就是最小化狀態變化機的排序[7]；BTSP也可以應用于最小化機器的最大變化狀態[8]；Ahmed[10]應用遺傳算法求解BTSP.元啟發式算法的其他相關應用：Liao等人[11]將蟻群算法應用于混合可變的優化問題；Ferreira[12]將一種蟻群優化方法用于眼底圖像滲出物的分割.

本文在提出CBTSP模型的基礎上，進一步設計了一種混合粒子群優化算法(particle swarm optimi-zation, PSO)、模擬退火(simulated annealing, SA)和遺傳算法(genetic algorithm, GA)的啟發式算法PSGA來求解該問題.PSGA算法首先用二重染色體編碼來產生問題的解，然后用遺傳算法的交叉操作來更新該解，在求解過程中，活動強度控制著交叉長度，并受粒子半徑和環境溫度的影響.通過使用小尺度到大尺度CBTSP實例的實驗與分析表明：PSGA求解CBTSP是有效的，該算法的求解質量要好于對比算法.

本文的創新點包括：給出了一種問題模型CBTSP，該模型可對一類組合優化問題進行建模，并將該模型的規模擴展到大尺度；提出了一種新的混合算法PSGA用于求解CBTSP問題.

1 著色瓶頸旅行商問題

1.1 CBTSP

著色旅行商問題(CTSP)[1]有m個旅行者和n個城市，其中m,n∈+={1,2,…},且m

在CTSP中，共享數據U有相關定義，U可被多個旅行者訪問，?a∈U，c(a)=Zm.如果di=0∈U且?a∈U，c(a)=Zm，該對應的CTSP整數編碼模型：變量xijk(i≠j,i,j∈V,k∈Zm)表示第k個旅行者是否經過城市i到j，變量uik(i∈V,k∈Zm)表示第k個旅行者從起點到城市i的城市數目[1].

CBTSP的定義與CTSP的定義基本相同，但兩者有不同的目標函數.

CTSP的目標函數[1]：

(1)

CBTSP的目標是最小化旅行回路的最大邊，其表達式為minf=maxedge(i,j)(i,j=0,1,…,n-1).

CBTSP與CTSP雖有不同的目標函數，但有相同的限制條件.

CBTSP的限制條件為

(2)

(3)

(4)

(5)

(6)

式(2)表示每個旅行者開始和結束該起點(起始點 0)；式(3)表示旅行者k不能訪問其他旅行者的獨享城市，另外，其他的旅行者也不能訪問旅行者k的獨享城市；式(4)表示旅行者l(l≠k)既不是開始于第k個城市數據集也不返回到該集；式(5)表示除了起始點0之外，每個城市必須只被訪問一次；式(6)表示每個旅行者必須同時訪問和退出共享城市集.

1.2 CBTSP與CTSP的比較

1) 相同點.CBTSP和CTSP一樣，都有獨享城市集和共享城市集、多個旅行者和多個城市，其中獨享城市只能被指定的旅行者訪問，共享城市可以被所有的旅行者訪問，每個城市只能被訪問一次.2個問題都是NP難問題，兩者可共用相同的數據集.

2) 不同點.2個問題的目標函數不同，CBTSP的目標是使所有旅行回路的最大的邊最小化，CTSP的目標是使所有旅行回路的總的距離最小.

CTSP可以較容易地轉化為CBTSP，改變CTSP的目標函數即可.BTSP轉化為CBTSP，需要經過更為復雜的過程：1)由BTSP單個旅行商變為CBTSP的多個旅行商；2)將BTSP的數據劃分成獨享城市集和共享城市集；3)增加式(2)～(6)的限制條件.

1.3 CBTSP相關應用

CTSP可以應用在多機器工程系統MES的規劃問題.CBTSP可用來建模一類具有協作與單獨任務的人員和車輛的路線規劃優化問題，該類問題一般具有一些特點：它們通常有多個人或車輛，個體不僅需要執行共同區域的合作任務，而且需要執行個體所屬區域的獨立任務，在行動過程中，因某些需要或原因，個體需要被增員，為了能被更容易地找到，需最小化旅行回路的最大邊.該問題的目標、個人和任務可與CBTSP的目標、旅行者和任務相匹配，因此可以用CBTSP來建模.

最大離散旅行商問題(maximum scatter traveling salesperson problem, MSTSP)[13-15]的目標是最大化旅行回路的最小邊，MSTSP是根據制造業與醫學圖像等方面的應用需要而提出的一種模型.例如MSTSP的一個應用實例[13]：假如有一個被冤枉的犯人，沒有違法卻可能會面臨著指控，于是他從警察局越獄逃跑，為了避免被拒捕他開始在全國流竄，當他的一些通緝信息如名字與照片在電視上播放后，一般來說，他在一個地方逗留的時間越長，就會越容易引起當地人的懷疑，有的人可能會到警察局報案，在相信他無辜的全國范圍朋友們的幫助下，他會采取從一個安全的地方到另外一個安全的地方的方式來逃跑，并盡可能地遠離之前安全的地方.用更準確的術語來表達，他需要尋找一個穿過安全地方的旅行回路，即任意相鄰2點的最小距離需盡可能最大化，這個問題的模型稱之為MSTSP.與MSTSP相反，BTSP是最小化旅行回路的最大邊長，使最大的相鄰點距離盡可能最小.如果上面的被誤指控的犯人的實例用BTSP建模，該犯人將會更容易被警察抓到.

本文給出了CBTSP的一類應用，它可以應用在具有合作與單獨任務的人員和車輛路線的規劃.以軍隊(裝甲車)路線規劃為例，有1支軍隊需執行殲敵或其他任務，譬如它可被分成6組，每組包括幾個人(裝甲車)，一方面，這6組需要在自己獨有的區域內執行單獨的任務；另一方面他們需要在共同的或共享的區域內互相合作來完成合作或協作的任務.如果在一個區域內有太多敵人或任務，這6組需要在這個共享的區域相互合作來完成任務，如果在某些區域敵人或任務非常少，每個組將負責一個獨立的區域來獨立完成自己所屬區域的單獨任務.因此，這6組不僅要在共同區域內執行合作的任務，而且需要在獨有區域內執行獨立任務.在共享的區域與獨立區域，分別有許多地點，對于共享區域，其地點可被所有組訪問，并且一個地點只能被訪問一次，一個地點的任務被完成之后，該地點無需再次被訪問，對于獨有的區域，每個地點只能被指定的組訪問.在執行任務的過程中，因為某些原因，例如遇到太多的敵人或新的作戰任務，需要增員或加派人員到某個組，在這幾組執行任務的行程過程中，為了能被增援人員更容易找到，他們需要在整個旅行路線中尋找一個旅行回路，使得旅行路線的最大邊盡可能最小，這個問題就可以通過CBTSP來建模.CBTSP的應用不僅僅局限于上面實例的路線規劃問題，它還可應用在一類具有協作與單獨任務的人員和車輛的路線規劃問題，該問題具有一些特征：它們一般包括多個人員與車輛，每個個體(人員或車輛)不僅需要執行共同區域或共享區域的協作任務，而且需要執行自己區域的獨立任務，在執行任務的行程過程當中，因為某些原因，單獨的個體需要增加人員或支援，為了能被其他人更容易地找到，需最小化旅行路線的最大邊.

1.4 CBTSP理論證明

定理1. CBTSP是NP難問題.

證明. 根據CBTSP的定義可知，CBTSP可由CTSP變化而來，改變CTSP的目標函數，也就是使所有旅行回路中最大的邊最小化，即可轉化為CBTSP.CBTSP可理解成有相同起始點的含部分重合區域(共享城市集)的多個BTSP問題的組合模型，因為BTSP是NP難問題，所以CBTSP也屬于NP難問題.從另一角度來講，CBTSP可由CTSP通過改變目標函數變化而來，因為CTSP屬于NP難問題[1]，而改變目標函數操作不會引起CBTSP的時間復雜度的變化，因此，CBTSP也屬于NP難問題.

2 混合算法求解CBTSP

2.1 解的表示

Fig. 1 Example of dual-chromosome coding圖1 二重染色體編碼的實例

文獻[1]使用二重染色體對遺傳算法求解CTSP的解編碼，本文也采用該編碼方法對CBTSP進行編碼.如圖1所示，城市染色體是n-1個城市的排列，而旅行商染色體是城市染色體的城市所對應的旅行商的排列.

我們給出一個實例來說明該編碼方法，如圖1所示，例如一個11個城市和2個旅行商的CBTSP，城市1～4為旅行商1的獨享城市集，城市5～8為旅行商2的獨享城市集，城市9～11為2個旅行商的共享城市集.如果起始點城市0也加入在內，則旅行商1的訪問路徑為0→2→1→9→3→4→0，旅行商2的訪問路徑為0→10→6→7→5→8→11→0.

2.2 算法步驟

本文提出的算法是一種基于伊藤過程的混合算法，粒子的活動強度受粒子半徑和環境溫度的影響，半徑越小、溫度越高則粒子的運動強度越劇烈，反之，運動強度越弱.粒子的運動強度控制著交叉操作中的交叉長度，強度越大，交叉長度越長.在求解過程中，將第i個粒子與最優解進行交叉，產生新的解.

PSGA求解CBTSP的步驟如算法1所示.

算法1. PSGA求解CBTSP的步驟.

① 參數初始化；

② 初始化M個粒子和環境溫度T；

③ 讀取CBTSP數據；

④ 計算粒子的適應度值，并保存最優值；

⑤ 用式(8)計算粒子半徑；

⑥ 通過式(9)計算環境溫度；

⑦ 用式(10)計算活動強度；

⑧ 根據適應度值，對粒子進行分類；

⑨ 用式(11)計算交叉長度；

⑩ while判斷終止條件是否滿足do

在算法1中，步驟①為算法參數的初始化；步驟②表示初始化種群M和初始溫度T；步驟③讀取CBTSP數據；步驟④計算粒子的適應度值，并保存最優值；步驟⑤～⑦計算粒子半徑、環境溫度和活動強度；步驟⑧為根據適應度值對粒子進行分類；步驟⑨根據式(11)計算粒子的交叉長度；步驟⑩～執行交叉操作；步驟對環境溫度、交叉長度更新.

2.3 混合算法

PSGA算法有4個部分：粒子半徑、環境溫度、活動強度、交叉算子.該混合算法具體設計如下：

Fig. 2 Example of crossover operator圖2 交叉算子的實例

1) 粒子半徑

根據愛因斯坦等人的大分子的運動定律，粒子半徑越小，環境溫度越高，粒子的隨機性運動就越激烈；相反，當粒子半徑越大和溫度越低，隨機運動會越弱.對于優化算法，擁有差的適應度值的粒子，其運動越劇烈；反之，有好的適應度值的粒子運動強度越弱.適應度值越好，越會降低運動強度保持當前解，因此，可以增加大的半徑；相反，可以使半徑變小，考慮到粒子的優缺點，動態調整粒子半徑可計算為[16-19]：

r(x)=g(f(x)),

(7)

其中，x表示當前種群的一個粒子，f(x)代表適應度值，g(x)表示單調函數.因為CBTSP是求最小化的問題，所以g(x)是一個單調遞減函數，該函數有許多選擇，包括線性轉換函數和基于分類的方法等.本文采用分類的方法來計算粒子半徑，首先，N個當前種群的粒子根據適應度值的等級按照最好到最差的順序分類，結果用x1,x2,…,xN來表示.粒子半徑計算為

(8)

其中，rmax和rmin分別表示最大和最小的半徑，所有的半徑被均勻分布在rmax和rmin之間.如果是缺省值，rmax=1，rmin=0.

2) 環境溫度

宏觀世界中，分子運動的劇烈程度是受外界環境溫度影響.當溫度越高，粒子運動越劇烈;當溫度越低，粒子運動越弱.模擬退火算法，在迭代過程中，環境溫度會被逐漸地降低，因此，環境溫度可計算為

Ti=ρ×Ti-1,

(9)

其中，Ti表示第i個規劃時間的溫度；ρ表示退火系數，它能控制溫度降低的速度，缺省值的情況下ρ=0.9.

3) 活動強度

活動強度控制著粒子的運動強度，根據愛因斯坦的分子運動和粒子熱力學運動定律，該強度與粒子半徑成反比、與溫度成正比，因此，當前種群微粒xi的活動強度δ計算為[16-19]

(10)

其中，δ表示粒子的活動強度，ri表示粒子xi的半徑，T表示溫度.

4) 交叉算子

因為交叉操作，粒子可以在解的空間中隨機變化，在環境影響下，將粒子與最優解進行交叉并產生新的解，從而避免陷入局部最優，保持解的多樣性.交叉操作有2部分：交叉長度和交叉過程.交叉長度由活動強度決定，通過環境溫度和粒子半徑計算，即：

li=γ×δi×n,

(11)

其中，li表示第i個粒子的交叉長度；γ為隨機服從均勻分布數，取值在0～1之間;n為二重染色體長度.

如圖2所示，混合算法PSGA的交叉過程的計算步驟如下：

步驟1. 在城市染色體的[1,n-li]之間隨機選擇一個起始位置，將灰色區域的連續li位置與最優解交叉.

步驟2. 用最優解中對應的城市替換灰色區域相應的城市.

例如城市替換的對應關系為：2—3，7—1，1—4，9—7，5—8，3—6，8—2.

步驟3. 根據交換規則，對第i個粒子的灰色區域之外的冗余城市進行替換，直到沒有冗余城市為止.

例如步驟2的灰色區域外的城市染色體4，灰色區域內也存在4，即為冗余，所以要對灰色區域外的4進行替換，根據交換關系，最優解中的4對應第i個粒子的1，于是4被1替換，但1仍然與灰色區域的城市染色體重復，繼續根據對應關系進行替換，最優解中的1對應著第i個粒子的7，但仍不滿足條件，7繼續被9替換，結果滿足條件，因此該操作的最終值為9.

步驟4. 對旅行商染色體錯誤的賦值進行校正.例如，在步驟3中，城市染色體中的1，7，6，2分別對應錯誤的旅行商染色體，因此需要校正，修正之后的結果，如步驟4所示.

通過以上4個步驟，粒子完成交叉過程，實現了對問題的解更新.

3 實驗與分析

實驗的運行環境為：AMD AthlonTMⅡ X4 640 3.01 GHz處理器和3.25 GB內存.實驗是用Java進行設計與開發.

表1為CBTSP的小尺度實驗數據.

在表1中，n的取值范圍是21～76，m的取值范圍是2～6，s表示共享城市的數量，e表示獨有城市的數量.

表1和表2中，城市數量n的取值為21～101的數據來自CTSP的論文[1].

表2和表3為中尺度和大尺度CBTSP的實驗數據.在表2中，n的取值為101～665，m的取值為4～33;在表3中，n的取值為1 379～2 361，m的取值為3～50.表3是根據CBTSP的要求，由TSPLIB TSP數據制作而成，TSPLIB TSP標準數據集已在網上共享.

Table 1 Small Scale Experiment Data for CBTSP表1 CBTSP的小尺度實驗數據

Table 2 Medium Scale Experiment Data for CBTSP表2 CBTSP的中尺度實驗數據

Table 3 Large Scale Experiment Data for CBTSP表3 CBTSP的大尺度實驗數據

混合算法參數的設置情況為：種群M=150，初始溫度T=1 000.遺傳算法和改進遺傳算法的參數設置情況為：種群大小為30，交叉概率0.7，變異概率0.1；模擬退火的遺傳算法參數情況，初始溫度100，總的冷卻時間60，每個溫度的步長為30，模擬系數0.9.以上算法的運行時間相同，即為算法的終止條件，其中算法求解小尺度、中尺度和大尺度CTSP的運行時間分別為20 s，40 s，60 s.

圖3所示為當n=31與m=4時，算法GAG，HCGA，SAGA，PSGA求解CBTSP的運行界面.

Fig. 3 Optimal routes of GAG, HCGA, SAGA and PSGA for CBTSP with n=31 and m=4圖3 當n=31和m=4時GAG，HCGA，SAGA，PSGA求解CBTSP最優路線圖

圖4所示為當n=51與m=5時，算法GAG，HCGA，SAGA，PSGA求解CBTSP的運行界面.GAG:最優解24.0，平均解27.8，迭代次數13 278；HCGA：最優解17.0，平均解23.0，平均求解最優解的迭代次數18 198；SAGA：最優解20.0，平均求解最優解的迭代次數22.4，平均求解最優解的迭代次數1 386；PSGA：最優解20.0，平均解21.8，迭代次數6 088.

圖5為當n=101與m=7時，4個算法求解CBTSP的運行界面.GAG：最優解25.0，最差解35.0，平均解31.1，平均求解最優解的迭代次數14 349；HCGA：最優解21.0，最差解34.0，平均解是27.7，平均求解最優解的迭代次數19 842；SAGA：最優解22.0，最差解27.0，平均解24.9，平均求解最優解的迭代次數2 378；PSGA：最優解21.0，最差解25.0，平均解22.5，平均求解最優解的迭代次數5 720.從該案例數據分析可看出，PSGA的求解質量最好.

表4～6為GA，GAG，HCGA這3種算法求解小尺度、中尺度和大尺度CBTSP的實驗結果對比數據.

表4～6中,n表示城市數量;m表示旅行者數目;best，worst，average為算法求解該問題運行10次的最優解、最差解和平均解;iterations表示算法運行10次的平均求解最優解的迭代次數.從表4～6中可看出，算法HCGA求解CBTSP的求解質量要好于算法GAG和GA.

Fig. 4 Optimal routes of GAG, HCGA, SAGA and PSGA for CBTSP with n=51 and m=5圖4 當n=51和m=5時GAG，HCGA，SAGA，PSGA求解CBTSP最優路線圖

InstancesnmGASolution Quality∕kmBestWorstAverageIterationsGAGSolution Quality∕kmBestWorstAverageIterationsHCGASolution Quality∕kmBestWorstAverageIterations121210.014.012.05711910.011.010.44114710.011.010.432428221311.016.012.14937211.012.011.52123911.014.011.515663331215.021.018.34784615.018.016.24051014.018.015.629989431315.018.016.94680215.021.016.73926815.018.016.418475531416.019.016.93048715.017.015.62437515.016.015.38066641220.027.023.23021116.022.018.95259716.021.018.045525741322.030.023.72887617.030.020.85280617.022.020.524136841419.026.022.93076722.030.023.8572817.023.019.425535951323.032.028.01535519.027.024.92857019.023.021.6221251051527.031.029.2885724.031.027.81327817.027.023.0181981176329.039.035.4936422.036.026.01483621.026.022.7242601276430.037.033.5817823.037.027.71263618.033.024.7136481376527.035.031.61067323.034.028.61360721.028.026.198621476629.040.032.91286626.038.030.91903621.034.026.217373

Notes： Bold number stands for the best values of best solution or average solution in the table.

Fig. 5 Optimal routes of GAG, HCGA, SAGA and PSGA for CBTSP with n=101 and m=7圖5 當n=101與m=7時 GAG，HCGA，SAGA，PSGA求解CBTSP最優路線圖

Notes： Bold number stands for the best values of best solution or average solution in the table.

Table 6 Experimental Results of GA, GAG and HCGA for Large Scale CBTSP表6 GA，GAG，HCGA求解大尺度CBTSP的實驗結果數據

Notes： Bold number stands for the best values of best solution or average solution in the table.

表7～9為算法SAGA和PSGA求解CBTSP的數據.表7～9中，best，worst，average分別為算法求解CBTSP運行10次的最優解、最差解和平均解;iterations為算法運行10次的平均求解最優解的迭代次數.從表7～9中數據可以看出，PSGA求解該問題的求解質量好于SAGA.

圖6和圖7為5種算法求解該問題的平均求解質量對比圖，數據來自于表4～9.圖6和圖7中X軸表示實例的序號，每個序號對應著具體實例數據；Y軸表示算法的平均求解質量.從圖6～7可以看出，PSGA求解CBTSP的解的質量要好于GA，GAG，HCGA，SAGA.

Table 7 Experimental Results of SAGA and PSGA for Small Scale CBTSP表7 SAGA與PSGA求解小尺度CBTSP的實驗結果數據

Notes： Bold number stands for the best values of best solution or average solution in the table.

Table 8 Experimental Results of SAGA and PSGA for Medium Scale CBTSP表8 SAGA與PSGA求解中尺度CBTSP的實驗結果數據

Notes： Bold number stands for the best values of best solution or average solution in the table.

Table 9 Experimental Results of SAGA and PSGA for Large Scale CBTSP表9 SAGA與PSGA求解大尺度CBTSP的實驗結果數據

Notes： Bold number stands for the best values of best solution or average solution in the table.

Fig. 6 Average solution quality of the algorithms for medium scale CBTSP圖6 算法求解中尺度CBTSP的平均求解質量對比圖

Fig. 7 Average solution quality of the algorithms for large scale CBTSP圖7 算法求解大尺度CBTSP的平均求解質量對比圖

為了驗證算法的有效性，引入文獻[20]的百分偏差對算法的求解質量進行對比分析.表10為5種算法求解3種不同尺度的CBTSP的百分偏差,PDbest表示最優解百分偏差，PDav表示平均解百分偏差.從表10可以看出，PSGA的PDbest與PDav的值在5種算法中最小，表明PSGA求解該問題的最優解和平均解要優于SAGA，HCGA，GAG，GA.

Table 10 Percentage Deviation of the Algorithms for Different Scale CBTSP Problem表10 算法求解不同尺度CBTSP的百分偏差數據

Notes： Bold number stands for the best values of best solution or average solution in the table.

PSGA是先進的群體智能算法，在求解組合優化等問題方面可表現出較好的性能，從本節分析可看出，PSGA求解CBTSP要好于其他算法.PSGA用二重染色體構建問題的解，并用交叉算子對問題的解動態更新，具有較強的求解問題的能力.為驗證本文所提混合算法的性能，使用小尺度、中尺度和大尺度的CBTSP實例數據進行實驗與對比分析，結果表明PSGA求解CBTSP是有效的，其求解質量要優于SAGA，HCGA，GAG，GA算法.

4 結語與展望

針對實際應用需求，本文給出一種CBTSP模型，可應用在具有部分重合工作區域的人員與車輛的路線規劃問題，為拓展該模型的應用領域，開發與研究高性能的求解算法具有一定的意義.本文提出一種新的混合算法，并將其應用在求解CBTSP之中，通過3種不同尺度的實驗表明：PSGA算法求解CBTSP解的質量要優于相關的對比算法.

今后的可能工作包括2方面：1)繼續研究與探索CBTSP的應用領域，開發更先進的算法和模型來求解該問題，使求解質量與速度得到進一步的提升；2)本文研究CBTSP的規模或尺度有限，以后的工作可集中在探討各種算法在求解更大尺度CBTSP的特點和規律.