Python代碼實現粒子群算法圖文詳解
1.引言
粒子群優化算法起源於對鳥群覓食活動的分析。鳥群在覓食的時候通常會毫無征兆的聚攏，分散，以及改變飛行的軌跡，但是在不同個體之間會十分默契的保持距離。所以粒子群優化算法模擬鳥類覓食的過程，將待求解問題的搜索空間看作是鳥類飛行的空間，將每隻鳥抽象成一個沒有質量和大小的粒子，用這個粒子來表示待求解問題的一個可行解。所以，尋找最優解的過程就相當於鳥類覓食的過程。
粒子群算法也是基於種群以及進化的概念，通過個體間的競爭與協作，實現復雜空間最優解的求解。但是與遺傳算法不同的是，他不會對每個個體進行“交叉”，“變異”等操作，而實以一定的規則，更新每個粒子的速度以及位置，使得每一個粒子向自身歷史最佳位置以及全局歷史最佳位置進行移動，從而實現整個種群向著最優的方向進化。
2.算法的具體描述：
2.1原理
在粒子群優化算法中，粒子之間通過信息共享機制，獲得其它粒子的發現與飛行經歷。粒子群算法中的信息共享機制實際上是一種合作共生的行為，在搜索最優解的過程中，每個粒子能夠對自己經過的最佳的歷史位置進行記憶，同時，每個粒子的行為有會受到群體中其他例子的影響，所以在搜索最優解的過程中，粒子的行為既受其他粒子的影響，有受到自身經驗的指導。
粒子群優化算法對於鳥群的模擬是按照如下的模式進行的：假設一群鳥在空中搜索食物，所有鳥知道自己當前距離食物有多遠（這裡的遠近會用一個值來衡量，適應度值），那麼每隻鳥最簡單的搜索策略就是尋找距離目前距離食物最近的鳥的周圍空間。因此，在粒子群算法中，每個粒子都相當於一隻鳥，每個粒子有一個適應度值，還有一個速度決定他們的飛行的距離與方向。所有的粒子追隨當前最優的粒子在解空間中搜索。每搜索一次，最優的粒子會發生變化，其他的粒子又會追隨新的最優粒子進行搜索，如此反復迭代。
在迭代開始的時候，每個粒子通過隨機的方式初始化在空間中的速度和位置，然後在迭代過程中，粒子通過跟蹤兩個極值來自己在解空間中的位置和速度，一個極值是單個粒子自身在迭代的過程中的最優位置（就是最優適應度值所對應的空間解），這個稱之為粒子的個體極值。另一個極值是種群中所有的粒子在迭代過程中所找到的最優位置，這個成為全局極值。如果粒子隻是跟蹤一個極值的話，則算法稱為局部粒子群算法或者全局粒子群算法。
PSO從這種模型中得到啟示並用於解決優化問題。PSO 中，每個優化問題的潛在解都是搜索空間中的一隻鳥，稱之為粒子。所有的粒子都有一個由被優化的函數決定的適值( fitness value) ，每個粒子還有一個速度決定它們飛翔的方向和距離。然後粒子們就追隨當前的最優粒子在解空間中搜索。
PSO初始化為一群隨機粒子(隨機解)，然後通過迭代找到最優解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自己；第一個就是粒子本身所找到的最優解，這個解稱為個體極值；另一個極值是整個種群目前找到的最優解，這個極值是全局極值。另外也可以不用整個種群而隻是用其中一部分作為粒子的鄰居，那麼在所有鄰居中的極值就是局部極值。
圖解：
2.2標準粒子群算法流程
算法的流程如下：
Step1:種群初始化:可以進行隨機初始化或根據被優化的問題設計特定的初始化方法,包括群體規模，每個粒子的位置 X i X_{i} Xi 和速度 V i V_i Vi ，然後計算每個粒子的適應度值，從而選擇出個體的局部最優位置向量和種群的全局最優位置向量。
Step2:迭代設置:設置迭代次數 g m a x g_{max} gmax ，令當前迭代次數g=1。
Step3:根據公式更新每個粒子的速度向量V。
Step4:根據公式更新每個粒子的位置向量X。
Step5:局部位置向量和全局位置向量更新:更新每個粒子的Pbest,和種群的Gbest。
Step6:終止條件判斷:判斷迭代次數時都達到 g m a x g_{max} gmax 或誤差是否足夠小，如果滿足則輸出Gbest.否則繼續進行迭代，跳轉至步驟(3)。
對於粒子群優化算法的實際應用，因為主要是對速度和位置向量迭代算子的設計計,選代算子是否合理將決定整個PSO算法性能的優劣.，所以如何設計 t pso的迭代算子是算法應用的研究重點和難點。
3.代碼案例
3.1問題
求解f(x,y)的最小值點
3.2繪圖
import numpy as np import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D # 生成X和Y的數據 X=np.arange(-5,5,0.1) Y=np.arange(-5,5,0.1) X,Y=np.meshgrid(X,Y) # 目標函數 Z=X**2+Y**2+X # 繪圖 fig=plt.figure() ax=Axes3D(fig) surf=ax.plot_surface(X,Y,Z,cmap=cm.coolwarm) plt.show()
3.3計算適應度
# 速度
# Vi+1 = w*Vi + c1 * r1 * (pbest_i - Xi) + c2 * r2 * (gbest_i - Xi)
# 位置
# Xi+1 = Xi + Vi+1
# vi, xi 分別表示粒子第i維的速度和位置
# pbest_i, gbest_i 分別表示某個粒子最好位置第i維的值、整個種群最好位置第i維的值
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認字體
mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題
def fitness_func(X):
"""計算粒子的的適應度值，也就是目標函數值，X 的維度是 size * 2 """
A = 10
pi = np.pi
x = X[:, 0]
y = X[:, 1]
return x**2+y**2+x
3.4更新速度
def velocity_update(V, X, pbest, gbest, c1, c2, w, max_val):
"""
根據速度更新公式更新每個粒子的速度
:param V: 粒子當前的速度矩陣，20*2 的矩陣
:param X: 粒子當前的位置矩陣，20*2 的矩陣
:param pbest: 每個粒子歷史最優位置，20*2 的矩陣
:param gbest: 種群歷史最優位置，1*2 的矩陣
"""
size = X.shape[0]
r1 = np.random.random((size, 1))
r2 = np.random.random((size, 1))
V = w*V+c1*r1*(pbest-X)+c2*r2*(gbest-X)
# 防止越界處理
V[V < -max_val] = -max_val
V[V > -max_val] = max_val
return V
3.5更新粒子位置
def position_update(X, V):
"""
根據公式更新粒子的位置
:param X: 粒子當前的位置矩陣，維度是 20*2
:param V: 粒子當前的速度舉著，維度是 20*2
"""
return X+V
3.6主要算法過程
def pos():
w = 1
c1 = 2
c2 = 2
r1 = None
r2 = None
dim = 2
size = 20
iter_num = 1000
max_val = 0.5
best_fitness = float(9e10)
fitness_val_list = []
# 初始化種群各個粒子的位置
X = np.random.uniform(-5, 5, size=(size, dim))
# 初始化各個粒子的速度
V = np.random.uniform(-0.5, 0.5, size=(size, dim))
# print(X)
p_fitness = fitness_func(X)
g_fitness = p_fitness.min()
fitness_val_list.append(g_fitness)
# 初始化的個體最優位置和種群最優位置
pbest = X
gbest = X[p_fitness.argmin()]
# 迭代計算
for i in range(1, iter_num):
V = velocity_update(V, X, pbest, gbest, c1, c2, w, max_val)
X = position_update(X, V)
p_fitness2 = fitness_func(X)
g_fitness2 = p_fitness2.min()
# 更新每個粒子的歷史最優位置
for j in range(size):
if p_fitness[j] > p_fitness2[j]:
pbest[j] = X[j]
p_fitness[j] = p_fitness2[j]
# 更新群體的最優位置
if g_fitness > g_fitness2:
gbest = X[p_fitness2.argmin()]
g_fitness = g_fitness2
# 記錄最優迭代記錄
fitness_val_list.append(g_fitness)
i += 1
# 輸出迭代結果
print("最優值是：%.5f" % fitness_val_list[-1])
print("最優解是：x=%.5f,y=%.5f" % (gbest[0], gbest[1]))
# 繪圖
plt.plot(fitness_val_list, color='r')
plt.title('迭代過程')
plt.show()
pos()
結果
最優值是：-0.23696
最優解是：x=-0.54359,y=-0.10555
參考：
蘇振裕.《Python最優化實戰》[M].北京大學出版社
總結
本篇文章就到這裡瞭，希望能給你帶來幫助，也希望您能夠多多關註WalkonNet的更多內容！