粒子群优化算法

一、概述

粒子群优化算法（Particle Swarm Optimization，PSO）的思想来源于对鸟捕食行为的模仿，最初，Reynolds.Heppner 等科学家研究的是鸟类飞行的美学和那些能使鸟群同时突然改变方向，分散，聚集的定律上，这些都依赖于鸟的努力来维持群体中个体间最佳距离来实现同步。而社会生物学家 E.O.Wilson 参考鱼群的社会行为认为从理论上说，在搜寻食物的过程中，尽管食物的分配不可知，群中的个体可以从群中其它个体的发现以及以往的经验中获益。

粒子群从这种模型中得到启发并用于解决优化问题。如果我们把一个优化问题看作是在空中觅食的鸟群，那么粒子群中每个优化问题的潜在解都是搜索空间的一只鸟，称之为“粒子”（Particle），“食物”就是优化问题的最优解。每个粒子都有一个由优化问题决定的适应度用来评价粒子的“好坏”程度，每个粒子还有一个速度决定它们飞翔的方向和距离，它根据自己的飞行经验和同伴的飞行经验来调整自己的飞行。粒子群初始化为一群随机粒子（随机解），然后通过迭代的方式寻找最优解，在每一次的迭代中，粒子通过跟踪两个“极值”来更新自己，第一个是粒子本身所经历过的最好位置，称为个体极值即$Pbest$；另一个是整个群体经历过的最好位置称为全局极值$Gbest$。每个粒子通过上述的两个极值不断更新自己，从而产生新一代的群体。

二、粒子群算法

算法的描述如下：

假设搜索空间是$L$维，并且群体中有$N$个粒子。那么群体中的第$i$个粒子可以表示为一个$L$维的向量，$X_i=(x_{i1},x_{i2},\cdots ,x_{iL}),i=1,2,\cdots ,N$，即第$i$个粒子在$L$维的搜索空间的位置是$X_i$，它所经历的“最好”位置记作$Pbes{t}_i=(p_{i1},p_{i2},\cdots ,p_{iL}),i=1,2,\cdots ,N$。粒子的每个位置代表要求的一个潜在解，把它代入目标函数就可以得到它的适应度值，用来评判粒子的“好坏”程度。整个群体迄今为止搜索到的最优位置记作$Gbes{t}_g=(p_{g1},p_{g2},\cdots ,p_{gL})$，$g$是最优粒子位置的索引。

$$V_i^{t+1}=\omega V_i^t+c_1{r}_1(Pbes{t}_i^t-X_i^t)+c_2{r}_2(Gbes{t}_g^t-X_i^t)\text{\hspace{1em}(1)}$$

$$X_i^{t+1}=X_i^t+V_i^{t+1}\text{\hspace{1em}(2)}$$

$\omega$为惯性权重（inertia weight），$Pbes{t}_i^t$为第$i$个粒子到第$t$代为止搜索到的历史最优解，$Gbes{t}_g^t$为整个粒子群到目前为止搜索到的最优解，$X_i^t$，$V_i^t$分别是第$i$个粒子当前的位置和飞行速度，$c_1,c_2$为非负的常数，称为加速度因子，$r_1,r_2$是$[0,1]$之间的随机数。

公式由三部分组成，第一部分是粒子当前的速度，表明了粒子当前的状态；第二部分是认知部分（Cognition Modal），表示粒子本身的思考（$c_1$也称为自身认知系数）；第三部分是社会认知部分（Social Modal），表示粒子间的信息共享（$c_2$为社会认知系数）。

参数的选择：

粒子数目一般取30~50，参数$c_1,c_2$一般取2。适应度函数、粒子的维数和取值范围要视具体问题而定。问题解的编码方式通常可以采用实数编码。

算法的主要步骤如下：

第一步：对粒子群的随机位置和速度进行初始设定，同时设定迭代次数。

第二步：计算每个粒子的适应度值。

第三步：对每个粒子，将其适应度值与所经历的最好位置$Pbes{t}_i$的适应度值进行比较，若较好，则将其作为当前的个体最优位置。

第四步：对每个粒子，将其适应度值与全局所经历的最好位置$Gbes{t}_g$的适应度值进行比较，若较好，则将其作为当前的全局最优位置。

第五步：根据公式(1),(2)对粒子的速度和位置进行优化，从而更新粒子位置。

第六步：如未达到结束条件（通常为最大循环数或最小误差要求），则返回第二步。

三、基于粒子群算法的非线性函数寻优

本案例寻优的非线性函数为

$$y=-c\times \exp (-0.2\sqrt{\frac{1}{n}\sum _{i=1}^n{x}_i^2})-\exp (\frac{1}{n}\sum _{i=1}^n\cos 2\pi x_i)+c+e$$

当$c=20$，$n=2$时，该函数为Ackley函数，函数图形如图1所示。

从函数图形可以看出，该函数有很多局部极小值，在$(0,0)$处取到全局最小值0。

本案例群体的粒子数为20，每个粒子的维数为2，算法迭代进化次数为100。加速度因子$c_1=1.4$，$c_2=1.5$，惯性权重$\omega =0.1$。

引入的模块如下：

import numpy as np
import matplotlib.pyplot as plt# 解决中文乱码和负号问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

适应度函数代码如下：

# 函数用于计算粒子适应度值
def funAckley(x):c = 20y = -1.0 * c * np.exp(-0.2 * np.sqrt((x[0] ** 2 + x[1] ** 2) / 2)) \- np.exp((np.cos(2 * np.pi * x[0]) + np.cos(2 * np.pi * x[1])) / 2) \+ c + np.exp(1)return y

PSO算法代码如下：

c1, c2 = 1.4, 1.5  # 加速度因子
maxgen = 100  # 进化次数
sizepop = 20  # 群体规模
w = 0.1  # 惯性权重
Vmax, Vmin = 1, -1  # 速度最大值，最小值
popmax, popmin = 5, -5  # 个体最大值，最小值
pop = np.zeros([sizepop, 2])  # 种群
V = np.zeros([sizepop, 2])  # 速度
fitness = np.zeros(sizepop)  # 适应度
trace = np.zeros(maxgen)  # 结果# 随机产生一个群体,初始粒子和速度
for i in range(sizepop):pop[i] = 5 * ((-2 * np.random.rand(1, 2) + 1))V[i] = (-2 * np.random.rand(1, 2) + 1)fitness[i] = funAckley(pop[i])# 个体极值和群体极值
bestfitness = np.min(fitness)
bestindex = np.argmin(fitness)
Gbest = pop[bestindex]  # 全局最佳
fitnessGbest = bestfitness  # 全局最佳适应度值
Pbest = pop.copy()  # 个体最佳
fitnessPbest = fitness.copy()  # 个体最佳适应度值# 迭代寻优
for i in range(maxgen):for j in range(sizepop):V[j] = w * V[j] + c1 * (-2 * np.random.rand(1, 2) + 1) * (Pbest[j] - pop[j]) \+ c2 * (-2 * np.random.rand(1, 2) + 1) * (Gbest - pop[j])  # 速度更新V[j, V[j] > Vmax] = VmaxV[j, V[j] < Vmin] = Vminpop[j] = pop[j] + V[j]  # 群体更新pop[j, pop[j] > popmax] = popmaxpop[j, pop[j] < popmin] = popminfitness[j] = funAckley(pop[j])  # 适应度值if fitness[j] < fitnessPbest[j]:Pbest[j] = pop[j].copy()fitnessPbest[j] = fitness[j]if fitness[j] < fitnessGbest:Gbest = pop[j].copy()fitnessGbest = fitness[j]trace[i] = fitnessGbestprint(Gbest, fitnessGbest)# 结果分析
plt.plot(trace)
plt.title("最优个体适应度")
plt.xlabel("进化代数")
plt.ylabel("适应度")
plt.show()

算法结果如下：

最终得到的个体适应度值为$4.660641074893945\times 10^{-05}$，对应的粒子位置为$(5.49894607\times 10^{-06},1.55305207\times 10^{-05})$，PSO算法寻优得到最优值接近函数实际最优值，说明PSO算法具有较强的函数极值寻优能力。

四、基于自适应变异粒子群算法的非线性函数寻优

本案例寻优的非线性函数（Shubert函数）为

$$f(x_1,x_2)=\sum _{i=1}^{5}i\cos [(i+1)x_1+i]\times \sum _{i=1}^{5}icos[(i+1)x_2+i],-10\le x_1,x_2\le 10$$

该函数图形如图3所示：

从函数图形可以看出，该函数有很多局部极小值，很难用传统的梯度下降方法进行全局寻优。

自适应变异是借鉴遗传算法中的变异思想，在PSO算法中引入变异操作，即对某些变量以一定的概率重新初始化。变异操作拓展了在迭代中不断缩小的种群搜索空间，使粒子能够跳出先前搜索到的最优值位置，在更大的空间中开展搜索，同时保持了种群多样性，提高算法寻找更优值的可能性。因此，在普通粒子群算法的基础上引入简单变异算子，在粒子每次更新之后，以一定概率重新初始化粒子。

本案例群体的粒子数为50，每个粒子的维数为2，算法迭代进化次数为1000。加速度因子$c_1=1.4$，$c_2=1.5$，惯性权重$\omega =0.8$。

适应度函数代码如下：

def funShubert(x):h1 = 0h2 = 0for i in range(1, 6):h1 += i * np.cos((i + 1) * x[0] + i)h2 += i * np.cos((i + 1) * x[1] + i)return h1 * h2

自适应变异PSO算法代码如下：

c1, c2 = 1.4, 1.5  # 加速度因子
maxgen = 1000  # 进化次数
sizepop = 50  # 群体规模
w = 0.8  # 惯性权重
Vmax, Vmin = 5, -5  # 速度最大值，最小值
popmax, popmin = 10, -10  # 个体最大值，最小值
pop = np.zeros([sizepop, 2])  # 种群
V = np.zeros([sizepop, 2])  # 速度
fitness = np.zeros(sizepop)  # 适应度
trace = np.zeros(maxgen)  # 结果# 随机产生一个群体,初始粒子和速度
for i in range(sizepop):pop[i] = 10 * (-2 * np.random.rand(1, 2) + 1)V[i] = 5 * (-2 * np.random.rand(1, 2) + 1)fitness[i] = funShubert(pop[i])# 个体极值和群体极值
bestfitness = np.min(fitness)
bestindex = np.argmin(fitness)
Gbest = pop[bestindex]  # 全局最佳
fitnessGbest = bestfitness  # 全局最佳适应度值
Pbest = pop.copy()  # 个体最佳
fitnessPbest = fitness.copy()  # 个体最佳适应度值# 迭代寻优
for i in range(maxgen):for j in range(sizepop):V[j] = w * V[j] + c1 * (-2 * np.random.rand(1, 2) + 1) * (Pbest[j] - pop[j]) \+ c2 * (-2 * np.random.rand(1, 2) + 1) * (Gbest - pop[j])  # 速度更新V[j, V[j] > Vmax] = VmaxV[j, V[j] < Vmin] = Vminpop[j] = pop[j] + V[j]  # 群体更新pop[j, pop[j] > popmax] = popmaxpop[j, pop[j] < popmin] = popminif np.random.rand() > 0.9: # 自适应变异pop[j] = 10 * (-2 * np.random.rand(1, 2) + 1)fitness[j] = funShubert(pop[j])  # 适应度值if fitness[j] < fitnessPbest[j]:Pbest[j] = pop[j].copy()fitnessPbest[j] = fitness[j]if fitness[j] < fitnessGbest:Gbest = pop[j].copy()fitnessGbest = fitness[j]trace[i] = fitnessGbestprint(Gbest, fitnessGbest)# 结果分析
plt.plot(trace)
plt.title("最优个体适应度")
plt.xlabel("进化代数")
plt.ylabel("适应度")
plt.show()

算法结果如下：

最终得到的个体适应度值为-186.6831831611777，对应的例子位置为$(5.4788292,-1.4228858)$，自适应变异PSO算法寻优得到最优值接近函数实际最优值，说明该算法具有较强的函数极值寻优能力。

五、补充

惯性权重$\omega$体现的是粒子当前速度多大程度上继承先前的速度，Shi.Y最先将惯性权重$\omega$引入到PSO算法中，并分析指出一个较大的惯性权值有利于全局搜索，而一个较小的惯性权值则更有利于局部搜索。为了更好地平衡算法的全局搜索与局部搜索能力，其提出了线性递减惯性权重（Linear Decreasing Inertia Weight，LDIW），即

$$\omega (t)={\omega }_{start}-({\omega }_{start}-{\omega }_{end})\times \frac{t}{T_{max}}$$

式中，${\omega }_{start}$为初始惯性权重，${\omega }_{end}$为迭代至最大次数时的惯性权重，$t$为当前迭代次数，$T_{m}ax$为最大迭代次数。

一般来说，惯性权重取值为${\omega }_{start}=0.9$，${\omega }_{end}=0.4$时算法性能最好。这样，随着迭代的进行，惯性权重由0.9线性递减至0.4，迭代初期较大的惯性权重使算法保持了较强的全局搜索能力，而迭代后期较小的惯性权值有利于算法进行更精确的局部搜索。

线性惯性权重只是一种经验做法，常用的惯性权重的选择还包括以下几种：

$$\omega (t)={\omega }_{start}-({\omega }_{start}-{\omega }_{end})\times (\frac{t}{T_{max}}{)}^2$$

$$\omega (t)={\omega }_{start}-({\omega }_{start}-{\omega }_{end})\times (\frac{2t}{T_{max}}-(\frac{t}{T_{max}}{)}^2)$$

六、练习

求测试函数的最小值，以及最小值点。

1. Rastrigin function：

$$f(\mathbf{x})=A_n+\sum _{i=1}^n[x_i^2-A\cos (2\pi x_i)],-5.12\le x_i\le 5.12$$

当$A=10,n=2$时，如下图所示：

2. Sphere function：

$$f(\mathbf{x})=\sum _{i=1}^n{x}_i^2,-\infty \le x_i\le \infty$$

当$n=2$时，如下图所示：

3. Beale function：

$$f(x,y)=(1.5-x+xy{)}^2+(2.25-x+x{y}^2{)}^2+(2.625-x+x{y}^3{)}^2,-4.5\le x,y\le 4.5$$

4. Booth function：

$$f(x,y)=(x+2y-7{)}^2+(2x+y-5{)}^2,-10\le x,y\le 10$$

5. Bukin function：

$$f(x,y)=100\sqrt{|y-0.01x^2|}+0.01|x+10|,-15\le x\le -5,-3\le y\le 3$$

6. three-hump camel function：

$$f(x,y)=2x^2-1.05x^4+\frac{x^6}{6}+xy+y^2,-5\le x,y\le 5$$

7. Hölder table function：

$$f(x,y)=-|\sin x\cos y\exp (|1-\frac{\sqrt{x^2+y^2}}{\pi }|)|,-10\le x,y\le 10$$