% function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACO(C,D,s,e,NC_max,m,Alpha,Beta,Rho,Q)
% function [Shortest_Route,Shortest_Length]=ACOR(C,D,s,e,NC_max,m,Alpha,Beta,Rho,Q)
%%=========================================================================
%% ACO.m
%% Ant Colony Optimization Algorithm for Road Select Problem
%% LiLixin,ShenYang Insitute of Aeronautical engineering ,ShenYang,China
%% Email:myassist@http://www.51wendang.com
%% All rights reserved
%%-------------------------------------------------------------------------
%% 主要符号说明
%% C n个城市的坐标，n×2的矩阵
%% D 道路连通加权矩阵，n×n的矩阵
%% s 起点
%% e 终点
%% NC_max 最大迭代次数
%% m 蚂蚁个数
%% Alpha 表征信息素重要程度的参数
%% Beta 表征启发式因子重要程度的参数
%% Rho 信息素蒸发系数
%% Q 信息素增加强度系数
%% R_best 各代最佳路线
%% L_best 各代最佳路线的长度
%%=========================================================================
%
clc
clear
% 设置初始参数如下：
m=10;Alpha=1;Beta=5;Rho=0.1;NC_max=20;Q=100;
%设定起始点
s=1;e=50;
% 城市坐标为：
load C.txt;
%%第一步：变量初始化
n=size(C,1);%n表示问题的规模（城市个数）
hold on
XX=C';
plot(XX(1 , :) , XX(2 , :) , 'k+', 'markersize' , 5) %画十字架
for i=1:n
text(C(i,1)+0.5,C(i,2)-0.5,int2str(i)); %加标号
end
load D.txt;
Eta=1./D;%Eta为启发因子，这里设为距离的倒数， ./表示元素对元素除法
Tau=ones(n,n);%Tau为信息素矩阵
Tabu=zeros(m,n);%存储并记录路径的生成
NC=1;%迭代计数器
R_best=zeros(NC_max,n);%各代最佳路线
L_best=inf.*ones(NC_max,1);%各代最佳路线的长度 .*元素对元素乘法
lastmin=inf; %上代最小路径
thesameNum=0; %终止算法条件之一
while NC<=NC_max%停止条件之一：达到最大迭代次数
%%第二步：将m只蚂蚁放到n个城市上
Randpos=ones(1,m);
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); %cumsum，元素累加即求和
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，