【智能算法】粒子群寻优算法

1.理论基础

粒子群算法（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）源于对鸟类捕食行为的研究，鸟类捕食时，找到食物最简单有限的策略就是搜寻当前距离食物最近的鸟的周围。

举个通俗的例子：

一群鸟在寻找食物，在这个区域内有一个食物，所有的鸟都不知道食物在哪里，但食物发出了香味，鸟通过香味的浓烈能判断出自己的当前位置距离食物有多远，同时鸟群之间是可以交流并告知离食物最近的鸟的位置的。试想一下这个时候会发生什么？

鸟甲：哈哈哈，原来我离食物最近！

鸟乙、丙、丁…:我得赶紧往鸟甲那里过去看看！

同时各自鸟在位置不停变化的时候，离食物的距离也不断的变化，所以一定有过离食物最近的位置，这是他们的一个参考。

鸟某某：我刚才的位置好像离食物很近了，我得往那里再靠近点！

通过这样一个过程，不断的变化，最终鸟群会向食物方向聚集，达到目标。在这个变化的过程中，影响鸟的运动状态变化的有两个因素：

离食物最近的鸟的位置
自己之前达到过的离食物最近的位置

而鸟的每次的位置变化，除了考虑以上的两个因素还有一个因素就是惯性。而对位置的变化，通过变化量（速度）来表示。

通过以上的例子，我们可以把鸟抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子I 在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)．每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi．这个可以看作是粒子自己的飞行经验．除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)．这个可以看作是粒子同伴的经验．粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。归结为下面的两个公式：

$[{ m{V}}_{id}^{k + 1}{ m{ = }}omega { m{V}}_{id}^k + {c_1}{r_1}(P_{id}^k - X_{id}^k) + {c_2}{r_2}(P_{gd}^k + X_{id}^k)]$ (公式-1）

$[{ m{X}}_{id}^{k + 1} = X_{id}^k + {V_{k + 1}}_{id}]$ （公式-2）

其中， $[omega ]$ 为惯性权重因子，（什么是权重因子，这个后面会讲到，此处大概知道有这么一个参数即可）；d=1,2，…，D；i=1,2，…，n；k为当前迭代次数；V_id为粒子的速度；c₁和c₂是非负的常数，称为加速因子；r₁和r₂是分布于[0,1]区间的随机数。为了防止粒子的满目搜索，一般建议限制位置和速度在一定的区间。

2.粒子群寻优算法的案例与实现

假设有如下非线性函数，需要在一定范围内找出最大值和最大值的位置，本文用Python实现。

$[f(x) = frac{{sin sqrt {{x^2} + {y^2}} }}{{sqrt {{x^2} + {y^2}} }} + {e^{frac{{cos 2pi x + cos 2pi y}}{2}}} - 2.71289]$ （公式-3）

我们先通过python把上面函数使用matplotlib把它画处理，让我们有个直观的感受。

# -*- coding: utf-8 -*-
from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-2, 2, 0.05)
Y = np.arange(-2, 2, 0.05)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2+Y**2))/np.sqrt(X**2+Y**2)+np.exp((np.cos(2*np.pi*X)+np.cos(2*np.pi*Y))/2)-2.71289

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='hot')
plt.show()

此处使用mpl_toolkits.mplot3d画3D图，在代码中，分别使用numpy中的arange创建等差长度为0.05的数据点向量，区间为[-2,2]之间，然后使用meshgrid将向量转换给对应的矩阵，这样的话，后面的函数式中就可以使用np中的三角函数和数学函数了，如果不这样的话，还得一个一个点的根据公式去做循环，拼接Z的向量点。

运行出来的图如下：

从函数图形中可以看出，这个函数有很多的局部极大值，而真正的极限位置就是在（0,0），极大值为1.005。下面就通过粒子群算法来找寻这个极大值和极大值的位置。根据粒子群算法函数求极值的算法流程如下图所示：

在pycham作为python的IDE中，建立一个名为PSO.py的文件，文件中编写一个名PSO的类。

2.1 参数的设定

设置惯性参数为常数1，也就是公式-1中的 $[omega ]$ 为常数1；设置加速度因子为非负的两个常数，公式-1的中的c₁和c₂；设置迭代次数为300；设置粒子群的规模为20；为了防止粒子的盲目搜索，将位置和速度限制在区间[-2,2]、[-0.5,0.5]；

具体代码如下：

def __init__(self):
        self.w = self.getweight()
        self.lr = self.getlearningrate()
        self.maxgen = self.getmaxgen()
        self.sizepop = self.getsizepop()
        self.rangepop = self.getrangepop()
        self.rangespeed = self.getrangespeed()

    def getweight(self):
        # 惯性权重
        weight = 1
        return weight

    def getlearningrate(self):
        # 分别是粒子的个体和社会的学习因子，也称为加速常数
        lr = (0.49445,1.49445)
        return lr

    def getmaxgen(self):
        # 最大迭代次数
        maxgen = 300
        return maxgen

    def getsizepop(self):
        # 种群规模
        sizepop = 20
        return sizepop

    def getrangepop(self):
        # 粒子的位置的范围限制,x、y方向的限制相同
        rangepop = (-2,2)
        return rangepop

    def getrangespeed(self):
        # 粒子的速度范围限制
        rangespeed = (-0.5,0.5)
        return rangespeed

2.2 种群初始化

随机初始化粒子位置和粒子的速度，并根据适应度函数计算粒子适应度值，我们的适应度函数就是公式-3这个函数。

首先定义一个适应度函数

def func(self,x):
        # x输入粒子位置
        # y 粒子适应度值
        if (x[0]==0)&(x[1]==0):
            y = np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        else:
            y = np.sin(np.sqrt(x[0]**2+x[1]**2))/np.sqrt(x[0]**2+x[1]**2)+np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        return y

然后初始化粒子和计算适应度值

def initpopvfit(self,sizepop):
        pop = np.zeros((sizepop,2))
        v = np.zeros((sizepop,2))
        fitness = np.zeros(sizepop)

        for i in range(sizepop):
            pop[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            v[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            fitness[i] = self.func(pop[i])
        return pop,v,fitness

2.3 寻找初始化后的极值

def getinitbest(self,fitness,pop):
        # 群体最优的粒子位置及其适应度值
        gbestpop,gbestfitness = pop[fitness.argmax()].copy(),fitness.max()
        #个体最优的粒子位置及其适应度值,使用copy()使得对pop的改变不影响pbestpop，pbestfitness类似
        pbestpop,pbestfitness = pop.copy(),fitness.copy()

        return gbestpop,gbestfitness,pbestpop,pbestfitness

此处需要分别找到群体极值和个体粒子的极值。

2.4 迭代寻优

通过循环迭代，不断的更新粒子的位置和速度，根据新粒子的适应度值更新个体和群体的极值。这部分的代码如下：

def run(self):
        pop, v, fitness = self.initpopvfit(self.sizepop)
        gbestpop, gbestfitness, pbestpop, pbestfitness = self.getinitbest(fitness, pop)

        result = np.zeros(self.maxgen)
        for i in range(self.maxgen):
            # 速度更新
            for j in range(self.sizepop):
                v[j] += self.lr[0] * np.random.rand() * (pbestpop[j] - pop[j]) + self.lr[1] * np.random.rand() * (
                gbestpop - pop[j])
            v[v < self.rangespeed[0]] = self.rangespeed[0]
            v[v > self.rangespeed[1]] = self.rangespeed[1]

            # 粒子位置更新
            for j in range(self.sizepop):
                pop[j] += 0.5 * v[j]
            pop[pop < self.rangepop[0]] = self.rangepop[0]
            pop[pop > self.rangepop[1]] = self.rangepop[1]

            # 适应度更新
            for j in range(self.sizepop):
                fitness[j] = self.func(pop[j])

            for j in range(self.sizepop):
                if fitness[j] > pbestfitness[j]:
                    pbestfitness[j] = fitness[j]
                    pbestpop[j] = pop[j].copy()

            if pbestfitness.max() > gbestfitness:
                gbestfitness = pbestfitness.max()
                gbestpop = pop[pbestfitness.argmax()].copy()

            result[i] = gbestfitness
        return result

2.4 结果分析

通过在主函数文件中，实例化PSO类，调用run这个函数，运行上面的代码

from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from PSO import *

pso = PSO()
result = pso.run()
plt.figure()
plt.plot(result)
plt.show()

当迭代了300次之后，画出每代最优个体适应度值的曲线，结果如下：

最终得到的最优个体适应度值为：1.00538866933，对应的粒子位置为：[ 2.37513360e-05 3.41265823e-04]，通过比较，PSO算法能得到最优值，接近函数实际的最优值1.005。

通过以上，证明PSO算法具有较强的函数极值寻优能力。

3.延伸与算法改进

讲到这里，我们再来回顾一下之前粒子更新速度（公式-1）中的 $[omega ]$ 这个惯性权重因子，所谓惯性权重因子，体现的是粒子继承上一次迭代速度的能力。较大的一个惯性权重因子有利于全局的搜索，而一个小的惯性权重因子则更加利于局部的搜索。为了更好的平衡算法的全局搜索与局部搜索的能力，常给惯性权重因子改为随迭代次数变化的一个函数，一般常用的惯性权重因子函数有以下几种：

$[omega (k) = {omega _{start}}-({omega _{start}} - {omega _{end}})({T_{max }} - k)/{T_{max }}]$ （公式-4）

（公式-5）

（公式-6）

（公式-7）

这几种 $[omega ]$ 的变化图如下：

从变化图中，可以看到（公式-5）中，前期 $[omega ]$ 变化较慢，取值较大，维持了算法的全局搜索能力；后期 $[omega ]$ 变化较快，极大的提高了算法的局部寻优能力，通过这种 $[omega ]$ 权重的变化函数，可以取得比 $[omega ]$ 固定不变的那种情况要好的效果。