%%清空环境变量
clear all
clc
%%导入数据
load citys_data.mat
%%计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
%计算城市两两之间的距离
for i =1:n
    for j =1:n
        if i~=j
            D(i,j)=sqrt(sum( ( citys(i,:) - citys(j,:) ).^2 ));
        else
            D(i,j)= 1e-4;
        end
    end
end
%%初始化参数
m = 31;%蚂蚁
alpha =1;%信息素重要程度因子
beta = 5;%启发函数重要程度因子
rho = 0.1;%信息素挥发因子
Q=1;%常系数
Eta = 1./D;%启发函数
Tau = ones(n,n);%信息素矩阵
Table = zeros(m,n); %路径记录表
iter =1;
iter_max = 200; %最大迭代次数
route_best = zeros(iter_max,n);%各代最佳路径
length_best = zeros(iter_max,1);%各代最佳路径的长度
length_ave = zeros(iter_max,1);%各代最佳路径的平均长度
%%迭代寻找最优路径
while iter<=iter_max
    %随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i =1:m
        temp = randperm(n);%随机生成了一个路径
        start(i) = temp(1);%生成路径的起点
    end
    Table(:,1)=start;
    %构建解空间
    city_index = 1:n;
    %逐个蚂蚁路径选择
    for i =1:m
        %逐个城市选择
        for j= 2:n
            tabu = Table(i,1:(j-1));%已访问的城市集合（禁忌表）
            allow_index = ~ismember(city_index,tabu);
            allow = city_index(allow_index);%待访问的城市集合
            P = allow;
            %计算城市转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha*Eta(tabu(end),allow(k))^beta;
            end
            P = P/sum(P);
            %使用轮盘赌方法选择下一个访问城市
            Pc =cumsum(P);
            target_index = find(Pc>=rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    %计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j= 1:(n-1)
            Length(i)=Length(i)+D(Route(j),Route(j+1));
        end
        Length(i) = Length(i)+D(Route(n),Route(1));
    end
    %计算最短路径与平均距离
    if iter ==1
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter-1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
        else
            Route_best(iter,:)=Route_best((iter-1),:);
        end
    end
    %更新信息素
    Delta_Tau = zeros(n,n);
    for i=1:m
        %逐个城市计算
        for j=1:(n-1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1))+Q/Length(i); 
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1))+Q/Length(i);
    end
    Tau = (1-rho)*Tau+Delta_Tau;
    %迭代次数加1，清空历史记录
    iter = iter+1;
    Table = zeros(m,n);
end
%% 结果显示
[shortest_length,index] = min(Length_best);
shortest_Route = Route_best(index,:);
disp(['最短距离：' num2str(shortest_length)]);
disp(['最短路径：' num2str([shortest_Route shortest_Route(1)])]);

%% 绘图
figure(1)
plot([citys(shortest_Route,1);citys(shortest_Route(1),1)],...
    [ citys(shortest_Route,2);citys(shortest_Route(1),2) ],'o-')
grid on
for i = 1:size(citys,1)
    text(citys(i,1),citys(i,2),[' ' num2str(i)]);
end
text(citys(shortest_Route(1),1),citys(shortest_Route(1),2),'起点')
text(citys(shortest_Route(end),1),citys(shortest_Route(end),2),'起点')
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['蚁群算法优化路径（最短距离：' num2str(shortest_length) ')'])

figure(2)
plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r')
legend('最短距离','平均距离')

　　内容参考  matlab 智能计算30个案例分析