智能优化算法及其MATLAB实例（第二版）——蚁群算法5.1-5.2

例5.1 旅行商问题(TSP问题)。假设有一个旅行商人要拜访全国31个省会城市，他需要选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择要求是:所选路径的路程为所有路径之中的最小值。
全国31个省会城市的坐标为「1304 2312:3639 1315;41772244: 3712 1399:3488 1535: 3326 1556; 3238 1229:4196 1004;4312 790:4386 570:3007 1970:2562 1756:2788 1491;23811676: 1332 695: 3715 1678: 3918 2179:4061 2370:3780 2212:3676 2578: 4029 2838: 4263 2931:3429 1908: 3507 2367: 33942643: 3439 3201; 2935 3240:3140 3550: 2545 2357:2778 2826:2370 2975]。

解:仿真过程如下:
(1)初始化蚂蚁个数m50，信息素重要程度参数Alpha=1，启发式大子重要程度参数Beta=5，信息素蒸发系数Rho=0.1，最大迭代次数 G200，信息素增加强度系数Q100。
(2)将m个蚂蚁置于n个城市上，计算待选城市的概率分布，m只蚂蚁按概率函数选择下一座城市，完成各自的周游。
(3)记录本次迭代最佳路线，更新信息素，禁忌表清零。
(4)判断是否满足终止条件:若满足，则结束搜索过程，输出优化值;若不满足，则继续进行迭代优化。
优化后的路径和适应度进化曲线如结果图所示。 

m=50;%%%蚂蚁个数
Alpha=1;%%%信息素重要程度参数
Beta=5;%%%%启发式因子重要程度参数
Rho=0.1;%%%信息素蒸发系数
G=200;%%%%最大迭代次数
Q=100;%%信息素增加强度系数
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;3240 3550;2545 2357;2778 2826;...
    2370 2975];
%%%%第一步，变量初始化%%%
n=size(C,1);%%%n表示城市的个数
D=zeros(n,n);%%%D用来储存两个城市的间隔
for i =1:n
    for j =1:n
        if i~=j
            D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;%%%%C（i,1）表示第一列存的是x坐标，C(i,2)表示第二列存的是y坐标
        else
            D(i,j)=eps;%%%i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps表示
        end
        D(j,i)=D(i,j);%%%对称矩阵
    end
end
Eta=1./D;%%%%Eta为启发因子，这里设为距离的倒数
Tau=ones(n,n);%%%Tau为信息素矩阵
Tabu=zeros(m,n);%%%存储并记录路径的生成，表示边(i,j)上的信息素量，禁忌表
NC=1;%迭代计数器
R_best=zeros(G,n);%%%各代最佳路线
L_best=inf.*ones(G,1);%%%各代最佳路线的长度
figure(1);%%%优化解
while NC<=G
%%%%%%第二步，将m只蚂蚁放在n个城市上，此时是m只蚂蚁的初始位置%%%%%
Randpos=[];%%%随即存取,定义随机位置
for i =1:(ceil(m/n))
    Randpos=[Randpos,randperm(n)];%%%%randperm整数的随机排列
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 %%%开始时置0
                J(Jc)=k;
                Jc=Jc+1; %%%访问的城市个数自加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); %%%开始距离为0，m*1的矩阵
for i=1:m
    R=Tabu(i,:);
    for j=1:(n-1)
        L(i)=L(i)+D(R(j),R(j+1)); %%%原距离加上第j个城市到第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),:);%%%%此次迭代后的最佳路线
%%%%第五步，更新信息素
Delta_Tau=zeros(n,n); %%%开始时，信息素为n*n的0矩阵
for i =1:m
    for j=1:(n-1)
        Delta_Tau(Tabu(i,j),Tabu(i,j+1))=...
             Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);%%此次循环在路径(i,j)上的信息素增量
    end
    Delta_Tau(Tabu(i,n),Tabu(i,1))=...
             Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);%%此次循环在整个路径上的信息素增量
end
Tau=(1-Rho).*Tau+Delta_Tau; %%%考虑信息素挥发后，更新后的信息素
%%%%第六步，禁忌表清零%%%
Tabu=zeros(m,n);  %%%直到最大迭代次数
%%%%历代最优路线%%%%%%%%
for i =1:n-1
    plot([C(R_best(NC,i),1),C(R_best(NC,i+1),1)],...
        [C(R_best(NC,i),2),C(R_best(NC,i+1),2)],'bo-');
    hold on;
end
plot([C(R_best(NC,n),1),C(R_best(NC,1),1)],...
        [C(R_best(NC,n),2),C(R_best(NC,1),2)],'ro-');
title(['优化最短距离：',num2str(L_best(NC))]);
hold off
pause(0.005); %%%暂停执行0.005秒
NC=NC+1;
end
%%%第七步，输出结果%%%
Pos=find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:);%%%%最佳路线
Shortest_length=L_best(Pos(1));%%%最佳路线长度
figure(2),
plot(L_best)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
            

 运行结果

 

 

 

解:仿真过程如下:
(1)初始化蚂蚁个数mF20，最大迭代次数G-200，信息素蒸发系数R=0.9，转移概率常数P=0.2，局部搜索步长step=0.1。
(2)随机产生蚂蚁初始位置，计算适应度承数值，设为初始信息素，计算状态转移概率。 

(3)进行位置更新:当状态转移概率小于转移概率常数时，进行局部搜索;当状态转移概率大于转移概率常数时，进行全局搜索，产生新的蚂蚁位置，并利用边界吸收方式进行边界条件处理，将蚂蚁位置界定在取值范围内。
(4)计算新的蚂蚁位置的适应度值，判断蚂蚁是否移动，更新信息素。
(5)判断是否满足终止条件:若满足，则结克搜索过程，输出优化值;若不满足，则继续进行迭代优化。
适应度值进化曲线如图所示 

%%%%%%%%%%%%%%%%%%%%蚁群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;               %清图
clc;                     %清屏
m=20;                    %蚂蚁个数
G_max=200;               %最大迭代次数
Rho=0.9;                 %信息素蒸发系数
P0=0.2;                  %转移概率常数
XMAX= 5;                 %搜索变量x最大值
XMIN= -5;                %搜索变量x最小值
YMAX= 5;                 %搜索变量y最大值
YMIN= -5;                %搜索变量y最小值
%%%%%%%%%%%%%%%%%随机设置蚂蚁初始位置%%%%%%%%%%%%%%%%%%%%%%
for i=1:m
    X(i,1)=(XMIN+(XMAX-XMIN)*rand);
    X(i,2)=(YMIN+(YMAX-YMIN)*rand);
    Tau(i)=func(X(i,1),X(i,2));
end
step=0.1;                %局部搜索步长
for NC=1:G_max
    lamda=1/NC;
    [Tau_best,BestIndex]=min(Tau);
    %%%%%%%%%%%%%%%%%%计算状态转移概率%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        P(NC,i)=(Tau(BestIndex)-Tau(i))/Tau(BestIndex);
    end
    %%%%%%%%%%%%%%%%%%%%%%位置更新%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        %%%%%%%%%%%%%%%%%局部搜索%%%%%%%%%%%%%%%%%%%%%%
        if P(NC,i)<P0
            temp1=X(i,1)+(2*rand-1)*step*lamda;
            temp2=X(i,2)+(2*rand-1)*step*lamda;
        else
            %%%%%%%%%%%%%%%%全局搜索%%%%%%%%%%%%%%%%%%%%%%%
            temp1=X(i,1)+(XMAX-XMIN)*(rand-0.5);
            temp2=X(i,2)+(YMAX-YMIN)*(rand-0.5);
        end
        %%%%%%%%%%%%%%%%%%%%%边界处理%%%%%%%%%%%%%%%%%%%%%%%
        if temp1<XMIN
            temp1=XMIN;
        end
        if temp1>XMAX
            temp1=XMAX;
        end
        if temp2<YMIN
            temp2=YMIN;
        end
        if temp2>YMAX
            temp2=YMAX;
        end
        %%%%%%%%%%%%%%%%%%蚂蚁判断是否移动%%%%%%%%%%%%%%%%%%
        if func(temp1,temp2)<func(X(i,1),X(i,2))
            X(i,1)=temp1;
            X(i,2)=temp2;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%更新信息素%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        Tau(i)=(1-Rho)*Tau(i)+func(X(i,1),X(i,2));
    end
    [value,index]=min(Tau);
    trace(NC)=func(X(index,1),X(index,2));
end
[min_value,min_index]=min(Tau);
minX=X(min_index,1)                           %最优变量
minY=X(min_index,2)                          %最优变量
minValue=func(X(min_index,1),X(min_index,2))  %最优值
figure
plot(trace)
xlabel('搜索次数');
ylabel('适应度值');
title('适应度进化曲线')

%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%
function value=func(x,y)
value=20* (x^2-y^2)^2-(1-y)^2-3* (1+y)^2+0.3;
end

 运行结果