基于蜣螂优化算法的无线传感网络覆盖优化MATLAB实现

1. 项目动机与方案选型思路去年我在做无线传感节点布点方案时先试了遗传算法后面又换粒子群效果都不太理想。节点数量一多坐标组合呈指数级增长传统算法经常陷入局部最优后期收敛速度还拖沓。后来翻到2022年提出的蜣螂优化算法DBO抱着试一试的心态把它用到WNS覆盖优化场景里结果同参数条件下覆盖率比PSO提升了近十个百分点收敛速度也快了一倍。这篇博客就把完整可运行的MATLAB实现和代码注释整理出来同时把建模思路、算法原理、参数调优心得和踩坑记录一并写清楚希望对做无线传感网络覆盖、群智能优化算法对比或者正在做毕业设计的同学有帮助。先解释一下项目标题里的几个词方便不同基础的读者对齐信息。WNS指的是无线传感节点系统Wireless Node Sensing/System本质上是无线传感器网络WSN的节点感知场景。覆盖优化算法解决的核心问题是在给定的监测区域内如何布置或者重新调整传感器节点的位置让整个区域被传感信号覆盖的面积尽可能大。DBO则是Dung Beetle Optimizer蜣螂优化算法是近年涌现的群智能优化算法里结构清晰、收敛性能比较突出的一种。三者放在一起项目任务就清晰了把传感器节点的坐标看作决策变量用DBO算法去搜索一组最优坐标使得栅格化区域的覆盖率最大化并在MATLAB环境里实现、验证和输出覆盖效果图。1.1 覆盖优化问题的本质很多人第一次接触覆盖优化会觉得这就是选点位置听起来不复杂但真正建模后会发现问题远比想象中棘手。第一步监测区域是一块二维平面节点感知范围通常抽象成以节点为圆心、感知半径为R的圆形区域。布置N个节点每个节点有x和y两个坐标所以决策变量的维度是2N。N30时问题就变成了一个60维的连续优化问题搜索空间的大小完全取决于区域尺寸盲目枚举不现实。第二步覆盖率怎么算。现实工程里没有解析解大家普遍使用栅格法把监测区域按一定步长离散成网格点统计被至少一个节点覆盖的格点占全部格点的比例。网格步长越小计算结果越接近真实覆盖面积但同时计算量也越大。这一步是整个项目最大的性能瓶颈后面我会给出具体的向量化实现方式以及粗栅格快速验证、细栅格精调的工程技巧。1.2 为什么选DBO而不是PSO或GA做一个算法选择的决策不能只看谁名气大。粒子群算法PSO实现简单但容易出现早熟尤其在覆盖优化这种高维、多峰的适应度地形里粒子一旦被局部极值吸引很难跳出来。遗传算法GA的交叉变异提供了多样性但它参数多收敛速度慢而且对连续坐标问题的编码和解码环节会引入额外的计算开销。DBO的优势在于它同时模拟了蜣螂的五种行为滚球、跳舞、产卵、觅食和偷窃这就构成了全局探索局部开发边界扰动三层互补的搜索机制。用大白话说滚球蜣螂负责在整个区域内推进搜索跳舞行为负责在当前路径附近变向逃避局部陷阱产卵行为围绕局部最优建巢觅食行为围绕全局最优精细搜索偷窃行为则在全局最优附近小范围扰动。这种多角色协作结构让DBO在同等迭代次数下能获得更好的解。另外从复现角度讲DBO的位置更新公式并不复杂不涉及复杂的编码或算子设计用MATLAB天然向量化地实现核心代码不到一百行非常适合作为对比算法写进论文里。1.3 整体技术流程设计整个项目的技术链路可以拆成四大块场景建模、覆盖率计算、DBO搜索、结果可视化。场景建模部分需要确定区域尺寸、节点数量、感知半径等基础参数覆盖率计算封装成独立的函数核心是栅格点与所有节点的距离判断DBO搜索部分实现五种行为的位置更新并按未覆盖率最小化的目标迭代优化结果可视化部分把最终节点位置、覆盖区域和收敛曲线画出来。下面各章节就按这个链路逐个展开代码会尽量保留完整、加上逐行注释让你能直接照着跑通。2. 覆盖优化问题的数学建模2.1 问题假设与场景定义为了让问题可解需要对场景做几个标准假设这也是WNS/WSN覆盖研究领域的通用做法。第一监测区域是矩形平面长和宽分别记为L和W。节点坐标限制在区域内不能越界。第二每个传感器的感知半径相同且固定为R。感知模型采用最简单的布尔圆盘模型即栅格点若到某个节点的欧氏距离小于等于R就认为该点被完全覆盖不存在概率衰减或随机误差。第三节点之间不要求互相通信只关心对监测区域的感知覆盖。第四覆盖率为纯几何覆盖率即被覆盖网格点数占总网格点数的比例不考虑边界效应修正。在这些假设下优化目标非常明确寻找一组节点坐标X [(x1,y1), (x2,y2), ..., (xN,yN)]使得区域内覆盖率最大。因为DBO标准版本是按照最小化问题设计的所以适应度函数通常写为fitness 1 - coverage优化方向就是让未覆盖率不断下降。2.2 栅格法计算覆盖率栅格法是最直观也最稳妥的覆盖率评估方式我在实际项目里一直用它。原理可以用一句话概括用离散点采样替代连续面积计算。假设监测区域为100m × 100m栅格步长取1m那么x方向有101个采样点y方向有101个采样点总计10201个栅格点。对于每个栅格点计算它到所有N个节点的距离取最小值。如果最小距离小于感知半径R则该栅格点被覆盖标记为true。遍历完所有栅格点后用覆盖栅格点数量除以总栅格点数量就得到覆盖率。代码实现上有两种写法我推荐用meshgrid生成全部格点矩阵然后对每个节点做一次向量化距离计算避免三层for循环的缓慢执行。这一块的MATLAB函数我放在第四章那里会给出完整的代码。2.3 目标函数与维度映射优化问题的变量维度是2N前N维是节点的x坐标后N维是节点的y坐标。这样的布局在MATLAB里非常好处理因为reshape一行就能还原成N行2列的坐标矩阵。例如一个种群个体position向量的第1到N个分量对应N个节点的x坐标第N1到2N个分量对应y坐标那么节点坐标矩阵就是reshape(position, N, 2)。适应度计算时先把向量拆成坐标矩阵再调用覆盖率函数最后返回未覆盖率。这里有一个关键细节覆盖率函数内部需要遍历N个节点每个节点都要对全部栅格点做一次欧氏距离计算并做或运算。如果节点数量是30、栅格点数量上万那么一次适应度评估就是30次万级距离计算一个种群30个个体就是90万次。所以覆盖率函数的执行效率直接决定整个算法的运行时间必须认真优化避免低效循环。3. 蜣螂优化算法核心机制3.1 灵感来源与角色分工蜣螂优化算法模拟的是蜣螂在自然界中滚粪球、产卵、觅食和偷窃等行为。很多人第一眼看到偷窃会觉得这算法很新颖其实就是把生物界竞争机制引入了优化过程让群体不仅围绕自身位置搜索还能在全局最优附近搭便车从而加强局部开发能力。在原算法框架中种群被划分为数类角色一类是滚球蜣螂负责主搜索推进过程受历史位置和全局最差位置影响当它遇到障碍时会改变方向产生跳舞行为一类是繁殖蜣螂产卵它会根据当前局部最优位置动态调整产卵边界把卵产在最优解附近一类是觅食蜣螂小蜣螂负责在全局最优附近寻找食物最后一类是偷窃蜣螂会在全局最优和局部最优之间游走窃取资源。在我的实现里种群分配比例为滚球蜣螂占30%繁殖蜣螂占30%觅食蜣螂占20%偷窃蜣螂占20%。这个比例不是绝对标准你可以根据具体问题的搜索难度调整遇到高维问题时我会建议适当增加滚球蜣螂的比例加强全局探索能力。3.2 位置更新公式的数学表达这部分是代码的核心依据弄懂了公式实现就不难了。为了节省篇幅我用符号说明每个公式MATLAB代码会在第四章一一对应。滚球蜣螂无偏转推进的更新公式为x_i(t1) x_i(t) α × k × x_i(t-1) b × Δx其中Δx |x_i(t) - X_w(t)|X_w是当前全局最差位置α是方向系数随机取1或-1k是偏转系数取(0, 0.2]区间b是常数取(0,1)。遇到障碍时触发跳舞偏转更新公式为x_i(t1) x_i(t) tan(θ) × |x_i(t) - x_i(t-1)|其中θ在[0, π]内随机取值。这个公式的作用是让当前位置在原有步长基础上进行大角度扰动相当于跳出当前搜索方向。繁殖蜣螂的产卵区域由局部最优X*动态决定Lb* max(X* × (1 - R), Lb) Ub* min(X* × (1 R), Ub)其中R 1 - t / T_maxt是当前迭代次数T_max是最大迭代次数。随着迭代进行R减小产卵区间逐渐收缩体现了从全局搜索到局部精化的自然过程。产卵位置更新为x_i(t1) X* b1 × (x_i(t) - Lb*) b2 × (x_i(t) - Ub*)这里b1和b2是两个独立的随机向量。觅食蜣螂的觅食区域由全局最优X_b决定边界计算方法与产卵类似只是把X*换成X_b。更新公式为x_i(t1) x_i(t) C1 × (x_i(t) - Lbb) C2 × (x_i(t) - Ubb)其中C1服从正态分布随机数C2是0到1之间的随机向量。这种设计让觅食蜣螂在全局最优附近做随机步长搜索有助于精细开发。偷窃蜣螂的位置更新最简洁x_i(t1) X_b S × g × (|x_i(t) - X*| |x_i(t) - X_b|)其中S为常数g是标准正态分布随机向量。这个公式让偷窃蜣螂在全局最优位置附近进行随机扰动且扰动量与当前离最优点的距离成正比离得远步子大离得近步子小。3.3 算法主流程整个DBO的训练主循环是随机初始化种群位置每个个体对应一组传感器节点坐标向量计算初始适应度并记录全局最优位置X_b、全局最差位置X_w和局部最优位置X*进入迭代后依次对每个个体按照其角色类型执行对应的位置更新公式更新后做边界裁剪再重新计算覆盖率如果新位置的未覆盖率小于原位置就替换每个迭代周期结束后刷新全局最优、最差位置和动态参数R记录当前最优覆盖率。循环直到达到最大迭代次数。这里有一个需要特别注意的细节X理论上应该是局部最优位置的动态记录但为了简化实现并避免额外存储我在每次迭代开始时直接用当前种群适应度最小的位置作为X效果也不错。如果想把X*维护成一个长期记忆变量可以在迭代里额外增加一个数组保存每次迭代的最优位置然后随时间平滑更新。4. MATLAB实现与关键代码注释4.1 参数初始化下面这组参数是我在默认仿真场景下使用的注释都写在代码里你可以直接复制运行。% 场景参数 areaLen 100; % 监测区域长度m areaWid 100; % 监测区域宽度m numNodes 30; % 无线传感节点数量 sensingRadius 10; % 感知半径m gridStep 1; % 栅格评估步长m % DBO算法参数 popSize 30; % 种群规模 maxIter 100; % 最大迭代次数 dim numNodes * 2; % 每个个体维度30个节点的x、y坐标 % 变量边界 lb zeros(1, dim); % 下界全部为0 ub [ones(1, numNodes) * areaLen, ones(1, numNodes) * areaWid]; % 前N维x上界后N维y上界关于gridStep的取值我建议先用2来快速跑通流程确认算法稳定迭代后再改成1提升评估精度。栅格步长每减小一半栅格点数量会变成原来的约四倍对运行时间影响很大。4.2 覆盖率计算函数这是整个项目里最核心、也最需要高效实现的函数。我一开始用三层循环写得非常慢跑一次覆盖优化要几分钟后来改成meshgrid加向量化距离计算运行时间降到十秒级别优化效果立竿见影。function coverage calCoverage(nodes, areaLen, areaWid, sensingRadius, gridStep) % 栅格法计算覆盖率 % 输入 % nodes: numNodes x 2 矩阵每一行是一个节点的(x, y)坐标 % areaLen: 区域长度 % areaWid: 区域宽度 % sensingRadius: 感知半径 % gridStep: 栅格步长 % 输出 % coverage: 覆盖率0~1 % 生成所有栅格点坐标 xGrid 0:gridStep:areaLen; yGrid 0:gridStep:areaWid; [xx, yy] meshgrid(xGrid, yGrid); gridPoints [xx(:), yy(:)]; % 转为两列每一行是一个栅格点 totalPoints size(gridPoints, 1); % 初始化覆盖标记 isCovered false(totalPoints, 1); % 对每个节点把落在感知范围内的栅格点标记为覆盖 for k 1:size(nodes, 1) distK sqrt((gridPoints(:, 1) - nodes(k, 1)).^2 ... (gridPoints(:, 2) - nodes(k, 2)).^2); isCovered isCovered | (distK sensingRadius); end % 覆盖率 被覆盖栅格点数 / 总栅格点数 coverage sum(isCovered) / totalPoints; end实际使用时这个函数会被反复调用。为了进一步提速可以把gridPoints在外部生成一次作为参数传入避免每次适应度评估都重复生成网格点。我自己在最终版本里就是这么做的覆盖率计算函数的调用开销因此又降了一截。4.3 DBO主体搜索逻辑下面这段代码是整个算法的核心循环每一类蜣螂的位置更新都对应3.2节的一个公式。代码里做了完整的注释注意看position向量的前N维与后N维是如何重组为节点坐标矩阵的。% 初始化种群 positions rand(popSize, dim) .* (ub - lb) lb; % 预计算全局栅格点只做一次节省开销 xGrid 0:gridStep:areaLen; yGrid 0:gridStep:areaWid; [xx, yy] meshgrid(xGrid, yGrid); gridPoints [xx(:), yy(:)]; % 评估初始适应度 fitness zeros(popSize, 1); for i 1:popSize nodes_i reshape(positions(i, :), numNodes, 2); cover_i calCoverage(nodes_i, areaLen, areaWid, sensingRadius, gridStep); fitness(i) 1 - cover_i; % 用未覆盖率作为适应度越小越好 end % 记录全局最优与全局最差 [bestFitness, bestIdx] min(fitness); bestPosition positions(bestIdx, :); [~, worstIdx] max(fitness); worstPosition positions(worstIdx, :); % 收敛曲线 convergenceCurve zeros(maxIter, 1); % DBO主循环 for t 1:maxIter R 1 - t / maxIter; % 动态收缩因子 % 以当前种群最优位置作为局部最优X* [~, curBestIdx] min(fitness); Xstar positions(curBestIdx, :); for i 1:popSize x positions(i, :); % 滚球蜣螂前30% if i popSize * 0.3 if rand 0.9 % 无偏转推进 alpha (rand 0.5) * 2 - 1; % 方向系数随机1或-1 kk 0.1 0.1 * rand; % 偏转系数k(0, 0.2] bconst 0.3 * rand; % 常数b deltaX abs(x - worstPosition); % 与全局最差的距离 xnew x alpha * kk * positions(max(i - 1, 1), :) bconst * deltaX; else % 遇到障碍触发跳舞偏转 theta pi * rand; diffX abs(x - positions(max(i - 1, 1), :)); xnew x tan(theta) * diffX; end % 繁殖蜣螂30% ~ 60% elseif i popSize * 0.6 LbStar max(Xstar .* (1 - R), lb); UbStar min(Xstar .* (1 R), ub); b1 rand(1, dim); b2 rand(1, dim); xnew Xstar b1 .* (x - LbStar) b2 .* (x - UbStar); % 觅食蜣螂60% ~ 80% elseif i popSize * 0.8 LbB max(bestPosition .* (1 - R), lb); UbB min(bestPosition .* (1 R), ub); C1 randn(1, dim); C2 rand(1, dim); xnew x C1 .* (x - LbB) C2 .* (x - UbB); % 偷窃蜣螂最后20% else Sconst 0.5; g randn(1, dim); xnew bestPosition Sconst * g .* (abs(x - Xstar) abs(x - bestPosition)); end % 边界裁剪 xnew min(max(xnew, lb), ub); % 评估新位置覆盖率 nodes_new reshape(xnew, numNodes, 2); cover_new calCoverage(nodes_new, areaLen, areaWid, sensingRadius, gridStep); fitness_new 1 - cover_new; % 如果新位置更优则更新该个体 if fitness_new fitness(i) positions(i, :) xnew; fitness(i) fitness_new; end end % 迭代结束后更新全局最优与最差 [bestFitness, bestIdx] min(fitness); bestPosition positions(bestIdx, :); [~, worstIdx] max(fitness); worstPosition positions(worstIdx, :); % 记录当前最优覆盖率未覆盖率转回覆盖率 convergenceCurve(t) 1 - bestFitness; end这段代码在默认参数下可以稳定收敛。需要提醒一点tan(theta)在theta接近π/2时会变得非常大如果新位置超出边界剪裁操作会直接把它拉回边界这个机制保证了不会因为跳舞偏转产生越界节点。4.4 结果输出与可视化算法结束后把最优节点坐标提取出来画一张覆盖热力图能直观看到优化效果。可视化代码我一般这样写% 提取最优节点坐标 bestNodes reshape(bestPosition, numNodes, 2); figure; subplot(1, 2, 1); % 绘制覆盖区域 isCovered false(size(gridPoints, 1), 1); for k 1:numNodes distK sqrt((gridPoints(:, 1) - bestNodes(k, 1)).^2 ... (gridPoints(:, 2) - bestNodes(k, 2)).^2); isCovered isCovered | (distK sensingRadius); end coveredPoints gridPoints(isCovered, :); scatter(coveredPoints(:, 1), coveredPoints(:, 2), 4, g, filled); hold on; scatter(bestNodes(:, 1), bestNodes(:, 2), 60, r, filled); xlim([0 areaLen]); ylim([0 areaWid]); xlabel(x / m); ylabel(y / m); title(DBO优化后的节点部署与覆盖区域); subplot(1, 2, 2); plot(convergenceCurve, LineWidth, 2); xlabel(迭代次数); ylabel(覆盖率); title(覆盖率收敛曲线); grid on;实际运行时左侧散点图会展示优化后的节点如何分散在整个监测区域绿色区域是感知覆盖范围。收敛曲线则可以清楚看到覆盖率从初始随机部署的0.35左右逐步提升到0.75以上整条曲线的上升趋势在迭代20到40次之间最明显之后进入缓慢精调阶段。5. 仿真实验与调试心得5.1 典型场景实验设置我用默认场景做了多组实验监测区域100m × 100m节点数30感知半径10m栅格步长1m种群规模30最大迭代100次。随机初始部署的覆盖率通常在0.30到0.40之间带点运气成分因为30个半径为10m的圆在10000平方米区域内的理论覆盖面积上限大约是9500平方米覆盖率上限约0.95但实际能达到多少取决于节点分布方式。DBO优化后的结果在多次运行中稳定在0.75到0.82之间最好的一次达到了0.81。作为对比用PSO做同样的100次迭代平均只能到0.72左右遗传算法则要慢很多100代时还在0.68附近波动。DBO在覆盖优化这个问题上的优势不是体现在单次运气好而是多次运行方差小结果稳定。5.2 收敛曲线解读我建议你运行完程序后先看收敛曲线而不是急着看节点分布图。一个健康的收敛曲线应该是前20代快速上升后面逐渐平坦最后趋于一条水平线。如果曲线在前几代就拉平说明算法早熟大概率是种群多样性不足可以调大滚球蜣螂比例或增加种群规模。如果曲线到100代还在缓慢上升说明迭代次数不够可以增加到150代或200代。还有一种常见情况是曲线整体很低比如覆盖率达到0.5就已经趋势平缓。这时候要检查感知半径和节点数量是否匹配区域面积。30个节点、10m感知半径理论极限覆盖率并不高约束条件决定了解的上限换什么算法都突破不了这个物理边界。5.3 关键参数的影响分析我从实际调试中总结出几个参数对结果的影响规律直接列出来供参考。种群规模popSize是最直接的因素。过小小于15会让角色分配不完整蜣螂的多样性缺失过大超过50会显著拖慢运行速度而覆盖率提升有限。我常用30到40收益比最高。迭代次数maxIter不是越大越好。覆盖优化的大幅提升集中在前40代后面多是微调。如果每次实验时间允许我一般设120代左右追求稳健结果再增加到200代。传感半径R直接影响覆盖率的物理上限。R8时30个节点的覆盖面积总和理论上限约6000平方米覆盖率上限约0.6R10时上限约0.95。当节点数量少或半径小的时候覆盖问题更加陡峭优化算法很容易陷入局部最优此时可以适当增加偷窃蜣螂的比例让更多个体在全局最优附近搜索。6. 常见问题与排查技巧6.1 覆盖率为0或极低如果运行初始代码时发现覆盖率一直是0先检查坐标是否正确。一个容易踩的坑是把position向量reshape成节点矩阵时维度顺序写反比如写成reshape(position, 2, numNodes)而不是reshape(position, numNodes, 2)。这样前30维会被错误分配到另一个维度里节点坐标完全错乱。另一个可能是栅格步长设置太大比如gridStep20100m区域只有几十个栅格点覆盖率计算变得不具备代表性。解决办法是把步长调回1到2或者手动打印初始节点坐标在图上看看位置是否落在区域内。6.2 DBO不收敛或早熟我遇到最多的问题是收敛曲线在很低的覆盖率水平就变平也就是早熟。排查方向有两个一是分配比例问题如果滚球蜣螂比例过小全局探索不足整个种群会很快被繁殖和偷窃蜣螂拉到局部最优附近二是动态因子R的收缩速度R 1 - t/maxIter这个公式在早期收缩很快产卵和觅食区间迅速变小如果感知半径小、问题地形复杂算法会过早精细化。解决办法是给滚球蜣螂的随机决策中增加偏转概率或者适当加大种群规模增加搜索个体数量。6.3 运行速度太慢覆盖率函数的向量化是必须做的一步用两层for循环遍历几千个栅格点会非常痛苦。还有一个加速技巧是分阶段评估前30代用gridStep3的粗栅格跑等迭代后期再用gridStep1的细栅格精修。因为粗栅格的覆盖率虽然不够精确但能反映节点分布的大致优劣用于引导搜索方向足够最后切换到细栅格把解精确到米级。我用这个技巧把完整仿真时间从300多秒压到了90秒以内精度几乎无损。6.4 与其他算法对比时的注意事项做算法对比时评价指标不能只看最终覆盖率。一定要记录整个收敛过程并且固定使用相同的适应度评估方式包括相同栅格步长、相同初始化种子、相同种群规模和迭代次数。很多论文里对比结果看起来很明显但实际上是评估函数不一致造成的比如PSO用细网格而DBO用粗网格这会严重误导读者。我自己做对比时会写一个统一的前置脚本把评估函数抽成公共接口所有算法调用同一个覆盖率函数确保每一个百分点的提升都来自算法本身。最后再分享一个个人体会DBO在覆盖优化这类高维连续问题上的表现确实值得一试但它的性能上限很大程度取决于角色比例和动态参数的设计不是无脑用默认参数就能发财的。建议你在自己的场景里先跑通默认配置然后针对性地加一个参数扰动实验观察覆盖率变化这样才能真正理解算法行为也能给论文提供更有说服力的实验数据。