第5章〓禁忌搜索算法
禁忌搜索(tabu search，TS)的思想最早由Glover于1986年提出，它是局部邻域搜索的扩展，是一种全局逐步寻优算法，也是对人类智力过程的一种模拟。禁忌搜索算法通过引入一个灵活的存储结构和相应的禁忌准则来避免迂回搜索，并通过藐视准则来赦免一些被禁忌的优良状态，进而保证多样化的有效探索以最终实现全局优化。相较于模拟退火和遗传算法，禁忌搜索是另一种搜索特点不同的元启发式算法。

本章首先介绍了禁忌搜索算法思想及特点，然后阐述了禁忌搜索算法的设计原则，最后通过几个案例分析了禁忌搜索算法在实际问题中的应用。

5.1禁忌搜索算法思想及特点
5.1.1算法思想

禁忌搜索算法是局部邻域搜索算法的扩展，是人工智能在组合优化算法中的成功应用。禁忌搜索算法可以看作比模拟退火法更为一般的一种邻域搜索算法，它采用了类似爬山法的移动原理，将最近若干步内所得到的解存储在禁忌表(tabu list，TL)中，从而强制搜索避免再次重复表中的解。如果说遗传算法开创了在解空间中从多出发点搜索问题最优解的先河，禁忌搜索法则是首次在搜索过程中使用了记忆功能，它们在求解各种实际应用问题中都取得了相当大的成就。

邻域是给定点附近其他点的集合。在距离空间中，邻域一般定义为以给定点为圆心的一个圆。而在组合优化问题中，邻域一般定义为由给定转化规则对给定的问题域上每个节点进行转化所得到的问题域上节点的集合。

N:x∈D→N(x)∈2D


其中，x∈N(x)称为一个邻域映射，2D表示D的所有子集组成的集合。N(x)称为x的邻域，y∈N(x)称为x的一个邻居。

局部搜索算法基本步骤如下： 

步骤1： 选定一个初始可行解，记录当前最优解xbest←x0，T=N(xbest)。

步骤2： 当T＼{xbest}=Φ或满足其他停止运算准则时，输出计算结果，停止运算； 否则，从T中选一集合S，得到S中的最好解xnow。若f(xnow)<f(xbest)，则xbest←xnow，T=N(xbest)； 否则，T←T＼S； 重复步骤2。

在局部搜索算法中，步骤1的初始可行解可用随机的方法选择，也可以用一些经验方法或其他算法计算得到。步骤2中集合S的选取可以大到N(xbest)本身，也可以小到只有一个元素，例如在N(xbest)中随机选择一点。可以看出，S选取得小，每一步的计算量就少，但是可比较的范围就小； S选取得大，每一步计算的时间增加，比较的范围自然增加。两种情况的应用效果依赖于实际问题。S选取大到N（xbest）本身，称为最优改进(best improvement)，S选取只有一个元素，称为一步改进(first improvement)。步骤2中，其他停止准则是除步骤2中T＼{xbest}=Φ外的准则，这些准则往往取决于对算法的计算时间、计算结果的要求。






局部搜索算法的计算结果主要依赖起点的选取和邻域的结构。同一个起点，不同的邻域结果会得到不同的计算结果。同样，同一个邻域结构，不同的初始点会得到不同的计算结果。因此，在使用局部搜索算法时，为了得到好的解，可以比较不同的邻域结构和不同的初始点。一个非常直观的结论： 如果初始点的选择足够多，总可以计算出全局最优解。



禁忌搜索算法的基本思想是： 给定一个当前解(初始解)和一种邻域，然后在当前解的邻域中确定若干候选解； 若最佳候选解对应的目标值优于“best so far”状态则忽视其禁忌特性，用其替代当前解和“best so far”状态，并将相应的对象加入禁忌表，同时修改禁忌表中各对象的任期； 若不存在上述候选解，则选择在候选解中选择非禁忌的最佳状态为新的当前解，而无视它与当前解的优劣，同时将相应的对象加入禁忌表，并修改禁忌表中各对象的任期； 如此重复上述迭代搜索过程，直至满足停止准则。

其基本步骤可大致叙述如下： 


步骤1： 给定算法参数，随机产生初始解x，置禁忌表为空。

步骤2： 判断算法终止条件是否满足。若是，则结束算法并输出优化结果； 否则，继续以下步骤。

步骤3： 利用当前解x的邻域函数产生其所有(或若干)邻域解，并从中确定若干候选解。

步骤4： 判断候选解是否满足藐视准则。
若成立，则用满足藐视准则的最佳状态y替代x成为新的当前解，即x=y，并用与y对应的禁忌对象替换最早进入禁忌表的禁忌对象，同时用y替换“best so far”状态，然后转步骤6； 否则，继续以下步骤。

步骤5： 判断候选解对应的各对象的禁忌属性，选择候选解集中非禁忌对象对应的最佳状态为新的当前解，同时用与之对应的禁忌对象替换最早进入禁忌表的禁忌对象元素。

步骤6： 转步骤2。

5.1.2算法特点

禁忌搜索算法采用了禁忌技术，所谓禁忌就是禁止重复前面的工作。其重要思想是标记已得到的局部最优解或求解的过程，并在进一步的迭代中避开这些局部最优解或过程。因此不仅可以避免一些重复无效的迭代，节省算法所需时间，还有利于避免陷入局部最优的困境，从而寻求到全局最优解。为了回避局部邻域搜索陷入局部最优的主要不足，禁忌搜索算法用一个禁忌表记录下已经到达过的局部最优点或达到局部最优的一些过程，在下一次搜索中，利用禁忌表中的信息不再或有选择地搜索这些点或过程，以此来跳出局部最优点。

与传统的优化算法相比，禁忌搜索算法的主要特点如下： 

(1) 在搜索过程中可以接受劣解，因此具有较强的“爬山”能力。

(2) 新解不是在当前解的邻域中随机产生，而或是优于“best so far”的解，或是非禁忌的最佳解，因此选取优良解的概率远远大于其他解。

由于禁忌搜索算法具有灵活的记忆功能和藐视准则，并且在搜索过程中可以接受劣解，所以具有较强的“爬山”能力，搜索时能够跳出局部最优解，转向解空间的其他区域，从而增强获得更好的全局最优解的概率，所以禁忌搜索算法是一种局部搜索能力很强的全局迭代寻优算法。

但是禁忌搜索算法也有明显的不足。

(1) 对初始解有较强的依赖性，好的初始解可使TS在解空间中搜索到好的解，而较差的初始解则会降低TS的收敛速度。

(2) 迭代搜索是串行的，仅是单一状态的移动，而非并行搜索。

为了进一步改善禁忌搜索的性能，一方面可以对禁忌搜索算法本身的操作和参数选取进行改进，另一方面则可以与模拟退火、遗传算法、神经网络以及基于问题信息的局部搜索相结合。

5.2禁忌搜索算法设计原则

与遗传算法、模拟退火法类似，禁忌搜索算法的运行效果也在很大程度上受其有关参数的影响。禁忌搜索算法的特征由禁忌对象和长度、候选集和评价函数、停止规则和一些计算信息组成。禁忌表特别指禁忌对象及其被禁的长度。禁忌对象多选择造成解变化的状态。候选集中的元素依评价函数而确定，根据评价函数的优劣选择一个可能替代被禁对象的元素，同时还需考虑禁忌规则和其他一些特殊规则。

作为一种人工智能算法，禁忌搜索算法的实现技术是算法的关键。下面按照禁忌对象、候选集合的构成、评价函数的构造、特赦规则、记忆频率信息和终止规则等分别进行介绍。

1. 禁忌对象、长度和候选集


禁忌表中的两个主要指标是禁忌对象和禁忌长度。禁忌对象指禁忌表中被禁的那些变化元素。因此，首先需要了解状态是怎样变化的。我们将状态的变化分为解的简单变化、解向量分量的变化和目标值变化三种情况。在这三种变化的基础上，讨论禁忌对象。下面还将同时介绍禁忌长度和候选集确定的经验方法。

1) 解的简单变化

这种变化最为简单。假设x,y∈D，其中D为优化问题的定义域，则简单解变化如下： 

x→y

是从一个解变化到另一个解。这种变化在局部搜索算法中经常采用，这种变化将问题的解看成变化最基本因素。

2) 向量分量的变化

这种变化考虑得更为精细，以解向量的每一个分量为变化的最基本因素。设原有的解向量为(Xl…Xl－1，Xl+1…Xn)，用数学表达式来描述向量分量的最基本变化为 (Xl…Xi－1，Xi，Xl+1…Xn)→(Xl…Xi－1，Yi，Xl+1…Xn)，即只有第i个分量发生变化。向量的分量变化包含多个分量发生变化的情形。部分优化问题的解可以用一个向量形式x－(X1，X2，…，Xn)r∈{0，1}n表示。解与解之间的变化可以表示某些分量的变化，如用分量X1=0变化为X1=1或从XK=1变化为XK=0，或是两者的结合。

3) 目标值变化

在优化问题的求解过程中，我们非常关心目标值是否发生变化，是否接近最优目标值。这就产生一种观察状态变化的方式： 观察目标值或平均值的变化。就犹如等位线的道理一样，把处在同一等位线的解视为相同。

II(a)={x∈D|f(x)=a}

其中，f(x)为目标函数。其表面是两个目标值的变化，即从a→b，但隐含着两个解集合的各种变化x∈II(a)→y∈II(b)的可能。

4) 禁忌对象的选取

由上面关于状态变化三种形式的讨论，禁忌的对象可以是其中的任何一种。在实际应用中，应根据具体问题采用一种方法。由于解的简单变化比解分量的变化和目标值变化的受禁忌范围要小，因此可能造成计算时间的增加。解分量的变化和目标值变化的禁忌范围大，减少了计算时间，可能引发的问题是禁忌范围太大导致陷入局部最优解。


因为NPhard问题不可能奢望计算得到最优解，因此在算法的构造和计算过程中，要尽量少地占用机器内存，这就要求禁忌长度尽可能短、候选集合尽可能量小。需要注意的是，禁忌长度过短易造成搜索的循环，候选集合过小易造成过早陷入局部最优。

5) 禁忌长度的确定


禁忌长度是被禁对象不允许选取的迭代次数。一般是给被禁对象X一个数(禁忌长度)t，要求对象X在t步迭代内被禁，在禁忌表中采用tabu(X)=t记忆，每迭代一步，该项指标做运算tabu(X)=t－1，直到tabu(X)=0时解禁。于是，可将所有元素分成两类，被禁元素和自由元素。有关禁忌长度t的选取，可以归纳为下面几种情况。

(1) t为常数，如t=0、t=n，其中n为邻域解的个数。这种规则容易在算法中实现。

(2) t=［tmin，tmax］，此时t是可以变化的数，它的变化依据是被禁对象的目标值和邻域的结构。此时tmin、tmax是确定的，确定tmin、tmax的常用方法是根据问题的规模T，限定变化区间［αt，βt］(0<α<β)； 也可以用邻域中解的个数n确定变化区间［αn，βn］(0<α<β)。当给定了变化区间时，确定t的大小主要依据实际问题、实验和设计者的经验。直观上，当函数值下降较大时，可能谷较深，欲跳出局部最优，希望被禁的长度较大。

(3) tmin、tmax的动态选取。有的情况下，用tmin、tmax的变化能达到更好的解。它的基本思路与(2)类似。

禁忌长度的选取同实际问题、实验和设计者的经验有紧密的联系。同时，它决定了计算的复杂性，过短会造成循环的出现，过长又造成计算时间较长。

6) 候选集合的确定

候选集合由邻域中的解组成。常规方法是从邻域中选择若干目标值或平均值最佳的解入选。有时认为计算量还是太大，则不在邻域中的解中选择，而是在邻域中的一部分解中选择若干目标值或评价值最佳的解入选。也可以用随机选取的方法实现部分解的选取。

2. 评价函数

评价函数是候选集合元素选取的一个评价公式，候选集合的元素通过评价函数值来选取。以目标函数作为评价函数是比较容易理解的。目标值是一个非常直观的指标，但有时为了方便或易于理解，会采用其他函数来取代目标函数，可以将评价函数分为基于目标函数的评价函数和其他方法两类。

1) 基于目标函数的评价函数

这一类主要包括以目标函数的运算所得到的评价方法，如记评价函数为p(x)，目标函数为f(x)，则评价函数可以采用目标函数p(x)=f(x)； 目标函数值与xnow目标值的差值p(x)=f(x)－f(xnow)，其中xnow是上一次迭代计算的解； 目标函数值与当前最优解目标值的差值p(x)=f(x)－f(xbest)，其中xbest是目前计算中的最好解。

基于目标函数的评价函数的形成主要通过对目标函数进行简单的运算，它的变形有很多。

2) 其他方法

有时计算目标值比较复杂或耗时较多，解决这一问题的方法之一是采用替代的评价函数。替代的评价函数还应该反映原目标函数的一些特征，如原目标函数对应的最优点还应该是替代函数的最优点。构造替代函数的目标是减少计算的复杂性。具体问题的替代函数构造依问题而定。

3. 特赦规则

在禁忌搜索算法的迭代过程中，会出现候选集中的全部对象都被禁忌，或有一对象被禁，但若解禁则其目标值将有非常大的下降情况。在这种情况下，为了达到全局的最优，我们会让一些禁忌对象重新可选。这种方法称为特赦，相应的规则称为特赦规则。常用的特赦规则如下： 

1) 基于评价值的规则

在整个计算过程中，记忆已出现的最好解xbest。当候选集中出现一个解xnow，其评价值(可能是目标值)满足c(xbest)>c(xnow)时，虽然从xbest达到xnow的变化是被禁忌的，此时，解禁xnow使其自由。直观理解，我们得到一个更好的解。

2) 基于最小错误的规则

当候选集中所有的对象都被禁忌时，而1)的规则又无法使程序继续下去。为了得到更好的解，从候选集的所有元素中选一个评价值最小的状态解禁。

3) 基于影响力的规则

有些对象的变化对目标值影响很大，而有的变化对目标值变化较小。我们应该关注影响力大的变化。

4. 记忆频率信息

在计算的过程中，记忆一些信息对解决问题是有利的。如一个最好的目标值出现的频率很高，这使我们有理由推测： 现有参数的算法可能无法再得到更好的解，因为重复的次数过高，使我们认为可能出现了多次循环。根据解决问题的需要，我们可以记忆解集合、有序被禁对象组、目标值集合等的出现频率。一般可以根据状态的变化将频率信息分为静态和动态两类。

静态的频率信息主要是某些变化，例如解、对换或目标值在计算中出现的频率。求解它们的频率相对比较简单，如可以记录它们在计算中出现的次数，出现的次数与总的迭代数的比率，从一个状态出发再回到该状态的迭代次数等。这些信息有助于我们了解一些解、对换或目标值的重要性，是否出现循环和循环的次数。在禁忌搜索中，为了更充分地利用信息，一定要记忆目前最优解。

动态的频率信息主要是从一个解、对换或目标值另一个解、对换或目标值的变化趋势，如记忆一个解的序列的变化，或记忆一个解序列变化的若个点等。由于记录比较复杂，因此它提供的信息量也较大。在计算动态频率时，通常采用的方法如下： 

(1) 一个序列的长度，即序列中的元素个数。在记录若干关键点的序列中，按这些关键点的序列长度的变化进行计算。

(2) 从序列中的一个元素出发，再回到该序列该元素的迭代次数。

(3) 一个序列的平均目标(评价)值，从序列中一个元素到另一个元素目标(评价)值的变化情况。

(4) 该序列出现的频率。

5. 终止规则

无论如何，禁忌搜索算法是一个启发式算法。我们不可能让禁忌长度充分大，只希望在可接受的时间中给出一个满意的解。于是，很多直观、易于操作的原则包含在终止规则中。常见终止规则如下： 

1) 确定步数终止

给定一个充分大的数N，总的迭代次数不超过N步。即使算法中包含其他的终止原则，算法的总迭代次数有保证。这种原则的优点是易于操作和可控计算时间，但无法保证解的效果。在采用这个规则时，应记录当前最优解。

2) 频率控制终止

当某一个解、目标值或元素序列的频率超过一个给定的标准时，如果算法不做改进，只会造成频率的增加，此时的循环对解的改进已无作用，因此，终止计算。这一规则认为如果不改进算法，解不会再改进。

3) 目标值变化控制原则

在禁忌搜索算法中，提倡记忆当前最优解。如果在一个给定的步数内，目标值没有改变，同2)相同的观点，如果算法没有其他改进，解不会改进。此时，停止运算。

4) 目标值偏离程度原则

对一些问题可以简单地计算出它们的下界(目标为极小)，目标计算得到的解与最优值很接近，终止计算。


5.3禁忌搜索算法的应用
5.3.1禁忌搜索算法在旅行商问题中的应用

旅行商问题(TSP)是经典的组合优化问题。假设有若干城市，任意两个城市之间的距离已知，则旅行商问题可简单描述如下： 一个旅行商从任意一个城市出发，在访问完其余城市后(每个城市只被访问一次，不允许多次访问)，最后返回出发城市，找到旅行商所行走的最短路线。

相关MATLAB代码见目录二维码中的附录7。

5.3.2禁忌搜索算法在双层级医疗设施选址问题中的应用

目前，不论是在城市还是农村，都存在着这样一个群体，他们害怕看病，更害怕进医院，于是就小病撑，大病扛。而纵观全国各大医院，更是人满为患，门诊“一号难求”，住院“一床难求”，“看病难”已是亟待解决的问题之一。 “看病难”这一问题出现的很大一部分原因在于医疗资源不足，配置不均。

医疗设施合理的层级布局可以有效分流患者，使大医院人满为患及患者“看病难”问题得到一定解决，使小医院的优质资源与设备的利用率得到提高，并对完善“大病进医院、小病进社区”的医疗模式具有很好的推动作用。

为了完善不同层级医院之间的分层就诊制度，以及提高医疗资源的整体利用效率，将以双层级多样流的医疗服务设施选址问题为研究对象，构建考虑距离、容量等多方面因素的双层级医院设施选址模型，为每个服务对象指派两个层级的医疗设施，使患者的需求在可以自行选择的情况下得到满足。

1. 分级诊疗问题描述

结合分级诊疗政策，将基层社区卫生服务中心及区域医疗中心分别作为第一层级医疗设施和第二层级医疗设施。区域医疗中心承担的一般性门诊、康复和护理等诊疗服务可分流到社区卫生服务中心，同时区域医疗中心担负与社区卫生服务中心上下联动、协同共进的任务，提供必要的转诊服务。患者首先选择基层社区卫生服务中心，检查后有一定比例的患者需要转诊到区域医疗中心继续治疗。这样原来作为供给点的基层社区卫生服务中心成为了需求点，新的供给点为区域医疗中心。但是，现实生活中还存在其他两种情况。

(1) 有一定比例的患者病况严重，需要直接选择区域医疗中心进行治疗。

(2) 有一定比例的患者虽然仅需要社区卫生服务中心可以提供的医疗服务，但其直接选择区域医疗中心进行诊疗，至此形成了多样流嵌套型层级系统，该系统如图51所示。



图51多样流嵌套型层级系统


2. 模型构建


设I={1，2，…，e}为需求点的集合。M={1，2，…，m}表示第一层级候选设施集合，共m个候选点； N={1，2，…，n}表示第二层级候选设施集合，共n个候选点。选址模型目的为在M中选择p个第一层级设施和在N中选择q个第二层级设施作为服务设施。

定义5个变量，含义如下： xj代表01变量，其取值为1时，第一层级候选设施j被选中设立； yk代表01变量，其取值为1时，第二层级候选设施k被选中设立； uij代表需求点与第一层级设施之间的流动人数； vjk代表第一层级设施j与第二层级设施k之间的流动人数； zik代表需求点i与第二层级设施k之间的流动人数。

其他参数定义如下： wi代表需求点处的需求量； dij代表需求点i与第一层级设施j之间的单位流动费用(元/人)； hik代表需求点i与第二层级设施k之间的单位流动费用(元/人)； rik代表第一层级设施j与第二层级设施k之间的单位流动费用(元/人)； fj代表第一层级候选设施j的费用； gk代表第二层级候选设施k的费用； θ代表需求点需要基础医疗服务的患者首选第一层级设施的比例； β代表被第一层级设施服务后的患者需要前往第二层级设施的比例； μ代表需求点居民需要基础医疗服务的患病率； α代表需求点居民需要第二层级医疗服务的患病率； c1j代表第一层级设施j的容量； c2k代表第二层级设施k的容量。

模型如下： 


minZ=∑ei=1∑mj=1dijuij+∑mj=1∑ni=1rjkvjk+∑ei=1∑nk=1hikzik+∑mj=1fjxj+∑nk=1gkyk(51)
s.t.∑mj=1uij=θ·μ·wi，i∈I(52)
∑nk=1vjk=β∑ei=1uij，i∈M(53)
∑nk=1zik=αwi+(1－θ)·μwi，i∈I(54)
∑ei=1uij≤c1jxj，j∈M(55)
∑mj=1vjk+∑ei=1zik≤c2kyk，k∈N(56)
∑mj=1xj=p，j∈M(57)
∑nk=1yk=q，k∈N(58)
xj={0，1}，j=1，2，3，…，m(59)
yk={0，1}，k=1，2，3，…，n(510)
uij≥0且为整数，i∈I，j∈M(511)
vjk≥0且为整数，j∈M，k∈N(512)
zik≥0且为整数，i∈I，k∈N(513)


其中，式(51)表示最小化总成本，包括第一、二层级的建设成本和需求点与层级之间流动成本。式(52)表示第一层级设施为需求点提供服务的总人数。式(53)表示由第一层级设施j服务后的患者需要前往第二层级设施就诊的人数。式(54)表示需求点i直接前往第二层级设施就诊的患者数。式(55)和式(56)为两个层级设施的容量限制，并确保了只有开放的设施才能提供医疗服务。式(57)和式(58)确定了第一层级和第二层级所建立设施数量。式(59)~式(513)为5个决策变量的取值范围。

3. 算法设计

本案例利用禁忌搜索算法求解上述双层级医疗设施选址模型，现将算法中涉及的编码方式、初始解产生方式、邻域结构计算规则、禁忌列表与藐视准则、适应度值计算方法等关键操作介绍如下。



图52编码方式

1) 编码方式

根据双层级选址模型的特征，采用01编码，解的编码长度为m+n。例如，图52中第一、二层级候选设施个数分别为6和4，则解的编码长度为10，且所有第一层级设施都是在第二层级设施之前排列。当解某一个位置上的值为1，表示该候选点的设施开放，否则关闭。图52中的编码表示第一层级编号为1、2、3的三个候选设施开放和第二层级编号为4的候选设施开放。


2) 初始解产生方式

以往研究表明： 以设施的成本效益值为标准选择开放候选设施，能够有效提升初始解的质量。为此，将固定成本和流动成本考虑在一起，设施的成本效益值定义为设施的固定成本与流动成本之和与设施容量的比率，则第一层级设施j的成本效益值表示为式(514)； 第二层级设施k的成本效益值表示为式(515)。分别选择两个层级成本效益值最低的前p和q个设施开放，若开放设施容量之和不能满足总需求，则利用邻域解的产生方式产生邻域解，在产生的p(m－p)q(n－q)个邻域解中选取一个可行解作为初始解。


fi+∑ei=1dijc1j(514)
gk+∑mj=1rjk+∑ei=1hijc2k(515)


3) 邻域结构计算规则

对当前解进行邻域搜索： 在开放设施中依次选择一个，将其关闭，同时在未开放设施中依次选择一个，将其开放。由于双层级设施选址的特殊性，需要对当前解中的每一个层级分别执行以上操作。对第一层级来讲，每个当前解对应p(m－p)个邻域解； 对第二层级来讲，每个当前解对应q(n－q)个邻域解。将两个层级的邻域解组合，产生p(m－p)q(n－q)个不同的邻域解，并从中剔除总容量不满足总需求的解； 为保证邻域解的质量，选取总成本效益值最低的前若干个解，然后利用CPLEX求解适应度值。

4) 禁忌列表与藐视准则

在双层级设施选址问题中，当前解与其邻域解相比，每层级均有一个不同的开放设施，因此将其作为禁忌对象列入禁忌列表，且禁忌长度设置为偶数。例如，禁忌长度为4，禁忌任期为2，当前禁忌列表为［2，6，3，9］，若此时新的禁忌对象为［4，7］，则禁忌列表更新为［3，9，4，7］。在迭代过程中，当候选解集中的全部或某一对象被禁忌，此时需采用藐视准则将禁忌对象解禁。本案例采用基于适应度准则，若某个禁忌候选解的适应度值优于当前最佳解，则将对应对象解禁，并更新当前解及最佳解。

5) 适应度值计算方法

确定了变量xi、yk后，模型的约束条件及决策变量的个数都得到一定减少。本案例在禁忌搜索过程中，利用CPLEX求解缩减后的模型，从而得到整数变量uij、vjk、zik的值，并将原模型的目标函数值作为当前解的适应度值。适应度值越优(小)，则解的质量越高。

6) 算法步骤

步骤1： 初始化参数。输入模型参数dn、rik、hik、M、N、β、θ等的值，并初始化算法迭代次数、禁忌搜索最大迭代次数以及禁忌长度。

步骤2： 按照式(514)和式(515)分别计算两个层级候选设施的成本效益值。

步骤3： 产生初始解。采用01编码方式，生成一个满足总需求量的初始解Sol。

步骤4： 判断算法是否满足终止条件满足。若是，结束迭代并输出优化结果； 否则，转步骤5。

步骤5： 利用当前解的邻域函数产生满足容量限制及设施数量的候选解。

步骤6： 判断候选解是否满足藐视准则。若成立，则用满足藐视准则适应度值最佳的候选解作为新的当前解，并将与新的当前解对应的禁忌对象加入禁忌列表，更新禁忌列表，转步骤4； 否则，转步骤7。

步骤7： 判断候选解对应的各对象的禁忌属性，选择候选解集中非禁忌对象对应的最佳状态为新的当前解，同时将与新的当前解对应的禁忌对象加入禁忌列表，更新禁忌列表，并转步骤4。

4. 算例分析

1) 算例说明

选择上海市某区社区卫生服务机构以及区域医疗中心作为研究主体。根据《2017年上海市统计年鉴》及各医院官网数据，该区的各级医疗设施分布如图53所示。


由图53可知，该区的12个街道(镇)中，街道2、3、6、11医疗资源比较充足，都拥有三个医院且至少有一个二级综合医院； 街道8、9、10、12的医疗资源配置量相对较低，都拥有两个医院； 街道1、4、5、7资源数量少于其他街道。因此对该区域医疗设施选址布局进行科学规划，改善医疗设施分布状况，有助于提高其医疗资源的配置水平。



图53上海市某区各级医疗设施分布


把该区的12个街道(镇)看作12个需求点，几何中心视为人口中心。根据《上海市医疗机构设置“十三五”规划》中对床位的规定标准： 每千人配备5.2张床位。将每个街道的居住人口量换算成所需床位量并作为需求点的需求量，具体信息如表51所示。另外，考虑到居民出行情况，以欧氏距离作为需求点与设施点的距离不符合实际情况，因此根据百度地图规划得到居民的出行路线。


表51需求点及其需求量



需求点编号人口(万人)需求量需求点编号人口(万人)需求量


1
2.725
142
7
9.033
470

2
19.255
1002
8
7.019
365
3
17.899
931
9
9.538
496
4
14.909
776
10
8.587
447
5
9.250
481
11
12.495
650
6
10.561
550
12
10.048
623


选择该区基层社区卫生服务机构作为第一层级候选设施点，包括社区卫生服务中心、社区卫生服务站。通过实地考察分析，剔除了一些硬件设施远不及社区卫生服务中心标准的候选设施点； 考虑到街道1、4、5、7等的医疗资源相对匮乏，新增该区域的家庭医生工作站和社区卫生服务中心延伸点作为候选点，最终得到41个候选设施点，候选点分布如图55所示。将床位数作为候选点的容量，根据城市社区卫生服务中心建设标准，每个中心床位数设置为20~99张。候选点设施的固定成本根据其床位数确定，平均每个床位占地30m2，对应建造成本为2000元/m2，第一层级候选点的需求量如表52所示。

选择该区二级及以上综合医院作为第二层级候选设施点，最终得到10个满足硬件设施要求的综合医院候选点，其分布如图54所示。每个候选点的固定成本根据该设施的床位数确定，平均每个床位占地45m2，对应建造成本为2000元/m2。第二层级候选点的需求量如表53所示。



图54医疗设施候选点分布



表52社区卫生服务中心候选点容量(床位数)



序号容量序号容量序号容量序号容量


1
36
12
78
23
60
34
60

2
45
13
72
24
72
35
96
3
45
14
84
25
60
36
79
4
87
15
48
26
84
37
84
5
74
16
60
27
72
38
60
6
78
17
72
28
60
39
72
7
79
18
78
29
240
40
72
8
72
19
90
30
60
41
78
9
42
20
48
31
48


10
62
21
66
32
48


11
60
22
78
33
86








表53区域医疗中心候选点容量(床位数)



序号容量序号容量


1
600
6
640

2
880
7
600
3
550
8
980
4
650
9
600
5
500
10
1012


根据《上海市医疗机构设置“十三五”规划》（以下简称《规划》）关于基层医疗卫生机构设置的规定，社区卫生服务机构按每新增5~10万居住人口增设1家社区卫生服务中心或分中心。根据第六次人口普查，该区常住人口数为131.32万，应在该区设立15个社区卫生服务中心，以承担一般常见病、多发病和诊断明确慢性病的初级诊疗及转诊服务。《规划》明确指出，按每30~50万人口设立1家区域医疗中心。因此应在该区建立3个区域医疗中心，以承担一般疑难杂症和转诊患者的治疗。

2) 算例求解

利用MATLAB 2016b对TS算法进行编程实现，算法参数设置如下： 患者首选第一层级设施的比例θ为0.55； 需求点居民需要基础医疗服务的患病率μ为0.186； 需求点居民需要第二层级医疗服务的患病率α为0.13； 第一层级设施服务后的患者需要至第二层级设施的比例β为0.1； 最大迭代次数Gmax为50； 禁忌长度len为10。

禁忌搜索过程中，选择合适长度的禁忌列表有助于提高算法的搜索效率。本案例将不同禁忌长度的算法求解结果进行对比，其结果如图55所示，禁忌长度为10时算法的总体迭代次数最少，运行结果较优。



图55不同禁忌长度的算法迭代情况


5. 案例总结

合理布局不同等级的医疗设施对分流患者、提高医疗资源的利用率及完善“大病进医院、小病进社区”的医疗模式具有积极的推动作用。本案例以街道作为资源配置的基本单位，均衡配置区域医疗资源，建立了更加适合分级诊疗政策的双层级多样流选址模型，并针对模型的具体特点，设计了求解该问题的禁忌搜索算法。最后以上海市某区社区卫生服务中心为第一层级医疗设施和区域医疗中心为第二层级医疗设施，对其布局进行选址优化，计算结果表明，优化后的设施布局更加合理，能够有效分流患者，满足居民的就医需求，减少了医疗资源的浪费。

5.3.3禁忌搜索算法在机场外航服务人员班型生成问题中的应用

机场外航服务人员班型生成问题是实现机场人员智能化排班的重要组成部分，主要是在满足各类约束条件的情况下生成一个安排员工完成一天全部任务的班型方案。有效的班型方案能降低后续人员排班过程中的复杂度，保证排班结果的稳定性，且通过班型方案管理人员可以了解完成现有任务所需人力资源情况。班型生成问题是一个NPhard问题，大规模、多目标及复杂的劳动法规约束使得该问题的求解非常困难。

1. 问题描述

本案例以首都机场外航服务人员班型生成问题为例，对多任务层次资质场景下的班型生成过程进行了阐述，提出了相应的优化模型和算法，求得了最优的班型方案，为机场外航服务部的班型生成问题提供了理论依据。

1) 多任务

外航服务部主要为停靠在机场的所有国际航班提供包括值机、接机、送机和签证审查等在内的多种地面保障服务。为了保证服务质量，航空公司与外航服务部签订了服务水平协议，详细规定了不同保障任务要求的到岗时间、离岗时间及员工需求情况等。

2) 层次资质

由于教育水平及工作经验等使得员工在执行不同任务时具有不同层次资历，外航服务部员工资历分为3个层次等级，分别是组长、控制人员和普通人员。其中，组长资历等级最高，控制人员次之，普通人员最低，且资质存在向下兼容的层次关系，即高等级资历的员工能够向下兼容执行低等级资历的任务，反之则不能。

3) 班型

班型是指员工一天的工作任务，在多任务层次资质场景下班型需要给出员工的到岗时间和离岗时间、保障的任务集合、对保障任务需要具备的层次资质等信息。

面向多任务层次资质场景下的班型生成问题是： 给定一天内的多任务集合和员工层次资质集合，生成一个优化的班型方案，要求该班型方案能够覆盖全部任务且每个班型需满足约束条件，优化目标为最小化班型方案中的总工作时间。

2. 模型建立

1) 符号说明


集合K={1,2,…,k}表示外航服务保障的多种任务类型，S={1,2,…,s}表示航班保障涉及员工的层次资质种类，|·|表示集合维度，即集合中元素的个数。

多任务集合T={(ki,hi,fi,Mi)|i=1,2,…,n}，其中n表示任务总数，ki表示任务类型，hi表示任务到岗时间，fi表示任务离岗时间，Mi={ms|s=1,2,…,|S|}表示任务所需各层次资质的员工数量集合，ms∈N表示需类资质员工的数量。

员工对各种任务的层次资质集合P={(p(k,s)d)|d=1,2,…,m； k=1,2,…,|K|； s=1,2,…,|S|}，其中m表示员工总数，p(k,s)d∈{0,1}表示员工d对任务类型k具备s级资质信息，取值为1表示具有该类资质，否则取值为0，因员工资质存在向下兼容的层次关系，故当p(k,s)d=1时，p(k,s－1)d=p(k,s－2)d=…=p(k,1)d=1。

优化变量为班型集合W={sj,cj,Tj,Pj|j=1,2,…,t}，其中t表示班型总数，sj表示班型到岗时间，cj表示班型离岗时间，Pj={p(k,s)|k∈Kj,s∈S}表示保障该班型的员工对班型内各任务需具备的资质信息，p(k,s)d∈{0,1}，取值为1表示要求服务该班型的员工对任务k必须具有s类的资质，TjT表示班型保障的任务集合。规定sj为保障任务集合Tj中最早的到岗时间，即sj=minhi,i∈Tj，cj为保障任务集合Tj中最晚的离岗时间，即cj=maxfi,i∈Tj。根据班型集合W可以得到xij∈{0,1}，表示班型j是否服务任务i，取值为1表示服务。

由于员工对不同的任务类型具有不同的层次资质，且不同任务类型之间的资质不存在交集，因此在生成班型时需考虑能否找到满足班型资质要求的员工。用pjd表示员工d对班型j是否具有服务资质，pjd∈{0,1}，员工d满足班型j的资质要求，取值为1，且员工d称为班型j的候选员工。为了使得求得的班型更具有实际意义，规定班型的候选员工数量必须大于K。


2) 优化目标及约束

结合机场外航服务部的实际业务，建立面向多任务层次资质场景下的班型生成模型如下： 

minxij∈{0,1}∑tj=1(cj－sj)(516)

s.t.


∑tj=1xij·p(k,s)=∑s1Msi; i=1,2,…,n; s=1,2,…,|S|(517)
xij·fi－xi′j·h′i>0,j; i,i′∈Tj; i≠i′(518)
Imin≤cj－sj≤Imax,j∈W(519)
pjd=1,∑|S|s=1p(ki,s)d≥∑|S|s=1p(ki,s),i∈Tj0,≤其他(520)
∑Md=1pjd≥K,j∈W(521)


其中，式(516)表示最小化班型方案的总工作时间。式(517)表示班型方案需覆盖全部的任务集合，即生成的班型方案能满足每个任务要求的各层次资质的人数。式(518)表示同一班型内保障任务的工作时间不能冲突，即任意两个任务的到岗时间和离岗时间不能有重叠。式(519)表示班型的服务时长应在给定的合法区间内。式(520)是pjd的计算方式，即当员工对班型保障任务集合中的所有任务具有的资质都高于或等于班型要求的服务人员的资质时，pjd取值为1，否则为0。式(521)规定班型候选员工数量必须大于或等于K，以保证所生成班型在进行人员安排时能找到具备服务该班型资质的员工。

3. 模型求解

上面介绍的班型生成优化模型的求解是一个集合划分问题，属于NPhard问题。传统算法很难在多项式时间内求得较优的解。禁忌搜索算法具有灵活的记忆功能和藐视准则，在搜索过程中可以接受劣质解，表现出较强的“爬山”能力。另外，禁忌搜索算法能够跳出局部最优解，转向解空间的其他区域，增大获得全局最优解的概率。故本案例基于禁忌搜索算法对所建立的面向多任务层次资质场景下的班型生成模型进行求解。

1) 多任务复制转换

在面向多任务层次资质场景下的班型生成问题中，任务不同所需各层次资质的人数不同，增加了邻域移动设置的困难性。为简化邻域移动设置，本案例提出如下复制方法对多任务集合进行等价转换： 任务i需各层次资质服务总人数p，则将任务i复制p次，且设置复制后的任务需要的服务人数为1，即用p个所需服务人数为1的等价任务替换1个需要服务人数为p的任务，复制过程保证了排班任务量不变。如表54所示，任务GA891复制为表55所示的5个等价任务。


表54原始任务



星期航班号开始时间结束时间需组长人数需控制人员人数
需普通人员人数


1
GA891
5:40
8:50
1
1
3



表55复制以后的等价任务



星期航班号开始时间结束时间需组长人数需控制人员人数需普通人员人数


1GA8915∶408∶50100
1GA8915∶408∶50010
1GA8915∶408∶50001
1GA8915∶408∶50001
1GA8915∶408∶50001


基于复制以后的等价任务，设计禁忌搜索算法还需解决的关键问题包括初始解生成、邻域移动和适应值函数选择等。

2) 初始解生成

模型优化最终结果对初始解依赖性很大，一个好的初始解，不仅能加快算法收敛速度，也更容易找到全局最优的解。本案例基于贪婪算法生成初始解，具体步骤如下。


步骤1： 将集合T中的任务按照到岗时间hi先后进行排序。

步骤2： 从T中的第1个任务t1开始，找到满足式(518)和式(519)的任务，形成备选任务集合T′。

步骤3： 遍历备选任务集合T′，按式(516)计算目标值，并找到目标值最小的任务形成班型j并加入班型集合W。

步骤4： 在任务集合T中删除班型j包含的所有任务。

步骤5： 如果任务集合T为空，得到初始解班型集合W，结束算法，否则转步骤2。

需要说明的是，通过上述算法生成的初始解满足式(517)、式(518)、式(519)，不满足式(521)，主要基于以下考虑。

(1) 式(521)需要检查所有员工对不同任务的资历信息，计算量较大。

(2) 在满足式(54)的情况下，班型保障任务较少，很大概率满足式(521)。

综合(1)(2)，为了加快初始解的生成速度，在生成初始解时未考虑式(521)。由于概率极小时可能不满足式(521)，因此在随后的禁忌搜索算法设计中，加入了对该状况的巨大惩罚。


3) 邻域移动

采用交换移动和插入移动两种邻域移动方式生成候选解。

插入移动是将任务i插入班型j的保障任务集合Tj中，其过程如图56所示。



图56插入移动


交换移动是将班型j中的任务i与班型j′中的任务i′进行交换，其过程如图57所示。



图57交换移动


进行插入移动和交换移动时需满足上述模型中的约束条件。

4) 适应值函数选择

适应值函数用于对搜索结果进行评价，本案例的适应值函数包括目标函数和惩罚函数两部分。惩罚函数也由两部分组成。

(1) 为了增加模型在资源不足情况下的适应性，将式(518)的约束条件设置为软约束。结合实际设置了如下惩罚函数： 班型的服务时长在法定区间之内，惩罚值取0； 班型的服务时长超出法定区间，惩罚值为服务时长到法定区间的中心距离，距离越远惩罚值越大，记为

g(x)=0,Imin≤x≤Imax
α·|x－d|,其他
x=cj－sj,d=Imax－Imin2+Imin(522)

(2) 若班型不满足式(521)，则惩罚值为一个较大的惩罚常数β，记为

h(j)=β，班型j不满足式(521)
0，其他(523)

结合式(516)、式(522)、式(523)，本案例的适配值函数设定如下： 

f(W)=∑tj=1g(wj)+h(wj)+χ·(cj－sj)(524)

5) 其他设置

(1) 禁忌表。禁忌表用于防止搜索过程中出现死循环，避免陷入局部最优解。禁忌对象设为每次发生邻域移动时的任务。候选解大小设置为固定值，本案例采用前人证明过的禁忌长度α=A时算法收敛速度和求得的解是最好的方式。

(2) 特赦准则。选取基于评价值的规则作为藐视准则，即当前解的适应值优于最优值，则满足特赦准则，更新最优值为当前值。

(3) 终止条件。终止准则采用固定步长的方式，达到设定的最大步长M终止算法。

6) 基于禁忌搜索算法的模型求解

基于禁忌搜索算法对模型进行求解的流程如图58所示。



图58基于禁忌搜索算法对模型进行求解的流程


4. 实验及分析

1) 实验数据

实验数据选用首都机场外航服务部2018年1月15日至2018年1月21日累计1周的生产数据，共包括100个任务信息如表56所示，84名员工在10家航空公司的4种层次资质信息如表57所示。为描述方便，采用3表示组长，2表示控制人员，1表示普通员工，0表示没有服务资质，且组长资质最高，控制人员次之，普通最低。


表56排版任务信息数据表



星期航班号开始时间结束时间需组长人数需控制人员人数
需普通人员人数


1
GA891
5:40
8:50
1
1
3

1
KE856
9:00
14:15
1
1
3
1
JS152
9:30
12:55
1
1
2
…
…
…
…
…
…
…




表57员工资质信息



员工姓名
航空公司类型


AA
KE
SU
GA
PK
J2
IR
7C
VN
JS

李晓敏
0
2
0
1
0
0
0
2
2
2
曹曦月
0
2
0
2
0
0
0
2
1
2
李艳香
2
0
0
3
0
0
0
0
3
1
︙
︙
︙
︙
︙
︙
︙
︙
︙
︙
︙



2) 实验结果

结合机场外航服务部的实际业务算法参数设置如下： 班型的最长工作时长Imax=9h； 班型的最短工作时长Imin=6h； 候选比例值K=10； 惩罚值α=100，β=1000，χ=100； 最大迭代次数M=100； 候选解大小A=100； 禁忌表长度α=10。


实验部分结果数据如表58所示。表58中一行代表1个班型，规定了到岗时间、离岗时间、保障的任务(航班号)及其任务要求的层次资质。例如星期一的第1个班型要求7：50到岗，11：45离岗，保障任务KE880、SU2852，且要求保障人员对任务KE必须具备普通资质，对任务SU必须具备控制人员资质。


表58班型方案示例



星期到岗时间离岗时间保障任务集合(航班号)及资质要求


17：5011：45KE880/1SU2852/2
15：0010：10J2067/1AA186/1
15：2014：15PK852/2KE856/1
︙︙︙︙


3) 实验结果分析

(1) 班型分析。为了验证算法的有效性，将本案例实验得到的班型方案与首都机场外航服务部现有人工生成的班型方案从总服务时间和总服务人数两方面进行对比，如表59所示。


从表59可看出，本案例模型求得的每天的班型方案在优化目标班型总工作时间上相比较于人工班型方案都有不同程度的降低，且在减少班型的总服务时间的同时并没有增加班型的服务人数，有效地提高了资源利用率。同时由于约束复杂，任务量大，人工生成班型方案中存在部分班型不能满足硬约束条件的情况，而本案例模型生成的班型方案中所有班型完全满足全部硬约束条件。


表59与人工班型对比结果



星期
总服务时间(h)总服务人数(人)
人 工 方 案本案例模型人 工 方 案本案例模型


1462.13444.416362
2
338.71
318.33
54
44
3
415.36
361.50
54
50
4
315.47
299.41
52
45
5
424.67
395.00
60
47
6
375.49
345.16
47
44
7
352.31
313.08
51
45




为了分析班型方案对班型工作时间的优化情况，统计了本案例实验得到班型方案和人工生成班型方案中服务时长在［0,6)，［6,9］，(9,11］，(11,13］区间的班型百分比，如表510所示。

从表510可以看出，人工班型方案服务时长在各时间范围内的分布比较分散，26.55%班型服务时长超过9h，有3.14%的班型服务时长超过11h。本案例模型求得的班型方案，49.55%班型服务时长在合法区间内，服务时长在(9,11］区间的班型较少，最长班型服务时长不会超过11h，其结果更高效更具人性化。


表510班型服务时长统计(%)



算法
班型服务时长区间(h)

［0,6)［6,9］(9,11］(11,13］


人工班型方案25.5639.0626.553.14
本案例模型生成的班型方案
35.60
49.55
14.84
0


(2) 算法特性分析。从算法的收敛性和稳定性分析算法的有效性，其结果如图59所示。


图59(a)给出了算法执行过程中每一步迭代对应的适应值，可以发现，在搜索初期，目标函数收敛速度较快，经过大概40次迭代以后，目标函数值趋于稳定，仅在一个很小范围内进行波动，可见本案例设计的算法具有良好的收敛性，同时也验证了迭代次数设置为100是合理的。智能优化算法具有一定的随机性，为了分析本案例算法的稳定性，图59(b)给出了算法多次运行时，求得的班型方案的总工作时间的变化情况，可以看出，算法多次运行时求得的总工作时间相差不大，表明算法具有较好的稳定性。

5. 案例总结

针对首都机场外航服务部班型生成面临的任务种类多、员工对任务具有层次资质、人工生成班型困难且资源利用率不高等问题，研究构建了面向多任务层次资历场景下的班型生成模型，并设计了禁忌搜索算法对模型进行求解。在首都机场外航服务部的实际数据集上进行实验，实验结果表明，提出的模型和算法很好地实现了多任务层次资质场景下的班型生成过程全自动化，且比较于人工生成的班型方案，所需工作总时间和总人数都有降低，提高了机场资源的利用率和运行效率，为机场智慧决策提供了理论依据。



图59算法有效性分析


5.4本章小结

本章先后介绍了禁忌搜索算法的基本思想和设计原则，并详细阐述了禁忌搜索算法在旅行商问题、双层级医疗设施选址问题以及机场外航服务人员班型生成问题中的应用。迄今为止，禁忌搜索算法在组合优化、生产调度、机器学习、电路设计和神经网络等领域取得了很大的成功。近年来，禁忌搜索算法在函数全局优化方面得到了较多的研究，并大有发展的趋势。


5.5习题

1. 简述禁忌搜索算法的原理。

2. 简述禁忌搜索算法的特点。

3. 总结禁忌搜索算法的优缺点。

4. 禁忌搜索有哪些关键参数？该如何设置？

5. 简述评价函数的作用。

6. 简述局部搜索算法的基本步骤。

7. 简述禁忌搜索算法求解旅行商问题的基本步骤。

8. 禁忌长度该如何设置？