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

文章目錄

• 粒子群算法
• • 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