粒子群優(yōu)化算法

一、概述

粒子群優(yōu)化算法（Particle Swarm Optimization，PSO）的思想來源于對鳥捕食行為的模仿，最初，Reynolds.Heppner 等科學(xué)家研究的是鳥類飛行的美學(xué)和那些能使鳥群同時突然改變方向，分散，聚集的定律上，這些都依賴于鳥的努力來維持群體中個體間最佳距離來實現(xiàn)同步。而社會生物學(xué)家 E.O.Wilson 參考魚群的社會行為認(rèn)為從理論上說，在搜尋食物的過程中，盡管食物的分配不可知，群中的個體可以從群中其它個體的發(fā)現(xiàn)以及以往的經(jīng)驗中獲益。

粒子群從這種模型中得到啟發(fā)并用于解決優(yōu)化問題。如果我們把一個優(yōu)化問題看作是在空中覓食的鳥群，那么粒子群中每個優(yōu)化問題的潛在解都是搜索空間的一只鳥，稱之為“粒子”（Particle），“食物”就是優(yōu)化問題的最優(yōu)解。每個粒子都有一個由優(yōu)化問題決定的適應(yīng)度用來評價粒子的“好壞”程度，每個粒子還有一個速度決定它們飛翔的方向和距離，它根據(jù)自己的飛行經(jīng)驗和同伴的飛行經(jīng)驗來調(diào)整自己的飛行。粒子群初始化為一群隨機粒子（隨機解），然后通過迭代的方式尋找最優(yōu)解，在每一次的迭代中，粒子通過跟蹤兩個“極值”來更新自己，第一個是粒子本身所經(jīng)歷過的最好位置，稱為個體極值即；另一個是整個群體經(jīng)歷過的最好位置稱為全局極值。每個粒子通過上述的兩個極值不斷更新自己，從而產(chǎn)生新一代的群體。

二、粒子群算法

算法的描述如下：

假設(shè)搜索空間是維，并且群體中有個粒子。那么群體中的第個粒子可以表示為一個維的向量，，即第個粒子在維的搜索空間的位置是，它所經(jīng)歷的“最好”位置記作。粒子的每個位置代表要求的一個潛在解，把它代入目標(biāo)函數(shù)就可以得到它的適應(yīng)度值，用來評判粒子的“好壞”程度。整個群體迄今為止搜索到的最優(yōu)位置記作，是最優(yōu)粒子位置的索引。

為慣性權(quán)重（inertia weight），為第個粒子到第代為止搜索到的歷史最優(yōu)解，為整個粒子群到目前為止搜索到的最優(yōu)解，，分別是第個粒子當(dāng)前的位置和飛行速度，為非負(fù)的常數(shù)，稱為加速度因子，是之間的隨機數(shù)。

公式由三部分組成，第一部分是粒子當(dāng)前的速度，表明了粒子當(dāng)前的狀態(tài)；第二部分是認(rèn)知部分（Cognition Modal），表示粒子本身的思考（也稱為自身認(rèn)知系數(shù)）；第三部分是社會認(rèn)知部分（Social Modal），表示粒子間的信息共享（為社會認(rèn)知系數(shù)）。

參數(shù)的選擇：

粒子數(shù)目一般取30~50，參數(shù)一般取2。適應(yīng)度函數(shù)、粒子的維數(shù)和取值范圍要視具體問題而定。問題解的編碼方式通?？梢圆捎脤崝?shù)編碼。

算法的主要步驟如下：

第一步：對粒子群的隨機位置和速度進行初始設(shè)定，同時設(shè)定迭代次數(shù)。

第二步：計算每個粒子的適應(yīng)度值。

第三步：對每個粒子，將其適應(yīng)度值與所經(jīng)歷的最好位置的適應(yīng)度值進行比較，若較好，則將其作為當(dāng)前的個體最優(yōu)位置。

第四步：對每個粒子，將其適應(yīng)度值與全局所經(jīng)歷的最好位置的適應(yīng)度值進行比較，若較好，則將其作為當(dāng)前的全局最優(yōu)位置。

第五步：根據(jù)公式(1),(2)對粒子的速度和位置進行優(yōu)化，從而更新粒子位置。

第六步：如未達到結(jié)束條件（通常為最大循環(huán)數(shù)或最小誤差要求），則返回第二步。

三、基于粒子群算法的非線性函數(shù)尋優(yōu)

本案例尋優(yōu)的非線性函數(shù)為

當(dāng)，時，該函數(shù)為Ackley函數(shù)，函數(shù)圖形如圖1所示。

從函數(shù)圖形可以看出，該函數(shù)有很多局部極小值，在處取到全局最小值0。

本案例群體的粒子數(shù)為20，每個粒子的維數(shù)為2，算法迭代進化次數(shù)為100。加速度因子，，慣性權(quán)重。

引入的模塊如下：

import numpy as np
import matplotlib.pyplot as plt

# 解決中文亂碼和負(fù)號問題
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

適應(yīng)度函數(shù)代碼如下：

# 函數(shù)用于計算粒子適應(yīng)度值
def funAckley(x):
    c = 20
    y = -1.0 * c * np.exp(-0.2 * np.sqrt((x[0] ** 2 + x[1] ** 2) / 2)) \
        - np.exp((np.cos(2 * np.pi * x[0]) + np.cos(2 * np.pi * x[1])) / 2) \
        + c + np.exp(1)
    return y

PSO算法代碼如下：

c1, c2 = 1.4, 1.5  # 加速度因子
maxgen = 100  # 進化次數(shù)
sizepop = 20  # 群體規(guī)模
w = 0.1  # 慣性權(quán)重
Vmax, Vmin = 1, -1  # 速度最大值，最小值
popmax, popmin = 5, -5  # 個體最大值，最小值
pop = np.zeros([sizepop, 2])  # 種群
V = np.zeros([sizepop, 2])  # 速度
fitness = np.zeros(sizepop)  # 適應(yīng)度
trace = np.zeros(maxgen)  # 結(jié)果

# 隨機產(chǎn)生一個群體,初始粒子和速度
for i in range(sizepop):
    pop[i] = 5 * ((-2 * np.random.rand(1, 2) + 1))
    V[i] = (-2 * np.random.rand(1, 2) + 1)
    fitness[i] = funAckley(pop[i])

# 個體極值和群體極值
bestfitness = np.min(fitness)
bestindex = np.argmin(fitness)
Gbest = pop[bestindex]  # 全局最佳
fitnessGbest = bestfitness  # 全局最佳適應(yīng)度值
Pbest = pop.copy()  # 個體最佳
fitnessPbest = fitness.copy()  # 個體最佳適應(yīng)度值

# 迭代尋優(yōu)
for i in range(maxgen):
    for j in range(sizepop):
        V[j] = w * V[j] + c1 * (-2 * np.random.rand(1, 2) + 1) * (Pbest[j] - pop[j]) \
               + c2 * (-2 * np.random.rand(1, 2) + 1) * (Gbest - pop[j])  # 速度更新
        V[j, V[j] > Vmax] = Vmax
        V[j, V[j] < Vmin] = Vmin
        pop[j] = pop[j] + V[j]  # 群體更新
        pop[j, pop[j] > popmax] = popmax
        pop[j, pop[j] < popmin] = popmin
        fitness[j] = funAckley(pop[j])  # 適應(yīng)度值

        if fitness[j] < fitnessPbest[j]:
            Pbest[j] = pop[j].copy()
            fitnessPbest[j] = fitness[j]
        if fitness[j] < fitnessGbest:
            Gbest = pop[j].copy()
            fitnessGbest = fitness[j]

    trace[i] = fitnessGbest
    print(Gbest, fitnessGbest)

# 結(jié)果分析
plt.plot(trace)
plt.title("最優(yōu)個體適應(yīng)度")
plt.xlabel("進化代數(shù)")
plt.ylabel("適應(yīng)度")
plt.show()

算法結(jié)果如下：

最終得到的個體適應(yīng)度值為，對應(yīng)的粒子位置為，PSO算法尋優(yōu)得到最優(yōu)值接近函數(shù)實際最優(yōu)值，說明PSO算法具有較強的函數(shù)極值尋優(yōu)能力。

四、基于自適應(yīng)變異粒子群算法的非線性函數(shù)尋優(yōu)

本案例尋優(yōu)的非線性函數(shù)（Shubert函數(shù)）為

該函數(shù)圖形如圖3所示：

從函數(shù)圖形可以看出，該函數(shù)有很多局部極小值，很難用傳統(tǒng)的梯度下降方法進行全局尋優(yōu)。

自適應(yīng)變異是借鑒遺傳算法中的變異思想，在PSO算法中引入變異操作，即對某些變量以一定的概率重新初始化。變異操作拓展了在迭代中不斷縮小的種群搜索空間，使粒子能夠跳出先前搜索到的最優(yōu)值位置，在更大的空間中開展搜索，同時保持了種群多樣性，提高算法尋找更優(yōu)值的可能性。因此，在普通粒子群算法的基礎(chǔ)上引入簡單變異算子，在粒子每次更新之后，以一定概率重新初始化粒子。

本案例群體的粒子數(shù)為50，每個粒子的維數(shù)為2，算法迭代進化次數(shù)為1000。加速度因子，，慣性權(quán)重。

適應(yīng)度函數(shù)代碼如下：

def funShubert(x):
    h1 = 0
    h2 = 0
    for i in range(1, 6):
        h1 += i * np.cos((i + 1) * x[0] + i)
        h2 += i * np.cos((i + 1) * x[1] + i)
    return h1 * h2

自適應(yīng)變異PSO算法代碼如下：

c1, c2 = 1.4, 1.5  # 加速度因子
maxgen = 1000  # 進化次數(shù)
sizepop = 50  # 群體規(guī)模
w = 0.8  # 慣性權(quán)重
Vmax, Vmin = 5, -5  # 速度最大值，最小值
popmax, popmin = 10, -10  # 個體最大值，最小值
pop = np.zeros([sizepop, 2])  # 種群
V = np.zeros([sizepop, 2])  # 速度
fitness = np.zeros(sizepop)  # 適應(yīng)度
trace = np.zeros(maxgen)  # 結(jié)果

# 隨機產(chǎn)生一個群體,初始粒子和速度
for i in range(sizepop):
    pop[i] = 10 * (-2 * np.random.rand(1, 2) + 1)
    V[i] = 5 * (-2 * np.random.rand(1, 2) + 1)
    fitness[i] = funShubert(pop[i])

# 個體極值和群體極值
bestfitness = np.min(fitness)
bestindex = np.argmin(fitness)
Gbest = pop[bestindex]  # 全局最佳
fitnessGbest = bestfitness  # 全局最佳適應(yīng)度值
Pbest = pop.copy()  # 個體最佳
fitnessPbest = fitness.copy()  # 個體最佳適應(yīng)度值

# 迭代尋優(yōu)
for i in range(maxgen):
    for j in range(sizepop):
        V[j] = w * V[j] + c1 * (-2 * np.random.rand(1, 2) + 1) * (Pbest[j] - pop[j]) \
               + c2 * (-2 * np.random.rand(1, 2) + 1) * (Gbest - pop[j])  # 速度更新
        V[j, V[j] > Vmax] = Vmax
        V[j, V[j] < Vmin] = Vmin
        pop[j] = pop[j] + V[j]  # 群體更新
        pop[j, pop[j] > popmax] = popmax
        pop[j, pop[j] < popmin] = popmin
        if np.random.rand() > 0.9: # 自適應(yīng)變異
            pop[j] = 10 * (-2 * np.random.rand(1, 2) + 1)
        fitness[j] = funShubert(pop[j])  # 適應(yīng)度值
        
        if fitness[j] < fitnessPbest[j]:
            Pbest[j] = pop[j].copy()
            fitnessPbest[j] = fitness[j]
        if fitness[j] < fitnessGbest:
            Gbest = pop[j].copy()
            fitnessGbest = fitness[j]

    trace[i] = fitnessGbest
    print(Gbest, fitnessGbest)

# 結(jié)果分析
plt.plot(trace)
plt.title("最優(yōu)個體適應(yīng)度")
plt.xlabel("進化代數(shù)")
plt.ylabel("適應(yīng)度")
plt.show()

算法結(jié)果如下：

最終得到的個體適應(yīng)度值為-186.6831831611777，對應(yīng)的例子位置為，自適應(yīng)變異PSO算法尋優(yōu)得到最優(yōu)值接近函數(shù)實際最優(yōu)值，說明該算法具有較強的函數(shù)極值尋優(yōu)能力。

五、補充

慣性權(quán)重體現(xiàn)的是粒子當(dāng)前速度多大程度上繼承先前的速度，Shi.Y最先將慣性權(quán)重引入到PSO算法中，并分析指出一個較大的慣性權(quán)值有利于全局搜索，而一個較小的慣性權(quán)值則更有利于局部搜索。為了更好地平衡算法的全局搜索與局部搜索能力，其提出了線性遞減慣性權(quán)重（Linear Decreasing Inertia Weight，LDIW），即

式中，為初始慣性權(quán)重，為迭代至最大次數(shù)時的慣性權(quán)重，為當(dāng)前迭代次數(shù)，為最大迭代次數(shù)。

一般來說，慣性權(quán)重取值為，時算法性能最好。這樣，隨著迭代的進行，慣性權(quán)重由0.9線性遞減至0.4，迭代初期較大的慣性權(quán)重使算法保持了較強的全局搜索能力，而迭代后期較小的慣性權(quán)值有利于算法進行更精確的局部搜索。

線性慣性權(quán)重只是一種經(jīng)驗做法，常用的慣性權(quán)重的選擇還包括以下幾種：

六、練習(xí)

求測試函數(shù)的最小值，以及最小值點。

1. Rastrigin function：

當(dāng)時，如下圖所示：

2. Sphere function：

當(dāng)時，如下圖所示：

3. Beale function：

4. Booth function：

5. Bukin function：

6. three-hump camel function：

7. H?lder table function：