一：前言

粒子群算法是由Kennedy和Eberhart在1995年提出的一種基于群體智能的優化算法。該算法模擬鳥群覓食行為，通過粒子間的信息共享和協作來尋找最優解。

二：算法原理

2.1 核心思想

粒子群算法的核心思想可以通過鳥群找食物的故事?理解，

假設一群小鳥在森林里找食物（最優解）：
?（1）每只鳥不知道食物在哪，但能感知自己當前位置離食物有多遠（通過目標函數計算適應度）。

（2）?小鳥們會交流?：每只鳥記住自己飛過的最好位置（個體經驗），同時整個鳥群知道所有鳥中最好的位置（群體經驗）。

?（3）飛行策略?：每只鳥決定下一步飛哪里時，會綜合三個方向：

• 慣性方向?：保持原來的飛行方向（不想急轉彎）。
• ?個體經驗方向?：飛向自己曾經找到過食物的地方。
• ?群體經驗方向?：飛向鳥群中其他鳥找到食物最多的地方。

?（4）動態調整?：小鳥們一邊飛一邊更新信息，最終整個鳥群會逐漸聚集到食物最豐富的位置

2.2 PSO核心公式?

2.3 PSO算法流程圖

三：python實現 二維Rastrigin函數 最低點檢索例子

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
import time# 設置matplotlib支持中文顯示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = Falsedef fitness_function(position):"""目標函數：二維Rastrigin函數"""x, y = positionreturn 20 + (x**2 - 10 * np.cos(2 * np.pi * x)) + (y**2 - 10 * np.cos(2 * np.pi * y))def pso_optimization(func, bounds, n_particles=30, max_iter=50, w=0.7, c1=1.5, c2=1.5):"""粒子群優化算法實現"""dim = len(bounds)# 初始化粒子群particles = {'position': np.array([np.random.uniform(low, high, n_particles) for (low, high) in bounds]).T,'velocity': np.zeros((n_particles, dim)),'best_position': np.zeros((n_particles, dim)),'best_fitness': np.full(n_particles, np.inf),'history': []  # 存儲迭代歷史}# 初始全局最優global_best_pos = np.zeros(dim)global_best_fit = np.inf# 計算初始適應度for i in range(n_particles):fitness = func(particles['position'][i])particles['best_position'][i] = particles['position'][i].copy()particles['best_fitness'][i] = fitnessif fitness < global_best_fit:global_best_fit = fitnessglobal_best_pos = particles['position'][i].copy()particles['history'].append({'positions': particles['position'].copy(),'global_best': global_best_pos.copy()})# 主循環for t in range(max_iter):for i in range(n_particles):# 更新速度和位置r1, r2 = np.random.rand(), np.random.rand()particles['velocity'][i] = (w * particles['velocity'][i] +c1 * r1 * (particles['best_position'][i] - particles['position'][i]) +c2 * r2 * (global_best_pos - particles['position'][i]))# 更新位置particles['position'][i] += particles['velocity'][i]# 邊界處理for d in range(dim):low, high = bounds[d]particles['position'][i][d] = np.clip(particles['position'][i][d], low, high)# 計算新適應度new_fitness = func(particles['position'][i])# 更新個體最優if new_fitness < particles['best_fitness'][i]:particles['best_position'][i] = particles['position'][i].copy()particles['best_fitness'][i] = new_fitness# 更新全局最優if new_fitness < global_best_fit:global_best_fit = new_fitnessglobal_best_pos = particles['position'][i].copy()# 記錄本次迭代狀態particles['history'].append({'positions': particles['position'].copy(),'global_best': global_best_pos.copy()})return global_best_pos, global_best_fit, particles['history']def visualize_optimization(history, bounds, func):"""可視化優化過程"""fig = plt.figure()  # 增加圖形寬度以容納colorbar# 控制變量animation_state = {'paused': False,'stopped': False,'current_frame': 0}# 創建函數網格x = np.linspace(bounds[0][0], bounds[0][1], 100)y = np.linspace(bounds[1][0], bounds[1][1], 100)X, Y = np.meshgrid(x, y)Z = np.zeros_like(X)for i in range(X.shape[0]):for j in range(X.shape[1]):Z[i, j] = func([X[i, j], Y[i, j]])# 創建子圖ax1 = fig.add_subplot(121, projection='3d')ax2 = fig.add_subplot(122)# 初始化繪圖元素surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.3)particles_3d = ax1.scatter([], [], [], c='red', s=30)best_3d = ax1.scatter([], [], [], c='gold', s=100, marker='*')contour = ax2.contourf(X, Y, Z, 20, cmap='viridis', alpha=0.6)# 添加等高線標簽contour_lines = ax2.contour(X, Y, Z, 10, colors='black', alpha=0.4, linewidths=0.5)ax2.clabel(contour_lines, inline=True, fontsize=8, fmt='%.1f')# 添加顏色條cbar = plt.colorbar(contour, ax=ax2, shrink=0.8)cbar.set_label('適應度值', rotation=270, labelpad=15)particles_2d = ax2.scatter([], [], c='red', s=30, alpha=0.7)best_2d = ax2.scatter([], [], c='gold', s=100, marker='*')# 設置標簽和標題ax1.set_xlabel('X')ax1.set_ylabel('Y')ax1.set_zlabel('f(X,Y)')ax1.view_init(elev=30, azim=45)ax2.set_xlabel('X')ax2.set_ylabel('Y')ax2.set_xlim(bounds[0])ax2.set_ylim(bounds[1])# 創建動畫函數def on_key_press(event):"""鍵盤事件處理"""if event.key == ' ':  # 空格鍵暫停/繼續animation_state['paused'] = not animation_state['paused']status = "暫停" if animation_state['paused'] else "繼續"print(f"動畫{status}")def show_final_results():"""顯示最終結果并暫停"""final_frame = len(history) - 1final_positions = history[final_frame]['positions']final_best = history[final_frame]['global_best']# 更新到最終狀態particles_3d._offsets3d = (final_positions[:, 0], final_positions[:, 1], [func(p) for p in final_positions])best_3d._offsets3d = ([final_best[0]], [final_best[1]], [func(final_best)])particles_2d.set_offsets(final_positions)best_2d.set_offsets([final_best])# 更新標題顯示最終結果final_fitness = func(final_best)ax1.set_title(f'優化完成! 總迭代: {final_frame}\n最優值: {final_fitness:.6f}\n最優位置: [{final_best[0]:.3f}, {final_best[1]:.3f}]')ax2.set_title(f'最終結果 - 誤差: {abs(final_fitness):.6f}')plt.draw()print(f"\n=== 優化結果 ===")print(f"最優位置: [{final_best[0]:.6f}, {final_best[1]:.6f}]")print(f"最優適應度值: {final_fitness:.6f}")print(f"理論最優值: 0.0000 (在位置 [0, 0])")print(f"誤差: {abs(final_fitness):.6f}")print("\n按任意鍵關閉窗口...")# 等待用戶輸入input()def update(frame):# 檢查是否被停止if animation_state['stopped']:return particles_3d, best_3d, particles_2d, best_2d# 檢查是否暫停if animation_state['paused']:return particles_3d, best_3d, particles_2d, best_2danimation_state['current_frame'] = frame# 當前迭代數據positions = history[frame]['positions']global_best = history[frame]['global_best']# 更新3D散點圖particles_3d._offsets3d = (positions[:, 0], positions[:, 1], [func(p) for p in positions])best_3d._offsets3d = ([global_best[0]], [global_best[1]], [func(global_best)])# 更新2D散點圖particles_2d.set_offsets(positions)best_2d.set_offsets([global_best])# 更新標題current_fitness = func(global_best)ax1.set_title(f'迭代 {frame}/{len(history)-1}\n當前最優值: {current_fitness:.4f}\n[空格]暫停')ax2.set_title(f'等高線圖 - 迭代 {frame}')# 如果是最后一幀，自動顯示最終結果if frame == len(history) - 1:ani.event_source.stop()# 延遲顯示最終結果plt.pause(0.5)  # 短暫暫停show_final_results()return particles_3d, best_3d, particles_2d, best_2d# 連接鍵盤事件fig.canvas.mpl_connect('key_press_event', on_key_press)# 創建動畫ani = FuncAnimation(fig, update, frames=len(history), interval=200, blit=False, repeat=False)  # 改為不重復plt.tight_layout(pad=2.0)  # 增加邊距以確保文字完整顯示# 顯示操作提示print("\n=== 操作說明 ===")print("[空格鍵] - 暫停/繼續動畫")print("[關閉窗口] - 退出程序")print("==================\n")return ani, figdef main():# 設置函數搜索邊界bounds = [(-5.12, 5.12), (-5.12, 5.12)]print("正在運行粒子群優化算法...")# 運行PSO優化best_pos, best_fit, history = pso_optimization(fitness_function, bounds)print(f"\n算法執行完成!")print(f"最優解位置: [{best_pos[0]:.6f}, {best_pos[1]:.6f}]")print(f"最優值: {best_fit:.6f}")# 可視化優化過程ani, fig = visualize_optimization(history, bounds, fitness_function)# 顯示實時動畫print("正在啟動粒子群算法可視化...")try:plt.show()except KeyboardInterrupt:print("\n程序被用戶中斷")print("程序結束")if __name__ == "__main__":main()

參考

【1】 粒子群優化算法(Particle Swarm Optimization, PSO)的詳細解讀
【2】一篇文章搞懂PSO（粒子群算法）理論講解＋Python代碼示例講解