%符号说明
%SumOfant表示蚂蚁数量
%SumOfCity表示城市数量
%sight为距离的倒数，表示每条边的能见度
%Q表示信息素增强的系数
%nc_now表示当前的迭代次数
%nc_max表示自主设置的最大迭代次数  表示终止条件
%first_address表示测试数据中的城市坐标
%RouteOfBest表示各代的最佳路线  
%LengthOfBest表示每次迭代的最佳路线的长度
%a表示启发因子，表示信息素的相对重要程度
%b表示期望因子，表示能见度的重要性
%p表示信息素的蒸发系数，（1-p）表示信息素持久性系数
%length_address表示两两城市间的距离


for   n=1:size(first_address)
for m=1:size(first_address)
length_address(n,m)=[(first_address(n,1)-first_address(m,1))^2+(first_address(n,2)-first_address(m,2))^2]^1/2;
end
end

sight=1./length_address;%表示每条边的能见度
SumOfCity=20;
>> SumOfant=40;
>> sight=1./length_address;
info_pre=ones(SumOfCity,SumOfCity);     %info_pre为信息素矩阵,可以理解为在蚂蚁还没有被放入城市前，每条道路上就已经存在了一定含量的信息素（只是为了方便计算）
>> EachOfRoute=zeros(SumOfCity,SumOfCity);  
> %存储并记录每次迭代时每只蚂蚁经历的路径生成；
>> nc_now=1;  nc_max=60;              %迭代计数器，记录迭代次数,此处设置为60次迭代
>> RouteOfBest=zeros(nc_max,SumOfCity);       %各代最佳路线
Length_best=inf.*ones(nc_max,1);   %各代最佳路线的长度
LengthOfAverage=zeros(nc_max,1);        %各代路线的平均长度


while nc_now<=nc_max%表示循环终止条件，迭代终止器
%%第二步：将m只蚂蚁放到n个城市上
Randpos=[];   %随即存取

for i=1:(ceil(SumOfant/SumOfCity))
%ceil为取整函数，表示每个城市放几只蚂蚁
Randpos=[Randpos,randperm(SumOfCity)];%randperm表示1-20随机分配，也就是每只蚂蚁的位置
end
%表示每只蚂蚁一开始所被放置的位置%第三步：所有蚂蚁按概率函数选择下一座城市，完成各自的周游
for j=2:SumOfCity    %每只蚂蚁从第一个城市出发，开始依次访问
for i=1:SumOfant
CityOfVisited=EachOfRoute(i,1:(j-1)); %记录已访问的城市，避免重复访问
J=zeros(1,(SumOfCity-j+1));       %待访问的城市
P=J;                      %待访问城市的选择概率分布,用P表示
Jc=1;%记录已经访问的城市数目
for k=1:SumOfCity
if length(find(CityOfVisited==k))==0   %开始时置0
J(Jc)=k;
Jc=Jc+1;                         %访问的城市个数自加1
end
end


%下面计算待选城市的概率分布
for k=1:length(J)%对每只蚂蚁还没有访问的城市依次计算概率
P(k)=(info_pre(CityOfVisited(end),J(k))^0.7)*(sight(CityOfVisited(end),J(k))^3.8);%在这里设置启发因子=0.7,期望因子=3.8
end
P=P/(sum(P));%按概率原则选取下一个城市
Pcum=cumsum(P);     %cumsum，元素累加即求和
Select=find(Pcum>=rand); %若计算的概率大于原来的就选择这条路线
CityOfNextVisit=J(Select(1));
EachOfRoute(i,j)=CityOfNextVisit;
end
end

if nc_now>=2
EachOfRoute(1,:)=RouteOfBest(nc_now-1,:);%如果迭代次数大于2，则将前依次迭代的访问顺序存入EachOfRoute
end%%第四步：记录本次迭代最佳路线
L=zeros(SumOfant,1);     %开始距离为0，m*1的列向量
for i=1:SumOfant
R=EachOfRoute(i,:);
for j=1:(SumOfCity-1)%R表示访问的城市编号,然后依次在length_address找到两两城市间的距离，计算路程
L(i)=L(i)+length_address(R(j),R(j+1));
end
L(i)=L(i)+length_address(R(1),R(n));    
end
Length_best(nc_now)=min(L);           %最佳距离取最小
pos=find(L==Length_best(nc_now));
RouteOfBest(nc_now,:)=EachOfRoute(pos(1),:); %此轮迭代后的最佳路线
LengthOfAverage(nc_now)=mean(L);           %此轮迭代后的平均距离
nc_now=nc_now+1 ; 


%%第五步：更新信息素
Update_info=zeros(SumOfCity,SumOfCity);        %开始时信息素为n*n的0矩阵
for i=1:SumOfant
for j=1:(SumOfCity-1)
Update_info(EachOfRoute(i,j),EachOfRoute(i,j+1))=Update_info(EachOfRoute(i,j),EachOfRoute(i,j+1))+400/L(i); %400表示表示信息素增强的系数         
%此次循环在路径（i，j）上的信息素增量
end
Update_info(EachOfRoute(i,n),EachOfRoute(i,1))=Update_info(EachOfRoute(i,n),EachOfRoute(i,1))+400/L(i);
%此次循环在整个路径上的信息素增量
end
info_pre=(1-0.75).*info_pre+Update_info; %考虑信息素挥发，更新后的信息素
%%第六步：禁忌表清零
EachOfRoute=zeros(SumOfant,SumOfCity);
end

Pos=find(Length_best==min(Length_best)); %找到最佳路径（非0为真）
Shortest_Route=RouteOfBest(Pos(1),:) %最大迭代次数后最佳路径
Shortest_Length=Length_best(Pos(1)) %最大迭代次数后最短距离
subplot(1,2,1)                  %绘制第一个子图形
DrawRoute(C,Shortest_Route)     %画路线图的子函数
subplot(1,2,2)                  %绘制第二个子图形
plot(Length_best)
hold on                         %保持图形
plot(LengthOfAverage,'r')
title('平均距离和最短距离')     %标题