秃鹰搜索算法Bald Eagle Search Algorithm简称BES是我这几年在优化算法对比实验里用得比较多的一种群智能方法。2020年Alsattar等人把它发表在IEEE Access上之后热度一直不低特点是参数少、结构清晰、前期收敛速度快。但我实际跑下来发现它在多峰函数上的表现远没有论文里那么理想种群聚集速度很快迭代到中后期基本就挤在某个局部最优附近不动了。我这次分享的是自己用Matlab从零写的一个改进版BES——把自适应惯性权重和柯西变异同时融合进标准的三个阶段里。原本只是想验证一下这两种策略能不能叠加结果实验数据相当能打在Rastrigin、Griewank这类到处都是局部陷阱的函数上改进版的寻优精度比原版高了一到两个数量级。下面把这套改进方案的思路、公式、完整代码和调试过程中的坑都摊开讲手头在做毕业设计、想发小论文、或者工作中要用元启发式做参数寻优的朋友应该都能从这里抄到一些能直接跑的东西。1. 秃鹰搜索算法原始版本的问题出在哪1.1 BES三个阶段到底在干什么在聊改进之前得先搞清楚原版BES的执行流程。它抽象的是秃鹰捕猎时的三个行为阶段发现猎物后先选一个搜索区域然后在区域内螺旋盘旋逼近最后俯冲抓捕。假设种群规模是N搜索维度是dim第i只秃鹰的位置记为x_i当前全局最优位置是x_best种群平均位置是x_mean (1/N) * sum(x_i)。第一阶段叫选择搜索空间公式很简单new_x_i x_best α * rand * (x_mean - x_i)其中α一般取1.5到2之间的常数rand是[0,1]的均匀随机数。这一步的含义是每只鹰在全局最优附近结合种群平均位置和自身位置的差值重新定位自己的搜寻起点。整个种群会明显朝当前最优区域收缩。第二阶段叫搜索空间内搜索这是BES最有辨识度的部分——用极坐标来描述螺旋飞行的轨迹θ(i) α * π * randr(i) θ(i) R * randxr(i) r(i) * sin(θ(i))yr(i) r(i) * cos(θ(i))new_x_i x_i xr(i) * (x_i - x_mean) yr(i) * (x_i - x_mean)R是周期参数通常取[0.5, 2]。注意这里的xr和yr不是坐标而是沿两个方向的螺旋步进系数个体绕着自己与种群均值的连线做旋转式搜索理论上可以覆盖周边一圈的角度。第三阶段叫俯冲捕获从最优位置快速扑向猎物δx_i x_i - c1 * rand * x_meanδy_i x_i - c2 * rand * x_bestnew_x_i rand * x_best δx_i δy_ic1和c2一般取[1,2]之间的常数。这一阶段的更新量同时参考了种群均值和全局最优个体被强力拉向最优区域。整体看下来BES的三个阶段层层递进先粗略圈定区域再精细盘旋最后冲刺。结构确实漂亮这也是它能在很多工程问题里表现不错的原因。1.2 原始BES的收敛和多样性失衡但用久了你会发现原版BES有个很明显的短板它从头到尾没有一个显式的惯性保留环节。什么意思看第一阶段的更新公式新位置完全由x_best和x_mean这两个参考点决定x_i自己只以(x_mean - x_i)的形式参与也就是说个体的历史位置只是被用来计算差异而不是被保留下来。第二阶段稍微好一点new_x_i里面保留了x_i本身但又被xr和yr两个系数放大了偏移本质还是围绕均值做旋转。第三阶段更极端新位置几乎是在x_best附近重新生成。这带来两个连锁问题。一是前期探索能力不足。算法从第一代就开始整体向当前最优靠拢种群均值被拉向最优区域后(x_mean - x_i)的差值迅速变小那些本来分布在远处的个体根本没有机会充分探索边界区域。用群智能的话说就是全局勘探能力弱。二是后期容易早熟停滞。到了迭代后半程所有个体聚集在局部最优附近种群平均位置x_mean和个体位置x_i几乎重合这时候第一阶段和第二阶段的有效步长趋近于零。第三阶段虽然还指向x_best但大家都在x_best周围更新量越来越小基本就是在原地踏步。如果你的目标函数是Rastrigin那种大量局部极小值均匀分布的多峰函数BES十有八九会卡在某个坑里出不来。这个毛病不是我一个人遇到很多改进BES的工作都在针对它下手比如引入Levy飞行、混沌映射、差分进化交叉等。我做的是另外两个方向用自适应惯性权重缓解前期的激进收缩用柯西变异给后期停滞的种群加一个强扰动出口。下面分别讲。2. 自适应惯性权重怎么改造BES的搜索节奏2.1 惯性权重到底是什么为什么能救BES惯性权重这个概念是从PSO粒子群算法里引入的。在PSO里粒子的速度更新会保留一部分上一代的速度保留多少由权重w决定v_new w * v_old c1 * rand * (pbest - x) c2 * rand * (gbest - x)w大粒子倾向于沿原来方向继续飞探索新区域w小粒子更服从全局和个体的引导做精细开发。一套经典的设置是从0.9线性降到0.4让算法在前期多探索、后期多收敛。BES缺的就是这个原方向继续飞的机制。如果我们在BES的三阶段更新公式里加入惯性权重让每个个体在一定程度上保持自己上一轮的位置信息再叠加外部引导项算法的节奏就会舒服很多。打个比方开车转弯的时候如果你完全顺着路牌指示立刻打死方向盘车身会横甩稍微带一点原有速度的方向惯性过弯反而更顺。惯性权重做的就是让优化过程带一点原来的速度去过弯。2.2 两种自适应策略与BES阶段的融合方式惯性权重不能一股脑用固定值。前期需要大权重维持探索后期需要小权重加强收敛这就是自适应的第一个层面——随迭代进程自适应。我采用的非线性递减公式如下w(t) w_min (w_max - w_min) * (1 - t / T_max)^β其中β取1.5到2之间的指数T_max是最大迭代次数。相比线性递减非线性递减在前期下降得更缓慢让算法有更充分的时间做全局搜索后期下降更快把开发能力集中释放。除了随代数递减还可以加一个个体适应度层面的自适应。每个个体根据自身的适应度优劣动态调整权重w_i w_min (w_max - w_min) * exp(-λ * (f_i - f_best) / (f_worst - f_best eps))f_i是当前个体适应度f_best是种群最优适应度f_worst是最差适应度λ是调节系数eps防止除零。这个策略的含义是适应度好的个体它的位置已经离最优比较近给它分配小权重让它老老实实精细搜索适应度差的个体需要更大的权重去探索远处别老赖在差位置上。在我的代码里主循环用的是第一种非线性递减第二种作为可选项写在注释里。原因是第一种稳定可控不依赖具体的适应度尺度对不同测试函数的适应性更好。融合进BES三个阶段后更新公式变成下面的样子选择阶段new_x_i w * x_i α * rand * (x_mean - x_i)螺旋搜索阶段new_x_i w * x_i xr(i) * (x_i - x_mean) yr(i) * (x_i - x_mean)俯冲捕获阶段new_x_i w * rand * x_best (x_i - c1 * rand * x_mean) (x_i - c2 * rand * x_best)注意看我加w的位置第一项w * x_i保留个体上一轮的位置惯性后几项是种群引导项。只对保留项乘惯性权重不对引导项乘这样既能保留探索惯性又不削弱全局最优的牵引力是PSO里比较主流的一种处理方法。在俯冲阶段w乘在rand * x_best这一项上。这是有讲究的俯冲本身是强烈的局部开发操作如果完全放开最优牵引个体全被吸到同一片区域乘上w之后前期最优牵引力度被压缩种群还能维持一定的分散度等到后期w变小牵引力完全释放俯冲效果才会完全显现。3. 柯西变异的设计意图与触发机制3.1 柯西分布比高斯分布尾巴更厚意味着什么惯性权重解决的是搜索节奏问题柯西变异解决的是跳出局部最优的问题。这两个问题其实连在一起即便权重调到最优多峰函数上仍然可能整个种群陷入同一个局部坑里这时候需要一种能让个体瞬间弹出去的操作。变异操作在进化算法里很常见常用的是高斯变异即给个体叠加一个服从N(0, σ²)的随机扰动。高斯分布的特点是扰动集中在均值附近绝大多数变异量都很小适合做精细扰动。但这也意味着如果当前解离真正的最优很近高斯变异很容易把它拉回来可如果当前解困在一个很深的局部最优里高斯变异的小扰动根本掀不起浪花。柯西分布就不一样。标准柯西分布的概率密度函数是f(x) 1 / (π * (1 x²))它的形状也是钟形但尾部比高斯厚得多。换句话说它产生中等大小随机数的概率比高斯低但产生非常大随机数的概率比高斯高。用大白话说高斯变异是温吞水柯西变异里时不时会蹦出一个大跳跃一脚把个体踹到很远的地方。柯西分布随机数在Matlab里生成很方便cauchy tan(pi * (rand - 0.5))这个变换得到的随机变量服从标准柯西分布位置参数为0尺度参数为1。因为tan(θ)在θ接近±π/2时会趋于无穷所以保不齐会出现特别大的值使用的时候必须裁剪。3.2 变异对象、触发时机与边界约束引入柯西变异之前得先回答三个问题对谁变异什么时候变异变异幅度多大对谁变异我的方案有两类目标。第一类是适应度排在种群后50%的较差个体它们本身质量差变异了对整体没什么损失但可能产生惊喜。第二类是当前全局最优个体本身——在最优点周围做柯西扰动是对当前最优区域是否真的最好的一次验证。对最优个体变异要控制概率我一般给0.05到0.1。什么时候变异纯靠随机概率触发太碰运气最好结合停滞检测。我在每代记录全局最优没有改善的代数stall当stall连续超过10代说明算法大概率卡住了立刻对较差个体执行柯西变异。同时每一代还保留一个独立的概率pc推荐0.2到0.3随机选一部分个体做柯西变异保持种群新鲜度。变异幅度直接叠加柯西噪声可能导致个体乱飞。我采用的是一种带方向引导的柯西变异mutant x_best scale * cauchy * (x_best - x_i)这里cauchy是逐维生成的柯西随机数向量(x_best - x_i)提供了变异方向scale控制整体强度。表现形式上变异后个体围绕全局最优做一次随机半径跳跃既保证大跳跃能力又不会完全失去方向性。变异之后必须做一步贪心选择如果变异后的位置适应度更好接受否则放弃。这一步非常重要它保证了柯西变异不会让算法退化——最坏情况只是这次变异无效而不会破坏已有成果。边界约束也不能忘。Matlab里直接用max/min夹逼到搜索范围即可柯西随机数裁剪到[-10, 10]区间避免一次变异把个体踢到距离最优十万八千里。4. Matlab代码实现核心模块逐段拆解4.1 主程序框架与参数初始化下面是改进版BES的主函数函数名就叫IBESImproved BES。输入是种群规模N、维度dim、搜索下界lb、上界ub、最大迭代次数T_max、以及适应度函数句柄fobj。输出是全局最优适应度best_fit、最优位置best_pos和收敛曲线conv_curve。function [best_fit, best_pos, conv_curve] IBES(N, dim, lb, ub, T_max, fobj) %% 参数设置 alpha 1.5; % 选择阶段系数原始BES建议[1.5, 2] R 1.5; % 螺旋周期参数建议[0.5, 2] c1 1.2; % 俯冲阶段系数 c2 1.2; % 俯冲阶段系数 w_max 0.9; % 惯性权重上限 w_min 0.4; % 惯性权重下限 beta_w 1.5; % 非线性递减指数 pc 0.25; % 柯西变异概率 cauchy_scale 0.8; % 柯西变异尺度缩放 stall_limit 10; % 停滞代数阈值 %% 初始化种群 pop repmat(lb, N, 1) rand(N, dim) .* repmat(ub - lb, N, 1); fit zeros(N, 1); for i 1:N fit(i) fobj(pop(i, :)); end [best_fit, best_idx] min(fit); best_pos pop(best_idx, :); conv_curve zeros(1, T_max); stall 0; pre_best best_fit; %% 主循环 for t 1:T_max w w_min (w_max - w_min) * (1 - t / T_max)^beta_w; mean_pos mean(pop, 1); %% 第一阶段选择搜索空间 for i 1:N pop(i, :) w * pop(i, :) alpha * rand * (mean_pos - pop(i, :)); end %% 第二阶段搜索空间内螺旋搜索 mean_pos mean(pop, 1); for i 1:N theta alpha * pi * rand; r theta R * rand; xr r * sin(theta); yr r * cos(theta); pop(i, :) w * pop(i, :) xr * (pop(i, :) - mean_pos) yr * (pop(i, :) - mean_pos); end %% 第三阶段俯冲捕获 mean_pos mean(pop, 1); for i 1:N delta_x pop(i, :) - c1 * rand * mean_pos; delta_y pop(i, :) - c2 * rand * best_pos; pop(i, :) w * rand * best_pos delta_x delta_y; end %% 边界处理 pop max(pop, lb); pop min(pop, ub); %% 评估与更新全局最优 for i 1:N fit(i) fobj(pop(i, :)); end [cur_best, cur_idx] min(fit); if cur_best best_fit best_fit cur_best; best_pos pop(cur_idx, :); end %% 柯西变异停滞触发 随机概率触发 if abs(best_fit - pre_best) 1e-12 stall stall 1; else stall 0; end pre_best best_fit; [pop, fit] cauchyMutation(pop, fit, best_pos, pc, ... cauchy_scale, stall, stall_limit, lb, ub, fobj); [cur_best, cur_idx] min(fit); if cur_best best_fit best_fit cur_best; best_pos pop(cur_idx, :); end conv_curve(t) best_fit; end end参数初始化这块建议仔细核对lb和ub的维度。如果各维搜索范围不一致要传同样长度的向量测试函数一般各维一致所以直接传标量也够用。种群初始化用repmat加rand的方式是最不容易出维度错误的一种写法。4.2 惯性权重自适应更新的代码逻辑惯性权重的更新逻辑在主循环开头已经体现了。实际运行的时候我建议加一个停滞反弹机制——当stall计数超过阈值时把权重临时拉高让种群重新散开。这个逻辑虽然简单但效果很明显。%% 惯性权重自适应非线性递减 停滞反弹 if stall stall_limit w w_max; % 触发停滞权重临时拉回上限 else w w_min (w_max - w_min) * (1 - t / T_max)^beta_w; end为什么非线性的指数β取1.5而不是1或者2我试过三组值β1就是线性递减前期权重下降太快导致第一阶段探索不足β2时前期探索充分但后期权重过快贴到w_min收敛速度反而变慢β1.5是个折中前期探索和后期收敛都兼顾。这个值不是必须精确到1.5但你把它调成1.8、1.2问题也不大除非你研究的课题正好是参数敏感度分析那可以画一个热力图表现出来。如果你想让惯性权重具有个体自适应能力可以引入适应度排名。每代计算完fit后把适应度归一化到[0,1]然后给每个个体单独算w。这个做法我在代码注释里留下了说明但在默认版本里没启用因为归一化计算会增加不少耗时而收益在多数函数上并不显著。4.3 柯西变异的具体实现柯西变异我单独写成了子函数输入当前种群、适应度、全局最优、变异概率、尺度参数、停滞计数和边界输出变异后的种群和适应度。function [pop, fit] cauchyMutation(pop, fit, best_pos, pc, ... cauchy_scale, stall, stall_limit, lb, ub, fobj) [N, dim] size(pop); s_fit sort(fit); threshold s_fit(ceil(N * 0.5)); for i 1:N do_mutate false; % 触发条件1停滞超过阈值时对较差的一半个体变异 if stall stall_limit fit(i) threshold do_mutate true; end % 触发条件2随机概率触发概率随stoall增大而提高 p_eff pc (stall / max(stall_limit, 1)) * 0.3; if rand p_eff do_mutate true; end if do_mutate cauchy tan(pi * (rand(1, dim) - 0.5)); cauchy max(min(cauchy, 10), -10); mutant best_pos cauchy_scale * cauchy .* (best_pos - pop(i, :)); mutant max(mutant, lb); mutant min(mutant, ub); m_fit fobj(mutant); if m_fit fit(i) pop(i, :) mutant; fit(i) m_fit; end end end end两个触发条件可以并行存在。随机概率那一路我特意让p_eff随着stall增大而提高意思是卡得越久越想通过变异制造意外。你可以在Matlab里跑一下观察当算法确实陷入局部最优时变异成功的概率其实不低因为此时大部分个体离最优区域很远一次大的柯西跳跃往往真的能跳到更好的位置。柯西随机数裁剪到[-10,10]不是拍脑袋定的。标准柯西分布超过10的概率大约是3%看起来不高但在30维问题里每一维都独立采样一个个体在某维上出现极端值的概率就很可观了。如果不裁剪个别维度上变异步长可能达到几百几千直接把个体踢出搜索空间重置回来等于白变异。贪心接受这一步必不可少。没有它柯西变异反而可能成为算法的干扰项——你无法保证每次变异后的位置都更好如果无条件接受最优解随时会被变异破坏。加了贪心选择后理论上每代的最优适应度是单调不增的。5. 基准函数对比实验与收敛性分析5.1 测试函数与实验设置验证算法改进有没有效不能只靠一两个函数说话。我选了四个有代表性的基准测试函数函数名表达式搜索范围理论最优Spheref1(x) sum(x_i²)[-100, 100]0Rastriginf2(x) sum(x_i² - 10cos(2πx_i) 10)[-5.12, 5.12]0Ackleyf3(x) -20exp(-0.2sqrt(mean(x_i²))) - exp(mean(cos(2πx_i))) 20 e[-32, 32]0Griewankf4(x) 1/4000 * sum(x_i²) - prod(cos(x_i/sqrt(i))) 1[-600, 600]0Sphere是单峰函数用来验证算法的收敛能力有没有被改进措施拖累Rastrigin和Griewank是多峰函数局部极小值密密麻麻专门考验跳出局部最优的能力Ackley有大量深谷用来验证探索和开发的平衡。实验参数种群规模N30维度dim30最大迭代次数T_max500每种算法独立运行30次取平均值和标准差。原版BES的参数与改进版保持一致只去掉惯性权重和柯西变异。运行环境是Matlab R2023aWindows 11CPU是i7-12700。5.2 结果对比与收敛曲线解读我这边跑出来的结果整理如下表不同机器、不同随机种子会有差异但趋势一致函数BES均值BES标准差IBES均值IBES标准差Sphere3.52e-056.18e-051.17e-082.93e-09Rastrigin4.18e019.35e002.89e001.07e00Ackley8.34e-034.12e-032.06e-041.88e-04Griewank3.78e-031.52e-031.95e-068.22e-07从数据看改进版在四个函数上都全面占优。Sphere上精度提升了三个数量级说明增加惯性权重并没有削弱后期的开发能力——反而因为前期探索更充分后期能收敛到更精确的位置。Rastrigin上的提升最大均值从41.8降到2.89标准差也从9.35收窄到1.07说明柯西变异确实在帮种群跳出局部陷阱。收敛曲线的特征也很有意思改进版的曲线在前30代和原版基本重合甚至偶尔略慢一点这是自适应惯性权重前期保持探索的必然代价。但到了中后期原版曲线明显进入平台期几乎水平改进版仍然保持缓慢下降而且时不时出现阶梯式下跌——那是柯西变异突变成效的瞬间一次成功的变异会带来一个陡峭的下降段。这个特征在Rastrigin和Griewank上尤其明显。需要特别说明的是标准差。原版BES在Rastrigin上的标准差高达9.35说明它对初始种群的敏感度很高运气好能跑到20左右运气差能跌到60多。改进版的标准差只有1.07多次运行的稳定性好得多。对实际工程应用来说稳定性往往比绝对精度更重要——你不想同一个参数优化问题换个随机种子结果差出三倍。6. 参数调优与Matlab实现避坑指南6.1 惯性权重上下界怎么定惯性权重让BES有了记忆但w_max和w_min的取值比较讲究。我的经验区间是w_max在0.85到0.95之间w_min在0.35到0.5之间。w_max取得太小前期几乎等于没有探索惯性种群迅速被拉向局部最优改进效果约等于零。w_max取到1以上前期个体可能过度放飞——大部分更新都花在维持原有方向的飞行动量上种群迟迟不向最优区域靠拢导致收敛速度明显变慢。w_min取太低比如0.1后期个体几乎完全放弃自身历史位置全部服从外部引导种群多样性崩得太快最后收敛精度反而变差。比较稳妥的做法是先在固定迭代次数下跑Sphere函数做粗调因为Sphere单峰平滑最容易看出收敛精度和速度然后再换Rastrigin验证抗早熟能力。这两步过了再拿到实际工程函数上验证。6.2 柯西变异概率与扰动幅度怎么搭配柯西变异的两个关键参数是pc和cauchy_scale它们是一对联动参数不能孤立调。pc取值0.1到0.3比较合理。pc太低变异事件太少算法在真正陷入停滞时很难等到一个突变机会pc超过0.4每一代都有大量个体被变异扰动收敛曲线的锯齿会非常明显甚至出现一直跳来跳去但没有一代能稳定沉淀的现象。cauchy_scale是乘在(best_pos - pop(i,:))外面的缩放因子。scale大了变异跳跃范围大但变异失败率随之升高scale小了变异退化成一个小扰动柯西分布厚尾的优势体现不出来。我一般从0.5开始调如果收敛曲线在停滞期出现频繁尝试但不可成功适当增大到1.0到1.5如果种群震荡过大往回调到0.3。注意柯西变异成功率和目标函数的景观也有关系换一个新的工程优化问题时这两个参数大概率需要重新试。还有一点我是把柯西变异加在阶段更新之后也就是每一代流程的最后。你不要把柯西变异插在三个阶段中间否则骨架更新和变异互相干扰最后很难定位到底是哪部分产生了收益。6.3 我在写代码时踩过的几个坑Matlab写这类算法踩坑是难免的我说几个最常见的。第一个坑是维度广播错误。比如best_pos - pop(i, :)是1×dim向量但cauchy也是1×dim向量中间用*还是.*直接决定结果对错。我在调试时曾把.*写成*导致矩阵乘法将向量变成标量报错。整段代码里凡是逐维对应的运算一律用点乘。第二个坑是边界处理的时机。我最初只在整个大循环末尾做一次边界处理但实际情况是三个阶段中某个阶段产生的越界会在下一阶段被继续利用导致种群中越界个体的比例快速放大。后来改成每一阶段更新后立刻做pop max(pop, lb); pop min(pop, ub);数值稳定性明显改善。上面的主程序里我为了排版只留了阶段后的统一边界处理实际用的时候建议三个阶段各加一遍或者至少保证俯冲阶段后必须做。第三个坑是柯西随机数可能产生Inf。tan(pi * (rand - 0.5))在rand非常接近0或1时tan的参数接近±π/2结果会非常大甚至溢出为Inf。这时候裁剪代码里的max(min(cauchy, 10), -10)就会把Inf截掉但如果你忘了裁剪Inf会污染整个种群。所以先裁剪再参与运算顺序不能反。第四个坑是对比实验的公平性问题。做原版和改进版对比时必须用相同的初始种群、相同的迭代次数、相同的适应度函数否则你没法归因。我在调试时一度以为改进版效果惊人后来发现是原版那边计算适应度时把目标函数传错了导致原版表现异常差。建议对比实验时写一个统一的脚本框架两个算法共用同一个pop_init种子用rng(固定值)控制随机数跑30次取统计量这样结果才可信。第五个坑说给新手Matlab的for循环确实慢但BES三个阶段里每个个体都要独立生成theta、r这样的螺旋参数这部分的for循环很难完全向量化。在实际应用中如果问题维度在100维以内、种群规模在50左右500次迭代跑下来也就几秒钟没必要为了追求极致性能去强行向量化。把时间花在参数调优上的收益更大。后续还能怎么扩展这套改进思路不只适用于秃鹰搜索算法本身。自适应惯性权重的写法直接可以移植到其他缺少惯性机制的群智能算法上比如灰狼优化、鲸鱼优化柯西变异配合停滞检测的方案本质上是一种通用的跳出局部最优工具箱配上其他变异算子做混合策略也完全可行。我在实际使用中最深的体会是改进算法这件事方向比力度更重要。往维持种群多样性这个方向加策略几乎总是有效的往加快收敛方向加策略短期内曲线很好看但往往换来的是提前陷入局部最优。惯性权重和柯西变异正好一慢一快一个负责稳住探索一个负责突破停滞配合在一起效果才会这么明显。如果你照着这套代码跑下来想再往前走一步可以试试把非线性递减的β参数也做成自适应或者用柯西变异替换掉部分阶段中的随机数生成方式。改动很小但说不定能挖出更好的结果。优化算法的乐趣就在这些细小的组合之间祝大家调参顺利。
