《粒子群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