粒子群優(yōu)化算法（PSO）python實現

　　粒子群算法屬于智慧算法的一類，與該類算法類似的還有蟻群算法，遺傳算法等。大家可以將這幾種算法進行比較。

　　粒子群優(yōu)化算法(Particle Swarm Optimization，PSO)屬于進化算法的一種，是通過模擬鳥群捕食行為設計的。從隨機解出發(fā)，通過迭代尋找最優(yōu)解，通過適應度來評價解的品質。在這里，我們舉一個例子來深入理解一下該算法：假設有一鳥群，在一座島上某個地方放有食物，但是鳥群并不知道食物在哪，只是知道與食物之間的距離。這時候鳥兒們很聰明，它們在目前距離食物最近的鳥的周圍搜索。

　　所有的粒子具有一下兩個屬性：

　　位置v

　　速度x

　　PSO初始化為一群隨機粒子(隨機解)，然后通過迭代找到最優(yōu)解。在每一次迭代種，粒子通過跟蹤兩個“極值”來跟新自己。第一個就是粒子本身所找到的最優(yōu)解pbest;另一個極值是整個種群目前找到的最優(yōu)解，及全局極值gbest。粒子通過下面的公式來更新自己速度和位置：

　　速度變換公式：

　　vi+1=w?vi+c1?rand1?(pbesti?xi)+c2?rand2?(gbesti?xi)v_{i+1}=w*v_i+c_1*rand_1*(pbest_i-x_i)+c_2*rand_2*(gbest_i-x_i)vi+1=w?vi+c1?rand1?(pbesti?xi)+c2?rand2?(gbesti?xi)

　　位置變換公式：

　　xi=xi+vi+1x_i=x_i+v_{i+1}xi=xi+vi+1

　　以上式子中：

　　w為慣性因子，一般取1

　　c1,c2c_1,c_2c1,c2為學習因子，一般取2

　　rand1,rand2rand_1,rand_2rand1,rand2為(0，1)之間的隨機數

　　vi和xiv_i和x_ivi和xi分別表示粒子第i維的速度和位置

　　pbesti,,gbestipbest_i,,gbest_ipbesti,,gbesti分別表示某個粒子最好位置第i維的值，整個種群最好位置第i維的值

　　注意：以上兩式是針對粒子的某一個維度進行跟新的。對粒子的每一維，都要用上述的式子進行更新。

　　流程圖如下：

　　# coding: utf-8

　　import numpy as np

　　import random

　　import matplotlib.pyplot as plt

　　class PSO():

　　# PSO參數設置

　　def __init__(self, pN, dim, max_iter):

　　self.w = 0.8

　　self.c1 = 2

　　self.c2 = 2

　　self.r1 = 0.6

　　self.r2 = 0.3

　　self.pN = pN # 粒子數量

　　self.dim = dim # 搜索維度

　　self.max_iter = max_iter # 迭代次數

　　self.X = np.zeros((self.pN, self.dim)) # 所有粒子的位置和速度

　　self.V = np.zeros((self.pN, self.dim))

　　self.pbest = np.zeros((self.pN, self.dim)) # 個體經歷的最佳位置和全局最佳位置

　　self.gbest = np.zeros((1, self.dim))

　　self.p_fit = np.zeros(self.pN) # 每個個體的歷史最佳適應值

　　self.fit = 1e10 # 全局最佳適應值

　　#目標函數Sphere函數

　　def function(self, X):

　　return X**4-2*X+3

　　#初始化種群

　　def init_Population(self):

　　for i in range(self.pN):

　　for j in range(self.dim):

　　self.X[i][j] = random.uniform(0, 1)

　　self.V[i][j] = random.uniform(0, 1)

　　self.pbest[i] = self.X[i]

　　tmp = self.function(self.X[i])

　　self.p_fit[i] = tmp

　　if tmp < self.fit:

　　self.fit = tmp

　　self.gbest = self.X[i]

　　# 更新粒子位置

　　def iterator(self):

　　fitness = []無錫×××醫(yī)院 https://yyk.familydoctor.com.cn/20612/

　　for t in range(self.max_iter):

　　for i in range(self.pN): # 更新gbest\pbest

　　temp = self.function(self.X[i])

　　if temp < self.p_fit[i]: # 更新個體最優(yōu)

　　self.p_fit[i] = temp

　　self.pbest[i] = self.X[i]

　　if self.p_fit[i] < self.fit: # 更新全局最優(yōu)

　　self.gbest = self.X[i]

　　self.fit = self.p_fit[i]

　　for i in range(self.pN):

　　self.V[i] = self.w * self.V[i] + self.c1 * self.r1 * (self.pbest[i] - self.X[i]) + \

　　self.c2 * self.r2 * (self.gbest - self.X[i])

　　self.X[i] = self.X[i] + self.V[i]

　　fitness.append(self.fit)

　　print(self.X[0], end=" ")

　　print(self.fit) # 輸出最優(yōu)值

　　return fitness

　　加入下列代碼驗證結果：

　　#程序

　　my_pso = PSO(pN=30, dim=1, max_iter=100)

　　my_pso.init_Population()

　　fitness = my_pso.iterator()

　　# 畫圖

　　plt.figure(1)

　　plt.title("Figure1")

　　plt.xlabel("iterators", size=14)

　　plt.ylabel("fitness", size=14)

　　t = np.array([t for t in range(0, 100)])

　　fitness = np.array(fitness)

　　plt.plot(t, fitness, color='b', linewidth=3)

　　plt.show()

　　在上面我們計算X4?2X+3X^4-2X+3X4?2X+3的優(yōu)化函數。

　　得出結果：

　　1.迭代圖像

　　2.迭代最后的結果可以看出來在x取0.1937左右，X4?2X+3X^4-2X+3X4?2X+3取得較優(yōu)值