一、粒子群算法的歷史
粒子群算法源於復雜適應系統（Complex Adaptive System,CAS）。CAS理論於1994年正式提出，CAS中的成員稱為主體。比方研究鳥群系統，每一個鳥在這個系統中就稱為主體。主體有適應性，它能夠與環境及其它的主體進行交流，而且依據交流的過程“學習”或“積累經驗”改變自身結構與行為。整個系統的演變或進化包括：新層次的產生（小鳥的出生）；分化和多樣性的出現（鳥群中的鳥分成很多小的群）；新的主題的出現（鳥尋找食物過程中，不斷發現新的食物）。
所以CAS系統中的主體具有4個基本特點（這些特點是粒子群算法發展變化的依據）：
首先，主體是主動的、活動的。
主體與環境及其它主體是相互影響、相互作用的，這樣的影響是系統發展變化的主要動力。
環境的影響是宏觀的，主體之間的影響是微觀的，宏觀與微觀要有機結合。
最后，整個系統可能還要受一些隨機因素的影響。
粒子群算法就是對一個CAS系統－－－鳥群社會系統的研究得出的。
粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy於1995年提出，它的基本概念源於對鳥群覓食行為的研究。設想這樣一個場景:一群鳥在隨機搜尋食物，在這個區域里僅僅有一塊食物，全部的鳥都不知道食物在哪里，可是它們知道當前的位置離食物還有多遠。那么找到食物的最優策略是什么呢?最簡單有效的就是搜尋眼下離食物近期的鳥的周圍區域。
PSO算法就從這樣的生物種群行為特性中得到啟示並用於求解優化問題。在PSO中，每一個優化問題的潛在解都能夠想象成d維搜索空間上的一個點，我們稱之為“粒子”（Particle），全部的粒子都有一個被目標函數決定的適應值(Fitness Value )，每一個粒子另一個速度決定他們飛翔的方向和距離，然后粒子們就追隨當前的最優粒子在解空間中搜索。Reynolds對鳥群飛行的研究發現。鳥僅僅是追蹤它有限數量的鄰居但終於的總體結果是整個鳥群好像在一個中心的控制之下.即復雜的全局行為是由簡單規則的相互作用引起的。
二、粒子群算法的詳細表述
上面羅嗦了半天，那些都是科研工作者寫論文的語氣，只是，PSO的歷史就像上面說的那樣。以下通俗的解釋PSO算法。
PSO算法就是模擬一群鳥尋找食物的過程，每一個鳥就是PSO中的粒子，也就是我們須要求解問題的可能解，這些鳥在尋找食物的過程中，不停改變自己在空中飛行的位置與速度。大家也能夠觀察一下，鳥群在尋找食物的過程中，開始鳥群比較分散，逐漸這些鳥就會聚成一群，這個群忽高忽低、忽左忽右，直到最后找到食物。這個過程我們轉化為一個數學問題。尋找函數 y=1-cos(3*x)*exp(-x)的在[0,4]最大值。該函數的圖形例如以下：
當x=0.9350-0.9450，達到最大值y=1.3706。為了得到該函數的最大值，我們在[0，4]之間隨機的灑一些點，為了演示，我們放置兩個點，而且計算這兩個點的函數值，同一時候給這兩個點設置在[0，4]之間的一個速度。以下這些點就會依照一定的公式更改自己的位置，到達新位置后，再計算這兩個點的值，然后再依照一定的公式更新自己的位置。直到最后在y=1.3706這個點停止自己的更新。這個過程與粒子群算法作為對比例如以下：
這兩個點就是粒子群算法中的粒子。
該函數的最大值就是鳥群中的食物
計算兩個點函數值就是粒子群算法中的適應值，計算用的函數就是粒子群算法中的適應度函數。
更新自己位置的一定公式就是粒子群算法中的位置速度更新公式。
以下演示一下這個算法執行一次的大概過程：
第一次初始化
第一次更新位置
第二次更新位置
第21次更新
最后的結果（30次迭代）
最后全部的點都集中在最大值的地方。
呵呵，如今粒子群算法的大概思想就講到這里。下節介紹標准的粒子群算法。
粒子群算法(2)----標准的粒子群算法
在上一節的敘述中，唯一沒有給大家介紹的就是函數的這些隨機的點（粒子）是怎樣運動的，僅僅是說依照一定的公式更新。這個公式就是粒子群算法中的位置速度更新公式。以下就介紹這個公式是什么。在上一節中我們求取函數y=1-cos(3*x)*exp(-x)的在[0,4]最大值。並在[0,4]之間放置了兩個隨機的點，這些點的坐標如果為x1=1.5； x2=2.5；這里的點是一個標量，可是我們常常遇到的問題可能是更一般的情況－－x為一個矢量的情況，比方二維的情況 z=2*x1+3*x22的情況。這個時候我們的每一個粒子為二維，記粒子P1＝(x11,x12),P2=(x21,x22),P3=(x31,x32)，......Pn=(xn1,xn2)。這里n為粒子群群體的規模，也就是這個群中粒子的個數，每一個粒子的維數為2。更一般的是粒子的維數為q，這樣在這個種群中有n個粒子，每一個粒子為q 維。
由n個粒子組成的群體對Q維（就是每一個粒子的維數）空間進行搜索。每一個粒子表示為：xi＝（xi1,xi2,xi3,...,xiQ），每一個粒子相應的速度能夠表示為vi=(vi1,vi2,vi3,....,viQ)，每一個粒子在搜索時要考慮兩個因素：
1。自己搜索到的歷史最優值 pi ，pi=(pi1,pi2,....,piQ)，i=1,2,3,....,n。
2。全部粒子搜索到的最優值pg，pg=(pg1,pg2,....,pgQ)，注意這里的pg僅僅有一個。
以下給出粒子群算法的位置速度更新公式：
這里有幾個重要的參數須要大家記憶，由於在以后的解說中將會常常常使用到：
它們是：
是保持原來速度的系數，所以叫做慣性權重。
是粒子跟蹤自己歷史最優值的權重系數，它表示粒子自身的認識，所以叫“認知”。通常設置為2。
是粒子跟蹤群體最優值的權重系數，它表示粒子對整個群體知識的認識，所以叫做“社會知識”，常常叫做“社會”。通常設置為2。
是[0,1]區間內均勻分布的隨機數。
是對位置更新的時候，在速度前面加的一個系數，這個系數我們叫做約束因子。通常設置為1。
這樣一個標准的粒子群算法就結束了。
以下對整個主要的粒子群的過程給一個簡單的圖形表示：
推斷終止條件可是設置適應值到達一定的數值或者循環一定的次數。
注意：這里的粒子是同一時候跟蹤自己的歷史最優值與全局（群體）最優值來改變自己的位置預速度的，所以又叫做全局版本號的標准粒子群優化算法。
粒子群算法(3)----標准的粒子群算法(局部版本號)
在全局版的標准粒子群算法中，每一個粒子的速度的更新是依據兩個因素來變化的，這兩個因素是：1. 粒子自己歷史最優值pi。2. 粒子群體的全局最優值pg。如果改變粒子速度更新公式，讓每一個粒子的速度的更新依據以下兩個因素更新，A. 粒子自己歷史最優值pi。B. 粒子鄰域內粒子的最優值pnk。其余保持跟全局版的標准粒子群算法一樣，這個算法就變為局部版的粒子群算法。
一般一個粒子i 的鄰域隨着迭代次數的添加而逐漸添加，開始第一次迭代，它的鄰域為0，隨着迭代次數鄰域線性變大，最后鄰域擴展到整個粒子群，這時就變成全局版本號的粒子群算法了。經過實踐證明：全局版本號的粒子群算法收斂速度快，可是easy陷入局部最優。局部版本號的粒子群算法收斂速度慢，可是非常難陷入局部最優。如今的粒子群算法大都在收斂速度與擺脫局部最優這兩個方面下功夫。事實上這兩個方面是矛盾的。看怎樣更好的折中了。
依據取鄰域的方式的不同，局部版本號的粒子群算法有非常多不同的實現方法。
第一種方法：依照粒子的編號取粒子的鄰域，取法有四種：1，環形取法 2，隨機環形取法 3，輪形取法 4，隨機輪形取法。
1 環形2 隨機環形
3 輪形 4隨機輪形
由於后面有以環形取法實現的算法，對環形取法在這里做一點點說明：以粒子1為例，當鄰域是0的時候，鄰域是它本身，當鄰域是1時，鄰域為2，8；當鄰域是2時，鄰域是2，3，7，8；......，以此類推，一直到鄰域為4，這個時候，鄰域擴展到整個樣例群體。據文獻介紹（國外的文獻），採用輪形拓撲結構，PSO的效果非常好。
另外一種方法：依照粒子的歐式距離取粒子的鄰域
在第一種方法中，依照粒子的編號來得到粒子的鄰域，可是這些粒子事實上可能在實際位置上並不相鄰，於是Suganthan提出基於空間距離的划分方案，在迭代中計算每一個粒子與群中其它粒子的距離。記錄不論什么2個粒子間的的最大距離為dm。對每一粒子依照||xa-xb||/dm計算一個比值。當中||xa-xb||是當前粒子a到b的距離。而選擇閾值frac依據迭代次數而變化。當另一粒子b滿足||xa-xb||/dm<frac時，覺得b成為當前粒子的鄰域。
這樣的辦法經過實驗，取得較好的應用效果，可是由於要計算全部粒子之間的距離，計算量大，且須要非常大的存儲空間，所以，該方法一般不常常使用。
粒子群算法(4)----粒子群算法分類
粒子群算法主要分為4個大的分支：
（1）標准粒子群算法的變形
在這個分支中，主要是對標准粒子群算法的慣性因子、收斂因子（約束因子）、“認知”部分的c1，“社會”部分的c2進行變化與調節，希望獲得好的效果。
慣性因子的原始版本號是保持不變的，后來有人提出隨着算法迭代的進行，慣性因子須要逐漸減小的思想。算法開始階段，大的慣性因子能夠是算法不easy陷入局部最優，到算法的后期，小的慣性因子能夠使收斂速度加快，使收斂更加平穩，不至於出現振盪現象。經過本人測試，動態的減小慣性因子w，的確能夠使算法更加穩定，效果比較好。可是遞減慣性因子採用什么樣的方法呢？人們首先想到的是線型遞減，這樣的策略的確非常好，可是是不是最優的呢？於是有人對遞減的策略作了研究，研究結果指出：線型函數的遞減優於凸函數的遞減策略，可是凹函數的遞減策略又優於線型的遞減，經過本人測試，實驗結果基本符合這個結論，可是效果不是非常明顯。
對於收斂因子，經過證明如果收斂因子取0.729,能夠確保算法的收斂，可是不能保證算法收斂到全局最優，經過本人測試，取收斂因子為0.729效果較好。對於社會與認知的系數c2,c1也有人提出：c1先大后小，而c2先小后大的思想，由於在算法執行初期，每一個鳥要有大的自己的認知部分而又比較小的社會部分，這個與我們自己一群人找東西的情形比較接近，由於在我們找東西的初期，我們基本依靠自己的知識取尋找，而后來，我們積累的經驗越來越豐富，於是大家開始逐漸達成共識（社會知識），這樣我們就開始依靠社會知識來尋找東西了。
2007年希臘的兩位學者提出將收斂速度比較快的全局版本號的粒子群算法與不easy陷入局部最優的局部版本號的粒子群算法相結合的辦法，利用的公式是
v＝n*v(全局版本號)＋（1－n）*v(局部版本號) 速度更新公式，v代表速度
w（k＋1）＝w(k)＋v 位置更新公式
該算法在文獻中討論了系數n取各種不同情況的情況，而且執行來了20000次來分析各種系數的結果。
（2）粒子群算法的混合
這個分支主要是將粒子群算法與各種算法相混合，有人將它與模擬退火算法相混合，有些人將它與單純形方法相混合。可是最多的是將它與遺傳算法的混合。依據遺傳算法的三種不同算子能夠生成3中不同的混合算法。
粒子群算法與選擇算子的結合，這里相混合的思想是：在原來的粒子群算法中，我們選擇粒子群群體的最優值作為pg，可是相結合的版本號是依據全部粒子的適應度的大小給每一個粒子賦予一個被選中的概率，然后依據概率對這些粒子進行選擇，被選中的粒子作為pg，其它的情況都不變。這樣的算法能夠在算法執行過程中保持粒子群的多樣性，可是致命的缺點是收斂速度緩慢。
粒子群算法與雜交算子的結合，結合的思想與遺傳算法的基本一樣，在算法執行過程中依據適應度的大小，粒子之間能夠兩兩雜交，比方用一個非常easy的公式
w（新）＝n×w1＋（1－n）×w2；
w1與w2就是這個新粒子的父輩粒子。這樣的算法能夠在算法的執行過程中引入新的粒子，可是算法一旦陷入局部最優，那么粒子群算法將非常難擺脫局部最優。
粒子群算法與變異算子的結合，結合的思想：測試全部粒子與當前最優的距離，當距離小於一定的數值的時候，能夠拿出全部粒子的一個百分比（如10％）的粒子進行隨機初始化，讓這些粒子又一次尋找最優值。
（3）二進制粒子群算法
最初的PSO是從解決連續優化問題發展起來的.Eberhart等又提出了PSO的離散二進制版.用來解決project實際中的組合優化問題。他們在提出的模型中將粒子的每一維及粒子本身的歷史最優、全局最優限制為1或0，而速度不作這樣的限制。用速度更新位置時，設定一個閾值，當速度高於該閾值時，粒子的位置取1，否則取0。二進制PSO與遺傳算法在形式上非常類似，但實驗結果顯示，在大多數測試函數中，二進制PSO比遺傳算法速度快，尤其在問題的維數添加時
（4）協同粒子群算法
協同PSO，該方法將粒子的D維分到D個粒子群中，每一個粒子群優化一維向量，評價適應度時將這些分量合並為一個完整的向量。比如第i個粒子群，除第i個分量外，其它D-1個分量都設為最優值，不斷用第i個粒子群中的粒子替換第i個分量，直到得到第i維的最優值，其它維同樣。為將有聯系的分量划分在一個群，可將D維向量分配到m個粒子群優化，則前D mod m個粒子群的維數是D/m的向上取整。后m－(D mod m)個粒子群的維數是D/m的向下取整。協同PSO在某些問題上有更快的收斂速度，但該算法easy被欺騙。
主要的粒子群算法的分支就着4個，大部分的粒子群算法都環繞着這4個分支在變化，當中粒子群算法的變形居多，從根本上來說，差點兒沒有什么新的思想的提出。
粒子群算法(5)-----標准粒子群算法的實現
標准粒子群算法的實現思想基本依照粒子群算法(2)----標准的粒子群算法的講述實現。主要分為3個函數。第一個函數為粒子群初始化函數
InitSwarm(SwarmSize......AdaptFunc)其主要作用是初始化粒子群的粒子，並設定粒子的速度、位置在一定的范圍內。本函數所採用的數據結構例如以下所看到的：
表ParSwarm記錄的是粒子的位置、速度與當前的適應度值，我們用W來表示位置，用V來代表速度，用F來代表當前的適應度值。在這里我們如果粒子個數為N，每一個粒子的維數為D。
W1,1 |
W1,2 |
... |
W1,D |
V1,1 |
V1,2 |
... |
V1,D-1 |
V1,D |
F1 |
第1個粒子 |
W2,1 |
W2,2 |
... |
W2,D |
V2,1 |
V2,2 |
... |
V2,D-1 |
V2,D |
F2 |
第2個粒子 |
... |
... |
... |
... |
... |
... |
... |
... |
... |
... |
....... |
WN-1,1 |
WN-1,2 |
... |
WN-1,D-1 |
VN-1,1 |
VN-1,2 |
... |
VN-1,D-1 |
VN-1,D |
FN-1 |
第N-1個粒子 |
WN,1 |
WN,2 |
... |
WN,D |
VN,1 |
VN,2 |
... |
VN,D-1 |
VN,D |
FN |
第N個粒子 |
表OptSwarm記錄每一個粒子的歷史最優解（粒子歷史最好的適應度）以及全部粒子搜索到的全局最優解。用Wg代表全局最優解，W.,1代表每一個粒子的歷史最優解。粒子群初始化階段表OptSwarm的前N行與表ParSwarm中的同樣，而Wg的值為表ParSwarm中適應度值的最大值相應的行。
Wj,1 |
Wj,2 |
... |
Wj,D-1 |
Wj,D |
第1個粒子的歷史最優解 |
Wk,1 |
Wk,2 |
... |
Wk,D-1 |
Wk,D |
第2個粒子的歷史最優解 |
... |
... |
... |
... |
... |
... |
Wl,1 |
Wl,2 |
... |
Wl,D-1 |
Wl,D |
第N-1個粒子的歷史最優解 |
Wm,1 |
Wm,2 |
... |
Wm,D-1 |
Wm,D |
第N個粒子的歷史最優解 |
Wg,1 |
Wg,2 |
... |
Wg,D-1 |
Wg,D |
全局粒子的歷史最優解 |
依據這樣的思想MATLAB代碼例如以下：
function [ParSwarm,OptSwarm]=InitSwarm(SwarmSize,ParticleSize,ParticleScope,AdaptFunc)
%功能描寫敘述：初始化粒子群，限定粒子群的位置以及速度在指定的范圍內
%[ParSwarm,OptSwarm,BadSwarm]=InitSwarm(SwarmSize,ParticleSize,ParticleScope,AdaptFunc)
%
%輸入參數：SwarmSize:種群大小的個數
%輸入參數：ParticleSize：一個粒子的維數
%輸入參數：ParticleScope:一個粒子在運算中各維的范圍；
% ParticleScope格式:
% 3維粒子的ParticleScope格式:
% [x1Min,x1Max
% x2Min,x2Max
% x3Min,x3Max]
%
%輸入參數：AdaptFunc：適應度函數
%
%輸出：ParSwarm初始化的粒子群
%輸出：OptSwarm粒子群當前最優解與全局最優解
%
%使用方法[ParSwarm,OptSwarm,BadSwarm]=InitSwarm(SwarmSize,ParticleSize,ParticleScope,AdaptFunc);
%
%異常：首先保證該文件在Matlab的搜索路徑中，然后查看相關的提示信息。
%
%編制人：XXX
%編制時間：2007.3.26
%參考文獻：無
%
%容錯控制
if nargin~=4
error('輸入的參數個數錯誤。')
end
if nargout<2
error('輸出的參數的個數太少，不能保證以后的執行。');
end
[row,colum]=size(ParticleSize);
if row>1|colum>1
error('輸入的粒子的維數錯誤，是一個1行1列的數據。');
end
[row,colum]=size(ParticleScope);
if row~=ParticleSize|colum~=2
error('輸入的粒子的維數范圍錯誤。');
end
%初始化粒子群矩陣
%初始化粒子群矩陣，全部設為[0-1]隨機數
%rand('state',0);
ParSwarm=rand(SwarmSize,2*ParticleSize+1);
%對粒子群中位置,速度的范圍進行調節
for k=1:ParticleSize
ParSwarm(:,k)=ParSwarm(:,k)*(ParticleScope(k,2)-ParticleScope(k,1))+ParticleScope(k,1);
%調節速度，使速度與位置的范圍一致
ParSwarm(:,ParticleSize+k)=ParSwarm(:,ParticleSize+k)*(ParticleScope(k,2)-ParticleScope(k,1))+ParticleScope(k,1);
end
%對每一個粒子計算其適應度函數的值
for k=1:SwarmSize
ParSwarm(k,2*ParticleSize+1)=AdaptFunc(ParSwarm(k,1:ParticleSize));
end
%初始化粒子群最優解矩陣
OptSwarm=zeros(SwarmSize+1,ParticleSize);
%粒子群最優解矩陣全部設為零
[maxValue,row]=max(ParSwarm(:,2*ParticleSize+1));
%尋找適應度函數值最大的解在矩陣中的位置(行數)
OptSwarm=ParSwarm(1:SwarmSize,1:ParticleSize);
OptSwarm(SwarmSize+1,:)=ParSwarm(row,1:ParticleSize);
以下的函數BaseStepPso實現了標准全局版粒子群算法的單步更新位置速度的功能
function [ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%功能描寫敘述：全局版本號：主要的粒子群算法的單步更新位置,速度的算法
%
%[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%
%輸入參數：ParSwarm:粒子群矩陣，包括粒子的位置，速度與當前的目標函數值
%輸入參數：OptSwarm：包括粒子群個體最優解與全局最優解的矩陣
%輸入參數：ParticleScope:一個粒子在運算中各維的范圍；
%輸入參數：AdaptFunc：適應度函數
%輸入參數：LoopCount：迭代的總次數
%輸入參數：CurCount：當前迭代的次數
%
%返回值：含意同輸入的同名參數
%
%使用方法：[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%
%異常：首先保證該文件在Matlab的搜索路徑中，然后查看相關的提示信息。
%
%編制人：XXX
%編制時間：2007.3.26
%參考文獻：XXX
%參考文獻：XXX
%
%改動記錄
%----------------------------------------------------------------
%2007.3.27
%改動人：XXX
% 加入2*unifrnd(0,1).*SubTract1(row,:)中的unifrnd(0,1)隨機數，使性能大為提高
%參照基於MATLAB的粒子群優化算法程序設計
%
% 總體評價：使用這個版本號的調節系數，效果比較好
%
%容錯控制
if nargin~=8
error('輸入的參數個數錯誤。')
end
if nargout~=2
error('輸出的個數太少，不能保證循環迭代。')
end
%開始單步更新的操作
%*********************************************
%*****更改以下的代碼，能夠更改慣性因子的變化*****
%---------------------------------------------------------------------
%線形遞減策略
w=MaxW-CurCount*((MaxW-MinW)/LoopCount);
%---------------------------------------------------------------------
%w固定不變策略
%w=0.7;
%---------------------------------------------------------------------
%參考文獻：陳貴敏，賈建援，韓琪，粒子群優化算法的慣性權值遞減策略研究，西安交通大學學報，2006，1
%w非線形遞減，以凹函數遞減
%w=(MaxW-MinW)*(CurCount/LoopCount)^2+(MinW-MaxW)*(2*CurCount/LoopCount)+MaxW;
%---------------------------------------------------------------------
%w非線形遞減，以凹函數遞減
%w=MinW*(MaxW/MinW)^(1/(1+10*CurCount/LoopCount));
%*****更改上面的代碼，能夠更改慣性因子的變化*****
%*********************************************
%得到粒子群群體大小以及一個粒子維數的信息
[ParRow,ParCol]=size(ParSwarm);
%得到粒子的維數
ParCol=(ParCol-1)/2;
SubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);
%*********************************************
%*****更改以下的代碼，能夠更改c1,c2的變化*****
c1=2;
c2=2;
%---------------------------------------------------------------------
%con=1;
%c1=4-exp(-con*abs(mean(ParSwarm(:,2*ParCol+1))-AdaptFunc(OptSwarm(ParRow+1,:))));
%c2=4-c1;
%----------------------------------------------------------------------
%*****更改上面的代碼，能夠更改c1,c2的變化*****
%*********************************************
for row=1:ParRow
SubTract2=OptSwarm(ParRow+1,:)-ParSwarm(row,1:ParCol);
TempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+2*unifrnd(0,1).*SubTract1(row,:)+2*unifrnd(0,1).*SubTract2;
%限制速度的代碼
for h=1:ParCol
if TempV(:,h)>ParticleScope(h,2)
TempV(:,h)=ParticleScope(h,2);
end
if TempV(:,h)<-ParticleScope(h,2)
TempV(:,h)=-ParticleScope(h,2)+1e-10;
%加1e-10防止適應度函數被零除
end
end
%更新速度
ParSwarm(row,ParCol+1:2*ParCol)=TempV;
%*********************************************
%*****更改以下的代碼，能夠更改約束因子的變化*****
%---------------------------------------------------------------------
%a=1;
%---------------------------------------------------------------------
a=0.729;
%*****更改上面的代碼，能夠更改約束因子的變化*****
%*********************************************
%限制位置的范圍
TempPos=ParSwarm(row,1:ParCol)+a*TempV;
for h=1:ParCol
if TempPos(:,h)>ParticleScope(h,2)
TempPos(:,h)=ParticleScope(h,2);
end
if TempPos(:,h)<=ParticleScope(h,1)
TempPos(:,h)=ParticleScope(h,1)+1e-10;
end
end
%更新位置
ParSwarm(row,1:ParCol)=TempPos;
%計算每一個粒子的新的適應度值
ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));
if ParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))
OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);
end
end
%for循環結束
%尋找適應度函數值最大的解在矩陣中的位置(行數)，進行全局最優的改變
[maxValue,row]=max(ParSwarm(:,2*ParCol+1));
if AdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow+1,:))
OptSwarm(ParRow+1,:)=ParSwarm(row,1:ParCol);
end
這兩個函數給出以后，須要一個函數來把這兩個函數組裝起來，以此實現一個完整的粒子群算法，這個函數就是PsoProcess
代碼例如以下：
function [Result,OnLine,OffLine,MinMaxMeanAdapt]=PsoProcess(SwarmSize,ParticleSize,ParticleScope,InitFunc,StepFindFunc,AdaptFunc,IsStep,IsDraw,LoopCount,IsPlot)
%功能描寫敘述：一個循環n次的PSO算法完整過程，返回這次執行的最小與最大的平均適應度,以及在線性能與離線性能
%[Result,OnLine,OffLine,MinMaxMeanAdapt]=PsoProcess(SwarmSize,ParticleSize,ParticleScope,InitFunc,StepFindFunc,AdaptFunc,IsStep,IsDraw,LoopCount,IsPlot)
%輸入參數：SwarmSize:種群大小的個數
%輸入參數：ParticleSize：一個粒子的維數
%輸入參數：ParticleScope:一個粒子在運算中各維的范圍；
% ParticleScope格式:
% 3維粒子的ParticleScope格式:
% [x1Min,x1Max
% x2Min,x2Max
% x3Min,x3Max]
%
%輸入參數:InitFunc:初始化粒子群函數
%輸入參數:StepFindFunc:單步更新速度，位置函數
%輸入參數：AdaptFunc：適應度函數
%輸入參數：IsStep：是否每次迭代暫停；IsStep＝0，不暫停，否則暫停。缺省不暫停
%輸入參數：IsDraw：是否圖形化迭代過程；IsDraw＝0，不圖形化迭代過程，否則，圖形化表示。缺省不圖形化表示
%輸入參數：LoopCount：迭代的次數；缺省迭代100次
%輸入參數：IsPlot：控制是否繪制在線性能與離線性能的圖形表示；IsPlot=0,不顯示；
% IsPlot=1；顯示圖形結果。缺省IsPlot=1
%
%返回值：Result為經過迭代后得到的最優解
%返回值：OnLine為在線性能的數據
%返回值：OffLine為離線性能的數據
%返回值：MinMaxMeanAdapt為本次完整迭代得到的最小與最大的平均適應度
%
%使用方法[Result,OnLine,OffLine,MinMaxMeanAdapt]=PsoProcess(SwarmSize,ParticleSize,ParticleScope,InitFunc,StepFindFunc,AdaptFunc,IsStep,IsDraw,LoopCount,IsPlot);
%
%異常：首先保證該文件在Matlab的搜索路徑中，然后查看相關的提示信息。
%
%編制人：XXX
%編制時間：2007.3.26
%參考文獻：XXXXX%
%改動記錄：
%加入MinMaxMeanAdapt，以得到性能評估數據
%改動人：XXX
%改動時間：2007.3.27
%參考文獻：XXX.
%容錯控制
if nargin<4
error('輸入的參數個數錯誤。')
end
[row,colum]=size(ParticleSize);
if row>1|colum>1
error('輸入的粒子的維數錯誤，是一個1行1列的數據。');
end
[row,colum]=size(ParticleScope);
if row~=ParticleSize|colum~=2
error('輸入的粒子的維數范圍錯誤。');
end
%設置缺省值
if nargin<7
IsPlot=1;
LoopCount=100;
IsStep=0;
IsDraw=0;
end
if nargin<8
IsPlot=1;
IsDraw=0;
LoopCount=100;
end
if nargin<9
LoopCount=100;
IsPlot=1;
end
if nargin<10
IsPlot=1;
end
%控制是否顯示2維以下粒子維數的尋找最優的過程
if IsDraw~=0
DrawObjGraphic(ParticleSize,ParticleScope,AdaptFunc);
end
%初始化種群
[ParSwarm,OptSwarm]=InitFunc(SwarmSize,ParticleSize,ParticleScope,AdaptFunc)
%在測試函數圖形上繪制初始化群的位置
if IsDraw~=0
if 1==ParticleSize
for ParSwarmRow=1:SwarmSize
plot([ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,1)],[ParSwarm(ParSwarmRow,3),0],'r*-','markersize',8);
text(ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,3),num2str(ParSwarmRow));
end
end
if 2==ParticleSize
for ParSwarmRow=1:SwarmSize
stem3(ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,2),ParSwarm(ParSwarmRow,5),'r.','markersize',8);
end
end
end
%暫停讓抓圖
if IsStep~=0
disp('開始迭代，按隨意鍵：')
pause
end
%開始更新算法的調用
for k=1:LoopCount
%顯示迭代的次數：
disp('----------------------------------------------------------')
TempStr=sprintf('第 %g 此迭代',k);
disp(TempStr);
disp('----------------------------------------------------------')
%調用一步迭代的算法
[ParSwarm,OptSwarm]=StepFindFunc(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,0.95,0.4,LoopCount,k)
%在目標函數的圖形上繪制2維以下的粒子的新位置
if IsDraw~=0
if 1==ParticleSize
for ParSwarmRow=1:SwarmSize
plot([ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,1)],[ParSwarm(ParSwarmRow,3),0],'r*-','markersize',8);
text(ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,3),num2str(ParSwarmRow));
end
end
if 2==ParticleSize
for ParSwarmRow=1:SwarmSize
stem3(ParSwarm(ParSwarmRow,1),ParSwarm(ParSwarmRow,2),ParSwarm(ParSwarmRow,5),'r.','markersize',8);
end
end
end
XResult=OptSwarm(SwarmSize+1,1:ParticleSize);
YResult=AdaptFunc(XResult);
if IsStep~=0
XResult=OptSwarm(SwarmSize+1,1:ParticleSize);
YResult=AdaptFunc(XResult);
str=sprintf('%g步迭代的最優目標函數值%g',k,YResult);
disp(str);
disp('下次迭代，按隨意鍵繼續');
pause
end
%記錄每一步的平均適應度
MeanAdapt(1,k)=mean(ParSwarm(:,2*ParticleSize+1));
end
%for循環結束標志
%記錄最小與最大的平均適應度
MinMaxMeanAdapt=[min(MeanAdapt),max(MeanAdapt)];
%計算離線與在線性能
for k=1:LoopCount
OnLine(1,k)=sum(MeanAdapt(1,1:k))/k;
OffLine(1,k)=max(MeanAdapt(1,1:k));
end
for k=1:LoopCount
OffLine(1,k)=sum(OffLine(1,1:k))/k;
end
%繪制離線性能與在線性能曲線
if 1==IsPlot
figure
hold on
title('離線性能曲線圖')
xlabel('迭代次數');
ylabel('離線性能');
grid on
plot(OffLine);
figure
hold on
title('在線性能曲線圖')
xlabel('迭代次數');
ylabel('在線性能');
grid on
plot(OnLine);
end
%記錄本次迭代得到的最優結果
XResult=OptSwarm(SwarmSize+1,1:ParticleSize);
YResult=AdaptFunc(XResult);
Result=[XResult,YResult];
這里給出一個使用的樣例代碼，並分別解釋各參數的含義：
%打開計時器
tic;
%
Scope=[-50 50
-50 50
-50 50
-50 50
-50 50
-50 50
-50 50
-50 50
-50 50
-50 50];
[v,on,off,minmax]=PsoProcess(20,10,Scope,@initswarm,@BasestepPSO,@Griewank,0,0,4000,0);
toc
在上面的代碼中函數PsoProcess中的20代表粒子群的規模為20個，10代表每一個粒子的維數為10，Scope是粒子的每一維的范圍，同一時候也是速度的范圍，@initswarm是初始化函數的句柄，@BasestepPSO是單步更新的函數句柄，@Griewank是適應度評價函數的句柄，4000代表真個算法循環4000次終止，其它參數參見說明文檔。
粒子群算法(6)-----幾個適應度評價函數
以下給出幾個適應度評價函數，並給出圖形表示
頭幾天機子種了病毒，又一次安裝了系統，不小心把程序全部格式化了，痛哭！！！沒辦法，好多程序不見了，如今把這幾個典型的函數又一次編寫了，把他們給出來，就算粒子群算法的一個結束吧！痛恨病毒！！！！
第一個函數：Griewank函數，圖形例如以下所看到的：
適應度函數例如以下：（為了求最大值，我去了全部函數值的相反數）
function y=Griewank(x)
%Griewan函數
%輸入x,給出相應的y值,在x=(0,0,…,0)處有全局極小點0.
%編制人：
%編制日期：
[row,col]=size(x);
if row>1
error('輸入的參數錯誤');
end
y1=1/4000*sum(x.^2);
y2=1;
for h=1:col
y2=y2*cos(x(h)/sqrt(h));
end
y=y1-y2+1;
y=-y;
繪制函數圖像的代碼例如以下：
function DrawGriewank()
%繪制Griewank函數圖形
x=[-8:0.1:8];
y=x;
[X,Y]=meshgrid(x,y);
[row,col]=size(X);
for l=1:col
for h=1:row
z(h,l)=Griewank([X(h,l),Y(h,l)]);
end
end
surf(X,Y,z);
shading interp
第二個函數：Rastrigin函數，圖形例如以下所看到的：
適應度函數例如以下：（為了求最大值，我去了全部函數值的相反數）
function y=Rastrigin(x)
%Rastrigin函數
%輸入x,給出相應的y值,在x=(0,0,…,0)處有全局極小點0.
%編制人：
%編制日期：
[row,col]=size(x);
if row>1
error('輸入的參數錯誤');
end
y=sum(x.^2-10*cos(2*pi*x)+10);
y=-y;
繪制函數圖像的代碼例如以下：
function DrawRastrigin()
%繪制Rastrigin函數圖形
x=[-5:0.05:5];
y=x;
[X,Y]=meshgrid(x,y);
[row,col]=size(X);
for l=1:col
for h=1:row
z(h,l)=Rastrigin([X(h,l),Y(h,l)]);
end
end
surf(X,Y,z);
shading interp
這樣粒子群算法不得不草草收場。