禁忌搜索（Tabu Search or Taboo Search,TS)算法解决3DTSP问题

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

clc;
clear all;
file_path = '..\chapter2\data.csv'; %读取文件
Location = csvread(file_path, 1, 1);
N=size(Location,1);                    %TSP问题的规模
D=zeros(N);                     %任意两个加工区域的距离间隔矩阵
%%%%%%%%%%%%%%%%%%%%%求任意两个城市距离间隔矩阵%%%%%%%%%%%%%%%%%%%%%
for i=1:N
    for j=1:N
        D(i,j)=((Location(i,1)-Location(j,1))^2+(Location(i,2)-Location(j,2))^2+(Location(i,3)-Location(j,3)))^0.5;
    end
end
Tabu=zeros(N);                      %禁忌表
TabuL=round((N*(N-1)/2)^0.5);       %禁忌长度
Ca=400;                             %候选集的个数(全部邻域解个数)
CaNum=zeros(Ca,N);                  %候选解集合
%S0=[22	78	122	101	49	131	80	21	100	97	103	99	96	37	12	93	31	102	98	115	41	44	19	20	109	107	111	112	64	35	46	42	75	113	16	71	74	72	69	70	73	68	67	63	29	76	120	50	9	105	62	130	126	124	128	129	125	127	123	84	34	33	66	60	116	59	57	53	54	58	56	52	13	94	25	121	117	26	11	132	85	87	91	88	89	92	86	90	28	15	36	38	17	48	134	39	2	114	61	135	82	95	10	65	77	1	5	3	8	7	4	6	55	79	24	51	104	133	14	23	18	27	106	81	110	108	30	118	45	47	40	136	83	32	119	43];  %ACO算法的结果设置为初始解
S0=randperm(N);               %随机产生初始解
bestsofar=S0;                       %当前最佳解
BestL=Inf;                          %当前最佳解距离
figure(1);
p=1;
Gmax=2000;                          %最大迭代次数   
%%%%%%%%%%%%%%%%%%%%%%%%%%%禁忌搜索循环%%%%%%%%%%%%%%%%%%%%%%%%%%
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
        plot3([Location(bestsofar(i),1),Location(bestsofar(i+1),1)],...
            [Location(bestsofar(i),2),Location(bestsofar(i+1),2)], ...
            [Location(bestsofar(i),3),Location(bestsofar(i+1),3)],'bo-');
        hold on;
    end
    plot3([Location(bestsofar(N),1),Location(bestsofar(1),1)],...
        [Location(bestsofar(N),2),Location(bestsofar(1),2)], ...
        [Location(bestsofar(N),3),Location(bestsofar(1),3)],'ro-');
    title(['优化最短距离:',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;