禁忌搜索演算法求解帶時間窗的車輛路徑規劃問題詳解（附Java程式碼）
- 2020 年 1 月 3 日
- 筆記
大家好呀！
眼看這9102年都快要過去了，小編也是越來越感覺著急了：
為什麼感覺自己今年還這麼蔡！
所以趕緊趁考試周來臨前，碼出了這篇禁忌搜索演算法解決VRPTW的文章，臨時抱佛腳，假裝自己今年學了一點東西。
本文附帶Java程式碼詳解，是根據過去學長寫的C++程式碼修改而來的：
乾貨 | 十分鐘掌握禁忌搜索演算法求解帶時間窗的車輛路徑問題(附C++程式碼和詳細程式碼注釋)
新的程式碼加入了原先忘加的藐視準則，將一些冗餘程式碼改為函數調用，同時加入大規模注釋，很適合沒嘗試過VRPTW的同學學習（沒錯就是我自己）。
算例是用上文的格式，所以建議在仔細閱讀本文程式碼之前先瀏覽一下上文。
下面就開始今天的分享吧！
VRPTW簡介
VRPTW問題可描述為：假設一個配送中心為周圍若干個位於不同地理位置、且對貨物送達時間有不相同要求的客戶點提供配送服務。其中，配送中心用於運行的車輛都是同一型號的（即擁有相同的容量、速度）；配送中心對車輛出入的時間有限制。我們的任務是找出使所有車輛行使路徑總和最小的路線。
VRPTW的更多詳細介紹可以參考之前的推文：
乾貨|十分鐘快速掌握CPLEX求解VRPTW數學模型（附JAVA程式碼及CPLEX安裝流程）
為了保持文章的獨立型，同時方便後續講解，這裡給出建模實例（參考文獻在文末標註）：
所求的所有車輛路線需滿足以下要求：
在此基礎上求出每輛車輛的總時間最短（由於車輛速度相同，時間最短相當於路程最短）的路線。（允許不使用某些車輛）
Tabu Search簡介
禁忌搜索演算法（Tabu Search Algorithm，簡稱TS）起源於對於人類記憶功能的模仿，是一種亞啟發式演算法(meta-heuristics)。它從一個初始可行解(initial feasible solution)出發，試探一系列的特定搜索方向(移動)，選擇讓特定的目標函數值提升最多的移動。為了避免陷入局部最優解，禁忌搜索對已經歷過的搜索過程資訊進行記錄，從而指導下一步的搜索方向。
禁忌搜索是人工智慧的一種體現，是局部搜索的一種擴展。禁忌搜索是在鄰域搜索(local search)的基礎上，通過設置禁忌表（tabu list）來禁忌一些曾經執行過的操作，並利用藐視準則來解禁一些優秀的解。
有關禁忌搜索演算法的具體內容可以參考往期推文：
TS+VRPTW
對鄰域搜索類演算法而言，採取的搜索運算元和評價函數至關重要。下面詳細介紹程式碼中針對VRPTW的插入運算元和評價函數。
插入運算元：
評價函數：
演算法概述
Java程式碼詳解
程式碼主要分為以下幾個類：
Main，主函數；
CustomerType，存放客戶節點的資訊；
RouteType，存放車輛路線資訊；
Parameter，存放全局變數；
EvaluateRoute，處理路線方法；
InitAndPrint，初始化與輸出對應方法；
TS，禁忌搜索方法。
接下來分別介紹。
Main
public class Main { public static void main (String arg[]) { long begintime = System.nanoTime(); ReadIn(); Construction(); TabuSearchnewnew(); Output(); long endtime = System.nanoTime(); double usedTime= (endtime - begintime)/(1e9); System.out.println(); System.out.println("程式耗時："+usedTime+"s"); } }
程式的入口。
CustomerType
public class CustomerType { int Number;//節點自身編號 int R;//節點所屬車輛路徑編號 double X, Y;//節點橫縱坐標 double Begin, End, Service;//節點被訪問的最早時間，最晚時間以及服務時長 double Demand;//節點的需求容量 public CustomerType() { this.Number=0; this.R=0; this.Begin =0; this.End=0; this.Service=0; this.X=0; this.Y=0; this.Demand=0; } public CustomerType(CustomerType c1) { this.Number=c1.Number; this.R=c1.R; this.Begin =c1.Begin; this.End=c1.End; this.Service=c1.Service; this.X=c1.X; this.Y=c1.Y; this.Demand=c1.Demand; } }
客戶類，對圖中的每一個客戶，分別構建客戶類，存放自身編號，所屬車輛路線，坐標位置，訪問時間窗，服務所需時長、需求。
RouteType
import java.util.ArrayList; public class RouteType { public double Load;//單條路徑承載量 public double SubT;//單條路徑違反各節點時間窗約束時長總和 public double Dis;//單條路徑總長度 public ArrayList <CustomerType> V=new ArrayList<>(); //單條路徑上顧客節點序列。在route中，第0個、最後一個都為depot，第k個為第k位。 }
路線類，記錄該路線的總承載量，總長度，對時間窗約束的總違反量，以及單條路徑上的客戶節點序列。
Parameter
public class Parameter { public static double INF=Double.MAX_VALUE; public static int CustomerNumber=100;//算例中除倉庫以外的顧客節點個數 public static int VehicleNumber = 25; public static int Capacity=200;//車輛最大容量 public static int IterMax=2000;//搜索迭代次數 public static int TabuTenure=20;//禁忌步長 public static int[][] Tabu=new int[CustomerNumber + 10][VehicleNumber + 10];//禁忌表用于禁忌節點插入操作:[i][j]將節點i插入路徑j中 public static int[] TabuCreate=new int[CustomerNumber + 10];//禁忌表用於緊急拓展新路徑或使用新車輛 public static double Ans;//最優解距離 public static double Alpha = 1, Beta = 1, Sita = 0.5;//Alpha，Beta為係數，計算目標函數值；Sita控制係數改變的速度 public static double[][] Graph=new double[CustomerNumber + 10][CustomerNumber + 10];//記錄圖 public static CustomerType[] customers=new CustomerType[CustomerNumber+10];//存儲客戶數據 public static RouteType[] routes=new RouteType[CustomerNumber+10];//存儲當前解路線數據 public static RouteType[] Route_Ans=new RouteType[CustomerNumber+10];//存儲最優解路線數據 }
參數類，有關VRPTW和TS的變數都存儲在這裡，在這裡修改數據。
EvaluateRoute
public static boolean Check ( RouteType R[] ) {//檢驗解R是否滿足所有約束 double Q = 0; double T = 0; //檢查是否滿足容量約束 for ( int i = 1; i <= VehicleNumber; ++i ) if ( R[i].V.size() > 2 && R[i].Load > Capacity )//對有客戶且超過容量約束的路徑，記錄超過部分 Q = Q + R[i].Load - Capacity; //檢查是否滿足時間窗約束 for ( int i = 1; i <= VehicleNumber; ++i ) T += R[i].SubT; //分別根據約束滿足的情況和控制係數Sita更新Alpha和Beta值 //新路徑滿足條件，懲罰係數減小， //新路徑違反條件，懲罰係數加大。 if ( Q == 0 && Alpha >= 0.001 ) Alpha /= ( 1 + Sita ); else if ( Q != 0 && Alpha <= 2000 ) Alpha *= ( 1 + Sita ); if ( T == 0 && Beta >= 0.001 ) Beta /= ( 1 + Sita ); else if ( T != 0 && Beta <= 2000 ) Beta *= ( 1 + Sita ); if ( T == 0 && Q == 0 ) return true; else return false; }
check函數是對產生解的檢驗。
由於插入運算元產生的解並不都滿足所有約束條件，對局部搜索產生的較優解需要判斷是否滿足時間窗約束和容量約束後，再決定是否為可行解。
在check局部最優解的過程中，修改懲罰係數Alpha、Beta的值。
public static void UpdateSubT(RouteType r) {//更新路徑r對時間窗的違反量 double ArriveTime =0; for ( int j = 1; j < r.V.size(); ++j ) {//對每一個節點分別計算超出時間窗的部分 ArriveTime = ArriveTime + r.V.get(j-1).Service //服務時間 + Graph[r.V.get(j-1).Number][r.V.get(j).Number];//路途經過時間 if ( ArriveTime > r.V.get(j).End )//超過，記錄 r.SubT = r.SubT + ArriveTime - r.V.get(j).End; else if ( ArriveTime < r.V.get(j).Begin )//未達到，等待 ArriveTime = r.V.get(j).Begin; } }
UpdateSubT函數更新一條車輛路線中在每一個客戶點的時間窗違反量。通過遍歷整條路線累加得到結果。
//計算路徑規劃R的目標函數值，通過該目標函數判斷解是否較優 public static double Calculation ( RouteType R[], int Cus, int NewR ) { //目標函數主要由三個部分組成：D路徑總長度（優化目標），Q超出容量約束總量，T超出時間窗約束總量 //目標函數結構為 f(R) = D + Alpha * Q + Beta * T, 第一項為問題最小化目標，後兩項為懲罰部分 //其中Alpha與Beta為變數，分別根據當前解是否滿足兩個約束進行變化，根據每輪迭代得到的解在Check函數中更新 double Q = 0; double T = 0; double D = 0; //計算單條路徑超出容量約束的總量 for ( int i = 1; i <= VehicleNumber; ++i ) if ( R[i].V.size() > 2 && R[i].Load > Capacity ) Q = Q + R[i].Load - Capacity; //計算總超出時間 for ( int i = 1; i <= VehicleNumber; ++i ) T += R[i].SubT; //計算路徑總長度 for ( int i = 1; i <= VehicleNumber; ++i ) D += R[i].Dis; return (D + Alpha * Q + Beta * T);//返回目標函數值 }
Calculate函數計算目標函數值，懲罰部分累加後乘以懲罰係數。
InitAndPrint
//計算圖上各節點間的距離 private static double Distance ( CustomerType C1, CustomerType C2 ) { return sqrt ( ( C1.X - C2.X ) * ( C1.X - C2.X ) + ( C1.Y - C2.Y ) * ( C1.Y - C2.Y ) ); }
根據計算距離。
//讀取數據 public static void ReadIn(){ for(int i=0;i<CustomerNumber+10;i++) { customers[i]=new CustomerType(); routes[i]=new RouteType(); Route_Ans[i]=new RouteType(); } try { Scanner in = new Scanner(new FileReader("c101.txt")); for ( int i = 1; i <= CustomerNumber + 1; ++i ) { customers[i].Number=in.nextInt()+1;//算例文件有點小問題，別介意。。 customers[i].X=in.nextDouble(); customers[i].Y=in.nextDouble(); customers[i].Demand=in.nextDouble(); customers[i].Begin=in.nextDouble(); customers[i].End=in.nextDouble(); customers[i].Service=in.nextDouble(); } in.close(); }catch (FileNotFoundException e) { // File not found System.out.println("File not found!"); System.exit(-1); } for ( int i = 1; i <= VehicleNumber; ++i ) { if ( routes[i].V.size()!=0 ) routes[i].V.clear(); routes[i].V.add ( new CustomerType (customers[1]) );//嘗試往這裡加入一個複製，後面也都要改。 routes[i].V.add ( new CustomerType (customers[1]) ); routes[i].V.get(0).End=routes[i].V.get(0).Begin;//起點 routes[i].V.get(1).Begin=routes[i].V.get(1).End;//終點 //算例中給出節點0有起始時間0和終止時間，所以如上賦值。 routes[i].Load = 0; } Ans = INF; for ( int i = 1; i <= CustomerNumber + 1; ++i ) for ( int j = 1; j <= CustomerNumber + 1; ++j ) Graph[i][j] = Distance ( customers[i], customers[j] ); }
從文件中讀取算例（在此修改算例，記得同時修改Parameter類中的參數），並對當前解routes[ ]的每條路線進行初始化，起終點都為配送中心。
記錄客戶間距離，存儲在Graph數組中。
//構造初始解 public static void Construction() { int[] Customer_Set=new int[CustomerNumber + 10]; for ( int i = 1; i <= CustomerNumber; ++i ) Customer_Set[i] = i + 1; int Sizeof_Customer_Set = CustomerNumber; int Current_Route = 1; //以滿足容量約束為目的的隨機初始化 //即隨機挑選一個節點插入到第m條路徑中，若超過容量約束，則插入第m+1條路徑 //且插入路徑的位置由該路徑上已存在的各節點的最早時間決定 while ( Sizeof_Customer_Set > 0 ) { int K = (int) (random() * Sizeof_Customer_Set + 1); int C = Customer_Set[K]; Customer_Set[K] = Customer_Set[Sizeof_Customer_Set]; Sizeof_Customer_Set--;//將當前服務過的節點賦值為最末節點值,數組容量減1 //隨機提取出一個節點，類似產生亂序隨機序列的程式碼 if ( routes[Current_Route].Load + customers[C].Demand > Capacity ) Current_Route++; //不滿足容量約束，下一條車輛路線 for ( int i = 0; i < routes[Current_Route].V.size() - 1; i++ )//對路徑中每一對節點查找，看是否能插入新節點 if ( ( routes[Current_Route].V.get(i).Begin <= customers[C].Begin ) && ( customers[C].Begin <= routes[Current_Route].V.get(i + 1).Begin ) ) { routes[Current_Route].V.add ( i + 1, new CustomerType (customers[C]) ); //判斷時間窗開始部分，滿足，則加入該節點。 routes[Current_Route].Load += customers[C].Demand; customers[C].R = Current_Route; //更新路徑容量，節點類。 break; } } //初始化計算超過時間窗約束的總量 for ( int i = 1; i <= VehicleNumber; ++i ) { routes[i].SubT = 0; routes[i].Dis = 0; for(int j = 1; j < routes[i].V.size(); ++j) { routes[i].Dis += Graph[routes[i].V.get(j-1).Number][routes[i].V.get(j).Number]; } UpdateSubT(routes[i]); } }
Construction構造初始解。根據前文偽程式碼構造初始解，每次隨機選擇節點（類似打亂有序數列）。
針對該節點找到符合容量約束，同時時間窗開啟時間符合要求的位置，插入該節點。記得在插入節點時同時更新該節點所屬的路徑。
對時間窗違反量進行初始化。
public static void Output () {//結果輸出 System.out.println("************************************************************"); System.out.println("The Minimum Total Distance = "+ Ans); System.out.println("Concrete Schedule of Each Route as Following : "); int M = 0; for ( int i = 1; i <= VehicleNumber; ++i ) if ( Route_Ans[i].V.size() > 2 ) { M++; System.out.print("No." + M + " : "); for ( int j = 0; j < Route_Ans[i].V.size() - 1; ++j ) System.out.print( Route_Ans[i].V.get(j).Number + " -> "); System.out.println( Route_Ans[i].V.get(Route_Ans[i].V.size() - 1).Number); } System.out.println("************************************************************"); }
輸出結果，不再贅述。
public static void CheckAns() { //檢驗距離計算是否正確 double Check_Ans = 0; for ( int i = 1; i <= VehicleNumber; ++i ) for ( int j = 1; j < Route_Ans[i].V.size(); ++j ) Check_Ans += Graph[Route_Ans[i].V.get(j-1).Number][Route_Ans[i].V.get(j).Number]; System.out.println("Check_Ans="+Check_Ans ); //檢驗是否滿足時間窗約束 boolean flag=true; for (int i=1;i<=VehicleNumber;i++){ UpdateSubT(Route_Ans[i]); if( Route_Ans[i].SubT>0 ) flag=false; } if (flag) System.out.println("Solution satisfies time windows construction"); else System.out.println("Solution not satisfies time windows construction"); }
最後加入一個CheckAns函數，檢驗一下輸出解是否滿足時間窗約束，計算的距離是否正確，有備無患~
TS
private static void addnode(int r,int pos,int Cus) {//節點加入的路徑routes[r],節點customer[Cus],節點加入路徑的位置pos //更新在路徑r中加上節點Cus的需求 routes[r].Load += customers[Cus].Demand; //更新在路徑r中插入節點Cus後所組成的路徑距離 routes[r].Dis = routes[r].Dis - Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos).Number] + Graph[routes[r].V.get(pos-1).Number][customers[Cus].Number] + Graph[routes[r].V.get(pos).Number][customers[Cus].Number]; //在路徑r中插入節點Cus routes[r].V.add (pos ,new CustomerType (customers[Cus]) );//插入i到下標為l處 } private static void removenode(int r,int pos,int Cus) {//節點去除的路徑routes[r],節點customer[cus],節點所在路徑的位置pos //更新在路徑r中去除節點Cus的需求 routes[r].Load -= customers[Cus].Demand; //更新在路徑r中去除節點Cus後所組成的路徑的距離 routes[r].Dis = routes[r].Dis - Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos).Number] - Graph[routes[r].V.get(pos).Number][routes[r].V.get(pos+1).Number] + Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos+1).Number]; //從路徑r中去除節點Cus routes[r].V.remove ( pos ); }
先是兩個輔助函數，addnode和removenode，他們是插入運算元的執行部分。
public static void TabuSearch() { //禁忌搜索 //採取插入運算元，即從一條路徑中選擇一點插入到另一條路徑中 //在該操作下形成的鄰域中選取使目標函數最小化的解 double Temp1; double Temp2; //初始化禁忌表 for ( int i = 2; i <= CustomerNumber + 1; ++i ) { for ( int j = 1; j <= VehicleNumber; ++j ) Tabu[i][j] = 0; TabuCreate[i] = 0; } int Iteration = 0; while ( Iteration < IterMax ) { int BestC = 0; int BestR = 0; int BestP = 0; int P=0; double BestV = INF; for ( int i = 2; i <= CustomerNumber + 1; ++i ) {//對每一個客戶節點 for ( int j = 1; j < routes[customers[i].R].V.size(); ++j ) {//對其所在路徑中的每一個節點 if ( routes[customers[i].R].V.get(j).Number == i ) {//找到節點i在其路徑中所處的位置j P = j;//標記位置 break; } } removenode(customers[i].R,P,i);//將客戶i從原路徑的第P個位置中移除 //找到一條路徑插入刪去的節點 for ( int j = 1; j <= VehicleNumber; ++j ) for ( int l = 1; l < routes[j].V.size(); ++l )//分別枚舉每一個節點所在位置 if ( customers[i].R != j ) { addnode(j,l,i);//將客戶l插入路徑j的第i個位置 Temp1 = routes[customers[i].R].SubT; //記錄原先所在路徑的時間窗違反總和 Temp2 = routes[j].SubT; //記錄插入的路徑時間窗違反總和 //更新i節點移出的路徑： routes[customers[i].R].SubT = 0; UpdateSubT(routes[customers[i].R]); //更新i節點移入的路徑j： routes[j].SubT = 0; UpdateSubT(routes[j]); double TempV = Calculation ( routes, i, j );//計算目標函數值 if((TempV < Ans)|| //藐視準則，如果優於全局最優解 (TempV < BestV && //或者為局部最優解，且未被禁忌 ( routes[j].V.size() > 2 && Tabu[i][j] <= Iteration ) || ( routes[j].V.size() == 2 && TabuCreate[i] <= Iteration ))) //禁忌插入操作，前者為常規禁忌表，禁忌插入運算元；後者為特殊禁忌表，禁忌使用新的車輛 //路徑中節點數超過2，判斷是否禁忌插入運算元；路徑中只有起點、終點，判斷是否禁忌使用新車輛。 if ( TempV < BestV ) { //記錄局部最優情況 BestV = TempV; //best vehicle 所屬車輛 BestC = i; //best customer客戶 BestR = j; //best route 所屬路徑 BestP = l; //best position所在位置 } //節點新路徑復原 routes[customers[i].R].SubT = Temp1; routes[j].SubT = Temp2; removenode(j,l,i); } //節點原路徑復原 addnode(customers[i].R,P,i); } //更新車輛禁忌表 if ( routes[BestR].V.size() == 2 ) TabuCreate[BestC] = Iteration + 2 * TabuTenure + (int)(random() * 10); //更新禁忌表 Tabu[BestC][customers[BestC].R] = Iteration + TabuTenure + (int)(random() * 10); //如果全局最優的節點正好屬於當前路徑，過 for ( int i = 1; i < routes[customers[BestC].R].V.size(); ++i ) if ( routes[customers[BestC].R].V.get(i).Number == BestC ) { P = i; break; } //依據上述循環中挑選的結果，生成新的總體路徑規劃 //依次更新改變過的路徑的：載重，距離長度，超出時間窗重量 //更新原路徑 removenode(customers[BestC].R,P,BestC); //更新新路徑 addnode(BestR,BestP,BestC); //更新超出時間 routes[BestR].SubT = 0; UpdateSubT(routes[BestR]); routes[customers[BestC].R].SubT = 0; UpdateSubT(routes[customers[BestC].R]); //更新被操作的節點所屬路徑編號 customers[BestC].R = BestR; //如果當前解合法且較優則更新存儲結果 if ( ( Check ( routes ) == true ) && ( Ans > BestV ) ) { for ( int i = 1; i <= VehicleNumber; ++i ) { Route_Ans[i].Load = routes[i].Load; Route_Ans[i].V.clear(); for ( int j = 0; j < routes[i].V.size(); ++j ) Route_Ans[i].V.add ( routes[i].V.get(j) ); } Ans = BestV; } Iteration++; } }
終！於！到了最後一步了！
現在萬事俱備，只欠東風，只需要按照禁忌搜索的套路，將所有工具整合在一起，搭建出程式碼框架就ok啦。
由於我們採用routes[ ]數組存儲當前解，因此在進行插入操作之前要存儲部分數據，在計算完目標函數之後要進行復原操作。
在更新禁忌表時，對禁忌步長的計算公式可以靈活改變。
記得對局部最優解進行判斷，再選取為可行的全局最優解。
算例展示
我們採用標準solomon測試數據c101.txt進行測試。（算例可在留言區下載獲取）
VehicleNumber = 25；
Capacity=200；
分別測試節點數25，50，100的情況。精確解分別為：
CustomerNumber=25：
CustomerNumber=50：
CustomerNumber=100：
可見我們的程式碼精確度還是很可以滴~~
當然不排除運氣不好，得出很差解的情況。不相信自己的人品可以手動調整迭代次數IterMax。
本期的內容到這裡就差不多結束了！開心！
在這裡提醒大家一下，在針對啟發式演算法的學習過程中，編寫程式碼的能力是很重要的。VRPTW是一個很好的載體，建議有時間的讀者盡量將學到的演算法知識運用到實踐中去。小編會和你們一起學習進步的！
下次再見ヾ(￣▽￣)Bye~Bye~
在此提前祝螢幕前的你元旦、聖誕快樂，
新的一年裡學業有成，萬事如意！
參考文獻：
Cordeau, J. F. , Laporte, G. , & Mercier, A. . (2001). A unified tabu search heuristic for vehicle routing problems with time windows. Journal of the Operational Research Society, 52(8), 928-936.
程式碼參考：
乾貨 | 十分鐘掌握禁忌搜索演算法求解帶時間窗的車輛路徑問題(附C++程式碼和詳細程式碼注釋)
【程式碼及參考資料見留言區】
贊 賞
長按下方二維碼打賞
感謝您，
支援學生們的原創熱情！
鄭重承諾
打賞是對工作的認可
所有打賞所得
都將作為酬勞支付給辛勤工作的學生
指導老師不取一文
精彩文章推薦
運籌學教學|Benders decomposition （二）應用實例及演算法實現（附源程式碼及詳細的程式碼注釋）
運籌學教學 | 十分鐘快速掌握最大流演算法（附C++程式碼及算例）
運籌學教學|列生成（Column Generation）演算法（附程式碼及詳細注釋）
運籌學教學 | 十分鐘教你求解分配問題(assignment problem)
-The End-
文案 && 編輯&&程式碼：周航
（華中科技大學管理學院本科一年級：[email protected]）
指導老師 / 秦虎 華中科技大學管理學院
掃一掃，獲取數據和模型
還有更多演算法學習課件分享喲
