一、簡介
粒子群優(yōu)化(PSO)是一種基于群體智能的數(shù)值優(yōu)化算法，由社會心理學(xué)家James Kennedy和電氣工程師Russell Eberhart于1995年提出。自PSO誕生以來，它在許多方面都得到了改進(jìn)，這一部分將介紹基本的粒子群優(yōu)化算法原理和過程。
1.1 粒子群優(yōu)化
粒子群優(yōu)化(PSO)是一種群智能算法，其靈感來自于鳥類的群集或魚群學(xué)習(xí)，用于解決許多科學(xué)和工程領(lǐng)域中出現(xiàn)的非線性、非凸性或組合優(yōu)化問題。

1.1.1 算法思想
許多鳥類都是群居性的，并由各種原因形成不同的鳥群。鳥群可能大小不同，出現(xiàn)在不同的季節(jié)，甚至可能由群體中可以很好合作的不同物種組成。更多的眼睛和耳朵意味著有更多的及時(shí)發(fā)現(xiàn)食物和捕食者的機(jī)會。鳥群在許多方面對其成員的生存總是有益的：
覓食：社會生物學(xué)家E.O.Wilson說，至少在理論上，群體中的個(gè)體成員可以從其他成員在尋找食物過程中的發(fā)現(xiàn)和先前的經(jīng)驗(yàn)中獲益[1]。如果一群鳥的食物來源是相同的，那么某些種類的鳥就會以一種非競爭的方式聚集在一起。這樣，更多的鳥類就能利用其他鳥類對食物位置的發(fā)現(xiàn)。
抵御捕食者：鳥群在保護(hù)自己免受捕食者侵害方面有很多優(yōu)勢。
更多的耳朵和眼睛意味著更多的機(jī)會發(fā)現(xiàn)捕食者或任何其他潛在的危險(xiǎn)；
一群鳥可能會通過圍攻或敏捷的飛行來迷惑或壓制捕食者；
在群體中，互相間的警告可以減少任何一只鳥的危險(xiǎn)。
空氣動力學(xué)：當(dāng)鳥類成群飛行時(shí)，它們經(jīng)常把自己排成特定的形狀或隊(duì)形。鳥群中鳥的數(shù)量不同，每只鳥煽動翅膀時(shí)產(chǎn)生不同的氣流，這都會導(dǎo)致變化的風(fēng)型，這些隊(duì)形會充分利用不同的分型，從而使得飛行中的鳥類能夠以最節(jié)能的方式利用周圍的空氣。
粒子群算法的發(fā)展需要模擬鳥群的一些優(yōu)點(diǎn)，然而，為了了解群體智能和粒子群優(yōu)化的一個(gè)重要性質(zhì)，值得提一下是鳥群的一些缺點(diǎn)。當(dāng)鳥類成群結(jié)隊(duì)時(shí)，也會給它們帶來一些風(fēng)險(xiǎn)。更多的耳朵和眼睛意味著更多的翅膀和嘴，這導(dǎo)致更多的噪音和運(yùn)動。在這種情況下，更多的捕食者可以定位鳥群，對鳥類造成持續(xù)的威脅。一個(gè)更大的群體也會需要更多的食物，這導(dǎo)致更多食物競爭，從而可能淘汰群體中一些較弱的鳥類。這里需要指出的是，PSO并沒有模擬鳥類群體行為的缺點(diǎn)，因此，在搜索過程中不允許殺死任何個(gè)體，而在遺傳算法中，一些較弱的個(gè)體會消亡。在PSO中，所有的個(gè)體都將存活，并在整個(gè)搜索過程中努力讓自己變得更強(qiáng)大。在粒子群算法中，潛在解的改進(jìn)是合作的結(jié)果，而在進(jìn)化算法中則是因?yàn)楦偁?。這個(gè)概念使得群體智能不同于進(jìn)化算法。簡而言之，在進(jìn)化算法中，每一次迭代都有一個(gè)新的種群進(jìn)化，而在群智能算法中，每一代都有個(gè)體使自己變得更好。個(gè)體的身份不會隨著迭代而改變。Mataric[2]給出了以下鳥群規(guī)則：

安全漫游：鳥類飛行時(shí)，不存在相互間或與障礙物間的碰撞；
分散：每只鳥都會與其他鳥保持一個(gè)最小的距離；
聚合：每只鳥也會與其他鳥保持一個(gè)最大的距離；
歸巢：所有的鳥類都有可能找到食物來源或巢穴。
在設(shè)計(jì)粒子群算法時(shí)，并沒有采用這四種規(guī)則來模擬鳥類的群體行為。在Kennedy和Eberhart開發(fā)的基本粒子群優(yōu)化模型中，對agent的運(yùn)動不遵循安全漫游和分散規(guī)則。換句話說，在基本粒子群優(yōu)化算法的運(yùn)動過程中，允許粒子群優(yōu)化算法中的代理盡可能地靠近彼此。而聚合和歸巢在粒子群優(yōu)化模型中是有效的。在粒子群算法中，代理必須在特定的區(qū)域內(nèi)飛行，以便與任何其他代理保持最大距離。這就相當(dāng)于在整個(gè)過程中，搜索始終停留在搜索空間的邊界內(nèi)或邊界處。第四個(gè)規(guī)則，歸巢意味著組中的任何代理都可以達(dá)到全局最優(yōu)。

在PSO模型的發(fā)展過程中，Kennedy和Eberhart提出了五個(gè)判斷一組代理是否是群體的基本原則：

就近原則：代理群體應(yīng)該能夠進(jìn)行簡單的空間和時(shí)間計(jì)算；
質(zhì)量原則：代理群體能夠?qū)Νh(huán)境中的質(zhì)量因素作出反應(yīng)；
多響應(yīng)原則：代理群體不應(yīng)在過于狹窄的通道從事活動；
穩(wěn)定性原則：代理群體不能每次環(huán)境變化時(shí)就改變其行為模式；
適應(yīng)性原則：計(jì)算代價(jià)不大時(shí)，代理群體可以改變其行為模式。

1.1.2 粒子群優(yōu)化過程
考慮到這五個(gè)原則，Kennedy和Eberhart開發(fā)了一個(gè)用于函數(shù)優(yōu)化的PSO模型。在粒子群算法中，采用隨機(jī)搜索的方法，利用群體智能進(jìn)行求解。換句話說，粒子群算法是一種群智能搜索算法。這個(gè)搜索是由一組隨機(jī)生成的可能解來完成的。這種可能解的集合稱為群，每個(gè)可能解都稱為粒子。
在粒子群優(yōu)化算法中，粒子的搜索受到兩種學(xué)習(xí)方式的影響。每一個(gè)粒子都在向其他粒子學(xué)習(xí)，同時(shí)也在運(yùn)動過程中學(xué)習(xí)自己的經(jīng)驗(yàn)。向他人學(xué)習(xí)可以稱為社會學(xué)習(xí)，而從自身經(jīng)驗(yàn)中學(xué)習(xí)可以稱為認(rèn)知學(xué)習(xí)。由于社會學(xué)習(xí)的結(jié)果，粒子在它的記憶中存儲了群中所有粒子訪問的最佳解，我們稱之為gbest。通過認(rèn)知學(xué)習(xí)，粒子在它的記憶中儲存了迄今為止它自己訪問過的最佳解，稱為pbest。

任何粒子的方向和大小的變化都是由一個(gè)叫做速度的因素決定的，速度是位置相對于時(shí)間的變化率。對于PSO，迭代的是時(shí)間。這樣，對于粒子群算法，速度可以定義為位置相對于迭代的變化率。由于迭代計(jì)數(shù)器單位增加，速度v的維數(shù)與位置x相同。

對于D維搜索空間，在時(shí)間步t下群體中的第ith個(gè)粒子由D維向量x i t = ( x i 1 t , ? ??, x i D t ) T x_i^t = {(x_{i1}^t, \cdots ,x_{iD}t)T}xit=(xi1t,?,xiDt)T來表示，其速度由另一個(gè)D維向量v i t = ( v i 1 t , ? ??, v i D t ) T v_i^t = {(v_{i1}^t, \cdots ,v_{iD}t)T}vit=(vi1t,?,viDt)T表示。第ith個(gè)粒子訪問過的最優(yōu)解位置用p i t = ( p i 1 t , ? ??, p i D t ) T p_i^t = {\left( {p_{i1}^t, \cdots ,p_{iD}^t} \right)^T}pit=(pi1t,?,piDt)T表示，群體中最優(yōu)粒子的索引為“g”。第ith個(gè)粒子的速度和位置分別由下式進(jìn)行更新：
v i d t + 1 = v i d t + c 1 r 1 ( p i d t ? x i d t ) + c 2 r 2 ( p g d t ? x i d t ) (1) v_{id}^{t + 1} = v_{id}^t + {c_1}{r_1}\left( {p_{id}^t - x_{id}^t} \right) + {c_2}{r_2}\left( {p_{gd}^t - x_{id}^t} \right)\tag 1vidt+1=vidt+c1r1(pidt?xidt)+c2r2(pgdt?xidt)(1)

x i d t + 1 = x i d t + v i d t + 1 (2) x_{id}^{t + 1} = x_{id}^t + v_{id}^{t + 1}\tag 2xidt+1=xidt+vidt+1(2)

其中d=1,2,…,D為維度，i=1,2,…,S為粒子索引，S是群體大小。c1和c2為常數(shù)，分別稱為認(rèn)知和社交縮放參數(shù)，或簡單地稱為加速系數(shù)。r1和r2是滿足均勻分布[0,1]之間的隨機(jī)數(shù)。上面兩個(gè)式子均是對每個(gè)粒子的每個(gè)維度進(jìn)行單獨(dú)更新，問題空間中不同維度之間唯一的聯(lián)系是通過目標(biāo)函數(shù)引入的，也就是當(dāng)前所找到的最好位置gbest和pbest[3]。PSO的算法流程如下：
1.1.3 解讀更新等式
速度更新等式（1）的右側(cè)包括三部分3：
前一時(shí)間的速度v，可以認(rèn)為是一動量項(xiàng)，用于存儲之前的運(yùn)動方向，其目的是防止粒子劇烈地改變方向。
第二項(xiàng)是認(rèn)知或自我部分，通過這一項(xiàng)，粒子的當(dāng)前位置會向其自己的最好位置移動，這樣在整個(gè)搜索過程中，粒子會記住自己的最佳位置，從而避免自己四處游蕩。這里需要注意的是，pidt-xidt是一個(gè)方向從xidt到pidt的向量，從而將當(dāng)前位置向粒子的最佳位置吸引，兩者的順序不能改變，否則當(dāng)前位置會遠(yuǎn)離最佳位置。
第三項(xiàng)是社交部分，負(fù)責(zé)通過群體共享信息。通過該項(xiàng)，粒子向群體中最優(yōu)的個(gè)體移動，即每個(gè)個(gè)體向群體中的其他個(gè)體學(xué)習(xí)。同樣兩者應(yīng)該是pgbt-xidt。
可以看出，認(rèn)知尺度參數(shù)c1調(diào)節(jié)的是粒子在其最佳位置方向上的最大步長，而社交尺度參數(shù)c2調(diào)節(jié)的是全局最優(yōu)粒子方向上的最大步長。圖2給出了粒子在二維空間中運(yùn)動的典型幾何圖形。

圖2 粒子群優(yōu)化過程中粒子移動的幾何說明

從更新方程可以看出，Kennedy和Eberhart的PSO設(shè)計(jì)遵循了PSO的五個(gè)基本原則。在粒子群優(yōu)化過程中，在d維空間中對一系列時(shí)間步進(jìn)行計(jì)算。在任何時(shí)間步，種群都遵循gbest和pbest的指導(dǎo)方向，即種群對質(zhì)量因素作出反應(yīng)，從而遵循質(zhì)量原則。由于速度更新方程中有均布隨機(jī)數(shù)r1和r2，在pbest和gbest之間的當(dāng)前位置隨機(jī)分配，這證明了響應(yīng)原理的多樣性。在粒子群優(yōu)化過程中，只有當(dāng)粒子群從gbest中接收到較好的信息時(shí)，才會發(fā)生隨機(jī)運(yùn)動，從而證明了粒子群優(yōu)化過程的穩(wěn)定性原則。種群在gbest變化時(shí)發(fā)生變化，因此遵循適應(yīng)性原則。
1.2 粒子群優(yōu)化中的參數(shù)
任何基于種群的算法的收斂速度和尋優(yōu)能力都受其參數(shù)選擇的影響。通常，由于這些算法的參數(shù)高度依賴于問題參數(shù)，因此不可能對這些算法的參數(shù)設(shè)置給出一般性的建議。但是，已有的理論和/或?qū)嶒?yàn)研究，給出了參數(shù)值的一般范圍。與其他基于種群的搜索算法類似，由于在搜索過程中存在隨機(jī)因素r1和r2，因此通用PSO的參數(shù)調(diào)整一直是一項(xiàng)具有挑戰(zhàn)性的任務(wù)。PSO的基礎(chǔ)版本只需要很少的參數(shù)。本章只討論了[4]中介紹的PSO基礎(chǔ)版本的參數(shù)。

一個(gè)基本的參數(shù)是群體規(guī)模，它通常是根據(jù)問題中決策變量的數(shù)量和問題的復(fù)雜性經(jīng)驗(yàn)地設(shè)置的。一般建議20-50個(gè)粒子。

另一個(gè)參數(shù)是縮放因子c1和c2。如前所述，這些參數(shù)決定了下一個(gè)迭代中粒子的步長。也就是說，c1和c2決定了粒子的速度。在PSO的基礎(chǔ)版本中，選擇c1=c2=2。在這種情況下，粒子s速度的增加是不受控制的，這有利于更快的收斂速度，但不利于更好地利用搜索空間。如果我們令c1=c2>0，那么粒子會吸引到pbest和gbest的平均值。c1>c2設(shè)置有利于多模態(tài)問題，而c2>c1有利于單模態(tài)問題。在搜索過程中，c1和c2的值越小，粒子軌跡越平滑，而c1和c2的值越大，粒子運(yùn)動越劇烈，加速度越大。研究人員也提出了自適應(yīng)加速系數(shù)[5]。
停止準(zhǔn)則不僅是粒子群算法的參數(shù)，也是任何基于種群的元啟發(fā)式算法的參數(shù)。常用的停止準(zhǔn)則通?；诤瘮?shù)評估或迭代的最大次數(shù)，該次數(shù)與算法所花費(fèi)的時(shí)間成正比。一個(gè)更有效的停止準(zhǔn)則是基于算法的搜索能力，如果一個(gè)算法在一定的迭代次數(shù)內(nèi)沒有顯著地改進(jìn)解，那么應(yīng)該停止搜索。

二、源代碼

%========================================================================== clear all; %這里是參數(shù)的預(yù)備部分，調(diào)試時(shí)就調(diào)整這個(gè)里面的參數(shù)即可preparation itmax=400;%進(jìn)化代數(shù)，就是預(yù)設(shè)的迭代次數(shù)。 W(1)=0.729;% 粒子先前速度保持。慣性權(quán)重 a(1)=0.316;% 用于計(jì)算W。 c1=2; %1.45449; %認(rèn)知部分 加速系數(shù) c2=2; %1.45449;% 社會部分 加速系數(shù) m=3; % 待選擇DC的個(gè)數(shù)（一定要知道，這是基數(shù)，不是序數(shù)，即不是需求點(diǎn)的標(biāo)號） M=[13 13 13];%各DC的容量均為13個(gè)單位，M是m維向量，此處就是3維了。%要×1000的。 R=1000000000; ?%罰系數(shù)，要求充分大即可 h=[ 0 ?1 ?6 ?7 ?4 ?3 ?4 ?6 ?6 ?9 ?8 ? 9; ? ?1 ?0 ?5 ?6 ?5 ?4 ?5 ?7 ?7 10 ?9 ?10; ? ?6 ?5 ?0 ?3 ?6 ?9 10 12 12 15 14 ?15; ? ?7 ?6 ?3 ?0 ?3 10 11 13 13 16 15 ?12; ? ?4 ?5 ?6 ?3 ?0 ?7 ?8 10 10 13 12 ? 9; ? ?3 ?4 ?9 10 ?7 ?0 ?6 ?4 ?9 10 ?6 ? 6; ? ?4 ?5 10 11 ?8 ?6 ?0 ?2 ?9 ?5 ?4 ? 9; ? ?6 ?7 12 13 10 ?4 ?2 ?0 10 ?6 ?2 ? 7; ? ?6 ?7 12 13 10 ?9 ?9 10 ?0 ?4 ?8 ?13; ? ?9 10 15 16 13 10 ?5 ?6 ?4 ?0 ?4 ? 9; ? ?8 ?9 14 15 12 ?6 ?4 ?2 ?8 ?4 ?0 ? 5; ? ?9 10 15 12 ?9 ?6 ?9 ?7 13 ?9 ?5 ? 0]; ?% h必須是對稱矩陣。 %h矩陣表示單位運(yùn)輸量的成本，第（i,j）個(gè)元素指從第i個(gè)需求點(diǎn)到第j個(gè)需求點(diǎn)的單位運(yùn)輸量的成本。 dem=[5 4 2 3 2 4 3 5 4 3 2 2];% dem 原來是demand ，需求量 F=[11 16 14 14 15 13 18 12 11 14 16 11]; %固定費(fèi)用 num_particle=100; % number of particle, should be more than D=length(h); % 需求點(diǎn)個(gè)數(shù) %********************************************************** %下面生成隨機(jī)初始值x,假設(shè)初始的服從[1,3]均勻分布,即[1,m]. %初始v 在[-2,2]之間的整數(shù)，即[-(m-1),m-1]. xmax=m; xmin=1; IndexIteration=1; tic; particle=zeros(2*num_particle,D,itmax); particle(:,:,1)=xmin+(xmax-xmin)*rand(2*num_particle,D); V(:,:,1)=round((xmin-xmax)+2*(xmax-xmin)*rand(2*num_particle,D)); %********************************************************** fit=zeros(num_particle,itmax);% 用于存儲粒子的適應(yīng)值 pbest=zeros(2*num_particle,D,itmax); % 用于存儲粒子的位置 num=numbering(m,D,2*num_particle,particle(:,:,1));% 粒子的編號及對應(yīng)關(guān)系 fit(:,1)=fitness(h,F,m,D,2*num_particle,dem,num,M,R); %********************************************************* pbest(:,:,1)=particle(:,:,1); pbest_value(:,1)=fit(:,1); ?%個(gè)體最優(yōu)值 [C,I]=min(pbest_value(:,1)); gbest_value(1)=C; % 群最優(yōu)值 for i=1:num_particle gbest(2*i-1:2*i,:,1)=particle(2*I-1:2*I,:,1); ?% 群最優(yōu)粒子位置 end %第一次迭代中最優(yōu)就是本身的結(jié)果，因?yàn)榫鸵粋€(gè)嘛，沒有選擇 %****************************************************** for IndexIteration=2:itmax %V(:,:,IndexIteration)=W(IndexIteration-1)*V(:,:,IndexIteration-1)+c1*rand*(pbest(:,:,IndexIteration-1)... % ? -particle(:,:,IndexIteration-1))+... %c2*rand*(gbest(:,:,IndexIteration-1)-particle(:,:,IndexIteration-1)); V(:,:,IndexIteration)=0.729*V(:,:,IndexIteration-1)+c1*rand*(pbest(:,:,IndexIteration-1)... ? ?-particle(:,:,IndexIteration-1))+... c2*rand*(gbest(:,:,IndexIteration-1)-particle(:,:,IndexIteration-1)); V(:,:,IndexIteration)=min(V(:,:,IndexIteration),xmax-xmin);% V(:,:,IndexIteration)=max(V(:,:,IndexIteration),xmin-xmax);% 所有粒子的速度均在[-(m-1),m-1]內(nèi) particle(:,:,IndexIteration)=particle(:,:,IndexIteration-1)+V(:,:,IndexIteration); particle(:,:,IndexIteration)=min(particle(:,:,IndexIteration),xmax);% particle(:,:,IndexIteration)=max(particle(:,:,IndexIteration),xmin);% 所有粒子均在[1,m]內(nèi) %for i=1:num_particle %V(2*i,:,IndexIteration)=min(V(2*i,:,IndexIteration),xmax-xmin);% %V(2*i,:,IndexIteration)=max(V(2*i,:,IndexIteration),xmin-xmax);% 所有粒子的速度均在[-(m-1),m-1]內(nèi) %end %particle(:,:,IndexIteration)=particle(:,:,IndexIteration-1)+V(:,:,IndexIteration); %for i=1:num_particle %particle(2*i,:,IndexIteration)=min(particle(2*i,:,IndexIteration),xmax);% %particle(2*i,:,IndexIteration)=max(particle(2*i,:,IndexIteration),xmin);% %所有粒子均在[1,m]內(nèi) num=numbering(m,D,2*num_particle,particle(:,:,IndexIteration));% 粒子的編號及對應(yīng)關(guān)系 fit(:,IndexIteration)=fitness(h,F,m,D,2*num_particle,dem,num,M,R); %下面更新 pbest and pbest_value pbest_value(:,IndexIteration)=min(pbest_value(:,IndexIteration-1),fit(:,IndexIteration)); for i=1:num_particle ? ?if pbest_value(i,IndexIteration)==fit(i,IndexIteration) ? ? ? ? ?pbest(2*i-1:2*i,:,IndexIteration)=particle(2*i-1:2*i,:,IndexIteration); ? else ? ? ? pbest(2*i-1:2*i,:,IndexIteration)=pbest(2*i-1:2*i,:,IndexIteration-1); ? ?end end function y=numbering(m,D,n,x); [xs,xid]=sort(x'); xz=floor(x')+1; xx=zeros(D,n); for i=1:D ? ?for j=1:2:n ? ? ? ?if i<=D-m ? ? ? ?xx(xid(i,j),j)=0; ? ? ? ?elseif i>D-m ? ? ?xx(xid(i,j),j)=i-(D-m); ? ? ? ?end ?% x是按列看的，最后要轉(zhuǎn)置的。 ? ?end ? ?for j=2:2:n ? ? ? if xz(i,j)>m ? ? ? ? ? xx(i,j)=m; ? ? ? elseif xz(i,j)<1 ? ? ? ? ? xx(i,j)=1; ? ? ? else ? ? ? ? ? ?xx(i,j)=xz(i,j); ? ? ? end ? ?end end

三、運(yùn)行結(jié)果

?

?