1.TS 介绍

TS是local search的拓展，是全局逐步寻优的全局性邻域搜索算法。TS模仿人类的记忆功能，在搜索过程红标记已经找到的局部最优解及求解过程，并于之后的搜索中避开它们。

算法通过禁忌策略实现记忆功能，通过破禁准则继承LS的强局部搜索能力，防止算法在优化中陷入局部最优。

2.Tabu Search 主要构成要素

（1）评价函数（Evaluation Function）:评价函数是用来评价淋雨中的邻居、判断其优劣的衡量指标。大多数情况下，评价函数为目标函数。

（2）邻域移动（Move Operatior）：邻域移动是进行解转移的关键，又称“算子”，影响整个算法的搜索速度。

（3）禁忌表（Tabu Table）：禁忌表记录被禁止的变化，以防出现搜索循环、陷入局部最优。设计的最关键因素是禁忌对象（禁忌表限定的对象）和禁忌步长（对象的禁忌在多少次迭代后失效）

（4）邻居选择策略（Neighbor selection strategy）:选择最佳邻域移动的规则。目前最广泛采用的是“最好解有限策略”和“第一个改进解优先策略”。前者需比较所有邻域，耗时较久，但解的收敛更有效；后者在发现第一个改进解救进行转移，耗时少，但收敛效率若于前者，对于邻域解空间较大的问题往往比较适合。

（5）破禁准则（Aspiration Criterion）：破禁准则是对于禁忌表的适度放松。

（6）停止规则（Stop Criterion）：停止规则如最大迭代数，算法运行时间，给定数目的迭代内不能改进解惑组合策略等

3.禁忌搜索解决VRPTW

3.1 禁忌搜索算法基本步骤：

1.初始化：利用贪婪算法等局部搜索算法生成一个初始解，清空禁忌表，设置禁忌长度

2.邻域搜索产生候选解，步骤1产生初始解通过搜索算子产生候选解，并计算各个候选解的适应值。

3.选择最好的候选解,从步骤2产生的所有候选解中选出适应值最好的候选解，将其于当前最好的解（即搜索算法开始到现在找到的最好的解）进行比较，如果优于当前最好的解，不考虑是否被禁忌，用这个最好的解来更新当前最好的解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表，

4.判断终止条件一般终止条件为是否到达一定的迭代次数。

TSP问题为例

下图，有5个城市，任何两个城市之间的距离都是确定的

1.编码和解码，构造初始解：

初始解对应的适应值为：F=5+7+8+12+9=41

2.搜索算子

（1）relocation 算子

（2）swap算子

**3.邻域：**从当前解对一系列的搜索方向进行一次试探（通过算子搜索一次）

3.2车辆路径问题c++代码（VRPTW）

采用插入算子，总迭代次数2000，禁忌长度40

输入的形式为表格表示

我们设置三个Alpha Beta Sita 1 1 0.5

目标函数设置为f(x)=D + AlphaQ+BetaT 此时Alpha 和Beta为可变参数，分别根据当前解是否满足两个约束来进行变化（在Check函数中更新，由于Check针对每轮迭代得到的解）

代码创造两个结构体

struct 1 : 记录各个消费者的坐标 编号 服务时间 时间窗等一些信息

struct 2： 记录的是一个路径信息，包含总长度及相关节点和超时总时长

对于到达时间有一个常识的公式

本节点的arrive time = 上一节点arrive time + serve time+路程时间

代码重点注意

（1）算子的插入与重组过程中，对应的路径记录数组，如何更新，如erase函数的运用来去除路径的中间节点：

v = [0,1,2,3,4];//想要删除第3个元素，对应的索引为2p=2;//删除第三个P=2v.erase(v.begin()+p);

（2）禁忌表的制定

对于禁忌表 tabu[customer][route]记录对应的禁忌

#include<iostream>#include<cstdio>#include<cstdlib>//伪随机数生成#include<fstream>#include<cmath>#include <ctime>#include <vector>using namespace std;#define cin fin // 重新定义了关键词，读写到文件中#define cout fout#define INF 0x3ffffff#define Customer_Number 100 //算例中除仓库以外的顾客节点个数#define Capacity 200 //车辆的容量#define Iter_Max 2000 //搜索最大迭代次数#define Tabu_tenure 20 //禁忌时长ifstream fin("R101.txt");ofstream fout("R101_Output.txt");struct Customer_Type {int Number;//节点自身编号int R;//节点所属车辆路径编号double X, Y;//节点横纵坐标double Begin, End, Service;//节点被访问的最早时间，最晚时间，服务时长double Demand;//节点需求量} Customer[Customer_Number + 10];//仓库节点编号为1，顾客节点编号为2-101struct Route_Type {double Load;//单条路径装载量double SubT;//单条路径违反各节点时间窗约束时长总和double Dis;//单条路径总长度vector<Customer_Type>V;//单条路径上顾客节点序列}Route[Customer_Number + 10], Route_Ans[Customer_Number + 10];//车辆路径及搜索到的最优的车辆路径int Vehicle_Number = Customer_Number;//由于无车辆数量限制，因此将上限设为顾客总数int Tabu[Customer_Number + 10][Customer_Number + 10];//禁忌表用于禁忌节点插入插入操作int TabuCreate[Customer_Number + 10];//禁忌表用于禁忌拓展的新路径或使用新车辆double Ans;double Alpha = 1, Beta = 1, Sita = 0.5;double Graph[Customer_Number + 10][Customer_Number + 10];//记录顾客与顾客之间 的距离//************************************************************double Distance(Customer_Type C1, Customer_Type C2) {//计算图上各节点间的距离return sqrt((C1.X - C2.X) * (C1.X - C2.X) + (C1.Y - C2.Y) * (C1.Y - C2.Y));}//************************************************************double Calculation(Route_Type R[], int Cus, int NewR) {//计算路径规划R的目标函数值//目标函数主要由三个部分组成：D路径总长度（优化目标），Q超出容量约束总量，T超出时间窗约束总量//目标函数结构为 f(R) = D + Alpha * Q + Beta * T, 第一项为问题最小化目标，后两项为惩罚部分//其中Alpha与Beta为可变参数，分别根据当前解是否满足两个约束来进行变化（在Check函数中更新，由于Check针对每轮迭代得到的解）double Q = 0;double T = 0;double D = 0;for (int i = 1; i <= Vehicle_Number; ++i)if (R[i].V.size() > 2 && R[i].Load > Capacity)Q = Q + R[i].Load > Capacity;//计算单条路径上各个节点超出时间创约束的总量（仅更新进行一处和插入节点操作的两条路径）double ArriveTime = 0;R[Customer[Cus].R].SubT = 0;for (int j = 1; j < R[Customer[Cus].R].V.size(); ++j) {ArriveTime = ArriveTime + R[Customer[Cus].R].V[j - 1].Service + Graph[R[Customer[Cus].R].V[j - 1].Number][R[Customer[Cus].R].V[j].Number];if (ArriveTime > R[Customer[Cus].R].V[j].End)R[Customer[Cus].R].SubT = R[Customer[Cus].R].SubT + ArriveTime - R[Customer[Cus].R].V[j].End;else if (ArriveTime > R[Customer[Cus].R].V[j].Begin)ArriveTime = R[Customer[Cus].R].V[j].Begin;}ArriveTime = 0;R[NewR].SubT = 0;for (int j = 1; j < R[NewR].V.size(); ++j) {ArriveTime = ArriveTime + R[NewR].V[j - 1].Service + Graph[R[NewR].V[j - 1].Number][R[NewR].V[j].Number];if (ArriveTime > R[NewR].V[j].End)R[NewR].SubT = R[NewR].SubT + ArriveTime - R[NewR].V[j].End;else if (ArriveTime < R[NewR].V[j].Begin)ArriveTime = R[NewR].V[j].Begin;}for (int i = 1; i < Vehicle_Number; ++i)T += R[i].SubT;//计算路径总长度for (int i = 1; i <= Vehicle_Number; ++i)D += R[i].Dis;return D + Alpha * Q + Beta * T;}//************************************************************bool Check(Route_Type R[]) {//检验路径规划R是否满足所有约束double Q = 0;double T = 0;double D = 0;//检查是否满足容量约束for(int i = 1;i<= Vehicle_Number;++i)if(R[i].V.size() > 2 && R[i].Load > Capacity)Q = Q + R[i].Load - Capacity;//检查是否满足时间窗约束for (int i = 1; i <= Vehicle_Number; ++i)T += R[i].SubT;//分别根据约束满足的情况更新Alpha和Beta值if (Q == 0 && Alpha >= 0.001)//满足容量约束，缩小alphaAlpha /= (1 + Sita);else if (Q != 0 && Alpha <= 2000)//不满足容量约束，增加系数值，寻找满足容量约束的解Alpha *= (1 + Sita);if (T == 0 && Beta >= 0.001)//beta 系数和alpha系数思路相同Beta /= (1 + Sita);else if (T != 0 && Beta <= 2000)Beta *= (1 + Sita);if (T == 0 && Q == 0)return true;elsereturn false;}//************************************************************void Copy_Route() {//将路径规划Route的内容复制给路径规划Route_Ansfor (int i = 1; i<= Vehicle_Number; ++i) {Route_Ans[i].Load = Route[i].Load;Route_Ans[i].V.clear();for (int j = 0; j < Route[i].V.size(); ++j)Route_Ans[i].V.push_back(Route[i].V[j]);}}//************************************************************void Output(Route_Type R[]) {//结果输出cout << "************************************************************" << endl;cout << "The Minimum Total Distance = " << Ans << endl;cout << "Concrete Schedule of Each Route as Following : " << endl;int M = 0;for(int i = 1;i<= Vehicle_Number;++i)if (R[i].V.size() > 2) {M++;cout << "No." << M << " : ";for (int j = 0; j < R[i].V.size() - 1; ++j)cout << R[i].V[j].Number << " -> ";cout << R[i].V[R[i].V.size() - 1].Number << endl;}//检验距离计算是否正确double Check_Ans = 0;for (int i = 1; i <= Vehicle_Number; ++i)for (int j = 1; j < R[i].V.size(); ++j)Check_Ans += Graph[R[i].V[j - 1].Number][R[i].V[j].Number];cout << "Check_Ans = " << Check_Ans << endl;cout << "************************************************************" << endl;}//************************************************************void Readin_and_Initialization() {//数据读入及初始化for (int i = 1; i <= Customer_Number + 1; ++i)cin >> Customer[i].Number >> Customer[i].X >> Customer[i].Demand >> Customer[i].Begin >> Customer[i].End >> Customer[i].Service;//初始化每条路径，默认路径收尾为仓库，且首仓库最早最晚时间均为原仓库最早时间，尾仓库则均为原仓库最晚时间Customer[1].R = -1;for (int i = 1; i <= Vehicle_Number; ++i) {if(!Route[i].V.empty())Route[i].V.clear();Route[i].V.push_back(Customer[1]);Route[i].V.push_back(Customer[1]);Route[i].V[0].End = Route[i].V[0].Begin;Route[i].V[1].Begin = Route[i].V[1].End;Route[i].Load = 0;}Ans = INF;for(int i = 1;i<= Customer_Number +1;++i)for(int j = 1;j <= Customer_Number + 1;++j)Graph[i][j] = Distance(Customer[i], Customer[j]);}//************************************************************void Construction() {// 构造初始路径int Customer_Set[Customer_Number + 10];for(int i = 1;i <= Customer_Number;++i)Customer_Set[i] = i + 1;int Sizeof_Customer_Set = Customer_Number;int Current_Route = 1;//以满足容量约束为目的的随机初始化//即随机挑选一个节点插入到第m条路径中，若超过容量约束，则插入第m+1条路径//且插入路径的位置由该路径上已存在的各节点最早时间的升序决定while (Sizeof_Customer_Set > 0) {int K = rand() % Sizeof_Customer_Set + 1;//产生一个1-customer number的随机数int C = Customer_Set[K];Customer_Set[K] = Customer_Set[Sizeof_Customer_Set];Sizeof_Customer_Set--;//将当前服务过的节点赋值为最末节点值,数组容量减1if (Route[Current_Route].Load + Customer[C].Demand > Capacity)Current_Route++;//即随机挑选一个节点插入到第m条路径中，若超过容量约束，则插入第m+1条路径for (int i = 0; i < Route[Current_Route].V.size() - 1; i++)if ((Route[Current_Route].V[i].Begin <= Customer[C].Begin) && (Customer[C].Begin <= Route[Current_Route].V[i + 1].Begin)) {Route[Current_Route].V.insert(Route[Current_Route].V.begin() + i + 1, Customer[C]);Route[Current_Route].Load += Customer[C].Demand;Customer[C].R = Current_Route;break;}}//初始化计算超过容量约束的总量和超过时间窗约束的总量for (int i = 1; i <= Vehicle_Number; ++i) {double ArriveTime = Route[i].V[0].Begin;Route[i].SubT = 0;Route[i].Dis = 0;for (int j = 1; j < Route[i].V.size(); ++j) {ArriveTime = ArriveTime + Route[i].V[j - 1].Service + Graph[Route[i].V[j - 1].Number][Route[i].V[j].Number];Route[i].Dis += Graph[Route[i].V[j - 1].Number][Route[i].V[j].Number];if (ArriveTime > Route[i].V[j].End)Route[i].SubT = Route[i].SubT + ArriveTime - Route[i].V[j].End;else if (ArriveTime < Route[i].V[j].Begin)ArriveTime = Route[i].V[j].Begin;}}}//************************************************************void Tabu_Search() {//禁忌搜索//禁忌搜索采取插入算子，即从一条路径中选择一点插入到另一条路径中//在该操作下形成的邻域中选取使目标函数最小的非禁忌解或者因满足藐视法则而被解禁的解double Temp1;double Temp2;for (int i = 2; i <= Customer_Number + 1; ++i) {for (int j = 1; j <= Vehicle_Number; ++j)Tabu[i][j] = 0;TabuCreate[i] = 0;}int Iteration = 0;//当前迭代次数while (Iteration < Iter_Max) {Iteration++;int BestC = 0, BestR = 0, BestP = 0, P;double BestV = INF;for (int i = 2; i <= Customer_Number + 1; ++i) {for (int j = 1; j < Route[Customer[i].R].V.size(); ++j)if (Route[Customer[i].R].V[j].Number == i) {//找到自身编号P = j;//记录该customer属于该路径第几个节点break;}//从节点原路径中去除该节点的需求Route[Customer[i].R].Load -= Customer[i].Demand;//从节点原路径中去除该节点所组成的路径并重组Route[Customer[i].R].Dis = Route[Customer[i].R].Dis - Graph[Route[Customer[i].R].V[P - 1].Number][Route[Customer[i].R].V[P].Number]- Graph[Route[Customer[i].R].V[P].Number][Route[Customer[i].R].V[P + 1].Number] + Graph[Route[Customer[i].R].V[P - 1].Number][Route[Customer[i].R].V[P + 1].Number];//从节点原路径中去除节点Route[Customer[i].R].V.erase(Route[Customer[i].R].V.begin() + P);for (int j = 1; j <= Vehicle_Number; ++j)//禁忌插入操作，后者为禁止使用新的车辆if ((Route[j].V.size() > 2 && Tabu[i][j] <= Iteration) || (Route[j].V.size() == 2 && TabuCreate[i] <= Iteration)) {for (int l = 1; l < Route[j].V.size(); ++l)if (Customer[i].R != j) {//在节点新路径中加上该节点的需求Route[j].Load += Customer[i].Demand;//在节点新路径中加上该节点插入后所组成的路径并断开原路径Route[j].Dis = Route[j].Dis - Graph[Route[j].V[l - 1].Number][Route[j].V[l].Number]+ Graph[Route[j].V[l - 1].Number][Customer[i].Number] + Graph[Route[j].V[l].Number][Customer[i].Number];//在节点新路径中插入节点Route[j].V.insert(Route[j].V.begin() + l, Customer[i]);Temp1 = Route[Customer[i].R].SubT;Temp2 = Route[j].SubT;double TempV = Calculation(Route, i, j);if (TempV < BestV) {BestV = TempV;BestC = i, BestR = j, BestP = l;}//节点新路径复原Route[Customer[i].R].SubT = Temp1;Route[j].SubT = Temp2;Route[j].V.erase(Route[j].V.begin() + l);Route[j].Load -= Customer[i].Demand;Route[j].Dis = Route[j].Dis + Graph[Route[j].V[l - 1].Number][Route[j].V[l].Number]- Graph[Route[j].V[l - 1].Number][Customer[i].Number] - Graph[Route[j].V[l].Number][Customer[i].Number];}}//节点原路径复原Route[Customer[i].R].V.insert(Route[Customer[i].R].V.begin() + P, Customer[i]);Route[Customer[i].R].Load += Customer[i].Demand;Route[Customer[i].R].Dis = Route[Customer[i].R].Dis + Graph[Route[Customer[i].R].V[P - 1].Number][Route[Customer[i].R].V[P].Number]+ Graph[Route[Customer[i].R].V[P].Number][Route[Customer[i].R].V[P + 1].Number] - Graph[Route[Customer[i].R].V[P - 1].Number][Route[Customer[i].R].V[P + 1].Number];}if (Route[BestR].V.size() == 2)TabuCreate[BestC] = Iteration + 2 * Tabu_tenure + rand() % 10;Tabu[BestC][Customer[BestC].R] = Iteration + Tabu_tenure + rand() % 10;for (int i = 1; i < Route[Customer[BestC].R].V.size(); ++i)if (Route[Customer[BestC].R].V[i].Number == BestC) {P = i;break;}//依据上述循环中挑选的结果，生成新的总体路径规划//更新改变过的各单条路径的载重，距离长度，超出时间窗的总量Route[Customer[BestC].R].Load -= Customer[BestC].Demand;Route[Customer[BestC].R].Dis = Route[Customer[BestC].R].Dis - Graph[Route[Customer[BestC].R].V[P - 1].Number][Route[Customer[BestC].R].V[P].Number]- Graph[Route[Customer[BestC].R].V[P].Number][Route[Customer[BestC].R].V[P + 1].Number] + Graph[Route[Customer[BestC].R].V[P - 1].Number][Route[Customer[BestC].R].V[P + 1].Number];Route[Customer[BestC].R].V.erase(Route[Customer[BestC].R].V.begin() + P);Route[BestR].Dis = Route[BestR].Dis - Graph[Route[BestR].V[BestP - 1].Number][Route[BestR].V[BestP].Number]+ Graph[Route[BestR].V[BestP - 1].Number][Customer[BestC].Number] + Graph[Route[BestR].V[BestP].Number][Customer[BestC].Number];Route[BestR].Load += Customer[BestC].Demand;Route[BestR].V.insert(Route[BestR].V.begin() + BestP, Customer[BestC]);double ArriveTime = 0;Route[BestR].SubT = 0;for (int j = 1; j < Route[BestR].V.size(); ++j) {ArriveTime = ArriveTime + Route[BestR].V[j - 1].Service + Graph[Route[BestR].V[j - 1].Number][Route[BestR].V[j].Number];if (ArriveTime > Route[BestR].V[j].End)Route[BestR].SubT = Route[BestR].SubT + ArriveTime - Route[BestR].V[j].End;else if (ArriveTime < Route[BestR].V[j].Begin)ArriveTime = Route[BestR].V[j].Begin;}ArriveTime = 0;Route[Customer[BestC].R].SubT = 0;for (int j = 1; j < Route[Customer[BestC].R].V.size(); ++j) {ArriveTime = ArriveTime + Route[Customer[BestC].R].V[j - 1].Service + Graph[Route[Customer[BestC].R].V[j - 1].Number][Route[Customer[BestC].R].V[j].Number];if (ArriveTime > Route[Customer[BestC].R].V[j].End)Route[Customer[BestC].R].SubT = Route[Customer[BestC].R].SubT + ArriveTime - Route[Customer[BestC].R].V[j].End;else if (ArriveTime < Route[Customer[BestC].R].V[j].Begin)ArriveTime = Route[Customer[BestC].R].V[j].Begin;}//更新被操作的节点所属路径编号Customer[BestC].R = BestR;//如果当前解合法且较优则更新存储结果if ((Check(Route) == true) && (Ans > BestV)) {Copy_Route();Ans = BestV;}}}//************************************************************int main() {clock_t Start, Finish;//Start = clock();srand((unsigned)time(NULL));Readin_and_Initialization();Construction();Tabu_Search();Output(Route_Ans);//Finish = clock();//cout << "Total Running Time = " << ( Finish - Start ) / 1000.0 << endl;return 0;}//************************************************************