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   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="" i="" .number="in.nextint()+1;" .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是一個很好的載體，建議有時間的讀者盡量將學到的算法知識運用到實踐中去。小編會和你們一起學習進步的！

參考文獻：

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.