% 灰狼优化算法GWO （求TSP问题）
clc;
clear;
%% 初始化参数
N=100;            % 灰狼个数
Max_iter=1000;    % 最大迭代次数
load('city.mat');       % 导入城市数据
City=city;              % 保存城市位置用于计算距离
M=size(city,1);         % 得到TSP问题的规模,即城市数目
dim=M;      % 维度
 
figure('name','灰狼优化算法');
plot(city(:,1),city(:,2),'o');  %描点   city--（30，2）
for i=1:M
    text(city(i,1)+0.5,city(i,2),num2str(i));       %标号  
end
title('城市初始位置');
 
% 初始化alpha,beta,delta
Alpha_pos=zeros(1,dim);    
Alpha_score=inf; 
Beta_pos=zeros(1,dim);
Beta_score=inf; 
Delta_pos=zeros(1,dim);
Delta_score=inf;
 
% 初始化种群位置
Positions=zeros(N,dim);
for i=1:N
    Positions(i,:)=randperm(dim);          %随机排列，比如[2 4 5 6 1 3]
end
 
Length_best = zeros(1, Max_iter);       % 定义每次迭代的最短距离
Length_ave = zeros(1, Max_iter);        % 定义每次迭代的平均距离
l = 1;   % 迭代计数器
%% 迭代寻优
while l < Max_iter+1
    for i = 1:N
        % 按行升序排列产生城市序列
        [~, sol] = sort(Positions, 2);             
        % 计算目标函数值(即路径距离)
        fitness = Fun(sol(i, :),City,M);
        Length_ave(l) = Length_ave(l)+fitness;
        % 更新Alpha, Beta, and Delta
        if fitness < Alpha_score 
            Alpha_score = fitness; 
            Alpha_pos = Positions(i, :);
        end     
        if fitness > Alpha_score && fitness < Beta_score 
            Beta_score = fitness; 
            Beta_pos = Positions(i, :);
        end  
        if fitness > Alpha_score && fitness > Beta_score && fitness < Delta_score 
            Delta_score = fitness; 
            Delta_pos = Positions(i, :);
        end
    end      
    a = 2-l*((2)/Max_iter); % a从2线性减到0    
    
    % 更新所有个体位置
    for i = 1:N
        for j = 1:dim                  
            r1 = rand(); % r1 is a random number in [0,1]
            r2 = rand(); % r2 is a random number in [0,1]   
            A1 = 2*a*r1-a;       
            C1 = 2*r2;              
            D_alpha = abs(C1*Alpha_pos(j)-Positions(i, j));
            X1 = Alpha_pos(j)-A1*D_alpha; 
                       
            r1 = rand();
            r2 = rand();     
            A2 = 2*a*r1-a;   
            C2 = 2*r2;         
            D_beta = abs(C2*Beta_pos(j)-Positions(i, j));  
            X2 = Beta_pos(j)-A2*D_beta;                              
            
            r1 = rand();
            r2 = rand();      
            A3 = 2*a*r1-a;     
            C3 = 2*r2;            
            D_delta = abs(C3*Delta_pos(j)-Positions(i, j));     
            X3 = Delta_pos(j)-A3*D_delta;                                       
            
            Positions(i, j) = (X1+X2+X3)/3;       
        end
    end  
    Length_best(l) = Alpha_score;               % 最短距离
    Length_ave(l) = Length_ave(l)/dim;          % 平均距离
    disp(['Iteration ' num2str(l) ': Best Fitness = ' num2str(Alpha_score)]);  % 命令行窗口展示
    l = l + 1;      % 更新迭代次数
    [~, BestSol] = sort(Alpha_pos);         % 得到最优解
end
 
figure(2);
for i=1:M-1
    plot([City(BestSol(i),1),City(BestSol(i+1),1)],[City(BestSol(i),2),City(BestSol(i+1),2)],'bo-');
    hold on;
end
plot([City(BestSol(M),1),City(BestSol(1),1)],[City(BestSol(M),2),City(BestSol(1),2)],'ro-');
for i=1:M
    text(City(i,1)+0.5,City(i,2),num2str(i));       %标号  
end
text(City(BestSol(1),1),City(BestSol(1),2),'    起点');
text(City(BestSol(M),1),City(BestSol(M),2),'    终点');
title('最终路线图');
 
figure(3);
t = 1:Max_iter;
plot(t, Length_best, 'r', t, Length_ave, 'b--', 'LineWidth', 2);
xlabel('Iteration');
ylabel('Best Cost');
legend('最短距离','平均距离')
xlabel('迭代次数')
ylabel('距离')
title('各代最短距离与平均距离对比')
P=zeros(1, 30);
 
function len=Fun(sol,City,M)  % 距离计算函数
len=0;
for k=1:M-1
    len=len+sqrt(sum((City(sol(1,k),:)-City(sol(1,k+1),:)).^2));
end
len=len+sqrt(sum((City(sol(1,M),:)-City(sol(1,1),:)).^2));
end