简介:群体智能优化算法通过模拟生物群体行为解决复杂连续优化问题,其中勘探与开发的平衡往往决定算法性能。斑马算法(ZBO)作为新型群体智能方法,模仿斑马的迁徙与反捕食策略,利用双阶段更新机制实现全局探索与局部精修。它不依赖梯度信息,在多峰、不可导问题上具备独特优势,且MATLAB实现简洁,仅需数十行主循环即可运行。然而,算法效果受参数设置与边界处理策略影响显著,常见的随机重置越界方法与线性衰减系数能有效避免早熟收敛。本文从算法原理出发,逐步讲解MATLAB代码实现,并给出基于基准函数的实验设计、统计对比与结果解读方法,帮助读者科学评估算法性能,为实际工程优化提供可靠参考。
1. 斑马算法不是“换了马甲的遗传算法”:先搞清楚它在解什么题
斑马算法(Zebra Optimization Algorithm,ZBO)是2022年提出的一种群体智能算法,它把斑马群的生存策略拆成了两件事:发现捕食者时跟着领队跑,狮子扑过来时各自沿“之”字形路线逃。这个双阶段模型放在连续优化问题上,正好对应群体智能里最核心的“勘探”与“开发”博弈,而且用MATLAB实现只需要几十行主循环。如果你遇到的是多峰、不可导、无法用梯度求解器的目标函数,又不想在matlab优化工具箱和fmincon的局部解上反复试初值,斑马算法是个值得评测的备选。文章会从策略公式讲起,给一版可以直接跑的MATLAB代码,再把“测试效果”这件事拆解成实验设计、统计指标和坑点规避,确保你跑出的收敛曲线真的能说明问题。
2. 斑马算法的机理拆解:从“反捕食行为”到可编码的迭代公式
斑马算法与遗传算法、粒子群(PSO)最大的区别,不在于“有没有种群”,而在于位置更新不是通过交叉变异或速度惯性完成的,而是用一次随机数直接把当前个体切换到两套截然不同的更新策略上。一套偏向全局侦察,一套偏向局部逃离,中间用一个概率门控切换。理解了这个门控,后面读代码时就不会觉得if-else是随手写的,调参数时也能判断到底该动哪一个变量。
2.1 勘探与开发两阶段:领队引导公式和“之”字逃跑公式怎么来的
勘探阶段模拟的是侦察者发现捕食者后的群体迁徙行为。群体中有一匹斑马在觅食时识别到危险,其他斑马会朝当前全局最优个体(领队)所在方向移动,移动的同时保留一定过冲量。简化后的更新式可以写成:
x_new = x_old + r * (x_best - c * x_old)
这里的r是(0,1)区间均匀分布的随机数,c是学习系数,它直接控制过冲量。c > 1时,x_best - c * x_old 这个差值会让当前个体越过最优位置,从而在全局最优附近形成更大的搜索半径;c < 1时,差值被压缩,个体向最优位置收缩,局部搜索特性更强。常见的做法是让c从2.0线性递减到1.0,这样前期勘探步长大,后期逐渐转入精修。
开发阶段对应狮子发动攻击后斑马群的逃跑行为。原始策略里斑马会向与捕食者进攻方向垂直的路线逃散,这样可以拉大与攻击线的距离。但在高维连续优化中,计算“垂直方向”需要额外构造空间向量,成本不低且维度越高收益越不稳定。工程实现里更常见的是把逃跑目标设为一个随机点:
x_new = x_old + r * (x_rand - x_old)
其中x_rand是搜索空间内的随机位置。它同样能实现“远离当前陷阱区域”的效果,而且不会引入越界修正的额外计算。不要觉得这个简化离论文太远,群体智能算法的工程价值本来就在“机制可复现、效果可统计”,而不是逐行复刻论文公式。
2.2 关键参数配置:种群规模、寻优维数、迭代次数怎么配
参数选择直接决定算法是否能在有限预算内收敛。给出我平时使用的参考范围:
| 参数 | 符号 | 参考值 | 作用 |
|---|---|---|---|
| 种群规模 | N | 20 ~ 50 | 维度高时建议取 2d ~ 4d |
| 最大迭代次数 | T | 500 ~ 2000 | 与目标函数评测预算挂钩 |
| 勘探概率 | P | 1.0 线性衰减至 0.3 | 前期多全局搜索、后期局部精修 |
| 学习系数 | c | 2.0 递减至 1.0 | 控制过冲量和步长 |
需要注意,zigbee优化算法这类种群算法每代都会调用 N 次目标函数,总评测次数约等于 N*T。假如你只有5000次函数评测预算,就不要机械地设N=50、T=1000,而应该把预算消耗控制在同一个数量级再和其他算法对比。
2.3 常见误用:把斑马算法当全局黑箱,靠灵感改代码
最容易踩的坑有三个。第一,边界处理随意裁剪,导致大量个体堆在边界上,种群多样性快速丧失;第二,P和c设置成固定值,前期勘探不足或后期震荡不收敛,然后误判为“算法不行”;第三,拿一次运行的适应度值直接下结论。斑马算法在两阶段门控下带有明显的随机性,单次运行的最优值可能相差几个数量级,正确做法是多次独立重复并做统计检验,这部分在第4章会展开。
3. 用MATLAB写一版可复现的斑马算法:函数文件、主循环与边界处理
3.1 目标函数接口设计:让算法函数不绑定具体问题
我把目标函数约定为“接收一行向量、返回一个标量”的格式。这样无论是Sphere、Rastrigin还是实际工程里套了仿真器的目标函数,都能通过函数句柄接入。标准结构如下:
function y = sphere(x) % x 为 1 x dim 行向量 y = sum(x.^2); end这种接口设计的好处是,算法主函数不需要关心目标函数内部逻辑。你可以在脚本里用@(x)定义匿名函数,免去反复创建文件的麻烦,例如:
fitFun = @(x) sum(x.^2);3.2 斑马算法主循环实现:勘探-开发两个阶段的更新逻辑
下面代码是完整可运行的斑马算法函数。整体结构包含初始化、主循环、边界修正三个阶段。勘探阶段加入“随机牺牲”算子模拟被捕食个体被随机重置,能有效防止种群过早集中。
function [gbest, gfit, curve] = zbo_opt(fitFun, dim, lb, ub, N, T) % zbo_opt: Zebra Optimization Algorithm % 输入: % fitFun - 适应度函数句柄,接收1*dim向量,返回标量 % dim - 搜索维度 % lb, ub - 搜索空间下界与上界,可传标量或向量 % N - 种群规模 % T - 最大迭代次数 % 输出: % gbest - 全局最优个体 % gfit - 全局最优适应度 % curve - 收敛曲线,T x 1 if numel(lb) == 1 lb = repmat(lb, 1, dim); ub = repmat(ub, 1, dim); else lb = lb(:)'; ub = ub(:)'; end % 初始化种群 X = rand(N, dim) .* (ub - lb) + lb; fit = zeros(N, 1); for i = 1:N fit(i) = fitFun(X(i, :)); end [gfit, idx] = min(fit); gbest = X(idx, :); curve = zeros(T, 1); for t = 1:T % P 线性衰减: 从 1.0 降到 0.3 P = 1.0 - (t / T) * 0.7; % c 线性衰减: 从 2.0 降到 1.0 c = 2.0 - (t / T); for i = 1:N if rand() < P % 勘探阶段: 向全局最优靠拢,c控制过冲 X_new = X(i, :) + rand() * (gbest - c * X(i, :)); else % 开发阶段: 向随机目标点移动,保持搜索多样性 target = lb + (ub - lb) .* rand(1, dim); X_new = X(i, :) + rand() * (target - X(i, :)); end % 随机牺牲算子: 个体直接被重置在搜索空间内 if rand() < 0.1 X_new = lb + (ub - lb) .* rand(1, dim); end % 越界处理: 越界分量随机重置,而不是裁剪 out_low = X_new < lb; out_up = X_new > ub; X_new(out_low) = lb(out_low) + rand(1, sum(out_low)) .* (ub(out_low) - lb(out_low)); X_new(out_up) = lb(out_up) + rand(1, sum(out_up)) .* (ub(out_up) - lb(out_up)); fnew = fitFun(X_new); if fnew < fit(i) X(i, :) = X_new; fit(i) = fnew; end if fnew < gfit gfit = fnew; gbest = X_new; end end curve(t) = gfit; end end这段代码需要重点看三个位置。第一是P的衰减公式P = 1.0 - (t / T) * 0.7,这意味着在整个迭代过程中,勘探阶段的执行概率从100%逐步降到30%,算法不会在后期还频繁进行大范围跳变。第二是越界处理中没有使用X_new = min(max(X_new, lb), ub)这种裁剪写法,原因是裁剪会把大量个体推到边界上,在多峰函数中很容易让种群失去多样性。第三是随机牺牲算子固定在10%概率,这个值不要给太大,否则算法会退化成随机搜索,收敛曲线呈现锯齿状。
3.3 MATLAB运行与测速:从脚本跑通到观察收敛行为
在脚本中调用上述函数,检查是否能够正常收敛:
lb = -5.12; ub = 5.12; dim = 30; N = 50; T = 500; % Sphere 函数,理论最优值 0 fitFun = @(x) sum(x.^2); [gbest, gfit, curve] = zbo_opt(fitFun, dim, lb, ub, N, T); fprintf('Best fitness: %.6e\n', gfit); % 绘制收敛曲线 figure; semilogy(1:T, curve, 'LineWidth', 1.5); grid on; xlabel('Iteration'); ylabel('Best Fitness'); title('ZBO Convergence on Sphere');这个流程在MATLAB R2023b到目前的新版本里都不需要额外安装工具箱,直接保存函数文件后运行脚本即可。如果你的目标是快速验证算法是否存在遍历性不足的问题,可以把N调大到80、T减少到200,对比收敛曲线中是否存在长时间的平台期。需要提醒的是,MATLAB循环中逐行调用fitFun在N和T上升到一定量级后会产生可感知的耗时,工程应用时可以先把目标函数向量化,再对X矩阵按行批量计算适应度。
4. 测试效果不是“画张收敛图”:实验设计、指标与对比基线
4.1 基准函数与实验配置:Sphere、Rastrigin、Griewank怎么选
选择基准函数的核心原则是覆盖不同类型的优化难度。单峰函数用于验证算法基本收敛能力,多峰函数用于检验全局勘探能力,不可导函数用于确认算法不依赖梯度信息也能下降。表里是四个常用函数及其配置:
| 函数 | 表达式 | 搜索范围 | 理论最优值 |
|---|---|---|---|
| Sphere | sum(x_i^2) | [-100, 100]^d | 0 |
| Rastrigin | 10d + sum(x_i^2 - 10cos(2πx_i)) | [-5.12, 5.12]^d | 0 |
| Ackley | 见代码 | [-32, 32]^d | 0 |
| Griewank | 1 + sum(x_i^2/4000) - prod(cos(x_i)/sqrt(i)) | [-600, 600]^d | 0 |
Ackley函数的MATLAB实现如下:
function y = ackley(x) d = numel(x); sum1 = sum(x.^2); sum2 = sum(cos(2 * pi * x)); y = -20 * exp(-0.2 * sqrt(sum1 / d)) - exp(sum2 / d) + 20 + exp(1); end如果你只是想确认代码没有逻辑错误,只跑Sphere就够;但想说明算法“测试效果”,Rastrigin更容易暴露问题,因为它的局部最优密度高,P和c衰减过快时曲线会在半途拉平。Griewank的特点是高维下函数表面相对平坦且带有周期扰动,边界处理策略的不同会带来明显差异,适合用来考察越界重置与边界吸收的效果差别。
4.2 与PSO、GWO、WOA做公平对比:预算、重复次数与Wilcoxon检验
算法效果对比最容易犯的错误是“各自用各自的迭代次数”。正确做法是把总预算对齐到相同的目标函数评测次数,例如统一N=40、T=500,所有算法使用相同初始边界,并且各自独立运行30次。下面脚本骨架展示了对比流程:
N = 40; T = 500; dim = 30; runs = 30; % 存放各算法每次运行的最优值 bestZBO = zeros(1, runs); bestGWO = zeros(1, runs); for r = 1:runs rng(r, 'twister'); % 复现实验 [~, vZBO, ~] = zbo_opt(@(x) sum(x.^2), dim, -100, 100, N, T); [~, vGWO, ~] = gwo_opt(@(x) sum(x.^2), dim, -100, 100, N, T); bestZBO(r) = vZBO; bestGWO(r) = vGWO; end % Wilcoxon 秩和检验 p = ranksum(bestZBO, bestGWO); fprintf('p-value: %.4f\n', p);代码中rng(r, 'twister')确保第r次运行时随机数流可复现,这对排查“某次效果特别好”的偶然因素很重要。ranksum是MATLAB基础函数,不需要统计工具箱,它会返回两个分布存在显著差异的置信度。只有当p值小于0.05时,才可以说斑马算法在该问题上与对比算法存在统计意义上的差别,否则所谓“效果更好”只是采样噪声。
4.3 结果解读:收敛曲线、箱型图与数字表格里到底该看什么
多峰函数上不同算法的最终结果可能相差多个数量级,需要灵活使用多种指标来展示。重复运行后建议输出三样东西:收敛曲线的中位数曲线,最终适应度分布的箱型图,以及均值、标准差、最优值、最差值四列数字。箱线图的作用是直接展示30次运行的离散程度,单纯画收敛曲线只能看到一条线,看不出稳定性。绘制箱型图可以直接使用MATLAB内置的boxplot,无需额外工具箱:
figure; boxplot([bestZBO', bestGWO'], {'ZBO', 'GWO'}); ylabel('Best Fitness'); grid on;这里给一个经验性的预期结果,不是固定结论:在Sphere这类单峰函数上,ZBO的收敛曲线通常会在迭代初期就压到1e-50以下,这其实是双精度浮点表示的下限,已经不需要继续比较精度;而在Rastrigin函数上,ZBO的前期勘探能力决定了能否跳出第一个局部陷阱,P如果线性衰减过快,曲线会呈现“台阶状”,每一段平台对应一个局部最优。遇到这种曲线不需要立刻调小P,更多情况下是c值过高导致个体反复越过目标区域,把c的初始值从2.0降为1.5再看差异。
5. 让算法在实际问题里“不翻车”的四个参数技巧
5.1 越界处理策略:随机重置与边界吸收的实际差异
边界裁剪写法简洁,但在边界附近会形成大量不可逆的“粘附”,尤其当最优解恰好位于搜索空间边缘时,裁剪策略会让种群丧失向内部探索的能力。推荐的方法是对越界分量的每个坐标独立随机重置,在代码里用X_new(out_low) = lb(out_low) + rand(1, sum(out_low)) .* (ub(out_low) - lb(out_low))处理。如果问题本身要求解必须严格落在边界上,再考虑裁剪方案。
5.2 停滞检测与局部逃逸:遇到多模态函数的正确应对
连续几十代最优适应度没有下降,基本可以判定陷入了局部最优。这时不必整篇重启算法,可以把当前种群按最优个体为中心重新分布:
if t > 50 && curve(t) == curve(t - 10) X = gbest + (rand(N, dim) - 0.5) * (ub - lb) * 0.2; X = max(min(X, ub), lb); end这段代码把种群重新铺在最优解周围一个较小的邻域内,保留已找到的解信息,同时给个体一个逃逸机会。这个方法在Rastrigin这类谷底密集的函数上有效,但在平滑单峰函数上会打乱正常的收敛节奏,使用前需要先判断目标函数的形态。
5.3 把固定学习系数改成自适应,收敛曲线不再“拉链化”
固定c=1.5时,迭代后期全局最优附近会出现频繁的震荡,适应度在小范围内上下跳动,曲线看起来像拉链。把c改成随时间递减后,后期步长自动缩小,曲线会平滑得多。我常用的替换形式是:
c = 1.5 - 0.5 * (t / T);它从1.5降到1.0,既保留了中期过冲能力,又保证后期能稳定落在最优附近。这个改动对高维问题尤其明显,维度越高,过冲带来的震荡代价越大。
5.4 验证实现的三个小手段:单峰函数收敛测试、重复运行统计、维度压力测试
第一,先用Sphere函数验证算法逻辑是否有问题,看能否在100代内降到1e-20以下;第二,把同一个测试函数重复运行30次,采集分布而不是单点值;第三,把维度从10逐渐升到30、50,观察收敛曲线变化,如果维度升高后曲线反而变平滑,往往说明前期勘探不足,需要调大P值的初始区间或增加随机牺牲概率。这样验证过后再上真实工程问题,才算是把斑马算法的测试效果做完整了。
本文还有配套的精品资源,点击获取