【群体智能优化算法系列 】一 粒子算法
- 一:前言
- 二:算法原理
- 2.1 核心思想
- 2.2 PSO核心公式
- 2.3 PSO算法流程图
- 三:python实现 二维Rastrigin函数 最低点检索例子
- 参考
一:前言
粒子群算法是由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代码示例讲解