一、鲸魚算法及LSSVM簡介
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
2 LSSVM簡介
2.1 最小二乘支援向量機LSSVM基本原理
最小二乘支援向量機是支援向量機的一種改進，它是将傳統支援向量機中的不等式限制改為等式限制, 且将誤差平方和(SumSquaresError)損失函數作為訓練集的經驗損失,這樣就把解二次規劃問題轉化為求解線性方程組問題, 提高求解問題的速度和收斂精度。
常用的核函數種類：
2.2 LSSVM工具箱的使用方法
2.1 最小二乘支援向量機Matlab工具箱下載下傳連結https://www.esat.kuleuven.be/sista/lssvmlab/（毫無疑問下載下傳最新版本）
2.2 将LS-SVM檔案添加到matlan使用路徑中，便可直接使用。
具體使用步驟：
（1）導入訓練資料：load 讀取mat檔案和ASCII檔案；xlsread讀取.xls檔案；csvread讀取.csv檔案。
（2）資料預處理：效果是加快訓練速度。
方法有：歸一化處理(把每組資料都變為 - 1～ +1之間的數, 所涉及到的函數有premnmx, post mnmx, tramnmx)
标準化處理(把每組資料都化為均值為 0, 方差為 1的一組資料, 所涉及到的函數有 prestd,poatstd, trastd)
主成分分析 (進行正交處理, 減少輸入資料的維數, 所涉及到的函數有 prepca, trapca)
（3）LS-SVM lab用于函數回歸主要用到 3個函數, trainlssvm函數用來訓練建立模型, simlssvm函數用于預估模型, plotlssvm函數是 LS-SVM lab工具箱的專用繪圖函數。
（4） 參數說明：
A =csvread(′traindata. csv′);
Ptrain0=A(:, [ 1:13] );Ttrain0=A(:, [ 14:16);
[ Ptrain, meanptrain, stdptrain] = prestd(Ptrain0′);
[ Ttrain, meant , stdt] = prestd(T train0′);
Prestd()是資料歸一化函數, 其中 meanptrain是未歸一化資料之前的向量平均值 stdptrain是未歸一化資料之前的向量标準差。
gam =10;sig2=0. 5;type=′function estimation′;
LS-SVM 要求調的參數就兩個。 gam 和 sig2是最小二乘支援向量機的參數, 其中 gam 是正則化參數, 決定了适應誤差的最小化和平滑程度, sig2是 RBF 函數的參數。 在工具箱中有一個函數 gridsearch可以在一定的範圍内用來尋找最優的參數範圍。 type有兩種類型, 一種是 classfication, 用于分類, 一種是 function estimation, 用于函數回歸。
[ alpha, b] =trainlssvm({Ptrain′, Ttrain′, type, gam, sig2,′RBF_kernel′, ′preprocess′});
alpha是支援向量, b是門檻值. 。 preprocess是表明資料已經進行歸一化, 也可以是′original ′, 表明資料沒有進行歸一化, 預設時是′preprocess′。
plotlssvm ({P, T, type, gam, sig2, ′RBF _ kernel ′,′preprocess′}, {alpha, b})plotlssvm函數是 LS-SVM 工具箱特有的繪圖函數, 和 plot函數原理相近。
simlssvm函數也是 LS-SVM 工具箱的重要函數, 其中的參數如上述所示, 原理類似于神經網絡工具箱中的 sim 函數。
通過調用 trainlssvm函數和 si m lssvm 函數我們可以看到最小二乘支援向量機和神經網絡的結構有很多共同之處。
與神經網絡進行對比：
神經網絡建立的模型要比 LS-SVM 好, 但是在預估上, LS-SVM 要優于神經網絡,具有較好的泛化能力, 而且訓練速度要比神經網絡快。
二、部分源代碼
%=====================================================================
%初始化
clc
close all
clear
format long
tic
%==============================================================
%%導入資料
data=xlsread('數值.xlsx','Sheet1','A2:E41');%訓練
data1=xlsread('數值.xlsx','Sheet1','G2:J31');%測試
[row,col]=size(data);
train_x=data(:,1:col-1);
train_y=data(:,col);
test_x=data(:,1:col-1);
% test_y=data(:,col);
train_x=train_x';
train_y=train_y';
test_x=test_x';
% test_y=test_y';
%%資料歸一化
[train_x,minx,maxx, train_yy,miny,maxy] =premnmx(train_x,train_y);
test_x=tramnmx(test_x,minx,maxx);
train_x=train_x';
train_yy=train_yy';
train_y=train_y';
test_x=test_x';
% test_y=test_y';
%% 參數初始化
eps = 10^(-6);
%%定義lssvm相關參數
type='f';
kernel = 'RBF_kernel';
proprecess='proprecess';
lb=[0.01 0.02];%參數c、g的變化的下限
ub=[1000 100];%參數c、g的變化的上限
dim=2;%次元，即一個優化參數
SearchAgents_no=20; % Number of search agents
Max_iter=50; % Maximum numbef of iterations
% 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);
Positions(:,1)=ceil(rand(SearchAgents_no,1).*(ub(1)-lb(1))+lb(1));
Positions(:,2)=ceil(rand(SearchAgents_no,1).*(ub(2)-lb(2))+lb(2));
Convergence_curve=zeros(1,Max_iter);
t=0;% Loop counter
% Main loop
woa1;
%% 結果分析
plot( Convergence_curve,'LineWidth',2);
title(['鲸魚優化算法适應度曲線','(參數c1=',num2str(Leader_pos(1)),',c2=',num2str(Leader_pos(2)),',終止代數=',num2str(Max_iter),')'],'FontSize',13);
xlabel('進化代數');ylabel('誤差适應度');
bestc = Leader_pos(1);
bestg = Leader_pos(2);
gam=bestc;
sig2=bestg;
model=initlssvm(train_x,train_yy,type,gam,sig2,kernel,proprecess);%原來是顯示
model=trainlssvm(model);%原來是顯示
%求出訓練集和測試集的預測值
[train_predict_y,zt,model]=simlssvm(model,train_x);
[test_predict_y,zt,model]=simlssvm(model,test_x);
%預測資料反歸一化
train_predict=postmnmx(train_predict_y,miny,maxy);%預測輸出
test_predict=postmnmx(test_predict_y,miny,maxy);
figure
plot(train_predict,':og')
hold on
plot(train_y,'- *')
legend('預測輸出','期望輸出')
title('鲸魚優化svm網絡預測輸出','fontsize',12)
ylabel('函數輸出','fontsize',12)
xlabel('樣本','fontsize',12)
disp(['預測輸出'])
YPred_best
toc %計算時間
三、運作結果
四、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.