粒子群算法 PSO【1】- 基本原理

算法簡介

粒子群算法，Particle Swarm optimization，簡稱PSO，是由Eberhart博士和kennedy博士發明的一種啟發式算法，

其思想來源於 群體間互相協作和信息共享，使得群體行為形成從無序到有序的演化，進而達到目的 的自然現象；

通過模擬鳥群覓食行為而發展起來的一種基於群體協作的隨機搜索算法。

通常認為它是群集智能 (Swarm intelligence, SI) 的一種。

 

算法原理

形象描述

粒子群算法源於鳥群覓食行為；　　【或者蟲兒覓食，蜂群覓食，都差不多】

假設有一群鳥（種群），在某片森林（解空間）尋找食物 [的位置]（最優解）；

起初，鳥兒都不知道食物在哪里，各自在隨意飛翔（種群初始化），

但他們有 一種超能力據說是嗅覺 能判斷自己離食物的距離，並且能記錄自己在飛行過程中離食物最近時的位置（個人經驗）；　　【位置對應解，距離代表適應度/目標函數】

同時，鳥群有個特異功能，類似有個微信群，能夠得知其他鳥經過的離食物最近的位置（群體經驗）；

鳥兒為了得到食物必須做出決策：在當前位置和速度狀態下，下一步應該往哪飛呢？（更新解）

實際上，鳥兒會結合自己的經驗和群體經驗做出最終決策；　　【決策是下一步的速度，包括大小和方向，更新解更新的是位置】

然后重復以上過程，即可找到食物

 

數學描述

粒子群算法涉及兩個最重要的公式

　　速度更新公式　　

　　　　位置更新公式

i 代表粒子編號，d 代表解維度編號；

v 代表速度，x 代表位置；

w 代表慣性因子，即當前飛行速度的慣性，

c1，c2為經驗權重，c1也叫個體學習因子，c2 也叫社會學習因子，r1，r2代表折扣因子，一般取值0-1；

α 代表約束因子，控制速度的權重；

pi 代表粒子 i 的個體最優經驗；

pg 代表群體最優經驗；

 

補充說明

部分內容屬於摘抄，不一定對，僅供參考

1. 越界處理：為了限制某個粒子的飛翔范圍過大，我們為粒子的每一個維度設置一個最大飛翔速度Vmax，

如果粒子i在d維上的飛翔速度大於Vmax，則將這一維上的飛翔速度設置為Vmax，如果小於-Vmax,則將其設置為-Vmax。

Vmax的取值一般為對應維度所代表的決策變量取值范圍的百分之十到百分之二十。即 0.2*(Vmax-Vmin)。 

2. 慣性因子：用來控制繼承多少粒子當前的速度的，越大則對於當前速度的繼承程度越小，越小則對於當前速度的繼承程度越大。有些同學可能會產生疑問，是不是說反了。其實不是，從公式中可以明確看出，其值越大，則速度的改變幅度就越大，則對於粒子的當前速度繼承越小；反之，速度的改變幅度越小，則對於粒子當前速度繼承越大。

因此值越大，則解的搜索范圍越大，可以提高算法的全局搜索能力，但也損失了局部搜索能力，有可能錯失最優解；

反之值越小，則解的搜索范圍也就越小，算法的全局搜索能力也就越小，容易陷入局部最優。

如果是變量，則其值應該隨着迭代次數的增加而減小（類似於梯度下降當中的學習率）。如果為定值，則建議在0.6到0.75之間進行選取。

3. 加速常數c1、c2：加速常數控制着速度的計算是更加看重自身經驗還是群體經驗。取值在學術界分歧很大主要有如下幾種情況：

學者	c1和c2
Clerc	c1=c2=2.05
Carlisle	c1=2.8,c2=1.3
Trelea	w=0.6,c1=c2=1.7
Eberhart	w=0.729,c1=c2=1.494

 

 

 

 

 

 

算法步驟

注意對迭代過程中不可行解的處理　　　　【可參考我的遺傳算法相關博客】

1. 懲罰

2. 人工修復

 

示例代碼

'''
以一個簡單的例子來介紹算法的流程。
假設我們需要尋找函數 y=x1^2+x2^2+x3^3+x4^4 在[1,30]之間的最大值。
我們很容易就知道，當x1=x2=x3=x4=30時，該函數能取到最大值。
'''

import random
import numpy as np
import matplotlib.pyplot as plt


class PSO:
    def __init__(self, NGEN, popsize, low, up):
        # 初始化
        self.NGEN = NGEN  # 迭代的代數
        self.pop_size = popsize  # 種群大小
        self.var_num = len(low)  # 變量維度
        self.bound = [low, up]

        # 初始化第0代初始全局最優解
        # self.pop_x = np.random.uniform(low, up, size=(self.pop_size, self.var_num))  # 所有粒子的位置
        # self.pop_x = np.random.randint(low, up, size=(self.pop_size, self.var_num))
        self.pop_x = np.random.randint(12, 21, size=(self.pop_size, self.var_num))
        print(self.pop_x)
        self.pop_v = np.random.rand(self.pop_size, self.var_num)  # 所有粒子的速度
        # print(self.pop_v)
        self.p_best = self.pop_x  # 每個粒子最優的位置
        temp = -1
        for i in range(self.pop_size):
            fit = self.fitness(self.p_best[i])
            if fit > temp:
                self.g_best = self.p_best[i]
                temp = fit

    def fitness(self, ind_var):
        """
        個體適應值計算
        """
        x1, x2, x3, x4 = ind_var
        y = x1 ** 2 + x2 ** 2 + x3 ** 3 + x4 ** 4
        return y

    def update_operator(self, pop_size):
        """
        更新算子：更新下一時刻的位置和速度
        """
        c1 = 0.8  # 學習因子，一般為2
        c2 = 0.99
        w = 0.4  # 自身權重因子
        for i in range(pop_size):
            # 更新速度
            self.pop_v[i] = w * self.pop_v[i] + \
                            c1 * random.uniform(0, 1) * (self.p_best[i] - self.pop_x[i]) + \
                            c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i])
            # 更新位置
            self.pop_x[i] = self.pop_x[i] + self.pop_v[i]

            # 越界保護  -- 不可行解的人工修復
            for j in range(self.var_num):
                if self.pop_x[i][j] < self.bound[0][j]:
                    self.pop_x[i][j] = self.bound[0][j]
                if self.pop_x[i][j] > self.bound[1][j]:
                    self.pop_x[i][j] = self.bound[1][j]

            # 更新p_best和g_best
            fit_x = self.fitness(self.pop_x[i])
            if fit_x > self.fitness(self.p_best[i]):  # 個體最優
                self.p_best[i] = self.pop_x[i]
            if fit_x > self.fitness(self.g_best):     # 群體最優
                self.g_best = self.pop_x[i]

    def main(self):
        popobj = []             # 記錄全局最優粒子的適應度
        self.ng_best = np.zeros((1, self.var_num))[0]     # 記錄全局最優粒子
        print(self.ng_best)
        for gen in range(self.NGEN):
            self.update_operator(self.pop_size)
            popobj.append(self.fitness(self.g_best))
            print('############ Generation {} ############'.format(str(gen + 1)))
            if self.fitness(self.g_best) > self.fitness(self.ng_best):
                self.ng_best = self.g_best.copy()
            print('最好的位置：{}'.format(self.ng_best))
            print('最大的函數值：{}'.format(self.fitness(self.ng_best)))
        print("---- End of (successful) Searching ----")

        plt.figure()
        plt.title("Figure1")
        plt.xlabel("iterators", size=14)
        plt.ylabel("fitness", size=14)
        plt.plot(range(self.NGEN), popobj, color='b', linewidth=2)
        plt.show()


if __name__ == '__main__':
    NGEN = 500
    popsize = 100
    low = [1, 10, 1, 1]
    up = [30, 50, 30, 30]
    pso = PSO(NGEN, popsize, low, up)
    pso.main()

 

 

 

 

參考資料：

https://blog.csdn.net/yy2050645/article/details/80740641　　粒子群算法

https://www.zhihu.com/question/23103725/answer/365298309　　如何直觀形象地理解粒子群算法？

https://www.cnblogs.com/yncxzdy/p/4280207.html　　標准粒子群算法（PSO）

 

https://finthon.com/python-pso/　　Python手把手構建粒子群算法（PSO）實現最優化搜索

https://www.jb51.net/article/128113.htm　　Python編程實現粒子群算法(PSO)詳解