一、鯨魚演算法及LSTM簡介

1 鯨魚優化演算法（Whale Optimization Algorithm，WOA）簡介
鯨魚優化演算法(WOA)，該演算法模擬了座頭鯨的社會行為，並引入了氣泡網狩獵策略。

1.1 靈感
鯨魚被認為是世界上最大的哺乳動物。一頭成年鯨可以長達 30 米，重 180 噸。這種巨型哺乳動物有 7 種不同的主要物種，如虎鯨，小鬚鯨，鰮鯨，座頭鯨，露脊鯨，長鬚鯨和藍鯨等。鯨通常被認為是食肉動物，它們從不睡覺，因為它們必須到海洋表面進行呼吸，但事實上，鯨魚有一半的大腦都處於睡眠狀態。
鯨魚在大腦的某些區域有與人類相似的細胞，這些細胞被稱為紡錘形細胞（spindle cells）。這些細胞負責人類的判斷、情感和社會行為。換句話說，紡錘形細胞使我們人類有別於其他生物。鯨魚的這些細胞數量是成年人的兩倍，這是它們具有高度智慧和更富情感的主要原因。已經證明，鯨魚可以像人類一樣思考、學習、判斷、交流，甚至變得情緒化，但顯然，這都只是在一個很低的智慧水平上。據觀察，鯨魚(主要是虎鯨)也能發展自己的方言。
另一個有趣的點是關於鯨魚的社會行為，它們可獨居也可群居，但我們觀察到的大多數仍然是群居。它們中的一些物種(例如虎鯨)可以在整個生命週期中生活在一個家族中。最大的鬚鯨之一是座頭鯨，一頭成年座頭鯨幾乎和一輛校車一樣大。它們最喜歡的獵物是磷蝦和小魚群。圖1顯示的就是這種哺乳動物。

圖1 座頭鯨的氣泡網進食行為
關於座頭鯨最有趣的事情是它們特殊的捕獵方法了。這種覓食行為被稱為氣泡網覓食法（bubble-net feeding method）。座頭鯨喜歡在接近海面的地方捕食磷蝦或小魚。據觀察，這種覓食是通過在圓形或類似數字“9”形路徑上製造獨特的氣泡來完成的，如圖 1 所示。在 2011 年之前，這一行為僅僅是基於海面觀測的。然而，有研究者利用標籤感測器研究了這種行為。他們捕獲了9頭座頭鯨身上300個由標籤得到的氣泡網進食事件。他們發現了兩種與氣泡有關的策略，並將它們命名為上升螺旋（upward-spirals）和雙螺旋（doubleloops）。在前一種策略中，座頭鯨會潛到水下 12 米左右，然後開始在獵物周圍製造一個螺旋形的泡泡，並遊向水面；後一種策略包括三個不同的階段:珊瑚迴圈,用尾葉拍打水面以及捕獲迴圈。這裡不展開詳細描述。
但是氣泡網捕食是隻有座頭鯨獨有的一種特殊行為，而鯨魚優化演算法就是模擬了螺旋氣泡網進食策略達到優化的目的。

1.2 數學建模和優化演算法
1.2.1 包圍捕食（Encircling prey）
座頭鯨可以識別獵物的位置並將其包圍。由於最優設計在搜尋空間中的位置不是先驗已知的，WOA 演算法假設當前的最佳候選解是目標獵物或接近最優解。在定義了最佳搜尋代理之後，其他搜尋代理將因此嘗試向最佳搜尋代理更新它們的位置。這種行為由下列方程表示：




  圖 2a 描述了等式（2）針對2D問題的基本原理，搜尋代理的位置( X , Y )可以根據當前最優解的位置( X ∗ , Y ∗ )進行更新，通過調整向量 A ⃗ 和C的值，可以找到相對於當前位置下一時刻最優代理附近的不同地方。在 3D 空間中搜索代理可能的更新位置如圖 2b。通過定義隨機向量 r ，可以到達圖 2 中所示關鍵點之間的搜尋空間內任何位置，因此等式（2）允許任何搜尋代理在當前最優解的鄰域內更新其位置，從而模擬了鯨魚的包圍捕食。相似的概念也可以擴充套件到 n 維搜尋空間。注意圖2中的兩幅圖均是在a=1和C=1情況下的。


圖2 2D和3D位置向量及其可能的下一個位置
1.2.2 氣泡網攻擊方式（Bubble-net attacking method）（利用階段）

共設計了兩種方法來對座頭鯨的氣泡網行為進行建模：
收縮包圍機制：通過降低式（3）中 a 的值實現。注意 A的波動範圍也通過 a降低，換句話說，A 是一個區間[-a,a]內的隨機值，a 隨著迭代進行從 2 降為 0。設定 A中的隨機值在[-1,1]之間，搜尋代理的新位置可以定義為代理原始位置與當前最優代理位置之間的任意位置。圖 3a 顯示了 2D 空間中當 0 ≤ A ≤ 1 0 時從 ( X , Y )靠近 ( X ∗ , Y ∗ ) 所有可能的位置。這種機制本質上就是包圍捕食。
螺旋更新位置。如圖 3b，該方法首先計算鯨魚位置 ( X , Y ) 與獵物位置 ( X ∗ , Y ∗ ) 之間的距離，然後在鯨魚與獵物位置之間建立一個螺旋等式，來模仿座頭鯨的螺旋狀移動：



(a)收縮包圍機制

(b)螺旋更新位置
圖3 WOA中實現的氣泡網搜尋機制
值得注意的是，座頭鯨在一個不斷縮小的圓圈內繞著獵物遊動，同時沿著螺旋形路徑遊動。為了對這種同時發生的行為進行建模，假設有 50%的可能性在收縮包圍機制和螺旋模型之間進行選擇，以便在優化過程中更新鯨魚的位置，數學模型如下：

其中 p pp 為[0,1]之間的隨機數。
1.2.3搜尋獵物（Search for prey）(exploration phase)
除了泡泡網方法，座頭鯨還會隨機尋找獵物，同樣基於可變 A向量，事實上，座頭鯨會根據彼此的位置進行隨機搜尋，因此使用隨機值大於1或小於-1的 A ⃗ 來迫使搜尋代理遠離參考鯨魚。與利用階段相反，這裡將根據隨機選擇的搜尋代理來更新搜尋代理在探索階段的位置，而不是根據目前為止最優的搜尋代理。該機制和 ∣ A ⃗ ∣ > 1 強調了探索，並允許WOA演算法執行全域性搜尋。數學模型如下：


其中 X → r a n d 為從當前種群中選擇的隨機位置向量（表示一頭隨機鯨魚）。
特定解附近滿足 A ⃗ > 1的一些可能解如圖 4 所示。

圖4 WOA中的探索機制(X*是一個隨機選擇的搜尋代理)
WOA演算法首先隨機初始化一組解,在每次迭代中，搜尋代理根據隨機選擇的搜尋代理或到目前為止獲得的最優解更新它們的位置。將 a aa 引數由 2 隨迭代次數降為 0，從而由探索逐步到利用。當 ∣ A ⃗ ∣ > 1 時選擇隨機搜尋代理，∣ A ⃗ ∣ < 1 時選擇最優解更新搜尋代理位置。根據 p pp 的值，WOA可以在螺旋運動和圓環運動之間進行切換。最後，通過滿足終止準則來終止WOA演算法。WOA演算法的虛擬碼如圖5所示。

圖5 WOA演算法虛擬碼
1.3 程式碼分析
只要明白了原理的基本流程，其實程式碼就沒有說明困難了，咱們主要介紹一下如何實現上述分析的幾個重要原理，所要優化的問題的是三十個數的平方和最小
（1）引數初始化。初始時主要設定代理數量和最大迭代次數即可，其他演算法相關的引數因為和當前迭代次數相關，需要在迭代中設定。

SearchAgents_no=30; % 搜尋代理數量
Max_iteration=500; % 最大迭代次數
``
**（2） 種群初始化**。隨機初始化所有代理各個維度上的位置值，需要保證在取值範圍內。

```c
Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;

（3）種群評估。評估種群中每個代理的目標值，如有某個代理由於當前最優解，則將其設為最優解。

for i=1:size(Positions,1)
    % 計算每個代理的目標值
    fitness=fobj(Positions(i,:));
    % 更新最優解
    if fitness<Leader_score % 如果是最大化問題，這裡就是">"
        Leader_score=fitness; 
        Leader_pos=Positions(i,:);
    end
end

（4）設定和迭代次數相關的演算法引數。

a=2-t*((2)/Max_iter); % 等式（3）中a隨迭代次數從2線性下降至0 
%a2從-1線性下降至-2，計算l時會用到
a2=-1+t*((-1)/Max_iter);

（5）對每個代理的每一維度進行位置更新。

% Update the Position of search agents 
for i=1:size(Positions,1)
    r1=rand(); % r1為[0,1]之間的隨機數
    r2=rand(); % r2為[0,1]之間的隨機數

    A=2*a*r1-a;  % 等式（3）
    C=2*r2;      % 等式（4）

    b=1;               %  等式（5）中的常數b
    l=(a2-1)*rand+1;   %  等式（5）中的隨機數l
    p = rand();        %  等式（6）中的概率p
    for j=1:size(Positions,2)
        if p<0.5   
            if abs(A)>=1
                rand_leader_index = floor(SearchAgents_no*rand()+1);
                X_rand = Positions(rand_leader_index, :);
                D_X_rand=abs(C*X_rand(j)-Positions(i,j)); % 等式（7）
                Positions(i,j)=X_rand(j)-A*D_X_rand;      % 等式（8）

            elseif abs(A)<1
                D_Leader=abs(C*Leader_pos(j)-Positions(i,j)); % 等式（1）
                Positions(i,j)=Leader_pos(j)-A*D_Leader;      % 等式（2）
            end
        elseif p>=0.5

            distance2Leader=abs(Leader_pos(j)-Positions(i,j));
            % 等式（5）
            Positions(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j);
        end
    end
end

二、部分原始碼

clear all 
clc

SearchAgents_no=30; % Number of search agents

Function_name='F1'; % Name of the test function that can be from F1 to F23 (Table 1,2,3 in the paper)

Max_iteration=500; % Maximum numbef of iterations

% Load details of the selected benchmark function
[lb,ub,dim,fobj]=Get_Functions_details(Function_name);

[Best_score,Best_pos,WOA_cg_curve]=WOA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);

figure('Position',[269   240   660   290])
%Draw search space
subplot(1,2,1);
func_plot(Function_name);
title('Parameter space')
xlabel('x_1');
ylabel('x_2');
zlabel([Function_name,'( x_1 , x_2 )'])

%Draw objective space
subplot(1,2,2);
semilogy(WOA_cg_curve,'Color','r')
title('Objective space')
xlabel('Iteration');
ylabel('Best score obtained so far');

axis tight
grid on
box on
legend('WOA')

display(['The best solution obtained by WOA is : ', num2str(Best_pos)]);
display(['The best optimal value of the objective funciton found by WOA is : ', num2str(Best_score)]);

    function [Leader_score,Leader_pos,Convergence_curve]=WOA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

% initialize position vector and score for the leader
Leader_pos=zeros(1,dim);
Leader_score=inf; %change this to -inf for maximization problems


%Initialize the positions of search agents
Positions=initialization(SearchAgents_no,dim,ub,lb);

Convergence_curve=zeros(1,Max_iter);

t=0;% Loop counter

% Main loop
while t<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 the leader
        if fitness<Leader_score % Change this to > for maximization problem
            Leader_score=fitness; % Update alpha
            Leader_pos=Positions(i,:);
        end
        
    end
    
    a=2-t*((2)/Max_iter); % a decreases linearly fron 2 to 0 in Eq. (2.3)
    
    % a2 linearly dicreases from -1 to -2 to calculate t in Eq. (3.12)
    a2=-1+t*((-1)/Max_iter);
    
    % Update the Position of search agents 
    for i=1:size(Positions,1)
        r1=rand(); % r1 is a random number in [0,1]
        r2=rand(); % r2 is a random number in [0,1]
        
        A=2*a*r1-a;  % Eq. (2.3) in the paper
          % Eq. (2.4) in the paper
        
        
                %  parameters in Eq. (2.5)
        l=(a2-1)*rand+1;   %  parameters in Eq. (2.5)
        
        p = rand();        % p in Eq. (2.6)
        
        for j=1:size(Positions,2)
            
            if p<0.5   
                if abs(A)>=1
                    rand_leader_index = floor(SearchAgents_no*rand()+1);
                    X_rand = Positions(rand_leader_index, :);
                    D_X_rand=abs(C*X_rand(j)-Positions(i,j)); % Eq. (2.7)
                    Positions(i,j)=X_rand(j)-A*D_X_rand;      % Eq. (2.8)
                    
                elseif abs(A)<1
                    D_Leader=abs(C*Leader_pos(j)-Positions(i,j)); % Eq. (2.1)
                    Positions(i,j)=Leader_pos(j)-A*D_Leader;      % Eq. (2.2)
                end
                
            elseif p>=0.5
              
              
                % Eq. (2.5)
                Positions(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j);
                
            end
            
        end
    end    

三、執行結果

四、matlab版本及參考文獻

1 matlab版本
2014a

2 參考文獻
[1] 包子陽,餘繼周,楊杉.智慧優化演算法及其MATLAB例項（第2版）[M].電子工業出版社，2016.
[2]張巖,吳水根.MATLAB優化演算法原始碼[M].清華大學出版社，2017.
[3]周品.MATLAB 神經網路設計與應用[M].清華大學出版社，2013.
[4]陳明.MATLAB神經網路原理與例項精解[M].清華大學出版社，2013.
[5]方清城.MATLAB R2016a神經網路設計與應用28個案例分析[M].清華大學出版社，2018.