新的代码加入了原先忘加的藐视准则，将一些冗余代码改为函数调用，同时加入大规模注释，很适合没尝试过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路线类，记录该路线的总承载量，总长度，对时间窗约束的总违反量，以及单条路径上的客户节点序列。
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~
在此提前祝屏幕前的你元旦、圣诞快乐，
新的一年里学业有成，万事如意！
