【演化(进化)算法】粒子群算法

文章目录

• 粒子群算法
• • 1 粒子群算法概述
  • 2 相关变量
  • 3 固定的参数
  • 4 粒子群算法求解优化问题
  • 5 实例
  • 6 python实现
  • 7 特点

粒子群算法

1 粒子群算法概述

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

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

2 相关变量

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

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

X i k = ( x i 1 k , x i 2 k , … , x i D k ) X_{i}^k=(x_{i1}^k,x_{i2}^k, \ldots ,x_{iD}^k) Xik​=(xi1k​,xi2k​,…,xiDk​)

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

P i k = ( p i 1 k , p i 2 k , … , p i D k ) P_{i}^k=(p_{i1}^k,p_{i2}^k, \ldots ,p_{iD}^k) Pik​=(pi1k​,pi2k​,…,piDk​)

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

P b e s t k = ( p b e s t 1 k , p b e s t 2 k , … , p b e s t D k ) P_{best}^k=(p_{best1}^k,p_{best2}^k, \ldots ,p_{bestD}^k) Pbestk​=(pbest1k​,pbest2k​,…,pbestDk​)

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

V i k = ( v i 1 k , v i 2 k , … , v i D k ) V_{i}^k=(v_{i1}^k,v_{i2}^k, \ldots ,v_{iD}^k) Vik​=(vi1k​,vi2k​,…,viDk​)

V i k + 1 = w ⋅ V i k + c 1 r 1 ( P i k − X i k ) + c 2 r 2 ( P b e s t k − X i k ) V^{k+1}_{i}=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}) Vik+1​=w⋅Vik​+c1​r1​(Pik​−Xik​)+c2​r2​(Pbestk​−Xik​)

其中 w w w为惯性权重； r 1 r_1 r1​和 r 2 r_2 r2​为0到1的随机数， c 1 c_1 c1​为局部学习因子， c 2 c_2 c2​为全局学习因子。

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

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

f i f_i fi​——个体历史最优适应值

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

f b e s t f_{best} fbest​——群体历史最优适应值

3 固定的参数

PSO算法中固定的参数有：

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

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

•　学习因子： c 1 c_1 c1​为局部学习因子， c 2 c_2 c2​为全局学习因子， c 1 c_1 c1​和 c 2 c_2 c2​通常可设定为2.0。

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

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

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

优化问题建模：

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

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

粒子群算法流程：

（3）设置粒子数 N N N，最大速度 V m a x {V_{max}} Vmax​，局部学习因子 c 1 c_1 c1​，全局学习因子 c 2 c_2 c2​，惯性权重 w w w，最大进化代数 G G G等。确定适应度函数。

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

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

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

V i k + 1 = w ( k ) ⋅ V i k + c 1 r 1 ( P i k − X i k ) + c 2 r 2 ( P b e s t k − X i k ) V^{k+1}_{i}=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}) Vik+1​=w(k)⋅Vik​+c1​r1​(Pik​−Xik​)+c2​r2​(Pbestk​−Xik​)

位置更新公式为：

X i k + 1 = X i k + V i k + 1 X_{i}^{k+1}=X_{i}^k+V_i^{k+1} Xik+1​=Xik​+Vik+1​

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

（8）更新群体全局最优：比较粒子当前适应值 f i K = f i f_i^K=f_i fiK​=fi​与群体最优值 f b e s t f_{best} fbest​，如果 f i K f_i^K fiK​优于 f b e s t f_{best} fbest​，则置全局最优(位置) P b e s t P_{best} Pbest​为当前值 X i K X_i^K XiK​，全局最优适应值 f b e s t f_{best} fbest​为当前值 f i K f_i^K fiK​。

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

5 实例

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

F ( x , y ) = 3 ( − x + 1 ) 2 e − x 2 − ( y + 1 ) 2 − ( − 10 x 3 + 2 x − 10 y 5 ) e − x 2 − y 2 − 3 − e − y 2 − ( x + 1 ) 2 F(x,y)=3 \left(- x + 1\right)^{2} e^{- x^{2} - \left(y + 1\right)^{2}} - \left(- 10 x^{3} + 2 x - 10 y^{5}\right) e^{- x^{2} - y^{2}} - 3^{- e^{- y^{2} - \left(x + 1\right)^{2}}} F(x,y)=3(−x+1)2e−x2−(y+1)2−(−10x3+2x−10y5)e−x2−y2−3−e−y2−(x+1)2

（2）确定目标函数的变量参数以及各种约束条件，即确定出遗传算法中的个体和问题的解空间，得到优化模型。

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

得到优化模型：

max ⁡ F ( x , y ) = 3 ( − x + 1 ) 2 e − x 2 − ( y + 1 ) 2 − ( − 10 x 3 + 2 x − 10 y 5 ) e − x 2 − y 2 − 3 − e − y 2 − ( x + 1 ) 2 \max F(x,y)=3 \left(- x + 1\right)^{2} e^{- x^{2} - \left(y + 1\right)^{2}} - \left(- 10 x^{3} + 2 x - 10 y^{5}\right) e^{- x^{2} - y^{2}} - 3^{- e^{- y^{2} - \left(x + 1\right)^{2}}} maxF(x,y)=3(−x+1)2e−x2−(y+1)2−(−10x3+2x−10y5)e−x2−y2−3−e−y2−(x+1)2

s.t. x , y ∈ [ − 3 ， 3 ] \text{s.t.}\qquad x,y\in[-3，3] s.t.x,y∈[−3，3]

待优化模型可视化如下：

粒子群算法流程：

（3）设置粒子数 N = 40 N=40 N=40，最大速度 V m a x = 3 {V_{max}}=3 Vmax​=3，局部学习因子 c 1 = 2 c_1=2 c1​=2，全局学习因子 c 2 = 2 c_2=2 c2​=2，惯性权重 w = 0.8 w=0.8 w=0.8，最大进化代数 G = 100 G=100 G=100。

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

F i t n e s s = F ( x , y ) Fitness=F(x,y) Fitness=F(x,y)

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

F i t n e s s = − F ( x , y ) Fitness=-F(x,y) Fitness=−F(x,y)

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

各个粒子的初始位置坐标取 [ − 3 , 3 ] [-3,3] [−3,3]的随机值，速度取 [ − V m a x , V m a x ] [-V_{max},V_{max}] [−Vmax​,Vmax​]的随机数

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

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

速度更新公式为：

V i k + 1 = w ( k ) ⋅ V i k + c 1 r 1 ( P i k − X i k ) + c 2 r 2 ( P b e s t k − X i k ) V^{k+1}_{i}=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}) Vik+1​=w(k)⋅Vik​+c1​r1​(Pik​−Xik​)+c2​r2​(Pbestk​−Xik​)

位置更新公式为：

X i k + 1 = X i k + V i k + 1 X_{i}^{k+1}=X_{i}^k+V_i^{k+1} Xik+1​=Xik​+Vik+1​

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

（8）更新群体全局最优：比较粒子当前适应值 f i K = f i f_i^K=f_i fiK​=fi​与群体最优值 f b e s t f_{best} fbest​，如果 f i K f_i^K fiK​优于 f b e s t f_{best} fbest​，则置全局最优(位置) P b e s t P_{best} Pbest​为当前值 X i K X_i^K XiK​，全局最优适应值 f b e s t f_{best} fbest​为当前值 f i K f_i^K fiK​。

（9）检查是否达到最大迭代代数。若满足，则结束寻优；否则 k = k ＋ 1 k=k＋1 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()
           

求最大：

求最小：

7 特点

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

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

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

e i π + 1 = 0 e^{i\pi}+1=0 eiπ+1=0