【演化(进化)算法】粒子群算法及其Python实现_演化算法相关公式

粒子群算法

1 粒子群算法概述

粒子群算法，也称粒子群优化（Particle SwarmOptimization，PSO）算法，源于对鸟群捕食的行为研究。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。那么，找到食物的最优策略就是搜寻目前离食物最近的鸟的周围区域。粒子群算法算法从这种模型中得到启示并用于解决优化问题。

PSO算法中，优化问题的每个解向量都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的目标函数确定的适应度值，适应度值越大越好。每个粒子还有一个速度向量决定它们飞行的方向和距离，粒子们追随当前的最优粒子在解空间中搜索。PSO算法首先将鸟群初始化为一群随机粒子（随机解向量），然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个“最优（适应度最大）位置向量”来更新自己的位置：第一个最优向量是每个粒子本身所找到的历史最优解，这个解向量叫作个体最优。另一个最优向量是整个群体找到的历史最优解向量，这个解向量称为全局最优。

2 相关变量

PSO算法中，假设在$D$维搜索空间中（目标函数的变量个数为$D$），有 $N$个粒子，每个粒子代表一个解，对于$i=1,2,...,N$，有以下变量：

第$k$次迭代，第$i$个粒子的位置向量（目标函数的变量组成的向量）为：

$$X_i^k=(x_{i1}^k,x_{i2}^k,\dots ,x_{iD}^k)$$

第$k$次迭代，第$i$个粒子搜索到的最优位置向量（个体最优解）为：

$$P_i^k=(p_{i1}^k,p_{i2}^k,\dots ,p_{iD}^k)$$

第$k$次迭代，整个群体搜索到的最优位置向量（全局最优解）为：

$$P_{best}^k=(p_{best1}^k,p_{best2}^k,\dots ,p_{bestD}^k)$$

第$k$次迭代，第$i$个粒子的速度向量（粒子移动的方向和大小）为：

$$V_i^k=(v_{i1}^k,v_{i2}^k,\dots ,v_{iD}^k)$$
$$V_i^{k+1}=w\cdot V_i^k+c_1{r}_1(P_i^k-X_i^k)+c_2{r}_2(P_{best}^k-X_i^k)$$

其中$w$为惯性权重；$r_1$和$r_2$为0到1的随机数，$c_1$为局部学习因子，$c_2$为全局学习因子。

例如，在二维空间中，各个向量表示如下：

第$k$次迭代，第$i$个粒子搜索到的最优位置的适应值（优化目标函数的值）为：

$f_i$——个体历史最优适应值

第$k$次迭代， 群体搜索到的最优位置的适应值为：

$f_{best}$——群体历史最优适应值

---

3 固定的参数

PSO算法中固定的参数有：

•　粒子数$N$：一般取20~40，对于比较难的问题，粒子数可以取到100或200。

•　最大速度$V_{max}$：决定粒子在一个迭代中最大的移动距离。如果是较大的$V_{max}$，可以保证粒子群体的全局搜索能力，如果是较小的$V_{max}$，则粒子群体的局部搜索能力加强。

•　学习因子：$c_1$为局部学习因子，$c_2$为全局学习因子，$c_1$和$c_2$通常可设定为2.0。

•　惯性权重$w$ ：一个大的惯性权重有利于展开全局寻优，而一个小的惯性权值有利于局部寻优。当粒子的最大速度$V_{max}$很小时，应使用接近于1的惯性权重。可选择使用递减权重。

•　终止条件：最大进化代数$G$或最小误差要求。

4 粒子群算法求解优化问题

优化问题建模：

（1）确定待优化的目标函数。

（2）确定目标函数的变量参数以及各种约束条件。

粒子群算法流程：

（3）设置粒子数$N$，最大速度$V_{max}$，局部学习因子$c_1$，全局学习因子$c_2$，惯性权重$w$，最大进化代数$G$等。确定适应度函数。

（4）初始化粒子群各个粒子的位置、速度。

（5）将各个粒子初始位置作为个体最优，计算群体中各个粒子的初始适应值$f_i^0$，并求出群体最优位置。

（6）更新粒子的速度和位置，产生新群体，并对粒子的速度和位置进行越界检查，速度向量公式为：

$$V_i^{k+1}=w(k)\cdot V_i^k+c_1{r}_1(P_i^k-X_i^k)+c_2{r}_2(P_{best}^k-X_i^k)$$

位置更新公式为：

$$X_i^{k+1}=X_i^k+V_i^{k+1}$$

（7）更新个体最优：比较粒子的当前适应值$f_i^K$和自身历史最优值$f_i$，如果$f_i^K$优于$f_i$，则置个体最优(位置)$P_i$为当前值$X_i^K$，个体最优适应值$f_i$为当前值$f_i^K$。

（8）更新群体全局最优：比较粒子当前适应值$f_i^K=f_i$与群体最优值$f_{best}$，如果$f_i^K$优于$f_{best}$，则置全局最优(位置)$P_{best}$为当前值$X_i^K$，全局最优适应值$f_{best}$为当前值$f_i^K$。

（9）检查结束条件，若满足，则结束寻优；否则$k=k＋1$，转至（6）。结束条件为寻优达到最大进化代数，或评价值小于给定精度。

5 实例

（1）确定待优化的目标函数。

$$F(x,y)=3{\left(-x+1\right)}^2{e}^{-x^2-{\left(y+1\right)}^2}-\left(-10x^3+2x-10y^5\right)e^{-x^2-y^2}-3^{-e^{-y^2-{\left(x+1\right)}^2}}$$
（2）确定目标函数的变量参数以及各种约束条件，即确定出遗传算法中的个体和问题的解空间，得到优化模型。

个体为$(x,y)$，约束条件为$x,y\in [-3，3]$。

得到优化模型：

$$\max F(x,y)=3{\left(-x+1\right)}^2{e}^{-x^2-{\left(y+1\right)}^2}-\left(-10x^3+2x-10y^5\right)e^{-x^2-y^2}-3^{-e^{-y^2-{\left(x+1\right)}^2}}$$
$$\text{s.t.}x,y\in [-3，3]$$

待优化模型可视化如下：

粒子群算法流程：

（3）设置粒子数$N=40$，最大速度$V_{max}=3$，局部学习因子$c_1=2$，全局学习因子$c_2=2$，惯性权重$w=0.8$，最大进化代数$G=100$。

然后确定适应度函数：粒子群算法没有选择操作，所以可以直接以目标函数作为适应度函数：

$$Fitness=F(x,y)$$

如果我们要求目标函数的最小值，仅需要取上式的负数：

$$Fitness=-F(x,y)$$

（4）初始化粒子群各个粒子的位置、速度。

各个粒子的初始位置坐标取$[-3,3]$的随机值，速度取$[-V_{max},V_{max}]$的随机数

（5）将各个粒子初始位置作为个体最优，计算群体中各个粒子的初始适应值$f_i^0$，并求出群体最优位置。

（6）更新粒子的速度和位置，产生新群体，并对粒子的速度和位置进行越界检查。速度大于$V_{max}$则置为$V_{max}$，小于$-V_{max}$则置为$-V_{max}$；同样地，位置坐标超过$-3$则置为$-3$，超过$3$则置为$3$。

速度更新公式为：

$$V_i^{k+1}=w(k)\cdot V_i^k+c_1{r}_1(P_i^k-X_i^k)+c_2{r}_2(P_{best}^k-X_i^k)$$

位置更新公式为：
$$X_i^{k+1}=X_i^k+V_i^{k+1}$$

（7）更新个体最优：比较粒子的当前适应值$f_i^K$和自身历史最优值$f_i$，如果$f_i^K$优于$f_i$，则置个体最优(位置)$P_i$为当前值$X_i^K$，个体最优适应值$f_i$为当前值$f_i^K$。

（8）更新群体全局最优：比较粒子当前适应值$f_i^K=f_i$与群体最优值$f_{best}$，如果$f_i^K$优于$f_{best}$，则置全局最优(位置)$P_{best}$为当前值$X_i^K$，全局最优适应值$f_{best}$为当前值$f_i^K$。

（9）检查是否达到最大迭代代数。若满足，则结束寻优；否则$k=k＋1$，转至（6）

6 python实现

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


def F(X):
    F = 3 * (1 - X[0]) ** 2 * np.exp(-(X[0] ** 2) - (X[1] + 1) ** 2) - 10 * (
            X[0] / 5 - X[0] ** 3 - X[1] ** 5) * np.exp(
        -X[0] ** 2 - X[1] ** 2) - 1 / 3 ** np.exp(-(X[0] + 1) ** 2 - X[1] ** 2)
    return F


def fitness(X):
    fit_value = F(X)  # 求最大
    #fit_value = -F(X)  # 求最小
    return fit_value


def plot_3d(ax):
    X = np.linspace(-3, 3, 100)
    Y = np.linspace(-3, 3, 100)
    X, Y = np.meshgrid(X, Y)
    Z = F([X, Y])
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm)
    ax.set_zlim(-10, 10)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    plt.pause(3)
    plt.show()


def main():
    fig1 = plt.figure()
    ax = Axes3D(fig1)
    plt.ion()  # 将画图模式改为交互模式，程序遇到plt.show不会暂停，而是继续执行
    plot_3d(ax)
    # 参数初始化
    p_num = 40         # 粒子数
    max_iter = 100
    w = 0.8           # 惯性权重
    c1 = 2            # 局部学习因子
    c2 = 2            # 全局学习因子
    dim = 2           # 搜索维度
    X = np.zeros((p_num, dim))  # 位置
    V = np.zeros((p_num, dim))  # 速度
    p = np.zeros((p_num, dim))  # 个体最优(位置)
    p_best = np.zeros((1, dim))  # 全局最优(位置)
    f = np.zeros(p_num)  # 个体最优
    f_best = -np.inf     # 全局最优
    V_max = 0.2 
    # 初始化
    for i in range(p_num):
        for j in range(dim):
            X[i][j] = np.random.random(1) * 6 - 3
            V[i][j] = np.random.random(1) * V_max * 2 - V_max
        p[i] = X[i]
        temp = fitness(X[i])  # 个体最优为初始位置
        f[i] = temp
        if temp > f_best:  # 全局最优为适应度最大的个体的位置
            f_best = temp  
            p_best = X[i]
    # 进入粒子群迭代
    for t in range(max_iter):
        F_value = []  # 储存目标函数值，用于可视化
        for i in range(p_num):
            r1 = np.random.random(1)
            r2 = np.random.random(1)
            # 更新速度
            V[i] = w * V[i] + c1 * r1 * (p[i] - X[i]) + c2 * r2 * (p_best - X[i])
            # 更新位置
            X[i] = X[i] + V[i]
            # 速度限制
            V[i][V[i] > V_max] = V_max
            V[i][V[i] < -V_max] = -V_max
            # 位置限制
            X[i][X[i] > 3] = 3
            X[i][X[i] < -3] = -3
            F_value.append(F(X[i]))
        # 更新个体最优和全局最优
        for i in range(p_num):
            temp = fitness(X[i])
            if temp > f[i]:  # 更新个体最优
                f[i] = temp
                p[i] = X[i]
                if temp > f_best:  # 更新全局最优
                    p_best = X[i]
                    f_best = temp

        if 'sca' in locals():
            sca.remove()  # 去除图像中上一个种群的点
        sca = ax.scatter(X[:, 0], X[:, 1], F_value, c='black', marker='o')
        plt.show()
        plt.pause(0.1)
    plt.ioff()
    plot_3d(ax)
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

求最大：

求最小：

7 特点

可以看到，和遗传算法相似，PSO算法也是从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质，但它比遗传算法规则更为简单，没有遗传算法的交叉和变异操作，通过追随当前搜索到的最优值来寻找全局最优。这种算法具有实现容易、精度高、收敛快等优点。

[1] 《智能控制：理论基础、算法设计与应用作者》 作者：刘金琨

十分感谢三连支持！最靓的公式送给最靓的你：

$$e^{i\pi }+1=0$$