粒子群算法优化不同维数的连续函数以及离散函数的最小值问题.

引言 (2)

一、问题描述 (3)

1.1 函数优化问题 (3)

1.2 粒子群算法基本原理 (3)

二、算法设计 (5)

2.1算法流程框图 (5)

2.2 算法实现 (5)

2.3 算法的构成要素 (6)

2.4 算法的改进 (7)

三、算例设计 (8)

3.1 测试函数介绍 (8)

3.2 优化函数特点 (8)

四、仿真实验设计 (10)

4.1 实验参数设计 (10)

4.2 基本粒子群算法在测试函数中应用 (11)

五、仿真实验结果分析 (12)

5.1 实验结果汇总 (12)

5.2 实验结果分析 (13)

六、总结与展望 (14)

6.1总结 (14)

6.2展望 (14)

附录一 (15)

附录二 (17)

引言

本文主要利用粒子群算法解决连续函数以及离散函数的最小值问题，粒子群优化是一种新兴的基于群体智能的启发式全局搜索算法，粒子群优化算法通过粒子间的竞争和协作以实现在复杂搜索空间中寻找全局最优点。它具有易理解、易实现、全局搜索能力强等特点，倍受科学与工程领域的广泛关注，已经成为发展最快的智能优化算法之一。

惯性权重是PSO标准版本中非常重要的参数，可以用来控制算法的开发(exploitation)和探索(exploration)能力。惯性权重的大小决定了对粒子当前速度继承的多少。较大的惯性权重将使粒子具有较大的速度，从而有较强的探索能力；较小的惯性权重将使粒子具有较强的开发能力。关于惯性权重的选择一般有常数和时变两种。算法的执行效果很大程度上取决于惯性权重的选取。

本文介绍了粒子群优化算法的基本原理，分析了其特点，并将其应用于函数优化问题求解。此外，本文根据惯性权重对粒子群优化算法性能影响的研究，提出了三种不同的惯性权重。通过仿真实验，验证了各自的收敛性．同时也说明了惯性权重在粒子群优化算法中有很大的自由度。

一、问题描述

1.1 函数优化问题

目标优化问题可以描述为:

f(1)S

max x

)

(

x∈

或：)

m i n x

f(2)S

(

x∈

这里S→Rn称为搜索空间，f(x):S→Rn称为目标函数。

(1)式描述的优化问题称为极大化问题，(2)式描述的称为极小化问题。

当把f(x)看成是一序列的函数时，上述的问题就转变为多目标优化问题。

对很多实际问题进行数学建模后，可将其抽象为一个数值函数的优化问题。由于问题种类的繁多、影响因素的复杂，这些数学函数会呈现出不同的数学特征，比如连续的、离散的、凸的、凹的、单峰值的、多峰值的函数等等，经常遇到的函数还有这些不同数学特征的组合，除了在函数是连续、可求导、低阶的简单情况下可解析地求出其最优解外，大部分情况需要通过数值计算方法来进行近似优化计算。尽管人们对这个问题研究了很多年，但至今仍无一种既能处理各种不同的复杂函数、又具有良好求解结果的数值计算方法。特别是当问题的规模比较大时，优化计算时的搜索空间急剧扩大，人们认识到要严格地求出其最优解不太现实。所以需要研究出一种能够在可接受的时间和可接受的精度范围内求出数值函数近似最优解的方法或通用算法。粒子群优化由于其算法的简单，易于实现，无需梯度信息，参数少等特点在连续优化问题和离散优化问题中都表现出了良好的效果，特别是因为其天然的实数编码特点适合于处理实优化问题。近年来成为国际上智能优化领域研究的热门。

1.2 粒子群算法基本原理

粒子群优化算法PSO(Particle Swarm Optimization)是一种基于群体的自适应的搜索优化方法。是由James Kennedy和Eberhart在1995年提出的。PSO 中，

每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的函数决定的适值( fitness value) ，每个粒子还有一个速度决定它们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。

PSO 初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解，这个极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。

假设在一个D 维的目标搜索空间中，有N 个粒子组成一个群落，其中第i 个粒子表示为一个D 维的向量

)

,,,(21iD i i i x x x X =，N i ,,2,1 =。

第i 个粒子的“飞行 ”速度也是一个D 维的向量，记为

)

,,21i iD i i v v v V ，（=，3,2,1 =i 。

第i 个粒子迄今为止搜索到的最优位置称为个体极值，记为

)

,,,(21iD i i best p p p p =，N i ,,2,1 =。

整个粒子群迄今为止搜索到的最优位置为全局极值，记为

)

,,,(21gD g g best p p p g =

在找到这两个最优值时，粒子根据如下的公式(1.1)和( 1.2)来更新自己的速度和位置：

())

(2211id gd id id id id x p r c x p r c v w v -+-+*= (1.1)

id

id id v x x += (1. 2)

其中：1c 和2c 为学习因子，也称加速常数(acceleration constant)，1r 和2r 为[0，1]范围内的均匀随机数。式(1.1)右边由三部分组成，第一部分为“惯性(inertia)”或“动量(momentum )”部分，反映了粒子的运动“习惯(habit)”，代表粒子有维持自己先前速度的趋势；第二部分为“认知(cognition)”部分，反映了粒子对自身历史经验的记忆(memory)或回忆(remembrance)，代表粒子有向自身历史最佳位置逼近的趋

势；第三部分为“社会(social)”部分，反映了粒子间协同合作与知识共享的群体历史经验。

二、算法设计

这部分内容主要是针对本文主要研究问题的类型确定粒子群算法具体实现过程和一些参数的选择。

2.1算法流程框图

2.2 算法实现

算法的流程如下：

① 初始化粒子群，包括群体规模N ，每个粒子的位置i x 和速度i V ② 计算每个粒子的适应度值][i F it ；

③ 对每个粒子，用它的适应度值][i F it 和个体极值）（i p best 比较，如果)(][i p i F best it > ，则用][i Fit 替换掉）（i best p ；

④ 对每个粒子，用它的适应度值][i Fit 和全局极值b e s t g 比较，如果

)(][i p i F b e s t it >则用][i F it 替best g ；

⑤ 根据公式（1.1），（1.2）更新粒子的速度i v 和位置i x ；

⑥ 如果满足结束条件(误差足够好或到达最大循环次数)退出，否则返回②。

2.3 算法的构成要素

对于构成粒子群算法的各个参数进行设定。本算法中主要的参数变量为惯性权值w ，学习因子1c ，2c ，群体的大小 N ，迭代次数 M ，粒子维数 D 。 （1）群体大小

通常，群体太小则不能提供足够的采样点，以致算法性能很差，容易陷入局部最优；群体太大尽管可以增加优化信息，阻止早熟收敛的发生，但无疑会增加计算量，造成收敛时间太长，表现为收敛速度缓慢。本文对函数的优化选择种群规模为500。

（2）最大迭代次数

迭代次数越多能保证解的收敛性，但是影响运算速度，本文对函数的优化选最大迭代次数为1000次。 （3）惯性权值

惯性权重w 表示在多大程度上保留原来的速度。w 较大，全局收敛能力强，局部收敛能力弱；w 较小，局部收敛能力强，全局收敛能力弱。 （4）学习因子

加速常数1c 和2c 分别用于控制粒子指向自身或邻域最佳位置的运动。建议

0.421≤+=c c φ，并通常取221==c c 。本文中取221==c c 。

（5）粒子维数

粒子维数取决于待优化函数的维数，例如本文主要有三个函数主要:第一个函数是2维的，第二个函数是10维的，第三个函数是30维的。

（6）粒子空间的初始化

较好的选择粒子初始化空间，将大大缩短收敛的时间。在本文中我们主要是选用随机对粒子进行初始化。

2.4 算法的改进

对于函数的优化我们主要选择的是对于惯性权重的优化。惯性权重是粒子优化算法的重要参数，算法的成败很大程度上就取决于参数的选取和调节，本文采用固定权重、时变权重和随机权重三种权重。 （1）固定权重

即赋予惯性权重以一个常数值，一般来说，该值在0和1之间。固定的惯性权重使粒子在飞行中始终具有相同的探索和开发能力。显然对于不同的问题，获得最好优化效果的常数是不同的，要找到这个值需要大量的实验。通过实验发现：种群规模越小，需要的惯性权重越大，因为此时种群需要更好的探索能力来弥补粒子数量的不足，否则粒子极易收敛；种群规模越大，需要的惯性权重越小，因为每个例子可以更专注于搜索自己附近的区域。 （2）时变权重

希望粒子群在飞行开始的时候具有较好的探索能力，随着迭代次数的增加，特别是在飞行后期，希望有较好的开发能力。所以使用动态调节惯性权重。可以通过时变的惯性权重来实现。设惯性权重的取值范围为：],[max min ωω，最大迭代次数为max _Iter ，则第i 次迭代时的惯性权重可以为：

i Iter i ?--

=max

_min

max max ωωωω

（3）随机权重

随机权重是在一定范围内随机取值，在本文中我们采用的是：

2

5.0Random

+

=ω

其中Random 为0-1之间的随机数。这样惯性指数将在0.5-1之间随机变化。均值为0.75。对于动态优化问题来说，不能够预测在给定的时间粒子群于要更好的探索能力还是更好的开发能力。所以，可以使惯性权重在一定范围内随机变化。 需要说明的是，本文的程序允许改变除惯性权重以外的其他参数，因为本文编写的程序参照MATLAB 工具箱，留给用户解决这类问题一个接口函数，上述的各个参数正是接口函数的参数，因此允许改变。另外对于c 也可采用变参数法，即随迭代次数增加，利用经验公式使它们动态调整，本文采用固定值。

三、算例设计

3.1 测试函数介绍

本文主要选取三个函数：一个2维连续函数、一个10离散函数、一个30维离散函数，利用MATLAB 编写粒子群算法程序来优化它们。本文选取了三个函数，分别如下：

()()()

2

22

212

22

2

11001.015.0sin

5.0x x

x

x f ++-+-

=

)10)2cos(10(210

1

2+-=∑=i i i x x f π

∑∏==+-=

301

30

1

231)cos(

4000

1i i i i i

x x f

其中f1求最大值，f2和f3为求最小值。

3.2 优化函数特点

(1)Schaffer 函数：()()()

2

22

212

22

211001.015.0sin

5.0x x

x

x f ++-+-

=求其最大值。

目标函数的效果图如图3.1下：

图3.1 Schaffer 函数的效果图

由图知此函数是个二维函数，常用于测试粒子群算法性能的测试函数，全局在（0，0）处取得最大值，具有强烈震荡的状态，而在Xi （-3.14,3.14）范围内，有无限个次全局最大点。

（2）Rastrigrin 函数：)10)2cos(10(210

12+-=∑=i i i x x f π

目标函数的效果图如图3.2所示：

图3.2 Rastrigrin 函数的效果图

由图知此函数是10维多峰值函数，存在大量按正弦拐点排列的、很深的局部最优点。其在（0,0,…0）处取得全局最小值，在Xi ∈（-5.12,5.12）范围内大约有10个局部极小点，不难优化查找到全局最优值。

（3）Girewank 函数：∑∏==+-=

30

1

30

1

231)cos(

4000

1i i i i

i x x f

目标函数的效果图如图3.3所示：

图3.3 Girewank 函数的效果图

由图知此函数为一个多峰值函数，为30维函数，变量之间有相互关系，该函数有很多局部最优点，其全局最优点全局最小值在（X 1，X 2，……Xn ）=（0,0,……0）取得，目标函数最优为0。

四、仿真实验设计

4.1 实验参数设计

（1）权重参数设计

本文对每个测试函数将使用三种不同的惯性权重策略进行实验，分别为固定

权重：6.0=ω；随机权重：2

5.0Random +=ω；时变权重：

i Iter i ?--

=max

_min

max max ωωωω，其中2.09.0min max ==ωω，。

（2）其他参数设计

对于Schaffer 函数,种群规模N=500，最大迭代次数M=1000，学习因子

221==c c ，1f 函数是2维的，粒子维数取21

=D 。对于Rastrigrin 函数,种群规模

N=500，最大迭代次数M=1000，学习因子221==c c ，2f 函数是10维的，粒子维数。102=D 对于Girewank 函数，种群规模N=500，最大迭代次数M=1000，学习因子221==c c ，3f 函数是30维，粒子维数303=D 。

4.2 基本粒子群算法在测试函数中应用

以Schaffer 函数、Rastrigrin 函数、Girewank 函数为例来说明基本粒子群算法在函数优化中问题中的效果。

步骤1：首先依据基本粒子群算法的流程图编写粒子群算法的MATLAB 程序（见附录一）；

步骤2：分别编写各个函数的适应值函数程序（见附录二）；

步骤3：在MATLAB 环境中用基本粒子群算法分别调用三个函数的适应值程序，得到收敛效果图。

Schaffer 函数、Rastrigrin 函数、Girewank 函数的收敛效果图分别如图4.1、4.2、4.3所示：

0.9988

0.9992

0.9994

0.9996

0.9998

图4.1 Schaffer 函数收敛效果图

图4.2 Rastrigrin 函数收敛效果图

0.01

0.02

0.03

0.04

0.05

0.06

0.07

图4.3 Girewank 函数收敛效果图

五、仿真实验结果分析

5.1 实验结果汇总

表1统计了固定迭代次数为1000次时，三个测试函数分别在固定权重、随机权重和时变权重三种惯性权重策略下求解的平均值和最差的的一次求解结果。

表1 算法运行100次搜索到的解的平均值和最差值

惯性权重 Schaffer Griewank Rastrigrin 平均值

最差值

平均值 最差值

平均值

最差值

固定权重 6.0=ω 1

1

8.26

16.91 1.82e-006 2.04e-005

随机权重

25.0Random

+

=ω 1 1 6.21 15.92 2.14e-006 1.74e-005

时变权重

i Iter i ?--

=max

_min

max max ωωωω

1 1 7.7

2 18.9 3.20e-009 2.16e-008

5.2 实验结果分析

惯性权重是粒子群优化算法标准版本的重要参数，算法的成败很大程度上取决于该参数的选取和调节。本文重点研究了惯性权重对优化效果的影响。从表1可以看出，（1）在其他参数选择适当的条件下，利用三种不同惯性权重策略的PSO 算法得到的平均值和最差值相差不太，较为稳定。（2）固定权重、随机权重以及时变权重对Schaffer 函数都具有较好的优化效果，并能较快地迭代到最优值。（2）对于Rastrigrin 函数，三种惯性权重策略下都未能搜索到目标函数的理想最优点，但相比之下利用随机权重进行函数优化具有较好的效果，时变权重并未取得比固定权重好的优化效果。（3）对于多峰函数Girewank 函数，惯性权重采用时变权重优化效果明显比固定权重、时变权重好，取得了不错的收敛结果。

六、总结与展望

6.1总结

本文主要用粒子群算法优化不同维数的连续函数以及离散函数的最小值问题，主要有以下几个方面：

⑴首先介绍了粒子群的算法在智能优化中的地位，也介绍了粒子群算法的主要特点，并通过对粒子群算法的学习和了解为下面粒子群改进算法在对不同维数的函数优化问题打下了很好的基础。

⑵其次对粒子群算法解决最优化问题的统一框架进行了分析，在此基础上提出了粒子群优化算法的设计步骤；又对粒子群优化算法的原理进行了分析，从而对粒子群算法有了更深刻的了解。

⑶最后将基本的粒子群算法与改进的粒子群算法分别对函数的优化问题在MATLAB中进行了仿真，从而将基本算法和改进算法在函数优化问题中的仿真结果进行了对比，从而验证了改进算法的相对优越性，并且验证了改进算法的实际可行性。

6.2展望

本文中主要利用改变权重的方法来对函数进行优化，实际用粒子群算法对函数进行优化还有很多方法：改变邻域拓扑结构、改变学习因子，对于离散函数的优化可以采用二进制编码和顺序编码来实现同时也可以使用基于遗传策略和梯度信息的集中改进方法例如：基于选择的改进算法。基于交叉的改进算法、基于变异的改进算法、带有梯度加速的的改进算法，同时智能优化方法处理约束的一般性策略都可以借鉴到粒子群算法中，也可以根据粒子群优化的特性来设计专门的约束处理方式，可以通过以上的优化方法来解决函数优化问题。

附录一

主函数

%fitness-是要优化的目标函数，N-种群数，c1,c2-学习因子，w-惯性权重，M-迭代次数，D-粒子维数。

format long;

%初始化种群

%c1=2;

%c2=2;

%N=500;

%w=0.6;

%M=1000;

c1=2;

c2=2;

N=500;

%w=0.6;

M=1000;

D=30;

x=randn(N,D); %随机初始化位子y=randn(N,D);

v=randn(N,D); %随机初始化速度

p=rand(N,1);

%先计算各个粒子的适应度，并初始化pi-粒子个体极值和pg-全局极值

for i=1:N

p(i)=fitness(x(i,:));

y(i,:)=x(i,:);

end

pg = x(N,:); %pg为全局极值

for i=1:(N-1)

if fitness(x(i,:))<(pg)

pg=x(i,:);

end

end

%进入粒子群算法主要循环

for t=1:M

%wmin=0.2;wmax=0.9;w=wmax-(wmax-wmin)/M*t

for i=1:N

w=0.5*randn(0,1)/2;

v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));

x(i,:)=x(i,:)+v(i,:);

if fitness(x(i,:))

p(i)=fitness(x(i,:));

y(i,:)=x(i,:);

end

if p(i)

pg=y(i,:);

end

end

Pbest(t)=fitness(pg);

end

xm = pg';

fv = fitness(pg);

t=[1:M];

plot(t,Pbest)

%------最后给出计算结果

disp('*************************************************************') %Solution=pg';

disp(pg')

disp(fitness(pg))

%Result=fitness(pg,D);

disp('*************************************************************') %------算法结束---DreamSun GL & HF-----------------------------------

附录二

三个测试函数的适应值函数程序

function f=fitness(x)

%f=sum(x.^2);

%f=0.5-(sin(sqrt(sum(x.^2)))^2-0.5)/(1+0.001*(sum(x.^2)))^2; %f=sum(x.^2-10*cos(pi*2*x)+10);

f1=1;

for i=1:30

f1=f1*cos(x(i)/sqrt(i));

end

f=1/4000*sum(x.^2)-f1+1;

end