粒子群优化算法(PSO)python实践
1 算法介紹和原理
1.1 算法原理
強烈推薦知乎大佬的這篇文章:粒子群優化算法(Particle Swarm Optimization, PSO)的詳細解讀 - 知乎 (zhihu.com)。該文章詳細介紹了算法的原理、算法流程、參數解釋和一些Tips,這里就不過多贅述了。
粒子群優化算法(PSO, Particle Swarm Optimization),屬于啟發式算法中的一種,常用于多目標優化,尋找全局最優解,具有收斂速度快、參數少、算法簡單的優點。
算法流程圖如下(圖片來自這篇文章):
1.2 更新公式
1.2.1 速度更新公式
vidk+1=ωvidk+c1r1(pid,pbest?k?xidk)+c2r2(pd,gbest?k?xidk)v_{i d}^{k+1}=\omega v_{i d}^k+c_1 r_1\left(p_{i d, \text { pbest }}^k-x_{i d}^k\right)+c_2 r_2\left(p_{d, \text { gbest }}^k-x_{i d}^k\right) vidk+1?=ωvidk?+c1?r1?(pid,?pbest?k??xidk?)+c2?r2?(pd,?gbest?k??xidk?)
vidk+1v_{i d}^{k+1}vidk+1? —— 粒子 iii 在第 kkk 次迭代中第 ddd 維的速度向量。
pid,pbest?kp_{i d, \text { pbest }}^kpid,?pbest?k? —— 粒子 iii 在第 kkk 次迭代中第 ddd 維的歷史最優位置。
速度可以看作一個向量,具有大小和方向。即是粒子下一輪迭代移動的距離和方向。公式分為三部分,第一部分為慣性項,由該粒子的當前速度和慣性權重 ω\omegaω 組成。第二部分為認知項,即是粒子當前位置和自身歷史最優位置間的距離和方向。 第三部分為社會項,即是粒子當前位置和群體歷史最優位置間的距離和方向。
對于更新速度的方向,等于三部分向量和向量的方向。
1.2.2 位置更新公式
xidk+1=xidk+vidk+1x_{i d}^{k+1}=x_{i d}^{k}+v_{i d}^{k+1} xidk+1?=xidk?+vidk+1?
點加向量等于點
大致掌握算法原理后,直接上手代碼。
2 代碼實現
示例問題:
求解如下函數的極小值
y=x1ex2+x3sinx2+x4x5y=x_1e^{x_2}+x_3sinx_2+x_4x_5 y=x1?ex2?+x3?sinx2?+x4?x5?
每個變量的取值都在(1,25)。
首先是定義一個求解類及其初始化方法。
class PSO:def __init__(self, D, N, M, p_low, p_up, v_low, v_high, w = 1., c1 = 2., c2 = 2.):self.w = w # 慣性權值self.c1 = c1 # 個體學習因子self.c2 = c2 # 群體學習因子self.D = D # 粒子維度self.N = N # 粒子群規模,初始化種群個數self.M = M # 最大迭代次數self.p_range = [p_low, p_up] # 粒子位置的約束范圍self.v_range = [v_low, v_high] # 粒子速度的約束范圍self.x = np.zeros((self.N, self.D)) # 所有粒子的位置self.v = np.zeros((self.N, self.D)) # 所有粒子的速度self.p_best = np.zeros((self.N, self.D)) # 每個粒子的最優位置self.g_best = np.zeros((1, self.D))[0] # 種群(全局)的最優位置self.p_bestFit = np.zeros(self.N) # 每個粒子的最優適應值self.g_bestFit = float('Inf') # float('-Inf'),始化種群(全局)的最優適應值,由于求極小值,故初始值給大,向下收斂,這里默認優化問題中只有一個全局最優解# 初始化所有個體和全局信息for i in range(self.N):for j in range(self.D):self.x[i][j] = random.uniform(self.p_range[0][j], self.p_range[1][j])self.v[i][j] = random.uniform(self.v_range[0], self.v_range[1])self.p_best[i] = self.x[i] # 保存個體歷史最優位置,初始默認第0代為最優fit = self.fitness(self.p_best[i])self.p_bestFit[i] = fit # 保存個體歷史最優適應值if fit < self.g_bestFit: # 尋找并保存全局最優位置和適應值self.g_best = self.p_best[i]self.g_bestFit = fit然后定義適應度計算函數,也就是我們要尋優的對象。
def fitness(x):"""根據粒子位置計算適應值,可根據問題情況自定義"""return x[0] * np.exp(x[1]) + x[2] * np.sin(x[1]) + x[3] * x[4]定義每次迭代的更新函數。
def update(self):for i in range(self.N):# 更新速度(核心公式)self.v[i] = self.w * self.v[i] + self.c1 * random.uniform(0, 1) * (self.p_best[i] - self.x[i]) + self.c2 * random.uniform(0, 1) * (self.g_best - self.x[i])# 速度限制for j in range(self.D):if self.v[i][j] < self.v_range[0]:self.v[i][j] = self.v_range[0]if self.v[i][j] > self.v_range[1]:self.v[i][j] = self.v_range[1]# 更新位置self.x[i] = self.x[i] + self.v[i]# 位置限制for j in range(self.D):if self.x[i][j] < self.p_range[0][j]:self.x[i][j] = self.p_range[0][j]if self.x[i][j] > self.p_range[1][j]:self.x[i][j] = self.p_range[1][j]# 更新個體和全局歷史最優位置及適應值_fit = self.fitness(self.x[i])if _fit < self.p_bestFit[i]:self.p_best[i] = self.x[i]self.p_bestFit[i] = _fitif _fit < self.g_bestFit:self.g_best = self.x[i]self.g_bestFit = _fit其中主要完成每輪迭代中單個粒子位置和速度,歷史最優位置和最優適應度的更新,以及群體(全局)的最優位置和最優適應度的更新。
最后,便是主要函數的實現。
def pso(self, draw = 1):best_fit = [] # 記錄每輪迭代的最佳適應度,用于繪圖w_range = Noneif isinstance(self.w, tuple):w_range = self.w[1] - self.w[0]self.w = self.w[1]time_start = time.time() # 記錄迭代尋優開始時間for i in range(self.M):self.update() # 更新主要參數和信息if w_range:self.w -= w_range / self.M # 慣性權重線性遞減print("\rIter: {:d}/{:d} fitness: {:.4f} ".format(i, self.M, self.g_bestFit, end = '\n'))best_fit.append(self.g_bestFit.copy())time_end = time.time() # 記錄迭代尋優結束時間print(f'Algorithm takes {time_end - time_start} seconds') # 打印算法總運行時間,單位為秒/sif draw:plt.figure()plt.plot([i for i in range(self.M)], best_fit)plt.xlabel("iter")plt.ylabel("fitness")plt.title("Iter process")plt.show()測試代碼如下。
if __name__ == '__main__':low = [1, 1, 1, 1, 1]up = [25, 25, 25, 25, 25]pso = PSO(5, 100, 50, low, up, -1, 1, w = 0.9)pso.pso()測試結果如下圖所示。
... Iter: 47/50 fitness: 4.5598 Iter: 48/50 fitness: 4.5598 Iter: 49/50 fitness: 4.5598 Algorithm takes 0.1444549560546875 seconds可以看到在第30輪就已經完全收斂了,且函數在求解空間中的極小值為4.5598。
3 總結
-
動態的慣性權重[1]^{[1]}[1]
w_range = self.w[1] - self.w[0] self.w = self.w[1] self.w -= w_range / self.M # 慣性權重線性遞減 -
fitness變化邏輯
fitness是適應度函數值,通常問題是尋找解空間內的粒子,使得該粒子所代表的解的fitness向下或向上收斂于某一定值。對于不同收斂方向,個體和全局最優fitness一般初始化賦值無窮大或者無窮小float('Inf')/float('-Inf')。并且在判斷更新最優適應值時也應當注意大小于符號。
-
程序復用
對于上面的PSO類代碼,不同多元尋優問題均可通過重寫類中的fitness函數實現。或者定義self.fitness_function屬性進行外部函數名傳參賦值。
參考
[1] 粒子群優化算法(Particle Swarm Optimization, PSO)的詳細解讀 - 知乎 (zhihu.com)
[2] 粒子群算法(PSO)的Python實現(求解多元函數的極值)_Cyril_KI的博客-CSDN博客_pso算法python
總結
以上是生活随笔為你收集整理的粒子群优化算法(PSO)python实践的全部內容,希望文章能夠幫你解決所遇到的問題。