用Python實現粒子群算法（PSO）

粒子群算法是一種基於鳥類覓食開發出來的優化算法，它是從隨機解出發，通過迭代尋找最優解，通過適應度來評價解的品質。

 

 From 《An Improved PSO Algorithm to Optimize BP Neural Network》

 

PSO算法的搜索性能取決於其全局探索和局部細化的平衡，這在很大程度上依賴於算法的控制參數，包括粒子群初始化、慣性因子w、最大飛翔速度和加速常數與等。

PSO算法具有以下優點：

不依賴於問題信息，采用實數求解，算法通用性強。

需要調整的參數少，原理簡單，容易實現，這是PSO算法的最大優點。

協同搜索，同時利用個體局部信息和群體全局信息指導搜索。

收斂速度快， 算法對計算機內存和CPU要求不高。

更容易飛越局部最優信息。對於目標函數僅能提供極少搜索最優值的信息，在其他算法無法辨別搜索方向的情況下，PSO算法的粒子具有飛越性的特點使其能夠跨過搜索平面上信息嚴重不足的障礙，飛抵全局最優目標值。比如Generalized Rosenbrock函數全局最小值在原占附近.但是此函數全局最優值與可到達的局部最優值之間右一條獨長的山路，曲面山谷中點的最速下降方向幾乎與到函數最小值的最佳方向垂直，找到全局最小值的可能性微乎其微， 但是PSO算法完全有可能找到全局最優值。

同時， PSO算法的缺點也是顯而易見的：

算法局部搜索能力較差，搜索精度不夠高。

算法不能絕對保證搜索到全局最優解。

 

PSO算法設計的具體步驟如下：

•  初始化粒子群（速度和位置）、慣性因子、加速常數、最大迭代次數、算法終止的最小允許誤差。
•  評價每個粒子的初始適應值。
•  將初始適應值作為當前每個粒子的局部最優值，並將各適應值對應的位置作為每個粒子的局部最優值所在的位置。
•  將最佳初始適應值作為當前全局最優值，並將最佳適應值對應的位置作為全局最優值所在的位置。
•  依據公式更新每個粒子當前的飛翔速度。
•  對每個粒子的飛翔速度進行限幅處理，使之不能超過設定的最大飛翔速度。
•  依據公式更新每個粒子當前所在的位置。
•  比較當前每個粒子的適應值是否比歷史局部最優值好，如果好，則將當前粒子適應值作為粒子的局部最優值，其對應的位置作為每個粒子的局部最優值所在的位置。
•  在當前群中找出全局最優值，並將當前全局最優值對應的位置作為粒子群的全局最優值所在的位置。
•  重復步驟（5）~（9），直到滿足設定的最小誤差或最大迭代次數
•  輸出粒子群的全局最優值和其對應的位置以及每個粒子的局部最優值和其對應的位置。

 

本文中我們假設要求解一個維度為10的向量，這里的適應度函數采用簡單的線性誤差求和。

#基本粒子群算法
#vi+1 = w*vi+c1*r1*(pi-xi)+c2*r2*(pg-xi)   速度更新公式
#xi+1 = xi + a*vi+1    位置更新公式（一般a=1）
#w = wmax -(wmax-wmin)*iter/Iter  權重更新公式
#iter當前迭代次數 Iter最大迭代次數 c1、c2學習因子  r1、r2隨機數 pi粒子當前最優位置  pg粒子群全局最優
#初始化 wmax=0.9 wmin=0.4 通常c1=c2=2 Iter對於小規模問題（10，20）對於大規模（100，200）
#算法優劣取決於w、c1和c2，迭代結束的條件是適應度函數的值符合具體問題的要求
#初始化粒子群，包括尺寸、速度和位置
#本算法假設想要的輸出是長度為10的矩陣，y=[1.7]*10,適應度函數f（x）= |x-y| <=0.001符合要求

import numpy as np

swarmsize = 500
partlen = 10
wmax,wmin = 0.9,0.4
c1 = c2 = 2
Iter = 400

def getwgh(iter):
    w = wmax - (wmax-wmin)*iter/Iter
    return w

def getrange():
    randompv = (np.random.rand()-0.5)*2
    return randompv

def initswarm():
    vswarm,pswarm = np.zeros((swarmsize,partlen)),np.zeros((swarmsize,partlen))
    for i in range(swarmsize):
        for j in range(partlen):
            vswarm[i][j] = getrange()
            pswarm[i][j] = getrange()
    return vswarm,pswarm
    
def getfitness(pswarm):
    pbest = np.zeros(partlen)
    fitness = np.zeros(swarmsize)
    for i in range(partlen):
        pbest[i] = 1.7
    
    for i in range(swarmsize):
        yloss = pswarm[i] - pbest
        for j in range(partlen):
            fitness[i] += abs(yloss[j])
    return fitness

def getpgfit(fitness,pswarm):
    pgfitness = fitness.min()
    pg = pswarm[fitness.argmin()].copy()
    return pg,pgfitness

vswarm,pswarm = initswarm()
fitness = getfitness(pswarm)
pg,pgfit = getpgfit(fitness,pswarm)
pi,pifit = pswarm.copy(),fitness.copy()

for iter in range(Iter):
    if pgfit <= 0.001:
        break
    #更新速度和位置
    weight = getwgh(iter)
    for i in range(swarmsize):
        for j in range(partlen):
            vswarm[i][j] = weight*vswarm[i][j] + c1*np.random.rand()*(pi[i][j]-pswarm[i][j]) + c2*np.random.rand()*(pg[j]-pswarm[i][j])
            pswarm[i][j] = pswarm[i][j] + vswarm[i][j]
    #更新適應值
    fitness = getfitness(pswarm)
    #更新全局最優粒子
    pg,pgfit = getpgfit(fitness,pswarm)
    #更新局部最優粒子
    for i in range(swarmsize):
        if fitness[i] < pifit[i]:
            pifit[i] = fitness[i].copy()
            pi[i] = pswarm[i].copy()

for j in range(swarmsize):
    if pifit[j] < pgfit:
        pgfit = pifit[j].copy()
        pg = pi[j].copy()
print(pg)
print(pgfit)

下面的結果分別是迭代300次和400次的結果。

 

 可以看到400次迭代雖然適應度沒有達到預期，得到的向量已經很接近期望的結果了。

 

寫在最后：粒子群算法最重要的參數就是慣性權重和學習因子，針對這兩個參數有了新的優化粒子群算法（IPSO）。還有初始化粒子群時速度和位置范圍的確定，包括種群的大小和迭代次數的選擇，這些都是‘摸着石頭過河’，沒有標准答案。