#============導(dǎo)入相關(guān)庫=================
import numpy as np
from scipy import spatial
import pandas as pd
import matplotlib.pyplot as plt
from sko.ACA import ACA_TSP
 
num_points = 25
 
points_coordinate = np.random.rand(num_points, 2)  # 生成點(diǎn)的坐標(biāo)
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')#函數(shù)用于計(jì)算兩個(gè)輸入集合的距離
 
 
def cal_total_distance(routine):
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
 
 
#=============ACA_TSP解決==================================
 
aca = ACA_TSP(func=cal_total_distance, n_dim=num_points,
              size_pop=50, max_iter=200,
              distance_matrix=distance_matrix)
 
best_x, best_y = aca.run()
 
#=============可視化=======================
 
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_x, [best_x[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
pd.DataFrame(aca.y_best_history).cummin().plot(ax=ax[1])
plt.show()

%蟻群算法尋找最短路徑程序
%% 清空環(huán)境變量
clear
clc
%% 導(dǎo)入數(shù)據(jù)，導(dǎo)入方式，看個(gè)人習(xí)慣
zuobiao=[1304  2312
3639  1315
4177  2244
3712  1399
3488  1535
3326  1556
3238  1229
4196  1004
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  2367
3394  2643
3439  3201
2935  3240
3140  3550
2545  2357
2778  2826
2370  2975];
%% 計(jì)算城市間相互距離
n = size(zuobiao,1);%城市個(gè)數(shù)
jl = zeros(n,n);%首先求得各個(gè)坐標(biāo)點(diǎn)的距離，這里是矩陣初始化
for i = 1:n
    for j = 1:n
        if i ~= j  %~=是不等于的意思，zuobiao矩陣中每行都有個(gè)坐標(biāo)，坐標(biāo)相減用i和j區(qū)分不同的坐標(biāo)點(diǎn)，然后求兩點(diǎn)距離
            jl(i,j) = sqrt(sum((zuobiao(i,:) - zuobiao(j,:)).^2));
%上式運(yùn)算如a=[2,2;1,1]=>a(1,:)-a(2,:)=>ans =1 1，然后1?+1?=2，最后開根號(hào)
        else
            jl(i,j) = 1e-4;%相等的點(diǎn)相減準(zhǔn)確說是等于0的，這里設(shè)置成了一個(gè)很小的數(shù)，是為了避免后面程序運(yùn)算出錯(cuò)
        end
    end    
end
%% 初始化參數(shù)
m = 50;         % 螞蟻數(shù)量，視情況而定，坐標(biāo)點(diǎn)多的話可以適當(dāng)增加螞蟻數(shù)量
a= 1;           % 信息素重要程度因子
b= 5;           % 啟發(fā)函數(shù)重要程度因子
r = 0.1;        % 信息素?fù)]發(fā)因子
Q = 1;          % 常數(shù)
qfhs = 1./jl;    % 啟發(fā)函數(shù),將jl矩陣中每個(gè)元素轉(zhuǎn)化為倒數(shù)
xxsjz = ones(n,n);       % 信息素矩陣初始化
ljjl = zeros(m,n);       % 路徑記錄表矩陣初始化
ddcs = 1;                % 迭代次數(shù)初值
ddcs_max = 200;          % 最大迭代次數(shù) 
Lujin_best = zeros(ddcs_max,n);      % 各代最佳路徑       
L_best = zeros(ddcs_max,1);     % 各代最佳路徑的長度  
L_ave = zeros(ddcs_max,1);      % 各代路徑的平均長度  
%% 迭代尋找最佳路徑
while ddcs <= ddcs_max%在ddcs小于ddcs_max前，一直循環(huán)
%% 隨機(jī)產(chǎn)生各個(gè)螞蟻的起點(diǎn)
      start = zeros(m,1);
      for i = 1:m
          temp = randperm(n);%功能是隨機(jī)打亂一個(gè)數(shù)字序列,也就是現(xiàn)將坐標(biāo)點(diǎn)排號(hào)再打亂，相當(dāng)于將螞蟻隨機(jī)分布在各個(gè)地點(diǎn)
          start(i) = temp(1);
      end
      ljjl(:,1) = start; 
%% 構(gòu)建解空間
      zuobiao_index = 1:n;
      % 逐個(gè)螞蟻路徑選擇
      for i = 1:m
          % 逐個(gè)地點(diǎn)路徑選擇
         for j = 2:n
             yfw = ljjl(i,1:(j - 1));           % 已訪問的地點(diǎn)集合(禁忌表)
             allow_index = ~ismember(zuobiao_index,yfw);%ismember用于判斷矩陣某個(gè)元素是否存在，用法詳見后文函數(shù)講解
             allow = zuobiao_index(allow_index);  % 待訪問的城市集合
             P = allow;
             % 計(jì)算城市間轉(zhuǎn)移概率
             for k = 1:length(allow)
                 P(k) = xxsjz(yfw(end),allow(k))^a * qfhs(yfw(end),allow(k))^b;%見文中公式
             end
             P = P/sum(P);
             % 選擇下一個(gè)訪問城市
             Plj = cumsum(P);     %cumsum函數(shù)用于累加，具體用法詳見后文函數(shù)講解
             yidong_index = find(Plj >= rand); 
             yidong = allow(yidong_index(1));
             ljjl(i,j) = yidong;
         end
      end
      % 計(jì)算各個(gè)螞蟻的路徑距離
      L = zeros(m,1);
      for i = 1:m
          Lujin = ljjl(i,:);
          for j = 1:(n - 1)
              L(i) = L(i) + jl(Lujin(j),Lujin(j + 1));
          end
          L(i) = L(i) + jl(Lujin(n),Lujin(1));
      end
      % 計(jì)算最短路徑距離及平均距離
      if ddcs == 1
          [min_L,min_index] = min(L);
          L_best(ddcs) = min_L;  
          L_ave(ddcs) = mean(L);
          Lujin_best(ddcs,:) = ljjl(min_index,:);
      else
          [min_L,min_index] = min(L);
          L_best(ddcs) = min(L_best(ddcs - 1),min_L);
          L_ave(ddcs) = mean(L);
          if L_best(ddcs) == min_L
              Lujin_best(ddcs,:) = ljjl(min_index,:);
          else
              Lujin_best(ddcs,:) = Lujin_best((ddcs-1),:);
          end
      end
%% 更新信息素
      S = zeros(n,n);
      % 逐個(gè)螞蟻計(jì)算
      for i = 1:m
          % 逐個(gè)城市計(jì)算
          for j = 1:(n - 1)
              S(ljjl(i,j),ljjl(i,j+1)) = S(ljjl(i,j),ljjl(i,j+1)) + Q/L(i);
          end
          S(ljjl(i,n),ljjl(i,1)) = S(ljjl(i,n),ljjl(i,1)) + Q/L(i);
      end
      xxsjz = (1-r) * xxsjz + S;
    % 迭代次數(shù)加1，清空路徑記錄表
    ddcs = ddcs + 1;
    ljjl = zeros(m,n);
end
%% 結(jié)果顯示
[Shortest_L,index] = min(L_best);
Shortest_Lujin = Lujin_best(index,:);
disp(['最短距離:' num2str(Shortest_L)]);
disp(['最短路徑:' num2str([Shortest_Lujin Shortest_Lujin(1)])]);
%% 繪圖
figure(1)
plot([zuobiao(Shortest_Lujin,1);zuobiao(Shortest_Lujin(1),1)],...
     [zuobiao(Shortest_Lujin,2);zuobiao(Shortest_Lujin(1),2)],'o-');
grid on
for i = 1:size(zuobiao,1)
    text(zuobiao(i,1),zuobiao(i,2),['   ' num2str(i)]);
end
text(zuobiao(Shortest_Lujin(1),1),zuobiao(Shortest_Lujin(1),2),'       起點(diǎn)');
text(zuobiao(Shortest_Lujin(end),1),zuobiao(Shortest_Lujin(end),2),'       終點(diǎn)');
xlabel('城市位置橫坐標(biāo)')
ylabel('城市位置縱坐標(biāo)')
title(['蟻群算法優(yōu)化路徑(最短距離:' num2str(Shortest_L) ')'])
figure(2)
plot(1:ddcs_max,L_best,'b',1:ddcs_max,L_ave,'r')
legend('最短距離','平均距離')
xlabel('迭代次數(shù)')
ylabel('距離')
title('各代最短距離與平均距離對(duì)比')