粒子群优化算法（PSO算法）

  粒子群优化算法（Particle Swarm Optimization, PSO算法）是一种进化计算技术，由美国电气工程师Eberhart博士和社会心理学家Kennedy博士发明，源于对鸟群捕食的行为研究。PSO算法同遗传算法类似，是一种基于迭代的优化工具。系统初始化为一组随机解，通过迭代搜寻最优值。但是并没有遗传算法用的交叉以及变异，而是粒子在解空间追随最优的粒子进行搜寻。

粒子群算法的发展

  粒子群算法的发展始于1995年Eberhart和Kennedy¹提出的基本粒子群算法。其中基本PSO的参数是固定的，在对某些函数优化上精度较差，且值得注意的是，在基本PSO中，并没有惯性权重$w$。
  后来Shi等²在1998年提出了惯性权重$w$（惯性因子）线性递减的改进算法，使算法在搜索初期具有较大搜索能力，而在后期又能够得到较准确的结果，此改进大大提高了基本PSO算法的性能。逐渐地大家都默认这个改进粒子群算法为标准的粒子群算法。2001年Shi又提出了自适应模糊调节$w$ 的PSO，在对单峰函数的处理中取得了良好的效果。(史玉回教授简介）
  Kennedy等³研究粒子群的拓扑结构，分析粒子间的信息流，提出了一系列的拓扑结构。
  Angeline⁴将选择算子引入到PSO中，选择每次迭代后的较好地粒子并复制到下一代，以保证每次迭代的粒子群都具有较好地性能。
  Higashi等⁵分别提出了自己的变异PSO算法，基本思想均是希望通过引入变异算子跳出局部极值点的吸引，从而提高算法的全局搜索能力，得到较高的搜索成功率。
  Al-Kazemi⁶所提出的Multi-Phase PSO 在粒子群中随机选取部分个体向全局最优飞，而其它个体向反方向飞，以扩大搜索空间。
  除以上的混合算法之外，还出现了量子PSO、模拟退火PSO、耗散PSO、自适应PSO等混合改进算法，也有采取PSO与基于梯度的优化方法相结合的办法等。粒子群算法概念简单，具有很强的发现较好解的能力，并不容易陷入局部最优解。为了进一步完善粒子群算法，许多学者已经提出了很多改进方法，综合来看，其改进模式主要可以分为以下3类：

1. 算法的具体参数设定和调整策略：如Shi和Eberhart在速度项添加的惯性权重，和后来提出的模糊自适应调整惯性权重的思想。
2. 算法的总体结构：粒子群的总体结构和组织模式，如协同粒子群算法和共同进化粒子群算法等。此外，群体的拓扑结构与粒子群算法的性能也备受关注。群体的拓扑结构主要有星形结构、环形结构和车轮形结构等等。不同的拓扑结构对应于不同的信息交换方式，从而进一步影响了算法的寻优能力。
3. 混合算法：粒子群算法和其它智能算法的结合。如基于自然选择的算法，基于杂交的算法，基于模拟退火的算法。

粒子群算法的概述

  如前文所述，目前被大家所广为接受的是Shi等人在1998年所改进的粒子群算法，即含有惯性权重的粒子群算法。故我将以带有惯性权重的粒子群算法作为基本粒子群算法，下文不再说明。
  在PSO算法中，每个优化问题的解都是搜索空间中的一只鸟，被抽象为没有质量和体积的微粒，并将其延伸到$N$维空间里的位置表示为一个向量，每个粒子的飞行速度也表示为一个向量。
  所有的粒子都有一个由被优化的函数决定的适应值（fitness），每个粒子还有一个速度决定它们飞翔的方向和距离。假设粒子们知道自己到目前为止发现的最好位置（pbest)和现在的位置，这个可以看做是粒子自己的飞行经验。除此之外，每个粒子还知道目前为止整个群体中所有粒子发现的最好位置（gbest,gbest是pbest中的最好值），这个可以看做是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
  PSO算法首先初始化一群随机粒子（随机解)，然后粒子们就追随当前的最优粒子在解空间中搜索，即通过迭代找到最优解。假设$d$维搜索空间（即待优化函数有$d$个自变量）中的第$i$个粒子的位置和速度分别为$X_i=(x_{i,1},x_{i,2},...,x_{i,d})$和$V_i=(v_{i,1},v_{i,2},...,v_{i,d})$，在每一次迭代过程中，粒子通过追踪两个最优解来更新自己，第一个就是粒子本身所找到的最优解，即个体极值pbest，$P_i=(p_{i,1},p_{i,2},...,p_{i,d})$；另一个就是整个种群目前所找到的最优解，即全局最优解gbest，$P_g$。在找到这两个最优值时，粒子根据如下的公式来更新自己的速度和位置。
v i , j ( t + 1 ) = w v i , j ( t ) + c 1 r 1 [ p i , j − x i , j ( t ) ] + c 2 r 2 [ p g , j − x i , j ( t ) ] x i , j ( t + 1 ) = x i , j ( t ) + v i , j ( t + 1 ) , j = 1 , 2 , . . . , d

$$\begin{array}{lcr} v_{i,j}(t+1)=wv_{i,j}(t)+c_1r_1[p_{i,j}-x_{i,j}(t)]+c_2r_2[p_{g,j}-x_{i,j}(t)]\\ x_{i,j}(t+1)=x_{i,j}(t)+v_{i,j}(t+1), j=1,2,...,d \end{array}$$ vi,j​(t+1)=wvi,j​(t)+c1​r1​[pi,j​−xi,j​(t)]+c2​r2​[pg,j​−xi,j​(t)]xi,j​(t+1)=xi,j​(t)+vi,j​(t+1),j=1,2,...,d​其中，$w$为惯性权重，最开始PSO算法是没有这个参数的，$c_1,c_2$为正的学习因子，$r_1,r_2$为0到1之间均匀分布的随机数，$t$为第$t$轮迭代。
  粒子群算法的性能很大程度上取决于算法的控制参数，例如粒子数、最大速度、学习因子、惯性权重等，各个参数的选取原则如下。

• 粒子数：粒子数的多少根据问题的复杂程度自行决定。对于一般的优化问题，20至40个粒子就完全可以得到很好的结果；对于比较简单的问题，10个粒子已经足够可以取得较好地结果；对于比较复杂的问题或者特定类别的问题，粒子数可以取到100以上。在MATLAB所带的库函数（particleswarm）中，关于粒子个数的选取，默认是$min(100,10\ast d)$,这里的$d$是粒子的维度。
• 粒子的维度：这是由优化问题决定的，就是问题解的维度。
• 粒子的范围：由优化问题决定，每一维可设定不同的范围。注：如果需要对解的范围做限定，可以在这里进行入手。
• 最大速度$V_{max}$：决定粒子在一个循环中最大的移动距离，通常设定为粒子的范围宽度。
• 学习因子：学习因子是粒子具有自我总结和向群体中优秀个体学习的能力，从而向群体内或邻域内最优点靠近，通常取$c_1$和$c_2$为2，但也有其他的取值，一般$c_1$等于$c_2$，且范围在0至4之间。
• 惯性权重：决定了对粒子当前速度继承的多少，合适的选择可以使粒子具有均衡的探索能力和开发能力，惯性权重的取法一般有常数法、线性递减法、自适应法等。

粒子群算法的步骤与实现

算法步骤

1. 初始化一群微粒（群体规模为N），包括随机位置和速度；
2. 评价每个微粒的适应度，将当前各微粒的位置和适应度值存储在各微粒的pbest中，将所有pbest中适应度值最优个体的位置和适应度值存储于gbest中；
3. 更新粒子的速度和位置：
   $$\begin{array}{l}{v}_{i,j}(t+1)=w{v}_{i,j}(t)+c_1{r}_1[p_{i,j}-x_{i,j}(t)]+c_2{r}_2[p_{g,j}-x_{i,j}(t)]\\ x_{i,j}(t+1)=x_{i,j}(t)+v_{i,j}(t+1),j=1,2,...,d\end{array}$$
4. 对每个微粒，将其适应值与其经过的最好位置作比较，如果较好，则将其作为当前的最好位置；
5. 比较当前所有pbest和gbest的值，更新gbest;
6. 若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则返回3继续搜索。

算法的流程图如下图所示：

算法的实现（MATLAB）

function [xm,fv] = PSO(fitness,N,c1,c2,w,M,D)
%待优化的目标函数：fitness
%粒子数目：N
%学习因子1、2：c1,c2
%惯性权重：w
%最大迭代次数：M
%自变量的个数：D
%目标函数取最小值时的自变量值：xm
%目标函数的最小值：fv

format long;
for i=1:N
	for j=1:D
		x(i,j)=randn;	%随机初始化位置
		v(i,j)=randn;	%随机初始化速度
	end
end
%粒子有N个，有D维。
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,:))<fitness(pg)
		pg=x(i,:);
	end
end

for t=1:M
	for i=1:N 	%速度、位移更新
		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)
			p(i)=fitness(x(i,:));
			y(i,:)=x(i,:);
		end
		if p(i)<fitness(pg)
			pg=y(i,:);
		end
	end
	Pbest(t)=fitness(pg);
end
xm=pg';
fv=fitness(pg);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

调用示例

• 编写待优化函数

function F = fitness(x)
F = 0;
for i=1:30
	F=F+x(i)^2;
end
1
2
3
4
5

• 在MATLAB命令窗口中输入

[xm,fv] = PSO(@fitness,40,2,2,0.5,5000,30)
1

关于惯性权重⁷

  为了改善基本粒子群算法的收敛性能，Shi和Eberhart在1998年的IEEE国际进化计算学术会议上发表了题为“A modified particle swarm optimizer”的论文，引入惯性权重，逐渐地大家都默认这个改进粒子群算法为标准的粒子群算法。
  在基本粒子群算法的速度公式上可见右边项包括了三部分：第一部分是粒子之前的速度；第二部分和第三部分是粒子对速度的调整。如果没有后面两部分，粒子将会保持相同的速度朝一个方向飞行，直到到达边界，这样粒子很大可能会找不到最优解，除非最优解在粒子飞行的轨迹上，但这种情况是很少的。此外，如果没有第一部分，粒子的飞行速度将仅由它们当前位置和历史最好位置决定，则速度自身是无记忆的。假定刚开始粒子$i$处于较优位置，那么粒子的飞行速度将会是0，即它会保持静止状态，知道其它粒子找到比粒子$i$所处位置还要好的解，从而替代了全局最优。此时，每个粒子将会飞向它自身最好位置和群体全局最好位置的权重中心。所以可以想象到如果没有第一部分，粒子群算法的搜索空间将会随着进化而收缩。此时只有当全局最优在初始搜索区间时，粒子群算法才可能找到解。所以最后解非常依赖于初始群体。当没有第一部分时，此算法更像是局部最优算法。
  对于不同的问题，局部最优能力和全局最优能力的权衡也不一样。考虑到这个问题，并结合以上的讨论，Shi和Eberhart添加了一个惯性权重到速度更新公式，即
$$v_{i,j}(t+1)=w{v}_{i,j}(t)+c_1{r}_1[p_{i,j}-x_{i,j}(t)]+c_2{r}_2[p_{g,j}-x_{i,j}(t)]$$
惯性权重$w$起着权衡局部最优能力和全局最优能力的作用。
  为了观察这个惯性权重对粒子群算法性能的影响,Shi 和 Eberhart 把此算法应用到 Schaffer’s f6 函数中（该函数是比较著名的评价优化算法的基准函数）。他们改变惯性权重的大小，通过大量的实验得到一些结论。当惯性权重较小时(<0.8)，如果粒子群算法能找到全局最优的话，那么它所经历的搜索时间是很短的，即所有的粒子趋向于快速汇集在一起。如果优解是在初始搜索空间内，粒子群算法将会很容易找到全局最优，否则它会找不到全局最优。当惯性权重较大时(>1.2)，粒子群算法更像全局搜索方法，且它总是探索新的区域。当然，这时的粒子群算法会需要更多的迭代来达到全局最优，且更有可能找不到全局最优当惯性权重适中时，粒子群算法将会有更大的机会找到全局最优，但迭代次数也会比第一种情况要多。
  根据这些分析，他们不是把惯性权重设为定值，而是设为一个随时间线性减少的函数（线性递减权重法），惯性权重的函数形式通常为
$$w=w_{max}-\frac{w_{max}-w_{min}}{t_{max}}\ast t$$
其中，$w_{max}$为初始权重；$w_{min}$为最终权重（即让惯性权重从初始权重-最大值，线性减小到最终权重-最小值）。$t$表示当前迭代步数，$t_{max}$表示最大迭代步数，通常取$w_{max}=0.9$，$w_{min}=0.4$。
  这个函数使得粒子群算法在刚开始的时候倾向于开掘(exploitation)，然后渐转向于开拓(exploration)，从而在局部区域调整解。这些改进使得粒子群算法的性能得到很大的提高。此外，针对惯性权重的改进，还有如下方法。

自适应权重法

  为了平衡PSO算法的全局搜索能力和局部改良能力，还可采用非线性的动态惯性权重系数公式，其表达式如下：
w = { w min ⁡ − ( w max ⁡ − w min ⁡ ) ∗ ( f − f min ⁡ ) ( f a v g − f min ⁡ ) , f ≤ f a v g w max ⁡ , f > f a v g

$$\begin{aligned} w=\left\{\begin{matrix}w_{\min}-\frac{\left(w_{\max}-w_{\min}\right)*\left(f-f_{\min}\right)}{\left(f_{a v g}-f_{\min}\right)},f\leq f_{avg}\\ w_{\max},f>f_{avg} \end{matrix}\right. \end{aligned}$$ w={wmin​−(favg​−fmin​)(wmax​−wmin​)∗(f−fmin​)​,f≤favg​wmax​,f>favg​​​
其中，$w_{\max },w_{\min }$分别表示$w$的最大值和最小值，$f$表示粒子当前的目标函数值，$f_{avg}$和$f_{\min }$分别表示当前所有微粒的平均目标值和最小目标值。在上式中，惯性权重随着微粒的目标函数值而自动改变，因此称为自适应权重。
  当各微粒的目标值趋于一致或者趋于局部最优时，将使惯性权重增加，而各微粒的目标值比较分散时，将使惯性权重减小，同时对于目标函数值优于平均目标值的微粒，其对应的惯性权重因子较小，从而保护了该微粒，反之对于目标函数值差于平均目标值的微粒，其对应的惯性权重因子较大，使得该微粒向较好的搜索区域靠拢。

随机权重法

  将标准 PSO算法中设定 $w$ 为服从某种随机分布的随机数，这样一定程度上可从两方面来克服 $w$ 的线性递减所带来的不足。
  首先，如果在进化初期接近最好点，随机 $w$ 可能产生相对小的 $w$ 值，加快算法的收敛速度，另外，如果在算法初期找不到最好点，$w$ 的线性递减，使得算法最终收敛不到此最好点，而 $w$ 的随机生成可以克服这种局限。
  $w$ 的计算公式如下:
{ w = μ + σ ∗ N ( 0 , 1 ) μ = μ min ⁡ + ( μ max ⁡ − μ min ⁡ ) ∗ r a n d ( 0 , 1 ) \left\{

$$\begin{matrix}w=\mu+\sigma^{*}N\left(0,1\right)\\ \mu=\mu_{\min}+\left(\mu_{\max}-\mu_{\min}\right)*r a n d\left(0,1\right)\end{matrix}$$ \right. {w=μ+σ∗N(0,1)μ=μmin​+(μmax​−μmin​)∗rand(0,1)​
其中N(0,1)表示标准正态分布的随机数，rand(0,1)表示0到1之的随机数，${\mu }_{\max }$为随机权重平均值的最大值，${\mu }_{\min }$为随机权重平均值的最小值，$\sigma$为随机权重的方差。一般情况下，可取${\mu }_{\max }=0.8,{\mu }_{\min }=0.4,\sigma =0.2$。

关于学习因子⁸

  学习因子是粒子具有自我总结和向群体中优秀个体学习的能力，从而向群体内或邻域内最优点靠近。学习因子一般固定为常数，并且取值为2。但是在实际的应用中，也有一些其他的取值方式，常见的有同步变化和异步变化的学习因子。

同步变化的学习因子

  同步变化的学习因子指的是将学习因子$c_1$和$c_2$的取值范围设定为$[c_{min},c_{max}]$，第$t$次迭代时的学习因子取值公式为：
$$c_1=c_2=c_{max}-\frac{c_{max}-c_{min}}{t_{max}}\ast t$$
其中，$t$为当前迭代次数，$t_{max}$为总的迭代次数。

异步变化的学习因子

  两个学习因子在优化过程中随时间进行不同的变化称之为异步变化的学习因子，这样使得在优化的初始阶段，粒子具有较大的自我学习能力和较小的社会学习能力，加强全局搜索能力。而在优化的后期，粒子具有较大的社会学习能力和较小的自我学习能力，有利于收敛到全局最优解。学习因子的变化公式为:
$$c_1=c_{1,ini}+\frac{c_{1,fin}-c_{1,ini}}{t_{max}}\ast t,c_2=c_{2,ini}+\frac{c_{2,fin}-c_{2,ini}}{t_{max}}\ast t$$
其中，$c_{1,ini}、c_{2,ini}$分别代表$c_1$和$c_2$的初始值，$c_{1,fin}、c_{2,fin}$代表$c_1$和$c_2$的迭代终值。对于大多数情况下，采用如下的参数设置效果较好：
$$c_{1,ini}=2.5，c_{1,fin}=0.5,c_{2,ini}=0.5，c_{2,fin}=2.5$$

参考文献

---

1. Kennedy J, Eberhart R. Particle swarm optimization[C]//Proceedings of ICNN’95-international conference on neural networks. IEEE, 1995, 4: 1942-1948. ↩︎

2. Shi Y, Eberhart R. A modified particle swarm optimizer[C]//1998 IEEE international conference on evolutionary computation proceedings. IEEE world congress on computational intelligence (Cat. No. 98TH8360). IEEE, 1998: 69-73. ↩︎

3. Mendes R, Kennedy J, Neves J. The fully informed particle swarm: simpler, maybe better[J]. IEEE transactions on evolutionary computation, 2004, 8(3): 204-210. ↩︎

4. Angeline P J. Using selection to improve particle swarm optimization[C]//1998 IEEE International Conference on Evolutionary Computation Proceedings. IEEE World Congress on Computational Intelligence (Cat. No. 98TH8360). IEEE, 1998: 84-89. ↩︎

5. Higashi N, Iba H. Particle swarm optimization with Gaussian mutation[C]//Proceedings of the 2003 IEEE Swarm Intelligence Symposium. SIS’03 (Cat. No. 03EX706). IEEE, 2003: 72-79. ↩︎

6. Al-Kazemi B, Mohan C K. Multi-phase generalization of the particle swarm optimization algorithm[C]//Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No. 02TH8600). IEEE, 2002, 1: 489-494. ↩︎

7. 龚纯,王正林. 精通MATLAB最优化计算（第四版）[M]. 电子工业出版社,2016 ↩︎

8. 纪震,廖惠连,吴青华. 粒子群算法及应用[M]. 科学出版社,2009 ↩︎