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

粒子群算法原理及其在函数优化中的应用

1粒子群优化(PSO)算法基本原理

1.1标准粒子群算法

假设在一个D 维的目标搜索空间中，有 m 个代表问题潜在解的粒子组成一 个种群x [X i ,X 2,...,X m ],第i 个粒子的信息可用D 维向量表示为 X i [X ii , X i2,..., X iD ]T ，其速度为V i [V ii ,V i2,...,V iD ]T 。算法首先初始化m 个随机粒

子，然后通过迭代找到最优解。每一次迭代中，粒子通过跟踪2个极值进行信息 交流，一个是第i 个粒子本身找到的最优解，称之为个体极值，即 P i [P il , P i2,...,厢]丁 ；另一个是所有粒子目前找到的最优解，称之为群体极值， 即P g [P gi ,P g2,..., P gD 「。粒子在更新上述2个极值后，根据式(1)和式(2)更新自

己的速度和位置。

t 1 t t t t t\

V i

WV i C 1「1（P i X i ） C 2「2（P g X i ） 式中，t 代表当前迭代次数，「1,「2是在[0,1]之间服从均匀分布的随机数，C 1,C 2 称为学习因子，分别调节粒子向个体极值和群体极值方向飞行的步长， w 为惯性 权重，一般在0.1~0.9之间取值。在标准的PSO 算法中，惯性权重w 被设为常数, 通常取w 0.5。在实际应用中，x 需保证在一定的范围内，即x 的每一维的变化 范围均为[X min ,X max ],这在函数优化问题中相当丁自变量的定义域

1.2算法实现步骤

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

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

步骤3:通过循环迭代更新x 、p i 和p g :

① 确定惯性权重w 的取值(当w 不是常数时)。

② 根据式(1)更新粒子的速度V ：1,若速度中的某一维超过了 V max ，则取为 V

max - ③ 根据式(2)更新自变量x ,若x 的取值超过其定义域，则在其定义域内重新 初t 1 X i t t 1

X i V

i

始化。

④代入fitness（x）求得适应度值，通过比较更新个体极值P i和全局极值p g <

步骤4:判断是否满足终止条件（通常设为达到最大迭代次数或达到估计精度

要求），若不满足，则转入步骤（3）,若满足，则输出估计结果，算法结束。

2程序实现

2.1各种测试函数（适应度函数）

测试函数是用来测试算法性能的一些通用函数，下面先给出一些测试函数的三维图（自变量为两维，加上函数值共三维）如图1-图17所示。

图1测试函数1

4

3

2

1

-1

10

图2测试函数2

2500

2000

1500

1000

500

10

图4测试函数4

1

0.5

-0.5

5

-5 -5

f6

1

0.8

0.6

0.4

0.2

5

图8测试函数8

-100 -100 120 100 80 60 40 20 0 10 100 80 60 40 20 0 100

图10测试函数10

-10 -10 0 -2 -4 -6 -8 10 200 150 100 50 0 10

图12测试函数12

-10 -10 1 0.8 0.6 0.4 0.2 0 10 15 10 5 0 10

20

15

10

5

10

图14测试函数14

120

100

80

60

40

10

250

200

150

100

50

10

图16测试函数16

2.5

2

1.5

1

0.5

10

-10

-10

2.2程序实现

首先给出绘制测试函数的程序：

%%绘图测试函数 draw_fitness.m clear;clc;close all;

%%

[x,y]=meshgrid(-10:0.5:10);

z2 = x.A 2-cos(18*x)+y.A 2-cos(18*y);

figure(1); surf(x,y,z2); minz2=min(min(z2));title( 'z2 = xA2-cos(18*x)+yA2-cos(18*y)');

%%

z4 = 4*x.A2-2.1*x.A4+x.A6/3+x.*y-4*y.A2+y.A4;

figure(2); surf(x,y,z4); minz4=min(min(z4));title( 'z4 = 4*xA2-2.1*xA4+xA6/3+x*y-4*yA2+yA4'

);

%%

z5 = (y-5.1*x.A2/4/pi/pi+5*x/pi-6).A2+10*(1-1/8/pi)*cos(x)+10; figure(3); surf(x,y,z5);

minz5=min(min(z5));title( 'z5 = (y-5.1*xA2/4/\pi/\pi+5*x/\pi-6)A2+10*(1-1/8/\pi)*cos(x)+10' );

%% [x,y]=meshgrid(-5:0.5:5);

z7 = x.*exp(-x.A2-y.A2);

a=(cos(x)).A 4+(cos(y)).A 4;

b=-2*((cos(x)).人2.*(cos(y)).人2);

figure(4); surf(x,y,z7); minz7=min(min(z7));title( %%

[x,y]=meshgrid(-5:0.25:5);

z8 = 3*(1-x).A2.*exp(-(x.A2) - &+1)人2) ...

-10*(x/5 - x.A3 - y.A5).*exp(-x.A2-y.A2)

-1/3*exp(-(x+1).A2 - y.A2);

figure(5); surf(x,y,z8); minz8=min(min(z8)); title( %%

r=sqrt(x.A2+y.A2)+eps; z9=sin(r)./r;

figure(6); surf(x,y,z9); minz9=min(min(z9));title(

%% [x,y]=meshgrid(-5:0.25:5);

num=sin(sqrt(x.A2+y.A2)).A2 - 0.5;

den=(1.0+0.01*(x.A2+y.A2)).A2;

z10=0.5 +num./den;

figure(7); surf(x,y,z10); minz10=min(min(z10));title(

%% [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);

z7 = x*exp(-xA2-yA2)');

'z8 ');

r=sqrt(xA2+yA2)+eps; z9=sin(r)/r;');

'f6');