张晓寒杨凌霄唐鹏程——粒子群算法.ppt

粒子群算法及其MATLAB实现,60组 张晓寒 杨凌霄 唐鹏程,PSO算法的基本理论 PSO算法概述 PSO算法源程序案例 PSO算法应用案例,粒子群算法,粒子群算法（particle swarm optimization,简称PSO）是依托群鸟觅食模型（Boid）寻找最优值。 鸟觅食模型需要对鸟的位置和飞翔速度赋予初值，采用对称随机初始化策略，使其在寻优空间范围内随机初始化分布，然后鸟根据自己寻找食物的经验以及鸟群信息共享，不断向目标靠近，经过有限次位移，绝大部分鸟聚集在一起并且找到了食物或离食物非常近。,PSO算法的基本理论,PSO算法最早是由美国电气工程师Eberhart和社会心理学家Kennedy在1995年基于群鸟觅食提出来的。 群鸟觅食过程中，每只鸟的初始状态和飞行方向是随机的，每只鸟都不知道食物在哪里，随着时间推移，处于随机状态的鸟会通过相互学习和自身经验，自行组织成一个群落，向着食物目标前进。,每只鸟能够通过自身经验和外界信息估计当前所处位置对于能找到食物有多大价值，即适应值。 适应值函数：基于某种标准对个体进行评价的函数 每只鸟能够记住自己所找到的最好位置，即局部最优（pbest）；每只鸟还能够记住所有个体能找到的所有最好位置，即全局最优(gbest)。 整个鸟群的觅食中心将向全局最优移动，通过鸟群位置的移动，即迭代过程，可以使鸟群向目标食物逼近。,第i个粒子在D维搜索空间的位置为 粒子的位置可以作为潜在解，将 带入目标函数可以计算出其适应值，根据适应值的大小衡量优劣。 粒子个体经过的最好位置记为 =( , ) 整个群体的所有粒子经过的最好位置记为 () 粒子i的速度记为（) ;为【0,1】之间的随机数，这两个参数是用来保持群体的多样性。 在每次迭代中，粒子根据以下公式更新速度和位置： +( ； 约束因子:控制速度的权重,PSO算法的参数选取,1.粒子数 粒子数一般取值为2040，特殊的难题需要100200个粒子。 粒子数越多，搜索范围越大，越容易找到全局最优解，算法运行时间也越长。 2.惯性因子 惯性因子w对于粒子群算法的收敛性起到很大作用。w值越大，粒子飞翔幅度越大，容易错失局部寻优能力，而全局搜索能力较强;w值越小，则局部寻优能力增强，全局寻优能力减弱。,通过调整w的大小来控制历史的速度对当前速度的影响程度，使其成为兼顾全局搜索和局部搜索的一个折中。惯性因子w的大小决定了对粒子当前速度继承的多少。 如果惯性因子w是变量，通常在迭代开始时将惯性因子w设置的较大，然后在迭代过程中逐步减小。这样可以使粒子群在开始优化时搜索到较大的解空间，得到合适的粒子，然后在后期逐渐收缩到较好的区域进行更精细的搜索，以加快收敛速度和目标精度。,3.加速常数 对于简单的常规问题，一般情况下去2.0。 加速常数调整自身经验和社会经验在其运动中所起作用的权重。若则粒子没有自身经验，只有社会经验，它的收敛速度可能较快，容易陷入局部最优点；若，则粒子群没有共享信息，只有自身经验，得到最优解的几率非常小。若粒子将在没有任何经验和信息的情况下盲目地搜索到有限区域，很难找到最优解。,4.最大飞翔速度 为使粒子有效地进行搜索，需要使用参数对粒子运动速度进行限制。参数有利于防止搜索范围毫无意义的发散，防止粒子群由于飞翔速度过大错过最优目标值。 为了跳出局部最优，需要较大的寻优步长，而在接近最优值时，采用更小的步长会更好。如果的选择是固定不变的，通常设定为每维变化范围的10%20%。,PSO算法程序设计,PSO算法实现的步骤： 1.初始化粒子群（速度和位置）、惯性因子、加速常数、最大迭代次数和算法终止的最小允许误差。 2.评价每个粒子的初始适应值。 3.将初始适应值作为当前全局最优值，将最佳适应值对应的位置作为全局最优值的位置。 4.将最佳初始适应值作为当前全局最优值，将最佳适应值的位置作为全局最优值的位置。 5.根据更新公式更新每个粒子当前的飞行速度。,6.对每个粒子的飞翔速度进行限幅处理，使其不超过设定的最大飞翔速度。 7.根据更新公式更新每个粒子当前所在位置。 8.比较当前粒子的适应值是否比历史局部最优好，如果好，则当前粒子适应值作为粒子局部最优值，其对应的位置作为局部最优值的位置。 9.找出全局最优值，将当前全局最优值的位置作为粒子群的全局最优值所在位置。 10.重复步骤59，直到满足设定的最小误差或达到最大迭代次数。 11.输出粒子群全局最优值和其对应的位置以及每个粒子的局部最优值和其对应的位置。,程序设计流程图,否 是,开始,初始化各参数,计算各粒子适应值,找出个体和群体最优值,更新各个粒子的速度和位置,终止条件？,结束,PSO算法源程序范例,function main() clc; clear all; close all; tic; %程序运行计时 E0=0.001; %允许误差 MaxNum=100; %粒子最大迭代次数 narvs=1; %目标函数的自变量个数 particlesize=30; %粒子群规模 c1=2; %每个粒子的个体学习因子，也称为加速常数 c2=2; %每个粒子的社会学习因子，也成为加速常数 w=0.6; %惯性因子,求解maxf(x)=2.1,vmax=0.8; %粒子的最大飞翔速度 x=-5+10*rand(particlesize,narvs); %粒子所在的位置 v=2*rand(particlesize,narvs); %粒子的飞翔速度 %用inline定义适应度函数以便将子函数文件与主程序文件放在一起 %目标函数是：y=1+(2.1*(1-x+2*x.2).*exp(-x.2/2) %inline命令定义适应度函数如下 fitness=inline(1/(1+(2.1*(1-x+2*x.2).*exp(-x.2/2),x); %inline定义的适应度函数会使程序运行速度大大降低 for i=1:particlesize for j=1:narvs f(i)=fitness(x(i,j); %fitness（适应度函数）用于评价每个粒子的初始适应值 end end,personalbest_x=x; %将初始适应值作为当前全局最优值 personalbest_faval=f; %将各适应值对应的位置作为每个粒子的局部最优值所在的位置 globalbest_faval i=min(personalbest_faval); %将最佳初始适应值对应的位置作为全局最优值所在的位置 globalbest_x=personalbest_x(i,:); %将最佳初始适应值作为当前全局最优值 k=1; while k = MaxNum for i=1:particlesize for j=1:narvs f(i)=fitness(x(i,j) end if f(i)personalbest_faval(i) %判断当前位置是否是历史上最佳位置,personalbest_faval(i)=f(i); personalbest_x(i,:)=x(i,:); end end %反复迭代，计算出局部最优值及其对应的位置 globalbest_faval i=min(personalbest_faval); globalbest_x=personalbest_x(i,:); for i=1:particlesize %更新粒子群里每个个体的最新位置 v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:). +c2*rand*(globalbest_x-x(i,:); for j=1:narvs %判断粒子的飞翔速度是否超过了最大飞翔速度 if v(i,j)vmax; v(i,j)=vmax; elseif v(i,j)-vmax;,v(i,j)=-vmax; %对每一个粒子的飞翔速度进行限幅处理 end end x(i,:)=x(i,:)+v(i,:); %更新粒子当前所在的位置 end if abs(globalbest_faval)E0,break, end k=k+1; %重复以上步骤，直到满足设定的最小误差或者达到最大迭代次数 end Value1=1/globalbest_faval-1; Value1=num2str(Value1); %将value1的数值转成字符串（number to string） disp(strcat(the maximum value,=,Value1); Value2=globalbest_x;,Value2=num2str(Value2); disp(strcat(the corresponding coordinate,=,Value2); x=-5:0.01:5; y=2.1*(1-x+2*x.2).*exp(-x.2/2); plot(x,y,m-,linewidth,3); hold on; plot(globalbest_x,1/globalbest_faval-1,kp,linewidth,4); legend(目标函数,搜索到的最大值); xlabel(x); ylabel(y); grid on; toc;,PSO算法应用案例,基于PSO算法和BP算法训练神经网络 BP神经网络容易陷入局部极小值，PSO算法在无约束非线性函数优化方面性能优越，通常可以直接找到全局最优解。本案例实现利用PSO算法和BP算法共同训练神经网络，先将网络进行PSO算法训练，然后BP算法进行小范围精细搜索，PSO算法训练神经网络的本质是将输出误差函数即能量函数看成目标函数，PSO对能量函数进行全局寻找最小值。,神经网络,BP算法,在网络训练中调整网络权值的训练算法，是反向传播算法（即BP学习算法）。网络权值其意义是两个神经元之间的连接强度。 当一对学习样本提供给输入神经元之后，神经元的激活值（该层神经元输出值）从输入层经过各隐含层向输出层传播，在输出层的各神经元获得网络输入响应，然后按照减少网络输出与实际输出样本之间误差的方向，从输出层反向经过各隐含层回到输入层，从而逐步修正各连接权值，这种算法称为误差反向传播算法，即BP算法。,BP算法的核心是数学中的“负梯度下降”理论，即BP网络的误差调整方向总之沿着误差下降最快的方向进行。常规三层BP网络权值和阈值调整公式如下： + + E为网络输出与实际输出样本之间的误差平方和， 为网络的学习速率即权值调整幅度 t时刻输入层第i个神经元与隐含层第j个神经元的连 接权值； B为神经元的阈值；,PSO-BP神经网络的设计指导原则 （1）PSO算法惯性因子采用线性递减，本例设定惯性因子初始值取0.9，PSO算法结束时等于0.4. （2）网络学习输出样本添加噪声，促进其泛化能力。 （3）PSO算法对网络训练完毕以后，BP算法接力训练神经网络。 （4）BP算法中学习速率尽量小，以便执行精细的搜索行为，BP算法的训练次数尽量不超过1000次。 （5）网络学习时可能会发生震荡，程序中设计了类似相机的功能即时拍下震荡图片，并在图片上显显示何时震荡。 （6）在保证网络得到足够的学习前提下，神经元的数目越少越好。 （7）为了使网络训练过程可视化，采用误差动画跟踪显示技术。 （8）PSO算法训练完成以后交由BP算法接着训练，在交接瞬间网络输出误差函数会出现略微上升。因为PSO算法是根据个体学习和社会信息共享原则的引导下完成寻优过程的，BP算法是按照确定的负梯度下降轨道完成寻优过程，两个过程的行走路线不重合。,网络交由BP算法训练,PSO-BP神经网络的实现,本案例网络结构为三层，基于PSO算法和BP算法先后训练神经网络的权值和阈值，然后逼近函数 f(x)=1.1（1-x+2）exp()。,求解程序,function NewW1,NewB1,NewW2,NewB2=PSOTrain(SamIn,SamOut,HiddenUnitNum); Maxgeneration=700; E0=0.0001; Xmin=-10; Xmax=10; Vmin=-5; Vmax=5; M=100; c1=2.7; c2=1.3; w=0.9; R,SamNum=size(SamIn); S2,SamNum=size(SamOut); generation=1; Done=0; Pb1=zeros(HiddenUnitNum,R+S2+1,M); Pb2=zeros(S2,M); Pg1=zeros(HiddenUnitNum,R+S2+1); Pg2=zeros(S2,1); E=zeros(size(SamOut); rand(state,sum(100*clock); startP1=rand(HiddenUnitNum,R+S2+1,M)-0.5;,startP2=rand(S2,M)-0.5; startV1=rand(HiddenUnitNum,R+S2+1,M)-0.5; startV2=rand(S2,M)-0.5; endP1=zeros(HiddenUnitNum,R+S2+1,M); endP2=zeros(S2,M); endV1=zeros(HiddenUnitNum,R+S2+1,M); endV2=zeros(S2,M); startE=zeros(1,M); endE=zeros(1,M); for i=1:M W1=startP1(1:HiddenUnitNum,1:R,i); W2=startP1(1:HiddenUnitNum,R+1:R+S2,i); B1=startP1(1:HiddenUnitNum,R+S2+1,i); B2=startP2(1:S2,i); for q=1:SamNum TempOut=logsig(W1*SamIn(:,q)+B1); NetworkOut(1,q)=W2*TempOut+B2; end E=NetworkOut-SamOut; startE(1,i)=sumsqr(E)/(SamNum*S2); Pb1(:,:,i)=startP1(:,:,i); Pb2(:,i)=startP2(:,i); end,val,position=min(startE(1,:); Pg1=startP1(:,:,position); Pg2=startP2(:,position); Pgvalue=val; Pgvalue_last=Pgvalue; while(Done) for num=1:M endV1(:,:,num)=w*startV1(:,:,num)+c1*rand*(Pb1(:,:,num)-startP1(:,:,num)+c2*rand*(Pg1-startP1(:,:,num); endV2(:,num)=w*startV2(:,num)+c1*rand*(Pb2(:,num)-startP2(:,num)+c2*rand*(Pg2-startP2(:,num); for i=1:HiddenUnitNum for j=1:(R+S2+1) endV1(i,j,num)=endV1(i,j,num); if endV1(i,j,num)Vmax endV1(i,j,num)=Vmax; elseif endV1(i,j,num)Vmin endV1(i,j,num)=Vmin; end end end,for s2=1:S2 endV2(s2,num)=endV2(s2,num); if endV2(s2,num)Vmax endV2(s2,num)=Vmax; elseif endV2(s2,num)Xmax endP1(i,j,num)=Xmax; elseif endP1(i,j,num)Xmin endP1(i,j,num)=Xmin; end end end,for s2=1:S2 if endP2(s2,num)Xmax endP2(s2,num)=Xmax; elseif endP2(s2,num)Xmin endP2(s2,num)=Xmin; end end W1=endP1(1:HiddenUnitNum,1:R,num); W2=endP1(1:HiddenUnitNum,R+1:R+S2,num); B1=endP1(1:HiddenUnitNum,R+S2+1,num); B2=endP2(1:S2,num); for q=1:SamNum TempOut=logsig(W1*SamIn(:,q)+B1); NetworkOut(1,q)=W2*TempOut+B2; end E=NetworkOut-SamOut; SSE=sumsqr(E) %便于在命令窗口观察网络误差的变化情况 endE(1,num)=sumsqr(E)/(SamNum*S2); if endE(1,num)startE(1,num),Pb1(:,:,num)=endP1(:,:,num); Pb2(:,num)=endP2(:,num); startE(1,num)=endE(1,num); end end if (generation=Maxgeneration) Done=1; end if Pgv