2022年2022年粒子群算法原理及在函数优化中的应用 .pdf

作者：胡斌（ Author ：xiao5） 2012-5-5 1 粒子群算法原理及其在函数优化中的应用1 粒子群优化（ PSO）算法基本原理1.1 标准粒子群算法假设在一个D维的目标搜索空间中，有m个代表问题潜在解的粒子组成一个 种 群12,.,mxx xx， 第i个 粒 子 的 信 息 可 用D维 向 量 表 示 为12,.,TiiiiDxxxx，其速度为12,.,TiiiiDvvvv。算法首先初始化m个随机粒子，然后通过迭代找到最优解。 每一次迭代中， 粒子通过跟踪 2 个极值进行信息交 流 ， 一 个 是 第i个 粒 子 本 身 找 到 的 最 优 解 ， 称 之 为 个 体 极 值 ， 即12,.,TiiiiDpppp；另一个是所有粒子目前找到的最优解，称之为群体极值，即12,.,TggggDpppp。粒子在更新上述2 个极值后，根据式 (1)和式(2)更新自己的速度和位置。11 12 2()()ttttttiiiigiwc rc rvvpxpx(1) 11tttiiixxv(2) 式中，t代表当前迭代次数，12,r r 是在0,1之间服从均匀分布的随机数，12,c c称为学习因子， 分别调节粒子向个体极值和群体极值方向飞行的步长，w为惯性权重，一般在0.1 0.9之间取值。在标准的 PSO算法中，惯性权重 w被设为常数，通常取0.5w。在实际应用中，x需保证在一定的范围内，即x的每一维的变化范围均为minmax,XX，这在函数优化问题中相当于自变量的定义域。1.2 算法实现步骤步骤 1：表示出 PSO算法中的适应度函数( )fitness x； （编程时最好以函数的形式保存，便于多次调用。 ）步骤 2：初始化 PSO 算法中各个参数 (如粒子个数，惯性权重，学习因子，最大迭代次数等 )，在自变量x定义域内随机初始化x，代入( )fitness x求得适应度值，通过比较确定起始个体极值ip 和全局极值gp 。步骤 3：通过循环迭代更新x、ip 和gp ：确定惯性权重 w的取值（当 w不是常数时）。根据式 (1)更新粒子的速度1kiv，若速度中的某一维超过了maxV，则取为maxV。根据式 (2)更新自变量x，若x的取值超过其定义域，则在其定义域内重新名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 2 初始化。代入( )fitness x求得适应度值，通过比较更新个体极值ip和全局极值gp。步骤 4： 判断是否满足终止条件 (通常设为达到最大迭代次数或达到估计精度要求)，若不满足，则转入步骤(3)，若满足，则输出估计结果，算法结束。2 程序实现2.1 各种测试函数（适应度函数）测试函数是用来测试算法性能的一些通用函数，下面先给出一些测试函数的三维图（自变量为两维，加上函数值共三维）如图1-图 17 所示。-10-50510-10-50510-50050100150200250z2 = x2-cos(18*x)+y2-cos(18*y)图 1 测试函数1 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 3 -10-50510-10-50510-101234x 105z4 = 4*x2-2.1*x4+x6/3+x*y-4*y2+y4图 2 测试函数2 -10-50510-10-5051005001000150020002500z5 = (y-5.1*x2/4/ / +5*x/-6)2+10*(1-1/8/)*cos(x)+10图 3 测试函数3 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 4 -505-505-0.4-0.200.20.4z7 = x*exp(-x2-y2)图 4 测试函数4 -505-505-10-50510z8 图 5 测试函数5 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 5 -505-505-0.500.51 r=sqrt(x2+y2)+eps; z9=sin(r)/r;图 6 测试函数6 -505-50500.20.40.60.81f6图 7 测试函数7 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 6 -10-50510-10-50510020406080100120图 8 测试函数8 -100-50050100-100-50050100020406080100图 9 测试函数9 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 7 02468100510-8-6-4-20图 10 测试函数10 -10-50510-10-50510050100150200NDparabola图 11 测试函数 11 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 8 -10-50510-10-5051000.20.40.60.81f6图 12 测试函数12 -10-50510-10-50510051015x 105Rosenbrock图 13 测试函数13 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 9 -10-50510-10-5051005101520ackley图 14 测试函数14 -10-50510-10-50510406080100120tripod图 15 测试函数15 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 10 -10-50510-10-50510050100150200250Rastrigin图 16 测试函数16 -10-50510-10-5051000.511.522.5Griewank图 17 测试函数17 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 11 2.2 程序实现首先给出绘制测试函数的程序：% 绘图测试函数draw_fitness.mclear;clc;close all;% x,y=meshgrid(-10:0.5:10);z2 = x.2-cos(18*x)+y.2-cos(18*y);figure(1); surf(x,y,z2); minz2=min(min(z2);title(z2 = x2-cos(18*x)+y2-cos(18*y);% z4 = 4*x.2-2.1*x.4+x.6/3+x.*y-4*y.2+y.4;figure(2); surf(x,y,z4); minz4=min(min(z4);title(z4 = 4*x2-2.1*x4+x6/3+x*y-4*y2+y4);%z5 = (y-5.1*x.2/4/pi/pi+5*x/pi-6).2+10*(1-1/8/pi)*cos(x)+10;figure(3); surf(x,y,z5); minz5=min(min(z5);title(z5 = (y-5.1*x2/4/pi/pi+5*x/pi-6)2+10*(1-1/8/pi)*cos(x)+10);%x,y=meshgrid(-5:0.5:5);z7 = x.*exp(-x.2-y.2);figure(4); surf(x,y,z7); minz7=min(min(z7);title(z7 = x*exp(-x2-y2);%x,y=meshgrid(-5:0.25:5);z8 = 3*(1-x).2.*exp(-(x.2) - (y+1).2) .- 10*(x/5 - x.3 - y.5).*exp(-x.2-y.2) .- 1/3*exp(-(x+1).2 - y.2); figure(5); surf(x,y,z8); minz8=min(min(z8); title(z8 );%r=sqrt(x.2+y.2)+eps; z9=sin(r)./r;figure(6); surf(x,y,z9); minz9=min(min(z9);title( r=sqrt(x2+y2)+eps; z9=sin(r)/r;);%x,y=meshgrid(-5:0.25:5);num=sin(sqrt(x.2+y.2).2 - 0.5;den=(1.0+0.01*(x.2+y.2).2;z10=0.5 +num./den;figure(7); surf(x,y,z10); minz10=min(min(z10);title(f6);% x,y=meshgrid(-10:0.5:10);z12=abs(x)+abs(y)+abs(x).*abs(y);figure(8); surf(x,y,z12); minz12=min(min(z12);%x,y=meshgrid(-100:5:100);z13=max(abs(x),abs(y);figure(9); surf(x,y,z13); minz13=min(min(z13);%x,y=meshgrid(0:0.5:10);名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 12 a=(cos(x).4+(cos(y).4;b=-2*(cos(x).2.*(cos(y).2);c=sqrt(x.2+2*y.2);z14=-abs(a-b)./c);figure(10); surf(x,y,z14); minz14=min(min(z14);%x,y=meshgrid(-10:0.5:10);NDparabola=x.2+y.2;figure(11); surf(x,y,NDparabola); title( NDparabola );%x,y=meshgrid(-10:0.5:10);num=sin(sqrt(x.2+y.2).2 - 0.5;den=(1.0+0.01*(x.2+y.2).2;f6=0.5 +num./den;figure(12); surf(x,y,f6); title(f6);%Rosenbrock=100*(x.2 - y).2 + (1-x).2;figure(13); surf(x,y,Rosenbrock); title(Rosenbrock );%x,y=meshgrid(-10:0.5:10);ackley=20 + exp(1) -20*exp(-0.2*sqrt(1/2).*(x.2+y.2) .-exp(1/2).*(cos(2*pi*x)+cos(2*pi*y);figure(14); surf(x,y,ackley);title(ackley);%x,y=meshgrid(-10:0.5:10);px1=(x) = 0);px2=(y) = 0);tripod= px2.*(1+px1) + abs(x + 50*px2.*(1-2*px1)+ abs(y + 50*(1-2.*px2) ;figure(15); surf(x,y,tripod); title(tripod );%x,y=meshgrid(-10:0.5:10);Rastrigin=(x.2-10*cos(2*pi*x)+10)+(y.2-10*cos(2*pi*y)+10);figure(16); surf(x,y,Rastrigin); title(Rastrigin );%x,y=meshgrid(-10:0.3:10);Griewank=1/4000*(x.2+y.2)-cos(x).*cos(y./sqrt(2)+1;figure(17); surf(x,y,Griewank); title(Griewank );以上只是绘制测试函数的 m文件， 其目的在于对测试函数有一个直观的认识。但以上的程序在 PSO算法的实现中是用不着的， 下面给出几个典型测试函数的代码。需要注意的是，一次只能调用其中一个测试函数。% 典型测试函数的调用函数fitness.m %function out=NDparabola(in) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 13 out = sum(in.2, 2); % % function out=f6(in)% x=in(:,1);% y=in(:,2);% num=sin(sqrt(x.2+y.2).2 - 0.5;% den=(1.0+0.01*(x.2+y.2).2;% out=0.5 +num./den;% % function out=DeJong_f2(in)% x= in(:,1);% y= in(:,2);% out = 100*(x.2 - y).2 + (1-x).2;% % function out=ackley(in)% % dimension is # of columns of input, x1, x2, ., xn% n=length(in(1,:);% x=in;% e=exp(1);% out = (20 + e .% -20*exp(-0.2*sqrt(1/n).*sum(x.2,2) .% -exp(1/n).*sum(cos(2*pi*x),2);% return% % function out=tripod(in)% x1=in(:,1);% x2=in(:,2);% px1=(x1) = 0);% px2=(x2) = 0);% out= ( px2.*(1+px1) .% + abs(x1 + 50*px2.*(1-2*px1).% + abs(x2 + 50*(1-2.*px2) );% function F=Rastrigin(in)% F=sum(in.2-10*cos(2*pi*in)+10);% function y=Griewank(in)% row,col=size(in);% y1=1/4000*sum(in.2);% y2=1;% for h=1:col% y2=y2*cos(in(h)/sqrt(h);% end% y=y1-y2+1; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 13 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 14 下面以测试函数 ackley为例，给出基本粒子群算法的程序实现代码：A、 主程序（可运行的程序）drawPSO.m：%基于PSO算法的函数寻优收敛图绘制drawPSO.m clear;clc;close all; %清除变量，清除命令窗口，关闭所以绘图窗口% %函数调用，注意函数 fitness前要加 gbest,M,xmin,fmin=PSO(fitness,5,2,1.5,0.5,50,2); figure(1); k=1:M; plot(k,gbest,:rp,LineWidth,2); %绘图title(ackley函数收敛曲线 ); %图形标题说明xlabel(迭代次数 ),ylabel(适应度函数最小值 ); %横纵坐标说明B、 PSO算法实现程序 PSO.m（核心程序）下面给出的是基本粒子群算法程序，没有做任何改进，但却很有效。%基本粒子群优化算法函数PSO.m function gbest,M,xmin,fmin=PSO(fitness,N,c1,c2,w,M,D) % 待优化的目标函数：fitness % 粒子数目： N % 学习因子 1：c1 % 学习因子 2：c2 % 惯性权重： w % 最大迭代次数：M % 问题的维数： D % 目标函数取最小值时的自变量值：xmin % 目标函数的最小值：fmin % gbest 是一存储了所有目标函数的最小值的向量，结合M 便于作出收敛图format long; %设置数据格式for i=1:N for j=1:D x(i,j)=rand; %随机初始化位置v(i,j)=rand; %随机初始化速度end end for i=1:N %利用循环初始化每个粒子的适应度fp(i)=fitness(x(i,:); xp(i,:)=x(i,:); %每个粒子对应的位置end xg=x(1,:); %给全局最优函数值的位置xg 赋初值for i=2:N %通过循环比较找出全局最优函数值的位置xg if fitness(x(i,:)fitness(xg) xg=x(i,:); end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 15 end for t=1:M for i=1:N v(i,:)=w*v(i,:)+c1*rand*(xp(i,:)-x(i,:)+c2*rand*(xg-x(i,:); %更新粒子的速度x(i,:)=x(i,:)+v(i,:); %更新粒子的位置if fitness(x(i,:)fp(i) %通过比较找出个体最优函数值及其位置fp(i)=fitness(x(i,:); xp(i,:)=x(i,:); end if fp(i)fitness(xg) %通过比较找出群体最优函数值及其位置xg=xp(i,:); end end gbest(t)=fitness(xg); %将每个全局最优函数值存储在gbest 向量中end xmin=xg; %输出最优位置及函数值fmin=fitness(xg); C、 适应度函数程序fitness.m （需保存为独立的m 文件）% 适应度函数函数fitness.m function out=ackley(in) n=length(in(1,:); x=in; e=exp(1); out = (20 + e-20*exp(-0.2*sqrt(1/n).*sum(x.2,2) . -exp(1/n).*sum(cos(2*pi*x),2); return 2.3 主程序运行结果运行主程序 drawPSO.m，结果如图 18所示。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 16 页 - - - - - - - - - 作者：胡斌（ Author ：xiao5） 2012-5-5 16 0510152025303540455000.20.40.60.811.21.41.61.82ackley 函 数 收 敛 曲 线迭 代 次 数适应度函数最小值图 18 PSO 函数寻优收敛曲线名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 16 页 - - - - - - - - -