简介蜣螂优化算法DBO是一种模拟蜣螂滚粪球行为的群智能全局优化方法非常适合科研人员、算法工程师以及竞赛选手在处理高维多模态函数最优化、模型参数辨识、机器学习超参数调优等复杂问题时使用。算法核心模拟蜣螂滚粪球时的寻路行为通过调整粪球滚动方向和距离不断探索解空间同时模拟个体间的社会交互并依据适应度评估进行选择淘汰代码完整呈现了种群初始化、滚动更新、位置调整、社会交互和迭代终止等关键环节便于读者理解局部搜索与全局搜索的平衡策略。资源以RAR压缩包形式发布整体大小约3.52MB包含可运行的DBO算法MATLAB源码以及配套的理论说明文档便于读者随时对照算法流程进行调试、改写与二次开发。已有1202人学习下载特别适合群智能优化算法入门者深入理解算法机制也适合工程实践者快速移植到自己的项目中。1. 蜣螂优化算法 DBO一份能直接跑通的 MATLAB 群智能优化实现我经常在标准测试函数上做对比实验粒子群在 Rastrigin 的 30 维搜索空间里提前收敛换用蜣螂优化算法DBO之后同样的边界和迭代次数平均最优值明显下降了一个台阶。这份 MATLAB 工程是 DBO 的完整落地实现从种群初始化、滚球、跳舞、繁殖、觅食到偷窃的位置更新全部收敛在一个主函数文件里入口一行命令目标函数换成自己的就能做参数寻优。适合两类人一类是刚接触群智能优化算法、想找个新对比算法的同学另一类是手上有 SVM 调参、图像阈值分割这类实际寻优问题、又不想从论文手搓公式的人。2. 算法原理拆解四种蜣螂行为如何映射到位置更新公式2.1 滚球与跳舞全局探索和转向重定向的数学表达蜣螂滚粪球时有两个状态没有障碍物就沿直线滚这时它要根据太阳或风向来导航遇到障碍物就爬上粪球“跳舞”重新判断方向再滚。放到优化问题里滚球就是常规的位置更新跳舞则是防止算法在一个方向上走死的重定向机制。先看滚球的核心更新公式MATLAB 工程里通常写成这样% 滚球行为朝全局最优方向滚动同时保留探索方向 k 0.1; % 偏转系数 b 0.3; % 环境变化系数 offset X(i, :) - Gbest; % 当前个体与全局最优的差值 if rand 0.5 % 一半概率靠向全局最优类似局部开发 X_new X(i, :) k .* offset - b .* abs(offset); else % 另一半概率背离最优保持探索能力 X_new X(i, :) - k .* offset b .* abs(offset); endk是偏转系数工程里通常取 0.1控制每一步相对自身位置的偏移幅度b是环境系数取 0.3模拟太阳光照、风向这些外部因素带来的随机扰动。Gbest是当前全局最优位置相当于蜣螂眼中的太阳方向。这里两个方向分支是很多复现版里默认的写法目的是让滚球群体不会全部挤向同一个点。跳舞行为的代码更短但坑更深% 跳舞遇到障碍物后重新定向theta 是随机角度 theta rand * pi; X_new X(i, :) tan(theta) .* abs(X(i, :) - X(max(1, i - 1), :));关键在tan(theta)。theta 在 0 到 pi 之间随机取正切值可能是正无穷也可能是负无穷所以一旦触发跳舞蜣螂可能瞬间跳到离当前位置非常远的点。这是 DBO 能跳出局部最优的核心原因。但正因如此跳舞触发概率不能设太高我习惯保持在 0.1 左右超过 0.3 之后算法行为就接近随机游走了。2.2 繁殖与觅食动态搜索区域如何随迭代收缩繁殖和觅食共享同一个核心设计动态边界收缩。雌性蜣螂产卵时会挑一个安全区域这个区域不固定而是围绕局部最优位置生成并且随着迭代次数增加逐渐缩小。% R 是动态收缩系数从 1 线性递减到 0 R 1 - t / MaxIter; % 产卵区域的动态上下界 Lb2 max(Lbest .* (1 - R), lb); Ub2 min(Lbest .* (1 R), ub); % 子代蜣螂的位置更新 b1 rand(1, dim); b2 rand(1, dim); X_new Lbest b1 .* (X(i, :) - Lb2) b2 .* (X(i, :) - Ub2);这里Lbest是局部最优位置工程实现里通常取排序后种群的中位数个体而不是全局最优。这一点很容易踩坑繁殖行为模拟的是“在一个还不错但并非最好的区域产卵”如果直接拿全局最优去生成产卵区域整个子群都会挤在同一个点附近探索能力就没了。R 1 - t / MaxIter是全局收缩节奏。迭代初期 R 接近 1产卵区域很大允许小蜣螂到处探索迭代后期 R 接近 0产卵区域收缩到Lbest附近精度就靠这个过程磨出来。觅食行为和繁殖对称区别在于它围绕全局最优构造搜索区域Lb3 max(Gbest .* (1 - R), lb); Ub3 min(Gbest .* (1 R), ub); C1 randn(1, dim); % 正态分布随机数 C2 rand(1, dim); % 均匀分布随机数 X_new X(i, :) C1 .* (X(i, :) - Lb3) C2 .* (X(i, :) - Ub3);注意C1用的是randn而不是rand。正态分布允许更新步长偶尔出现较大的负值相当于让觅食个体在动态区域内有概率反向探索这和滚球的跳舞行为形成互补。2.3 偷窃行为最差个体向全局最优靠拢的捷径最后一类行为是偷窃。小偷蜣螂自己不干活专门跟着最优个体去抢粪球。在算法里这部分个体是种群中适应度最差的群体它们的作用是保证极端差个体不会在边界外空转尽快回到有希望的区域。S 0.5; % 偷窃步长系数 g randn(1, dim); % 正态随机向量 X_new Gbest S * g .* (abs(X(i, :) - Lbest) abs(X(i, :) - Gbest));S控制偷窃步长g是正态随机向量乘上abs差值的含义是小偷会落在全局最优与局部最优之间的某个随机位置距离越远跳得越猛。这四类行为合起来构成一次完整的迭代滚球负责全局探索跳舞负责跳出局部最优繁殖在局部最优周边精细搜索觅食围绕全局最优扩大范围偷窃在两类最优之间补充样本。3. MATLAB 工程结构主函数与测试函数的代码逐段拆解3.1 文件清单与运行入口这套工程拿到手之后先别急着打开代码逐行读把文件结构摸清楚再说。最核心的只有两个文件入口脚本和算法主函数。DBO_Matlab/ ├── main_DBO.m % 入口脚本定义测试问题并画收敛曲线 ├── DBO.m % 算法主函数包含全部位置更新逻辑 ├── initialization.m % 种群初始化工具函数 ├── Sphere.m % 单峰测试函数 ├── Rastrigin.m % 多峰测试函数 └── Rosenbrock.m % 非凸测试函数入口脚本写得非常短核心就一段调用%% DBO 运行入口 clear; clc; close all; N 30; % 种群规模 MaxIter 500; % 最大迭代次数 lb -100; % 变量下界 ub 100; % 变量上界 dim 30; % 问题维度 [Best_pos, Best_fitness, Convergence_curve] ... DBO(N, MaxIter, lb, ub, dim, Sphere); disp([最优适应度: , num2str(Best_fitness)]); figure; semilogy(Convergence_curve, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优适应度); title(DBO 收敛曲线); grid on;DBO的六个入参顺序是固定的种群数、最大迭代、下界、上界、维度、目标函数句柄。前五个参数决定搜索空间和计算量第六个Sphere是你想优化的问题换成Rastrigin或Rosenbrock就能一键切换测试函数。返回的三个值里Best_pos是最优解的位置Best_fitness是最优适应度Convergence_curve是一条长度等于MaxIter的向量记录了每一代的最优适应度画图用。3.2 种群初始化与适应度排序初始化函数没有太多玄学但要注意 lb 和 ub 可能是标量也可能是向量。很多工程问题里不同维度的边界不一样比如 SVM 调参时惩罚系数 C 的范围是[0.01, 100]核参数 gamma 的范围是[0.001, 10]这时候 lb 和 ub 都得传向量。function X initialization(N, dim, ub, lb) % 种群初始化支持 lb/ub 传标量或行向量 if numel(lb) 1 % 标量边界直接向量化生成 X rand(N, dim) .* (ub - lb) lb; else % 向量边界每一维分别生成防止维度错位 X zeros(N, dim); for j 1:dim X(:, j) rand(N, 1) .* (ub(j) - lb(j)) lb(j); end end end初始化完成之后DBO 主函数做的第一件事是计算适应度并按升序排序。这一步是算法的关键设计后面循环里直接用索引区间来区分角色排序之后前 20% 是滚球蜣螂接下来 20% 是繁殖群体中间 40% 是觅食群体最后 20% 是偷窃群体。如果不排序角色分配就乱了。fitness zeros(N, 1); for i 1:N fitness(i) fobj(X(i, :)); end [fitness, index] sort(fitness); X X(index, :); Gbest X(1, :); % 全局最优位置 GbestScore fitness(1); % 全局最优适应度 Lbest X(ceil(N / 2), :); % 局部最优位置sort默认升序所以排序后X(1,:)是当前适应度最好的个体fitness(1)是对应的最优值。Lbest取中位数个体的位置而不是次优个体这是 DBO 与很多群智能算法不一样的地方它刻意避开了“所有参考点都指向同一个最好个体”的做法。3.3 主循环四种角色在同一轮迭代里的完整逻辑主循环是整个工程的心脏代码不长但每一行都有对应的生物行为。我贴一个精简且能直接跑的版本function [Best_pos, Best_fitness, Convergence_curve] ... DBO(N, MaxIter, lb, ub, dim, fobj) X initialization(N, dim, ub, lb); fitness zeros(N, 1); for i 1:N fitness(i) fobj(X(i, :)); end [fitness, index] sort(fitness); X X(index, :); Gbest X(1, :); GbestScore fitness(1); Lbest X(ceil(N / 2), :); Convergence_curve zeros(1, MaxIter); % 角色边界前 20% 滚球次 20% 繁殖中 40% 觅食末 20% 偷窃 r1 ceil(0.2 * N); r2 ceil(0.4 * N); r3 ceil(0.8 * N); for t 1:MaxIter R 1 - t / MaxIter; % 动态收缩系数从 1 降到 0 for i 1:N if i r1 % 滚球蜣螂 if rand 0.9 k 0.1; b 0.3; offset X(i, :) - Gbest; if rand 0.5 X_new X(i, :) k .* offset - b .* abs(offset); else X_new X(i, :) - k .* offset b .* abs(offset); end else % 跳舞重定向 theta rand * pi; X_new X(i, :) tan(theta) .* ... abs(X(i, :) - X(max(1, i - 1), :)); end elseif i r2 % 繁殖蜣螂 Lb2 max(Lbest .* (1 - R), lb); Ub2 min(Lbest .* (1 R), ub); b1 rand(1, dim); b2 rand(1, dim); X_new Lbest b1 .* (X(i, :) - Lb2) b2 .* (X(i, :) - Ub2); elseif i r3 % 觅食蜣螂 Lb3 max(Gbest .* (1 - R), lb); Ub3 min(Gbest .* (1 R), ub); C1 randn(1, dim); C2 rand(1, dim); X_new X(i, :) C1 .* (X(i, :) - Lb3) C2 .* (X(i, :) - Ub3); else % 偷窃蜣螂 S 0.5; g randn(1, dim); X_new Gbest S * g .* ... (abs(X(i, :) - Lbest) abs(X(i, :) - Gbest)); end % 边界约束 X_new max(X_new, lb); X_new min(X_new, ub); X(i, :) X_new; end % 重算适应度并排序 for i 1:N fitness(i) fobj(X(i, :)); end [fitness, index] sort(fitness); X X(index, :); Gbest X(1, :); GbestScore fitness(1); Lbest X(ceil(N / 2), :); Convergence_curve(t) GbestScore; end Best_pos Gbest; Best_fitness GbestScore; end这段代码里有几个容易忽略的工程细节。第一r1/r2/r3在循环开始前就算好不要在循环里用i 0.2 * N这种浮点比较整数边界判断在N不是整十的时候会出错。第二跳舞蹈里X(i - 1, :)在i 1时会越界所以写了max(1, i - 1)。第三边界约束max/min必须放在位置更新之后而且是在X_new上执行不是放在更新公式内部否则tan(theta)产生的极端值会直接污染下一次迭代。这段代码不做事后贪心筛选位置更新后全部直接替换。这是 DBO 原版设计的取舍四类行为本身已经承担了分工全量更新让种群在迭代早期保持足够的差异度。如果改成“只有新位置更优才接受”收敛确实会快不少但更容易陷入局部最优后面我会讲到这个改动。4. 运行与调参实战在标准测试函数上复现收敛曲线4.1 环境准备与三步跑通这套代码不依赖额外工具箱基础 MATLAB 环境就能跑。需要注意一点DBO主函数里用了randn、max、min、sort这些基础函数R2016a 以上的版本都能正常运行。如果你的机器上装了统计工具箱后面做显著性检验会更方便但不影响算法本身跑通。跑通的步骤就三步。第一步把main_DBO.m、DBO.m、initialization.m和三个测试函数文件放在同一个目录下。第二步打开main_DBO.m确认当前调用的目标函数是Sphere这个函数的边界是[-100, 100]30 维。第三步直接点运行。跑完你会看到两样东西命令行输出一个最优适应度数值以及一张迭代次数与适应度的收敛曲线图。以 Sphere 函数为例真实最优值在原点取到适应度理论上是 0。工程代码跑 500 代之后Best_fitness通常能落到1e-20量级这个精度对大多数工程问题完全够用。4.2 关键参数对收敛效果的影响DBO 需要手动调的核心参数不多但每个参数对搜索行为的影响差别很大。参数常见取值对搜索行为的影响种群规模 N30太小容易早熟太大计算量显著上升最大迭代 MaxIter500单峰函数 300 代足够多峰建议 1000 代偏转系数 k0.1越大探索半径越大但过大容易震荡环境系数 b0.3控制外部扰动的力度跳舞触发概率0.1越高越容易跳出局部最优但过高变随机游走偷窃系数 S0.5控制最差个体的跳跃幅度角色分配比例20/20/40/20影响全局探索与局部开发的平衡以 Rastrigin 函数为例它的边界是[-5.12, 5.12]30 维。如果沿用 Sphere 的[-100, 100]边界不会报错但搜索范围扩大了几十倍收敛精度会明显变差所以切换测试函数时边界一定要跟着改。% 切换成 Rastrigin 时入口参数改成这样 N 30; MaxIter 1000; lb -5.12; ub 5.12; dim 30; [Best_pos, Best_fitness, Convergence_curve] ... DBO(N, MaxIter, lb, ub, dim, Rastrigin);Rastrigin 是多峰函数局部极小值非常多直接跑 500 代经常会有几次实验落在局部最优附近。我的经验是迭代数调到 1000然后观察收敛曲线的形状如果曲线在中后期还在明显下降说明迭代数不够继续加如果曲线前 200 代就完全平了说明种群已经收敛再多迭代也白搭。4.3 多轮独立实验与算法对比设置单次运行的结果不具备统计意义群智能算法每次跑结果都不一样。严谨的做法是固定随机种子跑 30 次记录每次的最优适应度然后计算均值和标准差。我在做算法对比时用的模板是下面这样numRuns 30; results zeros(numRuns, 1); for r 1:numRuns rng(r); % 固定随机种子保证每轮实验可复现 [~, bestFitness, ~] DBO(N, MaxIter, lb, ub, dim, Rastrigin); results(r) bestFitness; end avgFit mean(results); stdFit std(results); fprintf(30 次平均最优值: %.6e标准差: %.6e\n, avgFit, stdFit);这里rng(r)是很多人容易漏掉的一步。固定随机种子之后同一台机器上每次跑这套脚本得到的结果完全一致这既是复现实验的基础也是排查算法 bug 的有力工具。如果你的结果每次差异很大第一件事就该检查是不是没固定随机种子。对比 PSO 或 GWO 时务必保持边界、维度、种群数、迭代数完全一致。另外每个算法独立跑 30 次不要用同一个随机种子序列这样对比结论才有说服力。5. 避坑与常见问题排查DBO 复现最容易翻车的五个点5.1 位置越界导致适应度全部变成 NaN现象跑着跑着Best_fitness变成 NaN收敛曲线中途断掉后面几代全是空值。原因跳舞公式里的tan(theta)在 theta 接近 pi/2 时会产出极大值一旦没做边界约束个体坐标可能飞出[lb, ub]几个数量级。如果目标函数里有log、sqrt或除法这些越界坐标会直接触发 NaN。解决强制在所有位置更新计算完成之后、写入X(i, :)之前加一行边界约束。最稳妥的写法是X_new max(X_new, lb); X_new min(X_new, ub);不要只约束部分分支滚球、繁殖、觅食、偷窃四条路径全部要过这一关。另外建议把 theta 限制在[0.05, pi - 0.05]区间从根本上避免 tan 爆炸。5.2 收敛曲线不下降或水平震荡现象跑了 500 代曲线几乎水平最优适应度比随机搜索还差或者曲线在后期来回跳动。原因最常见的是R的收缩方向写反了。正确写法是R 1 - t / MaxIter迭代初期 R 大、搜索区域宽后期 R 小、区域收缩。有人会顺手写成R t / MaxIter结果产卵和觅食区域越来越大算法后期还在乱飞曲线自然下不来。另一个常见原因是排序后没有更新Gbest整个循环一直拿初始最优位置当参考点。解决在循环里加一行disp([t, GbestScore])观察数值变化。如果GbestScore长时间不变检查Gbest的赋值语句是否在排序之后。我的习惯是每轮迭代结束后立刻打印前 10 代的收敛值确认曲线在下降再放长迭代跑完整实验。5.3 结果每次运行差异巨大无法复现现象同样的配置连续跑 5 次最优适应度差几十倍有时同一个参数组合这次表现好、下次表现差。原因没固定随机种子。MATLAB 默认每次启动rand和randn的状态都不同算法里大量使用随机数结果自然每次都不一样。这个不算 bug但会让参数调优变成玄学。解决在脚本最前面加rng(1)或者在多轮实验里用rng(r)其中 r 是循环变量。固定种子之后如果改动代码导致结果变化你能准确判断是改动引起的而不是随机波动。5.4 角色分配的整数边界算错群体行为混乱现象算法整体表现不稳定有时像纯随机搜索有时又像所有个体都在做同一种更新。原因角色边界用了浮点数比较。比如if i 0.2 * N当 N 不能整除时会因为浮点精度问题漏判或多判。更隐蔽的问题是边界重叠比如先写if i ceil(0.2 * N)再写elseif i ceil(0.4 * N)中间逻辑没问题但如果你写成i r1和i r2左开右闭区间就会漏掉第 r1 个个体。解决把所有角色边界提取到循环外用统一的整数变量r1 ceil(0.2 * N); r2 ceil(0.4 * N); r3 ceil(0.8 * N);然后在循环里严格按i r1、i r2、i r3、else四段划分。这样无论 N 是 30 还是 29四个区间都不会重叠或漏人。5.5 跳舞恒概率触发后期收敛精度被拖垮现象前期收敛速度不错但到迭代后期最优曲线还在上下波动最终精度明显不如同配置的 PSO。原因跳舞触发概率在整个迭代过程中是恒定的。后期种群已经聚集在最优区域附近此时再频繁触发tan(theta)这种大幅跳跃会把好不容易收敛的个体又炸出去。解决把跳舞概率改成随迭代衰减的值后期几乎不触发把搜索精力让给繁殖和觅食的局部收缩dance_prob 0.15 * (1 - t / MaxIter); if rand dance_prob % 跳舞重定向 else % 正常滚球 end这个改动我在 30 维 Rastrigin 上做过对比收敛精度提升非常明显而且不会破坏算法前期的探索能力属于性价比最高的一个改进。6. 二次开发把 DBO 换成自己的目标函数并验证效果6.1 替换目标函数的固定签名DBO 对目标函数的要求只有一个入参是 1×dim 的行向量返回一个标量适应度。这是这套工程扩展性最好的地方不管你是做 SVM 调参、图像阈值分割还是电力负荷预测只要把目标函数打包成这个签名就能接进来。以 SVM 参数寻优为例假设你有一份训练数据Xtrain和标签ytrain想用 DBO 同时优化惩罚系数 C 和高斯核参数 gamma目标函数可以写成function loss svm_loss(x, Xtrain, ytrain) C x(1); gamma x(2); model fitcsvm(Xtrain, ytrain, ... KernelFunction, rbf, ... BoxConstraint, C, ... KernelScale, 1 / sqrt(gamma)); loss kfoldLoss(crossval(model, KFold, 5)); end注意这里的参数空间只有二维维度改小之后收敛会快很多。在main_DBO.m里用匿名函数把训练数据传进去fobj (x) svm_loss(x, Xtrain, ytrain); % 参数范围C 和 gamma 的范围差异很大用向量边界 lb [0.01, 0.001]; ub [100, 10]; dim 2; [Best_pos, Best_fitness, ~] DBO(20, 100, lb, ub, dim, fobj);这里lb和ub传的是向量分别对应 C 和 gamma 的上下界。如果参数范围差几个数量级比如 C 到 100、gamma 到 0.001建议先取对数再寻优否则数值范围大的维度会主导搜索方向这是我调试时踩过的一个坑。6.2 加入自适应参数改进基础版 DBO 跑通之后二次开发最值得改的地方有两个跳舞概率衰减和精英保留。跳舞概率衰减的写法在上一章已经说过原理是前期多探索、后期少扰动。精英保留的做法更简单在每轮排序更新之前保存上一代最优如果更新后的最优个体反而变差就把旧最优放回第一位prevGbest Gbest; % ... 更新位置、重算适应度、排序 ... Gbest X(1, :); GbestScore fitness(1); % 精英保留历史最优永远不回退 if GbestScore fobj(prevGbest) X(1, :) prevGbest; Gbest prevGbest; GbestScore fobj(prevGbest); end加了精英保留之后收敛曲线就不会出现“回升”的现象。严格来说原始 DBO 的Gbest是当前代最优而不是历史最优理论上有回退的可能工程里加上这条保护是合理的。6.3 用显著性检验验证改进是否有效改了参数之后不能只凭一两次实验说“有效”我的验证习惯是跑 30 轮然后和原始版做秩和检验。MATLAB 里直接用内置函数[p, h] ranksum(dbo_baseline_results, dbo_improved_results); if h 1 disp(两组结果存在显著差异改进有效); else disp(差异不显著建议检查参数设置或增加迭代次数); endh 1说明改动前后有统计显著差异p值越小差异越明显。这套验证方法同样适用于 DBO 和 PSO、GWO 的对比实验。我把这套模板固定下来之后每次做算法实验都强制走一遍“固定随机种子、跑 30 轮、画收敛曲线、做显著性检验”的流程看起来麻烦但省下了很多重复劳动。希望帮到你。本文还有配套的精品资源点击获取