算法背景

粒子群算法（particle swarm optimization，PSO）是计算智能领域中的一种生物启发式方法，属于群体智能优化算法的一种，常见的群体智能优化算法主要有如下几类：

• （1）蚁群算法（Ant Colony Optimization，简称ACO）[1992年提出]；
• （2）粒子群优化算法（Particle Swarm Optimization，简称PSO）[1995年提出]（简单易于实现，也是目前应用最为广泛的群体智能优化算法）；
• （3）菌群优化算法（Bacterial Foraging Optimization，简称BFO）[2002年提出]；
• （4）蛙跳算法（Shuffled Frog Leading Algorithm，简称SFLA）[2003年提出]；
• （5）人工蜂群算法（Artificial Bee Colony Algorithm，简称ABC）[2005年提出]；

除了上述几种常见的群体智能算法以外，还有一些并不是广泛应用的群体智能算法，比如萤火虫算法、布谷鸟算法、蝙蝠算法以及磷虾群算法等等。

而其中的粒子群优化算法（PSO）源于对鸟类捕食行为的研究，鸟类捕食时，找到食物最简单有限的策略就是搜寻当前距离食物最近的鸟的周围。

• 相较于传统算法计算速度非常快，全局搜索能力也很强；
• PSO对于种群大小不十分敏感，所以初始种群设为500-1000，速度影响也不大；
• 粒子群算法适用于连续函数极值问题，对于非线性、多峰问题均有较强的全局搜索能力。

设想这样一个场景：一群鸟在随机的搜索食物。在这个区域里只有一块食物，所有的鸟都不知道食物在哪。但是它们知道自己当前的位置距离食物还有多远。那么找到食物的最优策略是什么？最简单有效的就是搜寻目前离食物最近的鸟的周围区域。

算法实现

算法流程图

Step1:确定一个粒子的运动状态是利用位置和速度两个参数描述的，因此初始化的也是这两个参数；
Step2:每次搜寻的结果（函数值）即为粒子适应度，然后记录每个粒子的个体历史最优位置和群体的历史最优位置；
Step3:个体历史最优位置和群体的历史最优位置相当于产生了两个力，结合粒子本身的惯性共同影响粒子的运动状态，由此来更新粒子的位置和速度。

粒子群算法流程图

初始化种群

位置和速度的初始化即在位置和速度限制内随机生成一个N x d 的矩阵，而对于速度则不用考虑约束，一般直接在0~1内随机生成一个50x1的数据矩阵。

此处的位置约束也可以理解为位置限制，而速度限制是保证粒子步长不超限制的，一般设置速度限制为[-1,1]。

粒子群的另一个特点就是记录每个个体的历史最优和种群的历史最优，因此而二者对应的最优位置和最优值也需要初始化。其中每个个体的历史最优位置可以先初始化为当前位置,而种群的历史最优位置则可初始化为原点。对于最优值，如果求最大值则初始化为负无穷，相反地初始化为正无穷。

每次搜寻都需要将当前的适应度和最优解同历史的记录值进行对比，如果超过历史最优值，则更新个体和种群的历史最优位置和最优解。

速度与位置更新

速度和位置更新是粒子群算法的核心，其原理表达式和更新方式：

速度和位置更新表达式

上面公式中：i表示粒子编号；d表示时刻，反映在迭代次数上；w是惯性权重，一般设置在0.4左右；c1、c2为学习因子，一般取2；pid表示的是粒子i的经验最佳位置；gid代表的是全局粒子的最优位置；r1、r2是0到1之间的随机值

每次更新完速度和位置都需要考虑速度和位置的限制，需要将其限制在规定范围内，此处仅举出一个常规方法，即将超约束的数据约束到边界（当位置或者速度超出初始化限制时，将其拉回靠近的边界处）。当然，你不用担心他会停住不动，因为每个粒子还有惯性和其他两个参数的影响。

代码实现

源程序

粒子群算法求平方和函数最小值，由于没有特意指定函数自变量量纲，不进行数据归一化。

#!/usr/bin/env python 
# -*- coding: UTF-8 -*-

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
 
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 = 3  # 搜索空间的维度
        self.max_steps = max_steps  # 迭代次数
        self.x_bound = [-10, 10]  # 解空间范围
        self.x = np.random.uniform(self.x_bound[0], self.x_bound[1],size= (self.population_size, self.dim))  # 初始化粒子群位置
        self.v = np.random.rand(self.population_size, self.dim)  # 初始化粒子群速度
        fitness = self.function_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 function_fitness(self, x):
        return np.sum(np.square(x), axis=1)
#        return np.sin(np.sqrt(x[0]**2+x[1]**2))
 
    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.subplot(111, projection='3d') 
            plt.scatter(self.x[:, 0], self.x[:, 1], self.x[:, 2], 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.function_fitness(self.x)
            # 需要更新的个体
            update_bool = np.greater(self.individual_best_fitness, fitness) # np.greater(a,b)相当于>，返回一个bool值
            update_id = np.argwhere(update_bool == True)
            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(100, 100)
pso.evolve()
plt.show()

python库调用

np.greater(a,b)相当于a>b的判断语句，返回一个bool值

np.argwhere()返回符合条件的元素的地址

3D图形的绘制通过调用mpl_tookits下的mplot3d类进行，
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow') #绘面
ax.scatter(x, y, z, c='r') #绘点