Python&Matlab實現灰狼優化演演算法的範例程式碼

2022-03-21 16:01:01

1 灰狼優化演演算法基本思想

灰狼優化演演算法是一種群智慧優化演演算法，它的獨特之處在於一小部分擁有絕對話語權的灰狼帶領一群灰狼向獵物前進。在瞭解灰狼優化演演算法的特點之前，我們有必要了解灰狼群中的等級制度。

灰狼群一般分為4個等級：處於第一等級的灰狼用α表示，處於第二階級的灰狼用β表示，處於第三階段的灰狼用δ表示，處於第四等級的灰狼用ω表示。按照上述等級的劃分，灰狼α對灰狼β、δ和ω有絕對的支配權；灰狼ω對灰狼δ和ω有絕對的支配權；灰狼δ對灰狼ω有絕對的支配權。

2 灰狼捕食獵物過程

GWO 優化過程包含了灰狼的社會等級分層、跟蹤、包圍和攻擊獵物等步驟，其步驟具體情況如下所示。

2.1 社會等級分層

當設計 GWO 時，首先需構建灰狼社會等級層次模型。計算種群每個個體的適應度，將狼群中適應度最好的三匹灰狼依次標記為α、β、δ而剩下的灰狼標記為ω 。也就是說，灰狼群體中的社會等級從高往低排列依次為α、β、δ及 ω。GWO 的優化過程主要由每代種群中的最好三個解（即 α、β、δ）來指導完成。

2.2 包圍獵物

灰狼群體通過以下幾個公式逐漸接近幷包圍獵物:

式中，t是當前的迭代代數，A和C是係數向量，Xp和X分別是獵物的位置向量和灰狼的位置向量。A和C的計算公式如下:

式中，a是收斂因子，隨著迭代次數從2線性減小到0，r1和r 2服從[ 0，1]之間的均勻分佈。

2.3 狩獵

狼群中其他灰狼個體Xi根據α、β和百的位置Xa、XB和Xo來更新各自的位置：

式中,Da，Dβ和D6分別表示a，β和5與其他個體間的距離;Xa，Xβ和X6分別代表a，β和5的當前位置;C1，C2，C3是隨機向量，X是當前灰狼的位置。

灰狼個體的位置更新公式如下:

2.4 攻擊獵物

構建攻擊獵物模型的過程中，根據2)中的公式，ａ值的減少會引起 A 的值也隨之波動。換句話說，A 是一個在區間［-a，a］（備註：原作者的第一篇論文裡這裡是[-2a,2a]，後面論文裡糾正為[-a,a]）上的隨機向量，其中ａ在迭代過程中呈線性下降。當 A 在［-1，１］區間上時，則捜索代理（Search Agent）的下一時刻位置可以在當前灰狼與獵物之間的任何位置上。

2.5 尋找獵物

灰狼主要依賴α、β、δ的資訊來尋找獵物。它們開始分散地去搜尋獵物位置資訊，然後集中起來攻擊獵物。對於分散模型的建立，通過｜A｜＞１使其捜索代理遠離獵物，這種搜尋方式使 GWO 能進行全域性搜尋。GWO 演演算法中的另一個搜尋係數是Ｃ。從2.2中的公式可知，Ｃ向量是在區間範圍［０，２］上的隨機值構成的向量，此係數為獵物提供了隨機權重，以便増加（｜Ｃ｜＞１）或減少（｜Ｃ｜＜１）。這有助於 GWO 在優化過程中展示出隨機搜尋行為，以避免演演算法陷入區域性最優。值得注意的是，Ｃ並不是線性下降的，Ｃ在迭代過程中是隨機值，該係數有利於演演算法跳出區域性，特別是演演算法在迭代的後期顯得尤為重要。

3 實現步驟及程式框圖

3.1 步驟

Step1:種群初始化:包括種群數量N，最大迭代次數Maxlter，調控引數a,A，C.

Step2:根據變數的上下界來隨機初始化灰狼個體的位置X。

Step3:計算每一頭狼的適應度值，並將種群中適應度值最優的狼的位置資訊儲存X_α，將種群中適應度值次優的狼的位置資訊儲存為X_β，將種群中適應度第三優的灰狼的位置資訊儲存為X_γ。

Step4:更新灰狼個體X的位置。

step5:更新引數a，A和C。

Step6:計算每一頭灰狼的適應度值，並更新三匹頭狼的最優位置。

Step7:判斷是否到達最大迭代次數Maxlter，若滿足則演演算法停止並返回Xa的值作為最終得到的最優解，否則轉到Step4。

3.2 程式框圖

4 Python程式碼實現

 
#=======匯入線管庫======
import random
import numpy
 
#完整程式碼見微信公眾號：電力系統與演演算法之美
#輸入關鍵字：灰狼演演算法
 
def GWO(objf, lb, ub, dim, SearchAgents_no, Max_iter):
 
    #===初始化 alpha, beta, and delta_pos=======
    Alpha_pos = numpy.zeros(dim)  # 位置.形成30的列表
    Alpha_score = float("inf")  # 這個是表示「正負無窮」,所有數都比 +inf 小；正無窮：float("inf"); 負無窮：float("-inf")
 
    Beta_pos = numpy.zeros(dim)
    Beta_score = float("inf")
 
    Delta_pos = numpy.zeros(dim)
    Delta_score = float("inf")  # float() 函數用於將整數和字串轉換成浮點數。
 
    #====list列表型別=============
    if not isinstance(lb, list):  # 作用：來判斷一個物件是否是一個已知的型別。 其第一個引數（object）為物件，第二個引數（type）為型別名，若物件的型別與引數二的型別相同則返回True
        lb = [lb] * dim  # 生成[100，100，.....100]30個
    if not isinstance(ub, list):
        ub = [ub] * dim
 
    #========初始化所有狼的位置===================
    Positions = numpy.zeros((SearchAgents_no, dim))
    for i in range(dim):  # 形成5*30個數[-100，100)以內
        Positions[:, i] = numpy.random.uniform(0, 1, SearchAgents_no) * (ub[i] - lb[i]) + lb[
            i]  # 形成[5個0-1的數]*100-（-100）-100
    Convergence_curve = numpy.zeros(Max_iter)
 
    #========迭代尋優=====================
    for l in range(0, Max_iter):  # 迭代1000
        for i in range(0, SearchAgents_no):  # 5
            #====返回超出搜尋空間邊界的搜尋代理====
            for j in range(dim):  # 30
                Positions[i, j] = numpy.clip(Positions[i, j], lb[j], ub[
                    j])  # clip這個函數將將陣列中的元素限制在a_min(-100), a_max(100)之間，大於a_max的就使得它等於 a_max，小於a_min,的就使得它等於a_min。
 
        
 
        #===========以上的迴圈裡，Alpha、Beta、Delta===========
        a = 2 - l * ((2) / Max_iter);  #   a從2線性減少到0
 
        for i in range(0, SearchAgents_no):
            for j in range(0, dim):
                r1 = random.random()  # r1 is a random number in [0,1]主要生成一個0-1的隨機浮點數。
                r2 = random.random()  # r2 is a random number in [0,1]
 
                A1 = 2 * a * r1 - a;  # Equation (3.3)
                C1 = 2 * r2;  # Equation (3.4)
                # D_alpha表示候選狼與Alpha狼的距離
                D_alpha = abs(C1 * Alpha_pos[j] - Positions[
                    i, j]);  # abs() 函數返回數位的絕對值。Alpha_pos[j]表示Alpha位置，Positions[i,j])候選灰狼所在位置
                X1 = Alpha_pos[j] - A1 * D_alpha;  # X1表示根據alpha得出的下一代灰狼位置向量
 
                r1 = random.random()
                r2 = random.random()
 
                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 = random.random()
                r2 = random.random()
 
                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  # 候選狼的位置更新為根據Alpha、Beta、Delta得出的下一代灰狼地址。
 
        Convergence_curve[l] = Alpha_score;
 
        if (l % 1 == 0):
            print(['迭代次數為' + str(l) + ' 的迭代結果' + str(Alpha_score)]);  # 每一次的迭代結果
 
#========函數==========
def F1(x):
    s=numpy.sum(x**2);
    return s
 
#===========主程式================
func_details = ['F1', -100, 100, 30]
function_name = func_details[0]
Max_iter = 1000#迭代次數
lb = -100#下界
ub = 100#上屆
dim = 30#狼的尋值範圍
SearchAgents_no = 5#尋值的狼的數量
x = GWO(F1, lb, ub, dim, SearchAgents_no, Max_iter)

5 Matlab實現

% 主程式 GWO
clear
close all
clc
 
%%完整程式碼見微信公眾號：電力系統與演演算法之美
 
%輸入關鍵字：灰狼演演算法
 
SearchAgents_no = 30 ; %　種群規模
dim = 10 ; % 粒子維度
Max_iter = 1000 ; %　迭代次數
ub = 5 ;
lb = -5 ;
 
%% 初始化三匹頭狼的位置
Alpha_pos=zeros(1,dim);
Alpha_score=inf; 
 
Beta_pos=zeros(1,dim);
Beta_score=inf; 
 
Delta_pos=zeros(1,dim);
Delta_score=inf; 
 
 
 
Convergence_curve = zeros(Max_iter,1);
 
%% 開始迴圈
for l=1:Max_iter
    for i=1:size(Positions,1)  
        
       %% 返回超出搜尋空間邊界的搜尋代理
        Flag4ub=Positions(i,:)>ub;
        Flag4lb=Positions(i,:)<lb;
        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
        
        %% 計算每個搜尋代理的目標函數
        fitness=sum(Positions(i,:).^2);
        
        %% 更新 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
    
    %% 更新搜尋代理的位置，包括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
  
    Convergence_curve(l)=Alpha_score;
    disp(['Iteration = ' num2str(l)  ', Evaluations = ' num2str(Alpha_score)]);
 
end
%========視覺化==============
figure('unit','normalize','Position',[0.3,0.35,0.4,0.35],'color',[1 1 1],'toolbar','none')
%% 目標空間
subplot(1,2,1);
x = -5:0.1:5;y=x;
L=length(x);
f=zeros(L,L);
for i=1:L
    for j=1:L
       f(i,j) = x(i)^2+y(j)^2;
    end
end
surfc(x,y,f,'LineStyle','none');
xlabel('x_1');
ylabel('x_2');
zlabel('F')
title('Objective space')
%% 狼群演演算法 
subplot(1,2,2);
semilogy(Convergence_curve,'Color','r','linewidth',1.5)
title('Convergence_curve')
xlabel('Iteration');
ylabel('Best score obtained so far');
 
axis tight
grid on
box on
legend('GWO')
display(['The best solution obtained by GWO is : ', num2str(Alpha_pos)]);
display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Alpha_score)]);