前两篇链接：
基本的速度更新公式为：
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_i^d=wv_i^{d-1}+c_1r_1(pbest_i^d-x_i^d)+c_2r_2(gbest^d-x_i^d)
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
)
其中，下标
i
i
i
表示第i个粒子，上标
d
d
d
表示第d次循环，
v
i
d
v_i^d
v
i
d
表示第i个粒子第d次循环的速度，
w
w
w
为惯性权重，
c
1
,
c
2
c_1,c_2
c
1
,
c
2
分别表示个体学习因子和社会学习因子，一般取2，
p
b
e
s
t
i
d
pbest_i^d
p
b
e
s
t
i
d
表示第i个粒子在前d-1次循环中出现的最佳位置，
g
b
e
s
t
d
gbest^d
g
b
e
s
t
d
表示所有粒子在前d-1次循环中出现的最佳位置，
r
1
,
r
2
r_1,r_2
r
1
,
r
2
分别是两个位于
[
0
,
1
]
[0,1]
[
0
,
1
]
的随机数。
随机惯性权重公式1
-
题目
求函数
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_1x_2-10x_1-4x_2+60
y
=
x
1
2
+
x
2
2
−
x
1
x
2
−
1
0
x
1
−
4
x
2
+
6
0
在
x1
,
x
2
∈
[
−
15
,
15
]
x_1,x_2∈[-15,15]
x
1
,
x
2
∈
[
−
1
5
,
1
5
]
内的最小值。 -
流程图
-
数学公式：
vi
d
=
r
0
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_i^d=r_0v_i^{d-1}+c_1r_1(pbest_i^d-x_i^d)+c_2r_2(gbest^d-x_i^d)
v
i
d
=
r
0
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
)
其中
r0
r_0
r
0
是一个均匀分布在
[0
,
1
]
[0,1]
[
0
,
1
]
上的随机数，其他参数含义与上面相同。 - Python代码：
# 第一步，绘制函数图像
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
def func(x,y):
return x**2+y**2-x*y-10*x-4*y+60
x0 = np.linspace(-15,15,100)
y0 = np.linspace(-15,15,100)
x0,y0 = np.meshgrid(x0,y0)
z0 = func(x0,y0)
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x0,y0,z0,cmap=plt.cm.viridis,alpha=0.7)
#ax.plot_wireframe(x0,y0,z0) # 另一种绘图方式
ax.set_title('$y = x_1^2+x_2^2-x_1x_2-10x_1-4x_2+60$')
# 第二步，设置粒子群算法的参数
n = 30 # 粒子数量
narvs = 2 # 变量个数
c1 = 2 # 个体学习因子
c2 = 2 # 社会学习因子
w_max = 0.9 # 惯性权重
w_min = 0.4
K = 40 # 迭代次数
vxmax = np.array([(15-(-15))*0.2,(15-(-15))*0.2]) # 粒子在x方向的最大速度
x_lb = np.array([-15,-15]) # x和y的下界
x_ub = np.array([15,15]) # x和y的上界
# 第三步，初始化粒子
x = x_lb + (x_ub-x_lb)*np.random.rand(n,narvs)
v = -vxmax + 2*vxmax*np.random.rand(n,narvs)
# 第四步，计算适应度
fit = func(x[:,0],x[:,1]) # 计算每个粒子的适应度
pbest = x # 初始化这n个例子迄今为止找到的最佳位置
ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标
gbest = x[ind,:]
gbest_total = np.zeros(K)
# 第五步，更新粒子速度和位置
for j in range(K): # 外层循环，共K次
w = np.random.rand() # 随机数作为惯性权重
for p in range(n):
v[p,:] = w*v[p,:] + c1*np.random.rand(narvs)*(pbest[p,:]-x[p,:]) + c2*np.random.rand(narvs)*(gbest-x[p,:])
loc_v = np.where(v<-vxmax)
v[loc_v] = -vxmax[loc_v[1]] # 速度小于-vmax的元素赋值为-vmax
loc_v = np.where(v>vxmax)
v[loc_v] = vxmax[loc_v[1]] # 速度大于vmax的元素赋值为vmax
x = x + v # 更新第i个粒子的位置
loc_x = np.where(x<x_lb)
x[loc_x] = x_lb[loc_x[1]]
loc_x = np.where(x>x_ub)
x[loc_x] = x_ub[loc_x[1]]
# 第六步，重新计算适应度并找到最优粒子
fit = func(x[:,0],x[:,1]) # 重新计算n个粒子的适应度
for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置
if fit[k]<func(pbest[k,0],pbest[k,1]):
pbest[k,:] = x[k,:]
if np.min(fit)<func(gbest[0],gbest[1]): # 更新所有粒子迄今找到的最佳位置
gbest = x[np.argmin(fit),:]
gbest_total[j] = func(gbest[0],gbest[1])
ax.scatter(x[:,0],x[:,1],fit,c='r',marker='x')
fig,ax = plt.subplots()
ax.plot(np.arange(K),gbest_total)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来显示负号
ax.set_title('算法找到的最小值随优化次数的变化曲线')
ax.set_xlabel('循环次数')
ax.autoscale(axis='x',tight=True)
print('找到的最优解为：',func(gbest[0],gbest[1]))
- 绘制结果
-
打印输出最优解：
找到的最优解为： 8.000000066898039
可以看出，比原始算法的精度提高了很多，与非线性递减权重效果可比拟。 -
显示最终粒子所在的位置：
其中
红色叉号
表示最后一次循环结束后粒子所在的位置。可以看出，相比原始算法，粒子更多分布在最优点周围。 -
循环过程中最优解随循环次数的变化曲线：
可以看出，收敛速度明显高于原始算法的速度。
随机惯性权重公式2
-
数学公式：
w=
μ
m
i
n
+
(
μ
m
a
x
−
μ
m
i
n
)
×
r
a
n
d
(
)
+
σ
×
r
a
n
d
n
(
)
w=\mu_{min}+(\mu_{max}-\mu_{min})×rand()+\sigma×randn()
w
=
μ
m
i
n
+
(
μ
m
a
x
−
μ
m
i
n
)
×
r
a
n
d
(
)
+
σ
×
r
a
n
d
n
(
)
其中
μm
i
n
\mu_{min}
μ
m
i
n
是随机惯性权重的最小值，
μm
a
x
\mu_{max}
μ
m
a
x
是随机惯性权重的最大值，
ra
n
d
(
)
rand()
r
a
n
d
(
)
是在
[0
,
1
]
[0,1]
[
0
,
1
]
上均匀分布的随机数，
σ\sigma
σ
是标准差，一般取0.2~0.5之间的一个数，
ra
n
d
n
(
)
randn()
r
a
n
d
n
(
)
为正态分布的随机数。 - Python代码
# 第一步，绘制函数图像
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
def func(x,y):
return x**2+y**2-x*y-10*x-4*y+60
x0 = np.linspace(-15,15,100)
y0 = np.linspace(-15,15,100)
x0,y0 = np.meshgrid(x0,y0)
z0 = func(x0,y0)
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x0,y0,z0,cmap=plt.cm.viridis,alpha=0.7)
#ax.plot_wireframe(x0,y0,z0) # 另一种绘图方式
ax.set_title('$y = x_1^2+x_2^2-x_1x_2-10x_1-4x_2+60$')
# 第二步，设置粒子群算法的参数
n = 30 # 粒子数量
narvs = 2 # 变量个数
c1 = 2 # 个体学习因子
c2 = 2 # 社会学习因子
w_max = 0.9 # 惯性权重
w_min = 0.4
sigma = 0.3 # 标准差
K = 40 # 迭代次数
vxmax = np.array([(15-(-15))*0.2,(15-(-15))*0.2]) # 粒子在x方向的最大速度
x_lb = np.array([-15,-15]) # x和y的下界
x_ub = np.array([15,15]) # x和y的上界
# 第三步，初始化粒子
x = x_lb + (x_ub-x_lb)*np.random.rand(n,narvs)
v = -vxmax + 2*vxmax*np.random.rand(n,narvs)
# 第四步，计算适应度
fit = func(x[:,0],x[:,1]) # 计算每个粒子的适应度
pbest = x # 初始化这n个例子迄今为止找到的最佳位置
ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标
gbest = x[ind,:]
gbest_total = np.zeros(K)
# 第五步，更新粒子速度和位置
for j in range(K): # 外层循环，共K次
w = w_min + (w_max-w_min)*np.random.rand() + sigma*np.random.randn()
for p in range(n):
v[p,:] = w*v[p,:] + c1*np.random.rand(narvs)*(pbest[p,:]-x[p,:]) + c2*np.random.rand(narvs)*(gbest-x[p,:])
loc_v = np.where(v<-vxmax)
v[loc_v] = -vxmax[loc_v[1]] # 速度小于-vmax的元素赋值为-vmax
loc_v = np.where(v>vxmax)
v[loc_v] = vxmax[loc_v[1]] # 速度大于vmax的元素赋值为vmax
x = x + v # 更新第i个粒子的位置
loc_x = np.where(x<x_lb)
x[loc_x] = x_lb[loc_x[1]]
loc_x = np.where(x>x_ub)
x[loc_x] = x_ub[loc_x[1]]
# 第六步，重新计算适应度并找到最优粒子
fit = func(x[:,0],x[:,1]) # 重新计算n个粒子的适应度
for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置
if fit[k]<func(pbest[k,0],pbest[k,1]):
pbest[k,:] = x[k,:]
if np.min(fit)<func(gbest[0],gbest[1]): # 更新所有粒子迄今找到的最佳位置
gbest = x[np.argmin(fit),:]
gbest_total[j] = func(gbest[0],gbest[1])
ax.scatter(x[:,0],x[:,1],fit,c='r',marker='x')
fig,ax = plt.subplots()
ax.plot(np.arange(K),gbest_total)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来显示负号
ax.set_title('算法找到的最小值随优化次数的变化曲线')
ax.set_xlabel('循环次数')
ax.autoscale(axis='x',tight=True)
print('找到的最优解为：',func(gbest[0],gbest[1]))
- 绘制结果
-
打印输出最优解：
找到的最优解为： 8.000000074157867
可以看出，比原始算法的精度提高了很多，与非线性递减权重效果可比拟。 -
显示最终粒子所在的位置：
其中
红色叉号
表示最后一次循环结束后粒子所在的位置。可以看出，相比原始算法，粒子更多分布在最优点周围。 -
循环过程中最优解随循环次数的变化曲线：
可以看出，收敛速度明显高于原始算法的速度。