VRPTW簡介

VRPTW問題可描述為：假設(shè)一個配送中心為周圍若干個位于不同地理位置、且對貨物送達時間有不相同要求的客戶點提供配送服務(wù)。其中，配送中心用于運行的車輛都是同一型號的（即擁有相同的容量、速度）；配送中心對車輛出入的時間有限制。我們的任務(wù)是找出使所有車輛行使路徑總和最小的路線。

VRPTW的更多詳細介紹可以參考之前的推文：

干貨|十分鐘快速掌握CPLEX求解VRPTW數(shù)學模型（附JAVA代碼及CPLEX安裝流程）

為了保持文章的獨立型，同時方便后續(xù)講解，這里給出建模實例（參考文獻在文末標注）：

所求的所有車輛路線需滿足以下要求：

在此基礎(chǔ)上求出每輛車輛的總時間最短（由于車輛速度相同，時間最短相當于路程最短）的路線。（允許不使用某些車輛）

Tabu Search簡介

禁忌搜索算法（Tabu Search Algorithm，簡稱TS）起源于對于人類記憶功能的模仿，是一種亞啟發(fā)式算法(meta-heuristics)。它從一個初始可行解(initial feasible solution)出發(fā)，試探一系列的特定搜索方向(移動)，選擇讓特定的目標函數(shù)值提升最多的移動。為了避免陷入局部最優(yōu)解，禁忌搜索對已經(jīng)歷過的搜索過程信息進行記錄，從而指導下一步的搜索方向。

禁忌搜索是人工智能的一種體現(xiàn)，是局部搜索的一種擴展。禁忌搜索是在鄰域搜索(local search)的基礎(chǔ)上，通過設(shè)置禁忌表（tabu list）來禁忌一些曾經(jīng)執(zhí)行過的操作，并利用藐視準則來解禁一些優(yōu)秀的解。

TS+VRPTW

對鄰域搜索類算法而言，采取的搜索算子和評價函數(shù)至關(guān)重要。下面詳細介紹代碼中針對VRPTW的插入算子和評價函數(shù)。

插入算子：

評價函數(shù)：

算法概述

Java代碼詳解

代碼主要分為以下幾個類：

Main，主函數(shù)；

CustomerType，存放客戶節(jié)點的信息；

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;//節(jié)點自身編號    int R;//節(jié)點所屬車輛路徑編號    double X, Y;//節(jié)點橫縱坐標    double Begin, End, Service;//節(jié)點被訪問的最早時間，最晚時間以及服務(wù)時長    double Demand;//節(jié)點的需求容量    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;    }}

客戶類，對圖中的每一個客戶，分別構(gòu)建客戶類，存放自身編號，所屬車輛路線，坐標位置，訪問時間窗，服務(wù)所需時長、需求。

RouteType

import java.util.ArrayList;public class RouteType {    public double Load;//單條路徑承載量    public double SubT;//單條路徑違反各節(jié)點時間窗約束時長總和    public double Dis;//單條路徑總長度    public ArrayList   V=new ArrayList<>();    //單條路徑上顧客節(jié)點序列。在route中，第0個、最后一個都為depot，第k個為第k位。}

路線類，記錄該路線的總承載量，總長度，對時間窗約束的總違反量，以及單條路徑上的客戶節(jié)點序列。

Parameter

public class Parameter {    public static double INF=Double.MAX_VALUE;    public static int CustomerNumber=100;//算例中除倉庫以外的顧客節(jié)點個數(shù)    public static int VehicleNumber = 25;    public static int Capacity=200;//車輛最大容量    public static int IterMax=2000;//搜索迭代次數(shù)    public static int TabuTenure=20;//禁忌步長    public static int[][] Tabu=new int[CustomerNumber + 10][VehicleNumber + 10];//禁忌表用于禁忌節(jié)點插入操作:[i][j]將節(jié)點i插入路徑j(luò)中    public static int[] TabuCreate=new int[CustomerNumber + 10];//禁忌表用于緊急拓展新路徑或使用新車輛    public static double Ans;//最優(yōu)解距離    public static double Alpha = 1, Beta = 1, Sita = 0.5;//Alpha，Beta為系數(shù)，計算目標函數(shù)值；Sita控制系數(shù)改變的速度    public static double[][] Graph=new double[CustomerNumber + 10][CustomerNumber + 10];//記錄圖    public static CustomerType[] customers=new CustomerType[CustomerNumber+10];//存儲客戶數(shù)據(jù)    public static RouteType[] routes=new RouteType[CustomerNumber+10];//存儲當前解路線數(shù)據(jù)    public static RouteType[] Route_Ans=new RouteType[CustomerNumber+10];//存儲最優(yōu)解路線數(shù)據(jù)}

參數(shù)類，有關(guān)VRPTW和TS的變量都存儲在這里，在這里修改數(shù)據(jù)。

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;        //分別根據(jù)約束滿足的情況和控制系數(shù)Sita更新Alpha和Beta值        //新路徑滿足條件，懲罰系數(shù)減小，        //新路徑違反條件，懲罰系數(shù)加大。        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函數(shù)是對產(chǎn)生解的檢驗。

由于插入算子產(chǎn)生的解并不都滿足所有約束條件，對局部搜索產(chǎn)生的較優(yōu)解需要判斷是否滿足時間窗約束和容量約束后，再決定是否為可行解。

在check局部最優(yōu)解的過程中，修改懲罰系數(shù)Alpha、Beta的值。

public static void UpdateSubT(RouteType r) {//更新路徑r對時間窗的違反量        double ArriveTime =0;        for ( int j = 1; j < r.V.size(); ++j ) {//對每一個節(jié)點分別計算超出時間窗的部分            ArriveTime = ArriveTime                     + r.V.get(j-1).Service //服務(wù)時間                    + Graph[r.V.get(j-1).Number][r.V.get(j).Number];//路途經(jīng)過時間            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函數(shù)更新一條車輛路線中在每一個客戶點的時間窗違反量。通過遍歷整條路線累加得到結(jié)果。

//計算路徑規(guī)劃R的目標函數(shù)值，通過該目標函數(shù)判斷解是否較優(yōu)    public static double Calculation ( RouteType R[], int Cus, int NewR ) {        //目標函數(shù)主要由三個部分組成：D路徑總長度（優(yōu)化目標），Q超出容量約束總量，T超出時間窗約束總量        //目標函數(shù)結(jié)構(gòu)為 f(R) = D + Alpha * Q + Beta * T, 第一項為問題最小化目標，后兩項為懲罰部分        //其中Alpha與Beta為變量，分別根據(jù)當前解是否滿足兩個約束進行變化，根據(jù)每輪迭代得到的解在Check函數(shù)中更新                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);//返回目標函數(shù)值    }

Calculate函數(shù)計算目標函數(shù)值，懲罰部分累加后乘以懲罰系數(shù)。

InitAndPrint

//計算圖上各節(jié)點間的距離    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 ) );    }

根據(jù)計算距離。

//讀取數(shù)據(jù)    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;//終點            //算例中給出節(jié)點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類中的參數(shù)），并對當前解routes[ ]的每條路線進行初始化，起終點都為配送中心。

記錄客戶間距離，存儲在Graph數(shù)組中。

//構(gòu)造初始解    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;        //以滿足容量約束為目的的隨機初始化        //即隨機挑選一個節(jié)點插入到第m條路徑中，若超過容量約束，則插入第m+1條路徑        //且插入路徑的位置由該路徑上已存在的各節(jié)點的最早時間決定        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--;//將當前服務(wù)過的節(jié)點賦值為最末節(jié)點值,數(shù)組容量減1            //隨機提取出一個節(jié)點，類似產(chǎn)生亂序隨機序列的代碼            if ( routes[Current_Route].Load + customers[C].Demand > Capacity )                Current_Route++;            //不滿足容量約束，下一條車輛路線            for ( int i = 0; i < routes[Current_Route].V.size() - 1; i++ )//對路徑中每一對節(jié)點查找，看是否能插入新節(jié)點                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]) );                    //判斷時間窗開始部分，滿足，則加入該節(jié)點。                    routes[Current_Route].Load += customers[C].Demand;                    customers[C].R = Current_Route;                    //更新路徑容量，節(jié)點類。                    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構(gòu)造初始解。根據(jù)前文偽代碼構(gòu)造初始解，每次隨機選擇節(jié)點（類似打亂有序數(shù)列）。

針對該節(jié)點找到符合容量約束，同時時間窗開啟時間符合要求的位置，插入該節(jié)點。記得在插入節(jié)點時同時更新該節(jié)點所屬的路徑。

對時間窗違反量進行初始化。

public static void Output () {//結(jié)果輸出        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("************************************************************");    }

輸出結(jié)果，不再贅述。

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函數(shù)，檢驗一下輸出解是否滿足時間窗約束，計算的距離是否正確，有備無患~

TS

private static void addnode(int r,int pos,int Cus) {//節(jié)點加入的路徑routes[r],節(jié)點customer[Cus],節(jié)點加入路徑的位置pos        //更新在路徑r中加上節(jié)點Cus的需求        routes[r].Load += customers[Cus].Demand;        //更新在路徑r中插入節(jié)點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中插入節(jié)點Cus        routes[r].V.add (pos ,new CustomerType (customers[Cus]) );//插入i到下標為l處    }    private static void removenode(int r,int pos,int Cus) {//節(jié)點去除的路徑routes[r],節(jié)點customer[cus],節(jié)點所在路徑的位置pos         //更新在路徑r中去除節(jié)點Cus的需求        routes[r].Load -= customers[Cus].Demand;        //更新在路徑r中去除節(jié)點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中去除節(jié)點Cus        routes[r].V.remove ( pos );    }

先是兩個輔助函數(shù)，addnode和removenode，他們是插入算子的執(zhí)行部分。

public static void TabuSearch() {        //禁忌搜索        //采取插入算子，即從一條路徑中選擇一點插入到另一條路徑中        //在該操作下形成的鄰域中選取使目標函數(shù)最小化的解        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 ) {//對每一個客戶節(jié)點                for ( int j = 1; j < routes[customers[i].R].V.size(); ++j ) {//對其所在路徑中的每一個節(jié)點                    if ( routes[customers[i].R].V.get(j).Number == i ) {//找到節(jié)點i在其路徑中所處的位置j                        P = j;//標記位置                        break;                    }                }                removenode(customers[i].R,P,i);//將客戶i從原路徑的第P個位置中移除                //找到一條路徑插入刪去的節(jié)點                for ( int j = 1; j <= VehicleNumber; ++j )                      for ( int l = 1; l < routes[j].V.size(); ++l )//分別枚舉每一個節(jié)點所在位置                            if ( customers[i].R != j ) {                                addnode(j,l,i);//將客戶l插入路徑j(luò)的第i個位置                                Temp1 = routes[customers[i].R].SubT;  //記錄原先所在路徑的時間窗違反總和                                Temp2 = routes[j].SubT;               //記錄插入的路徑時間窗違反總和                                //更新i節(jié)點移出的路徑：                                routes[customers[i].R].SubT = 0;                                UpdateSubT(routes[customers[i].R]);                                //更新i節(jié)點移入的路徑j(luò)：                                routes[j].SubT = 0;                                UpdateSubT(routes[j]);                                double TempV = Calculation ( routes, i, j );//計算目標函數(shù)值                                if((TempV < Ans)|| //藐視準則，如果優(yōu)于全局最優(yōu)解                                        (TempV < BestV &&   //或者為局部最優(yōu)解，且未被禁忌                                           ( routes[j].V.size() > 2 && Tabu[i][j] <= Iteration ) || ( routes[j].V.size() == 2 && TabuCreate[i] <= Iteration )))                                    //禁忌插入操作，前者為常規(guī)禁忌表，禁忌插入算子；后者為特殊禁忌表，禁忌使用新的車輛                                    //路徑中節(jié)點數(shù)超過2，判斷是否禁忌插入算子；路徑中只有起點、終點，判斷是否禁忌使用新車輛。                                if ( TempV < BestV ) { //記錄局部最優(yōu)情況                                    BestV = TempV; //best vehicle 所屬車輛                                    BestC = i;     //best customer客戶                                    BestR = j;     //best route   所屬路徑                                    BestP = l;     //best position所在位置                                }                                //節(jié)點新路徑復原                                routes[customers[i].R].SubT = Temp1;                                routes[j].SubT = Temp2;                                removenode(j,l,i);                            }                //節(jié)點原路徑復原                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);            //如果全局最優(yōu)的節(jié)點正好屬于當前路徑，過            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;                }            //依據(jù)上述循環(huán)中挑選的結(jié)果，生成新的總體路徑規(guī)劃            //依次更新改變過的路徑的：載重，距離長度，超出時間窗重量            //更新原路徑            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]);            //更新被操作的節(jié)點所屬路徑編號            customers[BestC].R = BestR;            //如果當前解合法且較優(yōu)則更新存儲結(jié)果            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++;        }    }

終！于！到了最后一步了！

現(xiàn)在萬事俱備，只欠東風，只需要按照禁忌搜索的套路，將所有工具整合在一起，搭建出代碼框架就ok啦。

由于我們采用routes[ ]數(shù)組存儲當前解，因此在進行插入操作之前要存儲部分數(shù)據(jù)，在計算完目標函數(shù)之后要進行復原操作。

在更新禁忌表時，對禁忌步長的計算公式可以靈活改變。

記得對局部最優(yōu)解進行判斷，再選取為可行的全局最優(yōu)解。

算例展示

我們采用標準solomon測試數(shù)據(jù)c101.txt進行測試。（算例可在留言區(qū)下載獲取）

VehicleNumber = 25；

Capacity=200；

分別測試節(jié)點數(shù)25，50，100的情況。精確解分別為：

CustomerNumber=25：

CustomerNumber=50：

CustomerNumber=100：

可見我們的代碼精確度還是很可以滴~~

當然不排除運氣不好，得出很差解的情況。不相信自己的人品可以手動調(diào)整迭代次數(shù)IterMax。

本期的內(nèi)容到這里就差不多結(jié)束了！開心！

在這里提醒大家一下，在針對啟發(fā)式算法的學習過程中，編寫代碼的能力是很重要的。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.