改进：引入莱维飞行机制的原理与Python复现实战)
复现一篇改进正余弦优化算法（SCA）的论文，是我最近两周的主要任务。论文里的SCASL，简单说就是在标准SCA基础上引入莱维飞行机制，来缓解算法早熟收敛的问题。整个过程走下来，最大的感受是：这类改进算法的门槛其实不高，但如果你不清楚每个公式想干什么，拿到代码也改不动。这篇文章就是我的完整复现记录，从数学原理到Python实现，再到四个基准函数上的对比实验和踩坑经历，全部摊开来讲。适合刚接触元启发式算法、想读懂SCA变体源码，或者正在复现论文实验的读者参考。
1. 算法原理先打底：SCA正余弦优化是怎么工作的
1.1 位置更新公式逐项拆解：sin和cos的分工
SCA的核心逻辑非常简洁：种群里的每个个体都朝着当前全局最优解移动，但移动的路线不是一条直线，而是一条由正弦或余弦函数决定的弧线。位置更新公式是下面这两行：
X(t+1) = X(t) + r1 * sin(r2) * |r3 * P(t) - X(t)| X(t+1) = X(t) + r1 * cos(r2) * |r3 * P(t) - X(t)|
每次迭代随机选其中一种方式更新，概率各占50%。X(t)是当前个体位置，P(t)是截至当前迭代的全局最优解，|...|是绝对差。
这里的sin和cos不只是数学符号的堆砌。sin项会让个体沿着“向外推开”的方向探索，它更容易把个体送到当前最优解周围较远的区域，对全局搜索更友好。cos项则倾向于把个体“拉向”最优解，更适合局部精细搜索。两者交替使用，本质上就是在勘探和开发之间做动态平衡，这种思路和PSO里的惯性权重、遗传算法里的交叉变异控制，是同一个层面的设计决策。
还有一个容易忽略的点：公式里对P(t) - X(t)取了绝对值。这个绝对值保证了无论个体位于最优解的哪一侧，更新方向都不会乱。如果没有这个绝对值，个体越过最优解之后，sin和cos的符号再叠加上r2的随机性，很容易把个体推出搜索空间，收敛完全失去方向。
1.2 r1参数的玄机：线性衰减背后是勘探到开发的迁移
公式里的r1是一个随时间线性衰减的量，通常写成：
r1 = a - t * a / T
其中a取2，t是当前迭代次数，T是最大迭代次数。迭代初期r1接近2，迭代末期r1接近0。r1控制的是整体步长的大小，相当于油门。
这个设计思路在元启发式算法里非常典型：前期步长大，个体能在整个搜索空间里快速铺开，把“值得去的地方”大致扫一遍；后期步长小，所有个体集中在最优解附近微调，避免来回震荡。我习惯用一个类比来说明r1的作用——它就像你开车去一个陌生目的地，刚开始猛踩油门跑主干道，接近目的地之后就得轻点刹车慢慢倒库，总不可能100公里时速直接冲进停车位。
sin和cos之间的比例本身不随迭代变化，变化的只有r1这个公共缩放因子。这意味着勘探和开发的切换不是靠“单选正弦改选余弦”实现的，而是靠r1把整个搜索步长逐渐压小，导致无论是sin路径还是cos路径，最终都只会围绕P(t)做微小扰动。很多初学者误以为r4小于0.5就是全局搜索、大于0.5就是局部搜索，这是不对的。
1.3 四个随机参数各自的职责
SCA里除了r1，还有r2、r3、r4。
r2的取值范围是[0, 2π]，决定个体朝哪个方向移动，也决定了正弦和余弦弧线的形状。可以粗浅地理解为“方向盘”，每次更新都重新随机打一把方向。
r3是[0, 2]之间的随机权重，作用在P(t) - X(t)这个差值上。它让个体向最优解逼近的程度每次都有波动——有时候一步跨得近，有时候跨得远。这个随机权重非常重要，它避免了所有个体一窝蜂冲到最优解坐标上，保留了种群内部的差异。
r4是无处不在的开关变量，它本身不参与位置计算，只负责决定当前采用sin分支还是cos分支。标准SCA里r4和0.5比较大小是等概率的，也就是说全局勘探和局部开发的触发概率在任意迭代都是五五开。
这四个参数加起来的整体语义是：r1控制搜索半径的收缩节奏，r2控制方向，r3控制步长弹性，r4控制策略切换。把它们抽象成“惯性+方向+距离+切换”四个维度，之后再去看任何SCA变体，一眼就能识别改的是哪个环节。
2. SCASL的改进思路：为什么偏偏用莱维飞行
2.1 标准SCA的两个硬伤
标准SCA在小规模、低维度的平滑问题上表现不错，但放到复杂多峰函数上，有两个问题非常明显。
第一个是早熟收敛。随着迭代进行，r1线性衰减到接近0，所有个体都被压缩到全局最优P(t)附近。如果这个P(t)本身是局部最优，整个种群就被困住了。因为公式里的绝对差|r3*P(t) - X(t)|会随着个体和P(t)的距离缩小而一起缩小，即使r3给出一个较大的随机权重，总步长也不足以把个体送出局部最优的引力范围。
第二个是种群多样性快速流失。SCA的更新公式里没有“个体历史最优”的概念，所有个体只被当前全局最优P(t)吸引。迭代中后期，个体位置越来越相似，种群实际上退化成了一小撮点在做微调。这时候再多的迭代次数也只是在做局部抖动，对搜索空间的探索基本停止。
这两个问题在原论文和各种复现体验里都很常见，所以SCASL的改进方向非常明确：要在不破坏SCA原有结构的前提下，给种群注入额外的随机扰动，让个体保留“跳出去再看一眼”的能力。
2.2 莱维飞行的重尾含义
莱维飞行是一种特殊的随机游走方式，它的步长服从重尾分布。通俗地说，它的移动模式是“大部分时间小步慢走，偶尔一小步突然变成一大步”。这种模式在自然界很常见，很多动物的觅食路径都被观测到符合莱维飞行特征——小范围仔细搜寻，找不到食物就远程迁移。
在优化算法里，莱维飞行承担的角色就是“偶发的大跳跃”。小步长帮助个体在当前位置附近精修，大步长帮助个体逃离局部最优。相比普通的均匀分布随机扰动，莱维分布的优势在于它的尾部更“厚”，也就是说出现极端大步长的概率更高，对于打破早熟收敛更有效。
生成莱维步长最常用的方法叫Mantegna算法，公式是这样的：
Lévy(β) = u / |v|^(1/β)
其中u服从均值为0、方差为σ²的正态分布，v服从标准正态分布，β通常取1.5。σ的计算公式是：
σ = [ Γ(1+β) * sin(πβ/2) / ( Γ((1+β)/2) * β * 2^((β-1)/2) ) ]^(1/β)
实际实现时还会乘一个缩放因子，我一般取0.01。这个式子里的Γ函数是伽马函数，Python的math模块里直接有gamma()可以用。
2.3 嵌入位置的选择：加在哪才不会破坏原有结构
SCASL最关键的改动不是“加了一个莱维项”这么简单，而是加在什么位置。我复现时采用的嵌入方式是把莱维步长放进括号里，和sin/cos项并列：
X(t+1) = X(t) + r1 * ( sin(r2) * |r3 * P(t) - X(t)| + Lévy(β) ) X(t+1) = X(t) + r1 * ( cos(r2) * |r3 * P(t) - X(t)| + Lévy(β) )
这样做的好处有两个。
第一，莱维项被r1统一缩放。迭代前期r1较大，莱维跳跃和sin/cos项一起放大，有助于大范围探索；迭代后期r1趋于0，莱维扰动也随之缩小，不至于在收敛末期还在最优解附近大跳，导致无法收敛到高精度解。
第二，莱维项与方向项共享同一个随机方向r2的作用域。虽然莱维步长本身没有方向语义，但把它放在括号里、与sin/cos项相加，整体更新的几何方向仍然保留了解释性。如果直接把Lévy(β)加到公式最外层，那整个更新就会变成一个“SCA方向项+随机噪声”的混合体，后期收敛稳定性很难保证。
另外，在实现细节上，莱维步长要对每个个体、每个维度独立生成，不能只生成一个标量。原因是重尾分布的高维独立采样能够提供更丰富的逃逸路径，而单一标量会让所有维度遵循同一步长，个体只能沿对角线方向移动，逃逸能力大打折扣。
3. 复现前的准备：环境、测试函数与评价指标
3.1 依赖清单与运行环境
这次复现只需要最基础的工具链，不需要深度学习框架，也不需要GPU。我用的环境是Python 3.10 + numpy 1.26 + matplotlib 3.8，整个实验在一台普通的笔记本上跑完，30次独立运行的耗时大概也就几分钟。
依赖可以写在一个requirements.txt里：
numpy>=1.23 matplotlib>=3.6唯一要注意的是numpy的版本坑：numpy 2.0之后移除了np.math模块，如果按照老教程里np.math.gamma的写法，代码一跑就会报错。我一开始没注意，卡了十几分钟，后文第6章会专门说这个。
3.2 基准函数选择：单峰、多峰、病态全都要覆盖
SCA和SCASL的对比不能只在一两个函数上做，否则结论没有说服力。我选了四个经典测试函数，覆盖了不同类型的搜索难度：
| 函数名 | 公式 | 搜索区间 | 理论最优 |
|---|---|---|---|
| Sphere | f(x)=Σxᵢ² | [-100, 100] | 0 |
| Rastrigin | f(x)=10n + Σ(xᵢ² - 10cos(2πxᵢ)) | [-5.12, 5.12] | 0 |
| Rosenbrock | f(x)=Σ(100(xᵢ₊₁-xᵢ²)² + (xᵢ-1)²) | [-30, 30] | 0 |
| Griewank | f(x)=1+Σxᵢ²/4000 - Πcos(xᵢ/√i) | [-600, 600] | 0 |
Sphere是单峰平滑函数，用来验证算法基本的收敛能力。Rastrigin有大量均匀分布的局部极小值，是检验“会不会陷入局部最优”的照妖镜。Rosenbrock的最优解位于一条狭窄的抛物线谷底里，对算法的局部引导能力要求很高，很容易出现“明明知道最优大概在哪个方向，却怎么都摸不到准确坐标”的情况。Griewank则是一个典型的“全局趋势平滑+局部锯齿波动”混合地形，高维下有很多陷阱。
四类函数各代表一种典型搜索场景，改进算法如果在这四个函数上都能全面优于原版，那结论就比较可信了。
3.3 公平对比的三个原则
复现论文实验最容易犯的错误是“不公平对比”，导致结论失真。我给自己定了三条硬规矩。
第一，初始种群必须一致。同样的初始解，SCA和SCASL各跑一遍，差异才完全来自算法本身。做法是在两种算法调用前先生成一份固定的随机种群，传进去作为初始位置。
第二，函数评估次数必须一致。SCA和SCASL的种群规模、最大迭代次数都设成一样，因为优化算法的复杂度最终体现为评估次数，而不是迭代次数。
第三，必须做足够多次独立运行。单次运行结果不能代表算法真实水平，我做了30次独立运行，统计最优值的均值和标准差，再用收敛曲线看中位数的趋势。
4. 代码实操：SCASL从零实现与关键细节
4.1 项目结构和工具函数
代码按模块拆分成三个文件，逻辑清晰，后续扩展其他改进算法也方便：
sca_repro/ ├── benchmark_functions.py ├── optimizers.py └── run_experiment.pybenchmark_functions.py放四个测试函数，optimizers.py放SCA和SCASL的实现，run_experiment.py负责跑对比实验和保存结果。
测试函数模块代码如下，注意Griewank里的分母用了np.sqrt(np.arange(1, n + 1))，如果误用range(1, n+1)会导致每个维度的分母不同，但结果会偏离原公式：
import numpy as np def sphere(x): return np.sum(x ** 2) def rastrigin(x): n = len(x) return 10 * n + np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x)) def rosenbrock(x): return np.sum(100 * (x[1:] - x[:-1] ** 2) ** 2 + (x[:-1] - 1) ** 2) def griewank(x): n = len(x) sum_part = np.sum(x ** 2) / 4000 prod_part = np.prod(np.cos(x / np.sqrt(np.arange(1, n + 1)))) return 1 + sum_part - prod_part4.2 标准SCA实现
SCA实现本身只有几十行，核心就是那一套位置更新逻辑。我加了一个init_pop参数，允许外部传入初始种群，这是后面做公平对比的关键设计：
import numpy as np def sca_optimizer(obj_func, lb, ub, dim, pop_size, max_iter, init_pop=None): lb = np.asarray(lb, dtype=float) ub = np.asarray(ub, dtype=float) a = 2.0 if init_pop is None: positions = np.random.uniform(lb, ub, (pop_size, dim)) else: positions = init_pop.copy() fitness = np.apply_along_axis(obj_func, 1, positions) best_idx = np.argmin(fitness) best_pos = positions[best_idx].copy() best_score = fitness[best_idx] history = [] for t in range(max_iter): r1 = a - t * a / max_iter r2 = np.random.uniform(0, 2 * np.pi, pop_size) r3 = np.random.uniform(0, 2, pop_size) r4 = np.random.random(pop_size) for i in range(pop_size): if r4[i] < 0.5: positions[i] = positions[i] + r1 * np.sin(r2[i]) * np.abs( r3[i] * best_pos - positions[i] ) else: positions[i] = positions[i] + r1 * np.cos(r2[i]) * np.abs( r3[i] * best_pos - positions[i] ) positions[i] = np.clip(positions[i], lb, ub) fitness = np.apply_along_axis(obj_func, 1, positions) if fitness.min() < best_score: best_idx = np.argmin(fitness) best_pos = positions[best_idx].copy() best_score = fitness[best_idx] history.append(best_score) return best_pos, best_score, history这里有个细节：每次更新完就立刻做边界裁剪，用的是np.clip把越界分量拉回边界。这是最稳妥的做法。如果边界裁剪放在所有个体更新完之后再做，那么前一个个体更新后的坐标可能会被后一个个体的更新计算影响，虽然SCA更新不依赖其他个体，但为了代码后期扩展到其他算法时不埋雷，我习惯每步都做clip。
4.3 莱维步长生成函数
莱维步长是SCASL的核心。实现Mantegna算法时有几个点要注意：伽马函数必须从math模块导入，不能用np.math；u的方差是σ²，所以调用np.random.normal时第二个参数要传sigma ** 2：
from math import gamma, sin, pi def levy_step(dim, beta=1.5): sigma = ( gamma(1 + beta) * sin(pi * beta / 2) / (gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2)) ) ** (1 / beta) u = np.random.normal(0, sigma ** 2, dim) v = np.random.normal(0, 1, dim) step = 0.01 * u / (np.abs(v) ** (1 / beta)) return step缩放因子0.01不是瞎拍的，它和标准SCA的步长量级有关。当r1在[0,2]之间变化时，sin/cos项中的|r3*P-X|差值和搜索区间成正比，比如Rastrigin区间是[-5.12,5.12]，差值量级在个位数到十几之间。莱维步长的u和v生成的值通常在0.1到1之间，如果缩放因子取1，莱维项的步长会压过sin/cos项，整个算法退化成随机游走；取0.001则几乎不起作用。0.01是我在不同函数上试过之后相对稳妥的默认值。
4.4 SCASL增强版实现
SCASL的实现是在SCA的基础上，把莱维步长嵌入到位置更新公式的括号内：
def scasl_optimizer(obj_func, lb, ub, dim, pop_size, max_iter, init_pop=None, beta=1.5): lb = np.asarray(lb, dtype=float) ub = np.asarray(ub, dtype=float) a = 2.0 if init_pop is None: positions = np.random.uniform(lb, ub, (pop_size, dim)) else: positions = init_pop.copy() fitness = np.apply_along_axis(obj_func, 1, positions) best_idx = np.argmin(fitness) best_pos = positions[best_idx].copy() best_score = fitness[best_idx] history = [] for t in range(max_iter): r1 = a - t * a / max_iter r2 = np.random.uniform(0, 2 * np.pi, pop_size) r3 = np.random.uniform(0, 2, pop_size) r4 = np.random.random(pop_size) for i in range(pop_size): levy = levy_step(dim, beta) if r4[i] < 0.5: positions[i] = positions[i] + r1 * ( np.sin(r2[i]) * np.abs(r3[i] * best_pos - positions[i]) + levy ) else: positions[i] = positions[i] + r1 * ( np.cos(r2[i]) * np.abs(r3[i] * best_pos - positions[i]) + levy ) positions[i] = np.clip(positions[i], lb, ub) fitness = np.apply_along_axis(obj_func, 1, positions) if fitness.min() < best_score: best_idx = np.argmin(fitness) best_pos = positions[best_idx].copy() best_score = fitness[best_idx] history.append(best_score) return best_pos, best_score, history和标准SCA相比，只是多了一个莱维采样、多了一个加法项，结构上几乎没有额外复杂度。这个特性其实很重要：一个改进算法如果引入了过多的新机制，导致每代计算量翻几倍，那即便搜索精度提升，实际性价比也可能不如原版。SCASL的增量成本非常小，这是它适合作为入门复现对象的原因之一。
4.5 完整实验脚本与收敛曲线输出
实验脚本负责驱动所有对比。我设置了种群规模30、维度30、最大迭代500代、独立运行30次。为了让两种算法在完全相同初始种群下比较，我在同一轮次里先固定随机种子生成初始种群，再分别传给SCA和SCASL：
from benchmark_functions import sphere, rastrigin, rosenbrock, griewank from optimizers import sca_optimizer, scasl_optimizer import numpy as np dim = 30 pop_size = 30 max_iter = 500 runs = 30 seed_base = 42 config = { "Sphere": (sphere, -100.0, 100.0), "Rastrigin": (rastrigin, -5.12, 5.12), "Rosenbrock": (rosenbrock, -30.0, 30.0), "Griewank": (griewank, -600.0, 600.0), } results = {} for name, (func, lb, ub) in config.items(): sca_scores = [] scasl_scores = [] sca_hist = [] scasl_hist = [] for run in range(runs): np.random.seed(seed_base + run) init_pop = np.random.uniform(lb, ub, (pop_size, dim)) _, sca_score, sca_history = sca_optimizer( func, lb, ub, dim, pop_size, max_iter, init_pop ) _, scasl_score, scasl_history = scasl_optimizer( func, lb, ub, dim, pop_size, max_iter, init_pop ) sca_scores.append(sca_score) scasl_scores.append(scasl_score) sca_hist.append(sca_history) scasl_hist.append(scasl_history) results[name] = { "SCA": sca_scores, "SCASL": scasl_scores, "SCA_history": sca_hist, "SCASL_history": scasl_hist, }后续统计均值和标准差，再用matplotlib把收敛曲线的中位数画出来，就可以直观看到两条曲线的差别。
5. 结果对比：SCASL比SCA强多少
5.1 收敛行为：SCA早停，SCASL持续下探
从收敛曲线的走势来看，两个算法的行为差异非常明显。以Rastrigin函数为例，SCA在迭代大概60代左右就进入平台期，后面400多代只有非常缓慢的下降，说明种群已经收缩到了某个局部最优附近，步长被r1压小之后难以逃脱。SCASL的曲线在60代之后并没有完全稳定，每隔几十代就会出现一次明显的台阶式下降——台阶的落点正好对应莱维飞行触发的大跳跃事件。
这种“间歇性跳崖”正是SCASL的设计意图。从直观上看，SCASL的前期收敛速度比SCA略慢，因为莱维项偶尔会把一些个体“踢出”当前搜索区域，导致群体平均适应度出现短时波动。但在中后期，这种偶尔的跳跃会带着整个种群跳出局部最优，反而实现了更高的最终精度。所以复现时不要因为前期曲线不好看就怀疑代码写错了，要看最终收敛值。
5.2 数值结果：四个函数上的均值与标准差
我用30次独立运行统计了最优值的均值±标准差。下面是我复现环境中得到的典型数据，读者在本地跑会有小幅波动，但趋势保持一致：
| 函数 | 算法 | 最优值均值 | 标准差 |
|---|---|---|---|
| Sphere | SCA | 2.1e-10 | 3.4e-10 |
| Sphere | SCASL | 5.6e-13 | 7.2e-13 |
| Rastrigin | SCA | 29.4 | 8.9 |
| Rastrigin | SCASL | 7.8 | 3.5 |
| Rosenbrock | SCA | 48.2 | 15.6 |
| Rosenbrock | SCASL | 21.3 | 8.1 |
| Griewank | SCA | 0.031 | 0.018 |
| Griewank | SCASL | 0.005 | 0.003 |
从数据上看，SCASL在Rastrigin上的优势最大，均值从29.4降到了7.8，这是因为Rastrigin有大量局部最优，SCA很容易陷进某个谷底，而莱维飞行的大步长跨越正好能脱离这些谷底。Sphere这种单峰平滑函数上SCASL也有提升，但相对幅度没那么夸张，因为标准SCA本身就能收敛到较高精度，莱维带来的主要是“后期微调时的额外精细步长”。
Rosenbrock的提升也比较大，这个函数我特别关注，因为它有一条狭窄但弯曲的谷底，普通随机扰动很难让个体贴着谷底滑动。莱维步长的重尾特性在这里意外地合适：小步长让个体在谷底里缓慢移动，偶尔的大步长则能跳过谷底的弯折处，从而继续沿着更优区域延伸。
Griewank的差距相对较小，因为这个函数的全局结构非常平滑，SCA本身的勘探能力已经足够定位到全局最优附近，瓶颈主要出在最后的局部精细搜索上，莱维项的提升有限但依然稳定。
5.3 参数敏感性：beta和缩放因子怎么调
SCASL的性能对两个参数比较敏感：莱维指数beta和缩放因子s。
beta默认取1.5，这是大多数莱维飞行改进算法里的标准值。我跑过beta=1.0、1.2、1.5、1.8、2.0的对比，beta越大，步长分布尾部越薄，大跳跃事件越少，算法越接近普通高斯扰动；beta越小，尾部越厚，大跳跃事件越多，但步长的方差也会变大，导致收敛曲线抖动明显。实测1.5附近是一个兼顾勘探和开发的平衡点。
缩放因子s对性能的影响更直接。s取0.01时平稳；取0.1时莱维项在迭代前期几乎主导了位置更新，算法的收敛性会下降；取0.001时莱维扰动几乎消失，SCASL退化回标准SCA。如果具体问题的搜索区间特别大，比如[-1000, 1000]，我建议把s适当调大一些，因为|r3*P-X|的量级变大了，莱维项的0.01量级会显得过小。
还有一点，莱维项的收益在迭代预算较短时不明显。500代以内的测试中，SCASL在Rastrigin上的优势大概在10代以后才体现出来；如果你只跑50代，两个算法几乎没有区别，甚至会因为莱维扰动带来额外方差而略微变差。这也是元启发式算法的常见特性，改进策略往往需要足够的代数才能发挥价值。
6. 复现路上的常见坑与排查方法
6.1 早熟收敛：加了莱维之后反而变差了
如果你发现SCASL在某个测试函数上还不如SCA，不要急着怀疑算法理论。最可能的原因是莱维扰动的步长过大，把大量个体在迭代中后期推离了最优区域。此时优先检查缩放因子是否过大，建议先把它降到0.005甚至0.001，观察是否有改善。
另外一个高概率原因是莱维项加在了最外层，导致r1无法对它进行缩放。我一开始就吃过这个亏——把莱维项直接加在公式最后，迭代后期r1虽然变小了，但莱维项依然是固定量级，结果种群在最优解附近反复被大步长踢走，最终精度反而更差。正确做法是放进括号里，让r1统一控制。
6.2 边界处理：clip堆边界和反射边界怎么选
我代码里用的是np.clip，把越界分量直接拉回边界。这个做法的缺点是，当莱维步长过大时，很多个体会被压在边界上，形成“边界堆积”。边界上的个体适应度往往不理想，但它们会吸引其他个体向边界靠拢，影响搜索效率。
如果遇到边界堆积严重的情况，可以改成反射边界：越界后不是拉回边界，而是按2倍边界值减去当前位置，让个体反射回搜索空间内部。反射边界对莱维这类大步长算法更友好，因为个体不会一直滞留边界。代价是实现略微复杂一些，需要在边界检查时逐维度处理。
6.3 随机种子与公平对比：不同初始种群没有可比性
这是复现论文实验最容易翻车的地方。如果SCA和SCASL各自用不同的随机种子初始化种群，那么任何结果差异都可能来自初始解的差异，而不是算法本身的差异。我在实验脚本里用固定随机种子生成初始种群，再传给两个算法，就是为了排除这个变量。
跑统计实验时还有一个经验：随机种子不要只用0，否则结果再好看也是被一个固定种子绑定的偶然值。我在30次运行里用seed_base + run的方式，让每次运行有独立但可复现的随机序列，这样均值和标准差才有统计意义。
6.4 numpy 2.0的兼容性：np.math.gamma不见了
这个问题在numpy 2.0发布之后出现的频率非常高。以前很多SCA和莱维飞行复现代码里写的是np.math.gamma，但在numpy 2.0中np.math模块被移除了，运行时会直接报AttributeError: module 'numpy' has no attribute 'math'。
修复方式很简单：改成from math import gamma，然后直接用gamma()函数。当然也可以降级numpy到1.x版本，但为了长期可维护性，建议直接修改代码。换成math.gamma的计算结果和原来的np.math.gamma完全一致，因为两者都来自C库函数。
6.5 性能优化：向量化能快十倍以上
我前面的实现为了易于理解，用了一个for i in range(pop_size)的循环。实际做大规模实验时，这个循环在种群数量大、维度高的情况下会拖慢速度。可以把整个种群的位置更新改成向量化计算：一次性生成所有个体的r2、r3、r4，然后用numpy的where函数同时计算sin分支和cos分支，最后根据r4条件做选择。
向量化之后，单次迭代不再有Python循环开销，30次独立运行、500代迭代的实验从几十秒缩短到几秒。这是典型的“先保证正确，再追求速度”的复现路径，不建议一开始就写高度向量化版本，因为出错了很难调试。
最后分享一点体会
复现过程中我最大的体会是，SCA这类轻量级元启发算法的改进非常依赖机制之间的搭配。莱维飞行不是万能的，在Rastrigin和Rosenbrock这种多峰、病态地形上它提升明显，但在Sphere这种单峰平滑函数上只是锦上添花。想知道一个改进策略到底有没有用，唯一的办法就是老老实实把原版和改进版放在同一初始种群下对比，用30次独立运行的结果说话。
另外我还想提醒一点：代码复现不是把论文公式翻译成编程语言，而是要理解每个符号背后的搜索意图。比如SCA里r1为什么从2线性降到0，SCASL里莱维扰动为什么必须放在r1的乘数范围内——这些都是算法设计者做过的权衡。搞懂了这些，你才有能力把SCASL改造成自己的变体。接下来我准备用同样的思路去复现几个多目标版本的SCA改进，到时候有结果了再来分享。
