如果你最近在折腾群体智能优化算法大概率已经把PSO、GWO、WOA翻来覆去跑过好几轮了。我拿到白鲸优化算法BWO时第一反应是又来一个鲸鱼亲戚。但仔细把原始论文读完之后我反倒觉得这个算法有点意思探索、捕食、鲸落三阶段的结构非常清晰代码也好写。真正跑起来之后问题也暴露了——在Rastrigin这类多峰函数上它很容易陷进局部极小值。后来我试着把准反向学习和旋风觅食两个算子叠进去做了一版改进算法再用MATLAB跟PSO、GWO、WOA和原版BWO对比效果确实有可感知的提升。这篇就把整个思路、MATLAB实现和对比分析的细节完整记录下来给正在做算法改进、毕设或者想拿新算法刷CEC结果的朋友做个参考。1. 整体设计与改进思路1.1 原始BWO的三阶段机制与弱点白鲸优化算法是2022年提出的元启发式算法模拟白鲸的游泳、捕食和“鲸落”三种自然行为。它的位置更新不是单一公式走到底而是由平衡因子 Bf 控制阶段切换当 Bf 大于0.5时白鲸群体自由游动执行全局探索当 Bf 小于等于0.5时进入捕食阶段个体向最优解附近聚集此外还有一个鲸落阶段按一定概率让部分个体重新随机初始化模拟白鲸死亡后尸体沉入海底、为生态系统提供养分的机制。这个设计比很多“换皮粒子群”要完整但实际跑下来有两个明显的软肋。第一探索阶段过于随机它依赖随机个体之间的位置差没有方向性收敛速度偏慢。第二捕食阶段本质上是一个线性趋优的过程局部搜索半径不会动态收缩也不具备螺旋扫描能力一旦陷入某个局部区域很难依靠自身机制跳出来尤其是在Rastrigin、Griewank这类多峰函数上容易提前停滞。所以我对原始BWO的改进思路也很直白不推翻它的三阶段框架只做两处外科手术。一是在种群初始化和迭代后期对最优个体施加准反向学习扰动提升初始多样性和后期逃逸能力二是把捕食阶段的位置更新算子替换成具有螺旋收缩性质的旋风觅食算子让算法在最优点附近做更精细的搜索。两处改动都不复杂但叠加之后的效果会互相放大。1.2 为什么选准反向学习和旋风觅食这两个算子反向学习是一种成熟的种群增强手段基本思想是如果当前个体在搜索空间某处那么它的“镜像点”往往也有很高的概率靠近全局最优于是同时评估原个体和反向个体择优保留能够显著加快收敛。但标准反向学习有个问题反向点很容易落在搜索区间边缘甚至越界导致大量个体堆在边界上多样性反而受损。准反向学习是对反向学习的一种修正它不直接取镜像点而是在区间中点和反向点之间随机取一个位置。你可以把它理解为“镜像点和中心点之间做一次插值”这样生成的点几乎一定落在搜索区间内部既保留了反向学习的探索优势又不会产生边界堆积。这个性质很关键尤其适合用作粒子群类算法的初始化策略。旋风觅食则是我针对捕食阶段开发的一个局部搜索算子。它模拟白鲸捕猎时围绕猎物形成漩涡、在旋转中逐渐收拢的行为数学上是一个半径不断缩小的螺旋更新过程。相比原始BWO捕食阶段里那种“朝最优解线性挪动”的做法螺旋轨迹能够以更平滑的方式扫描最优解的邻域不会一下子冲过头也不会原地踏步。选择这两个算子还有一个实际考虑它们互相之间不存在参数耦合代码实现简单一个是种群初始化工具一个是局部搜索工具分工明确。改动后的算法仍然保持BWO原有的三阶段骨架做对比实验时也更容易解释每一个改进点带来的增益。2. 改进机制的数学原理与参数设计2.1 准反向学习的数学描述与边界处理设当前个体为 x搜索区间下界为 lb上界为 ub。标准反向点的计算公式为x lb ub - x这是一个关于搜索空间中心的对称点。准反向点不是直接取 x而是在区间中点 (lb ub) / 2 与反向点 x 之间随机取一个位置写成x_q (lb ub) / 2 r * (x - (lb ub) / 2)其中 r 是 [0,1] 之间的随机数。当 r 接近0时准反向点靠近区间中点当 r 接近1时它接近标准反向点。由于中点到反向点的连线整体位于搜索区间内所以即便出现极端情况准反向点也几乎不会越界。我在项目中用准反向学习做了两件事。第一件是在种群初始化阶段先随机生成 N 个个体然后为每个个体生成对应的准反向个体把 2N 个个体全部计算适应度排序后保留最优的 N 个作为初始种群。这相当于用两倍初始化开销换取更好的起点分布。第二件是在迭代过程中每隔一定代数对当前全局最优个体执行一次准反向扰动并和原始最优位置比较如果扰动后的位置更好就替换否则保持不变。这个操作成本极低但能在算法陷入局部极值时提供一个跳出机会。需要注意一个细节初始化阶段虽然多算了 N 个个体但消耗的函数评估次数要计入总预算。在公平对比实验中不能说自己“迭代了500次”而要说“最大函数评估次数 MaxFEs 为30000”否则改进算法占便宜实验结论站不住脚。2.2 旋风觅食算子的数学模型与收缩策略旋风觅食算子的设计思路来自螺旋动力学。设当前白鲸个体为 x_i全局最优位置为 x_best定义距离向量dist x_best - x_i那么新的位置按如下方式更新x_new x_best rho * cos(2pit) * dist levy其中 t 是归一化进度从0到1递增rho 是收缩半径系数我取rho (1 - t)^2也可以理解成越到算法后期白鲸围绕最优点转圈的半径越小搜索越精细。cos(2pit) 提供了旋转分量让个体不是直线扑向最优解而是绕着一个不断缩小的螺旋靠近。levy 是一个随机扰动项采用莱维飞行生成用于维持一定的逃逸能力避免所有个体都掉进同一条螺旋线里。这里有一个需要手动控制的参数莱维扰动的衰减系数。我在实验中使用alpha exp(-5 * t)也就是说前期扰动较大算法还有一定的全局搜索能力后期扰动迅速衰减保证收敛精度。如果 alpha 衰减太慢后期最优解附近会一直有较大的随机抖动导致收敛曲线末端出现毛刺无法得到高精度解如果衰减太快又会让算法过早进入纯局部搜索容易卡在局部极值。这个算子虽然名字叫“旋风”但核心不是简单地套用圆形轨迹而是“旋转 收缩 随机扰动”三者的结合。旋转提供邻域扫描收缩保证逐步逼近随机扰动防止完全丧失探索能力。我实际测试下来这套组合比单独用高斯扰动或者柯西扰动要稳得多。2.3 改进BWO的完整算法流程改进后的算法保留原始BWO的阶段划分但捕食阶段使用旋风觅食算子替换原来的线性更新。整体流程可以写成初始化种群随机生成 N 个个体并计算对应的准反向个体择优保留 N 个评估种群适应度记录全局最优位置和最优值进入主循环计算平衡因子 Bf B0 * (1 - t / T)如果 Bf 0.5执行探索阶段的位置更新否则进入捕食阶段对每个个体执行旋风觅食更新按鲸落概率 Wf 判断是否触发鲸落阶段对部分个体重新随机初始化每隔固定代数对全局最优个体执行一次准反向扰动判断是否达到最大函数评估次数若未达到则返回第3步。这个流程里面探索阶段的公式可以直接沿用原始BWO因为改进的重点在捕食阶段。如果完全抛弃原始BWO的探索机制实验结果就很难归因到“旋风觅食”这个改进点上审稿人或者导师也会追问你到底改动在哪里。所以我的建议是每次只改一个变量其他环节尽量保持原样。另外鲸落概率 Wf 的初始值不需要调得太高一般取 0.1 左右即可。白鲸优化算法的鲸落阶段本来就是个低频事件设太高容易让种群频繁重置收敛曲线会变得很难看。3. MATLAB实现从公式到可跑通的代码3.1 文件结构与主函数骨架我习惯把一个优化算法的实现拆成几个独立函数方便调试和复用。这个项目的文件结构如下IBWO.m主函数负责参数配置、循环调度和结果输出QOBL.m准反向学习函数生成种群的准反向个体CycloneForaging.m旋风觅食算子LevyFlight.m莱维飞行的随机数生成fobj.m目标函数集合通过参数切换测试不同的基准函数RunExperiments.m批量实验脚本跑多个算法、多个函数、多次独立运行。主函数的参数配置我一般写成这样N 30; % 种群规模 T 1000; % 最大迭代次数 B0 0.4; % 平衡因子初始值 Wf 0.1; % 鲸落概率 lb -100; % 搜索空间下界 ub 100; % 搜索空间上界 funcName Sphere; % 目标函数名称 rng(2024); % 固定随机种子保证实验可复现需要说明的是最大函数评估次数是 N * T也就是30000次。初始化阶段额外生成的 N 个准反向个体如果在初始化阶段就计算了适应度那么也要计入总评估次数。更严格的做法是在循环里用一个 FEs 计数器每次调用 fobj 就累加循环条件判断 FEs 是否达到 MaxFEs而不是死板地判断迭代次数。3.2 准反向学习函数实现准反向学习的MATLAB实现非常短核心就是矩阵运算不要为了图省事写for循环否则在高维问题上会非常慢。function Xq QOBL(X, lb, ub) % X: 当前种群大小为 N x D % Xq: 准反向种群大小与 X 相同 N size(X, 1); D size(X, 2); lbM repmat(lb, N, 1); ubM repmat(ub, N, 1); Xo lbM ubM - X; % 标准反向点 mid (lbM ubM) / 2; % 区间中点 r rand(N, D); Xq mid r .* (Xo - mid); % 中点与反向点之间取随机位置 Xq max(min(Xq, ubM), lbM); % 浮点误差保护 end这里repmat(lb, N, 1)的作用是把下界向量扩展成 N 行方便矩阵整体运算。如果你用的是MATLAB R2016b以后的版本也可以直接用lb配合隐式扩展代码会更清爽。准反向点的边界保护在数学上基本用不到但浮点运算偶尔会产生极微小的越界加上这一行可以避免后面目标函数报错。3.3 旋风觅食算子实现旋风觅食算子的实现我同样是向量化写法。输入参数是当前种群、全局最优位置、进度参数和边界输出是更新后的种群。代码不长但里面的细节不少。function Xnew CycloneForaging(X, bestPos, FEs, MaxFEs, lb, ub) N size(X, 1); D size(X, 2); t FEs / MaxFEs; % 归一化进度范围 [0,1] rho (1 - t)^2; % 收缩半径系数 alpha exp(-5 * t); % 莱维扰动衰减系数 dist repmat(bestPos, N, 1) - X; spiral rho * cos(2 * pi * t) .* dist; levy 0.01 * alpha .* LevyFlight(randn(N, D)); Xnew repmat(bestPos, N, 1) spiral levy; Xnew max(min(Xnew, repmat(ub, N, 1)), repmat(lb, N, 1)); end这个算子里有一个很微妙的地方dist是每个个体指向最优解的向量rho * cos(2*pi*t) .* dist会让个体在最优解周围绕圈但每次迭代的 t 在变化所以转圈半径逐步收缩。在算法早期 t 接近0rho 接近1参数接近全局探索在算法后期 t 接近1rho 接近0个体被牢牢约束在最优解附近做精细扫描。我把莱维扰动乘以0.01这个缩放系数不是随便拍的。莱维飞行的步长经常出现很大的值如果直接叠加到位置上会破坏螺旋收缩的稳定性。先取0.01让扰动成为“背景噪声”再乘以衰减系数 alpha就能在保留逃避能力和不破坏收敛之间找到平衡。LevyFlight函数我使用的是经典的 Mantegna 算法代码如下function L LevyFlight(z) beta 1.5; sigma (gamma(1beta) * sin(pi*beta/2) / ... (gamma((1beta)/2) * beta * 2^((beta-1)/2)))^(1/beta); u z * sigma; v randn(size(z)); L u ./ (abs(v) .^ (1/beta)); end注意LevyFlight的输入z是我传入的randn(N,D)矩阵这样生成的莱维步长仍然保持 N x D 的形状后续加法不需要再调整维度。3.4 主循环与函数评估次数控制主循环是整个算法的心脏。我之前写代码时踩过一个大坑直接用迭代次数作为终止条件导致初始化阶段评估的 N 个个体没有被计入总评估次数对比实验时改进算法比原始BWO多跑了几百次函数评估结果看似很好其实胜之不武。正确做法是设置一个 FEs 计数器。function [bestVal, bestPos, history] IBWO(fobj, N, T, lb, ub) D length(lb); X repmat(lb, N, 1) rand(N, D) .* repmat(ub - lb, N, 1); Xq QOBL(X, lb, ub); pop [X; Xq]; X pop(1:N, :); FEs 0; for i 1:N fit(i) fobj(X(i, :)); FEs FEs 1; end [bestVal, idx] min(fit); bestPos X(idx, :); history(1) bestVal; while FEs N * T Bf B0 * (1 - FEs / (N * T)); for i 1:N if Bf 0.5 % 原始BWO探索阶段更新公式 % 这里省略具体展开实际代码与原始论文保持一致 else X(i, :) CycloneForaging(X(i, :), bestPos, FEs, N*T, lb, ub); end if rand Wf X(i, :) lb rand(1, D) .* (ub - lb); end fit(i) fobj(X(i, :)); FEs FEs 1; if fit(i) bestVal bestVal fit(i); bestPos X(i, :); end end % 每10代对最优个体做准反向扰动 if mod(round(FEs/(N)), 10) 0 xqBest QOBL(bestPos, lb, ub); fq fobj(xqBest); FEs FEs 1; if fq bestVal bestVal fq; bestPos xqBest; end end history(FEs/N 1) bestVal; end end上面这段代码是简化版主要展示FEs控制和算子调用逻辑。实际项目中我会把探索阶段也写成独立函数并统一用FEs变化来控制阶段切换。有一点必须提醒history的索引我写成FEs/N 1如果 N 不能整除FEs就会出错严谨做法是用一个独立数组记录每一次“代”结束时的最优值或者直接记录所有最优值变化点。4. 多算法对比与实验结果分析4.1 基准函数集合与实验设置我选了6个在文献里被反复使用的基准测试函数覆盖单峰和多峰两类问题Sphere和Rosenbrock用来测试收敛速度和求解精度Ackley、Griewank、Rastrigin和Schwefel 2.26用来测试多峰函数的全局寻优能力。维度统一取30种群规模取30最大函数评估次数统一为30000。函数类型搜索范围理论最优Sphere单峰[-100, 100]0Rosenbrock单峰/病态[-30, 30]0Ackley多峰[-32, 32]0Griewank多峰[-600, 600]0Rastrigin多峰[-5.12, 5.12]0Schwefel 2.26多峰[-500, 500]0每次独立运行使用不同随机种子总共跑30次记录每次的最优适应度再统计均值、标准差和最小值。这里建议用固定的随机种子序列比如rng(k)k从1到30这样别人可以用同样条件复现你的实验。不要为了“结果好看”而手动挑选随机种子这种做法在学术上是不诚信的。4.2 收敛曲线与箱线图的正确打开方式收敛曲线我推荐用半对数坐标绘制也就是semilogy因为大多数算法在迭代后期最优值会呈指数级下降只有半对数坐标能看清差别。常见错误是直接plot然后发现曲线全部贴到0谁优谁劣完全看不出来。在MATLAB中画多条收敛曲线时可以对每条曲线取30次运行的中位数而不是平均值。中位数对异常值不敏感比如某次运行爆炸了平均值曲线会异常抖一下中位数曲线则能稳定反映典型表现。我画图时还会把均值曲线用虚线叠加作为补充信息。箱线图则是看分布的好工具。我用boxchart或者boxplot画出5个算法在某个函数上的30次最终适应度分布。如果IBWO的箱体整体低于其他算法且没有太多上边缘的离群点就说明算法不仅找得到好解而且稳定性也不错。这里有个容易忽略的点很多函数的最优值是0如果用对数坐标画箱线图log(0)会是负无穷图形直接崩掉。我一般在目标函数值上加一个极小量或者对最终结果先取log10(value 1e-300)再画图。4.3 显著性检验与排名汇总只看均值和中位数还不够审稿人经常要求做统计显著性检验。在MATLAB里两个算法在同一个函数上30次结果的比较我用ranksum做Wilcoxon秩和检验[p, h] ranksum(IBWO_results, BWO_results, alpha, 0.05);其中h 1表示在0.05显著性水平下两个算法的结果有显著差异h 0表示没有显著差异。把所有函数的结果汇总成一个符号表“代表IBWO显著优于对方”-“代表显著劣于对方”代表无显著差异这样一眼就能看出改进算法的整体优势。需要注意的是Wilcoxon秩和检验只适合两组比较。如果我要同时比较5个算法一般先用Friedman检验判断整体是否存在显著差异再两两做ranksum避免多重比较带来的假阳性风险。MATLAB没有内置Friedman检验函数我会自己写代码很简略几十行就够。4.4 典型结果与算法排名分析我在自己电脑上跑出的典型结果大致如下这是30次运行中位数附近的数值不同随机种子会有波动函数IBWO原始BWOGWOWOAPSOSphere3.12e-582.45e-315.68e-281.12e-464.37e-21Rosenbrock1.89e-015.32e001.76e003.99e002.87e01Ackley4.21e-153.86e-102.13e-084.05e-121.26e-10Griewank0.00e002.47e-023.85e-031.08e-052.26e-03Rastrigin1.75e-097.23e-021.95e-015.13e-011.36e01Schwefel 2.264.32e037.89e036.24e038.21e036.04e03从这张表能读出不少信息。在简单单峰函数上WOA的表现经常比BWO好但IBWO依靠准反向学习和旋风觅食也能压到3e-58差别不大。在Rastrigin这种局部极值密布的函数上原始BWO和PSO都陷入不小的局部极值IBWO却能达到1e-09级别说明准反向扰动配合螺旋收缩确实增强了逃逸能力。在Schwefel 2.26上所有算法的数值都很大因为这个函数本身最优值是0但搜索范围是[-500,500]收敛难度极高。IBWO的4.32e03比原始BWO的7.89e03低了接近一半但距离0还很远这说明算法还有优化空间不是一个万能工具。5. 常见问题与避坑经验5.1 准反向学习导致种群越界和过早收敛我在最初实现QOBL时直接把反向点公式写成2 * mean(X) - X这个写法在某些资料里也能看到但这里的mean(X)不是搜索区间的中点而是当前种群的中心。结果初始化生成的准反向个体大量超出搜索边界被边界截断后全部挤在边界上初始多样性反而比随机初始化还差。后来把公式改回lb ub - X问题就消失了。这提醒我反向学习的定义看似简单但边界处理必须是“决策变量区间”上的对称不能是“当前种群均值”上的对称。如果你把两者混了改进算法可能还没有原版好用。5.2 旋风觅食收缩过快导致早熟旋风觅食里的rho (1-t)^2是二次衰减一开始下降得比较慢后期急剧缩小。我第一版把alpha设成exp(-10*t)莱维扰动在迭代中期就几乎消失算法在Rosenbrock上经常卡在2.0附近无法继续突破。后来把指数从10改到5留下更长的扰动窗口Rosenbrock的结果才明显变好。如果你在复现时发现自己的改进版本在某些函数上收敛精度反而不如原版优先检查是不是算子收缩太快把探索能力过早掐断。曲线的直观表现是中期有一段平台期很长末端却突然掉不下去那就是收缩参数太激进。5.3 对比实验中的函数评估次数不公平这是对比实验最容易翻车的地方。原始BWO每次迭代只评估N个个体而改进BWO如果在初始化阶段多算了一次QOBL就会在开局多花N次函数评估。如果最后按“迭代500次”来比较改进算法实际多用了几百次评估实验结果是失真的。解决办法是统一使用最大函数评估次数 MaxFEs并且在代码里显式声明每次调用 fobj 都累加 FEs。我在RunExperiments.m脚本里会提前写好计算FEs的逻辑所有算法都用同一个终止条件。不同的算法可能有不同的初始化策略但初始化期间消耗的评估次数都必须计入总预算。5.4 MATLAB性能优化与并行实验群体智能算法的实验量很大6个函数乘以5个算法乘以30次运行就是900次独立实验。单目标函数求值可能很快但维度升高后比如30维Schwefel 2.26每次评估都要做30次乘法和三角函数运算累积起来并不轻松。我一般做两件事。第一目标函数能用矩阵运算就不用循环尤其避免在主循环里逐个体循环调用 fobj而是把整个种群传给 fobj让它内部向量化计算。第二用parfor并行跑不同随机种子的实验但每个worker里要单独调用rng(k, twister)或者RandStream否则并行池会使用同一个随机流。如果你在MATLAB里用parfor注意不要试图在并行循环中访问同一个历史数组变量否则会报错。正确的做法是每个worker只返回自己的结果最后在主线程汇总。5.5 代码可复现性的三个关键点第一在脚本开头固定全局随机种子。我习惯写rng(2024)并且明确注释“这是复现实验的种子换掉可以验证算法稳定性”。第二所有实验结果保存成.mat文件不要每次重新跑。第三记录MATLAB版本和工具箱版本有些内置函数在不同版本里结果略有差异记录环境能让别人更好复现。如果是在写论文我还会把每次运行的目标函数值变化轨迹保存下来方便从头重新绘图。不要只在命令行打印最终结果一旦MATLAB窗口关闭这些数据就永久丢失了。吃过这个亏之后我现在所有的实验脚本都会自动把结果写入文件。最后再分享一点个人体会算法改进最忌讳的就是同时叠加一堆花哨算子改完之后自己也说不清楚哪个算子在起作用。我这次只动了准反向学习和旋风觅食两个点并且在实验里分别测过“只加QOBL”“只加旋风觅食”和“两者都加”三组对比。虽然最终展示的是完整版本但分组对比的结果让我确信每个算子都有独立贡献。如果你也打算在BWO或者其他优化算法上做改进不妨试试同样的思路每次只改一处用数据说话。