粒子群_Matlab程序代码.doc

粒子群优化算法(panicle swarm optimization，PSO)是kennedy和Eberhart在研究鸟类和鱼类的群体行为基础上于1995年提出的一种群智能算法，其思想米源予人工生命和演化计算理论，模仿鸟群飞行觅食行为，通过鸟集体协作使群体达到最优。1.粒子群算法的原理PSO中，每个优化问题的解看作搜索空间中的一只鸟(即粒子)，所有的粒子都有一个被优化的函数决定的适应值，并且有一个速度决定它们飞翔的方向和速率，粒子们追随当前的最优粒子在解空间中搜索。算法首先初始化一群随机粒子，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个“极值”即个体极值和全局极值来更新自己的速度与位置。在D维目标搜索空间中，由种群数为m的粒子组成粒子群，其中第f个粒子在第d维的位置为Xid，其飞行速度为Vid，该粒子当前搜索到的最优位置为Pid(goodvalue)和整个粒子群当前的最优位置Pgd(bestvalue)。每维的速度与位置更新公式如下W为惯性权重，C1和C2为学习因子，rand（）0,1范围内变化的随机数。2.参数介绍与设置（1）ww是保持粒子运动惯性的参数，能使种群扩展搜索空间，获得较好的求解效果。较大的w有利于群体在更大的范围内进行搜索。而较小的w能够保证群体收敛到最优位置，所以w的选择及在迭代中的变化对搜索能力和跳出局优能力具有重要影响。一般将w设定为0.8左右。（1）加速因子c1和c2c1和c2用于调整粒子自身经验和社会经验在其运动中的作用，表示将每个粒子拉向pbest和gbest位置的随机加速项的权重，低的值允许粒子在被拉回前可以在目标区域外徘徊, 而高的值则导致粒子突然冲向或越过目标区域。如果c1=0, 则粒子没有认知能力。在粒子的相互作用下,能到达新的搜索空间, 但容易陷入局部极值点; 如果c2=0, 粒子间没有社会信息共享, 其算法变成一个多起点的随机搜索; 如果c1=c2=0, 粒子将一直以当前的速度飞行, 直到到达边界。恰当地选择c1与c2能较好的得到最优解，一般都设定为2。（2）最大速度vmaxvmax决定当前位置与最好位置之间的区域的分辨率（或精度）。如果vmax太高, 微粒可能会飞跃最优解; 如果vmax太小, 则微粒易陷入局部最优。引入惯性权重w可消除对vmax的需要, 因为两者的作用都是维护全局和局部搜索能力的平衡。当vmax增加时, 可通过减小w来平衡搜索, 而w的减小可使得所需的迭代次数变小。但常常将设定一个最大的vmax。较好地选择和调整参数能够加大搜索能力，避免PSO的早熟现象。3.算法步骤及设计PSO的算法框架如下所述，图1给出了PS0的算法流程。(1)初始化所有的个体(粒子)，初始化它们的速度和位置，并且将个体的历史最优goodvalue设为当前位置，而群体中最优的个体历史位置作为当前的全局最优bestvalue。(2)在当代的进化中，计算各个粒子的适应度函数值。(3)如果该粒子当前的适应度函数值比其历史最优值要好，那么历史最优将会被当前位置所替代。（4）如果该粒子的历史最优比全局最优要好，那么全局最优将会被该粒子的历史最优所替代。(5)对每个粒子按照公式(1)和公式(2)对速度和位置进行更新。(6)进化代数增加1，如果还没有到达结束条件，转到(2)步，否则输出全局最优bestvalue并结束。开始初始化粒子群计算每个粒子新位置的适应值根据粒子适应值更新个体最优值和全局最优值根据公式更新粒子速度和位置N是否满足终止条件Y输出最优解结束图一 PSO算法流程4.程序设计与仿真（1）程序设计说明采用MATLAB软件编程，程序粒子种群个数设定为20个，搜索范围为全局版本，终止条件是粒子都到达最优点。（2）函数，其图形如图2所示，其在x=0.5，y=3.处有最小值f=-7.25.运行程序大约迭代220步，就可以让所有的粒子到达最优点。图3中坐标图从左到右从上到下分别为初始化种群的位置、k=50的粒子位置、k=100的粒子位置及迭代结束粒子的位置.从图中可见能够搜索到最优点。图2 f1三维图形 图3 仿真粒子位置图（2）Rastrigin函数，当n=2时，为二维图形如图4所示，在x1=0，x2=0.处有最大值f=0；在程序中取相反数求其最小值，迭代k=670步左右将可以将所有的粒子推向一点最优点。取前十个粒子的位置如图5，可知可以快速地得到满意的最优解，且能够跳出局优。图4 f2函数n=2三维图形图5 粒子位置坐标（3）rosenbrock函数：当n=2时，函数为二维，当x1=x2=1时有最小值f=0，其图形如图6，可知在最优点附近，函数值下降很慢，利用将所有粒子推向一点终止条件，大约迭代k=1200次左右结束。可以从坐标图中（图7）看到粒子的在迭代过程中位置。图6 f3函数n=2三维图形图7 粒子坐标位置当n=5时，即是5维的rosenbrock函数，当x1=x2=x3=x4=x5=1时有最小值f=0.采用粒子推向一点作为终止条件，运行程序迭代10万次（人为设定）后，全局最优值已到了理想值（f=0）但此时粒子的位置还是相离较远，且运行时间较长（大约一分钟），考虑程序目的就是找到最优值，而不是让所有粒子到达一点，因此采用迭代次数作为程序终止条件。依次增大设定的迭代次数观察运行结果，如下表1表1 取不同迭代次数最优点迭代次数k全局最优点位置x1 x2 x3 x4 x5最优解10000.0794 0.2819 0.5254 0.7267 0.85300.841530000.3718 0.6126 0.7832 0.8851 0.94070.215150000.5751 0.7583 0.8710 0.9328 0.96450.0816100000.7064 0.8413 0.9165 0.9572 0.97840.0348500000.9937 0.9969 0.9984 0.9992 0.99960.00001000000.9999 1.0000 1.0000 1.0000 1.00000.0000由表中可知，随着维数的增加，虽然迭代过程也是收敛的，但寻找最优点很难，前已述及，PSO的搜索能力与参数设置有关，且全局版本相对局部版本有较差的跳出局优的能力。为此若要增加搜索能力、提高迭代次数需要考虑改进参数和更换局部版本搜索方式。（5）程序代码（1）主函数：clearclchold offpsize=20;%粒子个数的设置pd=5;%粒子的维数lz=zeros(psize,pd);for i=1:psize %粒子群随机生成。psize行pd列，pd维。 lz(i,:)=rand(1,pd).*20-10;%【-1010之间】endlv=lz;lzfigure=lz;%画初始图用goodvalue=lz;%每个粒子自己历史最好值初始化,psize行pd列vmax=20;%速度上限c1=2;c2=2;%学习因子w=0.729;%随机因子和惯性因子。ri=0;bestvalue=zeros(1,pd);%全局最好值初始化，1行pd列for j=1:pd bestvalue(1,j)=goodvalue(1,j);endfnew=zeros(1,psize);for j=1:psize fnew(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) fold(i)=fnew(i);%fold记录每个粒子的最好历史值 for j=1:pd goodvalue(i,j)=lz(i,j); end end end%*.2.* for i=1:psize%每个粒子历史最好值相比较，更新全局最好值(及各维) f1=fold(i); if f1f0 f0=f1; for j=1:pd bestvalue(1,j)=goodvalue(i,j); end end end %*粒子趋一点终止条件*% % flagstop0=max(abs(fold);%比较当次的所有粒子的适应值， % flagstop1=min(abs(fold);%若它们之间差别较小，则可停止。 % if (flagstop0-flagstop1)1000 %迭代次数过大，强制终止，以防陷入死循环 flagstop=1; end %* %*速度和位置更新* for j=1:pd for i=1:psize lv(i,j)=w.*lv(i,j)+(c1*rand(1).*(goodvalue(i,j)-lz(i,j). +(c2*rand(1).*(bestvalue(1,j)-lz(i,j);%更新速度 if lv(i,j)vmax%速度最大值约束 lv(i,j)=vmax; end lz(i,j)=lz(i,j)+lv(i,j);%更新粒子位置 end end %*% % ri=ri+1; % rivalue(ri)=f0; % if ri=5000 % ri=0; % flagstop0=max(abs(rivalue); % flagstop1=min(abs(rivalue); % if (flagstop0-flagstop1)0.001 % flagstop=1; % end% end %*画图*% if k=0 %画出初始粒子的位置 subplot(2,2,1); t=1:psize; plot(lzfigure(t,1),lzfigure(t,2),*); axis(-10,10,-10,10); end for j=2:3%分别画出迭代300次和600次粒子的位置 if k=300*j subplot(2,2,j); t=1:psize; % gtext(k=600); plot(lz(t,1),lz(t,2),*); axis(0,10,0,10); end end if flagstop=1 %画出终止时粒子的位置 subplot(2,2,4); t=1:psize; plot(lz(t,1),lz(t,2),*); axis(0,10,0,10); end%*% k=k+1;endfprintf(总迭代次数：k=%ld。n,k-1);fprintf(每个粒子的历史最好值：n);for i=1:psize fprintf(粒子%d： %6.4f,%6.4fn,i,goodvalue(i,1),goodvalue(i,2);endfprintf(全局最好值：%6.4f,%6.4f.n,bestvalue(1,1),bestvalue(1,2);fprintf(最优解(最大值)：%6.4f.n,fpso(bestvalue(1,:);for j=1:20 fprintf(粒子%d位置：,j) ; fprintf(x=%6.4f，y=%6.4f,lz(j,1),lz(j,2); fprintf(n) ;end（2）调用函数：function fvalue=fpso(x)%fvalue=(x(1)+2)*(x(1)-3)+(x(2)-4)*(x(2)-2);%f1pd=5;%粒子的维数%*f3*fvalue=0; for r=1:(pd-1) f0=100*(x(r)-x(r+1)2)2+(x(r+1)-1)2; fvalue=fvalue+f0; end %* %fvalue=(x(1)-2)*(x(1)-2.5)+(x(2)-4