1. 起源
粒子群优化算法（Particle Swarm Optimization，PSO）属于进化算法的一种，是通过模拟鸟群捕食行为设计的。从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
三个属性
- 每只鸟自身的飞行惯性 — 自身惯性贡献（速度与位置）
- 每只鸟的历史最优状态 — 自身认知贡献（pbest-个体最优值）
- 整个鸟群的历史最优状态 — 群体经验贡献（gbest-群体最优值）
此处以与食物的距离来判断飞行状态的优劣。根据此三类贡献，状态空间内每只鸟将逐渐调整自身的飞行速度（包括大小、方向），并向食物位置（即局部香味最浓的位置）汇聚。因此，在相关参数设置合理的前提下，粒子群算法的最终解应该对应于给定状态空间内的最值。
2. 更新公式
粒子群速度更新公式：
Vi(t+1)=ωVi(t)+c1r1(pbesti−Xi(t))+c2r2(gbest−Xi(t))
该式右端三项分别代表：自身惯性贡献、自身认知贡献以及群体经验贡献。其中，ω代表惯性因子，c1、c2代表学习因子，r1、r2为[0,1]之间的均匀随机数，pbesti为第i个粒子已知的历史最优状态或位置，gbest整个粒子群已知的历史最优状态或位置。Vi(t)与Xi(t)分别代表t时刻粒子i的速度与位置。
粒子群位置更新公式：
Xi(t+1)=Xi(t)+Vi(t+1)
3. 使用案例（一元函数）
计算 f = sin(x^2) * (x^2 - 8*x + 6) 在区[-3,3]最大值
- 定义粒子类
# 定义“粒子”类
class parti(object):
def __init__(self, v, x):
self.v = v # 粒子当前速度
self.x = x # 粒子当前位置
self.pbest = x # 粒子历史最优位置
- 设置参数
partisNum= 10 #粒子数
iterMax= 1000 #迭代次数
w= 10 #惯性因子
c1= 0.05 #学习因子
c2= 0.05 #学习因子
v_max = (interval[1] - interval[0]) * 0.1 # 设置最大迁移速度，interval为变量x区间
# r1,r2为两个随机数
- 初始化所有粒子
def initPartis(self, partisNum):
partis_list = list()
for i in range(partisNum):
v_seed = random.uniform(-self.v_max, self.v_max) #速度
x_seed = random.uniform(*self.interval) #位置
partis_list.append(parti(v_seed, x_seed))
temp = 'find_' + self.tab
gbest = find_max(partis_list)
else:
exit('>>>tab标签传参有误："min"|"max"<<<')
return partis_list, gbest
def find_max(self, partis_list):
parti = max(partis_list, key=lambda parti: self.func(parti.pbest)) # 按状态函数最大值找到粒子群初始化的历史最优位置
return parti.pbest
４. 更新粒子状态
def solve(self):
for i in range(self.iterMax):
for parti_c in self.partis_list: # 对每一个粒子进行更新
f1 = self.func(parti_c.x) # 获取x对应函数值
# 更新粒子速度，并限制在最大迁移速度之内
parti_c.v = self.w * parti_c.v + self.c1 * random.random() * (parti_c.pbest - parti_c.x) + self.c2 * random.random() * (self.gbest - parti_c.x)
if parti_c.v > self.v_max:
parti_c.v = self.v_max
elif parti_c.v < -self.v_max:
parti_c.v = -self.v_max
# 更新粒子位置，并限制在待解空间之内
if self.interval[0] <= parti_c.x + parti_c.v <=self.interval[1]: # 保证下一个位置和当前位置在区间内
parti_c.x = parti_c.x + parti_c.v
else:
parti_c.x = parti_c.x - parti_c.v
f2 = self.func(parti_c.x) # 待确定更新位置
deal_max(f1, f2, parti_c) # 更新粒子历史最优位置与群体历史最优位置
def deal_max(self, f1, f2, parti_):
if f2 > f1: # 更新粒子历史最优位置
parti_.pbest = parti_.x
if f2 > self.func(self.gbest):
self.gbest = parti_.x # 更新群体历史最优位置
５. 展示
def display(self):
print('solution: {}'.format(self.gbest))
plt.figure(figsize=(8, 4))
x = np.linspace(self.interval[0], self.interval[1], 300)
y = self.func(x)
plt.plot(x, y, 'g-', label='function')
plt.plot(self.x_seeds, self.func(self.x_seeds), 'b.', label='seeds')
plt.plot(self.gbest, self.func(self.gbest), 'r*', label='solution')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('solution = {}'.format(self.gbest))
plt.legend()
plt.savefig('OptAlgorithm.png', dpi=500)
plt.show()
plt.close()
６. 结果
solution: -1.3355895427301705
4. 使用案例（二元函数）
计算 f = x^2 + y^2 在区[-10,10]最大值，图像如下：
- 定义粒子类
# 定义“粒子”类
def __init__(self, dim, size, iter_num, x_max, max_vel, best_fitness_value=float('Inf'), C1=2, C2=2, W=1):
self.C1 = C1
self.C2 = C2
self.W = W
self.dim = dim # 粒子的维度
self.size = size # 粒子个数
self.iter_num = iter_num # 迭代次数
self.x_max = x_max
self.max_vel = max_vel # 粒子最大速度
self.best_fitness_value = best_fitness_value
self.best_position = [0.0 for i in range(dim)] # 种群最优位置
self.fitness_val_list = [] # 每次迭代最优适应值
self.color = ['b', 'b', 'g', 'k', 'm', 'w', 'y']
# 对种群进行初始化
self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]
- 设置参数
partisNum= 10 #粒子数
iterMax= 1000 #迭代次数
w= 10 #惯性因子
c1= 0.05 #学习因子
c2= 0.05 #学习因子
v_max = (interval[1] - interval[0]) * 0.1 # 设置最大迁移速度，interval为变量x区间
# r1,r2为两个随机数
dim = 2 # 纬度
size = 100 #粒子数
iter_num = 1000 #迭代次数
x_max = 10 # 设置粒子运动区间
max_vel = 0.05 # 设置最大迁移速度
- 初始化所有粒子
class Particle:
# 初始化
def __init__(self, x_max, max_vel, dim):
self.__pos = [random.uniform(-x_max, x_max) for i in range(dim)] # 粒子的位置
self.__vel = [random.uniform(-max_vel, max_vel) for i in range(dim)] # 粒子的速度
self.__bestPos = [0.0 for i in range(dim)] # 粒子最好的位置
self.__fitnessValue = fit_fun(self.__pos) # 适应度函数值
４. 更新粒子状态
def update(self):
for i in range(self.iter_num):
for part in self.Particle_list:
self.update_vel(part) # 更新速度
self.update_pos(part) # 更新位置
self.fitness_val_list.append(self.get_bestFitnessValue()) # 每次迭代完把当前的最优适应度存到列表
# 迭代一次，绘画一次所有粒子
for Part in self.Particle_list:
x = Part.get_pos()[0]
y = Part.get_pos()[1]
x, y = np.meshgrid(x, y)
z = fit_fun([x, y])
col = random.choice(self.color)
self.ax.scatter(x, y, z, c=col, marker='.', s=50, label='', alpha=0.5)
# self.ax.set_title('iter-' + str(i))
plt.pause(0.1)
self.ax.scatter(self.fitness_val_list[0], self.fitness_val_list[1], self.get_bestPosition()[-1], c='r', marker='x', s=50, label='')
return self.fitness_val_list, self.get_bestPosition()
# 更新速度
def update_vel(self, part):
for i in range(self.dim):
vel_value = self.W * part.get_vel()[i] + self.C1 * random.random() * (part.get_best_pos()[i] - part.get_pos()[i]) \
+ self.C2 * random.random() * (self.get_bestPosition()[i] - part.get_pos()[i])
if vel_value > self.max_vel:
vel_value = self.max_vel
elif vel_value < -self.max_vel:
vel_value = -self.max_vel
part.set_vel(i, vel_value)
# 更新位置
def update_pos(self, part):
for i in range(self.dim):
pos_value = part.get_pos()[i] + part.get_vel()[i]
part.set_pos(i, pos_value)
value = fit_fun(part.get_pos())
if value < part.get_fitness_value():
part.set_fitness_value(value)
for i in range(self.dim):
part.set_best_pos(i, part.get_pos()[i])
if value < self.get_bestFitnessValue():
self.set_bestFitnessValue(value)
for i in range(self.dim):
self.set_bestPosition(i, part.get_pos()[i])
５. 展示
for Part in self.Particle_list:
x = Part.get_pos()[0]
y = Part.get_pos()[1]
x, y = np.meshgrid(x, y)
z = fit_fun([x, y])
col = random.choice(self.color)
self.ax.scatter(x, y, z, c=col, marker='.', s=50, label='', alpha=0.5)
# self.ax.set_title('iter-' + str(i))
plt.pause(0.1)
６. 结果
PSO最优位置:[5.661373470883246 10^-6, -2.38440134194550810^-5]
PSO最优解:6.005881255239746 * 10^-10
扫一扫
980
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
