《粒子群Matlab程序代码》由会员分享，可在线阅读，更多相关《粒子群Matlab程序代码（9页珍藏版）》请在金锄头文库上搜索。
1、粒子群优化算法（panicle swarm optimization , PSO是 kennedy 和 Eberhart 在 研究鸟类和鱼类的群体行为基础上于1995年提出的一种群智能算法，其思想米源予人工生命和演化计算理论，模仿鸟群飞行觅食行为，通过鸟集体协作使群体 达到最优。1粒子群算法的原理PSO中，每个优化问题的解看作搜索空间中的一只鸟（即粒子），所有的粒子都有一个被优化的函数决定的适应值，并且有一个速度决定它们飞翔的方向和速 率，粒子们追随当前的最优粒子在解空间中搜索。算法首先初始化一群随机粒子，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个“极值”即个体 极值和全局极值来
2、更新自己的速度与位置。在D维目标搜索空间中，由种群数为m的粒子组成粒子群，其中第f个粒子在第d维的位置为Xid，其飞行速度为 Vid，该粒子当前搜索到的最优位置为 Pid（goodvalue和整个粒子群当前的最优位 置 Pgd（bestvalue。每维的速度与位置更新公式如下vid 1vid c1 rand （1） （pid xid） c2 rand （1） （gidXid 1 Xid Vid 1（2）W为惯性权重，C1和C2为学习因子，rand（） 0,1范围内变化的随机数。 2参数介绍与设置（1）ww是保持粒子运动惯性的参数，能使种群扩展搜索空间，获得较好的求解效 果。较大的w有利于群体在
3、更大的范围内进行搜索。而较小的w能够保证群体收敛到最优位置，所以 w的选择及在迭代中的变化对搜索能力和跳出局优能力 具有重要影响。一般将 w设定为0.8左右。（1）加速因子c1和c2c1和c2用于调整粒子自身经验和社会经验在其运动中的作用，表示将每个 粒子拉向pbest和gbest位置的随机加速项的权重，低的值允许粒子在被拉回前 可以在目标区域外徘徊，而高的值则导致粒子突然冲向或越过目标区域。如果 c1=0,则粒子没有认知能力。在粒子的相互作用下，能到达新的搜索空间，但容易 陷入局部极值点；如果c2=0,粒子间没有社会信息共享，其算法变成一个多起点 的随机搜索；如果c仁c2=0,粒子将一直以当
4、前的速度飞行，直到到达边界。恰当 地选择c1与c2能较好的得到最优解，一般都设定为 2。（2）最大速度vmaxvmax决定当前位置与最好位置之间的区域的分辨率（或精度） 。如果vmax 太高，微粒可能会飞跃最优解；如果vmax太小，则微粒易陷入局部最优。引入惯 性权重w可消除对vmax的需要，因为两者的作用都是维护全局和局部搜索能力 的平衡。当vmax增加时，可通过减小w来平衡搜索，而w的减小可使得所需的 迭代次数变小。但常常将设定一个最大的 vmax。较好地选择和调整参数能够加大搜索能力，避免PSO的早熟现象。3算法步骤及设计PSO的算法框架如下所述，图1给出了 PS0的算法流程。（1）初始
5、化所有的个体（粒子），初始化它们的速度和位置，并且将个体的历史 最优goodvalue设为当前位置，而群体中最优的个体历史位置作为当前的全局最优 bestvalue。(2) 在当代的进化中，计算各个粒子的适应度函数值。(3) 如果该粒子当前的适应度函数值比其历史最优值要好，那么历史最优将会 被当前位置所替代。(4) 如果该粒子的历史最优比全局最优要好，那么全局最优将会被该粒子 的历史最优所替代。(5) 对每个粒子按照公式(1)和公式(2)对速度和位置进行更新。(6) 进化代数增加1,如果还没有到达结束条件，转到(2)步，否则输出全局最 优bestvalue并结束。开始图一 PSO算法流程4程序
6、设计与仿真(1)程序设计说明采用MATLAB软件编程，程序粒子种群个数设定为20个，搜索范围为全局版本, 终止条件是粒子都到达最优点。(2)函数fi (x 2) (x 3) (y 4) (y 2)，其图形如图2所示，其在x=0.5, y=3处有最小值f=-7.25.运行程序大约迭代220步，就可以让所有的粒子到达最优 点。图3中坐标图从左到右从上到下分别为初始化种群的位置、k=50的粒子位置、k=100的粒子位置及迭代结束粒子的位置从图中可见能够搜索到最优点。50C -H -ip图2 fi三维图形图3仿真粒子位置图图4 f2函数n=2三维图形n(2)Rastrigin 函数 f2(刃2 10
7、cos(2 刈 10)，当 n=2 时，为二维图i 1形如图4所示，在x仁0，x2=0处有最大值f=0；在程序中取相反数求其最小值， 迭代k=670步左右将可以将所有的粒子推向一点一一最优点。取前十个粒子的位 置如图5,可知可以快速地得到满意的最优解，且能够跳出局优粒子x=O. 0083, y=0. 0038 粒于2 位置：x=-0. 0000? v-0. 0000 粒子3位置；x=OLOOOO, rLOOOO 粒子4位置；=-0. 0000, =0. 0000 粒子5位置* X-0. 00001 y=0. 0000 粒子6位置：x=0.0000： y=0. 0000 粒子了 位置：=0. 0
8、000* y=0. 0000 粒子R位皆：x=0. 0000 y=-0. 0000 粒于9位置:x=0. 0000： y-0. 0000粒子 10位置:x=-0. 0000, y=0. 0000图5粒子位置坐标n(3) rosenbrock 函数：f3(100(xi 1 x(xi 1)2),x 10,10.当 n=2 时，函数i 1为二维，当x仁x2=1时有最小值f=0,其图形如图6，可知在最优点附近，函数 值下降很慢，利用将所有粒子推向一点终止条件， 大约迭代k=1200次左右结束可以从坐标图中（图7）看到粒子的在迭代过程中位置图6 f3函数n=2三维图形10 0 & 10图7粒子坐标位置当
9、n=5时，即是5维的rosenbrock函数，当x1= x2=x3=x4=x5=1时有最小值f=0. 采用粒子推向一点作为终止条件，运行程序迭代 10万次（人为设定）后，全局 最优值已到了理想值（f=0）但此时粒子的位置还是相离较远，且运行时间较长（大约一分钟），考虑程序目的就是找到最优值，而不是让所有粒子到达一点， 因此采用迭代次数作为程序终止条件。依次增大设定的迭代次数观察运行结果， 如下表1表1取不同迭代次数最优点迭代次数k全局最优点位置x1 x2 x3 x4 x5最优解10000.07940.28190.52540.72670.85300.841530000.37180.61260.7
10、8320.88510.94070.215150000.57510.75830.87100.93280.96450.0816100000.70640.84130.91650.95720.97840.0348500000.99370.99690.99840.99920.99960.00001000000.99991.00001.00001.00001.00000.0000由表中可知，随着维数的增加，虽然迭代过程也是收敛的，但寻找最优点很难， 前已述及，PSO的搜索能力与参数设置有关，且全局版本相对局部版本有较差的 跳出局优的能力。为此若要增加搜索能力、提高迭代次数需要考虑改进参数和更换局部版本搜索
11、方式( 5)程序代码(1)主函数：clearclchold offpsize=20;%4子个数的设置pd=5;%粒子的维数lz=zeros(psize,pd);for i=1:psize%粒子群随机生成。psize行pd列，pd维。lz(i,:)=rand(1 ,pd).*20-1 0;%【-1 01 0之间】endlv=lz;lzfigure=lz;%画初始图用goodvalue=lz;%每个粒子自己历史最好值初始化，psize行pd列 vmax=20;%速度上限c1=2;c2=2;%学习因子w=0.729;%随机因子和惯性因子。ri=0;bestvalue=zeros(1,pd);%b局最
12、好值初始化，1 行 pd 列for j=1:pdbestvalue(1,j)=goodvalue(1,j);endfnew=zeros(1,psize);for j=1:psizefnew(j)=fpso(lz(1,:);endfold=fnew;flagstop=0;%终止标志k=0;%迭代次数记录f0=fpso(bestvalue);%适应值初始化(默认)while flagstop=0%*.1.*for i=1:psize%适应值比较，更新各自历史最好值(位置) fnew(i)=fpso(lz(i,:);%记录每次每个粒子的适应值，便于以后设置终止条 件if fnew(i)fold(i)
13、fold(i)=fnew(i);%fold 记录每个粒子的最好历史值for j=1:pdgoodvalue(i,j)=lz(i,j);endendend*.2.*for i=1:psize%每个粒子历史最好值相比较，更新全局最好值(及各维)f1=fold(i);if f1f0f0=f1;for j=1:pdbestvalue(1,j)=goodvalue(i,j);endendend%*粒子趋一点终止条件 *%flagstopO=max(abs(fold);%t匕较当次的所有粒子的适应值，%flagstop仁mi n(abs(fold);%若它们之间差别较小，则可停止。% if (flagst
14、opO-flagstop1)1OOO %迭代次数过大，强制终止，以防陷入死循环 flagstop=1;end速度和位置更新 *%*%*for j=1:pdfor i=1:psizelv(i,j)=w.*lv(i,j)+(c1*rand(1).*(goodvalue(i,j)-lz(i,j). +(c2*ra nd(1).*(bestvalue(1,j)-lz(i,j);%更新速度 if lv(i,j)vmax%速度最大值约束lv(i,j)=vmax;endlz(i,j)=lz(i,j)+lv(i,j);%更 新粒子位置endend*%*%ri=ri+1;%rivalue(ri)=fO;%if ri
