深圳幻海软件技术有限公司 欢迎您！
更多

灰狼优化算法（GWO）（解决TSP问题，代码完整免费）

2023-04-02

end city rand
算法背景灰狼优化算法(GWO)，由澳大利亚格里菲斯大学学者Mirjalili等人于2014年提出来的一种群智能优化算法。灵感来自于灰狼群体捕食行为。优点：较强的收敛性能，结构简单、需要调节的参数少，容易实现，存在能够自适应调整的收敛因子以及信息反馈机制，能够在局部寻优与全局搜索之间实现平衡，因此在对

算法背景

灰狼优化算法(GWO)，由澳大利亚格里菲斯大学学者 Mirjalili 等人于2014年提出来的一种群智能优化算法。灵感来自于灰狼群体捕食行为。

优点：较强的收敛性能，结构简单、需要调节的参数少，容易实现，存在能够自适应调整的收敛因子以及信息反馈机制，能够在局部寻优与全局搜索之间实现平衡，因此在对问题的求解精度和收敛速度方面都有良好的性能。

缺点：存在着易早熟收敛，面对复杂问题时收敛精度不高，收敛速度不够快

应用：车间调度、参数优化、图像分类、路径规划。

算法思想

灰狼群体中有严格的等级制度，一小部分拥有绝对话语权的灰狼带领一群灰狼向猎物前进。灰狼群一般分为4个等级：  （权利从大到小）模拟领导阶层。

集体狩猎是灰狼的一种社会行为，社会等级在集体狩猎过程中发挥着重要的作用，捕食的过程在α的带领下完成。主要包括三个步骤：

  • 跟踪和接近猎物
  • 骚扰、追捕和包围猎物，直到它停止移动
  • 攻击猎物

社会等级分层

构建灰狼社会等级层次模型，对灰狼的社会等级进行数学建模。

考虑一个搜索空间，现实中狼群由眼睛嗅觉等等知道猎物的位置，然而计算机中模拟无法做到这点，但是我们可以假设，猎物位置是由我们在搜索空间中发现的最佳解决方案提供，我们可以使用该解决方案去找到更好的解决方案，不断迭代优化，此最优解决方案会非常接近最好的解决方案。

将  作为最优解（个体的适应度最优），次优解   ，最佳解决方案  ，剩下的候选解命名为 ， 狩猎过程由  引导， 跟随这三只狼。即我们总是去找到三个最佳解决方案，然后围绕该区域进行搜索，目的是找到更好的解决方案然后更新 。

包围猎物

在狩猎过程中，将灰狼围捕猎物的行为定义如下：

个体与猎间的距离公式：

        （1）

灰狼位置更新公式：

        （2）

系数向量：

        （3）

        （4）

其中：t 为迭代次数，D 是个体与猎间的距离向量， ‘’ 不是点积，是乘法， 是猎物位置向量，X 为灰狼位置向量，  是收敛因子（随迭代次数从 2 线性减小到 0 ），r1 和 r2 是随机向量，模取【0-1】之间的随机数。

从公式可以看出，经过移动后灰狼群向  移动，移动方向由自身位置和随机向量 C 决定，移动步长由隔离距离灰狼距离和系数向量 A 决定，即  线性减小意味着随机性和运动步长幅度，步长随迭代次数而减小越来越接近最优解。

狩猎

灰狼能够识别猎物的位置并包围它们。当灰狼识别出猎物的位置后,  在  的带领下指导狼群包围猎物。灰狼个体跟踪猎物位置的数学模型描述如下：

                        

                （5）

                        

其中,  ， 和  分别表示   和  , 与其他个体间的距离； 分别代表   和  ,当前的位置； 是随机向量,  X 是当前灰狼的位置。

                        

                （6）

                        

                 （7）

式（6）分别定义了狼群中  个体朝向    和  ,  前进的步长和方向 ，式 (7) 定义了ω的最终位置。

攻击猎物

当猎物停止移动时,灰狼通过攻击来完成狩猎过程.为了模拟逼近猎物,  的值被逐渐减小，因此  A的波动范围也随之减小。换句话说，在迭代过程中,当   的值从２线性下降到０时，其对应的 A 的值也在区间 [ -,  ] 内变化。

如图所示，当  A 的值位于区间内时,灰狼的下一位置可以位于其当前位置和猎物位置之间的任意位置（由于随机），当 -1<A<1 时，狼群向猎物发起攻击(陷入局部最优)。当 A>1 或 A<-1 时，灰狼远离猎物，探索其他区域（为了找到全局最优解）。

GWO 另一个有利于全局探索的组件 C（随机权重） ，C 是[0,2]之间的随机向量。表示狼所在的位置对猎物影响的随机权重，C>1表示影响权重大，反之，表示影响权重小。在算法陷入了局部最优并且不易跳出时， C 的随机性在避免局部最优方面发挥了非常重要的作用。

算法流程

伪代码

  1. Initialize the grey wolf population Xi (i=1,2,,,,n)
  2. Initialize a,A,C,t=0
  3.  
  4. Calculate the fitness of each search agent
  5. X(alpha) = the best search agent
  6. X(bete) = the second best search agent
  7. X(delta) = the third best search agent
  8.  
  9. while (t< Max number of iterations)
  10.     for each search agent
  11.         Update the position of the current search agent
  12.     end for
  13.  
  14.     Update a,A,and C
  15.     Calculate the fitness of all search agents
  16.     Update X(alpha),X(bete),and X(delta)
  17.     t=t+1
  18. end while
  19.  
  20. return X(alpha)

算法测试

求函数最大值（一元）

  1. % 灰狼优化算法（求函数极值）
  2. clc;
  3. clear;
  4. close all;
  5. %% 目标函数
  6. f= @(x) - (x - 10) .^ 2 + x .* sin(x) .* cos(2 * x) - 5 * x .* sin(3 * x) ; % 适应度函数表达式(求这个函数的最大值)  
  7. figure(1);
  8. fplot(f, [0 20], 'b-');                 % 画出初始图像 
  9. title('初始图像');
  10. hold on;
  11. %% 初始化参数
  12. N=30;       % 灰狼个数
  13. dim=1;      % 维度
  14. Iter=50;   % 最大迭代次数
  15. a=2;        % 收敛因子
  16. ub=20;      % 最大值限制
  17. lb=0;       % 最小值限制
  18.  
  19. % 初始化alpha,beta,delta
  20. Alpha_pos=zeros(1,dim);    
  21. Alpha_score=-inf; %求最大值改为inf
  22. Beta_pos=zeros(1,dim);
  23. Beta_score=-inf; 
  24. Delta_pos=zeros(1,dim);
  25. Delta_score=-inf;
  26.  
  27. Positions=rand(N,dim).*(ub-lb)+lb;      % 初始化个体位置
  28. Convergence_curve=zeros(1,Iter);        % 收敛曲线
  29. l=0;        %循环次数记录
  30.  
  31. %% 迭代求解
  32. while l<Iter
  33.     for i=1:size(Positions,1)  
  34.         
  35.        % 超出边界处理
  36.         Flag4ub=Positions(i,:)>ub;
  37.         Flag4lb=Positions(i,:)<lb;
  38.         Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
  39.         
  40.         % 计算个体适应度函数
  41.         fitness=f(Positions(i,:));
  42.         
  43.         % 更新 Alpha, Beta, and Delta
  44.         if fitness>Alpha_score 
  45.             Alpha_score=fitness; 
  46.             Alpha_pos=Positions(i,:);
  47.         end    
  48.         if fitness<Alpha_score && fitness>Beta_score 
  49.             Beta_score=fitness; 
  50.             Beta_pos=Positions(i,:);
  51.         end     
  52.         if fitness<Alpha_score && fitness<Beta_score && fitness>Delta_score 
  53.             Delta_score=fitness; 
  54.             Delta_pos=Positions(i,:);
  55.         end
  56.     end
  57.        
  58.     a=2-l*((2)/Iter); % 收敛因子从2线性递减到0
  59.     
  60.     % 更新灰狼个体的位置
  61.     for i=1:size(Positions,1)
  62.         for j=1:size(Positions,2)     
  63.             
  64.             r1=rand(); % r1 is a random number in [0,1]
  65.             r2=rand(); % r2 is a random number in [0,1]   
  66.             A1=2*a*r1-a;
  67.             C1=2*r2;
  68.             D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); 
  69.             X1=Alpha_pos(j)-A1*D_alpha;
  70.                        
  71.             r1=rand();
  72.             r2=rand();          
  73.             A2=2*a*r1-a;
  74.             C2=2*r2;       
  75.             D_beta=abs(C2*Beta_pos(j)-Positions(i,j));
  76.             X2=Beta_pos(j)-A2*D_beta;
  77.             
  78.             r1=rand();
  79.             r2=rand();          
  80.             A3=2*a*r1-a; 
  81.             C3=2*r2;          
  82.             D_delta=abs(C3*Delta_pos(j)-Positions(i,j));
  83.             X3=Delta_pos(j)-A3*D_delta;
  84.             
  85.             Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)         
  86.         end
  87.     end
  88.     l=l+1;    
  89.     Convergence_curve(l)=Alpha_score;
  90. end
  91. plot(Alpha_pos, f(Alpha_pos), '*r');
  92.  
  93. figure(2);
  94. plot(Convergence_curve);
  95. title('收敛过程'); 
  96.  
  97. display(['The best solution obtained by GWO is : ', num2str(Alpha_pos)]);
  98. display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Alpha_score)]);

二元函数求最大值

  1. % 灰狼优化算法（求函数极值）
  2. clc;
  3. clear;
  4. close all;
  5. %% 目标函数
  6. [ZuoBiao_x, ZuoBiao_y] = meshgrid(-10:0.1:10,-5:0.1:5);   % 产生二维坐标
  7. ZuoBiao_z = f_xy(ZuoBiao_x, ZuoBiao_y);
  8. figure(1);
  9. s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z);                % 画出初始图像 
  10. title('初始图像');
  11. %% 初始化参数
  12. N=30;       % 灰狼个数
  13. dim=2;      % 维度
  14. Iter=200;   % 最大迭代次数
  15. a=2;        % 收敛因子
  16. ub=[10 5];      % 最大值限制
  17. lb=[-10 -5];       % 最小值限制
  18.  
  19. % 初始化alpha,beta,delta
  20. Alpha_pos=zeros(1,dim);    
  21. Alpha_score=-inf; %求最大值改为inf
  22. Beta_pos=zeros(1,dim);
  23. Beta_score=-inf; 
  24. Delta_pos=zeros(1,dim);
  25. Delta_score=-inf;
  26.  
  27. Positions=rand(N,dim).*(ub-lb)+lb;      % 初始化个体位置
  28. Convergence_curve=zeros(1,Iter);        % 收敛曲线
  29. l=0;        %循环次数记录
  30.  
  31. %% 迭代求解
  32. while l<Iter
  33.     for i=1:size(Positions,1)  
  34.         
  35.        % 超出边界处理
  36.         Flag4ub=Positions(i,:)>ub;
  37.         Flag4lb=Positions(i,:)<lb;
  38.         Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
  39.         
  40.         % 计算个体适应度函数
  41.         fitness=f_xy(Positions(i,1), Positions(i,2));
  42.         
  43.         % 更新 Alpha, Beta, and Delta
  44.         if fitness>Alpha_score 
  45.             Alpha_score=fitness; 
  46.             Alpha_pos=Positions(i,:);
  47.         end    
  48.         if fitness<Alpha_score && fitness>Beta_score 
  49.             Beta_score=fitness; 
  50.             Beta_pos=Positions(i,:);
  51.         end     
  52.         if fitness<Alpha_score && fitness<Beta_score && fitness>Delta_score 
  53.             Delta_score=fitness; 
  54.             Delta_pos=Positions(i,:);
  55.         end
  56.     end
  57.        
  58.     a=2-l*((2)/Iter); % 收敛因子从2线性递减到0
  59.     
  60.     % 更新灰狼个体的位置
  61.     for i=1:size(Positions,1)
  62.         for j=1:size(Positions,2)     
  63.             
  64.             r1=rand(); % r1 is a random number in [0,1]
  65.             r2=rand(); % r2 is a random number in [0,1]   
  66.             A1=2*a*r1-a;
  67.             C1=2*r2;
  68.             D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); 
  69.             X1=Alpha_pos(j)-A1*D_alpha;
  70.                        
  71.             r1=rand();
  72.             r2=rand();          
  73.             A2=2*a*r1-a;
  74.             C2=2*r2;       
  75.             D_beta=abs(C2*Beta_pos(j)-Positions(i,j));
  76.             X2=Beta_pos(j)-A2*D_beta;
  77.             
  78.             r1=rand();
  79.             r2=rand();          
  80.             A3=2*a*r1-a; 
  81.             C3=2*r2;          
  82.             D_delta=abs(C3*Delta_pos(j)-Positions(i,j));
  83.             X3=Delta_pos(j)-A3*D_delta;
  84.             
  85.             Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)         
  86.         end
  87.     end
  88.     l=l+1;    
  89.     Convergence_curve(l)=Alpha_score;
  90. end
  91.  
  92. figure(2);
  93. plot(Convergence_curve);
  94. title('收敛过程'); 
  95.  
  96. display(['The best solution obtained by GWO is : ', num2str(Alpha_pos)]);
  97. display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Alpha_score)]);
  98.  
  99.

多元函书求最小值（代码，24个测试函数，感兴趣下载）

链接：https://pan.baidu.com/s/1sQxYMXhwF1xhRqhXB1kJpQ 
提取码：1234 

GWO解TSP问题（网上找了好久，都不全，自己补充一下，简单注释一下，有问题私聊哦）

  1. % 灰狼优化算法GWO （求TSP问题）
  2. clc;
  3. clear;
  4. %% 初始化参数
  5. N=100;            % 灰狼个数
  6. Max_iter=1000;    % 最大迭代次数
  7. load('city.mat');       % 导入城市数据
  8. City=city;              % 保存城市位置用于计算距离
  9. M=size(city,1);         % 得到TSP问题的规模,即城市数目
  10. dim=M;      % 维度
  11.  
  12. figure('name','灰狼优化算法');
  13. plot(city(:,1),city(:,2),'o');  %描点   city--（30，2）
  14. for i=1:M
  15.     text(city(i,1)+0.5,city(i,2),num2str(i));       %标号  
  16. end
  17. title('城市初始位置');
  18.  
  19. % 初始化alpha,beta,delta
  20. Alpha_pos=zeros(1,dim);    
  21. Alpha_score=inf; 
  22. Beta_pos=zeros(1,dim);
  23. Beta_score=inf; 
  24. Delta_pos=zeros(1,dim);
  25. Delta_score=inf;
  26.  
  27. % 初始化种群位置
  28. Positions=zeros(N,dim);
  29. for i=1:N
  30.     Positions(i,:)=randperm(dim);          %随机排列，比如[2 4 5 6 1 3]
  31. end
  32.  
  33. Length_best = zeros(1, Max_iter);       % 定义每次迭代的最短距离
  34. Length_ave = zeros(1, Max_iter);        % 定义每次迭代的平均距离
  35. l = 1;   % 迭代计数器
  36. %% 迭代寻优
  37. while l < Max_iter+1
  38.     for i = 1:N
  39.         % 按行升序排列产生城市序列
  40.         [~, sol] = sort(Positions, 2);             
  41.         % 计算目标函数值(即路径距离)
  42.         fitness = Fun(sol(i, :),City,M);
  43.         Length_ave(l) = Length_ave(l)+fitness;
  44.         % 更新Alpha, Beta, and Delta
  45.         if fitness < Alpha_score 
  46.             Alpha_score = fitness; 
  47.             Alpha_pos = Positions(i, :);
  48.         end     
  49.         if fitness > Alpha_score && fitness < Beta_score 
  50.             Beta_score = fitness; 
  51.             Beta_pos = Positions(i, :);
  52.         end  
  53.         if fitness > Alpha_score && fitness > Beta_score && fitness < Delta_score 
  54.             Delta_score = fitness; 
  55.             Delta_pos = Positions(i, :);
  56.         end
  57.     end      
  58.     a = 2-l*((2)/Max_iter); % a从2线性减到0    
  59.     
  60.     % 更新所有个体位置
  61.     for i = 1:N
  62.         for j = 1:dim                  
  63.             r1 = rand(); % r1 is a random number in [0,1]
  64.             r2 = rand(); % r2 is a random number in [0,1]   
  65.             A1 = 2*a*r1-a;       
  66.             C1 = 2*r2;              
  67.             D_alpha = abs(C1*Alpha_pos(j)-Positions(i, j));
  68.             X1 = Alpha_pos(j)-A1*D_alpha; 
  69.                        
  70.             r1 = rand();
  71.             r2 = rand();     
  72.             A2 = 2*a*r1-a;   
  73.             C2 = 2*r2;         
  74.             D_beta = abs(C2*Beta_pos(j)-Positions(i, j));  
  75.             X2 = Beta_pos(j)-A2*D_beta;                              
  76.             
  77.             r1 = rand();
  78.             r2 = rand();      
  79.             A3 = 2*a*r1-a;     
  80.             C3 = 2*r2;            
  81.             D_delta = abs(C3*Delta_pos(j)-Positions(i, j));     
  82.             X3 = Delta_pos(j)-A3*D_delta;                                       
  83.             
  84.             Positions(i, j) = (X1+X2+X3)/3;       
  85.         end
  86.     end  
  87.     Length_best(l) = Alpha_score;               % 最短距离
  88.     Length_ave(l) = Length_ave(l)/dim;          % 平均距离
  89.     disp(['Iteration ' num2str(l) ': Best Fitness = ' num2str(Alpha_score)]);  % 命令行窗口展示
  90.     l = l + 1;      % 更新迭代次数
  91.     [~, BestSol] = sort(Alpha_pos);         % 得到最优解
  92. end
  93.  
  94. figure(2);
  95. for i=1:M-1
  96.     plot([City(BestSol(i),1),City(BestSol(i+1),1)],[City(BestSol(i),2),City(BestSol(i+1),2)],'bo-');
  97.     hold on;
  98. end
  99. plot([City(BestSol(M),1),City(BestSol(1),1)],[City(BestSol(M),2),City(BestSol(1),2)],'ro-');
  100. for i=1:M
  101.     text(City(i,1)+0.5,City(i,2),num2str(i));       %标号  
  102. end
  103. text(City(BestSol(1),1),City(BestSol(1),2),'    起点');
  104. text(City(BestSol(M),1),City(BestSol(M),2),'    终点');
  105. title('最终路线图');
  106.  
  107. figure(3);
  108. t = 1:Max_iter;
  109. plot(t, Length_best, 'r', t, Length_ave, 'b--', 'LineWidth', 2);
  110. xlabel('Iteration');
  111. ylabel('Best Cost');
  112. legend('最短距离','平均距离')
  113. xlabel('迭代次数')
  114. ylabel('距离')
  115. title('各代最短距离与平均距离对比')
  116. P=zeros(1, 30);
  117.  
  118. function len=Fun(sol,City,M)  % 距离计算函数
  119. len=0;
  120. for k=1:M-1
  121.     len=len+sqrt(sum((City(sol(1,k),:)-City(sol(1,k+1),:)).^2));
  122. end
  123. len=len+sqrt(sum((City(sol(1,M),:)-City(sol(1,1),:)).^2));
  124. end

测试数据依旧用我们以前用的：

链接：https://pan.baidu.com/s/1j7omLCzx1uIJS0VdBbR-eA 
提取码：1234 

测试效果：

[1] S. Mirjalili, S. M. Mirjalili, and A. Lewis, "Grey Wolf Optimizer," Advances in Engineering Software, vol. 69, pp. 46-61, 2014.

文章知识点与官方知识档案匹配，可进一步学习相关知识
算法技能树首页概览42750 人正在系统学习中
【算法系列篇】滑动窗口
走进人工智能|深度学习 算法的创世纪