禁忌搜索（Tabu Search）算法及matlab实现（非旅行商（TSP）例子）

本文为原创，未经同意，不得转载

最近在自学智能算法，学到禁忌搜索算法时，发现无论是寻找到的相关教材还是CSDN博客中，对于禁忌搜索问题大多是以旅行商城市旅行路径最短为例来讲述的，就产生一种禁忌搜索算法只能对TSP类型的问题进行优化，显而易见，这是一种粗浅的认识。

任何智能算法都是一种寻找目标最优的独特思想，也就是说任何一种智能算法在适当设置参数和编写代码后都可以解决各类问题，只是存在这种算法解决这类问题的难易罢了。

当你读到这篇文章时，如果你还完全不了解禁忌搜索算法，可以先看看相关禁忌搜索文献或者CSDN里其他博主关于禁忌搜索的博客，里面都很详细的禁忌搜索算法的基本原理，在本文中就不在赘述，笔者更想是谈一些有关禁忌搜索的自我理解。

好了，废话不多说，进入到本文的主体。

---

一、禁忌搜索算法概要

1.一些禁忌搜索算法的见解

禁忌搜索算法说白了就是一种通过在给定的邻域内，运用确定的规则进行搜索，并通过相应的适应度函数对所搜索到的点进行比较寻优，最终找到目标函数的最优解。

禁忌搜索算法比较依赖初始值和邻域生成规则，两者选取不当，可能会陷入到局部最优当中。

禁忌搜索算法的基本构成是：初始值、目标函数、邻域产生规则、邻域搜索规则、禁忌表、特赦规则以及适应度函数。

2.禁忌搜索算法的基本流程

①设置初始值inittial_point

②将初始值赋给now_point；同时计算初始值的适应度函数，并将算出的值当做当前的目标最优值now_best；

③判断是否达到终止条件，若没有达到进行下步，否则数据输出。

④运用邻域产生规则，利用now_point产生当前的搜索邻域now_area；

⑤运用邻域搜索规则与禁忌表，产生当前邻域now_best的搜索点序列now_point；

⑥对当前点序列进行适应度函数的相关计算后，得到此次相应点的适应度序列shiyi；同时将算完适应度后的当前点序列，通过特赦规则放入到禁忌表中。

⑦如果产生的适应度序列中有优于当前目标最优值now_best的值H，将H赋值给now_best，同时将H点的信息point_H赋值给now_point后，并返回到第③步

⑧如果产生的适应度序列没有优于当前目标最值now_best，返回到第⑤步，即继续在now_area剩余的没有搜索过的邻域中选点。最终会再次出现两种情况：A、搜索到适应度值优于now_best的点H，则返回到第⑦步；B、now_area中所有的点都搜索了，但仍旧没有找到适应度值优于now_best的点，此时退而求其次选择当前邻域now_area中适应度值最优的点P的信息point_p赋值给now_point后返回到第③步。

3.禁忌搜索算法的一种改进

由于禁忌搜索算法十分依赖于初始值的选取，初始值选择不当会陷入到局部最优值当中，同时每次产生的邻域的大小是该算法搜索能力大小的体现，产生的邻域大，则搜索能力强。因此，一种改进的方式是动态的邻域生成规则，即当程序得到的目标最优值在迭代一定代数后仍旧不发生改变时，扩大邻域的大小，再次进行搜索可以跳出局部最优值。

二、matlab代码实现

1.问题描述

确定非线性二元函数f(x,y)的最大值

f(x,y)=3*(1-x)^2*exp(-x^2-(y+1)^2)-10*(x/5-x^3-y^5)*exp(-x^2-y^2)-1/3^exp(-(x+1)^2-y^2)

函数图像如下

可以看到该函数图像具有三个谷峰，对于没有改进的禁忌搜索算法，初值的选取十分重要，一旦选取不当便会陷入到局部最优值。

2.matlab代码

算法中将该函数用作适应度函数。如下是算法中的目标函数代码及画图用到的目标函数代码

①算法中的目标函数

function h=f(x,y)
h=3*(1-x)^2*exp(-x^2-(y+1)^2)-10*(x/5-x^3-y^5)*exp(-x^2-y^2)-1/3^exp(-(x+1)^2-y^2);
end

②画图用到的目标函数

function h=f1(x,y)
h=3*(1-x).^2.*exp(-x.^2-(y+1).^2)-10*(x./5-x.^3-y.^5).*exp(-x.^2-y.^2)-1./3.^exp(-(x+1).^2-y.^2);
end

③邻域生成函数

function s=producelinyu(x,y,h,SL)
s=cell(SL,SL);
m=(SL+1)/2;
for i=1:SL
    for j=1:SL
        s{i,j}=[x+h*(i-m),y+h*(j-m)];
    end
end
s{m,m}=inf;
end

④邻域搜索选点，生成点序列函数

function n=selectpoint(linyu,tab,SL,N,L,xlim,ylim)
for i=1:SL
    for j=1:SL
        if linyu{i,j}(1)==inf||linyu{i,j}(1)<xlim(1)||linyu{i,j}(1)>xlim(2)||linyu{i,j}(2)<ylim(1)||linyu{i,j}(2)>ylim(2)
            linyu{i,j}=inf;
        end
    end
end
tp=1;tpmax=(SL+1)^2;
for i=1:N
    kp=randi(SL,1,2);
    while(linyu{kp(1),kp(2)}(1)==inf||panduan(linyu,kp(1),kp(2),tab,L)==0)
        kp=randi(SL,1,2);
        tp=tp+1;
        if tp>tpmax
            break;
        end
    end
    if tp>tpmax
        n{i}=inf;
        continue;
    end
    n(i)=linyu(kp(1),kp(2));
    linyu{kp(1),kp(2)}=inf;
    tp=1;
end
end

⑤禁忌表函数

function tab=tabutable(L,n,tab1)
tab=tab1;
N=length(n);
for i=1:N
    if n{i}==inf
        continue;
    end
    for j=1:L
        if tab{j,2}==0||tab{j,3}==L+1
            kj=j;
            break;
        end
    end
    tab(kj,1)=n(i);tab{kj,2}=1;
    for j=1:L
        if tab{j,2}==1;
            tab{j,3}=tab{j,3}+1;
        end
        if tab{j,3}==L+2
            tab{j,3}=2;
        end
    end
end
    

⑥在邻域中搜某一点时，判断是否在禁忌表函数

function op=panduan(linyu,ki,kj,tab,L)
num=0;
for i=1:L
    if tab{i,1}==inf
        num=num+0;
    else
        num=num+1;
    end
end

if num==0
    op=1;
else
    for i=1:L
        if tab{i,1}(1)==inf
            continue;
        end
        if linyu{ki,kj}(1)==tab{i,1}(1)&&linyu{ki,kj}(2)==tab{i,1}(2)
            op=0;
            break;
        else
            op=1;
        end
    end
end
end

⑦判断某一邻域是否全部搜索完函数

function op=ifexitpoint(n)
num=0;
for i=1:length(n)
    if n{i}==inf
        num=num+0;
    else
        num=num+1;
    end
end
op=num;
end

⑧邻域动态变化函数

function [sl,edgepoint]=changelinyu(beststatus,SL,E)
num=1;k=1;
for i=1:length(beststatus)-1
    if abs(beststatus(i)-beststatus(i+1))<=E
        num=num+1;
        knum=num;
        ki(k)=i+1;
    else
        num=1;
        continue;
    end
    if knum>9
        k=k+1;
    end
end
if knum>9&&knum<14
    sl=SL+6;
else
    sl=SL;
end
for i=1:length(ki)
    edgepoint(i)=beststatus(ki(i));
end
end

⑨主函数

clear all;
clc;
%设置禁忌搜索算法的初始参数
h=0.1;%步长，代表整个禁忌搜索的整个邻域的精度
SL=5;%搜索邻域的长度(本算法中邻域生成方法必须保证SL为奇数)
N=10;%禁忌候选解的个数
L=500;%禁忌表的长度
K=100;%迭代次数
E=1e-6;
%所求函数的理论图像
x=-3:0.1:3;xlim=[x(1),x(length(x))];
y=-3:0.1:3;ylim=[y(1),y(length(y))];
xL=max(x)-min(x);
yL=max(y)-min(y);
[X,Y]=meshgrid(x,y);
ff=f1(X,Y);
figure(1);
surf(X,Y,ff);

%初始化禁忌表
tabu=cell(L,3);
for i=1:L
    tabu{i,1}=inf;
    tabu{i,2}=0;
    tabu{i,3}=1;
end
tab=tabu;
%设置搜索初始点
initial_x=0.4;initial_y=0.4;

best_x=initial_x;
best_y=initial_y;
beststatus(1)=f(best_x,best_y);
now_x=initial_x;
now_y=initial_y;
for k=1:K
    linyu=producelinyu(now_x,now_y,h,SL);
    n=selectpoint(linyu,tab,SL,N,L,xlim,ylim);
    op=ifexitpoint(n);
    kstart=0;kop=1;
    %clear shiyin;
    partbest=-inf;
    shiyin=-inf*ones(floor(SL^2/N)+1,N,k);
    while(op~=0)
        for i=1:N
            if n{i}(1)==inf
                continue;
            end
            shiyin(kop,i,k)=f(n{i}(1),n{i}(2));
        end
        [next_max,tp]=max(shiyin(kop,:,k));
        next_x=n{tp}(1);next_y=n{tp}(2);
        if next_max>partbest
            partbest=next_max;
            partbest_x=next_x;
            partbest_y=next_y;
        end
        if partbest>=beststatus(k)
            beststatus(k+1)=next_max;
            best_x=next_x;
            best_y=next_y;
            now_x=next_x;
            now_y=next_y;
            n{tp}=inf;
            tab=tabutable(L,n,tab);
            kstart=1;
            break;
        else
            kop=kop+1;
            tab=tabutable(L,n,tab);
            n=selectpoint(linyu,tab,SL,N,L,xlim,ylim);
            op=ifexitpoint(n);
            if kop>=floor(SL^2/N)+1
                break;
            end
        end
    end
    if kstart==1
        continue;
    else
        now_x=partbest_x;
        now_y=partbest_y;
        beststatus(k+1)=beststatus(k);
    end
    [SL,edgepoint{k,1}]=changelinyu(beststatus,SL,E);
end
figure(2);
plot(1:K+1,beststatus);

3.结果显示

主函数中倒数第四行[SL,edgepoint{k,1}]=changelinyu(beststatus,SL,E);为将禁忌搜索算法改进的重要算法。将其注释掉，即变成未改进的禁忌搜索算法。

所有方法都迭代100代，选取下图的四个点分别运用两种禁忌算法计算，四个点的分布图如下。

上图是四点分布周测图，下图为俯视图（等高图）

①未改进的禁忌搜索算法

选取初始点（3,3）

选取初始点（2，-1）

选取初始点（-1.4，-1.3）

选取（0.4,0.4）

注：这里也可能收敛到最优值7.14但是几率很小，大概50次中有一两次。

③改进的禁忌搜索算法

选取点（3,3）

选取点（2，-1）

选取点（-1.4，-1.3）

选取点（0.4,0.4）

三、结论

显而易见，通过上面的对比，可以很清晰的看到未改进的禁忌搜索算法极易陷入到局部最优值当中，而改进的禁忌搜索法可以跳出局部最优值而找到全局最优值（但有时也无法跳出局部最优值，此时再次重新运行1次主程序即可）。

值的注意的是：虽然未改进的禁忌搜索算法依赖于初始值的选取并且极易陷入到局部最优值，但是花费的时间很少；然而改进的禁忌搜索算法虽然可以跳出局部最优值，但是程序花费的时间较长，此问题改进的禁忌算法运行一次大概2分钟左右。这仅仅是一个比较简单的目标优化，若是一个十分复杂的问题，那时间肯定呈指数倍的增长。并且相同的问题运用遗传算法进行编程，很快就能够算出全局最优解。

因此，禁忌搜索法（局部搜索算法）一般不用做全局搜索最优解，而是作为全局搜索算法（如遗传算法）的后续操作，即全局搜索算法得到最优值作为禁忌搜索法的初始值而在局部更加精确的搜索全局最优解。