粒子群算法

1、 粒子群算法概述

  粒子群算法（Particle Swarm Optimization，PSO）由Kennedy和Eberhart于1995年提出。该算法的思想来源于对鸟类捕食行为的研究，鸟之间通过集体的协作使得群体能够找到最多的食物，PSO便是通过模拟鸟群飞行觅食的行为，来寻找最优解的算法，这是一种基于群体智能（Swarm Intelligence）的优化方法。在粒子群算法中，我们将鸟群抽象成粒子群，用一个粒子来代表一只鸟。PSO目标就是：通过一群粒子在解空间中进行搜索，找到使得适应度函数（fitness function）取得最大值（或最小值）的解。

2、算法介绍

  在正式介绍算法前，先约定好一些符号。假设粒子在 d 维空间中进行搜索，种群大小为N（一共有N个粒子），适应度函数为f(⋅)。
   x_i^(k) ∈ R^d 表示第 i 个粒子在迭代 k 次后的位置
   v_i^(k) ∈ R^d 表示第 i 个粒子的在迭代 k 次后的速度
   pbest_i ∈ R^d 表示第 i 个粒子所经过的历史最优位置
   gbest ∈ R^d 表示整个种群所找到的历史最优位置

则第 i 个粒子位置更新公式为：

其中
   w 是惯性系数，w 越大说明粒子更倾向于保持原来的运动状态
   c₁ 是自我认知系数，c₁ 越大说明粒子更倾向于相信自己的经验（认知）
   c₂ 是全局认知系数，c₂ 越大说明粒子更倾向于相信整个群体的经验（认知）
   r₁, r₂ 是随机数，通常是 [0,1] 上均匀产生的随机数

从更新公式可以看出，粒子下一时刻的运动速度收到三个因素的影响：当前时刻的速度、当前自己找到的最优解的位置以及当前全局最优解的位置。可以用一个图表示粒子位置的更新过程：

3、算法步骤

基本算法描述如下：

   1、首先，人为设置粒子的种群大小 N ，适应度函数f(⋅)，惯性系数 ω ，以及 c₁, c₂；
   2、随机初始化每个粒子的位置 x_i⁽⁰⁾，初始速度v_i⁽⁰⁾ ；然后记录当前每个粒子的历史最优解和全局最优解。
   3、不断循环迭代，每循环一次就更新种群中每个粒子的状态，直至达到算法的停止条件。最后，输出所找到的最优解。

算法的停止条件：

   1、最简单的就是直接设置一个最大迭代上限，超出最大迭代次数后直接退出。
   2、设置一个计数器，并设置一个阈值 C 。如果，PSO在搜索过程中，连续 C 次循环所找到的全局最优解都没有发生变化，就停止。
   3、设置一个阈值 T ，如果PSO所找到的全局最优的解满足 f(gbest) ≥ T ，那么停止算法。

4、实例： 用PSO算法搜索ackley函数的最小值

   ackley函数是一个具有非常多个局部极小值的函数，具体表达式如下：

通常取 a=20,b=0.2,c=2π 。其中 d 表示空间的维数，即 x∈R^d。该函数具有全局最小值 f(x_∗) = 0, x_∗=0 。函数图像如下（三维空间）：

.

完整代码

import numpy as np


# v(k+1) = w*v(k) + c1*r1*(gbest-x) + c2*r2*(pbest-x)
# x(k+1) = x(k) + v(k+1)


class Particle:
    """
    st_x：位置约束，粒子每个维度上的坐标范围必须处于[st_x[0], st_x[1]]之间
    st_v：速度约束，粒子每个维度上的速度范围必须处于[st_v[0], st_v[1]]之间
    position：粒子的当前位置
    velocity：粒子的当前速度
    pbest：粒子自身历史记录的最佳位置
    pvalue：粒子自身历史记录的最佳值
    """
    st_x = None
    st_v = None

    def __init__(self, x, v, pbest, pvalue):
        self.position = x
        self.velocity = v
        self.pbest = pbest
        self.pvalue = pvalue

    def update_velocity(self, v):
        v[v < self.st_v[0]] = self.st_v[0]
        v[v > self.st_v[1]] = self.st_v[1]
        self.velocity = v

    def update_position(self, func):
        # 更新粒子自身的位置，以及判断是否要更新pbest
        self.position += self.velocity
        self.position[self.position < self.st_x[0]] = self.st_x[0]
        self.position[self.position > self.st_x[1]] = self.st_x[1]
        if self.pvalue > func(self.position):
            self.pvalue = func(self.position)
            self.pbest = self.position.copy()


class PSO:
    """
    gbest：粒子群历史记录的最佳位置
    gvalue：粒子群历史记录的最佳值
    """
    gbest = None
    gvalue = np.inf

    def __init__(self, n_dims, n_particles, st_x, st_v, w, c1, c2, num_iter, func):
        # 初始化空间维度
        self.n_dims = n_dims
        # 初始化粒子群数目
        self.n_particles = n_particles
        # 目标函数
        self.func = func
        # 粒子惯性权重
        self.w = w
        # 全局部分学习率
        self.c1 = c1
        # 自我认知部分学习率
        self.c2 = c2
        # 迭代次数
        self.num_iter = num_iter
        # 存放粒子群的列表（容器）
        self.particles = []
        # 初始化粒子的位置和速度约束
        Particle.st_x = st_x
        Particle.st_v = st_v
        # 初始化粒子群
        for _ in range(n_particles):
            # 初始化粒子的随机位置在 st_x[0]~st_x[1] 
            x = (st_x[1] - st_x[0]) * np.random.rand(n_dims) + st_x[0]
            # 计算当前评估值
            pvalue = func(x)
            # 初始化一个粒子
            self.particles.append(
                Particle(
                    x=x,
                    v=(st_v[1] - st_v[0]) * np.random.rand(n_dims) + st_v[0],
                    pbest=x.copy(),
                    pvalue=pvalue
                )
            )

            if self.gvalue > pvalue:
                self.gvalue = pvalue
                self.gbest = x.copy()

    def solve(self):
        # 开始迭代
        for index in range(1, self.num_iter + 1):
            for particle in self.particles:
                v = self.w * particle.velocity + self.c1 * np.random.rand() * (self.gbest - particle.position) + \
                    self.c2 * np.random.rand() * (particle.pbest - particle.position)
                particle.update_velocity(v)
                particle.update_position(self.func)
            for particle in self.particles:
                if particle.pvalue < self.gvalue:
                    self.gvalue = particle.pvalue
                    self.gbest = particle.pbest.copy()
        return self.gbest, self.gvalue


def ackley(x):
    return - 20 * np.exp(-0.2 * np.sqrt((x * x).mean())) - np.exp(np.cos(2 * np.pi * x).mean()) + 20 + np.exp(1)


if __name__ == "__main__":
    # 测试不同迭代次数搜索出来的最优解情况
    for i in range(11):
        pso = PSO(
            n_dims=10,
            n_particles=100,
            st_x=(-20, 20),
            st_v=(-1, 1),
            w=0.8,
            c1=2,
            c2=2,
            num_iter=50 * i,
            func=ackley
        )
        gbest, gvalue = pso.solve()
        print("number of interations:%d\tgvalue:%f" % (50 * i, gvalue))