粒子群算法（PSO）原理与Python实现：从方程求根到优化问题求解-编程实验室

1. 从“暴力搜索”到“智能寻优”：为什么我们需要粒子群算法

在工程优化、参数调优甚至日常的数学建模中，我们常常会遇到一个核心问题：如何找到一个函数的最优解？对于简单的一元或二元方程求根，我们有很多经典方法，比如牛顿迭代法、二分法。但现实世界的问题往往复杂得多。比如，一个成本函数可能有几十个变量，而且这个函数可能“坑坑洼洼”（非凸），导数难以计算，甚至根本不存在解析形式。这时候，传统的基于梯度的方法就“抓瞎”了。

这就好比在一片漆黑、地形复杂的山区里寻找海拔最低点。牛顿法这类方法像是一个拿着精密地图和指南针的探险家，要求地形必须光滑可导（有地图）。但如果地图是错的，或者根本没有路（不可导），他就可能掉进悬崖（局部最优）或者彻底迷路（发散）。而粒子群算法（PSO）则像是一群被放出去的无人机，它们之间通过简单的通信规则，分享各自发现的好位置，最终引导整个群体向真正的“盆地”聚集。它不依赖梯度，对函数性质要求极低，只要你能计算出每个点的函数值就行，这种“黑盒优化”的特性让它大放异彩。

我最初接触PSO是为了调一个机器学习模型的超参数，手动网格搜索太慢，随机搜索又像碰运气。PSO用下来，感觉它像是一个有“集体智慧”的搜索策略，代码实现简单，效果却经常出人意料的好。今天，我们就从最基础的部分入手，用Python实现基础PSO，并把它应用到一个看似简单、实则能很好理解其原理的问题上：求解一元和二元方程的根。你会发现，把求根问题转化为优化问题，再用PSO来解决，是一个非常直观且富有启发性的实践。

2. 粒子群算法的核心思想：鸟群觅食的数学抽象

粒子群优化算法的灵感来源于鸟群或鱼群的社会行为。想象一下，一群鸟在随机搜索一片区域的食物。每只鸟都不知道食物具体在哪，但它们会记住自己飞过的最好位置（个体最优），同时也会通过鸣叫等方式，得知整个鸟群目前发现的最好位置（全局最优）。每只鸟决定下一步怎么飞，就基于三个因素：自己当前的飞行惯性、飞向自己曾找到的最好位置的意愿，以及飞向群体发现的最好位置的意愿。

我们将这个生动的场景转化为数学模型，需要定义几个核心概念和公式。理解这些公式，是写出正确代码和进行有效调参的基础。

2.1 算法中的“粒子”与它的记忆

在PSO中，每一只“鸟”被称为一个“粒子”（Particle）。对于一个D维的优化问题（比如求二元方程根，D=2），每个粒子i在时刻t具有以下属性：

位置（Position）：一个D维向量，记为 ( x_i(t) = [x_{i1}, x_{i2}, ..., x_{iD}] )。这代表粒子在搜索空间中的当前位置。在我们的问题里，位置就是方程变量的候选解。
速度（Velocity）：一个D维向量，记为 ( v_i(t) = [v_{i1}, v_{i2}, ..., v_{iD}] )。这决定了粒子下一步移动的方向和距离。
个体历史最优位置（Personal Best, pbest）：记为 ( p_i )。这是粒子i自搜索开始以来，它所经过的所有位置中，使得目标函数值最优（对于最小化问题，就是函数值最小）的那个位置。粒子有“记忆”，记住自己的高光时刻。
全局历史最优位置（Global Best, gbest）：记为 ( g )。这是整个粒子群中，所有粒子发现的个体最优位置里，最好的那一个。这是群体的“共识”和努力方向。

2.2 驱动粒子飞行的更新公式

粒子如何从时刻t的状态更新到t+1时刻？核心在于速度和位置的更新公式，这是PSO的“发动机”。

速度更新公式：[ v_{id}(t+1) = w \cdot v_{id}(t) + c_1 \cdot r_1 \cdot (p_{id} - x_{id}(t)) + c_2 \cdot r_2 \cdot (g_d - x_{id}(t)) ] 这个公式是PSO的灵魂，它融合了之前提到的三个驱动因素：

惯性部分（( w \cdot v_{id}(t) )）：w称为惯性权重。它代表了粒子维持先前速度的趋势。w较大时，粒子探索能力强，飞行速度快，适合全局搜索；w较小时，粒子开发能力强，飞行精细，适合局部精细搜索。通常w取值在0.4到0.9之间，也可以随着迭代线性递减，实现“先粗搜，后细找”的策略。
认知部分（( c_1 \cdot r_1 \cdot (p_{id} - x_{id}(t)) )）：c1称为个体学习因子。r1是一个[0,1]之间的随机数。这部分驱使粒子飞向它自己曾找到的最好位置，代表了粒子对自身经验的信赖。
社会部分（( c_2 \cdot r_2 \cdot (g_d - x_{id}(t)) )）：c2称为社会学习因子。r2是另一个[0,1]之间的随机数。这部分驱使粒子飞向群体发现的最好位置，代表了粒子向群体学习、追随共识的倾向。

位置更新公式：[ x_{id}(t+1) = x_{id}(t) + v_{id}(t+1) ] 这个公式很简单，就是用更新后的速度来移动粒子，得到新的位置。

2.3 算法流程与关键参数

理解了核心公式，整个算法的流程就清晰了。一个标准的PSO流程如下：

初始化：在搜索空间内，随机初始化一群粒子（比如20-50个）的位置和速度。同时，将每个粒子的当前位置设为它的个体最优pbest，并从中找出全局最优gbest。
迭代优化：对于每一轮迭代（比如100-500轮）： a.评估适应度：计算每个粒子新位置的适应度值（即目标函数值）。 b.更新个体最优：如果某个粒子当前位置的适应度优于其pbest的适应度，则用当前位置更新它的pbest。 c.更新全局最优：检查所有更新后的pbest，找出适应度最优的那个，更新gbest。 d.更新速度和位置：对每个粒子的每一维，应用上面的速度更新和位置更新公式。
终止与输出：达到最大迭代次数或满足其他停止条件（如gbest连续多轮无变化）后，算法停止，输出最终的gbest作为找到的最优解。

关键参数总结：

粒子数量（n_particles）：粒子越多，搜索能力越强，但计算开销也越大。通常20-50是个不错的起点。
维度（dim）：由你的优化问题决定。求一元方程根，dim=1；求二元方程根，dim=2。
惯性权重（w）：控制搜索范围。常用策略是从0.9线性递减到0.4。
学习因子（c1, c2）：平衡个体经验和群体智慧。经典设置是c1 = c2 = 2.0。
速度限制（v_max）：为防止粒子飞离搜索空间，通常会对速度进行钳制，例如v = np.clip(v, -v_max, v_max)。
位置边界（bounds）：定义搜索空间的范围，粒子位置会被限制在此边界内。

3. 将方程求根转化为优化问题：定义“适应度函数”

PSO是一个优化算法，它默认是寻找目标函数的最小值。那么，如何用它来求解方程 ( f(x) = 0 ) 或 ( f(x, y) = 0 ) 的根呢？

窍门在于构造一个合适的“适应度函数”（Fitness Function）或“目标函数”。我们的目标是找到使方程成立的x或(x, y)。一个最直接的想法是：方程成立时，等式左边与右边的差（即 ( f(x) ) 或 ( f(x, y) ) ）的绝对值应该为0。因此，我们可以定义适应度函数为这个差值的绝对值（或平方）。

对于一元方程 ( f(x) = 0 )：适应度函数定义为：fitness(x) = abs(f(x))我们的目标就变成了寻找一个x，使得fitness(x)的值最小，理想情况下为0。这个最小值点对应的x就是方程的根。

对于二元方程组 ( f(x, y) = 0 ) 和 ( g(x, y) = 0 )：我们需要同时满足两个方程。一个常见的处理方法是构造一个综合的适应度函数，例如两个方程绝对值的和：适应度函数定义为：fitness(x, y) = abs(f(x, y)) + abs(g(x, y))同样，寻找使这个和最小（趋近于0）的(x, y)，就是方程组的解。

注意：使用绝对值之和是一种简单有效的方法，但它隐含了一个假设：两个方程同等重要。如果两个方程的量级差异很大（比如一个结果是1e-6，另一个是1e3），直接相加可能会导致优化过程被量级大的方程主导。在实际复杂问题中，可能需要对各个分量进行归一化处理，或者使用平方和f^2 + g^2（可微，但PSO不关心是否可微）。

举个例子：假设我们要求解一元方程 ( x^2 - 4 = 0 )，它的根是x=2和x=-2。我们的适应度函数就是fitness(x) = abs(x**2 - 4)。 PSO会尝试寻找使abs(x**2 - 4)最小的x。当找到的x接近2或-2时，适应度值接近0，我们就认为找到了根。

4. Python实现基础PSO求解器：从零搭建代码框架

理论铺垫完成，现在进入实战环节。我们将一步步用Python实现一个基础但完整的PSO类，并用它来求解具体方程。我会在代码中穿插大量注释，解释每一步的意图和潜在陷阱。

4.1 构建PSO优化器类

首先，我们构建一个通用的PSO优化器。它不关心具体是什么方程，只关心我们传给它的适应度函数。

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class BasicPSO: """ 基础粒子群优化算法实现类。 用于求解最小化问题：minimize fitness_func(position) """ def __init__(self, fitness_func, dim, bounds, n_particles=30, w=0.8, c1=2.0, c2=2.0, max_iter=100, v_max_ratio=0.2): """ 初始化PSO优化器。 参数： fitness_func: 适应度函数，接受一个位置向量（形状为(dim,)或(n_particles, dim)），返回适应度值（标量或一维数组）。 dim: 问题的维度。 bounds: 搜索边界，列表形式，例如 [(x1_min, x1_max), (x2_min, x2_max), ...]。 n_particles: 粒子数量。 w: 惯性权重。 c1: 个体学习因子。 c2: 社会学习因子。 max_iter: 最大迭代次数。 v_max_ratio: 速度最大值与搜索范围的比例，用于计算v_max。 """ self.fitness_func = fitness_func self.dim = dim self.bounds = np.array(bounds) # 转换为数组便于计算 self.n_particles = n_particles self.w = w self.c1 = c1 self.c2 = c2 self.max_iter = max_iter # 计算搜索范围并确定速度限制 self.range_per_dim = self.bounds[:, 1] - self.bounds[:, 0] self.v_max = v_max_ratio * self.range_per_dim # 速度限制为搜索范围的20% # 初始化粒子群 self.positions = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1], (self.n_particles, self.dim)) self.velocities = np.random.uniform(-self.v_max, self.v_max, (self.n_particles, self.dim)) # 初始化个体最优和全局最优 self.personal_best_positions = self.positions.copy() # 评估初始适应度 self.personal_best_scores = self.fitness_func(self.positions) # 这里fitness_func应能处理批量输入 self.global_best_position = self.personal_best_positions[np.argmin(self.personal_best_scores)] self.global_best_score = np.min(self.personal_best_scores) # 记录历史用于分析 self.history_best_score = [] self.history_best_position = [] def optimize(self): """ 执行优化过程。 """ print(f"开始PSO优化，维度：{self.dim}, 粒子数：{self.n_particles}, 迭代次数：{self.max_iter}") print(f"初始全局最优解：{self.global_best_position}, 适应度：{self.global_best_score:.6e}") for iter in range(self.max_iter): # 1. 更新速度 r1 = np.random.rand(self.n_particles, self.dim) r2 = np.random.rand(self.n_particles, self.dim) cognitive = self.c1 * r1 * (self.personal_best_positions - self.positions) social = self.c2 * r2 * (self.global_best_position - self.positions) self.velocities = self.w * self.velocities + cognitive + social # 2. 应用速度限制（防止粒子飞得太快） self.velocities = np.clip(self.velocities, -self.v_max, self.v_max) # 3. 更新位置 self.positions += self.velocities # 4. 应用位置边界（将飞出边界的粒子拉回） # 方法1：直接钳制到边界（简单，但可能导致粒子聚集在边界） # self.positions = np.clip(self.positions, self.bounds[:, 0], self.bounds[:, 1]) # 方法2：反弹边界（更符合物理直觉，我更喜欢用这个） for d in range(self.dim): lower_bound, upper_bound = self.bounds[d] # 对于下界越界 below_mask = self.positions[:, d] < lower_bound self.positions[below_mask, d] = lower_bound self.velocities[below_mask, d] *= -0.5 # 反弹并损失部分能量 # 对于上界越界 above_mask = self.positions[:, d] > upper_bound self.positions[above_mask, d] = upper_bound self.velocities[above_mask, d] *= -0.5 # 5. 评估新位置的适应度 current_scores = self.fitness_func(self.positions) # 6. 更新个体最优 improved_mask = current_scores < self.personal_best_scores self.personal_best_positions[improved_mask] = self.positions[improved_mask] self.personal_best_scores[improved_mask] = current_scores[improved_mask] # 7. 更新全局最优 current_best_idx = np.argmin(self.personal_best_scores) current_best_score = self.personal_best_scores[current_best_idx] if current_best_score < self.global_best_score: self.global_best_score = current_best_score self.global_best_position = self.personal_best_positions[current_best_idx].copy() # print(f"迭代 {iter+1}: 发现新的全局最优，适应度：{self.global_best_score:.6e}") # 8. 记录历史 self.history_best_score.append(self.global_best_score) self.history_best_position.append(self.global_best_position.copy()) # 可选：提前终止条件（如果适应度已经足够好） if self.global_best_score < 1e-12: print(f"迭代 {iter+1}: 适应度已达阈值，提前终止。") break print(f"优化结束。最终全局最优解：{self.global_best_position}") print(f"最终适应度值：{self.global_best_score:.6e}") return self.global_best_position, self.global_best_score def plot_convergence(self): """ 绘制适应度收敛曲线。 """ plt.figure(figsize=(10, 6)) plt.plot(self.history_best_score, linewidth=2) plt.yscale('log') # 对数坐标能更清晰地看到适应度下降过程 plt.xlabel('迭代次数', fontsize=12) plt.ylabel('最佳适应度 (log scale)', fontsize=12) plt.title('PSO收敛曲线', fontsize=14) plt.grid(True, which='both', linestyle='--', alpha=0.7) plt.tight_layout() plt.show()

代码关键点解析与避坑经验：

适应度函数的批量处理：在__init__和optimize中，我们对所有粒子的位置self.positions（形状为(n_particles, dim)）一次性计算适应度。这就要求我们传入的fitness_func必须能够处理这种批量输入，并返回一个长度为n_particles的一维数组。这是为了利用NumPy的向量化运算，比在循环中逐个计算快成百上千倍。如果您的函数原本只处理单个粒子，可以用np.vectorize包装，或者稍作修改。
边界处理策略：代码中我提供了两种边界处理方式。直接钳制（np.clip）最简单，但缺点是粒子一旦撞到边界，速度分量会被置零，可能导致大量粒子“粘”在边界上，影响搜索效率。我采用的“反弹”策略（将越界位置设到边界，并让对应速度反向并衰减）更符合物理直觉，能让粒子在边界附近继续探索，通常效果更好。
速度限制（v_max）的重要性：如果没有速度限制，当w、c1、c2设置不当时，粒子的速度可能会指数级增长，导致粒子瞬间飞离搜索空间，算法立即失效。设置v_max是保证算法稳定的关键。我通常将其设置为搜索范围（bounds之差）的10%到30%。
惯性权重w的选择：代码中使用了固定的w。但在实践中，线性递减的惯性权重（LDW）策略往往效果更佳。你可以在optimize循环的开始加入：w = self.w_start - (self.w_start - self.w_end) * (iter / self.max_iter)。例如从0.9递减到0.4，这样前期全局探索能力强，后期局部开发能力强。

4.2 实战案例一：求解一元方程 ( x^3 - 2x - 5 = 0 )

这是一个经典的一元三次方程，在区间[2, 3]内有一个实根（大约2.0946）。我们用PSO来找到它。

# 定义一元方程的适应度函数 def fitness_unary(position): """ 适应度函数：f(x) = x^3 - 2x - 5 输入position可以是单个值，也可以是一维数组（多个粒子）。 返回适应度值（绝对值）。 """ # 确保输入是numpy数组，并处理批量计算 x = np.asarray(position) # 如果x是标量，np.abs会返回标量；如果x是数组，则返回数组。 return np.abs(x**3 - 2*x - 5) # 设置PSO参数并求解 dim = 1 # 一元方程，维度为1 bounds = [(-10, 10)] # 搜索范围，我们设定一个较宽的范围看看PSO能否找到 pso_unary = BasicPSO(fitness_func=fitness_unary, dim=dim, bounds=bounds, n_particles=20, w=0.7, # 使用固定惯性权重 c1=1.5, c2=1.5, max_iter=100, v_max_ratio=0.15) best_solution, best_fitness = pso_unary.optimize() print(f"\n方程 x^3 - 2x - 5 = 0 的近似根为：x = {best_solution[0]:.6f}") print(f"代入验证：f({best_solution[0]:.6f}) = {best_solution[0]**3 - 2*best_solution[0] - 5:.6e}") # 绘制收敛过程 pso_unary.plot_convergence() # 可视化搜索过程（可选，绘制函数曲线和粒子位置） x_plot = np.linspace(bounds[0][0], bounds[0][1], 400) y_plot = fitness_unary(x_plot) plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) plt.plot(x_plot, y_plot, 'b-', label='|f(x)|', linewidth=2) plt.scatter(pso_unary.positions, np.zeros_like(pso_unary.positions), c='red', alpha=0.6, label='最终粒子位置') plt.axvline(x=best_solution[0], color='green', linestyle='--', label=f'找到的根 ~{best_solution[0]:.4f}') plt.xlabel('x') plt.ylabel('|f(x)|') plt.title('一元方程 |x^3 - 2x - 5| 曲线与粒子分布') plt.legend() plt.grid(True) plt.subplot(1, 2, 2) plt.plot(pso_unary.history_best_position, marker='o', markersize=3, linestyle='-', alpha=0.7) plt.xlabel('迭代次数') plt.ylabel('全局最优位置 (x)') plt.title('全局最优位置随迭代的变化') plt.grid(True) plt.tight_layout() plt.show()

运行结果与分析：运行上述代码，你可能会得到类似x ≈ 2.094551的结果，适应度值在1e-7量级或更低。收敛曲线会显示适应度值在前几十次迭代中迅速下降，之后趋于平稳。

实操心得：对于一元方程，PSO看起来有点“杀鸡用牛刀”，因为二分法、牛顿法可能更快更准。但这个例子的意义在于验证我们的PSO框架工作正常，并且直观地展示了粒子如何在搜索空间中移动并聚集到最优点（根）附近。你可以尝试改变bounds范围，比如设为[(-100, 100)]，PSO依然有很大概率找到根，这体现了其全局搜索能力。而牛顿法如果初始点选得不好，很容易发散。

4.3 实战案例二：求解二元方程组

我们求解一个简单的非线性方程组： [ \begin{cases} x^2 + y^2 - 4 = 0 \ e^x + y - 1 = 0 \end{cases} ] 这个方程组在(x, y)平面内有解（可以通过数值方法验证，大约在(0.537, 0.629)和(-1.996, 6.389)附近）。我们构造适应度函数为两个方程绝对值的和。

# 定义二元方程组的适应度函数 def fitness_binary(position): """ 适应度函数：|f1(x,y)| + |f2(x,y)| 其中 f1 = x^2 + y^2 - 4 f2 = exp(x) + y - 1 输入position形状为 (n_particles, 2) 或 (2,) """ pos = np.asarray(position) # 处理批量输入和单个输入 if pos.ndim == 1: x, y = pos[0], pos[1] else: x, y = pos[:, 0], pos[:, 1] f1 = x**2 + y**2 - 4 f2 = np.exp(x) + y - 1 return np.abs(f1) + np.abs(f2) # 设置PSO参数并求解 dim = 2 bounds = [(-3, 2), (-3, 8)] # 根据对方程的粗略分析设定搜索范围 pso_binary = BasicPSO(fitness_func=fitness_binary, dim=dim, bounds=bounds, n_particles=40, # 二元问题，适当增加粒子数 w=0.8, c1=1.8, c2=1.8, max_iter=150, v_max_ratio=0.15) best_solution_2d, best_fitness_2d = pso_binary.optimize() print(f"\n方程组的近似解为：x = {best_solution_2d[0]:.6f}, y = {best_solution_2d[1]:.6f}") print(f"适应度（误差和）：{best_fitness_2d:.6e}") print(f"验证：") print(f" f1 = x^2+y^2-4 = {best_solution_2d[0]**2 + best_solution_2d[1]**2 - 4:.6e}") print(f" f2 = e^x+y-1 = {np.exp(best_solution_2d[0]) + best_solution_2d[1] - 1:.6e}") # 绘制收敛曲线 pso_binary.plot_convergence() # 可视化搜索过程（二维等高线图与粒子运动） 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 = fitness_binary(np.stack([X, Y], axis=-1).reshape(-1, 2)).reshape(100, 100) plt.figure(figsize=(15, 5)) plt.subplot(1, 3, 1) contour = plt.contour(X, Y, Z, levels=30, cmap='viridis') plt.colorbar(contour, label='适应度值 |f1|+|f2|') plt.scatter(pso_binary.positions[:, 0], pso_binary.positions[:, 1], c='red', s=20, alpha=0.7, label='最终粒子群') plt.scatter(best_solution_2d[0], best_solution_2d[1], c='gold', s=200, marker='*', edgecolors='black', label='全局最优解') plt.xlabel('x') plt.ylabel('y') plt.title('适应度函数等高线图与粒子分布') plt.legend() plt.grid(True) plt.subplot(1, 3, 2) plt.plot([pos[0] for pos in pso_binary.history_best_position], label='x') plt.plot([pos[1] for pos in pso_binary.history_best_position], label='y') plt.xlabel('迭代次数') plt.ylabel('参数值') plt.title('全局最优解 (x, y) 随迭代的变化') plt.legend() plt.grid(True) plt.subplot(1, 3, 3) # 绘制两个方程为零的曲线，交点即为解 x_line = np.linspace(bounds[0][0], bounds[0][1], 400) y_line1 = np.sqrt(4 - x_line**2) # x^2+y^2=4 -> y = sqrt(4-x^2) 和 -sqrt(...) y_line2 = 1 - np.exp(x_line) # e^x+y=1 -> y = 1 - e^x plt.plot(x_line, y_line1, 'b-', label='x^2 + y^2 = 4') plt.plot(x_line, -y_line1, 'b-') plt.plot(x_line, y_line2, 'r-', label='e^x + y = 1') plt.scatter(best_solution_2d[0], best_solution_2d[1], c='gold', s=200, marker='*', edgecolors='black', label='PSO解') plt.xlim(bounds[0]) plt.ylim(bounds[1]) plt.xlabel('x') plt.ylabel('y') plt.title('方程曲线与PSO解的位置') plt.legend() plt.grid(True) plt.tight_layout() plt.show()

运行结果与分析：PSO有很大概率会找到(x≈0.537, y≈0.629)这个解。从等高线图中可以清晰看到，粒子群最终聚集在了适应度函数（误差和）的“谷底”（颜色最蓝的区域）。第二个子图展示了最优解的x和y坐标是如何随着迭代逐步稳定下来的。第三个子图将PSO找到的解标在了两个方程曲线的交点附近，提供了非常直观的验证。

避坑经验：对于二元或更高维的问题，搜索范围（bounds）的设置非常关键。如果范围设得太大，粒子群需要探索的空间呈指数增长，找到精确解所需的迭代次数和粒子数会急剧增加，甚至可能失败。如果范围设得太小，可能直接错过了真解。一个实用的技巧是：先对方程进行粗略的分析或绘图，大致确定解可能存在的区域，再以此设定bounds。例如，对于x^2+y^2=4，我们知道解肯定在半径为2的圆内，这可以帮助我们设定x和y的大致范围。

5. 参数调优与算法改进：让PSO更高效、更稳定

基础PSO虽然能工作，但其性能严重依赖于参数设置。此外，它也有一些固有的缺陷，比如容易早熟收敛（陷入局部最优）。下面分享一些调参经验和简单的改进思路。

5.1 关键参数的影响与调参指南

粒子数量（n_particles）：
- 影响：粒子越多，搜索能力越强，探索空间更充分，但每次迭代的计算成本也越高。
- 建议：对于简单问题（低维、单峰），20-30个粒子足够。对于复杂问题（高维、多峰），可能需要50-200个甚至更多。这是一个需要权衡的参数。我的经验是，先从30-50开始，如果收敛效果不好，再适当增加。
惯性权重（w）：
- 影响：控制粒子的“探索-开发”平衡。高w（如0.9）利于全局探索，低w（如0.4）利于局部开发。
- 建议：强烈推荐使用线性递减权重（LDW）。例如w_start=0.9,w_end=0.4。这模拟了搜索过程从“粗放”到“精细”的自然过渡，效果通常比固定权重好很多。可以在optimize方法中动态计算：current_w = w_start - (w_start - w_end) * (iteration / max_iter)。
学习因子（c1, c2）：
- 影响：c1控制粒子向自身历史最优学习的强度（个体认知）；c2控制粒子向群体最优学习的强度（社会认知）。
- 建议：经典设置是c1 = c2 = 2.0。如果你想增强全局探索，可以适当增大c1（如2.5）并减小c2（如1.5），让粒子更依赖自己的经验。反之，如果想加快收敛，可以增大c2。也可以尝试让它们随着迭代变化。
速度限制（v_max）：
- 影响：防止粒子失速或爆炸。v_max太小，粒子移动慢，收敛慢；v_max太大，粒子可能跳过最优区域。
- 建议：通常设置为每个维度搜索范围（bounds[1]-bounds[0]）的10%~30%。这是一个比较稳健的范围。

5.2 基础改进策略：防止早熟收敛

基础PSO的一个常见问题是所有粒子快速聚集到当前全局最优gbest附近，如果gbest是一个局部最优解，算法就“僵住”了，很难再跳出来。以下是一些简单有效的改进策略：

引入随机扰动：以一定的小概率，随机重置部分粒子的位置或速度，为种群注入新的多样性。这被称为“变异”操作。

# 在optimize循环的合适位置（比如每20代）加入 if iter % 20 == 0 and iter > 0: # 随机选择10%的粒子进行重置 reset_idx = np.random.choice(self.n_particles, size=int(self.n_particles*0.1), replace=False) self.positions[reset_idx] = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1], (len(reset_idx), self.dim)) self.velocities[reset_idx] = np.random.uniform(-self.v_max, self.v_max, (len(reset_idx), self.dim)) # 记得也要更新这些粒子的pbest self.personal_best_positions[reset_idx] = self.positions[reset_idx].copy() self.personal_best_scores[reset_idx] = self.fitness_func(self.positions[reset_idx])

收缩因子（Constriction Factor）：这是一种更数学化的速度更新方式，可以保证算法收敛。速度更新公式变为： [ v_{id}(t+1) = \chi [ v_{id}(t) + \phi_1 r_1 (p_{id}-x_{id}(t)) + \phi_2 r_2 (g_d-x_{id}(t)) ] ] 其中 (\chi) 是收缩因子，通常根据 (\phi_1) 和 (\phi_2) 计算（(\phi = \phi_1 + \phi_2 > 4)）。当 (\phi_1=\phi_2=2.05) 时，(\chi \approx 0.729)。使用收缩因子时，通常不再需要惯性权重w和速度限制v_max，算法行为更稳定。
多种群PSO：将整个粒子群分成几个子群，每个子群有自己的局部最优lbest，子群之间定期交换信息。这有助于维持种群的多样性，避免全体粒子过早收敛到同一个点。