粒子群演算法（主要針對連續型函式最佳化問題）

码头牛牛發表於2024-04-02
演算法函式

文章主要參考了以下博文：

https://zhuanlan.zhihu.com/p/564819718

1. 簡介

粒子群演算法是一種解決最最佳化問題的通用方法，其優點是求解速度快，數值相對穩定，演算法簡單。粒子群演算法分為連續型粒子群演算法和離散型粒子群演算法，分別用於解決連續型問題和離散型問題。

粒子群最佳化演算法源自對鳥群捕食行為的研究：一群鳥在區域中隨機搜尋食物，搜尋的策略就是搜尋目前離食物最近的鳥的周圍區域。粒子群演算法利用這種模型得到啟示並應用於解決最佳化問題：

  • 鳥類的飛行空間： 搜尋空間，也可以理解為約束

  • 鳥（粒子）： 可行解

  • 鳥類尋找到的食物源： 最優解

在粒子群演算法中，每個最佳化問題的潛在解都是搜尋空間中的一隻鳥，稱之為粒子（在連續性例子群演算法中，每一個粒子就代表一個可行解）。所有的粒子都有一個由被最佳化的函式決定的適應度值，每個粒子還有一個速度決定它們飛翔的方向和距離。然後，粒子們就追隨當前的最優粒子在解空間中搜尋。

粒子群演算法首先在給定的解空間中隨機初始化粒子群，待最佳化問題的變數數決定了解空間的維數（這裡表示每個粒子的維度，後面我們會用D表示維度）。每個粒子有了初始位置與初始速度（x表示當前位置，v表示速度，粒子群演算法最重要的是對速度公式的更新，因為速度決定了更新的方向，進而決定解的值），然後透過迭代尋優。在每一次迭代中，每個粒子透過跟蹤兩個“極值”來更新自己在解空間中的空間位置與飛行速度：一個極值就是單個粒子本身在迭代過程中找到的最優解粒子，這個粒子叫作個體極值；（這裡的意思是，每個粒子具有記憶性，每個粒子子記得兩個值，一個是全部的粒子距離目標的最優值，另一個值是這個粒子曾經歷史的最優值） 另一個極值是種群所有粒子在迭代過程中所找到的最優解粒子，這個粒子是全域性極值。

2. 連續型粒子群最佳化演算法

2.1 基本原理

假設在一個\(D維\)的目標搜尋空間中，有\(N\)個粒子組成一個群落，其中第\(i\)個粒子表示一個\(D維\)的向量：

\[X_i=(x_{i1},x_{i2},\dots,x_{iD}),\quad{i}=1,2,\dots,N\quad\quad(1) \]

第\(i\)個粒子的“飛行”速度也是一個\(D維\)向量：

\[V_i=(v_{i1},v_{i2},\dots,v_{iD}),\quad{i}=1,2,\dots,N\quad\quad\quad(2) \]

第\(i\)個粒子迄今為止搜到的最優位置稱為個體極值（比如：\((p_{11},p_{12},\dots,p_{1D})\)就是第一個粒子搜尋到的最優位置），記為：

\[p_{best}=(p_{i1},p_{i2},\dots,p_{iD}),\quad{i}=1,2,\dots,N\quad\quad(3) \]

整個粒子群迄今為止搜尋到的最優位置為全域性極值（我的理解是\(g_{best}\)就是本次迭代的區域性最優解），記為：

\[g_{best}=(g_{1},g_{2},\dots,g_{D})\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(4) \]

在找到這兩個最優值時，粒子根據如下的式(5)和式(6)分別來更新自己的速度和位置:

\[v_{ij}(t+1)=v_{ij}(t)+c_1r_1[p_{ij}(t)-x_{ij}(t)]+c_2r_2[p_{gj}(t)-x_{ij}(t)]\quad(5) \]

\[x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)\quad\quad\quad\quad\quad\quad\quad\quad(6) \]

式(5)和式(6)是整個粒子群演算法最核心的部分，對於式(5)，其中：

  • \(c_1\)和\(c_2\)為學習因子，也稱加速常數；

  • \(r_1\)和\(r_2\)為\([0，1]\)範圍內的均勻隨機數；\(r_1\)和\(r_2\)是介於0和1之間的隨機數，增加了粒子飛行的隨機性。

  • \(j=1，2，…，D\)；

  • \(v_{ij}\)是粒子的速度，\(v_{ij}∈［-v_{max}，v_{max}］\)，\(v_{max}\)是常數，由使用者設定來限制粒子的速度。

  • 式（5）右邊由三部分組成：

    • 第一部分為“慣性”或“動量”部分，反映了粒子的運動“習慣”，代表粒子有維持自己先前速度的趨勢（就是原有的速度）；

    • 第二部分為“認知”部分，反映了粒子對自身歷史經驗的記憶或回憶，代表粒子有向自身歷史最佳位置逼近的趨勢；

    • 第三部分為“社會”部分，反映了粒子間協同合作與知識共享的群體歷史經驗，代表粒子有向群體或鄰域歷史最佳位置逼近的趨勢（對於\(p_{gj}\)，我的理解是把\(p_{best}\)的最優值中的\(i\)下標的值記為\(g\)，或者換句話來說，\(p_{gj}(t)\)就是整個群體的歷史最佳位置）。

上面的是基本粒子群演算法，下面的標準粒子群演算法的更新是這樣的，就是在速度\(v_{ij}\)的前面加一個權重係數\(w\)。使得權重係數是一個動態調整的。保證在迭代前期保證足夠的全域性搜尋能力，迭代後期能夠專注於區域性最優解的搜尋能力。

\[v_{ij}(t+1)=w·v_{ij}(t)+c_1r_1[p_{ij}(t)-x_{ij}(t)]+c_2r_2[p_{gj}(t)-x_{ij}(t)]\quad(7) \]

\[x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)\quad\quad\quad\quad\quad\quad\quad\quad(8) \]

權重更新公式如下（這是其中一種，還有其他的更新方式）：

\[w=w_{max}-\frac{(w_{max}-w_{min})·t}{T_{max}}\quad\quad\quad\quad\quad\quad\quad\quad(9) \]

2.2 演算法流程

（1）初始化粒子群，包括群體規模\(N\)，每個粒子的位置\(x_i\)和速度\(v_i\)。

（2）計算每個粒子的適應度值\(fitness_i\)。

（3）對每個粒子，用它的適應度值\(fitness_i\)和個體極值\(p_{best}(i)\)比較。如果\(fitness_i＜p_{best}(i)\)，則用\(fitness_i\)替換掉\(p_{best}(i)\)。

（4）對每個粒子，用它的適應度值\(fitness_i\)和全域性極值\(g_{best}\)比較。如果\(fitness_i＜g_{best}\)，則用\(fitness_i\)替換\(g_{best}\)。

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

（6）進行邊界條件處理。

（7）判斷演算法終止條件是否滿足：若是，則結束演算法並輸出最佳化結果；否則返回步驟（2）。

2.3 演算法詳解

求函式\(f（x，y）=3cos（xy）+x+y^2\)的最小值，其中\(x\)的取值範圍為\(［-4，4］\)，y的取值範圍為\(［-4，4］\)。為了方便程式設計，在演算法中，\(x\)用\(x[0]\)表示，\(y\)用\(x[1]\)來表示。

2.3.1 繪圖

首先我們來看看\(f（x，y）=3cos（xy）+x+y^2\)函式長啥樣：

# 求函式f(x,y) = 3*cos(x * y) + x + y**2的最小值，其中-4 <= x <= 4, -4 <= y <= 4

# 首先繪製這個函式的三維影像
import numpy as np
import matplotlib.pyplot as plt

X = np.arange(-4 ,4 ,0.01)
Y = np.arange(-4 ,4 ,0.01)
x, y = np.meshgrid(X ,Y)
Z = 3*np.cos(x * y) + x + y**2

# 作圖
fig = plt.figure(figsize=(10,15))
ax3 = plt.axes(projection = "3d")
ax3.plot_surface(x,y,Z ,cmap = "rainbow")
# ax3.contour(x ,y ,Z ,zdim = "z" ,offset=-2 ,cmap = "rainbow")
plt.show()

2.3.2 連續型粒子群演算法求解

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 設定字型和設定負號
matplotlib.rc("font", family="KaiTi")
matplotlib.rcParams["axes.unicode_minus"] = False
# 初始化種群，群體規模，每個粒子的速度和規模
N = 100 # 種群數目
D = 2 # 維度
T = 200 # 最大迭代次數
c1 = c2 = 1.5 # 個體學習因子與群體學習因子
w_max = 0.8 # 權重係數最大值
w_min = 0.4 # 權重係數最小值
x_max = 4 # 每個維度最大取值範圍，如果每個維度不一樣，那麼可以寫一個陣列，下面程式碼依次需要改變
x_min = -4 # 同上
v_max = 1 # 每個維度粒子的最大速度
v_min = -1 # 每個維度粒子的最小速度


# 定義適應度函式
def func(x):
    return 3 * np.cos(x[0] * x[1]) + x[0] + x[1] ** 2


# 初始化種群個體
x = np.random.rand(N, D) * (x_max - x_min) + x_min # 初始化每個粒子的位置
v = np.random.rand(N, D) * (v_max - v_min) + v_min # 初始化每個粒子的速度

# 初始化個體最優位置和最優值
p = x # 用來儲存每一個粒子的歷史最優位置
p_best = np.ones((N, 1))  # 每行儲存的是最優值
for i in range(N): # 初始化每個粒子的最優值，此時就是把位置帶進去，把適應度值計算出來
    p_best[i] = func(x[i, :])

# 初始化全域性最優位置和全域性最優值
g_best = 100 #設定真的全域性最優值
gb = np.ones(T) # 用於記錄每一次迭代的全域性最優值
x_best = np.ones(D) # 用於儲存最優粒子的取值

# 按照公式依次迭代直到滿足精度或者迭代次數
for i in range(T):
    for j in range(N):
        # 個更新個體最優值和全域性最優值
        if p_best[j] > func(x[j,:]):
            p_best[j] = func(x[j,:])
            p[j,:] = x[j,:].copy()
        # p_best[j] = func(x[j,:]) if func(x[j,:]) < p_best[j] else p_best[j]
        # 更新全域性最優值
        if g_best > p_best[j]:
            g_best = p_best[j]
            x_best = x[j,:].copy()   # 一定要加copy，否則後面x[j,:]更新也會將x_best更新
        # 計算動態慣性權重
        w = w_max - (w_max - w_min) * i / T
        # 更新位置和速度
        v[j, :] = w * v[j, :] + c1 * np.random.rand(1) * (p[j, :] - x[j, :]) + c2 * np.random.rand(1) * (x_best - x[j, :])
        x[j, :] = x[j, :] + v[j, :]
        # 邊界條件處理
        for ii in range(D):
            if (v[j, ii] > v_max) or (v[j, ii] < v_min):
                v[j, ii] = v_min + np.random.rand(1) * (v_max - v_min)
            if (x[j, ii] > x_max) or (x[j, ii] < x_min):
                x[j, ii] = x_min + np.random.rand(1) * (x_max - x_min)
    # 記錄歷代全域性最優值
    gb[i] = g_best
print("最優值為", gb[T - 1],"最優位置為",x_best)
plt.plot(range(T),gb)
plt.xlabel("迭代次數")
plt.ylabel("適應度值")
plt.title("適應度進化曲線")
plt.show()

相關文章