一、粒子群算法及RBF簡介
1 粒子群算法簡介
1.1 引言
自然界中的鳥群和魚群的群體行為一直是科學家的研究興趣所在。生物學家Craig Reynolds在1987年提出了一個非常有影響的鳥群聚集模型，在他的仿真中，每一個個體都遵循：避免與鄰域個體相撞：比對鄰域個體的速度；飛向鳥群中心，且整個群體飛向目标。仿真中僅利用上面三條簡單的規則，就可以非常接近地模拟出鳥群飛行的現象。1990年， 生物學家Frank Heppner也提出了鳥類模型， 它的不同之處在于：鳥類被吸引飛到栖息地。在仿真中，一開始每一隻鳥都沒有特定的飛行目标，隻是使用簡單的規則确定自己的飛行方向和飛行速度，當有一隻鳥飛到栖息地時，它周圍的鳥也會跟着飛向栖息地，最終整個鳥群都會落在栖息地。
1995年， 美國社會心理學家James Kennedy和電氣工程師RussellEberhart共同提出了粒子群算法(ParticleS warm Optimization， PSO) ， 該算法的提出是受對鳥類群體行為進行模組化與仿真的研究結果的啟發。他們的模型和仿真算法主要對Frank Heppner的模型進行了修正，以使粒子飛向解空間并在最優解處降落。粒子群算法一經提出，由于其算法簡單，容易實作，立刻引起了進化計算領域學者們的廣泛關注， 形成一個研究熱點。2001年出版的J.Kennedy與R.Eberhart合著的《群體智能》将群體智能的影響進一步擴大[] ， 随後關于粒子群優化算法的研究報告和研究成果大量湧現，繼而掀起了國内外研究熱潮[2-7]。
粒子群優化算法來源于鳥類群體活動的規律性，進而利用群體智能建立一個簡化的模型。它模拟鳥類的覓食行為，将求解問題的搜尋空間比作鳥類的飛行空間，将每隻鳥抽象成一個沒有品質和體積的粒
子，用它來表征問題的一個可能解，将尋找問題最優解的過程看成鳥類尋找食物的過程，進而求解複雜的優化問題。粒子群優化算法與其他進化算法一樣，也是基于“種群”和“進化”的概念，通過個體間
的協作與競争，實作複雜空間最優解的搜尋。同時，它又不像其他進化算法那樣對個體進行交叉、變異、選擇等進化算子操作，而是将群體中的個體看作在l維搜尋空間中沒有品質和體積的粒子，每個粒子以一定的速度在解空間運動， 并向自身曆史最佳位置P best和鄰域曆史最佳位置g best聚集， 實作對候選解的進化。粒子群算法具有很好的生物社會背景而易于了解，由于參數少而容易實作，對非線性、多峰問題均具有較強的全局搜尋能力，在科學研究與工程實踐中得到了廣泛關注。目前，該算法已廣泛應用于函數優化、神經網絡訓練、模式分類、模糊控制等領域。
1.2 粒子群算法理論
1.2.1粒子群算法描述
鳥類在捕食過程中，鳥群成員可以通過個體之間的資訊交流與共享獲得其他成員的發現與飛行經曆。在食物源零星分布并且不可預測的條件下，這種協作機制所帶來的優勢是決定性的，遠遠大于對食物
的競争所引起的劣勢。粒子群算法受鳥類捕食行為的啟發并對這種行為進行模仿，将優化問題的搜尋空間類比于鳥類的飛行空間，将每隻鳥抽象為一個粒子，粒子無品質、無體積，用以表征問題的一個可行解，優化問題所要搜尋到的最優解則等同于鳥類尋找的食物源。粒子群算法為每個粒子制定了與鳥類運動類似的簡單行為規則，使整個粒子群的運動表現出與鳥類捕食相似的特性，進而可以求解複雜的優化問題。
粒子群算法的資訊共享機制可以解釋為一種共生合作的行為，即每個粒子都在不停地進行搜尋，并且其搜尋行為在不同程度上受到群體中其他個體的影響[8]，同時這些粒子還具備對所經曆最佳位置的記
憶能力，即其搜尋行為在受其他個體影響的同時還受到自身經驗的引導。基于獨特的搜尋機制，粒子群算法首先生成初始種群，即在可行解空間和速度空間随機初始化粒子的速度與位置，其中粒子的位置用于表征問題的可行解，然後通過種群間粒子個體的合作與競争來求解優化問題。
1.2.2粒子群算法模組化
粒子群優化算法源自對鳥群捕食行為的研究：一群鳥在區域中随機搜尋食物，所有鳥知道自己目前位置離食物多遠，那麼搜尋的最簡單有效的政策就是搜尋目前離食物最近的鳥的周圍區域。粒子群算法
利用這種模型得到啟示并應用于解決優化問題。在粒子群算法中，每個優化問題的潛在解都是搜尋空間中的一隻鳥，稱之為粒子。所有的粒子都有一個由被優化的函數決定的适應度值，每個粒子還有一個速度決定它們飛翔的方向和距離。然後，粒子們就追随目前的最優粒子在解空間中搜尋[9]。
粒子群算法首先在給定的解空間中随機初始化粒子群，待優化問題的變量數決定了解空間的維數。每個粒子有了初始位置與初始速度，然後通過疊代尋優。在每一次疊代中，每個粒子通過跟蹤兩個“極值”來更新自己在解空間中的空間位置與飛行速度：一個極值就是單個粒子本身在疊代過程中找到的最優解粒子，這個粒子叫作個體極值：另一個極值是種群所有粒子在疊代過程中所找到的最優解粒子，這個粒子是全局極值。上述的方法叫作全局粒子群算法。如果不用種群所有粒子而隻用其中一部分作為該粒子的鄰居粒子，那麼在所有鄰居粒子中的極值就是局部極值，該方法稱為局部粒子群算法。
1.2.3粒子群算法的特點
粒子群算法本質是一種随機搜尋算法，它是一種新興的智能優化技術。該算法能以較大機率收斂于全局最優解。實踐證明，它适合在動态、多目标優化環境中尋優，與傳統優化算法相比，具有較快的計
算速度和更好的全局搜尋能力。
(1)粒子群算法是基于群智能理論的優化算法，通過群體中粒子間的合作與競争産生的群體智能指導優化搜尋。與其他算法相比，粒子群算法是一種高效的并行搜尋算法。
(2)粒子群算法與遺傳算法都是随機初始化種群，使用适應值來評價個體的優劣程度和進行一定的随機搜尋。但粒子群算法根據自己的速度來決定搜尋，沒有遺傳算法的交叉與變異。與進化算法相比，粒子群算法保留了基于種群的全局搜尋政策，但是其采用的速度-位移模型操作簡單，避免了複雜的遺傳操作。
(3)由于每個粒子在算法結束時仍保持其個體極值，即粒子群算法除了可以找到問題的最優解外，還會得到若幹較好的次優解，是以将粒子群算法用于排程和決策問題可以給出多種有意義的方案。
(4)粒子群算法特有的記憶使其可以動态地跟蹤目前搜尋情況并調整其搜尋政策。另外，粒子群算法對種群的大小不敏感，即使種群數目下降時，性能下降也不是很大。
1.3 粒子群算法種類
1.3.1基本粒子群算法
1.3.2标準粒子群算法
引入研究粒子群算法經常用到的兩個概念：一是“探索”，指粒子在一定程度上離開原先的搜尋軌迹，向新的方向進行搜尋，展現了一種向未知區域開拓的能力，類似于全局搜尋；二是“開發”，指粒子在一定程度上繼續在原先的搜尋軌迹上進行更細一步的搜尋，主要指對探索過程中所搜尋到的區域進行更進一步的搜尋。探索是偏離原來的尋優軌迹去尋找一個更好的解，探索能力是一個算法的全局搜尋能力。開發是利用一個好的解，繼續原來的尋優軌迹去搜尋更好的解，它是算法的局部搜尋能力。如何确定局部搜尋能力和全局搜尋能力的比例， 對一個問題的求解過程很重要。1998年， Shi Yuhui等人提出了帶有慣性權重的改進粒子群算法[10]，由于該算法能夠保證較好的收斂效果，是以被預設為标準粒子群算法。其進化過程為：
在式(6.7)中，第一部分表示粒子先前的速度，用于保證算法的全局收斂性能；第二部分、第三部分則使算法具有局部收斂能力。可以看出，式(6.7)中慣性權重w表示在多大程度上保留原來的速度：W
較大，則全局收斂能力較強，局部收斂能力較弱；w較小，則局部收斂能力較強，全局收斂能力較弱。
當w=1時，式(6.7)與式(6.5)完全一樣，表明帶慣性權重的粒子群算法是基本粒子群算法的擴充。實驗結果表明：w在0.8~1.2之間時，粒子群算法有更快的收斂速度；而當w>1.2時，算法則容易陷入局部極值。
另外，在搜尋過程中可以對w進行動态調整：在算法開始時，可給w賦予較大正值，随着搜尋的進行，可以線性地使w逐漸減小，這樣可以保證在算法開始時，各粒子能夠以較大的速度步長在全局範圍内探
測到較好的區域；而在搜尋後期，較小的w值則保證粒子能夠在極值點周圍做精細的搜尋，進而使算法有較大的機率向全局最優解位置收斂。對w進行調整，可以權衡全局搜尋和局部搜尋能力。目前，采用較多的動态慣性權重值是Shi提出的線性遞減權值政策， 其表達式如下：
1.3.3壓縮因子粒子群算法
Clerc等人提出利用限制因子來控制系統行為的最終收斂[11] ， 該方法可以有效搜尋不同的區域，并且能得到高品質的解。壓縮因子法的速度更新公式為：
實驗結果表明：與使用慣性權重的粒子群優化算法相比，使用具
有限制因子的粒子群算法具有更快的收斂速度。
1.3.4離散粒子群算法
基本的粒子群算法是在連續域中搜尋函數極值的有力工具。繼基本粒子群算法之後， Kennedy和Eberhart又提出了一種離散二進制版的粒子群算法[12]。在此離散粒子群方法中，将離散問題空間映射到連續粒子運動空間，并适當修改粒子群算法來求解，在計算上仍保留經典粒子群算法速度-位置更新運算規則。粒子在狀态空間的取值和變化隻限于0和1兩個值， 而速度的每一維vi y代表位置每一位xi取值為1的可能性。是以， 在連續粒子群中的vij更新公式依然保持不變， 但是P best和：best隻在[0， 1] 内取值。其位置更新等式表示如下：
1.4 粒子群算法流程
粒子群算法基于“種群”和“進化”的概念，通過個體間的協作與競争，實作複雜空間最優解的搜尋[13]，其流程如下：
(1)初始化粒子群，包括群體規模N，每個粒子的位置x；和速度Vio
(2) 計算每個粒子的适應度值fit[i] 。
(3) 對每個粒子， 用它的适應度值fit[門和個體極值P best(i)比較。如果fit[i] <P best(i) ， 則用fit[i] 替換掉P best(i) 。
(4) 對每個粒子， 用它的适應度值fit[i] 和全局極值g best比較。如果fit[i] < 8 best， 則用fit[i] 替換g best。
(5)疊代更新粒子的速度v；和位置xj。
(6)進行邊界條件處理。
(7)判斷算法終止條件是否滿足：若是，則結束算法并輸出優化結果；否則傳回步驟(2)。
粒子群算法的運算流程如圖6.1所示。
1.5 關鍵參數說明
在粒子群優化算法中，控制參數的選擇能夠影響算法的性能和效率；如何選擇合适的控制參數使算法性能最佳，是一個複雜的優化問題。在實際的優化問題中，通常根據使用者的經驗來選取控制參數。
粒子群算法的控制參數主要包括：粒子種群規模N，慣性權重w，加速系數c和c， 最大速度Via x， 停止準則， 鄰域結構的設定， 邊界條件處理政策等[14]，
粒子種群規模N
粒子種群大小的選擇視具體問題而定，但是一般設定粒子數為20~50。對于大部分的問題10個粒子，已經可以取得很好的結果：不過對于比較難的問題或者特定類型的問題，粒子的數量可以取到100或
200。另外，粒子數目越大，算法搜尋的空間範圍就越大，也就更容易發現全局最優解；當然，算法運作的時間也越長。
慣性權重w
慣性權重w是标準粒子群算法中非常重要的控制參數，可以用來控制算法的開發和探索能力。慣性權重的大小表示了對粒子目前速度繼承的多少。當慣性權重值較大時，全局尋優能力較強，局部尋優能力
較弱：當慣性權重值較小時，全局尋優能力較弱，局部尋優能力較強。慣性權重的選擇通常有固定權重和時變權重。固定權重就是選擇常數作為慣性權重值，在進化過程中其值保持不變，一般取值為
[0.8，1.2]：時變權重則是設定某一變化區間，在進化過程中按照某種方式逐漸減小慣性權重。時變權重的選擇包括變化範圍和遞減率。固定的慣性權重可以使粒子保持相同的探索和開發能力，而時變權重可以使粒子在進化的不同階段擁有不同的探索和開發能力。
加速常數c1和c2
加速常數c和c 2分别調節向P best和g best方向飛行的最大步長， 它們分别決定粒子個體經驗和群體經驗對粒子運作軌迹的影響，反映粒子群之間的資訊交流。如果cr=c2=0，則粒子将以目前的飛行速度飛到邊界。此時，粒子僅能搜尋有限的區域，是以難以找到最優解。如果q=0，則為“社會”模型，粒子缺乏認知能力，而隻有群體經驗，它的收斂速度較快，但容易陷入局部最優；如果oy=0，則為“認知”模
型，沒有社會的共享資訊，個體之間沒有資訊的互動，是以找到最優解的機率較小，一個規模為D的群體等價于運作了N個各行其是的粒子。是以一般設定c1=C2，通常可以取c1=cg=1.5。這樣，個體經驗和群體經驗就有了同樣重要的影響力，使得最後的最優解更精确。
粒子的最大速度vmax
粒子的速度在空間中的每一維上都有一個最大速度限制值vd max，用來對粒子的速度進行鉗制， 使速度控制在範圍[-Vimax， +va max] 内，這決定問題空間搜尋的力度， 該值一般由使用者自己設定。Vmax是一個非常重要的參數，如果該值太大，則粒子們也許會飛過優秀區域：而如果該值太小，則粒子們可能無法對局部最優區域以外的區域進行充分的探測。它們可能會陷入局部最優，而無法移動足夠遠的距離而跳出局部最優， 達到空間中更佳的位置。研究者指出， 設定Vmax和調整慣性權重的作用是等效的， 是以!max一般用于對種群的初始化進行設定， 即将vmax設定為每維變量的變化範圍， 而不再對最大速度進行細緻的選擇和調節。
停止準則
最大疊代次數、計算精度或最優解的最大停滞步數▲t(或可以接受的滿意解)通常認為是停止準則，即算法的終止條件。根據具體的優化問題，停止準則的設定需同時兼顧算法的求解時間、優化品質和
搜尋效率等多方面性能。
鄰域結構的設定
全局版本的粒子群算法将整個群體作為粒子的鄰域，具有收斂速度快的優點，但有時算法會陷入局部最優。局部版本的粒子群算法将位置相近的個體作為粒子的鄰域，收斂速度較慢，不易陷入局部最優
值。實際應用中，可先采用全局粒子群算法尋找最優解的方向，即得到大緻的結果，然後采用局部粒子群算法在最優點附近進行精細搜尋。
邊界條件處理
當某一維或若幹維的位置或速度超過設定值時，采用邊界條件處理政策可将粒子的位置限制在可行搜尋空間内，這樣能避免種群的膨脹與發散，也能避免粒子大範圍地盲目搜尋，進而提高了搜尋效率。
具體的方法有很多種， 比如通過設定最大位置限制Xmax和最大速度限制Vmax， 當超過最大位置或最大速度時， 在範圍内随機産生一個數值代替，或者将其設定為最大值，即邊界吸收。
二、BP神經網絡簡介
1 BP神經網絡概述
BP（Back Propagation）神經網絡是1986年由Rumelhart和McCelland為首的科研小組提出，參見他們發表在Nature上的論文 Learning representations by back-propagating errors 。
BP神經網絡是一種按誤差逆傳播算法訓練的多層前饋網絡，是目前應用最廣泛的神經網絡模型之一。BP網絡能學習和存貯大量的 輸入-輸出模式映射關系，而無需事前揭示描述這種映射關系的數學方程。它的學習規則是使用最速下降法，通過反向傳播來不斷 調整網絡的權值和門檻值，使網絡的誤差平方和最小。
2 BP算法的基本思想
上一次我們說到，多層感覺器在如何擷取隐層的權值的問題上遇到了瓶頸。既然我們無法直接得到隐層的權值，能否先通過輸出層得到輸出結果和期望輸出的誤差來間接調整隐層的權值呢？BP算法就是采用這樣的思想設計出來的算法，它的基本思想是,學習過程由信号的正向傳播與誤差的反向傳播兩個過程組成。
正向傳播時，輸入樣本從輸入層傳入,經各隐層逐層處理後,傳向輸出層。若輸出層的實際輸出與期望的輸出(教師信号)不符,則轉入誤差的反向傳播階段。
反向傳播時，将輸出以某種形式通過隐層向輸入層逐層反傳,并将誤差分攤給各層的所有單元,進而獲得各層單元的誤差信号,此誤差信号即作為修正各單元權值的依據。這兩個過程的具體流程會在後文介紹。
BP算法的信号流向圖如下圖所示
3 BP網絡特性分析——BP三要素
我們分析一個ANN時，通常都是從它的三要素入手，即
1)網絡拓撲結構；
2)傳遞函數；
3)學習算法。
每一個要素的特性加起來就決定了這個ANN的功能特性。是以，我們也從這三要素入手對BP網絡的研究。
3.1 BP網絡的拓撲結構
上一次已經說了，BP網絡實際上就是多層感覺器，是以它的拓撲結構和多層感覺器的拓撲結構相同。由于單隐層（三層）感覺器已經能夠解決簡單的非線性問題，是以應用最為普遍。三層感覺器的拓撲結構如下圖所示。
一個最簡單的三層BP：
3.2 BP網絡的傳遞函數
BP網絡采用的傳遞函數是非線性變換函數——Sigmoid函數（又稱S函數）。其特點是函數本身及其導數都是連續的，因而在處理上十分友善。為什麼要選擇這個函數，等下在介紹BP網絡的學習算法的時候會進行進一步的介紹。
單極性S型函數曲線如下圖所示。
雙極性S型函數曲線如下圖所示。
3.3 BP網絡的學習算法
BP網絡的學習算法就是BP算法，又叫 δ 算法（在ANN的學習過程中我們會發現不少具有多個名稱的術語）， 以三層感覺器為例，當網絡輸出與期望輸出不等時，存在輸出誤差 E ，定義如下
下面我們會介紹BP網絡的學習訓練的具體過程。
4 BP網絡的訓練分解
訓練一個BP神經網絡，實際上就是調整網絡的權重和偏置這兩個參數，BP神經網絡的訓練過程分兩部分：
前向傳輸，逐層波浪式的傳遞輸出值；
逆向回報，反向逐層調整權重和偏置；
我們先來看前向傳輸。
前向傳輸（Feed-Forward前向回報）
在訓練網絡之前，我們需要随機初始化權重和偏置，對每一個權重取[ − 1 , 1 ] [-1,1][−1,1]的一個随機實數，每一個偏置取[ 0 , 1 ] [0,1][0,1]的一個随機實數，之後就開始進行前向傳輸。
神經網絡的訓練是由多趟疊代完成的，每一趟疊代都使用訓練集的所有記錄，而每一次訓練網絡隻使用一條記錄，抽象的描述如下：
while 終止條件未滿足：
for record:dataset:
trainModel(record)
- 1
- 2
- 3
- 4
4.1 逆向回報（Backpropagation）
4.2 訓練終止條件
每一輪訓練都使用資料集的所有記錄，但什麼時候停止，停止條件有下面兩種：
設定最大疊代次數，比如使用資料集疊代100次後停止訓練
計算訓練集在網絡上的預測準确率，達到一定門限值後停止訓練
5 BP網絡運作的具體流程
5.1 網絡結構
輸入層有n nn個神經元，隐含層有p pp個神經元，輸出層有q qq個神經元。
5.2 變量定義
第九步：判斷模型合理性
判斷網絡誤差是否滿足要求。
當誤差達到預設精度或者學習次數大于設計的最大次數，則結束算法。
否則，選取下一個學習樣本以及對應的輸出期望，傳回第三部，進入下一輪學習。
6 BP網絡的設計
在進行BP網絡的設計是，一般應從網絡的層數、每層中的神經元個數和激活函數、初始值以及學習速率等幾個方面來進行考慮，下面是一些選取的原則。
6.1 網絡的層數
理論已經證明，具有偏差和至少一個S型隐層加上一個線性輸出層的網絡，能夠逼近任何有理函數，增加層數可以進一步降低誤差，提高精度，但同時也是網絡 複雜化。另外不能用僅具有非線性激活函數的單層網絡來解決問題，因為能用單層網絡解決的問題，用自适應線性網絡也一定能解決，而且自适應線性網絡的 運算速度更快，而對于隻能用非線性函數解決的問題，單層精度又不夠高，也隻有增加層數才能達到期望的結果。
6.2 隐層神經元的個數
網絡訓練精度的提高，可以通過采用一個隐含層，而增加其神經元個數的方法來獲得，這在結構實作上要比增加網絡層數簡單得多。一般而言，我們用精度和 訓練網絡的時間來恒量一個神經網絡設計的好壞：
（1）神經元數太少時，網絡不能很好的學習，訓練疊代的次數也比較多，訓練精度也不高。
（2）神經元數太多時，網絡的功能越強大，精确度也更高，訓練疊代的次數也大，可能會出現過拟合(over fitting)現象。
由此，我們得到神經網絡隐層神經元個數的選取原則是：在能夠解決問題的前提下，再加上一兩個神經元，以加快誤差下降速度即可。
6.3 初始權值的選取
一般初始權值是取值在(−1,1)之間的随機數。另外威得羅等人在分析了兩層網絡是如何對一個函數進行訓練後，提出選擇初始權值量級為s√r的政策， 其中r為輸入個數，s為第一層神經元個數。
6.4 學習速率
學習速率一般選取為0.01−0.8，大的學習速率可能導緻系統的不穩定，但小的學習速率導緻收斂太慢，需要較長的訓練時間。對于較複雜的網絡， 在誤差曲面的不同位置可能需要不同的學習速率，為了減少尋找學習速率的訓練次數及時間，比較合适的方法是采用變化的自适應學習速率，使網絡在 不同的階段設定不同大小的學習速率。
6.5 期望誤差的選取
在設計網絡的過程中，期望誤內插補點也應當通過對比訓練後确定一個合适的值，這個合适的值是相對于所需要的隐層節點數來确定的。一般情況下，可以同時對兩個不同 的期望誤內插補點的網絡進行訓練，最後通過綜合因素來确定其中一個網絡。
7 BP網絡的局限性
BP網絡具有以下的幾個問題：
(1)需要較長的訓練時間：這主要是由于學習速率太小所造成的，可采用變化的或自适應的學習速率來加以改進。
(2)完全不能訓練：這主要表現在網絡的麻痹上，通常為了避免這種情況的産生，一是選取較小的初始權值，而是采用較小的學習速率。
(3)局部最小值：這裡采用的梯度下降法可能收斂到局部最小值，采用多層網絡或較多的神經元，有可能得到更好的結果。
8 BP網絡的改進
P算法改進的主要目标是加快訓練速度，避免陷入局部極小值等，常見的改進方法有帶動量因子算法、自适應學習速率、變化的學習速率以及作用函數後縮法等。 動量因子法的基本思想是在反向傳播的基礎上，在每一個權值的變化上加上一項正比于前次權值變化的值，并根據反向傳播法來産生新的權值變化。而自适應學習 速率的方法則是針對一些特定的問題的。改變學習速率的方法的原則是，若連續幾次疊代中，若目标函數對某個權倒數的符号相同，則這個權的學習速率增加， 反之若符号相反則減小它的學習速率。而作用函數後縮法則是将作用函數進行平移，即加上一個常數。
三、部分源代碼
%% 清空環境變量
clc
clear
%% 節點個數
inputnum=8;
hiddennum=8;
outputnum=1;
%% 訓練資料預測資料提取及歸一化
%裝載輸入輸出資料
load temperature
load humidity
tday = 350; %訓練資料時間
sday = 30; %預測資料時間
sdata = (tday+sday)*24; %總的采集資料個數
oddata = 24; %一天采集的資料個數
%在10天的資料中将10分鐘的資料提取出來，總計2880個
temperature = Temp(2:1:sdata,:);
humidity = Humidity(2:1:sdata,:);
% information=[temp humidity];
% 初始化資料
[m,n]=size(temperature);
n = 8;
% m = (m-sday*oddata/10)-oddata/10-n+1;
m = m-sday*oddata-oddata-n+1;%隔天預測
input = zeros(m,n);
output = zeros(m,1);
%找出訓練資料和預測資料
for i = 1:m
input_train(i,1:n/2) = temperature(i:i+3,:);
input_train(i,n/2+1:n) = humidity(i:i+3,:);
output_train(i,:) = temperature(i+oddata,:);
end
[m,n] = size(temperature);
n = 8;
m = (m-tday*oddata)-oddata-n+1;
for i = 1:m
input_test(i,1:n/2) = temperature(i+tday*oddata:i+3+tday*oddata,:);
input_test(i,n/2+1:n) = humidity(i+tday*oddata:i+3+tday*oddata,:);
output_test(i,:) = temperature(i+oddata+tday*oddata,:);
end
input_train=input_train';
output_train=output_train';
input_test = input_test';
output_test=output_test';
%選連樣本輸入輸出資料歸一化
[inputn,inputps]=mapminmax(input_train,0,1);
[outputn,outputps]=mapminmax(output_train,0,1);
%%建構網絡，開始進行PSO的優化
net=newff(inputn,outputn,hiddennum);
% 參數初始化
%粒子群算法中的兩個參數
c1 = 1.49445;
c2 = 1.49445;
maxgen=100; % 進化次數
sizepop=30; %種群規模
% %位置和速度區間
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
%%初始化種群位置、速度、适應度
for i=1:sizepop
pop(i,:)=5*rands(1,81);
V(i,:)=rands(1,81);%rands為生成1*81均與分布在（-1~+1）之間的僞随機數
fitness(i)=fun(pop(i,:),inputnum,hiddennum,outputnum,net,inputn,outputn);
end
%% 個體極值和群體極值
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %個體最佳
fitnessgbest=fitness; %個體最佳适應度值
fitnesszbest=bestfitness; %全局最佳适應度值
%% 疊代尋優
for i=1:maxgen
i;
for j=1:sizepop
%速度更新
V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
V(j,find(V(j,:)>Vmax))=Vmax;
V(j,find(V(j,:)<Vmin))=Vmin;
%種群更新
pop(j,:)=pop(j,:)+0.2*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
% %自适應變異
% pos=unidrnd(81);%unidrnd（N）為産生一個從1到N的離散随機數
% if rand>0.95
% pop(j,pos)=5*rands(1,1);
% end
%适應度值
fitness(j)=fun(pop(j,:),inputnum,hiddennum,outputnum,net,inputn,outputn);
end
for j=1:sizepop
%個體最優更新
if fitness(j) < fitnessgbest(j)
gbest(j,:) = pop(j,:);
fitnessgbest(j) = fitness(j);
end
%群體最優更新
if fitness(j) < fitnesszbest
zbest = pop(j,:);
fitnesszbest = fitness(j);
end
end
yy(i)=fitnesszbest;
end
%該函數功能為傳回預測與實際值之間的絕對值誤差和
function error = fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)
%該函數用來計算适應度值
%x input 個體
%inputnum input 輸入層節點數
%outputnum input 隐含層節點數
%net input 網絡
%inputn input 訓練輸入資料
%outputn input 訓練輸出資料
%error output 個體适應度值
%提取
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
%網絡進化參數
net.trainParam.epochs=50;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
subplot(2,3,5)
plot(tr.epoch,tr.perf','r',tr.epoch,tr.vperf',':g',tr.epoch,tr.tperf','-.b');
% legend(['終止代數＝' num2str(a)])
legend('訓練','确證','測試',-1)
title( '誤差曲線')
ylabel('MSE誤差平方','fontsize',12)
xlabel('進化代數','fontsize',12)
%% 決定系數R^2
fenmu = 0 ;fenzi = 0;
for i= 1:m
fenmu = fenmu + (output_test(:,i)-mean(output_test,2)).^2;
fenzi = fenzi + (output_test(:,i)-BPoutput(:,i)).^2;
end
R = 1 - fenzi/fenmu;
%% 均方根誤差RMSE
MSE = 0;
for i = 1:m
MSE = MSE + ((BPoutput(:,i)-output_test(:,i)).^2);
end
RMSE = sqrt(MSE/m);
%% 畫出各個評價函數值
subplot(2,3,6)
plot(BPoutput,output_test,'.');
% legend(['RMSE = ',num2str(RMSE)]);
title(['R^2 = ',num2str(R),';RMSE = ',num2str(RMSE)])
ylabel('序列真實值','fontsize',12)
xlabel('預測值','fontsize',12)
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