介绍
粒子群算法（Particle swarm optimization，PSO）是模拟群体智能所建立起来的一种优化算法，主要用于解决最优化问题（optimization problems）。1995年由 Eberhart和Kennedy 提出，是基于对鸟群觅食行为的研究和模拟而来的。
假设一群鸟在觅食，在觅食范围内，只在一个地方有食物，所有鸟儿都看不到食物（即不知道食物的具体位置。当然不知道了，知道了就不用觅食了），但是能闻到食物的味道（即能知道食物距离自己是远是近。鸟的嗅觉是很灵敏的）。
假设鸟与鸟之间能共享信息（即互相知道每个鸟离食物多远。这个是人工假定，实际上鸟们肯定不会也不愿意），那么最好的策略就是结合自己离食物最近的位置和鸟群中其他鸟距离食物最近的位置这2个因素综合考虑找到最好的搜索位置。
粒子群算法与《遗传算法》等进化算法有很多相似之处。也需要初始化种群，计算适应度值，通过进化进行迭代等。但是与遗传算法不同，它没有交叉，变异等进化操作。与遗传算法比较，PSO的优势在于很容易编码，需要调整的参数也很少。
一、基本概念
与遗传算法类似，PSO也有几个核心概念。
- 粒子（particle）：一只鸟。类似于遗传算法中的个体。
- 种群（population）：一群鸟。类似于遗传算法中的种群。
- 位置（position）：一个粒子（鸟）当前所在的位置。
- 经验（best）：一个粒子（鸟）自身曾经离食物最近的位置。
- 速度（velocity ）：一个粒子（鸟）飞行的速度。
- 适应度（fitness）：一个粒子（鸟）距离食物的远近。与遗传算法中的适应度类似。
二、粒子群算法的过程
可以看出，粒子群算法的过程比遗传算法还要简单。
1）根据问题需要，随机生成粒子，粒子的数量可自行控制。
2）将粒子组成一个种群。这前2个过程一般合并在一起。
3）计算粒子适应度值。
4）更新种群中每个粒子的位置和速度。
5）满足退出条件就退出，不满足就转向步骤3）。
三、核心——“速度更新”
从上面PSO的算法流程中可以看出，核心步骤是更新种群中每个粒子的位置和速度，而速度的更新又是核心中的核心。
下面直接给出速度更新公式：
new_v = w*v + c1 * rand() * (pbest - position) + c2 * rand() * (gbest-position)
v为粒子当前的速度，w为惯性因子（有速度就有运动惯性）。rand()为随机数生成函数，能够生成0-1之间的随机数。position为粒子当前的位置，pbest为本粒子历史上最好的位置，gbest为种群中所有粒子中当前最好的位置。c1和c2表示学习因子，分别向本粒子历史最好位置和种群中当前最好位置进行学习。
参数好像也有很多，需要设置的是3个，w，c1和c2，但实际上一般都设置c1=c2=1，w一般设在0.5左右。所以也没什么好设置的。
从物理原理上来解释这个速度更新公式，该公式由加号分割为3个部分：
第一部分是惯性保持部分，粒子沿着当前的速度和方向惯性飞行，不会偏移，直来直去。（牛顿运动学第一定理）。
第二部分是自我认知部分，粒子受到自身历史最好位置的吸引力，有回到自身历史最好位置的意愿。（牛顿运动学第二定理）。
第三部分是社会认知部分，粒子处在一个社会中（种群中），社会上有更好的粒子（成功人士），粒子受到成功人士的吸引力，有去社会中成功人士位置的意愿。（牛顿运动学第二定理）。
速度更新公式的意义就是粒子在自身惯性和2种外力作用下，速度和方向发生的改变。
注意这3部分都有重要含义。没有惯性部分，粒子们将很快向当前的自身最优位置和全局最优粒子位置靠拢，变成了一个局部算法了。有了惯性，不同粒子将有在空间中自由飞行的趋势，能够在整个搜索区域内寻找食物（最优解）。而没有自我认知部分，粒子们将向当前的全局最优粒子位置靠拢，容易陷入局部最优。没有社会认知部分，粒子们各自向自身最优位置靠拢，各自陷入自身最优，整个搜索过程都不收敛了。
最后，有了速度更新公式，位置更新就简单了：
new_position = position + new_v * t
t一般默认取1。
四、代码实现
（1）
# -*- coding: utf-8 -*-
"""
粒子群算法求解函数最大值（最小值）
f(x)= x + 10*sin5x + 7*cos4x
"""
import numpy as np
import matplotlib.pyplot as plt
#粒子（鸟）
class Particle:
def __init__(self):
self.p = 0 # 粒子当前位置
self.v = 0 # 粒子当前速度
self.pbest = 0 # 粒子历史最好位置
class PSO:
def __init__(self,N=20,iter_N=100):
self.w = 0.2 # 惯性因子
self.c1 = 1 # 自我认知学习因子
self.c2 = 2 # 社会认知学习因子
self.gbest = 0 # 种群当前最好位置
self.N = N # 种群中粒子数量
self.POP = [] # 种群
self.iter_N = iter_N # 迭代次数
# 适应度值计算函数
def fitness(self,x):
return x + 10 * np.sin(5 * x) + 7 * np.cos(4 * x)
# 找到全局最优解
def g_best(self,pop):
for bird in pop:
if bird.fitness > self.fitness(self.gbest):
self.gbest = bird.p
# 初始化种群
def initPopulation(self,pop,N):
for i in range(N):
bird = Particle()
bird.p = np.random.uniform(-10,10)
bird.fitness = self.fitness(bird.p)
bird.pbest = bird.fitness
pop.append(bird)
# 找到种群中的最优位置
self.g_best(pop)
# 更新速度和位置
def update(self,pop):
for bird in pop:
v = self.w * bird.v + self.c1 * np.random.random() * (bird.pbest - bird.p) + self.c2 * np.random.random() * (self.gbest - bird.p)
p = bird.p + v
if -10 < p < 10:
bird.p = p
bird.v = v
# 更新适应度
bird.fitness = self.fitness(bird.p)
# 是否需要更新本粒子历史最好位置
if bird.fitness > self.fitness(bird.pbest):
bird.pbest = bird.p
def implement(self):
# 初始化种群
self.initPopulation(self.POP,self.N)
def func(x):
return x + 10 * np.sin(5 * x) + 7 * np.cos(4 * x)
x = np.linspace(-10, 10, 1000)
y = func(x)
# 迭代
for i in range(self.iter_N):
# 更新速度和位置
self.update(self.POP)
# 更新种群中最好位置
self.g_best(self.POP)
# 绘制动画
plt.clf()
scatter_x = np.array([ind.p for ind in pso.POP])
scatter_y = np.array([ind.fitness for ind in pso.POP])
scatter_x1 = pso.gbest
scatter_y1 = pso.fitness(pso.gbest)
plt.plot(x, y)
plt.scatter(scatter_x, scatter_y, c='b')
plt.scatter(scatter_x1, scatter_y1, c='r')
plt.pause(0.01)
pso = PSO(N=20,iter_N=100)
pso.implement()
for ind in pso.POP:
print("x = ",ind.p,"f(x) = ",ind.fitness)
print("最优解 x = ",pso.gbest,"相应最大值 f(x) = ",pso.fitness(pso.gbest))
plt.show()（2）
# -*- coding: utf-8 -*-
"""
f(x1,x2) = x1**2 + x2**2, x1,x2 belongs to [-10,10],求Minf
"""
import matplotlib.pyplot as plt
import numpy as np
class PSO(object):
def __init__(self, population_size, max_steps):
self.w = 0.6 # 惯性权重
self.c1 = self.c2 = 2
self.population_size = population_size # 粒子群数量
self.dim = 2 # 搜索空间的维度
self.max_steps = max_steps # 迭代次数
self.x_bound = [-10, 10] # 解空间范围
self.x = np.random.uniform(self.x_bound[0], self.x_bound[1],
(self.population_size, self.dim)) # 初始化粒子群位置
self.v = np.random.rand(self.population_size, self.dim) # 初始化粒子群速度
fitness = self.calculate_fitness(self.x)
self.p = self.x # 个体的最佳位置
self.pg = self.x[np.argmin(fitness)] # 全局最佳位置
self.individual_best_fitness = fitness # 个体的最优适应度
self.global_best_fitness = np.max(fitness) # 全局最佳适应度
def calculate_fitness(self, x):
return np.sum(np.square(x), axis=1)
def evolve(self):
fig = plt.figure()
for step in range(self.max_steps):
r1 = np.random.rand(self.population_size, self.dim)
r2 = np.random.rand(self.population_size, self.dim)
# 更新速度和权重
self.v = self.w*self.v+self.c1*r1*(self.p-self.x)+self.c2*r2*(self.pg-self.x)
self.x = self.v + self.x
plt.clf()
plt.scatter(self.x[:, 0], self.x[:, 1], s=30, color='k')
plt.xlim(self.x_bound[0], self.x_bound[1])
plt.ylim(self.x_bound[0], self.x_bound[1])
plt.pause(0.01)
fitness = self.calculate_fitness(self.x)
# 需要更新的个体
update_id = np.greater(self.individual_best_fitness, fitness)
self.p[update_id] = self.x[update_id]
self.individual_best_fitness[update_id] = fitness[update_id]
# 新一代出现了更小的fitness，所以更新全局最优fitness和位置
if np.min(fitness) < self.global_best_fitness:
self.pg = self.x[np.argmin(fitness)]
self.global_best_fitness = np.min(fitness)
# print('best fitness: %.5f, mean fitness: %.5f' % (self.global_best_fitness, np.mean(fitness)))
pso = PSO(10, 100)
pso.evolve()
plt.show()
