基本灰狼算法的MATLAB實現

灰狼算法：

灰狼優化算法（Grey Wolf Optimizer，GWO）由澳大利亞格里菲斯大學學者 Mirjalili 等人於2014年提出來的一種群智能優化算法。該算法受到了灰狼捕食獵物活動的啟發而開發的一種優化搜索方法，它具有較強的收斂性能、參數少、易實現等特點。近年來受到了學者的廣泛關注，它己被成功地應用到了車間調度、參數優化、圖像分類等領域中。

理解算法思想方面主要參考了原始論文，以及以下的博客

https://www.pianshen.com/article/35241673990/

https://blog.csdn.net/haha0332/article/details/88805910

https://www.it610.com/article/1288128297732976640.htm

https://blog.csdn.net/hba646333407/article/details/107297509

我個人在學習過程中，主要的困惑在於原論文中位置向量xp的獲取：

 

 

                                             第一種                                                                                                            第二種

因為如果位置已知，那還為什么要優化，覺得很奇怪。

后面查閱第一篇博客才理解，兩種定義式中，第一種（左圖）是從想法出發（如果知道位置向量就好了），第二種（右圖）其實是回到現實，現實中不可能知道獵物的位置向量，因此這里引出了灰狼的一個核心想法：將xp用三只頭狼的信息估計出來，等效代替，理解了這一點，其他的細節都清楚了。

之后學習源碼，由於筆者的MATLAB基礎不深厚，因此順便總結了源代碼的一些MATLAB用法：

rand(m,n)%返回m行n列的矩陣，矩陣中每一個數都是（0，1）間的隨機數。

A(i,:)%返回矩陣A的第i行

A(:,i)%返回矩陣A的第i列

@f %函數句柄，表示調用此函數

size(A,1)%返回矩陣A的行數

size(A,2)%返回矩陣A的列數

axis tight%使邊框貼近曲線
grid on%畫網格線
box on%圖形邊框
legend('GWO') %標注字符串到圖中

另外一些值得總結的地方在於函數實現過程中的一些細節

1.遍歷Position矩陣的行，並確定當前的α、β、γ狼

2.初始化以及調用時行列中代表的含義

源代碼以及二者的示意圖如下

% Main loop
while l<Max_iter
    for i=1:size(Positions,1)  
        
       % Return back the search agents that go beyond the boundaries of the search space
        Flag4ub=Positions(i,:)>ub;
        Flag4lb=Positions(i,:)<lb;
        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
        
        % Calculate objective function for each search agent
        fitness=fobj(Positions(i,:));
        
        % Update Alpha, Beta, and Delta
        if fitness<Alpha_score 
            Alpha_score=fitness; % Update alpha
            Alpha_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness<Beta_score 
            Beta_score=fitness; % Update beta
            Beta_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score 
            Delta_score=fitness; % Update delta
            Delta_pos=Positions(i,:);
        end
    end
    
    
    a=2-l*((2)/Max_iter); % a decreases linearly fron 2 to 0
    
    % Update the Position of search agents including omegas
    for i=1:size(Positions,1)
        for j=1:size(Positions,2)     
                       
            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; % Equation (3.3)
            C1=2*r2; % Equation (3.4)
            
            D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
                       
            r1=rand();
            r2=rand();
            
            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
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2       
            
            r1=rand();
            r2=rand(); 
            
            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
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3             
            
            Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)
            
        end
    end
    l=l+1;    
    Convergence_curve(l)=Alpha_score;
end