粒子群算法原理及在函数优化中的应用(附程序)

粒子群算法原理及其在函数优化中的应用1 粒子群优化（PSO ）算法基本原理1.1 标准粒子群算法假设在一个D 维的目标搜索空间中，有m 个代表问题潜在解的粒子组成一个种群12[,,...,]m =x x x x ，第i 个粒子的信息可用D 维向量表示为12[,,...,]T i i i iD x x x =x ，其速度为12[,,...,]T i i i iD v v v =v 。

算法首先初始化m 个随机粒子，然后通过迭代找到最优解。

每一次迭代中，粒子通过跟踪2个极值进行信息交流，一个是第i 个粒子本身找到的最优解，称之为个体极值，即12[,,...,]T i i i iD p p p =p ；另一个是所有粒子目前找到的最优解，称之为群体极值，即12[,,...,]T g g g gD p p p =p 。

粒子在更新上述2个极值后，根据式(1)和式(2)更新自己的速度和位置。

11122()()t t t t t t i i i i g i w c r c r +=+-+-v v p x p x (1)11t t t i i i ++=+x x v (2)式中，t 代表当前迭代次数，12,r r 是在[0,1]之间服从均匀分布的随机数，12,c c 称为学习因子，分别调节粒子向个体极值和群体极值方向飞行的步长，w 为惯性权重，一般在0.1~0.9之间取值。

在标准的PSO 算法中，惯性权重w 被设为常数，通常取0.5w =。

在实际应用中，x 需保证在一定的围，即x 的每一维的变化围均为min max [,]X X ，这在函数优化问题中相当于自变量的定义域。

1.2 算法实现步骤步骤1：表示出PSO 算法中的适应度函数()fitness x ；（编程时最好以函数的形式保存，便于多次调用。

）步骤2：初始化PSO 算法中各个参数(如粒子个数，惯性权重，学习因子，最大迭代次数等)，在自变量x 定义域随机初始化x ，代入()fitness x 求得适应度值，通过比较确定起始个体极值i p 和全局极值g p 。

步骤3：通过循环迭代更新x 、i p 和g p ： ①确定惯性权重w 的取值（当w 不是常数时）。

②根据式(1)更新粒子的速度1k i +v ，若速度中的某一维超过了max V ，则取为maxV 。

③根据式(2)更新自变量x ，若x 的取值超过其定义域，则在其定义域重新初始化。

④代入()fitness x 求得适应度值，通过比较更新个体极值i p 和全局极值g p 。

步骤4：判断是否满足终止条件(通常设为达到最大迭代次数或达到估计精度要求)，若不满足，则转入步骤(3)，若满足，则输出估计结果，算法结束。

2 程序实现2.1 各种测试函数（适应度函数）测试函数是用来测试算法性能的一些通用函数，下面先给出一些测试函数的三维图（自变量为两维，加上函数值共三维）如图1-图17所示。

z2 = x 2-cos(18*x)+y 2-cos(18*y)图1 测试函数1z4 = 4*x2-2.1*x4+x6/3+x*y-4*y2+y4图2 测试函数2z5 = (y-5.1*x2/4/π/π+5*x/π-6)2+10*(1-1/8/π)*cos(x)+10-10图3 测试函数3z7 = x*exp(-x2-y2)图4 测试函数4z8-5图5 测试函数5r=sqrt(x2+y2)+eps; z9=sin(r)/r;图6 测试函数6f6-5图7 测试函数7图8 测试函数8-100图9 测试函数9图10 测试函数10NDparabola-10图11 测试函数11f6图12 测试函数12Rosenbrock-10图13 测试函数13图14 测试函数14tripod图15 测试函数15图16 测试函数16Griewank-10图17 测试函数172.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 = x^2-cos(18*x)+y^2-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*x^2-2.1*x^4+x^6/3+x*y-4*y^2+y^4'); %%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*x^2/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(-x^2-y^2)');%%[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(x^2+y^2)+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);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)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;下面以测试函数ackley为例，给出基本粒子群算法的程序实现代码：A、主程序（可运行的程序）drawPSO.m：%基于PSO算法的函数寻优收敛图绘制drawPSO.mclear;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（核心程序）下面给出的是基本粒子群算法程序，没有做任何改进，但却很有效。