​

一、 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)。

二、粒子群演算法

粒子群演算法是在1995年由Eberhart博士和Kennedy博士一起提出的，它源於對鳥群捕食行為的研究。它的基本核心是利用群體中的個體對資訊的共享從而使整個群體的運動在問題求解空間中產生從無序到有序的演化過程，從而獲得問題的最優解。設想這麼一個場景：一群鳥進行覓食，而遠處有一片玉米地，所有的鳥都不知道玉米地到底在哪裡，但是它們知道自己當前的位置距離玉米地有多遠。那麼找到玉米地的最佳策略，也是最簡單有效的策略就是搜尋目前距離玉米地最近的鳥群的周圍區域。

在PSO中，每個優化問題的解都是搜尋空間中的一隻鳥，稱之為"粒子"，而問題的最優解就對應於鳥群中尋找的"玉米地"。所有的粒子都具有一個位置向量（粒子在解空間的位置）和速度向量（決定下次飛行的方向和速度），並可以根據目標函式來計算當前的所在位置的適應值（fitness value）,可以將其理解為距離"玉米地"的距離。在每次的迭代中，種群中的例子除了根據自身的經驗（歷史位置）進行學習以外，還可以根據種群中最優粒子的"經驗"來學習，從而確定下一次迭代時需要如何調整和改變飛行的方向和速度。就這樣逐步迭代，最終整個種群的例子就會逐步趨於最優解。

上面的解釋可能還比較抽象，下面通過一個簡單的例子來進行說明

​

在一個湖中有兩個人他們之間可以通訊，並且可以探測到自己所在位置的最低點。初始位置如上圖所示，由於右邊比較深，因此左邊的人會往右邊移動一下小船。

​

現在左邊比較深，因此右邊的人會往左邊移動一下小船

一直重複該過程，最後兩個小船會相遇

​

得到一個區域性的最優解

​

​

將每個個體表示為粒子。每個個體在某一時刻的位置表示為，x(t),方向表示為v(t)

​

p（t）為在t時刻x個體的自己的最優解，g(t)為在t時刻所有個體的最優解，v(t)為個體在t時刻的方向，x(t)為個體在t時刻的位置

​​

下一個位置為上圖所示由x,p,g共同決定了

​

種群中的粒子通過不斷地向自身和種群的歷史資訊進行學習，從而可以找到問題的最優解。

但是，在後續的研究中表表明，上述原始的公式中存在一個問題：公式中V的更新太具有隨機性，從而使整個PSO演算法的全域性優化能力很強，但是區域性搜尋能力較差。而實際上，我們需要在演算法迭代初期PSO有著較強的全域性優化能力，而在演算法的後期，整個種群應該具有更強的區域性搜尋能力。所以根據上述的弊端，shi和Eberhart通過引入慣性權重修改了公式，從而提出了PSO的慣性權重模型：

每一個向量的分量表示如下

​

​

其中w稱為是PSO的慣性權重，它的取值介於【0,1】區間，一般應用中均採用自適應的取值方法，即一開始令w=0.9，使得PSO全域性優化能力較強，隨著迭代的深入，引數w進行遞減，從而使的PSO具有較強的區域性優化能力，當迭代結束時，w=0.1。引數c1和c2稱為學習因子，一般設定為1,4961；而r1和r2為介於[0,1]之間的隨機概率值。

整個粒子群優化演算法的演算法框架如下：

step1種群初始化，可以進行隨機初始化或者根據被優化的問題設計特定的初始化方法，然後計算個體的適應值，從而選擇出個體的區域性最優位置向量和種群的全域性最優位置向量。

step2 迭代設定：設定迭代次數，並令當前迭代次數為1

step3 速度更新：更新每個個體的速度向量

step4 位置更新：更新每個個體的位置向量

step5 區域性位置和全域性位置向量更新：更新每個個體的區域性最優解和種群的全域性最優解

step6 終止條件判斷：判斷迭代次數時都達到最大迭代次數，如果滿足，輸出全域性最優解，否則繼續進行迭代，跳轉至step 3。

　　對於粒子群優化演算法的運用，主要是對速度和位置向量迭代運算元的設計。迭代運算元是否有效將決定整個PSO演算法效能的優劣，所以如何設計PSO的迭代運算元是PSO演算法應用的研究重點和難點。

三、PSO優化BP神經網路的步驟

Step1：初始化BP神經網路的權值和閾值
Step2：計算粒子群優化演算法的決策變數長度，選取均方誤差作為優化的目標函式。
Step3：設定演算法停止準則，使用遺傳優化演算法優化神經網路的權值和閾值引數。
Step4：將優化得到的權值和閾值引數賦給BP神經網路。
Step5：優化後的BP神經網路訓練與測試，與優化前的BP神經網路進行誤差分析和精度對比。

四、演示程式碼

clc;clear;close all;

%% 初始化種群

N = 500; % 初始種群個數

d = 24; % 空間維數

ger = 300; % 最大迭代次數


% 設定位置引數限制(矩陣的形式可以多維)

vlimit = [-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;

-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;

-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;

-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;-0.5, 0.5;]; % 設定速度限制

c_1 = 0.8; % 慣性權重

c_2 = 0.5; % 自我學習因子

c_3 = 0.5; % 群體學習因子

for i = 1:d

x(:,i) = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, 1);%初始種群的位置

end

v = 0.5*rand(N, d); % 初始種群的速度

xm = x; % 每個個體的歷史最佳位置

ym = zeros(1, d); % 種群的歷史最佳位置

fxm = 100000*ones(N, 1); % 每個個體的歷史最佳適應度

fym = 10000; % 種群歷史最佳適應度

%% 粒子群工作

iter = 1;

times = 1;

record = zeros(ger, 1); % 記錄器

while iter <= ger

for i=1:N

fx(i) = calfit(x(i,:)) ; % 個體當前適應度

end

for i = 1:N

if fxm(i) > fx(i)

fxm(i) = fx(i); % 更新個體歷史最佳適應度

xm(i,:) = x(i,:); % 更新個體歷史最佳位置

end

end

if fym > min(fxm)

[fym, nmax] = min(fxm); % 更新群體歷史最佳適應度

ym = xm(nmax, :); % 更新群體歷史最佳位置

end

v = v * c_1 + c_2 * rand *(xm - x) + c_3 * rand *(repmat(ym, N, 1) - x);% 速度更新

% 邊界速度處理

for i=1:d

for j=1:N

if v(j,i)>vlimit(i,2)

v(j,i)=vlimit(i,2);

end

if v(j,i) < vlimit(i,1)

v(j,i)=vlimit(i,1);

end

end

end

x = x + v;% 位置更新

% 邊界位置處理

for i=1:d

for j=1:N

if x(j,i)>limit(i,2)

x(j,i)=limit(i,2);

end

if x(j,i) < limit(i,1)

x(j,i)=limit(i,1);

end

end

end

record(iter) = fym;%最大值記錄

iter = iter+1;

times=times+1;

end

disp(['最小值：',num2str(fym)]);

disp(['變數取值：',num2str(ym)]);

figure

plot(record)

xlabel('迭代次數');

ylabel('適應度值')

四、模擬結果

​

​

​

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

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

​

​