简介
1995年，美国学者Kennedy和Eberhart共同提出了粒子群算法，其基本思想源于对鸟类群体行为进行建模与仿真的研究结果的启发。
粒子群算法，其全称为粒子群优化算法(Particle Swarm Optimization,PSO) 。它是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的搜索算法。
符号说明
| 符号 | 含义 |
|---|---|
| n n n | 粒子个数 |
| c 1 c_{1} c1 | 粒子的个体学习因子，也称个体加速因子 |
| c 2 c_{2} c2 | 粒子的社会学习因子，也称社会加速因子 |
| w w w | 速度的惯性权重 |
| v i d v^{d}_{i} vid | 第d次迭代时，第i个粒子的速度 |
| x i d x^{d}_{i} xid | 第d次迭代时，第i个粒子所在位置 |
| f ( x ) f(x) f(x) | 在位置x时的适应度值（一般为目标函数值） |
| p b e s t i d pbest^{d}_{i} pbestid | 到第d次迭代为止，第i个粒子经过的最好的位置 |
| g b e s t d gbest^{d} gbestd | 到第d次迭代为止，所有粒子经过的最好的位置 |
核心思想
粒子群算法是通过模拟鸟群觅食行为一种基于群体协作的搜索算法。
有一群鸟需要觅食，我们对这群鸟作出几个假设：
（1）所有鸟都不知道食物最丰盛的地方在哪；
（2）所有鸟同时共享当前已知最为丰盛的食物所在地方；
（3）每只鸟有一定的“判断能力”。对于其中的某一只鸟它的“选择”（并非主观选择，由随机数控制）有三种：
- 独狼。不做出任何调整，沿着本来的方向飞行一步
- 随群。做出一个从众的调整，朝向目前最优的位置飞行一步
- 自我。做出一个盲目自信的调整，朝向自己去过的最优的位置飞行一步
第i只鸟（粒子）第d步的速度 = 上一步自身的惯性速度 + 自我认知部分 + 社会认知部分(r1、r2为[0, 1]上的随机数)： v i d = w v i d − 1 + c 1 ∗ r 1 ∗ ( p b e s t i d − x i d ) + c 2 ∗ r 2 ∗ ( g b e s t d − x i d ) v^{d}_{i} = wv^{d-1}_{i}+c_{1}*r_{1}*(pbest^{d}_{i} - x^{d}_{i}) + c_{2}*r_{2}*(gbest^{d} - x^{d}_{i}) vid=wvid−1+c1∗r1∗(pbestid−xid)+c2∗r2∗(gbestd−xid)第i只鸟（粒子）第d + 1步的位置 = 第d步所在的位置 + 第d步的速度 * 运动时间（运动时间一般取1）： x i d + 1 = x i d + v i d x^{d + 1}_{i} = x^{d}_{i} + v^{d}_{i} xid+1=xid+vid
让我们来对速度转移公式进行分析： w v i d − 1 wv^{d-1}_{i} wvid−1：每一个粒子都有自身的惯性速度，由惯性权重控制取值； c 1 ∗ r 1 ∗ ( p b e s t i d − x i d ) c_{1}*r_{1}*(pbest^{d}_{i} - x^{d}_{i}) c1∗r1∗(pbestid−xid)：自我认知部分，每个粒子渴望去往曾经自己去过最优的地方，由随机数和个体学习因子控制取值； c 2 ∗ r 2 ∗ ( g b e s t d − x i d ) c_{2}*r_{2}*(gbest^{d} - x^{d}_{i}) c2∗r2∗(gbestd−xid)：社会认知部分，每个粒子有去往所有粒子去过最优的地方，由随机数和社会学习因子控制取值。
可以直观得出几种运动状态：
（1）某个粒子距离全局最优的部分越近则社会认知部分取值越小，反之越大；
（2）某个粒子距离自身去过最优的部分越近则自我认知部分取值越小，反之越大；
（3）上一步的选择，会通过惯性速度保留影响。惯性速度越大，则粒子越活跃，全局搜索能力越强。
如此，粒子在反复调整下，搜索并寻找最优值。
流程图
文章使用到的测试函数
一元函数：
y
=
−
11
s
i
n
(
x
)
−
7
c
o
s
(
5
x
)
,
x
∈
[
−
3
,
3
]
y = -11sin(x)-7cos(5x), \ \ \ x \in [-3, 3]
y=−11sin(x)−7cos(5x), x∈[−3,3]参考最小值为：
m
i
n
(
y
)
=
−
17.4905
min(y) = -17.4905
min(y)=−17.4905
二元函数：
y
=
x
1
2
+
x
2
2
−
x
1
x
2
−
10
x
1
−
4
x
2
+
60
y = x_{1}^{2}+x_{2}^{2}-x_{1}x_{2}-10x_{1}-4x_{2}+60
y=x12+x22−x1x2−10x1−4x2+60参考最小值为：
m
i
n
(
y
)
=
8.0
min(y) = 8.0
min(y)=8.0
粒子群算法代码
Python版本：
from random import random
from copy import deepcopy
import math
def func(x):
# 一元函数测试
return -11 * math.sin(x[0]) - 7 * math.cos(5 * x[0])
# # 二元函数测试
# return x[0]**2 + x[1]**2 - x[0]*x[1] - 10 * x[0] - 4 * x[1] + 60
class PSO:
def __init__(self, func, n, var = 1, w = 0.9, iter = 50, lb = None, ub = None):
""" 默认寻找最小值，以及对应自变量取值。 若需要寻找最大值，对目标函数乘-1，并对最终结果乘-1即可
Args:
:param func: type: 函数， des: 所要求解优化问题的目标函数
:param n: type: int， des: 粒子群粒子的个数
:param var: type: int， des: 函数中自变量的个数，即：几元函数
:param lb: type: list(double)， des: 每一种自变量的下界，注意应和自变量一一对应
:param ub: type: list(double)， des: 每一种自变量的上界，注意应和自变量一一对应
:param iter: type: int， des: 迭代次数
"""
if lb is None:
lb = []
if ub is None:
ub = []
# 目标函数
self.func = func
# 粒子个数
self.n = n
# 迭代次数
self.iter = iter
# 函数中自变量个数
self.var = var
# 惯性权重
self.w = w
# 初始化每个自变量的范围 len(lb) = len(ub) = var
self.lb = lb
self.ub = ub
# 存储每一个粒子的位置
## 二维数组(n * m)：n表示n个粒子， m表示每一个粒子第m个自变量位置
self.x = [[0] * var for _ in range(n)]
## 初始化粒子位置
for i in range(n):
for j in range(var):
self.x[i][j] = lb[j] + random() * (ub[j] - lb[j])
# 存储每一个粒子每一个自变量的速度
self.v = [[0] * var for _ in range(n)]
## 粒子的最大速度通常为可行域的10% - 20%
## 控制所有速度在-max_v到max_v之间
for i in range(n):
for j in range(var):
max_v = (ub[j] - lb[j]) * 0.2
self.v[i][j] = -max_v + 2 * max_v * random()
# 全局最优解 pbest：最优解位置
# 先赋一个初始值， 待会儿更新
self.pbest = deepcopy(self.x[0])
# 每个粒子的最优解位置
# gbest: 每个粒子最优解位置
self.gbest = [[0] * var for _ in range(n)]
# 存储每一个粒子的适应度
self.fit = [0] * n
## 使用初始值计算函数值（适应度）
for i in range(n):
self.fit[i] = self.func(self.x[i])
# 对pbest赋初值
if self.fit[i] < self.func(self.pbest):
self.pbest = deepcopy(self.x[i])
# 以第一组计算值作为，每个粒子目前的最优值
self.gbest[i] = deepcopy(self.x[i])
def run(self):
# 迭代次数
for k in range(self.iter):
# n个粒子
for i in range(self.n):
# j个因变量
for j in range(self.var):
# 更新速度 这里c1 = c2 = 2
self.v[i][j] = self.w * self.v[i][j] + 2 * random() * (self.pbest[j] - self.x[i][j]) \
+ 2 * random() * (self.gbest[i][j] - self.x[i][j])
# 限制速度范围
max_v = (self.ub[j] - self.lb[j]) * 0.2
if self.v[i][j] < -max_v:
self.v[i][j] = -max_v
if self.v[i][j] > max_v:
self.v[i][j] = max_v
# 更新粒子位置
self.x[i][j] += self.v[i][j]
# 调整位置避免跑出范围
if self.x[i][j] < self.lb[j]:
self.x[i][j] = self.lb[j]
if self.x[i][j] > self.ub[j]:
self.x[i][j] = self.ub[j]
# 计算函数值
self.fit[i] = self.func(self.x[i])
# 更新全局最优位置
if self.fit[i] < self.func(self.pbest):
# 最优位置
self.pbest = deepcopy(self.x[i])
# 更新每个粒子的最优位置
if self.fit[i] < self.func(self.gbest[i]):
# 第i个粒子的最优位置
self.gbest[i] = deepcopy(self.x[i])
print(f'函数最小值为{self.func(self.pbest)}')
print(f'最小值对应的一组x取值为{self.pbest}')
"""一元函数测试"""
pso = PSO(func, 50, 1, 0.9, 50, [-3], [3])
"""二元函数测试"""
# pso = PSO(func, 50, 2, 0.9, 50, [-15, -15], [15, 15])
pso.run()
粒子群算法优化
速度转移公式： v i d = w v i d − 1 + c 1 ∗ r 1 ∗ ( p b e s t i d − x i d ) + c 2 ∗ r 2 ∗ ( g b e s t d − x i d ) v^{d}_{i} = wv^{d-1}_{i}+c_{1}*r_{1}*(pbest^{d}_{i} - x^{d}_{i}) + c_{2}*r_{2}*(gbest^{d} - x^{d}_{i}) vid=wvid−1+c1∗r1∗(pbestid−xid)+c2∗r2∗(gbestd−xid)
基于惯性权重的优化
惯性权重 w w w体现的是粒子继承先前的速度的能力，Shi,Y最先将惯性权重 w w w引入到粒子群算法中，并分析指出：一个较大的惯性权值有利于全局搜索，而一个较小的权值则更利于局部搜索。
递减惯性权重优化
为了更好地平衡算法的全局搜索以及局部搜索能力，Shi,Y提出了线性递减惯性权重LDIW(Linear Decreasing Inertia Weight),公式如下： w d = w s t a r t − ( w s t a r t − w e n d ) × d K w^{d} = w_{start} - (w_{start}-w_{end})\times\frac{d}{K} wd=wstart−(wstart−wend)×Kd其中： d d d为当前迭代次数， K K K为总迭代次数。 w s t a r t w_{start} wstart一般取0.9， w e n d w_{end} wend一般取0.4。
与原来相比，现在的惯性权重和迭代次数有关。早期权重较大，有利于全局搜索最优解；后期权重较小，有利于精确锁定最优解。
另外两种递减惯性优化：
w
d
=
w
s
t
a
r
t
−
(
w
s
t
a
r
t
−
w
e
n
d
)
×
(
d
K
)
2
w^{d} = w_{start} - (w_{start}-w_{end})\times (\frac{d}{K})^2
wd=wstart−(wstart−wend)×(Kd)2和
w
d
=
w
s
t
a
r
t
−
(
w
s
t
a
r
t
−
w
e
n
d
)
×
[
2
d
K
−
(
d
K
)
2
]
w^{d} = w_{start} - (w_{start}-w_{end})\times [\frac{2d}{K} - (\frac{d}{K})^2 ]
wd=wstart−(wstart−wend)×[K2d−(Kd)2]按照次序依次作出三种优化方式，权重
w
w
w随迭代次数的取值：
如图可以看出，第一种优化方式
w
w
w随迭代次数线性递减，即随着迭代次数增加逐步偏重局部搜索；第二种优化方式
w
w
w前期缓慢下降后期快速下降，即前期偏重全局搜索，后期偏重局部搜索；第三种优化方式基本上与第二种相反。
综合比较下，推荐使用第二种优化方式：
只需要将PSO类中速度转移部分的w进行修改即可：
# 更新速度 此处默认c1 = c2 = 2
# k + 1 由于k从0开始
w = 0.9 - 0.5 * ((k + 1) / self.iter) ** 2
self.v[i][j] = w * self.v[i][j] + 2 * random() * (self.pbest[j] - self.x[i][j]) \
+ 2 * random() * (self.gbest[i][j] - self.x[i][j])
输出比较（使用二元函数测试）：
# 未优化前
函数最小值为8.001278684912592
最小值对应的一组x取值为[8.028295198814478, 6.040190279868983]
# 使用方法二优化后
函数最小值为8.00000035536938
最小值对应的一组x取值为[7.99935395039432, 5.999882728234289]
效果非常显著（非常自然的测试，并没有刻意调整结果）！精度从 1 e − 3 1e^{-3} 1e−3进一步精确到了 1 e − 7 1e^{-7} 1e−7。
这里给出所修改部分的代码（只修改了PSO类中run函数）：
def run(self):
# 迭代次数
for k in range(self.iter):
# n个粒子
for i in range(self.n):
# j个因变量
for j in range(self.var):
# 更新速度 此处默认c1 = c2 = 2
# k + 1 由于k从0开始
w = 0.9 - 0.5 * ((k + 1) / self.iter) ** 2
self.v[i][j] = w * self.v[i][j] + 2 * random() * (self.pbest[j] - self.x[i][j]) \
+ 2 * random() * (self.gbest[i][j] - self.x[i][j])
# 限制速度范围
max_v = (self.ub[j] - self.lb[j]) * 0.2
if self.v[i][j] < -max_v:
self.v[i][j] = -max_v
if self.v[i][j] > max_v:
self.v[i][j] = max_v
# 更新粒子位置
self.x[i][j] += self.v[i][j]
# 调整位置避免跑出范围
if self.x[i][j] < self.lb[j]:
self.x[i][j] = self.lb[j]
if self.x[i][j] > self.ub[j]:
self.x[i][j] = self.ub[j]
# 计算函数值
self.fit[i] = self.func(self.x[i])
# 更新全局最优位置
if self.fit[i] < self.func(self.pbest):
# 最优位置
self.pbest = deepcopy(self.x[i])
# 更新每个粒子的最优位置
if self.fit[i] < self.func(self.gbest[i]):
# 第i个粒子的最优位置
self.gbest[i] = deepcopy(self.x[i])
print(f'函数最小值为{self.func(self.pbest)}')
print(f'最小值对应的一组x取值为{self.pbest}')
自适应惯性权重优化
我们可以根据当前位置计算出的适应度（函数值）来调整对应的权重取值，即通过当前的适应度判断需要偏重全局还是偏重局部搜索。
假设当前在求最小值问题，那么：
w
i
d
=
{
w
m
i
n
+
(
w
m
a
x
−
w
m
i
n
)
f
(
x
i
d
)
−
f
m
i
n
d
f
a
v
e
r
a
g
e
d
−
f
m
i
n
d
,
f
(
x
i
d
)
≤
f
a
v
e
r
a
g
e
d
w
m
a
x
,
f
(
x
i
d
)
>
f
a
v
e
r
a
g
e
d
w^{d}_{i} = \left\{\begin{matrix} w_{min} + (w_{max} - w_{min})\frac{f(x_{i}^{d})-f^{d}_{min}}{ f^{d}_{average} - f^{d}_{min} }, \ \ \ f(x_{i}^{d}) \le f^{d}_{average} \\ \\ w_{max}, \ \ \ \ \ \ f(x_{i}^{d}) > f^{d}_{average} \end{matrix}\right.
wid=⎩
⎨
⎧wmin+(wmax−wmin)faveraged−fmindf(xid)−fmind, f(xid)≤faveragedwmax, f(xid)>faveraged其中：
（1）
w
m
i
n
w_{min}
wmin和
w
m
a
x
w_{max}
wmax是我们预先给定的最小惯性系数和最大惯性系数，一般取0.4和0.9；
（2） f a v e r a g e d = ∑ i = 1 n f ( x i d ) / n f^{d}_{average} = \sum_{i=1}^{n} f(x^{d}_{i})/n faveraged=∑i=1nf(xid)/n，即第d次迭代时所有粒子的平均适应度；
（3） f m i n d = m a x ( f ( x 1 d ) , f ( x 2 d ) , . . . , f ( x n d ) ) f^{d}_{min} =max(f(x^{d}_{1}),f(x^{d}_{2}),...,f(x^{d}_{n})) fmind=max(f(x1d),f(x2d),...,f(xnd))，即第d次迭代时所有粒子的最小适应度。
适应度越小，说明距离最优解越近，此时更需要局部搜索；适应度越大，说明距离最优解越远，此时更需要全局搜索。
输出比较（使用二元函数测试）：
# 未优化前
函数最小值为8.001278684912592
最小值对应的一组x取值为[8.028295198814478, 6.040190279868983]
# 使用自适应惯性权重优化
函数最小值为8.00000000000297
最小值对应的一组x取值为[8.00000177487541, 6.000001662692137]
效果异常显著（非常自然的测试，并没有刻意调整结果）！精度从 1 e − 3 1e^{-3} 1e−3进一步精确到了 1 e − 12 1e^{-12} 1e−12！
这里给出所修改部分的代码（只修改了PSO类中run函数）：
def run(self):
# 迭代次数
for k in range(self.iter):
# 计算n个粒子的 平均适应度 和 最小适应度
faverage = 0
fmin = self.fit[0]
for i in range(self.n):
# 计算平均适应度
faverage += self.fit[i]
# 计算最小适应度
fmin = min(fmin, self.fit[i])
# 计算平均适应度
faverage /= self.n
# n个粒子
for i in range(self.n):
# j个因变量
for j in range(self.var):
# 更新速度 此处默认c1 = c2 = 2
# 计算自适应权重
w = 0
# wmax = 0.9 wmin = 0.4
if self.fit[i] > faverage: # 大于平均值说明更需要全局搜索
w = 0.9
else: # 否则，需要局部搜索，进一步定位最优值
w = 0.4 + 0.5 * (self.fit[i] - fmin) / (faverage - fmin)
self.v[i][j] = w * self.v[i][j] + 2 * random() * (self.pbest[j] - self.x[i][j]) \
+ 2 * random() * (self.gbest[i][j] - self.x[i][j])
# 限制速度范围
max_v = (self.ub[j] - self.lb[j]) * 0.2
if self.v[i][j] < -max_v:
self.v[i][j] = -max_v
if self.v[i][j] > max_v:
self.v[i][j] = max_v
# 更新粒子位置
self.x[i][j] += self.v[i][j]
# 调整位置避免跑出范围
if self.x[i][j] < self.lb[j]:
self.x[i][j] = self.lb[j]
if self.x[i][j] > self.ub[j]:
self.x[i][j] = self.ub[j]
# 计算函数值
self.fit[i] = self.func(self.x[i])
# 更新全局最优位置
if self.fit[i] < self.func(self.pbest):
# 最优位置
self.pbest = deepcopy(self.x[i])
# 更新每个粒子的最优位置
if self.fit[i] < self.func(self.gbest[i]):
# 第i个粒子的最优位置
self.gbest[i] = deepcopy(self.x[i])
print(f'函数最小值为{self.func(self.pbest)}')
print(f'最小值对应的一组x取值为{self.pbest}')
基于学习因子的优化
压缩因子法
基于前人的研究成果（参考粒子群算法相关论文），在之前的计算中个体学习因子 c 1 c_{1} c1和社会学习因子 c 2 c_{2} c2默认取2。
个体学习因子 c 1 c_{1} c1和社会(群体)学习因子 c 2 c_{2} c2决定了粒子本身经验信息和其他粒子的经验信息对粒子运行轨迹的影响，其反映了粒子群之间的信息交流。
设置 c 1 c_{1} c1较大的值，会使粒子过多地在自身的局部范围内搜索，而较大的 c 2 c_{2} c2的值，则又会促使粒子过早收敛到局部最优值。为了有效地控制粒子的飞行速度，使算法达到全局搜索与局部搜索两者间的有效平衡，M.Clerc构造了引入收缩因子的PSO模型，采用了压缩因子，这种调整方法通过合适选取参数，可确保PSO算法的收敛性，并可取消对速度的边界限制。
压缩因子法中 c 1 c_{1} c1和 c 2 c_{2} c2取2.05， C = c 1 + c 2 = 4.1 C = c_{1} + c_{2} = 4.1 C=c1+c2=4.1，收缩因子： Φ = 2 ∣ 2 − C − c 2 − 4 C ∣ = 0.7298 \Phi = \frac{2}{|2 - C - \sqrt{c^2-4C} |} = 0.7298 Φ=∣2−C−c2−4C∣2=0.7298原速度转移公式更新为： v i d = Φ ∗ [ w v i d − 1 + c 1 ∗ r 1 ∗ ( p b e s t i d − x i d ) + c 2 ∗ r 2 ∗ ( g b e s t d − x i d ) ] v^{d}_{i} = \Phi * [wv^{d-1}_{i}+c_{1}*r_{1}*(pbest^{d}_{i} - x^{d}_{i}) + c_{2}*r_{2}*(gbest^{d} - x^{d}_{i})] vid=Φ∗[wvid−1+c1∗r1∗(pbestid−xid)+c2∗r2∗(gbestd−xid)]注意： Φ \Phi Φ是一个准确值0.7298，代入原方程即可。
输出比较（使用二元函数测试）：
# 未优化前
函数最小值为8.001278684912592
最小值对应的一组x取值为[8.028295198814478, 6.040190279868983]
# 使用压缩因子法优化
函数最小值为8.000000000403176
最小值对应的一组x取值为[8.000017106999763, 6.000022107058988]
效果相对比较显著（非常自然的测试，并没有刻意调整结果）！精度从 1 e − 3 1e^{-3} 1e−3进一步精确到了 1 e − 10 1e^{-10} 1e−10！
建议使用混合优化方式，惯性权重 + 压缩因子效果异常异常显著！！！
这里给出所修改部分的代码（只修改了PSO类中run函数，注意：只使用了压缩因子法，因为要切题，建议读者自行修改混合优化版本，非常容易）：
def run(self):
# 迭代次数
for k in range(self.iter):
# n个粒子
for i in range(self.n):
# j个因变量
for j in range(self.var):
# 更新速度
self.v[i][j] = self.w * self.v[i][j] + 2.05 * random() * (self.pbest[j] - self.x[i][j]) \
+ 2.05 * random() * (self.gbest[i][j] - self.x[i][j])
# 压缩因子法
self.v[i][j] *= 0.7298
# 限制速度范围
max_v = (self.ub[j] - self.lb[j]) * 0.2
if self.v[i][j] < -max_v:
self.v[i][j] = -max_v
if self.v[i][j] > max_v:
self.v[i][j] = max_v
# 更新粒子位置
self.x[i][j] += self.v[i][j]
# 调整位置避免跑出范围
if self.x[i][j] < self.lb[j]:
self.x[i][j] = self.lb[j]
if self.x[i][j] > self.ub[j]:
self.x[i][j] = self.ub[j]
# 计算函数值
self.fit[i] = self.func(self.x[i])
# 更新全局最优位置
if self.fit[i] < self.func(self.pbest):
# 最优位置
self.pbest = deepcopy(self.x[i])
# 更新每个粒子的最优位置
if self.fit[i] < self.func(self.gbest[i]):
# 第i个粒子的最优位置
self.gbest[i] = deepcopy(self.x[i])
print(f'函数最小值为{self.func(self.pbest)}')
print(f'最小值对应的一组x取值为{self.pbest}')
扫一扫
1905
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
