Matlab实现蚁群算法求解旅行商优化问题（TSP）（理论+例子+程序）

一、蚁群算法

        蚁群算法由意大利学者Dorigo M等根据自然界蚂蚁觅食行为提岀。蚂蚁觅食行为表示大量蚂蚁组成的群体构成一个信息正反馈机制，在同一时间内路径越短蚂蚁分泌的信息就越多，蚂蚁选择该路径的概率就更大。

        蚁群算法的思想来源于自然界蚂蚁觅食，蚂蚁在寻找食物源时，会在路径上留下蚂蚁独有的路径标识——信息素，蚂蚁会感知其他蚂蚁在各条路径上留下的信息素，并根据各条路径上的信息素浓度来选择之后要走的路，路径上留有的信息浓度越高，则蚂蚁更倾向于选择该路径。在蚂蚁选择某条路径后也会在改路径上留下信息素吸引更多蚂蚁选择该路径，随着时间的推移，信息素浓度不断增大，蚂蚁选择路径的概率也随之增高，由此形成了正反馈机制。由于蚁群算法的正反馈性，因此蚁群算法也属于增强型学习算法的其中一种。

初始时刻，不妨将P kij (t)设为t时刻蚂蚁k从结点i转移到结点j的概率。“蚂蚁TSP”策略收到两方面的左右，首先是访问某结点的概率，这个概率的大小依赖于其他蚂蚁释放的信息素浓度。所以定义：

        式中，nkij(t)为启发函数，表示蚂蚁从结点i转移到结点j的概率；allowk为蚂蚁k下一步可转移结点的集合，随着时间的推移，allowk储存的元素数量会减小，最终会变为空集合。a 为信息素重要程度因子。

        与实际情况类似的一点是：随着时间的推移，残留在路径上的信息素会逐渐挥发，蚂蚁在经过路径时残留的信息素量也会逐渐等同于信息素挥发量，最终使信息素残留量趋于稳定。令α表示信息素挥发程度，那么所有蚂蚁遍历完所有结点之后，各路径上的信息素残留量的数学表达式如下：

        式中，ckij为第k只蚂蚁在连接结点i 与结点k的路径上释放信息素而增加的信息素浓度。Δckij为所有蚂蚁在结点i 与结点k 连接路径上释放信息素而增加的信息素浓度，通常情况下：

        式中，Q为路径信息素常量，I为第k 只蚂蚁所经过路径的总长度。

二、蚁群算法改进（自适应）

改进（自适应蚁群算法）：

1）每次循环结束后求出最优解将其保留。

2）自适应的改变 值。

1. 信息素挥发系数的存在，会让没有搜索到的信息素的量减小到接近于0，降低了算法的全局搜索能力。

2. 当过大，且解的信息量过大时，曾经搜索过的解被重新搜索的可能性会变大。

3. 减小提高算法的全局搜索能力，但这会让算法的收敛速度降低

实现方法：

1. 的初始值为

2. 当算法求得的最优解在N次循环内没有改进时，减小为

三、实现步骤

四、代码结果

改进前：

改进后：

五、改进后运行的数据

1. 迭代最后城市之间的信息素（部分）

2. 最佳路径记录（部分）

3. 城市初始化顺序

4. 最终城市顺序（最优结果）

5. 最终蚁群算法初始参数

sumnum=0;       %记录选择概率全为0的次数
m=100;         %蚂蚁个数           
Alpha=1;      %信息素重要程度参数                         
Beta=5;       %启发式因子重要程度参数            
Rho=1;        %信息素蒸发系数           
Rho_min=0.2;    %最小信息素蒸发系数   
num_G=0;         %迭代多少次最优值不变得次数
num_G_max=15;    %最大迭代多少次最优值不变得次数
G_max=200;    %最大迭代次数            
Q=100;        %信息素增加强度系数  

六、代码（改进前）

clear all;      %清除所有变量           
close all;      %清图        
clc;            %清屏%% 初始化
%蚂蚁个数：50
%信息素重要程度参数：1
%启发式因子重要程度参数：5
%信息素蒸发系数：0.1
%最大迭代次数：200
%信息素增加强度系数：100
m=100;         %蚂蚁个数           
Alpha=1;      %信息素重要程度参数                         
Beta=5;       %启发式因子重要程度参数            
Rho=0.15;      %信息素蒸发系数            
G_max=300;    %最大迭代次数            
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;3140 3550;2545 2357;2778 2826;
% 2370 2975];   %31 个省会城市坐标C=[6734	1453;2233	10;5530	1424;401	841;3082	1644;7608	4458;
7573	3716;7265	1268;6898	1885;1112	2049;5468	2606;5989	2873;
4706	2674;4612	2035;6347	2683;6107	669;7611	5184;7462	3590;
7732	4723;5900	3561;4483	3369;6101	1110;5199	2182;1633	2809
4307	2322;675	1006;7555	4819;7541	3981;3177	756;7352	4506;
7545	2801;3245	3305;6426	3173;4608	1198;23	2216;7248	3779;
7762	4595;7392	2244;3484	2829;6271	2135;4985	140;1916	1569;
7280	4899;7509	3239;10	2676;6807	2993;5185	3258;3023	1942];        %% 第一步：变量初始化n=size(C,1);          %n 表示问题的规模（城市个数）  
D=zeros(n,n);         %D 表示两个城市距离间隔矩阵  
for i=1:nfor j=1:nif i~=jD(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;    %计算两两城市之间的距离elseD(i,j)=eps;endD(j,i)=D(i,j);end
end
Eta=1./D;                       %Eta 为启发因子，这里设为距离的倒数
Tau=ones(n,n);                  %Tau 为信息素矩阵 
Tabu=zeros(m,n);                %存储并记录路径的生成
NC=1;                           %迭代计数器
R_best=zeros(G_max,n);          %各代最佳路线
L_best=inf.*ones(G_max,1);      %各代最佳路线的长度
figure(1);                      %优化解
%% 判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。
while NC<=G_max            
%% 第二步：将 m 只蚂蚁放到 n 个城市上Randpos=[];for i=1:(ceil(m/n))Randpos=[Randpos,randperm(n)];endTabu(:,1)=(Randpos(1,1:m))'; 
%% 第三步：m 只蚂蚁按概率函数选择下一座城市，完成各自的周游
%将 m 个蚂蚁置于n个城市上，计算待选城市的概率分布，m 只蚂蚁按概率函数选择下一座城市，完成各自的周游。for j=2:nfor i=1:mvisited=Tabu(i,1:(j-1));  %己访问的城市J=zeros(1,(n-j+1));       %待访问的城市P=J;                      %待访问城市的选择概率分布Jc=1;for k=1:nif length(find(visited==k))==0J(Jc)=k;Jc=Jc+1;endend
%计算待选城市的概率分布for k=1:length(J)P(k)=(Tau(visited(end),J(k))^Alpha)...*(Eta(visited(end),J(k))^Beta);endP=P/(sum(P));
%按概率原则选取下一个城市Pcum=cumsum(P);         % 如 P=[1 2 3 4]，则cumsum(P)=[1 3 6 10],要累加，轮盘赌法，依次看是否在转得的区域内Select=find(Pcum>=rand);%轮盘赌法随机选择to_visit=J(Select(1));  %待选择的城市Tabu(i,j)=to_visit;     %访问的城市endendif NC>=2Tabu(1,:)=R_best(NC-1,:);   %最佳路线end
%% 第四步：记录本次迭代最佳路线L=zeros(m,1);for i=1:mR=Tabu(i,:);  %第i只蚂蚁走过的城市for j=1:(n-1)L(i)=L(i)+D(R(j),R(j+1));  %计算i只蚂蚁走过的路径endL(i)=L(i)+D(R(1),R(n)); %加上初始位置的路径endL_best(NC)=min(L);      %获取路径最少的一只pos=find(L==L_best(NC));%在50只蚂蚁中寻找走最少的一只R_best(NC,:)=Tabu(pos(1),:);%记录最佳路径%% 第五步：更新信息素%蚁环算法更新信息素增量%离线更新——蚁群Delta_Tau=zeros(n,n);for i=1:mfor 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);%信息素增量的更新endDelta_Tau(Tabu(i,n),Tabu(i,1))=...Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);%加上起使信息素增量的更新endTau=(1-Rho).*Tau+Delta_Tau;  %更新公式
%% 第六步: 禁忌表清零Tabu=zeros(m,n);
%历代最优路线for i=1:n-1plot([ 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;endplot([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);NC=NC+1;   %迭代+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('适应度进化曲线')

七、代码（改进后）

clear all;      %清除所有变量           
close all;      %清图        
clc;            %清屏
%% 初始化sumnum=0;       %记录选择概率全为0的次数
m=100;         %蚂蚁个数           
Alpha=1;      %信息素重要程度参数                         
Beta=5;       %启发式因子重要程度参数            Rho=1;        %信息素蒸发系数           
Rho_min=0.2;    %最小信息素蒸发系数   
num_G=0;         %迭代多少次最优值不变得次数
num_G_max=3;    %最大迭代多少次最优值不变得次数G_max=200;    %最大迭代次数            
Q=100;        %信息素增加强度系数            %48 个省会城市坐标
C=[6734	1453;2233	10;5530	1424;401	841;3082	1644;7608	4458;
7573	3716;7265	1268;6898	1885;1112	2049;5468	2606;5989	2873;
4706	2674;4612	2035;6347	2683;6107	669;7611	5184;7462	3590;
7732	4723;5900	3561;4483	3369;6101	1110;5199	2182;1633	2809
4307	2322;675	1006;7555	4819;7541	3981;3177	756;7352	4506;
7545	2801;3245	3305;6426	3173;4608	1198;23	2216;7248	3779;
7762	4595;7392	2244;3484	2829;6271	2135;4985	140;1916	1569;
7280	4899;7509	3239;10	2676;6807	2993;5185	3258;3023	1942];        %% 第一步：变量初始化
aaaa=eps;
n=size(C,1);          %n 表示问题的规模（城市个数）  
D=zeros(n,n);         %D 表示两个城市距离间隔矩阵  
for i=1:nfor j=1:nif i~=jD(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;    %计算两两城市之间的距离elseD(i,j)=eps;%eps表示从 1.0 到下一个最大双精度数的距离endD(j,i)=D(i,j);end
end
Eta=1./D;                       %Eta 为启发因子，这里设为距离的倒数
Tau=ones(n,n);                  %Tau 为信息素矩阵 
Tabu=zeros(m,n);                %存储并记录路径的生成
NC=1;                           %迭代计数器
R_best=zeros(G_max,n);          %各代最佳路线
L_best=inf.*ones(G_max,1);      %各代最佳路线的长度
figure(1);                      %优化解
%% 判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。
while NC<=G_max            
%% 第二步：将 m 只蚂蚁放到 n 个城市上Randpos=[];for i=1:(ceil(m/n))%蚂蚁个数除以城市个数向上取整Randpos=[Randpos,randperm(n)];%生成ceil(m/n)个1*n的矩阵并合并endTabu(:,1)=(Randpos(1,1:m))'; %将Randpos的第一行前m个放到Tabu的第1列
%% 第三步：m 只蚂蚁按概率函数选择下一座城市，完成各自的周游
%将 m 个蚂蚁置于n个城市上，计算待选城市的概率分布，m 只蚂蚁按概率函数选择下一座城市，完成各自的周游。for j=2:n%第j个城市for i=1:m%第i个蚂蚁visited=Tabu(i,1:(j-1));  %己访问的城市J=zeros(1,(n-j+1));       %待访问的城市P=J;                      %待访问城市的选择概率分布Jc=1;for k=1:nif length(find(visited==k))==0%判断第k个城市有没有被访问J(Jc)=k;Jc=Jc+1;endend
%计算待选城市的概率分布for k=1:length(J)P(k)=(Tau(visited(end),J(k))^Alpha)...*(Eta(visited(end),J(k))^Beta);end%% 修改%开始if sum(P)==0to_visit=J(ceil(length(J)*rand)); %% 如果所选择的全部城市信息素为0将随机选择sumnum=sumnum+1;else%结束    P=P/(sum(P));
%按概率原则选取下一个城市Pcum=cumsum(P);         % 如 P=[1 2 3 4]，则cumsum(P)=[1 3 6 10],要累加，轮盘赌法，依次看是否在转得的区域内Select=find(Pcum>=rand);%轮盘赌法随机选择to_visit=J(Select(1));  %待选择的城市endTabu(i,j)=to_visit;     %访问的城市endendif NC>=2Tabu(1,:)=R_best(NC-1,:);   %最佳路线end
%% 第四步：记录本次迭代最佳路线L=zeros(m,1);%m*1的距离矩阵for i=1:mR=Tabu(i,:);  %第i只蚂蚁走过的城市for j=1:(n-1)L(i)=L(i)+D(R(j),R(j+1));  %计算i只蚂蚁走过的路径endL(i)=L(i)+D(R(1),R(n)); %加上初始位置的路径endL_best(NC)=min(L);      %获取路径最少的一只pos=find(L==L_best(NC));%在m只蚂蚁中寻找走最少的一只R_best(NC,:)=Tabu(pos(1),:);%记录最佳路径
%% 改进算法%开始
if NC>1
if L_best(NC)==L_best(NC-1)num_G=num_G+1;if num_G >= num_G_maxif Rho>=Rho_minRho=0.95*Rho;if Rho<=Rho_minRho=Rho_min;endnum_G=0;endend
end
end
Rho_list(NC)=Rho;
%结束%% 第五步：更新信息素%蚁环算法更新信息素增量%离线更新——蚁群Delta_Tau=zeros(n,n);for i=1:mfor 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);%信息素增量的更新endDelta_Tau(Tabu(i,n),Tabu(i,1))=...Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);%加上起使信息素增量的更新endTau=(1-Rho).*Tau+Delta_Tau;  %更新公式
%% 第六步: 禁忌表清零Tabu=zeros(m,n);
%历代最优路线for i=1:n-1plot([ 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;endplot([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);NC=NC+1;   %迭代+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('适应度进化曲线')
figure(3),
plot(Rho_list)
xlabel('迭代次数')
ylabel('挥发系数')
title('挥发系数自适应曲线')