最大最小目標的多旅行商問題求解①

袁 志

(廣州大學 華軟軟件學院, 廣州 510990)

1 概述

多旅行問題(Multiple Traveling Salesmen Problem,MTSP)是TSP的擴展. 給定一個中心城市和n個訪問城市, 將訪問城市分配給m個旅行商, 每個旅行商從中心城市出發巡游若干訪問城市后回到中心城市, 要求所有旅行商經過的總路程(total)盡量小且其中最長環路的長度(max)盡量小. 事實上, 無法保證total和max兩個目標同時達到最小, 本文尋求最小化max, 稱之為MinMax-MTSP. 這一類問題的算法可應用在工作均衡調度, 印刷機調度、衛星測量系統設計、機器人應急響應[1–3]等領域.

城市集合C ={c0,c1,c2,…,cn}, c0表示中心城市, c1,c2,…,cn表示訪問城市. 距離矩陣D={di,j}, di,j為ci到cj的距離. 每條環路用城市號序列來編碼, 例如cycle =(0,1,2)表示 c0→c1→c2→c0. 圖1顯示了 n=10, m=3 的MinMax-MTSP的一個解(非最優解), 其中的3條環路分別是(0,1,2), (0,5,4,3)和(0,6,7,8,9,10), 解表示為S=(0,1,2,0,5,3,4,0,6,7,8,9,10). MinMax-MTSP求解目標是尋找S的最佳序, 使max最小.

將進化算法與局部搜索算法結合來求解MTSP是近年來的主要方法. 一類是遺傳算法, Carter[4]提出了一種雙基因編碼(城市基因和組基因)的遺傳算法, 并提出了12個測試例; Brown[5]提出一種分組遺傳算法.Yuan[6]在分組遺傳算法中引入了一種新的交叉算子(TCX), 提出了3個新的測試例; Singh[7]改進了分組遺傳算法(GGA-SS), 用組的rk值(rk=路徑長度/城市數)衡量一個組屬于最優解的可能性, rk值較大的組優先保留在子代中, 其余一些零散的城市用貪心算法插入到各個組, 使用2-opt局部搜索算子進行組內優化.另一類是群體智能算法, Liu[8]用蟻群算法求解MTSP;VENKATESH[9]提出了兩種蜂群算法(ABCFC、ABCVC) 和一種雜草入侵算法 (IWO), 同樣用2-opt進行組內優化.

圖1 3個旅行商10個訪問城市的一個解

以上算法的一致性在于: 使用局部搜索算法來快速加速尋優過程, 使用群體進化算法不斷累計優化結果. 這些算法采用2-opt作為局部搜索算子, 而且僅將2-opt作用于單條環路的優化. 實際上, 在MTSP中, 可以將局部搜索從單條環路的優化擴展為兩條環路的重組優化, 加速算法的尋優過程(第2節討論). 這些算法主要依賴個體之間交換信息來完成迭代優化, 很少考慮到局部搜索算子自身的特點, 在進化算法中根據局部搜索算子自身的特點, 設計新的進化機制, 是本文關注的另一個重要內容(在第3節討論).

本文設計了一個新的局部搜索算子reverse/ move(轉置/移動), 該算子既能進行一條環路的優化, 也能重組優化兩條環路, 即便在一條環路上進行優化, 其能力也明顯好于2-opt; 在分析reverse/move算子特點的基礎上, 提出了“搜索-選優-變異-搜索”策略, 設計了競爭搜索算法(Competitive Search Algorithm, CSA), 在文獻[4]和文獻[6]的15個測試例上進行實驗, 與文獻[6–9]進行比較, CSA在計算結果上有明顯的改進.

2 局部搜索算子reverse/move

為了指導搜索過程, 需要對解有合適的評價方法.MinMax-MTSP的目標是最小化max, 但注意到, 將目標改為優先最小化max, 其次最小化total, 有助于最小化max. 因為后者要求每一條環路自身次序最優, 在大多數情況下, 這有利于最小化max.

以圖1為例, 算法運行到當前步驟, (0,6,7,8,9,10)是最長環路, 次序已經是最優; (0,5,3,4)是另外一條稍短的環路, 次序不是最優. 此時對(0,5,3,4)進行優化得到(0,3,4,5), 然后從(0,6,7,8,9,10)中移動“6”到(0,3,4,5)可以得到(0,3,4,5,6), 這是一條新的最長環路,而且比前一條最長環路更短.

因此, 將解的適應值設計為一個二元組(max,total), 規定: 解 S1優于解 S2, 當且僅當 (max1＜max2) 或(max1=max2且 total1＜total2).

本文局部搜索算子通過依次檢查S中的每一個位置來完成. 下文中, 將被檢查位置上的城市標記為c1,它在S中的后一個城市標記為c3, c1的鄰域N(c1)中的城市的標記為c2, c4和c5分別是c2在S中的后一個城市和前一個城市. 按照Lin-kernighan算法的建議,N(c1)取離c1最近的6個城市.

我們的局部搜索方法是: 依次檢查S的每一個位置, 對于該位置上的城市c1, 對N(c1)中每一個c2, 做兩種嘗試: 先嘗試轉置c2與c3之間(包含c2與c3)的次序, 得到 S′, 如果 S′優于 S, 則用 S′代替 S; 否則嘗試將c2移動到 c1和 c3之間, 得到 S′, 如果 S′優于 S, 則用S′代替S. S包含n+m個位置, 如果連續檢查n+m個位置都無法優化 S, 算子停止. 轉置和移動合稱“reverse/move”.

c2和c3可能位于S中的同一環路, 也可能位于兩條不同的環路, 以下分別進行討論.

當c2和c3在同一環路中, 轉置c2和c3之間的次序將刪除兩條邊并新增兩條邊, 效果如圖2(a), 例如:cycle=(0,1,2,3,4,5,6), c2=1, c3=5, reverse(cycle)=(0,5,4,3,2,1,6); 移動c2將刪除三條邊并新增三條邊, 效果如圖2(b), 例如: cycle=(0,1,2,3,4,5,6), c2=1,c3=5, move(cycle)=(0,2,3,4,1,5,6).

當c2和c3屬于分屬兩條環路. 轉置c2和c3之間的次序, 將對兩條環路進行重組, 效果如圖3(a), 例如:S=(0,1,2,0,5,3,4,0,6,7,8,9,10), c2=1, c3=7,reverse(S)=(0,7,6,0,4,3,5,0,2,1,8,9,10). 移動 c2, 相當于將c2從一條環路中移除, 插入到另一條環路, 效果如圖3(b), 例如: S=(0,1,2,0,5,3,4,0,6,7,8,9,10), c2=1, c3=7,move(s)=(0,2,0,5,3,4,0,6,1,7,8,9,10).

圖2 一條環路中的轉置/移動

圖3 跨兩條環路的轉置/移動

我們在經典TSP的公開測試數據TSPLIB上比較reverse/move與2-opt. 經典TSP問題相當于m=1的MTSP問題, 其上的測試結果能夠反映局部搜索算子的尋優能力. 對每個測試例, 分別產生800個隨機的初始解, 然后各用reverse/move與2-opt兩種局部搜索進行優化, 得到800個最終解. 從兩個方面進行比較: 1)環路長度(length)的平均值, 2) 檢查位置數(checkedpoints)的平均值, 結果如表1.

表1 reverse/move與2-opt的搜索能力比較

表1顯示, 對每一個測試例, 兩種局部搜索算法在停止時, reverse/move檢查的位置數略大于2-opt, 與此同時, 在搜索的結果上, reverse/move的明顯好于2-opt.

3 競爭搜索算法

當reverse/move算子停止搜索時, 對所得到的解S, 選擇其中一個片段, 轉置其次序, 然后再執行reverse/move, 可能得到更優的解. 以圖4為例, 圖4(a)是reverse/move得到的一個解, 圖4(b)是一次轉置的結果, 圖4(c)和圖4(d)是再次執行reverse/move的中間過程和結果.

圖4 轉置一個片段后再次執行reverse/move

用一對城市號(cstart,cend)來標記S中的一個片段,這樣的對共有(n+m)×(n+m-1)/2個, 構成S的候選對集. 用隨機方式來選擇一對城市, 轉置其所標記的片段,然后再次執行reverse/move, 可能有三種結果: 1) 得到更好的解, 2) 恢復到轉置之前的解, 3) 得到變差的解.每選出一對城市, 就從S的候選對集中將其刪除, 直到候選對集為空.

根據reverse/move的以上特點, 我們提出一種群體迭代策略: “搜索-選優-變異-搜索”, 用圖5 說明. 圖5中的橫坐標是解空間, 縱坐標是解的適應值. 如果一組解經過局部搜索得到相同的局部最優解, 則將這組解集中在一個區域, 該區域中適應值最小的解就是該區域的局部最優解. 選擇最好的若干個局部最優解, 對其進行變異(隨機的片段轉置), 得到的一部分新的解將到達新的區域, 經再次搜索得到新的局部最優解. 如此迭代, 逐步提高局部最優解的質量, 直至最好的若干個局部最優解的候選對集全部為空.

圖5 多個解經過局部搜索得到相同的局部最優解

基于這一策略, 我們設計了競爭搜索算法(CSA),如算法1.

算法1. CSA算法1) 設定種群規模p和選擇比例θ(θ建議取0.2), 隨機生成p個初始解, 對每個解執行reverse/move.2) 對所有解按適應值從小到大排序, 保留θ×p個排名靠前且互異的解.3) 對保留的解, 從其候選對集中隨機選(1–θ)/θ個城市對, 并從候選對集中刪除這些對; 對保留解, 分別轉置這些片段產生新的解.4) 在新的解上執行reverse/move, 然后加入種群.5) 循環執行2)到4), 直至所有保留解的候選集為空.6) 輸出群體中的最優解.

CSA算法與保留精英的遺傳算法相似, 其特異性在于: 用局部搜索reverse/move取代雜交操作; 用最優且互異的若干個體作為父代個體, 取代概率性選擇; 用固定的小尺度變異(轉置一個片段)替代概率性變異.

4 實驗

我們用文獻[4]的12個測試例和文獻[6]的3個測試例作為實驗數據. 文獻[4] 的12個測試例中的城市數據是二維平面坐標, 包括MTSP-51的3個問題(m=3,m=5,m=10), MTSP-100的4個問題(m=3,m=5,m=10,m=20)和MTSP-150的5個問題(m=3,m=5,m=10,m=20,m=30), 用第一個城市作為中心城市. 文獻[6]的3個測試例中的數據是128個城市的經緯度和距離矩陣, 包括sgb128(m=10,m=15,m=30), 用第一個城市作為出發城市.

我們在2.8 GHz, 2 Core, 4 G RAM的Windows 8.1系統上實現了CSA, 在實驗中種群規模p設置為50. 表2給出了CSA在15個測試例上得到的最小的最長環路長度(best max)、迭代次數(iterations)和計算時間(time).

表2 CSA算法的結果、迭代次數和時間

表2說明, CSA能在較短時間內終止, 滿足實際應用的需求.

對于文獻[4]的12個測試例, CSA所得的best max與文獻[6–9]的結果比較如表3.

表3 幾種算法在文獻[4]測試例上的best max比較

由表3, 對文獻[4]的12個問題, 在所有算法中,CSA的結果都優于或等于現有算法的最好結果, 其中兩個測試例MTSP-150 (m=3,m=5)的解展示在圖6中.

圖6 CSA在文獻[6]問題上的兩個解

對于文獻[6]的3個測試例, CSA所得的best max與文獻[6–9]的結果比較如表4.

由表4, 對文獻[6]的3個測試例, CSA的結果有大幅度改進. 由于文獻[6]提供的128城市的坐標用經緯度來表示, 在二維平面較難展示, 下面給出m=10,max=2748的解, 其中每一對括號包含的是一條環路.

表4 幾種算法在文獻[6]測試例上的best max比較

5 結論

為了求解最大最小目標的多旅行商問題, 在對現有文獻進行研究的基礎上, 提出了競爭搜索算法(CSA), 與近期文獻中的相比, 明顯提高了解的質量. 局部搜索算子和變異算子是CSA算法中的兩個關鍵算子. 我們曾采用多種不同的變異算子, 包括: 隨機移動一個城市、隨機轉置一個片段、隨機轉置兩個或多個片段以及這些操作的組合, 我們觀察到, 采用不同的變異算子, 在收斂速度和最終解質量上有明顯差異, 其中,隨機轉置一個片段的變異方法明顯好于其他方法. 改進變異方法, 可能進一步提高CSA算法性能.