算法笔记 · 粒子群优化 · 附完整可运行代码
先讲一个真实的故事。 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 全部精髓浓缩在两条更新规则里。每一代(每一轮迭代),每只鸟这样飞: 8 r- r, ^6 ?, ]
$$v_{k+1} = w \cdot v_k + c_1 r_1 \cdot (pbest - x_k) + c_2 r_2 \cdot (gbest - x_k)$$- Q5 |+ m) C" r9 |
3 }# J |5 c$ r, M$$x_{k+1} = x_k + v_{k+1}$$" T% `% W- ?/ F8 }
第一条是速度更新,即"下一刻往哪飞"。拆开看,一共三股力量在拉扯这只鸟: [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(满地局部最优,考跳出陷阱的能力)。 - # -*- coding: utf-8 -*-6 g4 R, D& h/ k8 z1 C( c |2 G
- """粒子群优化算法(PSO)完整演示,零依赖可直接运行"""
. \3 E6 \6 H7 a( S F4 F - import random
* H. w3 D5 U6 S2 i+ k2 e/ V - import math: p& @* X4 D' `7 s9 w+ Y
- 2 e& ^0 @* n( w7 i3 f6 t% ? b$ Y
- def sphere(x):
3 ~' l6 _, F$ |) O$ ~ - """Sphere 函数:f(0,0,...,0)=0,单峰"""\" d+ F# U3 A, H\" {
- return sum(xi ** 2 for xi in x)
' I6 D0 u9 l6 W. _7 a, Z - 3 F2 E5 Q1 p+ h: V! x
- def rastrigin(x):1 M5 r/ K6 W\" a
- """Rastrigin 函数:大量局部最优点,全局最优仍在原点"""
0 F; y4 {) `+ U - return sum(xi**2 - 10*math.cos(2*math.pi*xi) + 10 for xi in x)
7 X0 Q5 o5 u. t% i# V w0 r
) t) H' D' J$ i- class Particle:
0 `' k! P# w$ x - """一只『鸟』:有位置、速度,记得自己飞过的最好地方"""' ]6 L8 e4 d$ E. s* J
- def __init__(self, dim, bounds):9 _5 \- V9 ?' P
- self.pos = [random.uniform(bounds[i][0], bounds[i][1]) for i in range(dim)]9 G5 R: ^# ?1 J\" t( i. ^9 @8 j
- self.vel = [random.uniform(-1, 1) for _ in range(dim)]
, N5 N; ^/ X1 \$ Z0 x - self.best_pos = list(self.pos) # 个体极值 pbest
) c' u\" h4 d8 f6 \9 f, \2 F! K) ] - self.best_val = float("inf")4 r% R) S, X6 [: K8 i
- 0 M; K6 {8 e. I8 r$ ~- t/ B
- def update_velocity(self, gbest, w, c1, c2): K* e( L. Y5 L$ D2 \' _4 _\" O
- # 速度更新:v = w*v + c1*r1*(pbest-x) + c2*r2*(gbest-x)
! Y) v, A) s# I! @3 a( n - for i in range(len(self.vel)):
5 P* v& H; r8 F7 Q f6 L# g: V - r1, r2 = random.random(), random.random()
9 @3 J+ U, T0 a! j; f: B - cognitive = c1 * r1 * (self.best_pos[i] - self.pos[i]) # 自己的经验
) A) O5 s6 Q- v. S7 [' { - social = c2 * r2 * (gbest[i] - self.pos[i]) # 群体的情报4 P; y6 c( w6 @5 c I- w: \' w! u
- self.vel[i] = w * self.vel[i] + cognitive + social
, U/ _, ~9 V7 }$ ` - / [* {* R. d% @6 m- o X\" B
- def update_position(self, bounds):! z% P1 N. r# y; ]: C
- # 位置更新:x = x + v,并夹回边界
1 W3 R5 R- I8 }2 c/ Q$ F - for i in range(len(self.pos)):, S7 {- D& x: ]# u# U
- self.pos[i] += self.vel[i]
# x6 z$ }3 q1 A4 f1 n% [, K5 R - lo, hi = bounds[i]
$ K1 A* P' E( f E( r2 y - if self.pos[i] < lo: self.pos[i] = lo
% k\" p+ M5 K& Y/ s4 ]6 Z4 Q9 a - elif self.pos[i] > hi: self.pos[i] = hi- N1 k+ G* i- }8 B' b2 L6 F. d
7 `8 @\" Z6 E3 A. U h+ q* @ y: h4 c- def pso(func, dim=2, bounds=None, n_particles=30, max_iter=100,
% l2 c u8 B& N# N - w=0.7, c1=1.5, c2=1.5, verbose_every=10):
( x! C/ y' c+ Q - """标准 PSO 主流程"""1 D6 @: U+ Y. T, _\" J0 V& B8 b
- if bounds is None:
$ [- w. n9 ?2 ^, W, |: G+ J2 K - bounds = [(-5.12, 5.12)] * dim, z\" G0 W I- w) j5 x) H3 q: P
- particles = [Particle(dim, bounds) for _ in range(n_particles)]
) Q& Q# A4 |/ k# H6 V - gbest, gbest_val = list(particles[0].pos), float("inf")
+ d4 ?) q# M+ X* L, N V3 ^ - history = []) { ^5 ^8 f# n4 X; x. F9 m
/ I& X) w$ `\" _8 O. {; ^* j& ^- for it in range(max_iter):/ m' A6 F# [, n( K
- # 第 1 步:评估每只鸟,更新 pbest 和 gbest
9 a6 y) O+ O& F9 i4 z - for p in particles:
# D2 c5 V/ _# `\" V2 i3 h: l\" q - val = func(p.pos)7 f9 t+ r- e% y
- if val < p.best_val:
6 G: J& J! c% _- E, k, ^7 w - p.best_val, p.best_pos = val, list(p.pos)
% l+ e4 Z7 E\" N1 C - if val < gbest_val:! Y9 T1 A# H) ^
- gbest_val, gbest = val, list(p.pos)
4 ^$ ?- P. Q. D y - history.append(gbest_val)' _+ ?9 l t! z/ A J( B8 c4 @. _/ i
7 H7 w( Q: }, e) [4 o4 Y1 G- # 第 2 步:更新速度和位置
1 ]7 E8 B w1 }/ G/ B - for p in particles:4 N' A! z6 o; o; L8 a' Y8 I
- p.update_velocity(gbest, w, c1, c2)
# \5 O7 n\" x. {8 K8 g - p.update_position(bounds)
$ E; |1 ?4 |3 N' K8 n4 q& @- j* J
4 f! o- B2 U$ v Z+ N- if verbose_every and (it + 1) % verbose_every == 0:- a. L- w& k7 O8 s
- print(f" 第 {it+1:3d} 代 | 全局最优值 = {gbest_val:.6f}")0 W6 r& f# T' B3 N- i
- return gbest, gbest_val, history
- p: a: P! E. M - ) n6 T* `& I+ Q8 g
- if __name__ == "__main__":
8 O- x, i) A) l5 R& H3 `, Y2 K - # 实验 1:单峰函数,看收敛速度/ A! @4 k; `' _0 w
- print("[实验 1] Sphere 函数(理论最优值 = 0)")
) ?6 z1 [3 j- K5 D; I# m1 s - best_pos, best_val, _ = pso(sphere, dim=2, bounds=[(-10, 10)]*2,4 Y; B% C6 J2 F1 D* [; K
- n_particles=30, max_iter=100)
' b6 u0 q* b2 g `6 @# i- a& v m9 Y - print(f" >>> 最优位置: ({best_pos[0]:.4f}, {best_pos[1]:.4f})")+ W1 l* M, \' Q, p2 z6 M
- print(f" >>> 最优值: {best_val:.8f}")8 A+ _3 D$ }5 Q3 ^
) Z; q& X) [# T8 i4 o' W, _- # 实验 2:多峰函数,看跳出局部最优2 W5 b- D+ Z6 W( ^5 J8 t/ X0 e
- print("[实验 2] Rastrigin 函数(局部最优遍地都是)")7 m' z; ?3 Y0 H3 P
- best_pos, best_val, _ = pso(rastrigin, dim=2,/ q/ ^9 N8 A+ ^1 _+ O$ d' b
- n_particles=50, max_iter=200)0 @ p7 A9 p$ P! M. ^' m6 U X
- print(f" >>> 最优值: {best_val:.6f}")
1 {4 j+ C; k; w/ s Q) C: [ u - . ^; d% F& Q1 Y: L7 t& X, t d
- # 实验 3:惯性权重对比(10 维)% k! a\" b& _% |* b3 P# W! h _# A
- print("[实验 3] 惯性权重 w 的影响(10 维 Rastrigin)")- A5 c0 l# Z M% Z6 s3 u; e) H: |
- for w in (0.2, 0.7, 1.2):* V4 u; O; G1 _& Z9 o6 V: M1 k1 ^
- _, val, hist = pso(rastrigin, dim=10, bounds=[(-5.12, 5.12)]*10,
* X, ^/ d& r1 \6 z9 @- K/ @; D( O - n_particles=50, max_iter=200, w=w, verbose_every=0)
8 ?1 d' U; L5 }\" ~ - 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,没有淘汰,只有一群朴素的个体,各自带着一点点记忆和一点点从众,互相分享情报,就把一片茫茫的搜索空间梳了个遍。 有人说这是对"群体智慧"最好的计算演绎:没有一只鸟知道答案,但鸟群知道。 下次遇到解不出的优化题,不妨也撒一把粒子出去——记得先广搜,后精修。
如果这篇文章帮到了你,点个「在看」让更多备赛的朋友看到
- i6 Q2 [! [/ z* U- v- x1 L6 H: L4 g4 t
+ ?; T0 h; ?' C+ G- _. T- a& R2 Y) f% d* W) s! ?; P* z
3 _& ]! E) i; T& r. i! b; f+ V0 c8 R7 I# A" Z) P! _' g% {" g
|