aco_knight_3.m

三维广义骑士巡游问题的蚁群算法的程序

%%=========================================================================
%%  Ant Colony Algorithm(ACO) For Search the Generial Knight's Tour 
                      %By Delei Jiang.2007.9.28
%%  Ant Colony Algorithm for Numbers of the knight's tour
%%-------------------------------------------------------------------------
%%  主要符号说明
%%  C        广义棋盘上的格子的坐标矩阵，n×3的矩阵
%%  NC_max   最大迭代次数
%%  S_max    合法的骑士巡游矩阵的最大选取次数
%%  m        蚂蚁个数
%%  Alpha    表征信息素重要程度的参数
%%  Beta     表征启发式因子重要程度的参数
%%  Rho      信息素蒸发系数
%%  Q        信息素增加强度系数
%%  R_best   各代最佳路线
%%  L_best   各代最佳路线的长度
%%=========================================================================

clear;clc;
%设置初始参数如下：
m=300;Alpha=2;Beta=5;Rho=0.25;NC_max=10;Q=1.0;
% C=[1 1; 1 2;1 3; 1 4; 1 5;2 1;2 2 ;2 3;2 4; 2 5;3 1; 3 2;... 
% 3 3; 3 4; 3 5;4 1; 4 2;4 3; 4 4;4 5;5 1;5 2 ;5 3; 5 4;5 5;];
row=8;%input('Please input the Width of the chessboard:');
col=8;%input('Please input the Height of the chessboard:');
dim=3;%input('Please input the dimenion of the chessboard:');
% x0=input('Please input intial postion of x:');
% y0=input('Please input intial postion of y:');
% z0=input('Please input intial postion of z:');
%%*************************模块1S***************************%%

C=zeros(dim*row*col,3);%广义骑士巡游棋盘的各个位置的坐标
count=0;
for k=1:dim
    for i=1:row
        for j=1:col
            count=count+1;
            C(count,1)=k;  C(count,2)=i;C(count,3)=j;
        end
    end
end
   
%%*************************模块1E***************************%%

Knight=[1 2 2];%以[1 2 2]广义马为例进行说明
%%第一步：变量初始化
n=size(C,1);%n表示问题的规模（广义棋盘的格数）
 
%%*************************模块2S***************************%%
D=zeros(n,n);%D表示完全图的赋权邻接矩阵
for i=1:n
    for j=1:n
        if (sum((C(i,:)-C(j,:)).^2)==sum(Knight.^2))...
                &(C(i,1)~=C(j,1))&(C(i,2)~=C(j,2))&(C(i,3)~=C(j,3))
            D(i,j)=eps;
        else
            D(i,j)=10e6;
        end
        D(j,i)=D(i,j);
    end
end
%*************************模块2E***************************%%

Eta=1./D;%Eta为启发因子，这里设为距离的倒数
Tau=10.^(-1).*ones(n,n);%Tau为信息素矩阵
Tabu=zeros(m,n);%存储并记录路径的生成
NC=1;%迭代计数器
R_best=zeros(NC_max,n);%各代最佳路线
L_best=inf.*ones(NC_max,1);%各代最佳路线的长度
L_ave=zeros(NC_max,1);%各代路线的平均长度
S=1;S_max=2;%作为计数器来选取合法的骑士巡游矩阵

%%*************************模块3S***************************%%
while S<S_max
 while NC<=NC_max  %停止条件之一：达到最大迭代次数
    %%第二步：将m只蚂蚁放到n个棋盘格子上
    Randpos=[];
    for i=1:(ceil(m/n))
        %ZB=(x0-1)*col*row+(y0-1)*col+z0; %初始化的坐标
        Randpos=[Randpos,randperm(n)];    %蚂蚁所放置的位置(指随机的多个位置)
        %Randpos=[Randpos,ZB.*ones(1,n)]; %蚂蚁所放置的位置(指人为的同一位置)
    end
    Tabu(:,1)=(Randpos(1,1:m))';
   
    %%第三步：m只蚂蚁按概率函数选择下一个格子，完成各自的周游
    for j=2:n
        for i=1:m
            visited=Tabu(i,1:(j-1));%已访问的格子
            J=zeros(1,(n-j+1));%待访问的格子
            P=J;%待访问格子的选择概率分布
            Jc=1;
            for k=1:n
                if length(find(visited==k))==0
                    J(Jc)=k;
                    Jc=Jc+1;
                end
            end
            %下面计算待选格子的概率分布
            for k=1:length(J)
                P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);
            end
            P=P/(sum(P));
            %按概率原则选取下一个格子
            Pcum=cumsum(P);
            Select=find(Pcum>=rand);
            to_visit=J(Select(1));
            Tabu(i,j)=to_visit;
        end
    end
    if NC>=2
        Tabu(1,:)=R_best(NC-1,:);
    end
   
    %%第四步：记录本次迭代最佳路线
    L=zeros(m,1);
    for i=1:m
        R=Tabu(i,:);
        for j=1:(n-1)
            L(i)=L(i)+D(R(j),R(j+1));
        end
        L(i)=L(i)+D(R(1),R(n));
    end
    L_best(NC)=min(L);
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:);
    L_ave(NC)=mean(L);
    NC=NC+1;
   
    %%第五步：更新信息素  (关键点)
    Delta_Tau=zeros(n,n);
    for i=1:m
        for j=1:(n-1)
            %[r move]=find(Tabu(i,:)==j);  %<1.1>
            %move=sum(Tabu(i,1:j));  %<1.2>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q*(move-j)/(n-j);  %<1>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q*(visited(j)-j)/(n-j);    %<2>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q*(visited(end)-j)/(n-j);  %<3>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/D(j,j+1);  %<4>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q*(L(i)-j)/(n-j);  %<5>
            Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);  %<6>
            %Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q*(Tabu(i,j)-j)/(n-j);  %<7>
        end
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q*(move-n)/(n-j);   %<1>
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q*(visited(j)-n)/(n-j);  %<2>
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q*(visited(end)-n)/(n-j);   %<3>
        %(没有这个的话就是闭路，否则就是开路)
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/D(j,j+1);   %<4>
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q*(L(i)-n)/(n-j)  %<5>
        Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i); %<6>
        %Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q*(Tabu(i,n)-n)/(n-j);   %<7>
    end
    Tau=(1-Rho).*Tau+Delta_Tau;

    %%第六步：禁忌表清零
    Tabu=zeros(m,n);
 end
%%*************************模块3E***************************%%

%%*************************模块4S***************************%%
%%第七步：输出结果
 Pos=find(L_best==min(L_best));
 Shortest_Route=R_best(Pos(1),:);
 Shortest_Length=L_best(Pos(1));
%%*************************模块4E***************************%%


%%*************************模块5S***************************%%
 for i=1:n
    U(i)=ceil(Shortest_Route(i)/col/row);
    V(i)=ceil((Shortest_Route(i)-(U(i)-1)*col*row)/col);
    W(i)=mod((Shortest_Route(i)-(U(i)-1)*col*row),col);
 end
 V(V==0)=1;     %修正坐标V
 W(W==0)=col;   %修正坐标W
%%*************************模块5E***************************%%

if sum((diff(W).^2+diff(V).^2+diff(U).^2)==sum(Knight.^2).*ones(1,n-1))==n-1
    break;
 else
     S=S+1;
     %if S>=S_max
         %S=S_max-1;
     %end
end
 
end

%%*************************模块6S***************************%%
 chessboard=zeros(row,col,dim);
 for i=1:n
    chessboard(V(i),W(i),U(i))=i;  %产生骑士巡游矩阵
 end
%%*************************模块6E***************************%%

if S==S_max
    figure(1);
    axis([0.9 dim 0.9 row 0.9 col])
    text((1+dim)/2,(1+row)/2,(1+col)/2,'此次没有找到合法的骑士巡游路径','FontSize',18,'HorizontalAlignment','center' );
else
%%*************************模块7S***************************%%
 figure(1);
 plot(L_best,'r')
 hold on;
 plot(L_ave,'k.-')

 h = legend('各代最佳路线长度','各代路线平均长度',2);
 set(h,'Interpreter','none')

 figure(2);
 scatter3(C(:,1),C(:,2),C(:,3));
 hold on;
% whitebg('w')

for i=1:n
    if mod(i,2)==0
        plot3([U(i) U(i)],[V(i) V(i)],[W(i) W(i)],'o-','LineWidth',1,'MarkerFaceColor','c','MarkerSize',20);
    else
        plot3([U(i) U(i)],[V(i) V(i)],[W(i) W(i)],'o-','LineWidth',1,'MarkerFaceColor','r','MarkerSize',20);
    end
    text(U(i),V(i),W(i),num2str(i),'FontSize',15,'HorizontalAlignment','center' );
    hold on;
end

for i=1:n-1
    if mod(i,2)==0
        plot3([U(i) U(i+1)],[V(i) V(i+1)],[W(i) W(i+1)]);
    else
        plot3([U(i) U(i+1)],[V(i) V(i+1)],[W(i) W(i+1)]);
    end
    hold on;
end

if  (W(n)-W(1)).^2+(V(n)-V(1)).^2+(U(n)-U(1)).^2==sum(Knight.^2) %~(mod(row,2))|~(mod(col,2))
    plot([U(n) U(1)],[V(n) V(1)],[W(n) W(1)]);
end

 axis([0.9 dim 0.9 row 0.9 col])
%%*************************模块7E***************************%%
end