1. 为什么我决定复现这个“三件套”改进
蝴蝶优化算法（BOA）在群体智能算法里不算冷门，它靠“气味浓度”来引导个体位置更新的机制很特别，代码写起来也比粒子群简单。但这两年我陆陆续续看了不少BOA改进文章，发现一个普遍现状：直接在原始BOA上跑基准函数，结果往往不够看，收敛慢、后期不动、局部最优里出不来。这次我复现的这篇改进文章，标题里就把组合拳写在明面上：Circle混沌初始化种群、非线性因子w/p/r、融合正余弦算法。说白了就是分别在种群初始化、参数控制、位置更新三个层面各补一刀。这篇博文就把我从零复现这个组合算法的完整过程写出来，包括每个改进点背后的原理、关键公式、可运行的Python核心代码，以及我实际踩过的一堆坑。适合正在做群智能优化算法改进、打算复现论文又不清楚细节的人看。
1.1 原始BOA的收敛逻辑，以及它到底差在哪
先过一遍原始BOA的核心逻辑。每只蝴蝶的位置就是一个候选解，个体i会释放“香味”，香味浓度的计算公式是 f_i = c * I_i^a，c是感知因子，a是幂指数，I_i是刺激强度（通常取适应度相关的值）。搜索过程按概率p切换全局探索和局部开采。全局探索时，个体朝当前全局最优位置靠近，常见更新式是 x_i^{t+1} = x_i^t + f_i * (r^2 * g^* - x_i^t)，其中g^*是当前全局最优位置；局部开采时，随机选两个个体j和k，更新式为 x_i^{t+1} = x_i^t + f_i * (r^2 * x_j^t - x_k^t)。这里r是[0,1]的均匀随机数。
只要把这个公式拆开就能发现问题。c、a、p全是固定值，意味着整个迭代过程中步长性质和搜索偏好完全不变；初始种群如果只用普通均匀随机数生成，在高维空间里大概率出现局部扎堆、大片区域空白的情况。更尴尬的是，后期算法已经陷入局部最优了，步长还是以前那个量级，想跳跳不出来，想精细搜索又太粗。所以原始BOA在单峰函数上表现尚可，一到Rastrigin、Griewank这类多峰函数，经常跑出糟糕的最终精度。
1.2 三个改进点为什么能形成互补
看清楚原版的三个短板，标题里的三个改进就很好理解了。Circle混沌初始化是解决开局问题的：混沌序列有遍历性和低相关性，能让种群在搜索空间里分布更均匀，避免初始解过度集中在某一小块区域。非线性因子w/p/r是解决参数僵化问题的：w控制位置更新时的惯性保留比例，p控制全局搜索和局部搜索的切换节奏，r控制气味强度对适应度差异的敏感度，三个参数都随迭代次数动态变化，前期偏探索、后期偏开发。融合正余弦算法则是在位置更新算子这一层做增强：正弦余弦震荡项既能提供大范围的跳跃能力，又能在后期保留小步长的精细搜索，补上了BOA局部搜索多样性不足的毛病。
这三个模块作用在不同环节，不重复，组合起来逻辑是通的。但复现的时候不能一股脑全塞进去，我建议先单独跑通每个模块，再逐步叠加，否则出了问题你根本不知道是哪个模块引起的。后文我会按这个思路把每个模块拆开讲。
2. 改进点逐个拆解：从Circle混沌初始化到SCA融合
2.1 Circle混沌初始化种群：均匀撒网，避免开局“挤成一团”
混沌初始化的目的不是让代码看起来高级，而是真能改善初始种群在解空间中的覆盖程度。Circle混沌映射的迭代式很简单：
x_{k+1} = mod( x_k + b - (a / (2π)) * sin(2π * x_k), 1 )
常用参数a=0.5，b=0.2。输出的序列被mod到[0,1)区间，理论上具有较好的遍历性。和Logistic混沌、Tent混沌相比，Circle的实现只多了个正弦项，但分布形态更均匀，不容易出现局部过密的情况。
实际初始化时，我对每个维度各生成一条独立的混沌序列，而不是把所有维度拼成一条。这样做是为了避免维度之间的相关性问题，否则种群会落在一个低维子空间里，等于变相缩小了搜索范围。代码如下：
import numpy as np def circle_init(N, dim, lb, ub): chaos = np.zeros((N, dim)) for d in range(dim): x = np.random.rand() + 0.1 seq = [] for _ in range(N + 50): x = np.mod(x + 0.2 - (0.5 / (2 * np.pi)) * np.sin(2 * np.pi * x), 1.0) seq.append(x) # 丢弃前50个暂态点，避免混沌序列初期的异常波动 chaos[:, d] = np.array(seq[50:50 + N]) population = lb + chaos * (ub - lb) return population, chaos这里有几个细节一定要留意。第一，每个维度的混沌序列初值x最好不要取0、0.25、0.5这类特殊值，否则序列可能直接进入周期点；第二，生成序列时多产生50个点再丢掉，是为了绕过混沌映射的暂态过程；第三，混沌序列本身对浮点精度敏感，如果发现N比较大时序列出现重复，可以改用np.longdouble提升中间精度。我复现时一开始没丢暂态点，结果前10个个体几乎挤在同一条线上，后来加了这个处理才正常。
2.2 非线性因子w、p、r：给BOA的三个关键旋钮加上“时间轴”
原始BOA里的c、a、p都是固定常数，改进文章通常会把其中一部分参数改成随迭代变化的非线性函数。标题里说的w、p、r，在我复现的这篇文章背景下，对应的是惯性权重w、转换概率p、气味强度幂指数r（有的文章里也叫a）。
先说w。惯性权重的经典用法是放在位置更新式前面的系数，保证个体保留一部分原有运动信息。我们这里用的是抛物线递减：
w(t) = w_min + (w_max - w_min) * (1 - (t/T)^2)
w_max取0.9，w_min取0.4。前期w大，个体不容易被最优位置带偏，能维持较高的探索能力；后期w小，步长跟着收缩，更容易在当前最优附近做精细搜索。选择抛物线而不是线性递减，是因为前段下降更快、后段变化更缓和，更符合“前期快速探索、后期慢慢精修”的需求。
然后是p。p在BOA里是全局搜索与局部搜索的切换概率，原始版本固定为0.8，有点偏高，导致局部开采机会不足。改进后我用的公式是：
p(t) = p_max - (p_max - p_min) * (1 - cos(πt / (2T)))
p_max=0.8，p_min=0.2。前期p接近0.8，大多数个体做全局搜索；后期p接近0.2，把重心转向局部搜索。余弦形式的好处是中段过渡平滑，不像线性切换那样容易产生阶段突变。
最后是r。r直接影响气味浓度 f_i = c * I_i^r。r较小的时候，不同个体的气味强度差异被压缩，种群能保持多样；r较大的时候，最优个体的气味被放大，其他人更容易被吸过去。所以我的取值是：
r(t) = 0.1 + 0.7 * sin(πt / (2T))
前期接近0.1，后期慢慢涨到0.8。这个参数一旦取反或者变化太快，算法会在中期突然“抱团”，收敛曲线直接掉进平台期。调参的时候可以优先检查r的变化曲线有没有过大。
三个参数的变化趋势我整理成下面这张表，方便复现时对照：
| 参数 | 初始值 | 终止值 | 变化形式 | 主要作用 |
|---|---|---|---|---|
| w | 0.9 | 0.4 | 抛物线递减 | 控制步长和历史信息残留 |
| p | 0.8 | 0.2 | 余弦递减 | 从全局搜索切换到局部搜索 |
| r | 0.1 | 0.8 | 正弦递增 | 增强最优个体的气味吸引强度 |
2.3 正余弦算法融合：用振荡项补上局部搜索的多样性
正余弦算法（SCA）的核心思想是用正弦和余弦函数产生的振荡方向作为搜索步。标准更新式是这样：
X_new = X + a_sca * sin(r2) * |r3 * X_best - X|， 当 r4 < 0.5 X_new = X + a_sca * cos(r2) * |r3 * X_best - X|， 当 r4 >= 0.5
其中a_sca从2线性递减到0，r2是[0,2π]的随机角度，r3是[0,2]的随机系数，r4是[0,1]的随机数。正弦项负责向远离最优点的方向试探，余弦项负责向靠近最优点的方向细化，两者配合就能在局部区域形成一种震荡搜索。
融合BOA的时候，我采用的方法是：在全局搜索分支里同时生成两个候选位置。一个来自BOA原始公式加惯性权重，记为v_boa；另一个来自SCA更新式，记为v_sca。然后比较两者的适应度，哪个好就保留哪个。这种“贪婪选择”比直接加权平均更稳，因为它保证了每代更新至少不会比BOA原本的候选更差。
融合后的全局更新逻辑可以这样理解：BOA提供的是基于气味梯度的方向性搜索，SCA提供的是基于三角振荡的随机性试探。两者交替生效，能让种群在飞向全局最优的同时保留跳出局部最优的可能性。
3. 实操实现：从Circle种群初始化到主循环完整代码
3.1 实验配置与基准函数选择
复现这类算法，实验配置必须写清楚，否则后续对比没有意义。我的环境是Python 3.9 + NumPy 1.24，基准函数选了五个有代表性的：Sphere是单峰函数，用于检验收敛精度；Rastrigin有大量局部极小点，用于检验跳出局部最优的能力；Griewank维度间存在乘积耦合，综合性比较强；Ackley坡度平缓且多峰，是很多改进算法的“照妖镜”；Schwefel 2.26的全局最优离搜索子空间边缘很远，很容易陷入次优解。
参数配置我用的是30个种群个体、30维、500次迭代，每个函数独立跑30次，记录每一代的最优值和最终最优值，最后统计均值和标准差。为什么跑30次？因为群智能算法随机性很强，跑一次得出的结论和掷骰子没区别。具体参数如下：
| 参数 | 取值 |
|---|---|
| 种群规模N | 30 |
| 维度D | 30 |
| 最大迭代次数T | 500 |
| 独立运行次数 | 30 |
| 感知因子c | 0.01 |
| w | 0.9→0.4（抛物线递减） |
| p | 0.8→0.2（余弦递减） |
| r | 0.1→0.8（正弦递增） |
| a_sca | 2→0（线性递减） |
3.2 主循环与位置更新代码
整个算法的主循环实现如下。我把Circle初始化、刺激强度转换、动态参数、SCA融合和边界反射都整合在一个函数里，便于对照复现。
import numpy as np def reflect(pos, lb, ub): pos = np.array(pos, dtype=float) pos = np.where(pos < lb, 2 * lb - pos, pos) pos = np.where(pos > ub, 2 * ub - pos, pos) return np.clip(pos, lb, ub) def stimulus(fit_value): # 把最小化问题转换为“刺激强度越大越好” return 1.0 / (1.0 + fit_value) def boa_improved(obj_func, N=30, D=30, T=500, lb=-100, ub=100): pop, _ = circle_init(N, D, lb, ub) fitness = np.array([obj_func(p) for p in pop]) best_idx = np.argmin(fitness) g_best = pop[best_idx].copy() g_best_f = fitness[best_idx] c = 0.01 w_max, w_min = 0.9, 0.4 p_max, p_min = 0.8, 0.2 r_start, r_end = 0.1, 0.8 history = [] for t in range(T): w = w_min + (w_max - w_min) * (1 - (t / T) ** 2) p = p_max - (p_max - p_min) * (1 - np.cos(np.pi * t / (2 * T))) r = r_start + (r_end - r_start) * np.sin(np.pi * t / (2 * T)) a_sca = 2.0 * (1.0 - t / T) for i in range(N): f_i = c * (stimulus(fitness[i]) ** r) if np.random.rand() < p: # BOA 全局搜索候选 ra = np.random.rand() v_boa = w * pop[i] + f_i * (ra * ra * g_best - pop[i]) # SCA 候选 r2 = np.random.uniform(0, 2 * np.pi) r3 = np.random.uniform(0, 2) if np.random.rand() < 0.5: v_sca = pop[i] + a_sca * np.sin(r2) * np.abs(r3 * g_best - pop[i]) else: v_sca = pop[i] + a_sca * np.cos(r2) * np.abs(r3 * g_best - pop[i]) # 贪婪选择更优的候选 if obj_func(v_sca) < obj_func(v_boa): new_pos = v_sca else: new_pos = v_boa else: # BOA 局部搜索候选 choice = [idx for idx in range(N) if idx != i] j, k = np.random.choice(choice, 2, replace=False) ra = np.random.rand() new_pos = w * pop[i] + f_i * (ra * ra * pop[j] - pop[k]) new_pos = reflect(new_pos, lb, ub) new_fit = obj_func(new_pos) if new_fit < fitness[i]: pop[i] = new_pos fitness[i] = new_fit best_idx = np.argmin(fitness) g_best_f = fitness[best_idx] g_best = pop[best_idx].copy() history.append(g_best_f) return g_best, g_best_f, np.array(history)这里有个很关键的小经验：不要把最小化问题的原始函数值直接当成刺激强度I。BOA原始文献更适合最大化问题的语境，因为香味越浓的个体吸引力越强。复现时如果直接代入5000、6000这种量级的函数值，f_i会非常大，一步就把所有个体震到搜索空间外面去。所以代码里我用1/(1+fitness)做了转换，把刺激强度压到(0,1)区间，同时保持了“适应值越小刺激越强”的方向。
3.3 如何判断复现结果是否靠谱
跑完代码之后，不要只看最终一个数字，要记录下来每次独立运行的最优值和收敛曲线。我通常会保存30次运行的逐代最优值，然后做三件事。第一，画平均收敛曲线，Y轴用对数坐标，因为群智能算法的收敛曲线前期下降快、后期变化慢，线性坐标下后半段基本是一条直线，看不出差异。第二，输出最终最优值的均值和标准差，标准差太大说明算法不稳定。第三，把改进版和原始BOA的结果做Wilcoxon符号秩检验，看改善是否显著，这比“我肉眼看着图更好”有说服力得多。
还有一个容易忽略的点：做对比时，原版BOA里也要保证种群规模、迭代次数和评估次数完全一致，否则对比结果没有意义。如果发现改进版在个别函数上不如原版，这很正常，没有哪个算法能在所有函数上全面碾压，关键是整体胜出。
4. 复现路上的坑：常见问题与排查技巧
4.1 Circle混沌序列退化：为什么生成的点“不够混沌”
我复现时遇到过一个诡异现象：用Circle混沌初始化30个个体，画出来只有五六个位置。检查发现是混沌序列进入了短周期。原因之一是没有丢弃暂态点，另一个原因是初值x取到了特殊值。Circle映射对初值和参数比较敏感，一旦序列进入周期，整条链就废了。快速排查方法是对每个维度的序列统计唯一值数量：
_, chaos = circle_init(30, 30, -100, 100) for d in range(3): unique_cnt = len(np.unique(np.round(chaos[:, d], 8))) print(unique_cnt)如果唯一值数量明显小于种群规模，就说明该维度混沌序列退化。解决办法有三个：初值加一个随机偏移，丢弃前50个暂态点，必要时在生成序列时每隔一个点取一个样本，降低相邻点相关性。这三个操作都不改变混沌的本质，只是让实际使用的序列更稳定。
4.2 气味强度爆炸：最小化问题上f_i异常大
原版BOA公式里的f_i = c * I^a，在最小化基准函数上直接套用非常容易翻车。举个例子，Sphere函数在d=30、初始范围[-100,100]时，一个随机位置的函数值可能超过10000。如果r又随着迭代涨到0.8左右，c=0.01，那么f_i = 0.01 * 10000^0.8 ≈ 2.5，这还勉强能接受，但在Rastrigin或Ackley上初始值更大，f_i可能到几十上百，新位置直接就飞出天边。所以我在代码里用1/(1+fit_value)把刺激强度压缩到(0,1)，再加上reflect反射边界处理，这个问题才彻底消失。如果你不想改变气味强度语义，也可以把每个个体的刺激强度除以当前种群最大值做归一化，效果类似，但要注意这种归一化会引入额外的随机性波动。
4.3 正余弦扰动过大，收敛曲线不但不降反而震荡
加入SCA之后，最容易出现的问题就是收敛曲线在后期来回震荡。SCA的a_sca初值如果保持2.0，前期扰动确实不小，但到了后期如果a_sca衰减不够快速，正弦余弦项还会产生较大的位移，把已经接近最优的个体又弹开。解决这个问题有几条路：把a_sca的线性衰减改成指数衰减，也就是a_sca = 2 * exp(-3 * t / T)；或者只在迭代前80%的阶段使用SCA，最后20%完全交给BOA原始局部搜索；再或者像我代码里那样，对所有越界位置做反射而不是简单裁剪。反射比裁剪好在哪？裁剪会把大量个体压到边界上，形成一个虚假的边界聚集点，反射则能把越界的个体弹回搜索空间内部，保留了多样性。实测下来反射让Ackley函数上的最终精度高了一个数量级左右。
4.4 改进后仍然打不过原版BOA
如果你的复现结果里改进版并没有明显优于原始BOA，先别急着怀疑算法有问题。我踩过的最深一个坑是：原版BOA的随机初始化和改进版用同一个随机种子，但原版可能会偶然生成一组非常“走运”的初始解，导致它在某个函数上胜出。所以必须多次独立运行。其次要检查p和w的变化范围。p如果后期降得太低，比如降到0.05，几乎所有个体都去做局部搜索，一旦陷入局部最优就没有拯救机会了，因此p_min不要低于0.2。最后，建议做消融实验，单独测试Circle初始化、非线性因子、SCA融合各自的效果。我实际做下来发现，Circle初始化在大部分函数上的提升是稳定的，非线性因子则在收敛速度上收益明显，而SCA融合如果参数调得不好，反而是负优化。找到负优化点之后，针对性地调整那个模块的强度，比盲目加参数更有效。
4.5 问题速查表
| 症状 | 可能原因 | 处理办法 |
|---|---|---|
| 初始化位置大量重复 | 混沌序列退化或初值特殊 | 丢弃暂态点、加随机偏移、提升浮点精度 |
| 第一次迭代飞越边界 | f_i过大或步长不受控 | 刺激强度归一化、使用反射边界 |
| 后期收敛曲线震荡 | a_sca衰减太慢或SCA权重过大 | 指数衰减a_sca、限制SCA使用阶段 |
| 曲线停在一个平台不下降 | p后期过低或r变化过慢 | p_min调到0.2以上、加快r递增速度 |
| 改进效果不稳定 | 独立运行次数过少 | 至少30次运行、统计均值和方差 |
5. 复现之外的一些经验
整套流程跑下来，我最大的体会是：算法复现不是把公式抄进代码就完事，而是每个改进点背后的逻辑都要能对应到算法实际发生的环节上去。Circle混沌初始化管的是“起点分布”，非线性因子管的是“过程节奏”，正余弦融合管的是“位置更新方式”。三个模块之间如果参数衔接不好，很容易出现模块间互相干扰。
我个人的习惯是每加一个改进模块就做一次A/B测试，比如先用Circle初始化替换随机初始化，其他部分保持原版，记录一次结果；再加非线性因子，再记录一次；最后加SCA。这样一旦某个环节结果异常，立刻能定位。另一个小技巧是，把每次运行的随机种子、参数配置、最终适应度都存成CSV，别只记一个平均值，否则后面想回溯某个异常结果根本无从下手。复现算法这种事，功夫都花在实验设计和数据记录上，真不是代码本身有多难。希望这篇记录能让你在复现BOA改进时少走我走过的弯路。
