算法简介
粒子群算法（PSO）来源于对一个简化社会模型的模拟，该算法最初是为了图形化地模拟鸟群优美而不可预测的运动。在鸟群觅食过程中，我们假定鸟群中信息共享，每个个体都知道自己与其它个体的当前觅食地点和历史最佳觅食地点，因此个体下一运动状态趋势如下：
- 受飞行惯性影响，个体将沿原飞行方向运动
- 个体向自己的历史最佳觅食地点运动
- 个体向当前鸟群中最佳觅食地点运动
将这三个运动矢量合成，即可得知下一觅食地点。
这是一个迭代过程，我们认为鸟群中个体的运动是同时发生的，并将一次运动过程视为一代，当群体的规模较大且迭代次数较多时，鸟群中所有个体将向该有限空间内的最佳觅食地点靠拢。
在处理优化问题时，粒子群算法将自变量集合视为个体觅食地点的坐标，将优化目标函数视为觅食地点的评价方法，将自变量在其领域的偏移视为鸟群个体的运动过程，当鸟群聚集于最佳觅食点时，此时自变量的取值刚好满足优化目标的全局最优。
进化方程
粒子群算法将每个个体看作$D$维空间中的一个没有体积的微粒，在搜索空间中以一定速度飞行，这个速度根据它本身的飞行经验和同伴的飞行经验来动态调整。第$i$个微粒表示为$X_i=(x_{i1},x_{i2},…,x_{iD})$，它经历过的最好位置（有最好适应值）记为$P_i=(p_{i1},p_{i2},…,p_{iD})$.在群体所有微粒经过的最好位置的索引号用符号$g$表示，即$P_g$. 微粒$i$的速度为$V_i=(v_{i1},v_{i2},…,v_{iD})$，对于每一维度的第$t$代，它的第$t+1$代根据如下方程进行变化：
参数含义如下：
- $V_i^k$：第$i$个粒子第$k$代的速度，即第$k-1$代到第$k$代的运动矢量
- $X_i^k$：第$i$个粒子第$k$代的位置坐标
- $w$：自身权重系数，取$0.9$到$1.2$较合适
- $c_1$：自身认知常数，取$2.0$较合适
- $c_2$：社会认知常数，取$2.0$较合适
- $r_1,r_2$：区间$[0,1]$内随机数
- $P_i^k$：第$i$个粒子运动到第$k$代时的历史最好位置
- $P_{gi}^k$：所有微粒运动到第$k$代时的历史最好位置
除了以上在公式中出现的参数，还有其它参数需要注意，如种群规模即粒子数量，粒子位置的范围即自变量定义域，运动速度上下限等。
算法步骤
$step\ 1$ 初始化微粒（群体规模为$m$），微粒的初始位置和速度取随机值
$step\ 2$ 评价每个微粒的适应度（优化目标函数值）
$step\ 3$ 对每个微粒，将它的适应值和它经历过的最好位置pbest的作比较，如果较好，则将其作为该粒子历史最佳位置pbest
$step\ 4$ 对每个微粒，将它的适应值和全局所经历最好位置 gbest的作比较，如果较好，则将其作为全局最佳位置gbest
$step\ 5$ 依据进化方程更新微粒的速度与位置
$step\ 6$ 如未达到结束条件（通常为足够好的适应值或达到一个预设最大代数 $G_{max}$ ），回到$step\ 2$
算法案例
无约束条件
以Bukin函数为例，其函数形式如下：
已知该函数在$x^{*}=(-10,1)$时有最小值$0$。下面用粒子群算法求解最优值
该例中，粒子运动空间的维度为$2$，取种群规模为$100$，更新代数为$500$，速度范围为$[-1,1]，$自身权重因子$0.9$，自身认知常数$2.0$，社会认知常数$2.0$。程序如下：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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61import numpy as np
import random
import math
class PSO:
def __init__(self, size, nVar, xMin, xMax, vMin, vMax, times):
self.size = size # 种群规模
self.nVar = nVar # 空间维数
self.times = times # 迭代次数
self.x = np.zeros((size, nVar), dtype=float) # 粒子坐标
self.v = np.zeros((size, nVar), dtype=float) # 粒子速度
self.xMin = np.array(xMin, dtype = float) # 粒子位置下限
self.xMax = np.array(xMax, dtype = float) # 粒子位置上限
self.vMin = vMin # 粒子速度下限
self.vMax = vMax # 粒子速度上限
# 粒子坐标与速度随机值
for i in range(size):
for j in range(nVar):
self.x[i][j] = random.uniform(self.xMin[j], self.xMax[j])
self.v[i][j] = random.uniform(self.vMin, self.vMax)
# 单个粒子最有值与全局最优值索引
self.pBest = self.x.copy()
self.gBest = 0
def fitness(self, x):
x1 = x[0]
x2 = x[1]
f = 100 * np.sqrt(np.abs(x2 - 0.01*x1**2)) + 0.01 * np.abs(x1 + 10)
return f
def update(self, w, c1, c2):
for t in range(self.times):
for i in range (self.size):
# 更新粒子速度
self.v[i] = w * self.v[i] + random.uniform(0, 1) * c1 * (self.pBest[i] - self.x[i]) + random.uniform(0, 1) * c2 * (self.pBest[self.gBest] - self.x[i])
# 确保速度不超过上下限
self.v[i][self.v[i] < self.vMin] = self.vMin
self.v[i][self.v[i] > self.vMax] = self.vMax
# 更新粒子位置
self.x[i] += self.v[i]
# 确保位置不超过上下限
for j in range(self.nVar):
if self.x[i][j] > self.xMax[j]:
self.x[i][j] = self.xMax[j]
if self.x[i][j] < self.xMin[j]:
self.x[i][j] = self.xMin[j]
# 比较每个粒子的新的适应度与历史最佳适应度
if self.fitness(self.x[i]) < self.fitness(self.pBest[i]):
self.pBest[i] = self.x[i]
# 比较当前粒子的新适应度与全局最佳适应度
if self.fitness(self.x[i]) < self.fitness(self.pBest[self.gBest]):
self.gBest = i
def display(self):
print(self.pBest[self.gBest])
print(self.fitness(self.pBest[self.gBest]))
if __name__ == '__main__':
# 种群规模为100，更新500代
bukin = PSO(size=100, nVar=2, xMin=[-15, -3], xMax=[-5, 3], vMin=-1, vMax=1, times=500)
# 自身权重因子0.9，自身认知常数2.0，社会认知常数2.0
bukin.update(w = 0.9, c1 = 2.0, c2 = 2.0)
bukin.display()
筛选最优输出结果为1
2[-10.77631364 1.16128901]
0.0662286658094341
有约束条件
与模拟退火算法相同，在处理约束条件时我们依旧采用罚函数的思想，将约束条件转换为限制粒子飞行的因素。
对于该例子，利用惩罚函数法，将约束条件转化为惩罚项，定义$P_1=(\min(x_1+0.5x_2-0.4,0))^2$，$P_2=(\min(0.5x_1+x_2-0.5,0))^2$，此时构造新的增广目标函数$min\ F(X,M_t)=f(X)+M_t(P_1+P_2)$，其中$M_t$为惩罚因子。
粒子运动空间的维度为$2$，取种群规模为$100$，更新代数为$500$，速度范围为$[-1,1]，$自身权重因子$0.9$，自身认知常数$2.0$，社会认知常数$2.0$。该例自变量没有上限，但我们可以粗略估计在最优解时自变量$x_1$，$x_2$都落在$[0,1.0]$内，因此将粒子坐标的上限设为$1.0$。
1 | import numpy as np |
输出结果为1
2[0.34072776 0.32958779]
0.2455774890528109
1 | def update2d(self, w, c1, c2): |
1 | [-10.39337108 1.08022159] |
在全局算法中，一般地，希望算法前期有较高的搜索能力以得到合适的种子，而在后期有较高的开发能力以加快收敛速度。也就是希望算法前期全局搜索能力较强，后期剧不搜索能力较强。因此，PSO算法做如下更改：
$step\ 4$ 依据二阶振荡PSO进化方程更新微粒的速度与位置。当前代数$t<G_{max}/2$时，取$0<\xi_2<\frac{1+\xi_1}{2}$，$0<\xi_4<\frac{1+\xi_3}{2}$，$0<\xi_1<1$，$0<\xi_3<1$；当前代数$t\geq G_{max}/2$时，取$\xi_1<\xi_2-1$，$\xi_3<\xi_4-1$，$0<\xi_2<1$，$0<\xi_4<1$.
二阶振荡粒子群算法程序如下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
38def update2dvib(self, w, c1, c2):
for t in range(self.times):
# 上一代粒子的位置，初始为所有粒子的初始位置
xPrev = self.x.copy()
for i in range (self.size):
# 设置参数
phi1 = random.uniform(0, 1) * c1
phi2 = random.uniform(0, 1) * c2
# 振荡因子在前半段和后半段取值不同
if t < self.times/2:
xi1 = random.uniform(0, 1)
xi3 = random.uniform(0, 1)
xi2 = random.uniform(0, (1 + xi1)/2)
xi4 = random.uniform(0, (1 + xi3)/2)
else:
xi2 = random.uniform(0, 1)
xi4 = random.uniform(0, 1)
xi1 = random.uniform(0, xi2 - 1)
xi3 = random.uniform(0, xi4 - 1)
# 二阶振荡粒子群算法进化方程
self.v[i] = w * self.v[i] + phi1 * (self.pBest[i] - (1+xi1)*self.x[i] + xi2 * xPrev[i]) + phi2 * (self.pBest[self.gBest] - (1+xi3)*self.x[i] + xi4 * xPrev[i])
# 确保速度不超过上下限
self.v[i][self.v[i] < self.vMin] = self.vMin
self.v[i][self.v[i] > self.vMax] = self.vMax
# 更新粒子位置
self.x[i] += self.v[i]
# 确保位置不超过上下限
for j in range(self.nVar):
if self.x[i][j] > self.xMax[j]:
self.x[i][j] = self.xMax[j]
if self.x[i][j] < self.xMin[j]:
self.x[i][j] = self.xMin[j]
# 比较每个粒子的新的适应度与历史最佳适应度
if self.fitness(self.x[i]) < self.fitness(self.pBest[i]):
self.pBest[i] = self.x[i]
# 比较当前粒子的新适应度与全局最佳适应度
if self.fitness(self.x[i]) < self.fitness(self.pBest[self.gBest]):
self.gBest = i
依旧以Bukin函数为例，筛选出的最优结果为1
2[-10. 1.]
0.0
可以看到结果非常接近标准最优解，这说明二阶振荡PSO算法相较于二阶PSO算法与标准PSO算法精度得到很大提升。