3. 动手实现第一个PSO

前两节我们建立了PSO的直觉：每个粒子带着位置和速度在搜索空间中飞行，速度由“个体历史最优”和“群体历史最优”共同牵引。这一节要把这套直觉变成可运行的代码。我们只用 NumPy，不依赖任何优化库，从零写出一个基础 PSO，并在 Sphere 函数上跑起来，重点观察两件事：粒子轨迹如何被 pbest 与 gbest 拉扯，以及全局最优值如何随迭代下降。

一、先明确要实现的数学模型

设搜索空间维度为 D，粒子数为 N。第 i 个粒子在第 t 代的位置记为 x_i，速度记为 v_i，个体历史最优位置为 p_i，群体历史最优位置为 g。基础 PSO 的速度与位置更新为：

v_i(t+1) = w * v_i(t)
         + c1 * r1 * (p_i - x_i(t))
         + c2 * r2 * (g  - x_i(t))
x_i(t+1) = x_i(t) + v_i(t+1)

其中 w 是惯性权重，c1、c2 是学习因子，r1、r2 是每次更新时在 [0,1] 上独立采样的随机数。注意 r1、r2 是逐粒子、逐维度重新采样的，这一点在写代码时最容易写错——如果整代只采样一次，多样性会明显下降。本节示例参数取：N=30，D=2（便于画轨迹），w=0.7，c1=c2=1.5，最大速度 vmax=1.0，迭代 60 代。这些是教学用的合理取值，不是调优结果。

二、Sphere 函数与边界处理

Sphere 函数 f(x)=Σx_d²，全局最优在原点，最优值 0。它是个碗形单峰函数，适合第一次验证算法是否正确：如果连 Sphere 都收不拢，问题多半在更新公式或边界处理上。我们把搜索范围设为 [-5, 5]^D，粒子飞出边界时采用“截断 + 速度反向减半”的简单策略，既保证位置合法，又避免粒子贴边后彻底失去运动能力。

三、完整实现代码

import numpy as np

def sphere(x):
    return np.sum(x ** 2, axis=1)

def pso_sphere(n_particles=30, dim=2, iters=60,
               w=0.7, c1=1.5, c2=1.5,
               vmax=1.0, bound=5.0, seed=42):
    rng = np.random.default_rng(seed)
    # 初始化位置与速度
    x = rng.uniform(-bound, bound, (n_particles, dim))
    v = rng.uniform(-vmax, vmax, (n_particles, dim))

    pbest = x.copy()
    pbest_val = sphere(x)
    g_idx = np.argmin(pbest_val)
    gbest = pbest[g_idx].copy()
    gbest_val = pbest_val[g_idx]

    history = [gbest_val]
    traj = [x.copy()]          # 记录每代所有粒子位置

    for t in range(iters):
        r1 = rng.random((n_particles, dim))
        r2 = rng.random((n_particles, dim))
        v = (w * v
             + c1 * r1 * (pbest - x)
             + c2 * r2 * (gbest - x))
        # 速度限幅
        v = np.clip(v, -vmax, vmax)
        x = x + v
        # 边界处理：截断并反向减半
        out = np.abs(x) > bound
        x = np.clip(x, -bound, bound)
        v[out] *= -0.5

        val = sphere(x)
        better = val < pbest_val
        pbest[better] = x[better]
        pbest_val[better] = val[better]

        if pbest_val.min() < gbest_val:
            g_idx = np.argmin(pbest_val)
            gbest = pbest[g_idx].copy()
            gbest_val = pbest_val[g_idx]

        history.append(gbest_val)
        traj.append(x.copy())

    return gbest, gbest_val, np.array(history), np.array(traj)

gbest, gbest_val, hist, traj = pso_sphere()
print("最优位置:", gbest)
print("最优值:", gbest_val)
print("前5代最优值:", hist[:5])
print("末5代最优值:", hist[-5:])

代码里有几个关键点值得对照检查。第一，pbest 与 pbest_val 必须同步更新，只更新值不更新位置会让牵引方向错乱。第二，gbest 用 .copy() 保存，否则后续 pbest 被原地修改时会连带改掉 gbest。第三，速度限幅放在位置更新之前，边界处理放在位置更新之后，顺序不能颠倒。

四、观察收敛过程与粒子轨迹

运行后，hist 是一条长度 61 的曲线。在示例参数下，典型现象是：前 10 代下降很快，因为粒子散布在 [-5,5] 内，梯度信息强；20 代之后曲线变平，进入精细搜索；到 60 代时最优值通常已降到 10^-3 量级以下，但未必精确为 0。如果曲线在中途出现平台后又突然下降，往往说明有粒子被边界截断后重新获得速度，属于正常现象。

轨迹比曲线更直观。traj 的形状是 (61, 30, 2)，取任意一个粒子，把它的 61 个位置连成折线，就能看到它先朝自己的 pbest 方向摆动，再被 gbest 拉向原点。下图把这一过程抽象成“速度合成的三个分量”与“收敛曲线形态”的对照，帮助你在看实际轨迹时知道该关注什么。

看轨迹时建议关注三个特征。其一，粒子是否在早期出现“集体转向”——所有粒子几乎同时朝同一方向移动，说明 gbest 已经形成强吸引。其二，是否有粒子长时间在边界附近徘徊，这通常意味着 vmax 偏大或边界策略过于激进。其三，末代粒子是否聚集在原点附近但仍有小幅抖动，这是 w 与 c1、c2 共同作用下的正常精细搜索，不必强行压到零。

五、用几个小实验验证实现是否正确

第一次跑通不代表实现正确，建议做三组对照。第一组把 c2 设为 0，只保留个体认知，此时粒子应各自收敛到自己的初始附近，gbest 下降缓慢且不集中，这验证了社会项确实在起作用。第二组把 w 设为 0，粒子失去惯性，轨迹会变得“急转弯”，收敛更快但更容易早熟，这验证了惯性项的存在。第三组把 vmax 调到 0.05，粒子几乎挪不动，60 代后最优值仍较大，这验证了速度限幅确实生效。

到这里，你已经拥有一个可运行、可验证、可修改的基础 PSO。后续调参时，所有对 w、c1、c2、vmax 的讨论都可以直接落到这段代码上做对照实验，而不是停留在概念层面。