一、简介
粒子群算法源于复杂适应系统（Complex Adaptive System,CAS）。CAS理论于1994年正式提出，CAS中的成员称为主体。比如研究鸟群系统，每个鸟在这个系统中就称为主体。主体有适应性，它能够与环境及其他的主体进行交流，并且根据交流的过程“学习”或“积累经验”改变自身结构与行为。整个系统的演变或进化包括：新层次的产生（小鸟的出生）；分化和多样性的出现（鸟群中的鸟分成许多小的群）；新的主题的出现（鸟寻找食物过程中，不断发现新的食物）。
所以CAS系统中的主体具有4个基本特点（这些特点是粒子群算法发展变化的依据）：
首先，主体是主动的、活动的。
主体与环境及其他主体是相互影响、相互作用的，这种影响是系统发展变化的主要动力。
环境的影响是宏观的，主体之间的影响是微观的，宏观与微观要有机结合。
最后，整个系统可能还要受一些随机因素的影响。
粒子群算法就是对一个CAS系统－－－鸟群社会系统的研究得出的。
粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究。设想这样一个场景:一群鸟在随机搜寻食物，在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢?最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
PSO算法就从这种生物种群行为特性中得到启发并用于求解优化问题。在PSO中，每个优化问题的潜在解都可以想象成d维搜索空间上的一个点，我们称之为“粒子”（Particle），所有的粒子都有一个被目标函数决定的适应值(Fitness Value )，每个粒子还有一个速度决定他们飞翔的方向和距离，然后粒子们就追随当前的最优粒子在解空间中搜索。Reynolds对鸟群飞行的研究发现。鸟仅仅是追踪它有限数量的邻居但最终的整体结果是整个鸟群好像在一个中心的控制之下.即复杂的全局行为是由简单规则的相互作用引起的。
二、源代码
%% 基于粒子群工具箱的函数优化算法
%% 清空环境
clear
clc
%% 参数初始化
x_range=[-50,50]; %参数x变化范围
y_range=[-50,50]; %参数y变化范围
range = [x_range;y_range]; %参数变化范围(组成矩阵)
Max_V = 0.2*(range(:,2)-range(:,1)); %最大速度取变化范围的10%~20%
n=2; %待优化函数的维数，此例子中仅x、y两个自变量，故为2
PSOparams= [25 2000 24 2 2 0.9 0.4 1500 1e-25 250 NaN 0 0];
%% 粒子群寻优
pso_Trelea_vectorized('Rosenbrock',n,Max_V,range,0,PSOparams) %调用PSO核心模块
% goplotpso.m
% default plotting script used in PSO functions
%
% this script is not a function,
% it is a plugin for the main PSO routine (pso_Trelea_vectorized)
% so it shares all the same variables, be careful with variable names
% when making your own plugin
% setup figure, change this for your own machine
clf
set(gcf,'Position',[651 31 626 474]); % this is the computer dependent part
%set(gcf,'Position',[743 33 853 492]);
set(gcf,'Doublebuffer','on');
% particle plot, upper right
subplot('position',[.7,.6,.27,.32]);
set(gcf,'color','k')
plot3(pos(:,1),pos(:,D),out,'b.','Markersize',7)
hold on
plot3(pbest(:,1),pbest(:,D),pbestval,'g.','Markersize',7);
plot3(gbest(1),gbest(D),gbestval,'r.','Markersize',25);
% crosshairs
offx = max(abs(min(min(pbest(:,1)),min(pos(:,1)))),...
abs(max(max(pbest(:,1)),max(pos(:,1)))));
offy = max(abs(min(min(pbest(:,D)),min(pos(:,D)))),...
abs(min(max(pbest(:,D)),max(pos(:,D)))));
plot3([gbest(1)-offx;gbest(1)+offx],...
[gbest(D);gbest(D)],...
[gbestval;gbestval],...
'r-.');
plot3([gbest(1);gbest(1)],...
[gbest(D)-offy;gbest(D)+offy],...
[gbestval;gbestval],...
'r-.');
hold off
xlabel('Dimension 1','color','y')
ylabel(['Dimension ',num2str(D)],'color','y')
zlabel('Cost','color','y')
title('Particle Dynamics','color','w','fontweight','bold')
set(gca,'Xcolor','y')
set(gca,'Ycolor','y')
set(gca,'Zcolor','y')
set(gca,'color','k')
% camera control
view(2)
try
axis([gbest(1)-offx,gbest(1)+offx,gbest(D)-offy,gbest(D)+offy]);
catch
axis([VR(1,1),VR(1,2),VR(D,1),VR(D,2)]);
end
% error plot, left side
subplot('position',[0.1,0.1,.475,.825]);
semilogy(tr(find(~isnan(tr))),'color','m','linewidth',2)
%plot(tr(find(~isnan(tr))),'color','m','linewidth',2)
xlabel('epoch','color','y')
ylabel('gbest val.','color','y')
if D==1
titstr1=sprintf(['%11.6g = %s( [ %9.6g ] )'],...
gbestval,strrep(functname,'_','\_'),gbest(1));
elseif D==2
titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g ] )'],...
gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2));
elseif D==3
titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g, %9.6g ] )'],...
gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2),gbest(3));
else
titstr1=sprintf(['%11.6g = %s( [ %g inputs ] )'],...
gbestval,strrep(functname,'_','\_'),D);
end
title(titstr1,'color','m','fontweight','bold');
grid on
% axis tight
set(gca,'Xcolor','y')
set(gca,'Ycolor','y')
set(gca,'Zcolor','y')
set(gca,'color','k')
set(gca,'YMinorGrid','off')
% text box in lower right
% doing it this way so I can format each line any way I want
subplot('position',[.62,.1,.29,.4]);
clear titstr
if trelea==0
PSOtype = 'Common PSO';
xtraname = 'Inertia Weight : ';
xtraval = num2str(iwt(length(iwt)));
elseif trelea==2 | trelea==1
PSOtype = (['Trelea Type ',num2str(trelea)]);
xtraname = ' ';
xtraval = ' ';
elseif trelea==3
PSOtype = (['Clerc Type 1"']);
xtraname = '\chi value : ';
xtraval = num2str(chi);
end
if isnan(errgoal)
errgoalstr='Unconstrained';
else
errgoalstr=num2str(errgoal);
end
if minmax==1
minmaxstr = ['Maximize to : '];
elseif minmax==0
minmaxstr = ['Minimize to : '];
else
minmaxstr = ['Target to : '];
end
if rstflg==1
rststat1 = 'Environment Change';
rststat2 = ' ';
else
rststat1 = ' ';
rststat2 = ' ';
end
titstr={'PSO Model: ' ,PSOtype;...
'Dimensions : ' ,num2str(D);...
'# of particles : ',num2str(ps);...
minmaxstr ,errgoalstr;...
'Function : ' ,strrep(functname,'_','\_');...
xtraname ,xtraval;...
rststat1 ,rststat2};
text(.1,1,[titstr{1,1},titstr{1,2}],'color','g','fontweight','bold');
hold on
text(.1,.9,[titstr{2,1},titstr{2,2}],'color','m');
text(.1,.8,[titstr{3,1},titstr{3,2}],'color','m');
text(.1,.7,[titstr{4,1}],'color','w');
text(.55,.7,[titstr{4,2}],'color','m');
text(.1,.6,[titstr{5,1},titstr{5,2}],'color','m');
text(.1,.5,[titstr{6,1},titstr{6,2}],'color','w','fontweight','bold');
text(.1,.4,[titstr{7,1},titstr{7,2}],'color','r','fontweight','bold');
% if we are training a neural net, show a few more parameters
if strcmp('pso_neteval',functname)
% net is passed from trainpso to pso_Trelea_vectorized in case you are
% wondering where that structure comes from
hiddlyrstr = [];
for lyrcnt=1:length(net.layers)
TF{lyrcnt} = net.layers{lyrcnt}.transferFcn;
Sn(lyrcnt) = net.layers{lyrcnt}.dimensions;
hiddlyrstr = [hiddlyrstr,', ',TF{lyrcnt}];
end
hiddlyrstr = hiddlyrstr(3:end);
text(0.1,.35,['#neur/lyr = [ ',num2str(net.inputs{1}.size),' ',...
num2str(Sn),' ]'],'color','c','fontweight','normal',...
'fontsize',10);
text(0.1,.275,['Lyr Fcn: ',hiddlyrstr],...
'color','c','fontweight','normal','fontsize',9);
end
legstr = {'Green = Personal Bests';...
'Blue = Current Positions';...
'Red = Global Best'};
text(.1,0.025,legstr{1},'color','g');
text(.1,-.05,legstr{2},'color','b');
text(.1,-.125,legstr{3},'color','r');
hold off
set(gca,'color','k');
set(gca,'visible','off');
drawnow
三、运行结果
四、备注
版本：2014a