仿生智能算法系列之六—粒子群算法

粒子群算法(particle swarm optimization，pso)是一种新兴的演化算法, 由Kennedy和Eberhart于1995年用计算机模拟鸟群觅食这一简单的社会行为时，受到启发，简化之后而提出的。

这类算法的仿生基点是：群集动物(如蚂蚁、鸟、鱼等)通过群聚而有效的觅食和逃避追捕。在这类群体的动物中，每个个体的行为是建立在群体行为的基础之上的，即在整个群体中信息是共享的，而且在个体之间存在着信息的交换与协作。如在蚁群中，当每个个体发现食物之后，它将通过接触或化学信号来招募同伴，使整个群落找到食源；在鸟群的飞行中，每只鸟在初始状态下处于随机位置，且朝各个方向随机飞行，但随着时间推移，这些初始处于随机状态的鸟通过相互学习(相互跟踪)组织的聚集成一个个小的群落，并以相同的速度朝着相同的方向飞行，最终整个群落聚集在同一位置──食源。这些群集动物所表现的智能常称为“群体智能”，它可表述为：一组相互之间可以进行直接通讯或间接通讯(通过改变局部环境)的主体，能够通过合作对问题进行分布求解。换言之，一组无智能的主体通过合作表现出智能行为特征。粒子群算法就是以模拟鸟的群集智能为特征，以求解连续变量优化问题为背景的一种优化算法。因其概念简单、参数较少、易于实现等特点，自提出以来已经受到国内外研究者的高度重视并被广泛应用于许多领域。

在算法开始时，随机初始化粒子的位置和速度构成初始种群，初始种群在解空间中为均匀分布。其中第i个粒子在n维解空间的位置和速度可分别表示为Xi=（xi1,xi2,…,xid）和Vi=（vi1,vi2,…,vid），然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己的速度和位置。一个极值是粒子本身到目前为止所找到的最优解，这个极值称为个体极值Pbi=（Pbi1,Pbi2,…,Pbid）。另一个极值是该粒子的邻域到目前为止找到的最优解，这个极值称为整个邻域的最优粒子Nbesti=(Nbesti1,Nbesti2,…,Nbestid)。粒子根据如下的式(1)和式(2)来更新自己的速度和位置：

Vi=Vi+c1•rand()•(Pbesti-Xi)+c2•rand()•(Nbesti-Xi) (1)

Xi= Xi+ Vi (2)

式中c1和c2是加速常量，分别调节向全局最好粒子和个体最好粒子方向飞行的最大步长，若太小，则粒子可能远离目标区域,若太大则会导致突然向目标区域飞去，或飞过目标区域。合适的c1，c2可以加快收敛且不易陷入局部最优。rand()是0到1之间的随机数。粒子在每一维飞行的速度不能超过算法设定的最大速度Vmax。设置较大的Vmax可以保证粒子种群的全局搜索能力，Vmax较小则粒子种群优化算法的局部搜索能力加强。

粒子群优化算法的流程：

粒子群优化算法目前的应用较多，例如：各类连续问题和离散问题的参数优化（包括模糊控制器的设计、机器人路径规划和模式识别等）、 0-1规划问题及带有排序关系的优化问题（包括背包问题、电网机组控制、数据挖掘、tsp问题、vrp问题、job-shop及资源分配等）、系统设计、多目标优化、自动目标检测、时频分析、模式识别、调度、信号处理、决策、机器人应用等，如下图所示的极值优化问题。

粒子群算法的matlab实现

主函数：

function [Result,OnLine,OffLine,MinMaxMeanAdapt]=PSO_Stand(SwarmSize,ParticleSize,ParticleScope,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为本次完整迭代得到的最小与最大的平均适应度

[ParSwarm,OptSwarm]=InitSwarm(SwarmSize,ParticleSize,ParticleScope);

%初始化粒子群

if IsStep~=0

disp(‘开始迭代，按任意键：’)

pause

end

%开始更新算法的调用%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:LoopCount

%显示迭代的次数：

disp(‘———————————————————-‘)

TempStr=sprintf(‘第 %g 次迭代’,k);

disp(TempStr);

disp(‘———————————————————-‘)

%调用一步迭代的算法

[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,ParticleScope,0.95,0.4,LoopCount,k);

figure(1);

plot(ParSwarm(:,1),ParSwarm(:,2),’g*’,’markersize’,8);grid on;

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));

%mean函数为取有效值函数

End

初始化函数：

function [ParSwarm,OptSwarm]=InitSwarm(SwarmSize,ParticleSize,ParticleScope)

%初始化粒子群矩阵，全部设为[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);

%将适应度值最大的粒子的位置最为全局粒子的最优值

每步更新函数：

function [ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,ParticleScope,MaxW,MinW,LoopCount,CurCount)

%输入参数：ParSwarm:粒子群矩阵，包含粒子的位置，速度与当前的目标函数值

%输入参数：OptSwarm：包含粒子群个体最优解与全局最优解的矩阵

%输入参数：ParticleScope:一个粒子在运算中各维的范围；

%输入参数：AdaptFunc：适应度函数

%输入参数：AdaptFunc：适应度函数

%输入参数：MaxW MinW：惯性权重(系数)的最大值与最小值

%输入参数：CurCount：当前迭代的次数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%线形递减策略（惯性权值的变化）

w=MaxW-CurCount*((MaxW-MinW)/LoopCount);

%得到粒子群群体大小以及一个粒子维数的信息

[ParRow,ParCol]=size(ParSwarm);

%得到粒子的维数

ParCol=(ParCol-1)/2;

SubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%求解出历史最优值与当前位置的差值

%*****更改下面的代码，可以更改c1,c2的变化*****

c1=2;

c2=2;

%完成一次粒子位置 速度 最优值的更新迭代

for row=1:ParRow

SubTract2=OptSwarm(ParRow+1,:)-ParSwarm(row,1:ParCol);%计算出全局最优值与当前该粒子位置的差值

%速度更新公式

TempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*SubTract1(row,:)+c2*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=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;

%加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

%寻找适应度函数值最大的解在矩阵中的位置(行数)，进行全局最优值的改变

[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

适应度函数：

function y=AdaptFunc(x)

%针对不同的应用场合，需要自行编写适应度函数，例如最大、最小极值；

End

上述程序中%后为注释，非常详细，有兴趣欢迎交流。