如何用粒子群算法求解局部最优值
粒子群算法的原理和求解思路
粒子群算法基本思想是经验学习和群体学习。在核辐射探测中,有一种仪器叫做盖革计数器,切尔诺贝利事件中,施工人员多半会佩戴盖革计数器。盖革计数器响了,代表人所处点位的受辐射水平高;施工人员可以通过对讲机相互交流,团队中所有人都能知道谁所处的点位有超标辐射。这只是一个抽象的例子,但是可以基本了解粒子群算法的基本思路。
用鸟群捕食做更形象的描述。一个鸟群在一片区域中捕食一只兔子,每只鸟都不知道这只兔子在哪里,但是都知道自己离这只兔子的距离是多少。每只鸟都会记下自己的历史最优位置;此外,鸟群可以相互交流,大家都可以知道谁有最佳的历史最优点(即当前的全局历史最优点)。在调整位置的过程中,每只鸟都会参考自己的历史最优位置和当前的全局历史最优位置,从而挪动到新的位置。
因此,粒子群算法的基本思路是这样的式子:
v
i
=
ω
×
v
i
+
c
1
×
r
a
n
d
o
m
(
0
,
1
)
×
(
h
b
e
s
t
i
−
p
i
)
+
c
2
×
r
a
n
d
o
m
(
0
,
1
)
×
(
g
b
e
s
t
i
−
p
i
)
v_{i} = \omega\times v_{i}+c_{1}\times random(0, 1)\times(hbest_{i}-p_{i})+c_{2}\times random(0, 1)\times(gbest_{i}-p_{i})
vi=ω×vi+c1×random(0,1)×(hbesti−pi)+c2×random(0,1)×(gbesti−pi)
其中,
v
i
v_{i}
vi 是速度向量,具有大小和方向,
h
b
e
s
t
i
hbest_{i}
hbesti 是个体
i
i
i 的历史最优位置,
g
b
e
s
t
i
gbest_{i}
gbesti 是当前全局历史最优位置,
p
i
p_{i}
pi 是个体当前位置。
ω
\omega
ω 被称为“惯性系数”,形象地描述了改变当前速度的“积极程度”;
c
1
、
c
2
c_1、c_2
c1、c2 被称为加速系数,描述了两方向的分速度改变的快慢。经过不断迭代,每个点更新自己的坐标,最后得到优化后的坐标点解和最优值。
用粒子群算法求解局部最小值问题
求解 y = f ( x 1 , x 2 ) = x 1 2 + x 2 2 y=f(x_1, x_2)=x_1^2 + x_2^2 y=f(x1,x2)=x12+x22 的最小值,其中 − 10 ≤ x 1 , x 2 ≤ 10 -10\leq x_1,x_2\leq10 −10≤x1,x2≤10。
'''
全局版 PSO 算法
粒子数 n = 3
惯性权重 omega = 0.5
加速系数 c1 = c2 = 2
'''
import numpy as np
from tqdm import tqdm # 监视程序运行进程
def func(x):
return (x ** 2).sum(axis=0)
def Initialize(n, dimension, lowerBound=-10, upperBound=10): # 点数 维度 下界 上界
return np.random.uniform(lowerBound, upperBound, (n, dimension))
def scoreEvaluation(func, matrix): # 转化为评分
return np.apply_along_axis(func, 1, matrix)
def historyBestPoint(matrix, presentValue, historyBestValue, historyBestPosition): # 更新和记录历史最优点
historyBestPosition[presentValue < historyBestValue] = matrix[presentValue < historyBestValue]
historyBestValue[historyBestValue > presentValue] = presentValue[historyBestValue > presentValue]
return historyBestPosition, historyBestValue
def Diversion(matrix, historyBestPosition, globalBestPosition, omega, c1, c2, lowerBound=-10, upperBound=10):
matrix = omega * matrix + \
c1 * np.random.random() * (historyBestPosition - matrix) + \
c2 * np.random.random() * (globalBestPosition - matrix)
# 强制性调整范围防止越界
matrix[matrix < lowerBound] = lowerBound
matrix[matrix > upperBound] = upperBound
return matrix
在主函数中调用模型粒子群算法:
if __name__ == '__main__':
omega = 0.5 # 惯性系数
c1 = c2 = 2 # 加速系数
n = 3 # 粒子数量 粒子为二维粒子(x1, x2)
PointGroup = Initialize(n, 2) # 上下界已经有默认参数了
print('初始点集为:\n', PointGroup)
# 初始化个体历史最优解(第一次迭代前,当前最优解就是历史最优解和全局历史最优解)
historyBestValue = scoreEvaluation(func, PointGroup)
historyBestPosition = PointGroup
# 初始化全局历史最优解
globalBestValue = historyBestValue.min()
globalBestPosition = historyBestPosition[historyBestValue.argmin()]
# 迭代 10000 次
for i in tqdm(range(10000)):
# 计算现值
presentValue = scoreEvaluation(func, PointGroup)
# 决定是否替换全局最优位置
if globalBestValue > presentValue.min():
globalBestPosition = PointGroup[presentValue.argmin()]
globalBestValue = presentValue.min()
# 更新历史最优位置
historyBestPosition, historyBestValue = historyBestPoint(PointGroup, presentValue, historyBestValue, historyBestPosition)
# 迭代更新点位置
PointGroup = Diversion(PointGroup, historyBestPosition, globalBestPosition, omega, c1, c2)
# 输出最后的全局最优解
print(f'(近似)最优解为: ', globalBestPosition)
print(f'(近似)最优值为: ', globalBestValue)
最后,得到的结果并作图示意:
下图中红点即为用粒子群算法求解得到的全局最优点。