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: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; %计算两两城市之间的距离
else
D(i,j)=eps;
end
D(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)];
end
Tabu(:,1)=(Randpos(1,1:m))';
%% 第三步:m 只蚂蚁按概率函数选择下一座城市,完成各自的周游
%将 m 个蚂蚁置于n个城市上,计算待选城市的概率分布,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); % 如 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; %访问的城市
end
end
if NC>=2
Tabu(1,:)=R_best(NC-1,:); %最佳路线
end
%% 第四步:记录本次迭代最佳路线
L=zeros(m,1);
for i=1:m
R=Tabu(i,:); %第i只蚂蚁走过的城市
for j=1:(n-1)
L(i)=L(i)+D(R(j),R(j+1)); %计算i只蚂蚁走过的路径
end
L(i)=L(i)+D(R(1),R(n)); %加上初始位置的路径
end
L_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: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);%信息素增量的更新
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);
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: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; %计算两两城市之间的距离
else
D(i,j)=eps;%eps表示从 1.0 到下一个最大双精度数的距离
end
D(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的矩阵并合并
end
Tabu(:,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:n
if length(find(visited==k))==0%判断第k个城市有没有被访问
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
%% 修改
%开始
if sum(P)==0
to_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)); %待选择的城市
end
Tabu(i,j)=to_visit; %访问的城市
end
end
if NC>=2
Tabu(1,:)=R_best(NC-1,:); %最佳路线
end
%% 第四步:记录本次迭代最佳路线
L=zeros(m,1);%m*1的距离矩阵
for i=1:m
R=Tabu(i,:); %第i只蚂蚁走过的城市
for j=1:(n-1)
L(i)=L(i)+D(R(j),R(j+1)); %计算i只蚂蚁走过的路径
end
L(i)=L(i)+D(R(1),R(n)); %加上初始位置的路径
end
L_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_max
if Rho>=Rho_min
Rho=0.95*Rho;
if Rho<=Rho_min
Rho=Rho_min;
end
num_G=0;
end
end
end
end
Rho_list(NC)=Rho;
%结束
%% 第五步:更新信息素
%蚁环算法更新信息素增量
%离线更新——蚁群
Delta_Tau=zeros(n,n);
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);%信息素增量的更新
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);
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('挥发系数自适应曲线')