数学建模社区-数学中国

标题: 一群"鸟"帮你解数学题：聊聊粒子群算法 [打印本页]
作者: 2744557306    时间: 2026-8-23 17:38
标题: 一群"鸟"帮你解数学题：聊聊粒子群算法

算法笔记 · 粒子群优化 · 附完整可运行代码

先讲一个真实的故事。

1995 年，两位研究者在模拟鸟群觅食：一群鸟在一片区域里随机飞，谁也不知道食物在哪。但每只鸟都记得自己飞过的最好位置，也能看到整个鸟群目前找到的最好位置。于是每只鸟的飞行策略变得很简单——一部分凭自己的经验继续飞，一部分朝群里公认的好方向靠拢。

结果让人惊讶：不需要任何"领导"，不需要知道食物在哪，鸟群竟能自发地聚拢到食物周围。

这两位研究者一位是电子工程师 James Kennedy，一位是社会心理学家 Russell Eberhart。他们把这个模拟简化成几个公式，就成了今天要讲的主角——粒子群优化算法（Particle Swarm Optimization，PSO）。顺带一提，这个算法诞生 25 年后，他们拿到了 IEEE 的诺贝尔级别的荣誉——它至今仍是工程优化里最常用的算法之一。

如果你正在准备数学建模竞赛，这篇文章能帮你在"优化类"赛题里多一件趁手的兵器。

一、把优化问题想成"找食物"

先忘掉算法，想一个具体场景：

你要找一个函数的最小值，比如 $$f(x, y) = x² + y²$$，最小值显然在原点。但现实问题里，函数可能长得像一片丘陵——到处是小山包（局部最优），真正的谷底（全局最优）藏在某个角落。梯度下降这类方法从一点出发，很容易卡在最近的小山包里下不来。

PSO 的思路完全不同：与其派一个人去找，不如派一群鸟去找。

[td]
概念
在算法里是什么
一只鸟 / 粒子（Particle）一个候选解（一组变量值）
位置（Position）解本身，比如 (x, y) = (3.2, -1.5)
速度（Velocity）下一步往哪个方向飞、飞多远
食物的位置函数的最小值点
离食物多近适应度：把位置代入函数算出的值
自己的经验个体极值 pbest：我飞过的最好位置
群体的情报全局极值 gbest：全群找到的最好位置
二、灵魂公式：只有两行

PSO 全部精髓浓缩在两条更新规则里。每一代（每一轮迭代），每只鸟这样飞：


; y6 x; U( F1 o3 p$$v_{k+1} = w \cdot v_k + c_1 r_1 \cdot (pbest - x_k) + c_2 r_2 \cdot (gbest - x_k)$$' e( \, |$ S+ q/ W: G. U0 D: \

( `7 e7 h2 Z2 P: [: t& x7 E$$x_{k+1} = x_k + v_{k+1}$$) {3 N. q; P  C; ]$ x0 [* P

第一条是速度更新，即"下一刻往哪飞"。拆开看，一共三股力量在拉扯这只鸟：

[td]
部分
含义
生活翻译
w·vk惯性：保留原来的飞行方向"我正往这边飞呢，别让我急转弯"
c₁r₁·(pbest − x)认知项：朝自己最好位置飞"我以前在那儿吃过食，回去看看"
c₂r₂·(gbest − x)社会项：朝全群最好位置飞"老王说他找到吃的了，跟过去！"

第二条是位置更新，简单到不像话：新位置 = 旧位置 + 速度。就是"按刚才决定的方向飞过去"。

这里的 r₁、r₂ 是 0 到 1 之间的随机数——这点很关键。如果每次都精确地朝 pbest 和 gbest 飞，所有鸟很快会挤成一条线，搜索就死了。加上随机扰动，鸟群才能保持"散着搜索、又大致向好看的方向去"的状态。

一句话记住 PSO： 每只鸟被三股力拉着——惯性（接着飞）、自己的记忆（回自己找到过的好地方）、群体的情报（去大家公认的好地方）。飞到新位置后，更新记忆，再飞。如此反复，鸟群整体就向最优解聚拢。

三、三个参数怎么调[td]
参数
名字
常用值
调大调小会怎样
w惯性权重0.6 ~ 0.8大→搜得广但收敛慢；小→收敛快但易陷局部最优
c₁认知因子1.5 ~ 2大→各飞各的，多样性好
c₂社会因子1.5 ~ 2大→一窝蜂跟风，收敛快但容易早熟

工程上还有个进阶技巧：惯性权重线性递减——从 0.9 慢慢降到 0.4。前期惯性大，鸟群撒开网全局搜索；后期惯性小，大家稳稳地精修局部。这被称为"先广搜、后精修"，是 PSO 最常用的改良。

四、上代码：70 行实现完整 PSO

下面这段代码零依赖，保存成 pso_demo.py 直接运行。我们拿两个经典函数开刀：Sphere（单峰，考收敛速度）和 Rastrigin（满地局部最优，考跳出陷阱的能力）。

  1. # -*- coding: utf-8 -*-
    . X, T& ^# M: V
  2. """粒子群优化算法（PSO）完整演示，零依赖可直接运行"""; P* T, v/ S6 p2 b
  3. import random
    ) o" e& U( f: o
  4. import math
    + {% g# ~* g7 S& r
  5. ! g: U# J6 T) b8 F
  6. def sphere(x):  D8 u3 |- h& C0 m& J: s
  7.     """Sphere 函数：f(0,0,...,0)=0，单峰""": b. i- S  S1 W$ Z( l" k
  8.     return sum(xi ** 2 for xi in x)/ p0 O2 [2 w- _0 R' b/ _  z* n0 B

  9. 0 b) ?4 T0 p# w3 @; V
  10. def rastrigin(x):+ O% D  A% ^& s* |% f( K  g
  11.     """Rastrigin 函数：大量局部最优点，全局最优仍在原点"""6 U. |9 Z9 m2 g, w! x
  12.     return sum(xi**2 - 10*math.cos(2*math.pi*xi) + 10 for xi in x)6 K1 h2 [  Z. D8 `

  13. ! B/ s  D6 ?: d* O
  14. class Particle:
    4 q5 s% @$ X1 O. f
  15.     """一只『鸟』：有位置、速度，记得自己飞过的最好地方"""
    $ J9 \5 P1 t; u3 G' @
  16.     def __init__(self, dim, bounds):' k; h: N" p; {/ Z  ?' p2 N; k* A
  17.         self.pos = [random.uniform(bounds[i][0], bounds[i][1]) for i in range(dim)]$ }8 t, N) s: Y; R
  18.         self.vel = [random.uniform(-1, 1) for _ in range(dim)]
    % l& Y' Y- |8 q: W- G
  19.         self.best_pos = list(self.pos)      # 个体极值 pbest: e" p% I, ?' M4 _/ |' C
  20.         self.best_val = float("inf")
    4 ~* @' u) ]& l* ?

  21. ) [) `4 g% B, F/ S9 U
  22.     def update_velocity(self, gbest, w, c1, c2):2 `' A9 o9 ~8 Z/ V% _$ H# s8 k" Q
  23.         # 速度更新：v = w*v + c1*r1*(pbest-x) + c2*r2*(gbest-x)- L# n" _5 E1 i/ @: g5 l; V
  24.         for i in range(len(self.vel)):" E' `. g. r, r. i. ^- P; O2 X1 ~( p
  25.             r1, r2 = random.random(), random.random()
    2 X0 K# Z7 {( t
  26.             cognitive = c1 * r1 * (self.best_pos[i] - self.pos[i])  # 自己的经验
    ; E% a1 u* s. Q
  27.             social    = c2 * r2 * (gbest[i] - self.pos[i])          # 群体的情报1 ?) {; A7 ?! E  j
  28.             self.vel[i] = w * self.vel[i] + cognitive + social
    1 K( w$ ^; k1 {1 k
  29. $ h& ~: O4 o' C1 J) p- Y
  30.     def update_position(self, bounds):4 C! @  U! c' A6 S& a
  31.         # 位置更新：x = x + v，并夹回边界5 E9 D* y0 f' S* K2 n2 a6 a% O: ^
  32.         for i in range(len(self.pos)):
    ! c8 R* X, j$ L2 W( c# f, \& \
  33.             self.pos[i] += self.vel[i]
    * ^- d- l- ]- M8 k' g5 ?. j1 Y" n2 r
  34.             lo, hi = bounds[i]
    ' a# ], i/ a% j# o5 R
  35.             if self.pos[i] < lo:   self.pos[i] = lo' p% H* g, c# D
  36.             elif self.pos[i] > hi: self.pos[i] = hi
    " ?8 P$ n: g8 o5 h
  37. . u4 G5 ^- o* r- z
  38. def pso(func, dim=2, bounds=None, n_particles=30, max_iter=100,
    , V, g1 }/ n" b" \- b, ]( m( ]
  39.         w=0.7, c1=1.5, c2=1.5, verbose_every=10):# w! |& f5 {$ W& \8 X! {1 s
  40.     """标准 PSO 主流程"""
      x3 b- T' F* a/ h$ L
  41.     if bounds is None:
    4 _! f! Z2 u- x/ f  {4 G% i
  42.         bounds = [(-5.12, 5.12)] * dim. P4 t) ?5 S; a+ a3 Z
  43.     particles = [Particle(dim, bounds) for _ in range(n_particles)]
    ' |+ X: m0 u* [0 Z
  44.     gbest, gbest_val = list(particles[0].pos), float("inf")
    & W! l) m& Z5 }+ N) C* i
  45.     history = []
    : z3 a) _5 n* X, ~4 j

  46. 5 m! O. n1 e0 K$ }( b# C3 v
  47.     for it in range(max_iter):
    % s. q1 N1 F$ W% A2 I
  48.         # 第 1 步：评估每只鸟，更新 pbest 和 gbest
    # L3 J+ A, j6 M) [9 G
  49.         for p in particles:
    ( `, |( I" G9 N$ `
  50.             val = func(p.pos)
    & ]1 T) j$ ^2 E* }
  51.             if val < p.best_val:9 ?7 W0 W4 W; A3 y
  52.                 p.best_val, p.best_pos = val, list(p.pos)
    & L5 l$ c- J9 F- V1 X# {9 O
  53.             if val < gbest_val:
    7 f4 r+ W7 z; Y8 g0 ~
  54.                 gbest_val, gbest = val, list(p.pos)
    9 g7 |: v& _9 O+ ~4 x9 C2 m
  55.         history.append(gbest_val)
    / f0 r! O( g6 Z( \# C$ p' {

  56. ! i8 M6 B; m/ G! ~& \
  57.         # 第 2 步：更新速度和位置
    # a1 A7 G) [1 y; k8 L
  58.         for p in particles:0 q% w+ H/ i3 j
  59.             p.update_velocity(gbest, w, c1, c2)
    : i+ B8 U* J: K
  60.             p.update_position(bounds)  O% ]; Z: P! f3 ^
  61. 7 x* w  o$ j+ u$ [
  62.         if verbose_every and (it + 1) % verbose_every == 0:
    - F7 c0 W! D) g8 X% x$ C: @
  63.             print(f"    第 {it+1:3d} 代 | 全局最优值 = {gbest_val:.6f}")8 n* {) W) d9 ~. [/ w) ?+ a$ c: L
  64.     return gbest, gbest_val, history
    2 r: S4 O( k% ?' E* C9 k
  65. 6 c: C/ o$ q0 B- V
  66. if __name__ == "__main__":# c, F8 P0 c% H
  67.     # 实验 1：单峰函数，看收敛速度
    1 K$ ?0 B3 V5 C0 w1 ?/ M
  68.     print("[实验 1] Sphere 函数（理论最优值 = 0）")& T. G' W: I4 J4 Q; C) ]( ?$ J5 o
  69.     best_pos, best_val, _ = pso(sphere, dim=2, bounds=[(-10, 10)]*2,8 _' J) w2 _0 \$ P2 R% O
  70.                                 n_particles=30, max_iter=100)& ]% i7 e9 M4 H- m" [5 S( p4 Q
  71.     print(f"  >>> 最优位置: ({best_pos[0]:.4f}, {best_pos[1]:.4f})"): f3 F) l/ s9 q6 O8 N7 t
  72.     print(f"  >>> 最优值: {best_val:.8f}")$ s1 b$ u% E1 J/ X' \, l
  73. $ @5 V# ?$ Z7 }; E7 R0 _, u
  74.     # 实验 2：多峰函数，看跳出局部最优& t# E/ f4 _% q& v% M1 W" I
  75.     print("[实验 2] Rastrigin 函数（局部最优遍地都是）"); j2 L6 P' C. F) F; i
  76.     best_pos, best_val, _ = pso(rastrigin, dim=2,$ A! F4 p- E. w% k1 M- `! S
  77.                                 n_particles=50, max_iter=200)
    6 }* y" N. f; J0 I- _7 A
  78.     print(f"  >>> 最优值: {best_val:.6f}")
    7 m9 ?4 B$ h% w2 c* v3 W) a8 k5 P; u* t

  79. 0 U4 j% T1 K' B0 i5 S% h
  80.     # 实验 3：惯性权重对比（10 维）: c  Y/ _; |: ]& h+ E% b( ]/ g
  81.     print("[实验 3] 惯性权重 w 的影响（10 维 Rastrigin）")
    4 q2 Q/ Q7 q% h" c' A3 Q: F# q
  82.     for w in (0.2, 0.7, 1.2):7 |) h, j9 f' _, K/ H( j
  83.         _, val, hist = pso(rastrigin, dim=10, bounds=[(-5.12, 5.12)]*10,$ g$ ], E3 M3 L, l  }; H* R5 O
  84.                            n_particles=50, max_iter=200, w=w, verbose_every=0)
    : t% U! @* ?0 Y: I& J" Y1 e8 u8 }
  85.         print(f"  w = {w:<4} | 最终最优值 = {val:10.5f}")
复制代码

三个实验分别说明：

实验 1：单峰函数太简单，50 只鸟迭代 50 代就精确命中原点，误差小于百万分之一。

实验 2：Rastrigin 函数像个布满小坑的月球表面，梯度下降几乎必卡坑里，PSO 靠"一群鸟分散搜索"硬是跳了出来。

实验 3最有意思——升到 10 维后参数的差别立刻显现：w=1.2 时惯性太大，鸟群一直"飘着"落不下来（101.6，基本没搜到）；w=0.2 时收敛太快，早熟地卡在局部最优（24.9）；w=0.7 恰到好处（7.0）。这就是参数敏感性的现场教学。

五、什么时候用 PSO，什么时候别用适合的场景[td]
场景
举例
函数不可导 / 不光滑目标函数是查表、仿真得出的，没有解析式
多峰、怕陷入局部最优神经网络调参、TSP 路径优化、portfolio 投资组合
混合整数 / 离散问题稍加改造即可处理 0-1 变量（粒子位置取整）
需要"能跑就行"的快速方案数模竞赛限时三天，PSO 半小时就能调通
要小心的地方

1. 不是全局最优保证。 PSO 是启发式算法，理论上可能错过全局最优。论文里别写"证明了最优解"，要写"搜索到的近似最优解"。

2. 高维会退化。 维度上百时收敛变慢，可考虑增加粒子数或换差分进化（DE）。

3. 结果有随机性。 每次运行结果略有不同。写论文时要多次运行取平均/最优，并报告标准差——评委很看重这一点。

4. 约束处理要自己写。 常见做法是罚函数：越界就往目标函数上加一个很大的惩罚。

和其他算法的关系[td]
算法
思想
和 PSO 的区别
梯度下降沿坡往下走单点出发、需可导；PSO 群体出发、不需要导数
遗传算法 GA选择、交叉、变异有"淘汰"机制；PSO 不淘汰，所有粒子活到最后、共享信息
模拟退火 SA以一定概率接受坏解单点搜索；PSO 是群体协作
差分进化 DE向量差扰动与 PSO 各有胜负，高维时 DE 常更稳
六、写给数模竞赛的你

回到竞赛场景。PSO 在数模里的出场率非常高，典型的用法：

题型一：优化决策。"如何分配 X 使总成本最小/总收益最大"，目标函数复杂、变量多、不可导——PSO 直接上。

题型二：模型调参。 先用回归/机器学习建了模型，再有一堆超参数要调。可以用 PSO 包在外面做参数寻优，论文里就是个漂亮的"双层优化"结构。

题型三：拟合非线性模型。 最小二乘拟合非线性曲线时，初值选不好就拟合歪。用 PSO 先粗搜一轮找初值，再交给 LM 算法精修——这是很扎实的组合拳写法。

论文里怎么写才加分： 不要只丢一个结果。写清楚三件事——① 粒子数、迭代数、参数 w/c₁/c₂ 怎么定的（最好说明试过几组）；② 收敛曲线（画"迭代代数 vs 最优值"的图，评委一眼看到收敛行为）；③ 多次运行的结果稳定性（跑 10 次报告均值和标准差）。做到这三点，PSO 就从"套模型"升级成了"用模型"。

写在最后

PSO 之所以流传三十年，恰恰因为它简单得"不像话"——两条公式，没有 leaders，没有淘汰，只有一群朴素的个体，各自带着一点点记忆和一点点从众，互相分享情报，就把一片茫茫的搜索空间梳了个遍。

有人说这是对"群体智慧"最好的计算演绎：没有一只鸟知道答案，但鸟群知道。

下次遇到解不出的优化题，不妨也撒一把粒子出去——记得先广搜，后精修。

如果这篇文章帮到了你，点个「在看」让更多备赛的朋友看到


% s" k- `; k% @  k& a0 U6 g& `; ]* V7 V  r# s" H7 v
7 l5 |/ z$ R0 x8 b" |

9 G2 u2 X* c' T
$ Q8 |; \' h2 p  y) a, ]
6 n- h/ A" V( c! V$ m' y( {+ u9 ]! h



欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5