【智能算法】粒子羣尋優算法

時間  2019-12-09
標籤 智能算法 粒子 算法 简体版
原文   原文鏈接

1.理論基礎

粒子羣算法（particle swarm optimization，PSO）是計算智能領域中的一種生物啓發式方法，屬於羣體智能優化算法的一種，常見的羣體智能優化算法主要有以下幾類：
　　（1）蟻羣算法（Ant Colony Optimization，簡稱ACO）[1992年提出]；
　　（2）粒子羣優化算法（Particle Swarm Optimization，簡稱PSO）[1995年提出]（簡單易於實現，也是目前應用最爲普遍的羣體智能優化算法）；
　　（3）菌羣優化算法（Bacterial Foraging Optimization，簡稱BFO）[2002年提出]；
　　（4）蛙跳算法（Shuffled Frog Leading Algorithm，簡稱SFLA）[2003年提出]；
　　（5）人工蜂羣算法（Artificial Bee Colony Algorithm，簡稱ABC）[2005年提出]；
　　除了上述幾種常見的羣體智能算法之外，還有一些並非普遍應用的羣體智能算法，好比螢火蟲算法、布穀鳥算法、蝙蝠算法以及磷蝦羣算法等等。python

而其中的粒子羣優化算法（PSO）源於對鳥類捕食行爲的研究，鳥類捕食時，找到食物最簡單有限的策略就是搜尋當前距離食物最近的鳥的周圍。算法

舉個通俗的例子：dom

一羣鳥在尋找食物，在這個區域內有一個食物，全部的鳥都不知道食物在哪裏，但食物發出了香味，鳥經過香味的濃烈能判斷出本身的當前位置距離食物有多遠，同時鳥羣之間是能夠交流並告知離食物最近的鳥的位置的。試想一下這個時候會發生什麼？ide

鳥甲：哈哈哈，原來我離食物最近！函數

鳥乙、丙、丁…:我得趕忙往鳥甲那裏過去看看！學習

同時各自鳥在位置不停變化的時候，離食物的距離也不斷的變化，因此必定有過離食物最近的位置，這是他們的一個參考。優化

鳥某某：我剛纔的位置好像離食物很近了，我得往那裏再靠近點！spa

經過這樣一個過程，不斷的變化，最終鳥羣會向食物方向彙集，達到目標。在這個變化的過程當中，影響鳥的運動狀態變化的有兩個因素：3d

  • 離食物最近的鳥的位置
  • 本身以前達到過的離食物最近的位置

而鳥的每次的位置變化，除了考慮以上的兩個因素還有一個因素就是慣性。而對位置的變化，經過變化量（速度）來表示。code

經過以上的例子，咱們能夠把鳥抽象爲沒有質量和體積的微粒(點)，並延伸到N維空間，粒子I 在N維空間的位置表示爲矢量Xi＝(x1，x2，…，xN)，飛行速度表示爲矢量Vi＝(v1，v2，…，vN)．每一個粒子都有一個由目標函數決定的適應值(fitness value)，而且知道本身到目前爲止發現的最好位置(pbest)和如今的位置Xi．這個能夠看做是粒子本身的飛行經驗．除此以外，每一個粒子還知道到目前爲止整個羣體中全部粒子發現的最好位置(gbest)(gbest是pbest中的最好值)．這個能夠看做是粒子同伴的經驗．粒子就是經過本身的經驗和同伴中最好的經驗來決定下一步的運動。歸結爲下面的兩個公式：

                                                    (公式-1）

                                                                                              （公式-2）

其中，爲慣性權重因子，（什麼是權重因子，這個後面會講到，此處大概知道有這麼一個參數便可）；d=1,2，…，D；i=1,2，…，n；k爲當前迭代次數；Vid爲粒子的速度；c1和c2是非負的常數，稱爲加速因子；r1和r2是分佈於[0,1]區間的隨機數。爲了防止粒子的滿目搜索，通常建議限制位置和速度在必定的區間。

2.粒子羣尋優算法的案例與實現

假設有以下非線性函數，須要在必定範圍內找出最大值和最大值的位置，本文用Python實現。

                                                    （公式-3）

咱們先經過python把上面函數使用matplotlib把它畫處理，讓咱們有個直觀的感覺。

# -*- coding: utf-8 -*-
from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-2, 2, 0.05)
Y = np.arange(-2, 2, 0.05)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2+Y**2))/np.sqrt(X**2+Y**2)+np.exp((np.cos(2*np.pi*X)+np.cos(2*np.pi*Y))/2)-2.71289

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='hot')
plt.show()

此處使用mpl_toolkits.mplot3d畫3D圖，在代碼中，分別使用numpy中的arange建立等差長度爲0.05的數據點向量，區間爲[-2,2]之間，而後使用meshgrid將向量轉換給對應的矩陣，這樣的話，後面的函數式中就可使用np中的三角函數和數學函數了，若是不這樣的話，還得一個一個點的根據公式去作循環，拼接Z的向量點。

運行出來的圖以下：

figure_1

從函數圖形中能夠看出，這個函數有不少的局部極大值，而真正的極限位置就是在（0,0），極大值爲1.005。下面就經過粒子羣算法來找尋這個極大值和極大值的位置。根據粒子羣算法函數求極值的算法流程以下圖所示：

image

在pycham做爲python的IDE中，創建一個名爲PSO.py的文件，文件中編寫一個名PSO的類。

2.1 參數的設定

設置慣性參數爲常數1，也就是公式-1中的爲常數1；設置加速度因子爲非負的兩個常數，公式-1的中的c1和c2；設置迭代次數爲300；設置粒子羣的規模爲20；爲了防止粒子的盲目搜索，將位置和速度限制在區間[-2,2]、[-0.5,0.5]；

具體代碼以下：

def __init__(self):
        self.w = self.getweight()
        self.lr = self.getlearningrate()
        self.maxgen = self.getmaxgen()
        self.sizepop = self.getsizepop()
        self.rangepop = self.getrangepop()
        self.rangespeed = self.getrangespeed()

    def getweight(self):
        # 慣性權重
        weight = 1
        return weight

    def getlearningrate(self):
        # 分別是粒子的個體和社會的學習因子，也稱爲加速常數
        lr = (0.49445,1.49445)
        return lr

    def getmaxgen(self):
        # 最大迭代次數
        maxgen = 300
        return maxgen

    def getsizepop(self):
        # 種羣規模
        sizepop = 20
        return sizepop

    def getrangepop(self):
        # 粒子的位置的範圍限制,x、y方向的限制相同
        rangepop = (-2,2)
        return rangepop

    def getrangespeed(self):
        # 粒子的速度範圍限制
        rangespeed = (-0.5,0.5)
        return rangespeed

2.2 種羣初始化

隨機初始化粒子位置和粒子的速度，並根據適應度函數計算粒子適應度值，咱們的適應度函數就是公式-3這個函數。

首先定義一個適應度函數

def func(self,x):
        # x輸入粒子位置
        # y 粒子適應度值
        if (x[0]==0)&(x[1]==0):
            y = np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        else:
            y = np.sin(np.sqrt(x[0]**2+x[1]**2))/np.sqrt(x[0]**2+x[1]**2)+np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        return y

而後初始化粒子和計算適應度值

def initpopvfit(self,sizepop):
        pop = np.zeros((sizepop,2))
        v = np.zeros((sizepop,2))
        fitness = np.zeros(sizepop)

        for i in range(sizepop):
            pop[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            v[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            fitness[i] = self.func(pop[i])
        return pop,v,fitness

2.3 尋找初始化後的極值

def getinitbest(self,fitness,pop):
        # 羣體最優的粒子位置及其適應度值
        gbestpop,gbestfitness = pop[fitness.argmax()].copy(),fitness.max()
        #個體最優的粒子位置及其適應度值,使用copy()使得對pop的改變不影響pbestpop，pbestfitness相似
        pbestpop,pbestfitness = pop.copy(),fitness.copy()

        return gbestpop,gbestfitness,pbestpop,pbestfitness

此處須要分別找到羣體極值和個體粒子的極值。

2.4 迭代尋優

經過循環迭代，不斷的更新粒子的位置和速度，根據新粒子的適應度值更新個體和羣體的極值。這部分的代碼以下：

def run(self):
        pop, v, fitness = self.initpopvfit(self.sizepop)
        gbestpop, gbestfitness, pbestpop, pbestfitness = self.getinitbest(fitness, pop)

        result = np.zeros(self.maxgen)
        for i in range(self.maxgen):
            # 速度更新
            for j in range(self.sizepop):
                v[j] += self.lr[0] * np.random.rand() * (pbestpop[j] - pop[j]) + self.lr[1] * np.random.rand() * (
                gbestpop - pop[j])
            v[v < self.rangespeed[0]] = self.rangespeed[0]
            v[v > self.rangespeed[1]] = self.rangespeed[1]

            # 粒子位置更新
            for j in range(self.sizepop):
                pop[j] += 0.5 * v[j]
            pop[pop < self.rangepop[0]] = self.rangepop[0]
            pop[pop > self.rangepop[1]] = self.rangepop[1]

            # 適應度更新
            for j in range(self.sizepop):
                fitness[j] = self.func(pop[j])

            for j in range(self.sizepop):
                if fitness[j] > pbestfitness[j]:
                    pbestfitness[j] = fitness[j]
                    pbestpop[j] = pop[j].copy()

            if pbestfitness.max() > gbestfitness:
                gbestfitness = pbestfitness.max()
                gbestpop = pop[pbestfitness.argmax()].copy()

            result[i] = gbestfitness
        return result

2.4 結果分析

經過在主函數文件中，實例化PSO類，調用run這個函數，運行上面的代碼

from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from PSO import *

pso = PSO()
result = pso.run()
plt.figure()
plt.plot(result)
plt.show()

當迭代了300次以後，畫出每代最優個體適應度值的曲線，結果以下：

figure_2

最終獲得的最優個體適應度值爲：1.00538866933，對應的粒子位置爲：[  2.37513360e-05   3.41265823e-04]，經過比較，PSO算法能獲得最優值，接近函數實際的最優值1.005。

經過以上，證實PSO算法具備較強的函數極值尋優能力。

3.延伸與算法改進

講到這裏，咱們再來回顧一下以前粒子更新速度（公式-1）中的這個慣性權重因子，所謂慣性權重因子，體現的是粒子繼承上一次迭代速度的能力。較大的一個慣性權重因子有利於全局的搜索，而一個小的慣性權重因子則更加利於局部的搜索。爲了更好的平衡算法的全局搜索與局部搜索的能力，常給慣性權重因子改成隨迭代次數變化的一個函數，通常經常使用的慣性權重因子函數有如下幾種：

          （公式-4）

4                                   （公式-5）

2                       （公式-6）

3                 （公式-7）

 

這幾種的變化圖以下：

image

從變化圖中，能夠看到（公式-5）中，前期變化較慢，取值較大，維持了算法的全局搜索能力；後期變化較快，極大的提升了算法的局部尋優能力，經過這種權重的變化函數，能夠取得比固定不變的那種狀況要好的效果。

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程