​

一、 BP神經網路預測演算法簡介

說明：1.1節主要是概括和幫助理解考慮影響因素的BP神經網路演算法原理，即常規的BP模型訓練原理講解（可根據自身掌握的知識是否跳過）。1.2節開始講基於歷史值影響的BP神經網路預測模型。

使用BP神經網路進行預測時，從考慮的輸入指標角度，主要有兩類模型：

1.1 受相關指標影響的BP神經網路演算法原理

如圖一所示，使用MATLAB的newff函式訓練BP時，可以看到大部分情況是三層的神經網路（即輸入層，隱含層，輸出層）。這裡幫助理解下神經網路原理：
1）輸入層：相當於人的五官，五官獲取外部資訊，對應神經網路模型input埠接收輸入資料的過程。
2）隱含層：對應人的大腦，大腦對五官傳遞來的資料進行分析和思考，神經網路的隱含層hidden Layer對輸入層傳來的資料x進行對映，簡單理解為一個公式hiddenLayer_output=F(w*x+b)。其中，w、b叫做權重、閾值引數，F()為對映規則，也叫啟用函式，hiddenLayer_output是隱含層對於傳來的資料對映的輸出值。換句話說，隱含層對於輸入的影響因素資料x進行了對映，產生了對映值。
3）輸出層：可以對應為人的四肢，大腦對五官傳來的資訊經過思考（隱含層對映）之後，再控制四肢執行動作（向外部作出響應）。類似地，BP神經網路的輸出層對hiddenLayer_output再次進行對映，outputLayer_output=w *hiddenLayer_output+b。其中，w、b為權重、閾值引數，outputLayer_output是神經網路輸出層的輸出值（也叫模擬值、預測值）（理解為，人腦對外的執行動作，比如嬰兒拍打桌子）。
4）梯度下降演算法：通過計算outputLayer_output和神經網路模型傳入的y值之間的偏差，使用演算法來相應調整權重和閾值等引數。這個過程，可以理解為嬰兒拍打桌子，打偏了，根據偏離的距離遠近，來調整身體使得再次揮動的胳膊不斷靠近桌子，最終打中。

再舉個例子來加深理解：

圖一所示BP神經網路，具備輸入層、隱含層和輸出層。BP是如何通過這三層結構來實現輸出層的輸出值outputLayer_output，不斷逼近給定的y值，從而訓練得到一個精準的模型的呢？

從圖中串起來的埠，可以想到一個過程：坐地鐵，將圖一想象為一條地鐵線路。王某某坐地鐵回家的一天：在input起點站上車，中途經過了很多站（hiddenLayer），然後發現坐過頭了（outputLayer對應現在的位置），那麼王某某將會根據現在的位置離家（目標Target）的距離（誤差Error），返回到中途的地鐵站（hiddenLayer）重新坐地鐵（誤差反向傳遞，使用梯度下降演算法更新w和b），如果王某某又一次發生失誤，那麼將再次進行這個調整的過程。

從在嬰兒拍打桌子和王某某坐地鐵的例子中，思考問題：BP的完整訓練，需要先傳入資料給input，再經過隱含層的對映，輸出層得到BP模擬值，根據模擬值與目標值的誤差，來調整引數，使得模擬值不斷逼近目標值。比如（1）嬰兒受到了外界的干擾因素（x），從而作出反應拍桌（predict），大腦不斷的調整胳膊位置，控制四肢拍準（y、Target）。（2）王某某上車點（x），過站點（predict），不斷返回中途站來調整位置，到家（y、Target）。

在這些環節中，涉及了影響因素資料x，目標值資料y（Target）。根據x，y，使用BP演算法來尋求x與y之間存在的規律，實現由x來對映逼近y，這就是BP神經網路演算法的作用。再多說一句，上述講的過程，都是BP模型訓練，那麼最終得到的模型雖然訓練準確，但是找到的規律（bp network）是否準確與可靠呢。於是，我們再給x1到訓練好的bp network中，得到相應的BP輸出值（預測值）predict1，通過作圖，計算Mse，Mape，R方等指標，來對比predict1和y1的接近程度，就可以知道模型是否預測準確。這是BP模型的測試過程，即實現對資料的預測，並且對比實際值檢驗預測是否準確。
​
圖一 3層BP神經網路結構圖

1.2 基於歷史值影響的BP神經網路

以電力負荷預測問題為例，進行兩種模型的區分。在預測某個時間段內的電力負荷時：

一種做法，是考慮t時刻的氣候因素指標，比如該時刻的空氣溼度x1，溫度x2，以及節假日x3等的影響，對t時刻的負荷值進行預測。這是前面1.1所說的模型。

另一種做法，是認為電力負荷值的變化，與時間相關，比如認為t-1，t-2，t-3時刻的電力負荷值與t時刻的負荷值有關係，即滿足公式y(t)=F(y(t-1),y(t-2),y(t-3))。採用BP神經網路進行訓練模型時，則輸入到神經網路的影響因素值為歷史負荷值y(t-1),y(t-2),y(t-3)，特別地，3叫做自迴歸階數或者延遲。給到神經網路中的目標輸出值為y(t)。

二、人工蜂群演算法

受到蜜蜂群體的有組織的覓食過程的啟發，Karaboga提出了模擬蜜蜂群體覓食過程的人工蜂群(Artificial Bee Colony)演算法用於解決多維度多峰谷的優化問題。該演算法創始之初被用來尋找Sphere、Rosenbrock和Rastrigin函式的最小值。 首先對蜜蜂基於搖擺舞進行覓食的過程特徵進行介紹。在圖1中，存在兩個已發現的食物源A和B。初始時，潛在工蜂以非僱傭蜂的身份進行搜尋。它並不知道蜂房附近的任何蜜源的資訊。因此，它有以下兩個可能的選擇： （1）成為一個偵察蜂，秉著自身潛在動力或外在因素自發的搜尋蜂房附近的區域(見圖1中的S)； （2）在觀看擺尾舞后，成為一個被招募者，並開始搜尋蜜源(見圖1中的R)。 在定位蜜源之後，該蜜蜂能夠利用自身的能力來記住食物源的位置，並立刻對它進行探索。該蜜蜂現在成為了一個僱傭蜂。僱傭蜂採到蜂蜜後，從蜜源處返回蜂房並將蜂蜜解除安裝到蜜室中。在解除安裝完蜂蜜後，僱傭蜂有下列三個選擇： （1）放棄已經採集過的蜜源，成為一個受其他搖尾舞招募的跟隨者(UF)。 （2）施展搖尾舞技，招募蜂房內的同伴，再次回到原先採集過的食物源(EF1)。 （3）不招募其它的蜜蜂，繼續探索採集過的食物源(EF2)。​​

圖1 蜜蜂覓食行為圖

二、演算法流程

人工蜂群演算法由連續的四個階段組成，分別是初始化階段、引領(僱傭)蜂階段、跟隨蜂階端和偵察蜂階段。 人工蜂群演算法中將人工蜂群分為引領蜂、跟隨蜂和偵察蜂三類，每一次搜尋過程中，引領蜂和跟隨蜂是先後開採食物源，即尋找最優解，而偵察蜂是觀察是否陷入區域性最優，若陷入區域性最優則隨機地搜尋其它可能的食物源。每個食物源代表問題一個可能解，食物源的花蜜量對應相應解的質量（適應度值f i t fitfit）。 ABC演算法流程圖如圖2所示。​​

圖2 ABC演算法流程圖

1、初始化階段

​​

2、引領蜂階段

​​3、跟隨蜂階段

​​

4、偵察蜂階段

​​

5、食物源

​

三、部分程式碼

%

clc;
clear;
close all;

%% Problem Definition

CostFunction=@(x) Sphere(x);        % Cost Function

nVar=5;             % Number of Decision Variables

VarSize=[1 nVar];   % Decision Variables Matrix Size

VarMin=-10;         % Decision Variables Lower Bound
VarMax= 10;         % Decision Variables Upper Bound

%% ABC Settings

MaxIt=200;              % Maximum Number of Iterations

nPop=100;               % Population Size (Colony Size)

nOnlooker=nPop;         % Number of Onlooker Bees

L=round(0.6*nVar*nPop); % Abandonment Limit Parameter (Trial Limit)

a=1;                    % Acceleration Coefficient Upper Bound

%% Initialization

% Empty Bee Structure
empty_bee.Position=[];
empty_bee.Cost=[];

% Initialize Population Array
pop=repmat(empty_bee,nPop,1);

% Initialize Best Solution Ever Found
BestSol.Cost=inf;

% Create Initial Population
for i=1:nPop
    pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
    pop(i).Cost=CostFunction(pop(i).Position);
    if pop(i).Cost<=BestSol.Cost
        BestSol=pop(i);
    end
end

% Abandonment Counter
C=zeros(nPop,1);

% Array to Hold Best Cost Values
BestCost=zeros(MaxIt,1);

%% ABC Main Loop

for it=1:MaxIt
    
    % Recruited Bees
    for i=1:nPop
        
        % Choose k randomly, not equal to i
        K=[1:i-1 i+1:nPop];
        k=K(randi([1 numel(K)]));
        
        % Define Acceleration Coeff.
        phi=a*unifrnd(-1,+1,VarSize);
        
        % New Bee Position
        newbee.Position=pop(i).Position+phi.*(pop(i).Position-pop(k).Position);
        
        % Evaluation
        newbee.Cost=CostFunction(newbee.Position);
        
        % Comparision
        if newbee.Cost<=pop(i).Cost
            pop(i)=newbee;
        else
            C(i)=C(i)+1;
        end
        
    end
    
    % Calculate Fitness Values and Selection Probabilities
    F=zeros(nPop,1);
    MeanCost = mean([pop.Cost]);
    for i=1:nPop
        F(i) = exp(-pop(i).Cost/MeanCost); % Convert Cost to Fitness
    end
    P=F/sum(F);
    
    % Onlooker Bees
    for m=1:nOnlooker
        
        % Select Source Site
        i=RouletteWheelSelection(P);
        
        % Choose k randomly, not equal to i
        K=[1:i-1 i+1:nPop];
        k=K(randi([1 numel(K)]));
        
        % Define Acceleration Coeff.
        phi=a*unifrnd(-1,+1,VarSize);
        
        % New Bee Position
        newbee.Position=pop(i).Position+phi.*(pop(i).Position-pop(k).Position);
        
        % Evaluation
        newbee.Cost=CostFunction(newbee.Position);
        
        % Comparision
        if newbee.Cost<=pop(i).Cost
            pop(i)=newbee;
        else
            C(i)=C(i)+1;
        end
        
    end
    
    % Scout Bees
    for i=1:nPop
        if C(i)>=L
            pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
            pop(i).Cost=CostFunction(pop(i).Position);
            C(i)=0;
        end
    end
    
    % Update Best Solution Ever Found
    for i=1:nPop
        if pop(i).Cost<=BestSol.Cost
            BestSol=pop(i);
        end
    end
    
    % Store Best Cost Ever Found
    BestCost(it)=BestSol.Cost;
    
    % Display Iteration Information
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
end
    
%% Results

figure;
%plot(BestCost,'LineWidth',2);
semilogy(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;

四、模擬結果

​

圖2人工蜂群演算法收斂曲線

測試統計如下表所示

五、參考文獻及程式碼私信博主

《基於BP神經網路的寧夏水資源需求量預測》

​