一、簡介
粒子群優化(PSO)是一種基于群體智能的數值優化算法,由社會心理學家James Kennedy和電氣工程師Russell Eberhart于1995年提出。自PSO誕生以來,它在許多方面都得到了改進,這一部分将介紹基本的粒子群優化算法原理和過程。
1.1 粒子群優化
粒子群優化(PSO)是一種群智能算法,其靈感來自于鳥類的群集或魚群學習,用于解決許多科學和工程領域中出現的非線性、非凸性或組合優化問題。
1.1.1 算法思想
許多鳥類都是群居性的,并由各種原因形成不同的鳥群。鳥群可能大小不同,出現在不同的季節,甚至可能由群體中可以很好合作的不同物種組成。更多的眼睛和耳朵意味着有更多的及時發現食物和捕食者的機會。鳥群在許多方面對其成員的生存總是有益的:
覓食:社會生物學家E.O.Wilson說,至少在理論上,群體中的個體成員可以從其他成員在尋找食物過程中的發現和先前的經驗中獲益[1]。如果一群鳥的食物來源是相同的,那麼某些種類的鳥就會以一種非競争的方式聚集在一起。這樣,更多的鳥類就能利用其他鳥類對食物位置的發現。
抵禦捕食者:鳥群在保護自己免受捕食者侵害方面有很多優勢。
更多的耳朵和眼睛意味着更多的機會發現捕食者或任何其他潛在的危險;
一群鳥可能會通過圍攻或靈活的飛行來迷惑或壓制捕食者;
在群體中,互相間的警告可以減少任何一隻鳥的危險。
空氣動力學:當鳥類成群飛行時,它們經常把自己排成特定的形狀或隊形。鳥群中鳥的數量不同,每隻鳥煽動翅膀時産生不同的氣流,這都會導緻變化的風型,這些隊形會充分利用不同的分型,進而使得飛行中的鳥類能夠以最節能的方式利用周圍的空氣。
粒子群算法的發展需要模拟鳥群的一些優點,然而,為了了解群體智能和粒子群優化的一個重要性質,值得提一下是鳥群的一些缺點。當鳥類成群結隊時,也會給它們帶來一些風險。更多的耳朵和眼睛意味着更多的翅膀和嘴,這導緻更多的噪音和運動。在這種情況下,更多的捕食者可以定位鳥群,對鳥類造成持續的威脅。一個更大的群體也會需要更多的食物,這導緻更多食物競争,進而可能淘汰群體中一些較弱的鳥類。這裡需要指出的是,PSO并沒有模拟鳥類群體行為的缺點,是以,在搜尋過程中不允許殺死任何個體,而在遺傳算法中,一些較弱的個體會消亡。在PSO中,所有的個體都将存活,并在整個搜尋過程中努力讓自己變得更強大。在粒子群算法中,潛在解的改進是合作的結果,而在進化算法中則是因為競争。這個概念使得群體智能不同于進化算法。簡而言之,在進化算法中,每一次疊代都有一個新的種群進化,而在群智能算法中,每一代都有個體使自己變得更好。個體的身份不會随着疊代而改變。Mataric[2]給出了以下鳥群規則:
安全漫遊:鳥類飛行時,不存在互相間或與障礙物間的碰撞;
分散:每隻鳥都會與其他鳥保持一個最小的距離;
聚合:每隻鳥也會與其他鳥保持一個最大的距離;
歸巢:所有的鳥類都有可能找到食物來源或巢穴。
在設計粒子群算法時,并沒有采用這四種規則來模拟鳥類的群體行為。在Kennedy和Eberhart開發的基本粒子群優化模型中,對agent的運動不遵循安全漫遊和分散規則。換句話說,在基本粒子群優化算法的運動過程中,允許粒子群優化算法中的代理盡可能地靠近彼此。而聚合和歸巢在粒子群優化模型中是有效的。在粒子群算法中,代理必須在特定的區域内飛行,以便與任何其他代理保持最大距離。這就相當于在整個過程中,搜尋始終停留在搜尋空間的邊界内或邊界處。第四個規則,歸巢意味着組中的任何代理都可以達到全局最優。
在PSO模型的發展過程中,Kennedy和Eberhart提出了五個判斷一組代理是否是群體的基本原則:
就近原則:代理群體應該能夠進行簡單的空間和時間計算;
品質原則:代理群體能夠對環境中的品質因素作出反應;
多響應原則:代理群體不應在過于狹窄的通道從事活動;
穩定性原則:代理群體不能每次環境變化時就改變其行為模式;
适應性原則:計算代價不大時,代理群體可以改變其行為模式。
1.1.2 粒子群優化過程
考慮到這五個原則,Kennedy和Eberhart開發了一個用于函數優化的PSO模型。在粒子群算法中,采用随機搜尋的方法,利用群體智能進行求解。換句話說,粒子群算法是一種群智能搜尋算法。這個搜尋是由一組随機生成的可能解來完成的。這種可能解的集合稱為群,每個可能解都稱為粒子。
在粒子群優化算法中,粒子的搜尋受到兩種學習方式的影響。每一個粒子都在向其他粒子學習,同時也在運動過程中學習自己的經驗。向他人學習可以稱為社會學習,而從自身經驗中學習可以稱為認知學習。由于社會學習的結果,粒子在它的記憶中存儲了群中所有粒子通路的最佳解,我們稱之為gbest。通過認知學習,粒子在它的記憶中儲存了迄今為止它自己通路過的最佳解,稱為pbest。
任何粒子的方向和大小的變化都是由一個叫做速度的因素決定的,速度是位置相對于時間的變化率。對于PSO,疊代的是時間。這樣,對于粒子群算法,速度可以定義為位置相對于疊代的變化率。由于疊代計數器機關增加,速度v的維數與位置x相同。
對于D維搜尋空間,在時間步t下群體中的第ith個粒子由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來表示,其速度由另一個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個粒子通路過的最優解位置用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表示,群體中最優粒子的索引為“g”。第ith個粒子的速度和位置分别由下式進行更新:
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為常數,分别稱為認知和社交縮放參數,或簡單地稱為加速系數。r1和r2是滿足均勻分布[0,1]之間的随機數。上面兩個式子均是對每個粒子的每個次元進行單獨更新,問題空間中不同次元之間唯一的聯系是通過目标函數引入的,也就是目前所找到的最好位置gbest和pbest[3]。PSO的算法流程如下:
1.1.3 解讀更新等式
速度更新等式(1)的右側包括三部分3:
前一時間的速度v,可以認為是一動量項,用于存儲之前的運動方向,其目的是防止粒子劇烈地改變方向。
第二項是認知或自我部分,通過這一項,粒子的目前位置會向其自己的最好位置移動,這樣在整個搜尋過程中,粒子會記住自己的最佳位置,進而避免自己四處遊蕩。這裡需要注意的是,pidt-xidt是一個方向從xidt到pidt的向量,進而将目前位置向粒子的最佳位置吸引,兩者的順序不能改變,否則目前位置會遠離最佳位置。
第三項是社交部分,負責通過群體共享資訊。通過該項,粒子向群體中最優的個體移動,即每個個體向群體中的其他個體學習。同樣兩者應該是pgbt-xidt。
可以看出,認知尺度參數c1調節的是粒子在其最佳位置方向上的最大步長,而社交尺度參數c2調節的是全局最優粒子方向上的最大步長。圖2給出了粒子在二維空間中運動的典型幾何圖形。
圖2 粒子群優化過程中粒子移動的幾何說明
從更新方程可以看出,Kennedy和Eberhart的PSO設計遵循了PSO的五個基本原則。在粒子群優化過程中,在d維空間中對一系列時間步進行計算。在任何時間步,種群都遵循gbest和pbest的指導方向,即種群對品質因素作出反應,進而遵循品質原則。由于速度更新方程中有均布随機數r1和r2,在pbest和gbest之間的目前位置随機配置設定,這證明了響應原理的多樣性。在粒子群優化過程中,隻有當粒子群從gbest中接收到較好的資訊時,才會發生随機運動,進而證明了粒子群優化過程的穩定性原則。種群在gbest變化時發生變化,是以遵循适應性原則。
1.2 粒子群優化中的參數
任何基于種群的算法的收斂速度和尋優能力都受其參數選擇的影響。通常,由于這些算法的參數高度依賴于問題參數,是以不可能對這些算法的參數設定給出一般性的建議。但是,已有的理論和/或實驗研究,給出了參數值的一般範圍。與其他基于種群的搜尋算法類似,由于在搜尋過程中存在随機因素r1和r2,是以通用PSO的參數調整一直是一項具有挑戰性的任務。PSO的基礎版本隻需要很少的參數。本章隻讨論了[4]中介紹的PSO基礎版本的參數。
一個基本的參數是群體規模,它通常是根據問題中決策變量的數量和問題的複雜性經驗地設定的。一般建議20-50個粒子。
另一個參數是縮放因子c1和c2。如前所述,這些參數決定了下一個疊代中粒子的步長。也就是說,c1和c2決定了粒子的速度。在PSO的基礎版本中,選擇c1=c2=2。在這種情況下,粒子s速度的增加是不受控制的,這有利于更快的收斂速度,但不利于更好地利用搜尋空間。如果我們令c1=c2>0,那麼粒子會吸引到pbest和gbest的平均值。c1>c2設定有利于多模态問題,而c2>c1有利于單模态問題。在搜尋過程中,c1和c2的值越小,粒子軌迹越平滑,而c1和c2的值越大,粒子運動越劇烈,加速度越大。研究人員也提出了自适應加速系數[5]。
停止準則不僅是粒子群算法的參數,也是任何基于種群的元啟發式算法的參數。常用的停止準則通常基于函數評估或疊代的最大次數,該次數與算法所花費的時間成正比。一個更有效的停止準則是基于算法的搜尋能力,如果一個算法在一定的疊代次數内沒有顯著地改進解,那麼應該停止搜尋。
二、源代碼
% pso_Trelea_vectorized.m
% a generic particle swarm optimizer
% to find the minimum or maximum of any
% MISO matlab function
%
% Implements Common, Trelea type 1 and 2, and Clerc's class 1". It will
% also automatically try to track to a changing environment (with varied
% success - BKB 3/18/05)
%
% This vectorized version removes the for loop associated with particle
% number. It also *requires* that the cost function have a single input
% that represents all dimensions of search (i.e., for a function that has 2
% inputs then make a wrapper that passes a matrix of ps x 2 as a single
% variable)
%
% Usage:
% [optOUT]=PSO(functname,D)
% or:
% [optOUT,tr,te]=...
% PSO(functname,D,mv,VarRange,minmax,PSOparams,plotfcn,PSOseedValue)
%
% Inputs:
% functname - string of matlab function to optimize
% D - # of inputs to the function (dimension of problem)
%
% Optional Inputs:
% mv - max particle velocity, either a scalar or a vector of length D
% (this allows each component to have it's own max velocity),
% default = 4, set if not input or input as NaN
%
% VarRange - matrix of ranges for each input variable,
% default -100 to 100, of form:
% [ min1 max1
% min2 max2
% ...
% minD maxD ]
%
% minmax = 0, funct minimized (default)
% = 1, funct maximized
% = 2, funct is targeted to P(12) (minimizes distance to errgoal)
% PSOparams - PSO parameters
% P(1) - Epochs between updating display, default = 100. if 0,
% no display
% P(2) - Maximum number of iterations (epochs) to train, default = 2000.
% P(3) - population size, default = 24
%
% P(4) - acceleration const 1 (local best influence), default = 2
% P(5) - acceleration const 2 (global best influence), default = 2
% P(6) - Initial inertia weight, default = 0.9
% P(7) - Final inertia weight, default = 0.4
% P(8) - Epoch when inertial weight at final value, default = 1500
% P(9)- minimum global error gradient,
% if abs(Gbest(i+1)-Gbest(i)) < gradient over
% certain length of epochs, terminate run, default = 1e-25
% P(10)- epochs before error gradient criterion terminates run,
% default = 150, if the SSE does not change over 250 epochs
% then exit
% P(11)- error goal, if NaN then unconstrained min or max, default=NaN
% P(12)- type flag (which kind of PSO to use)
% 0 = Common PSO w/intertia (default)
% 1,2 = Trelea types 1,2
% 3 = Clerc's Constricted PSO, Type 1"
% P(13)- PSOseed, default=0
% = 0 for initial positions all random
% = 1 for initial particles as user input
%
% plotfcn - optional name of plotting function, default 'goplotpso',
% make your own and put here
%
% PSOseedValue - initial particle position, depends on P(13), must be
% set if P(13) is 1 or 2, not used for P(13)=0, needs to
% be nXm where n<=ps, and m<=D
% If n<ps and/or m<D then remaining values are set random
% on Varrange
% Outputs:
% optOUT - optimal inputs and associated min/max output of function, of form:
% [ bestin1
% bestin2
% ...
% bestinD
% bestOUT ]
%
% Optional Outputs:
% tr - Gbest at every iteration, traces flight of swarm
% te - epochs to train, returned as a vector 1:endepoch
%
% Example: out=pso_Trelea_vectorized('f6',2)
% Brian Birge
% Rev 3.3
% 2/18/06
function [OUT,varargout]=pso_Trelea_vectorized(functname,D,varargin)
rand('state',sum(100*clock));
if nargin < 2
error('Not enough arguments.');
end
% PSO PARAMETERS
if nargin == 2 % only specified functname and D
VRmin=ones(D,1)*-100;
VRmax=ones(D,1)*100;
VR=[VRmin,VRmax];
minmax = 0;
P = [];
mv = 4;
plotfcn='goplotpso';
elseif nargin == 3 % specified functname, D, and mv
VRmin=ones(D,1)*-100;
VRmax=ones(D,1)*100;
VR=[VRmin,VRmax];
minmax = 0;
mv=varargin{1};
if isnan(mv)
mv=4;
end
P = [];
plotfcn='goplotpso';
elseif nargin == 4 % specified functname, D, mv, Varrange
mv=varargin{1};
if isnan(mv)
mv=4;
end
VR=varargin{2};
minmax = 0;
P = [];
plotfcn='goplotpso';
elseif nargin == 5 % Functname, D, mv, Varrange, and minmax
mv=varargin{1};
if isnan(mv)
mv=4;
end
VR=varargin{2};
minmax=varargin{3};
P = [];
plotfcn='goplotpso';
elseif nargin == 6 % Functname, D, mv, Varrange, minmax, and psoparams
mv=varargin{1};
if isnan(mv)
mv=4;
end
VR=varargin{2};
minmax=varargin{3};
P = varargin{4}; % psoparams
plotfcn='goplotpso';
elseif nargin == 7 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn
mv=varargin{1};
if isnan(mv)
mv=4;
end
VR=varargin{2};
minmax=varargin{3};
P = varargin{4}; % psoparams
plotfcn = varargin{5};
elseif nargin == 8 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn, PSOseedValue
mv=varargin{1};
if isnan(mv)
mv=4;
end
VR=varargin{2};
minmax=varargin{3};
P = varargin{4}; % psoparams
plotfcn = varargin{5};
PSOseedValue = varargin{6};
else
error('Wrong # of input arguments.');
end
% sets up default pso params
Pdef = [100 2000 24 2 2 0.9 0.4 1500 1e-25 250 NaN 0 0];
Plen = length(P);
P = [P,Pdef(Plen+1:end)];
df = P(1);
me = P(2);
ps = P(3);
ac1 = P(4);
ac2 = P(5);
iw1 = P(6);
iw2 = P(7);
iwe = P(8);
ergrd = P(9);
ergrdep = P(10);
errgoal = P(11);
trelea = P(12);
PSOseed = P(13);
% used with trainpso, for neural net training
if strcmp(functname,'pso_neteval')
net = evalin('caller','net');
Pd = evalin('caller','Pd');
Tl = evalin('caller','Tl');
Ai = evalin('caller','Ai');
Q = evalin('caller','Q');
TS = evalin('caller','TS');
end
% error checking
if ((minmax==2) & isnan(errgoal))
error('minmax= 2, errgoal= NaN: choose an error goal or set minmax to 0 or 1');
end
if ( (PSOseed==1) & ~exist('PSOseedValue') )
error('PSOseed flag set but no PSOseedValue was input');
end
if exist('PSOseedValue')
tmpsz=size(PSOseedValue);
if D < tmpsz(2)
error('PSOseedValue column size must be D or less');
end
if ps < tmpsz(1)
error('PSOseedValue row length must be # of particles or less');
end
end
% set plotting flag
if (P(1))~=0
plotflg=1;
else
plotflg=0;
end
% preallocate variables for speed up
tr = ones(1,me)*NaN;
% take care of setting max velocity and position params here
if length(mv)==1
velmaskmin = -mv*ones(ps,D); % min vel, psXD matrix
velmaskmax = mv*ones(ps,D); % max vel
elseif length(mv)==D
velmaskmin = repmat(forcerow(-mv),ps,1); % min vel
velmaskmax = repmat(forcerow( mv),ps,1); % max vel
else
error('Max vel must be either a scalar or same length as prob dimension D');
end
posmaskmin = repmat(VR(1:D,1)',ps,1); % min pos, psXD matrix
posmaskmax = repmat(VR(1:D,2)',ps,1); % max pos
posmaskmeth = 3; % 3=bounce method (see comments below inside epoch loop)
% PLOTTING
message = sprintf('PSO: %%g/%g iterations, GBest = %%20.20g.\n',me);
% INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE
% initialize population of particles and their velocities at time zero,
% format of pos= (particle#, dimension)
% construct random population positions bounded by VR
pos(1:ps,1:D) = normmat(rand([ps,D]),VR',1);
if PSOseed == 1 % initial positions user input, see comments above
tmpsz = size(PSOseedValue);
pos(1:tmpsz(1),1:tmpsz(2)) = PSOseedValue;
end
% construct initial random velocities between -mv,mv
vel(1:ps,1:D) = normmat(rand([ps,D]),...
[forcecol(-mv),forcecol(mv)]',1);
% initial pbest positions vals
pbest = pos;
% VECTORIZE THIS, or at least vectorize cost funct call
out = feval(functname,pos); % returns column of cost values (1 for each particle)
%---------------------------
pbestval=out; % initially, pbest is same as pos
% assign initial gbest here also (gbest and gbestval)
if minmax==1
% this picks gbestval when we want to maximize the function
[gbestval,idx1] = max(pbestval);
elseif minmax==0
% this works for straight minimization
[gbestval,idx1] = min(pbestval);
elseif minmax==2
% this works when you know target but not direction you need to go
% good for a cost function that returns distance to target that can be either
% negative or positive (direction info)
[temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2);
gbestval = pbestval(idx1);
end
% preallocate a variable to keep track of gbest for all iters
bestpos = zeros(me,D+1)*NaN;
gbest = pbest(idx1,:); % this is gbest position
% used with trainpso, for neural net training
% assign gbest to net at each iteration, these interim assignments
% are for plotting mostly
if strcmp(functname,'pso_neteval')
net=setx(net,gbest);
end
%tr(1) = gbestval; % save for output
bestpos(1,1:D) = gbest;
% this part used for implementing Carlisle and Dozier's APSO idea
% slightly modified, this tracks the global best as the sentry whereas
% their's chooses a different point to act as sentry
% see "Tracking Changing Extremea with Adaptive Particle Swarm Optimizer",
% part of the WAC 2002 Proceedings, June 9-13, http://wacong.com
sentryval = gbestval;
sentry = gbest;
if (trelea == 3)
% calculate Clerc's constriction coefficient chi to use in his form
kappa = 1; % standard val = 1, change for more or less constriction
if ( (ac1+ac2) <=4 )
chi = kappa;
else
psi = ac1 + ac2;
chi_den = abs(2-psi-sqrt(psi^2 - 4*psi));
chi_num = 2*kappa;
chi = chi_num/chi_den;
end
end