一、粒子群算法簡介
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聚集， 實作對候選解的進化。粒子群算法具有很好的生物社會背景而易于了解，由于參數少而容易實作，對非線性、多峰問題均具有較強的全局搜尋能力，在科學研究與工程實踐中得到了廣泛關注。目前，該算法已廣泛應用于函數優化、神經網絡訓練、模式分類、模糊控制等領域。
2 粒子群算法理論
2.1粒子群算法描述
鳥類在捕食過程中，鳥群成員可以通過個體之間的資訊交流與共享獲得其他成員的發現與飛行經曆。在食物源零星分布并且不可預測的條件下，這種協作機制所帶來的優勢是決定性的，遠遠大于對食物
的競争所引起的劣勢。粒子群算法受鳥類捕食行為的啟發并對這種行為進行模仿，将優化問題的搜尋空間類比于鳥類的飛行空間，将每隻鳥抽象為一個粒子，粒子無品質、無體積，用以表征問題的一個可行解，優化問題所要搜尋到的最優解則等同于鳥類尋找的食物源。粒子群算法為每個粒子制定了與鳥類運動類似的簡單行為規則，使整個粒子群的運動表現出與鳥類捕食相似的特性，進而可以求解複雜的優化問題。
粒子群算法的資訊共享機制可以解釋為一種共生合作的行為，即每個粒子都在不停地進行搜尋，并且其搜尋行為在不同程度上受到群體中其他個體的影響[8]，同時這些粒子還具備對所經曆最佳位置的記
憶能力，即其搜尋行為在受其他個體影響的同時還受到自身經驗的引導。基于獨特的搜尋機制，粒子群算法首先生成初始種群，即在可行解空間和速度空間随機初始化粒子的速度與位置，其中粒子的位置用于表征問題的可行解，然後通過種群間粒子個體的合作與競争來求解優化問題。
2.2粒子群算法模組化
粒子群優化算法源自對鳥群捕食行為的研究：一群鳥在區域中随機搜尋食物，所有鳥知道自己目前位置離食物多遠，那麼搜尋的最簡單有效的政策就是搜尋目前離食物最近的鳥的周圍區域。粒子群算法
利用這種模型得到啟示并應用于解決優化問題。在粒子群算法中，每個優化問題的潛在解都是搜尋空間中的一隻鳥，稱之為粒子。所有的粒子都有一個由被優化的函數決定的适應度值，每個粒子還有一個速度決定它們飛翔的方向和距離。然後，粒子們就追随目前的最優粒子在解空間中搜尋[9]。
粒子群算法首先在給定的解空間中随機初始化粒子群，待優化問題的變量數決定了解空間的維數。每個粒子有了初始位置與初始速度，然後通過疊代尋優。在每一次疊代中，每個粒子通過跟蹤兩個“極值”來更新自己在解空間中的空間位置與飛行速度：一個極值就是單個粒子本身在疊代過程中找到的最優解粒子，這個粒子叫作個體極值：另一個極值是種群所有粒子在疊代過程中所找到的最優解粒子，這個粒子是全局極值。上述的方法叫作全局粒子群算法。如果不用種群所有粒子而隻用其中一部分作為該粒子的鄰居粒子，那麼在所有鄰居粒子中的極值就是局部極值，該方法稱為局部粒子群算法。
2.3粒子群算法的特點
粒子群算法本質是一種随機搜尋算法，它是一種新興的智能優化技術。該算法能以較大機率收斂于全局最優解。實踐證明，它适合在動态、多目标優化環境中尋優，與傳統優化算法相比，具有較快的計
算速度和更好的全局搜尋能力。
(1)粒子群算法是基于群智能理論的優化算法，通過群體中粒子間的合作與競争産生的群體智能指導優化搜尋。與其他算法相比，粒子群算法是一種高效的并行搜尋算法。
(2)粒子群算法與遺傳算法都是随機初始化種群，使用适應值來評價個體的優劣程度和進行一定的随機搜尋。但粒子群算法根據自己的速度來決定搜尋，沒有遺傳算法的交叉與變異。與進化算法相比，粒子群算法保留了基于種群的全局搜尋政策，但是其采用的速度-位移模型操作簡單，避免了複雜的遺傳操作。
(3)由于每個粒子在算法結束時仍保持其個體極值，即粒子群算法除了可以找到問題的最優解外，還會得到若幹較好的次優解，是以将粒子群算法用于排程和決策問題可以給出多種有意義的方案。
(4)粒子群算法特有的記憶使其可以動态地跟蹤目前搜尋情況并調整其搜尋政策。另外，粒子群算法對種群的大小不敏感，即使種群數目下降時，性能下降也不是很大。
3 粒子群算法種類
3.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提出的線性遞減權值政策， 其表達式如下：
3.3壓縮因子粒子群算法
Clerc等人提出利用限制因子來控制系統行為的最終收斂[11] ， 該方法可以有效搜尋不同的區域，并且能得到高品質的解。壓縮因子法的速度更新公式為：
實驗結果表明：與使用慣性權重的粒子群優化算法相比，使用具
有限制因子的粒子群算法具有更快的收斂速度。
3.4離散粒子群算法
基本的粒子群算法是在連續域中搜尋函數極值的有力工具。繼基本粒子群算法之後， Kennedy和Eberhart又提出了一種離散二進制版的粒子群算法[12]。在此離散粒子群方法中，将離散問題空間映射到連續粒子運動空間，并适當修改粒子群算法來求解，在計算上仍保留經典粒子群算法速度-位置更新運算規則。粒子在狀态空間的取值和變化隻限于0和1兩個值， 而速度的每一維vi y代表位置每一位xi取值為1的可能性。是以， 在連續粒子群中的vij更新公式依然保持不變， 但是P best和：best隻在[0， 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所示。
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， 當超過最大位置或最大速度時， 在範圍内随機産生一個數值代替，或者将其設定為最大值，即邊界吸收。
二、部分源代碼
clear all
clc
%opt計算相關量
M=20;%粒子個數
I=4;%水電站個數
T=24;%時間
K=400;%疊代次數
lamda=100;%罰函數懲罰因子
error=0.01;%誤差精度
wini=0.9;wend=0.4;c1=2;c2=2;%粒子群算法慣性因子和學習因子
x=zeros(I,T,M);%粒子位置
v=zeros(I,T,M);%粒子速度
vmax=zeros(I,1);%粒子最大速度
alpha=60;beta=-1.355;gama=0.0105;yeta=0.4968;deta=0.01925;%目标函數e的系數
e=zeros(M,1);%粒子适應值
pbest=zeros(I,T,M);%個體最好位置
epbest=zeros(M,1);%個體最好适應值
gbest=ones(1,1);%全局最好序号
as=5000;bs=19.2;cs=0.002;%煤耗成本系數
f=0;%煤耗成本
%負載相關量
Pd=[1370 1390 1360 1290 1290 1410 1650 2000 2240 2320 2230 2310 2230 2200 2130 2070 2130 2140 2240 2280 2240 2120 1850 1590];%24小時的負載功率
%火電站相關量
Ps=zeros(1,T);%火電輸出功率
Psmax=2500;%火電輸出最大功率
Psmin=500;%火電輸出最小功率
Psbest=zeros(1,T);%火電輸出功率最優解
%水電站相關量
Vh=zeros(I,T);%水電站庫容量
Vhmax=[150 120 240 160];%水電站庫容量最大值
Vhmin=[80 60 100 70];%水電站庫容量最小值
Vhini=[100 80 170 120];%水電站庫容量初始值
Vhend=[120 70 170 140];%水電站庫容量最終值
Vhx=ones(I,1);%調整發電流量Qh時，水電站庫容量計算值與最終值的內插補點
Vhbest=zeros(I,T);%水電站庫容量最優解
Qh=zeros(I,T);%發電流量
Qhmax=[15 15 30 25];%發電流量最大值
Qhmin=[5 6 10 13];%發電流量最小值
Qhbest=zeros(I,T);%發電流量最優解
Ph=zeros(I,T);%水電輸出功率
Phsum=zeros(1,T);%水電輸出功率之和
Phmax=500;%水電輸出功率最大值
Phmin=0;%水電輸出功率最小值
Phbest=zeros(I,T);%水電輸出功率最優解
Phsumbest=zeros(1,T);%水電輸出功率最優解之和
Ih=[10 9 8 7 6 7 8 9 10 11 12 10 11 12 11 10 9 8 7 6 7 8 9 10;
8 8 9 9 8 7 6 7 8 9 9 8 8 9 9 8 7 6 7 8 9 9 8 8;
8.1 8.2 4 2 3 4 3 2 1 1 1 2 4 3 3 2 2 2 1 1 2 2 1 0;
2.8 2.4 1.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%水電站徑流量
c=[-0.0042 -0.42 0.030 0.90 10.0 -50;
-0.0040 -0.30 0.015 1.14 9.5 -70;
-0.0016 -0.30 0.014 0.55 5.5 -40;
-0.0030 -0.31 0.027 1.44 14.0 -90];%水電站輸出功率系數
tao=[2 3 4 0];%水流時滞系數
function[Qh,Pst,Vh,Ph]=adjust(Qh)
%opt計算相關量
I=4;%水電站個數
T=24;%時間
Tadjust=round(rand()*(T-2))+1;%調整時間
error=[0.001 0.001 0.001 0.001];%精度誤差
s=0;%調整次數
smax=20;%調整次數最大值
%負載相關量
Pd=[750 780 700 650 670 800 950 1010 1090 1080 1100 1150 1110 1030 1010 1060 1050 1120 1070 1050 910 860 850 800];%24小時的負載功率
%火電站相關量
Pst=zeros(3,T);%火電輸出功率
Psmax=[175 300 500];%火電輸出最大功率
Psmin=[20 40 50];%火電輸出最小功率
%水電站相關量
Vh=zeros(I,T);%水電站庫容量
Vhini=[100 80 170 120];%水電站庫容量初始值
Vhend=[120 70 170 140];%水電站庫容量最終值
Vhx=ones(I,1);%調整發電流量Qh時，水電站庫容量計算值與最終值的內插補點
Qhmax=[15 15 30 25];%發電流量最大值
Qhmin=[5 6 10 13];%發電流量最小值
Ph=zeros(I,T);%水電輸出功率
Phmin=0;%水電輸出功率最小值
Ih=[10 9 8 7 6 7 8 9 10 11 12 10 11 12 11 10 9 8 7 6 7 8 9 10;
8 8 9 9 8 7 6 7 8 9 9 8 8 9 9 8 7 6 7 8 9 9 8 8;
8.1 8.2 4 2 3 4 3 2 1 1 1 2 4 3 3 2 2 2 1 1 2 2 1 0;
2.8 2.4 1.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%水電站徑流量
c=[-0.0042 -0.42 0.030 0.90 10.0 -50;
-0.0040 -0.30 0.015 1.14 9.5 -70;
-0.0016 -0.30 0.014 0.55 5.5 -40;
-0.0030 -0.31 0.027 1.44 14.0 -90];%水電站輸出功率系數
tao=[2 3 4 0];%水流時滞系數
%計算Vh
for i=1:I
for t=1:T
if t==1
Vh(i,t)=Vhini(i)+Ih(i,t)-Qh(i,t);
else
Vh(i,t)=Vh(i,t-1)+Ih(i,t)-Qh(i,t);
end
if i==3
if t>tao(1)
Vh(i,t)=Vh(i,t)+Qh(1,t-tao(1));
end
if t>tao(2)
Vh(i,t)=Vh(i,t)+Qh(2,t-tao(2));
end
end
if i==4
if t>tao(3)
Vh(i,t)=Vh(i,t)+Qh(3,t-tao(3));
end
end
end
Vhx(i)=Vhend(i)-Vh(i,T);
end
%調整Qh使得水量平衡
for i=1:I
while abs(Vhx(i))>error(i)
for t=Tadjust:T
Qh(i,t)=Qh(i,t)-Vhx(i)/(T-Tadjust)/50;
if Qh(i,t)>Qhmax(i)
Qh(i,t)=Qhmax(i);
end
if Qh(i,t)<Qhmin(i)
Qh(i,t)=Qhmin(i);
end
end
for t=Tadjust:T
if t==1
Vh(i,t)=Vhini(i)+Ih(i,t)-Qh(i,t);
else
Vh(i,t)=Vh(i,t-1)+Ih(i,t)-Qh(i,t);
end
if i==3
if t>tao(1)
Vh(i,t)=Vh(i,t)+Qh(1,t-tao(1));
end
if t>tao(2)
Vh(i,t)=Vh(i,t)+Qh(2,t-tao(2));
end
end
if i==4
if t>tao(3)
Vh(i,t)=Vh(i,t)+Qh(3,t-tao(3));
end
end
end
Vhx(i)=Vhend(i)-Vh(i,T);
s=s+1;
%重新選取調節時間
if s>smax
s=0;
Tadjust=round(rand()*(T-2))+1;
end
end
end
%計算Ph
Ph(i,t)=c(i,1)*Vh(i,t)*Vh(i,t)+c(i,2)*Qh(i,t)*Qh(i,t)+c(i,3)*Vh(i,t)*Qh(i,t)+c(i,4)*Vh(i,t)+c(i,5)*Qh(i,t)+c(i,6);
if Ph(i,t)<Phmin
Ph(i,t)=Phmin;
end
end
end
%Ph每小時求和
Phsum=zeros(1,T);
for t=1:T
for i=1:I
Phsum(1,t)=Phsum(1,t)+Ph(i,t);
end
end
%計算Ps
Psum=zeros(1,T);
for t=1:T
for i=1:2
Ps(t)=Pd(t)-Phsum(t);
Pst(i,t)=Psmin(i)+rand()*(Psmax(i)-Psmin(i));
Psum(1,t)=Psum(1,t)+Pst(i,t);
end
Pst(3,t)=Ps(t)-Psum(1,t);
if Pst<0
Pst(3,t)=0;
end
if Pst(3,t)<Psmin(3)
Pst(3,t)=0;
Pst(1,t)=Pst(1,t)+ Pst(3,t)./2;
Pst(2,t)=Pst(2,t)+ Pst(3,t)./2;
end
end
三、運作結果
四、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 包子陽,餘繼周,楊杉.智能優化算法及其MATLAB執行個體（第2版）[M].電子工業出版社，2016.
[2]張岩,吳水根.MATLAB優化算法源代碼[M].清華大學出版社，2017.