通过alpha, beta,delta狼的领导去寻找范围中的最优解。

个人认为理解为三个最佳猎物去吸引狼群移动更加贴切。

在初始化种群时，第一遍计算所有个体的适应度，找到最优的三个个体的位置，将这三个个体的位置定义为三个猎物的位置。后续每一次种群位置的更新都是在这三个猎物的吸引下共同作用的结果，每更新一代时重新评估所有个体适应度，如果存在个体的适应度的值比当前三个猎物更优秀，则需要对该猎物进行更新。

以求解 f(x,y) = ( x-1)^2 + ( y -3)^2为例

GWO.m代码如下

%___________________________________________________________________%
%  Grey Wolf Optimizer (GWO) source codes version 1.0               %
%                                                                   %
%  Developed in MATLAB R2011b(7.13)                                 %
%                                                                   %
%  Author and programmer: Seyedali Mirjalili                        %
%                                                                   %
%         e-Mail: ali.mirjalili@gmail.com                           %
%                 seyedali.mirjalili@griffithuni.edu.au             %
%                                                                   %
%       Homepage: http://www.alimirjalili.com                       %
%                                                                   %
%   Main paper: S. Mirjalili, S. M. Mirjalili, A. Lewis             %
%               Grey Wolf Optimizer, Advances in Engineering        %
%               Software , in press,                                %
%               DOI: 10.1016/j.advengsoft.2013.12.007               %
%                                                                   %
%___________________________________________________________________%

% 灰狼优化算法
% 返回值：
% Alpha_score：最终Alpha狼的适应度值(种群最终的最佳适应度值)
% Alpha_pos：最终Alpha狼的位置（种群最终的最优解）
% Convergence_curve ：种群适应度变化曲线（用于后期作图）
% 输入参数：
% SearchAgents_no:种群个体数量
% dim:求解问题的维度
% ub:不同维度的上边界矩阵
% lb:不同维度的下边界矩阵
function [Alpha_score,Alpha_pos,Convergence_curve]=GWO(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

% 初始化alpha, beta,delta狼，开辟空间
% 三狼的位置为1*dim矩阵，要求适应度函数最小值，所以三匹狼的初始适应度值为无穷大inf
Alpha_pos=zeros(1,dim);
Alpha_score=inf; % 将此值改为-inf求解适应度函数最大值问题

Beta_pos=zeros(1,dim);
Beta_score=inf; %change this to -inf for maximization problems

Delta_pos=zeros(1,dim);
Delta_score=inf; %change this to -inf for maximization problems

% 狼群种群个体位置初始化
Positions=initialization(SearchAgents_no,dim,ub,lb);

%用一个1*Max_iter矩阵记录最佳适应度值的变化趋势方便后续画图
Convergence_curve=zeros(1,Max_iter);

l=0;% 当前循环次数

% 主循环，每循环一次代表1次迭代，一共迭代Max_iter此
while l<Max_iter
    
    % 在该循环中依次对每一个个体的位置进行边界修正，并完成对三匹头狼的更新
    % size矩阵返回Positions矩阵的行数（种群个体数）    
    for i=1:size(Positions,1)  
        % 返回的Flag4ub与Flag4lb是1*dim的逻辑矩阵，对于满足条件的位置置1，否则为0
        Flag4ub=Positions(i,:)>ub;
        Flag4lb=Positions(i,:)<lb;
        % Flag4ub+Flag4lb：如果一个位置超出了上边界或下边界则这个位置的和为1
        % (~(Flag4ub+Flag4lb))：表示超出边界的位置为0，未超出边界的位置为1
        % Positions(i,:).*(~(Flag4ub+Flag4lb))：将超出边界的数据清零，正常数据不变
        % ub.*Flag4ub+lb.*Flag4lb：将超出上边界的数据的位置写入上边界最大值，下边界同理
        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
        
        % 分别计算每一个个体的适应度值
        fitness=fobj(Positions(i,:));
        
        % 判断新一代的群体中有没有找到更优的猎物，如果找到了更优的猎物
        % 则需要对Alpha狼、Beta狼、Delta狼(三个目前最优的猎物的位置)进行替换
        % Alphal狼的适应度函数值是最低的
        if fitness<Alpha_score
            % 如果当前个体的适应度值比Alphal狼更优秀，则当前个体成为Alphal狼
            Alpha_score=fitness;
            Alpha_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness<Beta_score 
            % 如果当前个体的适应度值大于Alphal狼但是小于此时的Beta狼，则当前个体成为Beta狼
            Beta_score=fitness;
            Beta_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score 
            % 如果当前个体比Alphal狼和Beta狼差但比当前Delta狼优秀，则当前个体成为Delta狼
            Delta_score=fitness; 
            Delta_pos=Positions(i,:);
        end
    end
    
    % 在迭代初期鼓励狼群去全局搜索更优的猎物，在迭代后期去围攻（收敛）当前最优的三个猎物
    a=2-l*((2)/Max_iter); % a逐渐递减
    
    % 更新狼群个体的位置
    % 在该循环中对种群中每一个个体都进行更新
    for i=1:size(Positions,1)
        % 在该循环中对每个个体的每个维度进行更新
        for j=1:size(Positions,2)     
            
            % r1与r2取0到1之间的随机数
            r1=rand(); 
            r2=rand(); 
            % A1的范围为[-a,a]，C1的范围为[0,2]
            A1=2*a*r1-a; % Equation (3.3)
            C1=2*r2;     % Equation (3.4)
            D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
            % 计算得出种群在Alpha狼（最优猎物位置）的影响下种群的下一代位置
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
                  
            % r1与r2取0到1之间的随机数
            r1=rand();
            r2=rand();
            % A2的范围为[-a,a]，C2的范围为[0,2]
            A2=2*a*r1-a; % Equation (3.3)
            C2=2*r2; % Equation (3.4)
            D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2
            % 计算得出种群在beta狼（次优猎物位置）的影响下种群的下一代位置
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2       
            
            % r1与r2取0到1之间的随机数
            r1=rand();
            r2=rand(); 
            % A3的范围为[-a,a]，C3的范围为[0,2]
            A3=2*a*r1-a; % Equation (3.3)
            C3=2*r2; % Equation (3.4)
            D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3
            % 计算得出种群在delta狼（次次优猎物位置）的影响下种群的下一代位置
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3             
            
            % 最终狼群下一代的真正位置是在三个最优解的吸引作用下的共同结果，即对上述位置求平均值
            Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)
            
        end
    end
    % 准备进行下一次迭代，记录下当前的Alpha狼（当前最优解）的适应度值（表示整个种群当前最优适应度）
    l=l+1;    
    Convergence_curve(l)=Alpha_score;
end

%% 初始化函数
function Positions=initialization(SearchAgents_no,dim,ub,lb)

% size用于获取矩阵或数组的维度信息
% size(ub,2)表示返回矩阵ub的列数；size(ub,1)表示返回矩阵ub的行数
Boundary_no= size(ub,2); % 边界的数量

% 如果所有维度上的变量的边界都是相同的或者用户只输入了一个数字
if Boundary_no==1
    Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end

% 如果对于不同维度的变量的边界也不同（边界数量>1）
if Boundary_no>1
    %分别处理两个维度的边界
    for i=1:dim
        ub_i=ub(i);
        lb_i=lb(i);
        Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
    end
end

适应度函数文件TestFun2D.m代码：

function [fitvalue]=TestFun2D(x)

fitvalue = ( x(1) -1)^2 + ( x(2) -3)^2;

end

主函数文件TestMain.m代码：

ObjFun=@TestFun2D;
%function [Alpha_score,Alpha_pos,Convergence_curve]=
% GWO(SearchAgents_no,Max_iter,lb,ub,dim,fobj)
[Alpha_score,Alpha_pos,Convergence_curve] = GWO(50,50,[0,3],[2,4],2,ObjFun);

display(Alpha_pos);
display(Alpha_score);