% I.清空环境变量
clear all
clc

% II. 符号说明
%{

iter_max                          -- 最大迭代次数
m                                 -- 蚂蚁数量,一般设置为目标数的1.5倍
D(i,j)                            -- 城市i和j之间的距离
Eta(i, j) = 1./ D(i, j)           -- 启发函数
alpha                             -- 表征信息素重要程度的参数,通常[1, 4]之间
beta                              -- 表征启发函数重要程度的参数,通常[0, 5]之间
rho                               -- 信息素挥发因子,通常[0.1,0.2, 0.5]之间
Q                                 -- 常系数
route_best                        -- 各代最佳的路线
length_best                       -- 各代最佳路线的长度
length_average                    -- 各代的平均长度

%}

% III.导入数据
load cities.mat

%  IV. 计算城市建的距离（欧式距离）
n = size(cities,1); %城市个数
D = zeros(n,n);
for i = 1:n
    for j = i+1:n
        %D(i,j) = sqrt(sum((cities(1,:)-cities(j,:)).^2));  
        D(i, j) = sqrt((cities(i, 1) - cities(j, 1))^2 + (cities(i, 2) - cities(j, 2))^2);
        D(j,i) = D(i,j);
    end
    D(i,i) = 1e-4;
end

% V. 初始化参数
m = 50;             % 蚂蚁数量
alpha = 1;          %信息素重要程度因子
beta = 5;           % 启发函数重要程度因子
rho = 0.1;          % 信息素挥发因子
Q = 20;            %常系数
Eta = 1./D;         % 启发函数
Tau = ones(n,n);    %信息素矩阵,Tau(i, j)表示边(i, j)的信息素量，一开始都为1
Table = zeros(m,n); %路径记录表
iter = 1;           %迭代次数初值
iter_max = 200;     %最大迭代次数
route_best = zeros(iter_max,n); %各代最佳路径
length_best = zeros(iter_max,1); %各代最佳路径长度
length_average = zeros(iter_max,1); %各代路径的平均长度

% VI.迭代寻找最佳路径
while iter <= iter_max
    
    %step 1,随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n);     %1-n的随机数
        start(i) = temp(1);
    end
    Table(:,1) = start;     %Table第一列即是所有蚂蚁的起点城市
    cities_index = 1:n;     %所有城市索引的一个集合
    
    %step 2,逐个蚂蚁路径选择
    for i = 1:m
        %逐个城市路径选择
        for j = 2:n
            visited = Table(i,1:(j-1));     % 蚂蚁i已访问的城市集合(禁忌表)
            allow_index = ~ismember(cities_index, visited);
            allow = cities_index(allow_index);   % 存放待访问的城市集合
            P = allow;
            
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(visited(end), allow(k))^alpha * Eta(visited(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
    
    % step 3,计算各个蚂蚁的路径距离
    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

    % step 4,计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        length_best(iter) = min_Length;
        length_average(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_average(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
    
    % step 5,更新信息素
    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;
    
    %step 6,迭代次数加1, 清空路径记录表
    iter = iter + 1;
    Table = zeros(m, n);
end

% VII. 结果显示
[Shortest_Length, index] = min(length_best);
Shortest_Route = route_best(index, :);
disp([\'最短距离\', num2str(Shortest_Length)]);
disp([\'最短路径\', num2str([Shortest_Route Shortest_Route(1)])]);

% VIII. 绘图
figure(1);
plot([cities(Shortest_Route,1);cities(Shortest_Route(1),1)],...
     [cities(Shortest_Route,2);cities(Shortest_Route(1),2)],\'o-\');
grid on;
for i = 1:size(cities,1)
    text(cities(i,1),cities(i,2),[\'   \' num2str(i)]);
end
text(cities(Shortest_Route(1),1),cities(Shortest_Route(1),2),\'       起点\');
text(cities(Shortest_Route(end),1),cities(Shortest_Route(end),2),\'       终点\');
xlabel(\'城市位置横坐标\')
ylabel(\'城市位置纵坐标\')
title([\'蚁群算法优化路径(最短距离:\' num2str(Shortest_Length) \')\'])
figure(2)
plot(1:iter_max,length_best,\'b\',1:iter_max,length_average,\'r:\')
legend(\'最短距离\',\'平均距离\')
xlabel(\'迭代次数\')
ylabel(\'距离\')
title(\'各代最短距离与平均距离对比\')