智能优化算法--粒子群算法

1、引言

自然界中的鸟群和鱼群的群体行为一直是科学家的研究兴趣所在。

1990年，生物学家FrankHeppner提出了鸟类模型。在仿真中，一开始每一只鸟都没有特定的飞行目标，只是使用简单的规则确定自己的飞行方向和飞行速度，当有一只鸟飞到栖息地时，它周围的鸟也会跟着飞向栖息地，最终整个鸟群都会落在栖息地。

1995年，美国社会心理学家James Kennedy和电气工程师Russell Eberhart共同提出了粒子群算法（Particle SwarmOptimization，PSO），该算法的提出是受对鸟类群体行为进行建模与仿真的研究结果的启发。他们的模型和仿真算法主要对Frank Heppner的模型进行了修正，以使粒子飞向解空间并在最优解处降落。

2、粒子群算法理论

2.1 算法描述

粒子群算法受鸟类捕食行为的启发并对这种行为进行模仿，将优化问题的搜索空间类比于鸟类的飞行空间，将每只鸟抽象为一个粒子，粒子无质量、无体积，用以表征问题的一个可行解，优化问题所要搜索到的最优解则等同于鸟类寻找的食物源。

2.2 算法建模

粒子群算法首先在给定的解空间中随机初始化粒子群，待优化问题的变量数决定了解空间的维数。每个粒子有了初始位置与初始速度，然后通过迭代寻优。在每一次迭代中，每个粒子通过跟踪两个“极值”来更新自己在解空间中的空间位置与飞行速度：一个极值就是单个粒子本身在迭代过程中找到的最优解粒子，这个粒子叫作个体极值；另一个极值是种群所有粒子在迭代过程中所找到的最优解粒子，这个粒子是全局极值。上述的方法叫作全局粒子群算法。如果不用种群所有粒子而只用其中一部分作为该粒子的邻居粒子，那么在所有邻居粒子中的极值就是局部极值，该方法称为局部粒子群算法。

2.3粒子群算法种类

2.3.1 基本粒子群算法

假设在一个D维的目标搜索空间中，有N个粒子组成一个群落，其中第i个粒子表示为一个D维的向量：

$${\mathbf{X}}_i=\left(x_{i1},x_{i2},\cdots ,x_{iD}\right),i=1,2,\cdots ,N$$

第i个粒子的“飞行”速度也是一个D维的向量，记为

$${\mathbf{V}}_i=\left(v_{i1},v_{i2},\cdots ,v_{iD}\right),i=1,2,\cdots ,N$$

第i个粒子迄今为止搜索到的最优位置称为个体极值，记为

$${\mathbf{p}}_{\text{best }}=\left(p_{i1},p_{i2},\cdots ,p_{iD}\right),i=1,2,\cdots ,N$$

整个粒子群迄今为止搜索到的最优位置为全局极值，记为

$${\mathbf{g}}_{\text{best }}=\left(g_1,g_2,\cdots ,g_D\right)$$

更新公式为

$$\begin{array}{c}{v}_{ij}(t+1)=v_{ij}(t)+c_1{r}_1(t)\left[p_{ij}(t)-x_{ij}(t)\right]+c_2{r}_2(t)\left[p_{gj}(t)-x_{ij}(t)\right]\\ x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)\end{array}$$
其中，
$c_1$ 和 $c_2$ 为学习因子，也称加速常数；

$r_1$ 和 $r_2$ 为［0,1］范围内的均匀随机数。

速度更新公式由三部分组成：第一部分为“惯性”或“动量”部分，反映了粒子的运动“习惯”，代表粒子有维持自己先前速度的趋势；第二部分为“认知”部分，反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为“社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势。

2.3.2 标准粒子群算法

1998年，Shi Yuhui等人提出了带有惯性权重的改进粒子群算法，
$$\begin{array}{c}{v}_{ij}(t+1)=w\cdot v_{ij}(t)+c_1{r}_1(t)\left[p_{ij}(t)-x_{ij}(t)\right]+c_2{r}_2(t)\left[p_{gi}(t)-x_{ij}(t)\right]\\ x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)\end{array}$$

第一部分表示粒子先前的速度，用于保证算法的全局收敛性能；第二部分、第三部分则使算法具有局部收敛能力。

实验结果表明：w在0.8～1.2之间时，粒子群算法有更快的收敛速度；而当w＞1.2时，算法则容易陷入局部极值。

另外，在搜索过程中可以对w进行动态调整：在算法开始时，可给w赋予较大正值，随着搜索的进行，可以线性地使w逐渐减小，这样可以保证在算法开始时，各粒子能够以较大的速度步长在全局范围内探测到较好的区域；而在搜索后期，较小的w值则保证粒子能够在极值点周围进行精细的搜索，从而使算法有较大的概率向全局最优解位置收敛。

$$w=w_{\max }-\frac{\left(w_{\max }-w_{\min }\right)\cdot t}{T_{\max }}$$
其中，$T_{max}$ 表示最大进化代数，$w_{min}$ 表示最小惯性权重，$w_{max}$ 表示最大惯性权重，t 表示当前迭代次数。在大多数的应用中，$w_{\max }=0.9,w_{\min }=0.4$。

2.3.3 压缩因子粒子群算法

2.3.4 离散粒子群算法

3、粒子群算法流程

（1）初始化粒子群，包括群体规模 N，每个粒子的位置 $x_i$ 和速度 $v_i$。

（2）计算每个粒子的适应度（即目标函数）值 fit[i]。

（3）对每个粒子，用它的适应度值 fit[i] 和个体极值 $p_{best}(i)$ 比较。如果 fit[i] < $p_{best}(i)$，则用 fit[i] 替换掉 $p_{best}(i)$。

（4）对每个粒子，用它的适应度值 fit[i] 和全局极值 $g_{best}$ 比较。如果fit[i] < $p_{best}$，则用 fit[i] 替换 $p_{best}$。

（5）迭代更新粒子的速度 $v_i$ 和位置 $x_i$ 。

（6）进行边界条件处理。

（7）判断算法终止条件是否满足：若是，则结束算法并输出优化结果；否则返回步骤（2）。

4、关键参数说明

（1）粒子种群规模 N

粒子种群大小的选择视具体问题而定，但是一般设置粒子数为20～50。粒子数目越大，算法搜索的空间范围就越大，也就更容易发现全局最优解；当然，算法运行的时间也越长。

（2）惯性权重 w

惯性权重的大小表示了对粒子当前速度继承的多少。当惯性权重值较大时，全局寻优能力较强，局部寻优能力较弱；当惯性权重值较小时，全局寻优能力较弱，局部寻优能力较强。一般取值范围为［0.8，1.2］。

（3）加速常数 $c_1$ 和 $c_2$

如果 $c_1=0$，则为“社会”模型，粒子缺乏认知能力，而只有群体经验，它的收敛速度较快，但容易陷入局部最优；如果 $c_2=0$，则为“认知”模型，没有社会的共享信息，个体之间没有信息的交互，所以找到最优解的概率较小。一般设置$c_2=c_2$，通常可以取 $c_1=c_2=1.5$。

5、仿真实例

计算 $f(x)=x^4-2x+3$ 的最小值

代码

# @Time : 2021/11/16 19:29
# @Author : xiao cong
# @Function : 粒子群算法

import numpy as np
import matplotlib.pyplot as plt


class PSO:
    def __init__(self, pN, dim, max_iter):
        self.w = 0.8
        self.c1 = 2
        self.c2 = 2
        self.r1 = 0.6
        self.r2 = 0.3
        self.pN = pN  # 粒子数量
        self.dim = dim  # 维度
        self.max_iter = max_iter  # 迭代次数
        self.X = np.zeros((self.pN, self.dim))  # 位置
        self.v = np.zeros((self.pN, self.dim))  # 速度
        self.pbest = np.zeros((self.pN, self.dim))  # 个体经历的最佳位置
        self.gbest = np.zeros((1, self.dim))  # 全局最佳位置
        self.p_fit = np.zeros(self.pN)  # 每个个体的历史最佳适应值
        self.g_fit = 0  # 全局最佳适应值

    # 定义目标函数
    def function(self, X):
        return np.float(X ** 4 - 2 * X + 3)

    # 初始化种群
    def init_population(self):
        for i in range(self.pN):
            for j in range(self.dim):
                self.X[i][j] = np.random.uniform(0, 1)
                self.v[i][j] = np.random.uniform(0, 1)

            self.pbest[i] = self.X[i]  # 初始化第i个粒子的位置
            self.p_fit[i] = self.function(self.X[i])  # 初始第i个粒子适应值，即个体最优值

        self.g_fit = np.max(self.p_fit)  # 初始全局最优值

    # 更新粒子位置
    def iterator(self):
        g_fits = []  # 记录全局最优值

        for t in range(self.max_iter):
            for i in range(self.pN):  # 更新gbest\pbest
                temp = self.function(self.X[i])
                if temp < self.p_fit[i]:  # 更新个体最优
                    self.p_fit[i] = temp
                    self.pbest[i] = self.X[i]

                    if temp < self.g_fit:  # 更新全局最优
                        self.g_fit = temp  # 最优值
                        self.gbest = self.X[i]  # 最优位置

            for i in range(self.pN):  # 更新公式
                self.v[i] = self.w * self.v[i] + self.c1 * self.r1 * (self.pbest[i] - self.X[i]) + \
                            self.c2 * self.r2 * (self.gbest - self.X[i])
                self.X[i] = self.X[i] + self.v[i]

            g_fits.append(self.g_fit)

        return g_fits


my_PSO = PSO(pN=30, dim=1, max_iter=100)
my_PSO.init_population()
g_fits = my_PSO.iterator()
print(g_fits)

# 绘图
t = np.arange(0, 100)
plt.plot(t, g_fits, color='b', linewidth=3)
plt.title("Figure1")
plt.xlabel("iterators", size=14)
plt.ylabel("g_fits", size=14)
plt.show()

最终收敛效果较好。