如果你最近在折腾群体智能优化算法,大概率已经把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 接近0,rho 接近1,参数接近全局探索;在算法后期 t 接近1,rho 接近0,个体被牢牢约束在最优解附近做精细扫描。
我把莱维扰动乘以0.01,这个缩放系数不是随便拍的。莱维飞行的步长经常出现很大的值,如果直接叠加到位置上,会破坏螺旋收缩的稳定性。先取0.01让扰动成为“背景噪声”,再乘以衰减系数 alpha,就能在保留逃避能力和不破坏收敛之间找到平衡。
LevyFlight函数我使用的是经典的 Mantegna 算法,代码如下:
function L = LevyFlight(z) beta = 1.5; sigma = (gamma(1+beta) * sin(pi*beta/2) / ... (gamma((1+beta)/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] | 0 |
| Rosenbrock | 单峰/病态 | [-30, 30] | 0 |
| Ackley | 多峰 | [-32, 32] | 0 |
| Griewank | 多峰 | [-600, 600] | 0 |
| Rastrigin | 多峰 | [-5.12, 5.12] | 0 |
| Schwefel 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 | 原始BWO | GWO | WOA | PSO |
|---|---|---|---|---|---|
| Sphere | 3.12e-58 | 2.45e-31 | 5.68e-28 | 1.12e-46 | 4.37e-21 |
| Rosenbrock | 1.89e-01 | 5.32e+00 | 1.76e+00 | 3.99e+00 | 2.87e+01 |
| Ackley | 4.21e-15 | 3.86e-10 | 2.13e-08 | 4.05e-12 | 1.26e-10 |
| Griewank | 0.00e+00 | 2.47e-02 | 3.85e-03 | 1.08e-05 | 2.26e-03 |
| Rastrigin | 1.75e-09 | 7.23e-02 | 1.95e-01 | 5.13e-01 | 1.36e+01 |
| Schwefel 2.26 | 4.32e+03 | 7.89e+03 | 6.24e+03 | 8.21e+03 | 6.04e+03 |
从这张表能读出不少信息。在简单单峰函数上,WOA的表现经常比BWO好,但IBWO依靠准反向学习和旋风觅食也能压到3e-58,差别不大。在Rastrigin这种局部极值密布的函数上,原始BWO和PSO都陷入不小的局部极值,IBWO却能达到1e-09级别,说明准反向扰动配合螺旋收缩确实增强了逃逸能力。
在Schwefel 2.26上,所有算法的数值都很大,因为这个函数本身最优值是0,但搜索范围是[-500,500],收敛难度极高。IBWO的4.32e+03比原始BWO的7.89e+03低了接近一半,但距离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或者其他优化算法上做改进,不妨试试同样的思路,每次只改一处,用数据说话。