禁忌搜索算法TS求解TSP问题

目录

一、局部邻域搜索

二、禁忌搜索

三、禁忌搜索算法流程

四、算法求解例题 

---

一、局部邻域搜索

局部邻域搜索是基于贪婪准则持续地在当前的邻域中进行搜索，虽然算法通用，易于实现，且容易理解，但其搜索性能完全依赖于邻域结构和初始解，尤其容易陷入局部极小值无法保证全局优化

算法可以描述为：

1）选定一个初始可行解：；记录当前最优解，，其中表示的邻域。

2）当（空集），或满足其他停止运算准则是，输出计算结果，停止运算，否则，继续步骤3）

3）从中选一个集合，得到中的最好解。若，则，；否则，，重复步骤2），继续搜索

这种邻域搜索的方法易于理解，易于实现，而且具有很好的通用性，但是搜索结果的好坏完全依赖于初始解和邻域的结构。

二、禁忌搜索

对于一个初始解，在一种邻域范围内对其进行一系列变化，从而得到许多候选解。从这些候选解中选出最优候选解，将候选解对应的目标值于best-so-far状态进行比较。若其目标值优于best-so-far状态，就将该候选解解禁，用来替代当前最优解及其best-so-far状态，然后将其加入禁忌表，再将禁忌表中相应对象的禁忌长度改变；如果候选解的目标值都不优于best-so-far，就从候选解中选出不属于禁忌对象的最佳状态，并将其作为新的当前解，不用与当前解比较，直接将其所对应的对象作为禁忌对象，并将禁忌表中相应对象的禁忌长度进行修改。

三、禁忌搜索算法流程

禁忌搜索算法基本思想是：给定一个当前解（初始解）和一种邻域，然后在当前解的邻域中确定若干候选解；若候选解对应的目标值优于best-so-far状态，并将相应的对象加入禁忌表，同时修改禁忌表中各对象的任期，若不存在上述候选解，则在候选解中选择非禁忌的最佳状态为新的当前解，而无视它与当前解的优劣，同时将相应的对象加入禁忌表，并修改禁忌表中各对象人气。如此重复，直到满足停止准则。算法步骤可描述如下：

1）给定禁忌搜索算法参数，随机产生初始解，置禁忌表为空。

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

3）利用当前解的邻域函数产生其所有邻域解，并从中确定若干候选解

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

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

6）判断算法终止条件是否满足

 

四、算法求解例题 

旅行商问题（TSP问题）。假设有一个旅行商人要拜访全国31个省会城市，他需要选择所要走的路径，路径的限制是每个城市只能拜访一次，二球要最后回到原来出发的城市。路径的选择要求是：所选的路径的路程之和中的最小。

    全国31个省会的坐标为[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;4196 1044;4312  790;4386  570;3007 1970;2562 1756;2788 1491;2381 1676;1332  695;3715 1678;3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]

解:仿真过程如下:

(1)初始化优化城市规模N= 31,，禁忌长度TbuL=22,候选集的个数Ca=200，最大迭代次数G=1000。

(2)计算任意两个城市的距离间隔矩阵D;随机产生组路径为初始解 S%,计算其适配值，并将其赋给当前最佳解bestsofar.

(3)定义初始解的邻域映射为2-opt形式，即初始解路径中的两个城市坐标进行对换。产生Ca个候选解，计算候选解的适配值，并保留前Ca/2个最好候选解。

(4)对候选解判断是否满足藐视准则:若满足，则用满足藐视准则的解替代初始解成为新的当前最佳解，并更新禁忌表Tabu和禁忌长度TabuL,然后转步骤(6); 否则，继续以下步骤。

(5)判断候选解对应的各对象的禁忌属性，选择候选解集中非禁忌对象所对应的最佳状态为新的当前解，同时更新禁忌表Tabu和禁忌长度TabuL。

(6)判断是否满足终止条件:若满足，则结束搜索过程，输出优化值:若不满足，则继续进行迭代优化。

优化后的路径如图所示

%%%%%%%%%%%%%%%%%%%%%禁忌搜索算法解决TSP问题%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                        %清除所有变量
close all;                        %清图 
clc;                              %清屏
C = [1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
    3238 1229;4196 1044;4312  790;4386  570;3007 1970;2562 1756;...
    2788 1491;2381 1676;1332  695;3715 1678;3918 2179;4061 2370;...
    3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
    3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
    2370 2975];                   %31个省会城市坐标
N=size(C,1);                      %TSP问题的规模,即城市数目
D=zeros(N);                       %任意两个城市距离间隔矩阵
%%%%%%%%%%%%%%%%%%%%%求任意两个城市距离间隔矩阵%%%%%%%%%%%%%%%%%%%%%
for i=1:N
    for j=1:N
        D(i,j)=((C(i,1)-C(j,1))^2+...
            (C(i,2)-C(j,2))^2)^0.5;
    end
end
Tabu=zeros(N);                      %禁忌表
TabuL=round((N*(N-1)/2)^0.5);       %禁忌长度
Ca=200;                             %候选集的个数(全部领域解个数)
CaNum=zeros(Ca,N);                  %候选解集合
S0=randperm(N);                     %随机产生初始解
bestsofar=S0;                       %当前最佳解
BestL=Inf;                          %当前最佳解距离
figure(1);
p=1;
Gmax=300;                          %最大迭代次数   
%%%%%%%%%%%%%%%%%%%%%%%%%%%禁忌搜索循环%%%%%%%%%%%%%%%%%%%%%%%%%%
while p<Gmax
    ALong(p)=func1(D,S0);            %当前解适配值
    %%%%%%%%%%%%%%%%%%%%%%%%%%%交换城市%%%%%%%%%%%%%%%%%%%%%%%%%%
    i=1;
    A=zeros(Ca,2);                   %解中交换的城市矩阵
    %%%%%%%%%%%%%%%%%求领域解中交换的城市矩阵%%%%%%%%%%%%%%%%%%%%%
    while i<=Ca
        M=N*rand(1,2);
        M=ceil(M);         
        if M(1)~=M(2)
            A(i,1)=max(M(1),M(2));
            A(i,2)=min(M(1),M(2));
            if i==1
                isa=0;
            else
                for j=1:i-1
                    if A(i,1)==A(j,1) && A(i,2)==A(j,2)
                        isa=1;
                        break;
                    else
                        isa=0;
                    end
                end
            end
            if ~isa
                i=i+1;
            else
            end
        else
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%产生领域解%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%保留前BestCaNum个最好候选解%%%%%%%%%%%%%%%%%%%
    BestCaNum=Ca/2;
    BestCa=Inf*ones(BestCaNum,4);
    F=zeros(1,Ca);
    for i=1:Ca
        CaNum(i,:)=S0;
        CaNum(i,[A(i,2),A(i,1)])=S0([A(i,1),A(i,2)]);
        F(i)=func1(D,CaNum(i,:));
        if i<=BestCaNum
            BestCa(i,2)=F(i);
            BestCa(i,1)=i;
            BestCa(i,3)=S0(A(i,1));
            BestCa(i,4)=S0(A(i,2));
        else
            for j=1:BestCaNum
                if F(i)<BestCa(j,2)
                    BestCa(j,2)=F(i);
                    BestCa(j,1)=i;
                    BestCa(j,3)=S0(A(i,1));
                    BestCa(j,4)=S0(A(i,2));
                    break;
                end
            end
        end
    end
    [JL,Index]=sort(BestCa(:,2));
    SBest=BestCa(Index,:);
    BestCa=SBest;
    %%%%%%%%%%%%%%%%%%%%%%%%藐视准则%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if BestCa(1,2)<BestL
        BestL=BestCa(1,2);
        S0=CaNum(BestCa(1,1),:);
        bestsofar=S0;
        for m=1:N
            for n=1:N
                if Tabu(m,n)~=0
                    Tabu(m,n)=Tabu(m,n)-1;    
                    %更新禁忌表
                end
            end
        end
        Tabu(BestCa(1,3),BestCa(1,4))=TabuL;
        %更新禁忌表
    else
        for  i=1:BestCaNum
            if  Tabu(BestCa(i,3),BestCa(i,4))==0
                S0=CaNum(BestCa(i,1),:);
                for m=1:N
                    for n=1:N
                        if Tabu(m,n)~=0
                            Tabu(m,n)=Tabu(m,n)-1;
                            %更新禁忌表
                        end
                    end
                end
                Tabu(BestCa(i,3),BestCa(i,4))=TabuL;
                %更新禁忌表
                break;
            end
        end
    end
    ArrBestL(p)=BestL;
    p=p+1;   
    for i=1:N-1
        plot([C(bestsofar(i),1),C(bestsofar(i+1),1)],...
            [C(bestsofar(i),2),C(bestsofar(i+1),2)],'bo-');
        hold on;
    end
    plot([C(bestsofar(N),1),C(bestsofar(1),1)],...
        [C(bestsofar(N),2),C(bestsofar(1),2)],'ro-');
    title(['迭代次数：',num2str(p),'   优化最短距离:',num2str(BestL)]);
    hold off;
    pause(0.005);
end
BestShortcut=bestsofar;            %最佳路线
theMinDistance=BestL;              %最佳路线长度
figure(2);
plot(ArrBestL); 
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
%%%%%%%%%%%%%%%%%%%%%%%%%适配值函数%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=func1(D,s)
DistanV=0;
n=size(s,2);
for i=1:(n-1)
    DistanV=DistanV+D(s(i),s(i+1));
end
DistanV=DistanV+D(s(n),s(1));
F=DistanV;
end