一、节约算法及禁忌搜索算法简介
1 节约算法简介
1.1节约算法定义
1.2 基本思想
行时通过这一条弧。
1.3 迭代步骤
2禁忌搜索算法理论
2.1局部邻域搜索
局部邻域搜索是基于贪婪准则持续地在当前的邻域中进行搜索，虽然其算法通用，易于实现，且容易理解，但其搜索性能完全依赖于邻域结构和初始解，尤其容易陷入局部极小值而无法保证全局优化。
局部搜索的算法可以描述为：
这种邻域搜索方法易于理解，易于实现，而且具有很好的通用性，但是搜索结果的好坏完全依赖于初始解和邻域的结构。若邻域结构设置不当，或初始解选择不合适，则搜索结果会很差，可能只会搜
索到局部最优解，即算法在搜索过程中容易陷入局部极小值。因此，若不在搜索策略上进行改进，要实现全局优化，局部邻域搜索算法采用的邻域函数就必须是“完全”的，即邻域函数将导致解的完全枚
举。而这在大多数情况下是无法实现的，而且穷举的方法对于大规模问题在搜索时间上也是不允许的。为了实现全局搜索，禁忌搜索采用允许接受劣质解的策略来避免局部最优解。
2.2禁忌搜索
禁忌搜索算法是模拟人的思维的一种智能搜索算法，即人们对已搜索的地方不会再立即去搜索，而是去对其他地方进行搜索，若没有找到，可再搜索已去过的地方。禁忌搜索算法从一个初始可行解出
发，选择一系列的特定搜索方向(或称为“移动”)作为试探，选择使目标函数值减小最多的移动。为了避免陷入局部最优解，禁忌搜索中采用了一种灵活的“记忆”技术，即对已经进行的优化过程进行记录，指导下一步的搜索方向，这就是禁忌表的建立。禁忌表中保存了最近若干次迭代过程中所实现的移动，凡是处于禁忌表中的移动，在当前迭代过程中是禁忌进行的，这样可以避免算法重新访问在最近若
干次迭代过程中已经访问过的解，从而防止了循环，帮助算法摆脱局部最优解。另外，为了尽可能不错过产生最优解的“移动”，禁忌搜索还采用“特赦准则”的策略。
对一个初始解，在一种邻域范围内对其进行一系列变化，从而得到许多候选解。从这些候选解中选出最优候选解，将候选解对应的目标值与“best so far”状态进行比较。若其目标值优于“best sofar”状态， 就将该候选解解禁， 用来替代当前最优解及其“best sofar”状态， 然后将其加入禁忌表， 再将禁忌表中相应对象的禁忌长度改变：如果所有的候选解中所对应的目标值都不存在优于“best sofar”状态， 就从这些候选解中选出不属于禁忌对象的最佳状态， 并将其作为新的当前解，不用与当前最优解进行比较，直接将其所对应的对象作为禁忌对象，并将禁忌表中相应对象的禁忌长度进行修改。
2.3禁忌搜索算法的特点
禁忌搜索算法是在邻域搜索的基础上，通过设置禁忌表来禁忌一些已经进行过的操作，并利用藐视准则来奖励一些优良状态，其中邻域结构、候选解、禁忌长度、禁忌对象、藐视准则、终止准则等是影
响禁忌搜索算法性能的关键。邻域函数沿用局部邻域搜索的思想，用于实现邻域搜索；禁忌表和禁忌对象的设置，体现了算法避免迂回搜索的特点：藐视准则，则是对优良状态的奖励，它是对禁忌策略的一种放松。
与传统的优化算法相比，禁忌搜索算法的主要特点是：
(1)禁忌搜索算法的新解不是在当前解的邻域中随机产生，它要么是优于“best so far”的解， 要么是非禁忌的最佳解， 因此选取优良解的概率远远大于其他劣质解的概率。
(2)由于禁忌搜索算法具有灵活的记忆功能和藐视准则，并且在搜索过程中可以接受劣质解，所以具有较强的“爬山”能力，搜索时能够跳出局部最优解，转向解空间的其他区域，从而增大获得更好的全局最优解的概率。因此，禁忌搜索算法是一种局部搜索能力很强的全局迭代寻优算法。
2.4禁忌搜索算法的改进方向 禁忌搜索是著名的启发式搜索算法，但是禁忌搜索也有明显的不足，即在以下方面需要改进： (1)对初始解有较强的依赖性，好的初始解可使禁忌搜索算法在解空间中搜索到好的解，而较差的初始解则会降低禁忌搜索的收敛速度。因此可以与遗传算法、模拟退火算法等优化算法结合，先产生较好的初始解，再用禁忌搜索算法进行搜索优化。 (2)迭代搜索过程是串行的，仅是单一状态的移动，而非并行搜索。为了进一步改善禁忌搜索的性能，一方面可以对禁忌搜索算法本身的操作和参数选取进行改进，对算法的初始化、参数设置等方面实 施并行策略，得到各种不同类型的并行禁忌搜索算法[9]：另一方面则可以与遗传算法、神经网络算法以及基于问题信息的局部搜索相结合。 (3)在集中性与多样性搜索并重的情况下，多样性不足。集中性搜索策略用于加强对当前搜索的优良解的邻域做进一步更为充分的搜索，以期找到全局最优解。多样性搜索策略则用于拓宽搜索区域，尤 其是未知区域，当搜索陷入局部最优时，多样性搜索可改变搜索方向，跳出局部最优，从而实现全局最优。增加多样性策略的简单处理手段是对算法的重新随机初始化，或者根据频率信息对一些已知对象进行惩罚。
3 禁忌搜索算法流程
简单禁忌搜索算法的基本思想是：给定一个当前解(初始解)和一种邻域，然后在当前解的邻域中确定若干候选解；若最佳候选解对应的目标值优于“best so far”状态， 则忽视其禁忌特性， 用它替代当前解和“best so far”状态， 并将相应的对象加入禁忌表， 同时修改禁忌表中各对象的任期：若不存在上述候选解，则在候选解中选择非禁忌的最佳状态为新的当前解，而无视它与当前解的优劣，同时将相应的对象加入禁忌表，并修改禁忌表中各对象的任期。如此重复上述迭代搜索过程，直至满足停止准则。其算法步骤可描述如下：
(1)给定禁忌搜索算法参数，随机产生初始解x，置禁忌表为空。
(2)判断算法终止条件是否满足：若是，则结束算法并输出优化结果：否则，继续以下步骤。
(3)利用当前解的邻域函数产生其所有(或若干)邻域解，并从中确定若干候选解。
(4)对候选解判断藐视准则是否满足：若满足，则用满足藐视准则的最佳状态y替代x成为新的当前解，即x=y，并用与y对应的禁忌对象替换最早进入禁忌表的禁忌对象， 同时用y替换“best so far”状态，然后转步骤(6)：否则，继续以下步骤。
(5)判断候选解对应的各对象的禁忌属性，选择候选解集中非禁忌对象对应的最佳状态为新的当前解，同时用与之对应的禁忌对象替换最早进入禁忌表的禁忌对象。
(6)判断算法终止条件是否满足：若是，则结束算法并输出优化结果：否则，转步骤(3)。
禁忌搜索算法的运算流程如图8.1所示。
4 关键参数说明
一般而言，要设计一个禁忌搜索算法，需要确定算法的以下环节：初始解、适配值函数、邻域结构、禁忌对象、候选解选择、禁忌表、禁忌长度、藐视准则、搜索策略、终止准则[10，11]。面对如此众
多的参数，针对不同邻域的具体问题，很难有一套比较完善的或非常严格的步骤来确定这些参数。
初始解
禁忌搜索算法可以随机给出初始解，也可以事先使用其他启发式算法等给出一个较好的初始解。由于禁忌搜索算法主要是基于邻域搜索的，初始解的好坏对搜索的性能影响很大。尤其是一些带有很复杂
约束的优化问题，如果随机给出的初始解很差，甚至通过多步搜索也很难找到一个可行解，这时应该针对特定的复杂约束，采用启发式方法或其他方法找出一个可行解作为初始解；再用禁忌搜索算法求解，以提高搜索的质量和效率。也可以采用一定的策略来降低禁忌搜索对初始解的敏感性。
适配值函数
禁忌搜索的适配值函数用于对搜索进行评价，进而结合禁忌准则和特赦准则来选取新的当前状态。目标函数值和它的任何变形都可以作为适配值函数。若目标函数的计算比较困难或耗时较长，此时可采
用反映问题目标的某些特征值来作为适配值，进而改善算法的时间性能。选取何种特征值要视具体问题而定，但必须保证特征值的最佳性与目标函数的最优性一致。适配值函数的选择主要考虑提高算法的效率、便于搜索的进行等因素。
邻域结构
所谓邻域结构，是指从一个解(当前解)通过“移动”产生另一个解(新解)的途径，它是保证搜索产生优良解和影响算法搜索速度的重要因素之一。邻域结构的设计通常与问题相关。邻域结构的设计方法很多，对不同的问题应采用不同的设计方法，常用设计方法包括互换、插值、逆序等。不同的“移动”方式将导致邻域解个数及其变化情况的不同，对搜索质量和效率有一定影响。
通过移动，目标函数值将产生变化，移动前后的目标函数值之差，称之为移动值。如果移动值是非负的，则称此移动为改进移动：否则，称之为非改进移动。最好的移动不一定是改进移动，也可能是 非改进移动，这一点能保证在搜索陷入局部最优时，禁忌搜索算法能自动把它跳出局部最优。
禁忌对象 所谓禁忌对象，就是被置入禁忌表中的那些变化元素。禁忌的目的则是为了尽量避免迂回搜索而多搜索一些解空间中的其他地方。归纳而言，禁忌对象通常可选取状态本身或状态分量等。
候选解选择 候选解通常在当前状态的邻域中择优选取，若选取过多将造成较大的计算量，而选取较少则容易“早熟”收敛，但要做到整个邻域的择优往往需要大量的计算，因此可以确定性地或随机性地在部分邻域中选取候选解，具体数据大小则可视问题特征和对算法的要求而定。 禁忌表 不允许恢复(即被禁止) 的性质称作禁忌(Tabu) 。禁忌表的主要目的是阻止搜索过程中出现循环和避免陷入局部最优，它通常记录前若干次的移动，禁止这些移动在近期内返回。在迭代固定次数后， 禁忌表释放这些移动，重新参加运算，因此它是一个循环表，每迭代一次，就将最近的一次移动放在禁忌表的末端，而它的最早的一个移动就从禁忌表中释放出来。 从数据结构上讲，禁忌表是具有一定长度的先进先出的队列。禁忌搜索算法使用禁忌表禁止搜索曾经访问过的解，从而禁止搜索中的局部循环。禁忌表可以使用两种记忆方式：明晰记忆和属性记忆。明 晰记忆是指禁忌表中的元素是一个完整的解，消耗较多的内存和时间：属性记忆是指禁忌表中的元素记录当前解移动的信息，如当前解移动的方向等。
禁忌长度 所谓禁忌长度，是指禁忌对象在不考虑特赦准则的情况下不允许被选取的最大次数。通俗地讲，禁忌长度可视为禁忌对象在禁忌表中的任期。禁忌对象只有当其任期为0时才能被解禁。在算法的设计和构 造过程中，一般要求计算量和存储量尽量小，这就要求禁忌长度尽量小。但是，禁忌长度过小将造成搜索的循环。禁忌长度的选取与问题特征相关，它在很大程度上决定了算法的计算复杂性。 一方面，禁忌长度可以是一个固定常数(如t=c，c为一常数)，或者固定为与问题规模相关的一个量(如t=√n，n为问题维数或规模)，如此实现起来方便、简单，也很有效：另一方面，禁忌长度也可以是动态变化的，如根据搜索性能和问题特征设定禁忌长度的变化区间，而禁忌长度则可按某种规则或公式在这个区间内变化。
藐视准则 在禁忌搜索算法中，可能会出现候选解全部被禁忌，或者存在一个优于“best so far”状态的禁忌候选解， 此时特赦准则将某些状态解禁，以实现更高效的优化性能。特赦准则的常用方式有： (1) 基于适配值的原则：某个禁忌候选解的适配值优于“bestso far”状态， 则解禁此候选解为当前状态和新的“best so far”状态。 (2)基于搜索方向的准则：若禁忌对象上次被禁忌时使得适配值有所改善，并且目前该禁忌对象对应的候选解的适配值优于当前解，则对该禁忌对象解禁。
搜索策略 搜索策略分为集中性搜索策略和多样性搜索策略。集中性搜索策略用于加强对优良解的邻域的进一步搜索。其简单的处理手段可以是在一定步数的迭代后基于最佳状态重新进行初始化，并对其邻域进行再次搜索。在大多数情况下，重新初始化后的邻域空间与上一次的邻域空间是不一样的，当然也就有一部分邻域空间可能是重叠的。多样性搜索策略则用于拓宽搜索区域，尤其是未知区域。其简单的处理手段可以是对算法的重新随机初始化，或者根据频率信息对一些已知对象进行惩罚。
终止准则 禁忌搜索算法需要一个终止准则来结束算法的搜索进程，而严格理论意义上的收敛条件，即在禁忌长度充分大的条件下实现状态空间的遍历，这显然是不可能实现的。因此，在实际设计算法时通常采用 近似的收敛准则。常用的方法有： (1)给定最大迭代步数。当禁忌搜索算法运行到指定的迭代步数之后，则终止搜索。 (2)设定某个对象的最大禁忌频率。若某个状态、适配值或对换等对象的禁忌频率超过某一阈值，或最佳适配值连续若干步保持不变，则终止算法。 (3)设定适配值的偏离阈值。首先估计问题的下界，一旦算法中最佳适配值与下界的偏离值小于某规定阈值，则终止搜索。
二、部分源代码
function [outcome1,outcome2,outcome3]=cw(Numberoffacilities,assignofpoint,ttimeu,timewindow,distMatrix,quantity,Qofcar,p1)
H=1;
outcome3=zeros(1,Numberoffacilities);%与设施点的伪编号是一一对应的
%%%%%%%找到实际分配的设施点，因为之前的chrom中为1的点，可能%%%%%%%%%%%%%
facilitypop=zeros(1:Numberoffacilities);
facilitypop(assignofpoint)=1;
trueSelectefacilities=find(facilitypop(1,:)==1); %被选择的设施编号
trueselectNumberoffacilities=size(trueSelectefacilities,2);
outcome1=cell(1,trueselectNumberoffacilities);
outcome2=cell(1,trueselectNumberoffacilities);
for i=1:trueselectNumberoffacilities%此时的设施为i
a=trueSelectefacilities(i);%设施的伪编号
pointofsubroute=find(assignofpoint==a);%属于该设施点的所有需求点的伪编号
Numberofpointsofsubroute=size(pointofsubroute,2);%该设施点的需求点数量
judge=zeros(1,Numberofpointsofsubroute); %判断需求点的位置情况
chrom1=zeros(1,Numberofpointsofsubroute);%最后输入的是伪编号
chrom2=sort(pointofsubroute); %需求点伪编号从小到大排序
originalchrom2=chrom2;
%%对关键节点先分配一个车辆给它？？？？？？？？？？？？？？？？？？？？？对不对？
%for j=1:Numberofpointsofsubroute
%if timewindow(chrom2(j),2)==timewindow(chrom2(j),3)
%chrom1(j)=H;
%H=H+1;
%end
%end
%%%%%%%%计算到达时间%%%%%%%%%%%%%%%%%%%%
arrivetime=zeros(1,Numberofpointsofsubroute);%需求点的到达时间,与当前的originalchrom2位置一一对应，与分配给该设施点的需求点一一对应
for j=1:Numberofpointsofsubroute
if ttimeu(chrom2(j)+Numberoffacilities,a)<timewindow(chrom2(j),2)
arrivetime(j)=timewindow(chrom2(j),2);
else
arrivetime(j)=ttimeu(chrom2(j)+Numberoffacilities,a);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%路径节约值列表%%%%%%%%%%%%%%%%%%%%%%
savingnumber=zeros(Numberofpointsofsubroute);%与分配给该设施的需求点的伪编号一一对应
for j=1:Numberofpointsofsubroute
for z=1:Numberofpointsofsubroute
if j~=z
savingnumber(j,z)=distMatrix(chrom2(j)+Numberoffacilities,a)+distMatrix(chrom2(z)+Numberoffacilities,a)-distMatrix(chrom2(j)+Numberoffacilities,chrom2(z)+Numberoffacilities);
end
end
end
[a b]=max(savingnumber);%a输入值，b中对应的数字为行，对应的列数为列
[c d]=max(a);%c输入具体的当前节约值，d为列数
e=b(d);%e为行
%%%在节约值列表中e行d列是当前最大的节约量
E=chrom2(e);%此时的E为节约量最大的需求点伪编号之一
D=chrom2(d);%此时的D为节约量最大的需求点伪编号之1
done=1;
while(done<2)
if judge(e)==0&&judge(d)==0
PD=0;
a=quantity(E,2)+quantity(D,2);
aa=quantity(E,3)+quantity(D,3);
pp1=normcdf(Qofcar,a,sqrt(aa));
a2=find(chrom2==D);
b2=find(chrom2==E);
if pp1>p1
EFj=arrivetime(e)+ttimeu(E+Numberoffacilities,D+Numberoffacilities)-arrivetime(d);
if EFj<0
aheadtime=arrivetime(d)-timewindow(D,2);
if aheadtime>=-EFj
PD=1;
end
elseif EFj==0
PD=1;
elseif EFj>0
delaytime=timewindow(D,3)-arrivetime(d);
if delaytime>=EFj
PD=1;
end
end
if PD==1
chrom1(e)=H;
chrom1(d)=H;%H为车辆编号
judge(e)=1;%与设施直接相连的起点为1
judge(d)=2;%与设施直接相连的终点为2
H=H+1;
arrivetime(d)=arrivetime(d)+EFj;
if b2>a2
chrom2(a2)=E;
chrom2(b2)=D;
end
end
end
savingnumber(e,d)=0;
elseif judge(e)==0&&judge(d)==1;
a2=find(chrom2==D);
b2=find(chrom2==E);
a=find(chrom1==chrom1(a2));
b=sum(quantity(chrom2(a),2))+quantity(E,2);
bb=sum(quantity(chrom2(a),3))+quantity(E,3);
pp1=normcdf(Qofcar,b,sqrt(bb));
c=size(a,2);
if pp1>p1
EFj=arrivetime(e)+ttimeu(E+Numberoffacilities,D+Numberoffacilities)-arrivetime(d);
if EFj<0
a1=10000;
for j=1:c
b1=find(originalchrom2==chrom2(a(j)));
Aheadtime=arrivetime(b1)-timewindow(chrom2(a(j)),2);
if Aheadtime<a1
aheadtime=Aheadtime;
a1=Aheadtime;
end
end
if aheadtime>=-EFj
PD=1;
end
elseif EFj==0
PD=1;
elseif EFj>0
a1=10000;
for j=1:c
b1=find(originalchrom2==chrom2(a(j)));
Delaytime=timewindow(chrom2(a(j)),3)-arrivetime(b1);
if Delaytime<a1
delaytime=Delaytime;
a1=Delaytime;
end
end
if delaytime>=EFj
PD=1;
end
end
if PD==1
chrom1(b2)=chrom1(a2);
a1=find(chrom1==chrom1(a2));
CHROM2=chrom2;
chrom2(a1(1))=CHROM2(b2);
k=2;
while(k<=c+1)
chrom2(a1(k))=CHROM2(a(k-1));
k=k+1;
end
judge(d)=10;%10为内点
judge(e)=1;
end
end
savingnumber(e,d)=0;
elseif judge(e)==0&&judge(d)==2
a2=find(chrom2==D);
b2=find(chrom2==E);
a=find(chrom1==chrom1(a2));
b=sum(quantity(chrom2(a),2))+quantity(E,2);
bb=sum(quantity(chrom2(a),3))+quantity(E,3);
pp1=normcdf(Qofcar,b,sqrt(bb));
c=size(a,2);
if pp1>p1
EFj=arrivetime(d)+ttimeu(E+Numberoffacilities,D+Numberoffacilities)-arrivetime(e);
if EFj<0
aheadtimee=arrivetime(e)-timewindow(E,2);
aheadtimed=arrivetime(d)-timewindow(D,2);
if aheadtimee<aheadtimed
aheadtime=aheadtimee;
else
aheadtime=aheadtimed;
end
if aheadtime>=-EFj
PD=1;
end
elseif EFj==0
PD=1;
elseif EFj>0
delaytimee=timewindow(E,3)-arrivetime(e);
delaytimed=timewindow(D,3)-arrivetime(d);
if delaytimee<delaytimed
delaytime=delaytimee;
else
delaytime=delaytimed;
end
if delaytime>=EFj
PD=1;
end
end
if PD==1
chrom1(b2)=chrom1(a2);
a1=find(chrom1==chrom1(a2));
CHROM2=chrom2;
k=1;
while(k<=c)
chrom2(a1(k))=CHROM2(a(k));
k=k+1;
end
chrom2(a1(k))=E;
judge(d)=10;
judge(e)=2;
end
end
savingnumber(e,d)=0;
elseif judge(e)==1&&judge(d)==0
a2=find(chrom2==D);
b2=find(chrom2==E);
a=find(chrom1==chrom1(b2));
b=sum(quantity(chrom2(a),2))+quantity(D,2);
bb=sum(quantity(chrom2(a),3))+quantity(D,3);
pp1=normcdf(Qofcar,b,sqrt(bb));
c=size(a,2);
if pp1>p1
EFj=arrivetime(e)+ttimeu(E+Numberoffacilities,D+Numberoffacilities)-arrivetime(d);
if EFj<0
a1=10000;
for j=1:c
b1=find(originalchrom2==chrom2(a(j)));
Aheadtime=arrivetime(b1)-timewindow(chrom2(a(j)),2);
if Aheadtime<a1
aheadtime=Aheadtime;
a1=Aheadtime;
end
end
三、运行结果
四、matlab版本及参考文献
1 matlab版本 2014a
2 参考文献 [1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016. [2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
