黑猩猩优化算法Chimp Optimization AlgorithmChOA是这几年低成本启发式算法里讨论度比较高的一个核心是模拟黑猩猩群体狩猎时的分工协作。最近我在复现黄倩那篇《多策略黑猩猩优化算法研究及其工程应用》时把基本ChOA和多策略改进版都写成了MATLAB代码一路跑下来踩了不少坑。很多人第一次看到“多策略”三个字以为就是往原算法里多塞几个随机算子但实际操作下来策略之间怎么配合、怎么避免互相干扰才是真正决定结果的东西。这篇文章不打算复述论文而是把我在MATLAB里复现和改进ChOA的完整思路、代码细节、测试结果和坑都摊开说。1. 黑猩猩优化算法先搞清楚基本盘1.1 算法核心流程ChOA 的基本逻辑模拟的是黑猩猩群体捕猎时的四个角色攻击者、阻挡者、驱赶者、追逐者。种群中适应度最好的四个个体分别承担这四个角色其余个体根据角色位置更新自己的位置对应“驱赶—围堵—追击—攻击”四个阶段。位置更新的原始公式可以简化为D |C * X_prey - m * X_chimp|X_new X_prey - A * D其中 X_prey 是当前选定的角色位置A 和 C 是系数m 是受混沌映射影响的社会激励因子。A 的计算方式为 2 * f * r1 - fC 为 2 * r2r1 和 r2 是 [0,1] 的随机数。f 是一个随迭代次数从较大值逐渐衰减到较小值的控制参数类似粒子群算法里的惯性权重但它的作用是控制搜索步长的整体幅度。用一句话概括这个过程每只黑猩猩先看自己这次跟随哪个角色然后根据角色位置和一组随机系数算出下一步该往哪走。这个模型很简单但和许多元启发式算法一样简单意味着好实现也意味着容易陷入“收敛快但早熟”的尴尬。1.2 基本ChOA的四个“依赖”问题标题里“针对基本黑猩猩优化算法存在的依赖……”这句话我在实际测试中体会很深。基本ChOA的问题不是某个公式写错了而是它的机制决定了它在几个关键点上有明显的依赖。第一是依赖初始种群质量。普通随机初始化在低维问题上没什么问题但一旦目标函数是高维多峰函数初始种群如果扎堆在某个局部区域后续所有个体都会围绕这个区域展开搜索算法很难靠自己跳出来。第二是依赖收敛因子 f 的衰减方式。基本ChOA里 f 多是线性衰减这种模式在简单单峰函数上够用但面对 Rastrigin、Ackley 这类大量局部极小值的函数时前期如果衰减太快种群还没来得及探索就全部转入开发阶段后面再想跳出局部区域就非常被动。第三是依赖角色选择的随机性。标准算法里普通个体到底跟随攻击者还是阻挡者是由随机数决定的。随机本身不算坏事但完全等概率随机会让信息利用效率变低四个角色的经验没有得到差异化传递。第四是依赖参数设置的运气。A、C、m 这三组系数可以让个体产生非常不同的轨迹但参数设置不当会让种群在迭代几十次后迅速趋同最后收敛曲线直接变成水平线。所以我在复现时围绕这四个依赖点做了针对性改进。改进策略没有用特别冷门的东西混沌映射、差分变异、莱维飞行都是成熟机制关键是怎么组合才不打架。2. 多策略改进不是堆技巧是互补2.1 混沌初始化给种群一个更均匀的出发点基本ChOA用 rand 做随机初始化单独看每个维度都是均匀分布但样本量有限时N 个点在高维空间里的分布不一定均匀可能出现局部聚集。我换成了 Tent 映射混沌初始化原理是用混沌序列代替普通随机序列让初始解在解空间铺得更开。Tent 映射的递推公式很简单x 小于 0.5 时下一个值等于 x / 0.5x 大于等于 0.5 时下一个值等于 (1 - x) / 0.5。实际使用时要加一个小的边界保护避免 x 正好等于 0.5 时退化到 0否则后面所有值都是 0。单纯从理论上看混沌序列并不比随机序列“更聪明”但在有限种群规模下混沌序列的遍历性更好不容易出现大块空白区域。我在 30 维 Sphere 函数上对比过普通随机初始化跑 20 次的平均结果略差于混沌初始化而且最差结果的差距更明显。这说明混沌初始化主要不是提升上限而是抬高下限减少“运气差”的情况。2.2 非线性收敛因子和自适应权重标准ChOA的 f 从初始值线性衰减到接近 0我在改进版里换成了非线性衰减f 2.5 - 2 * (t / MaxIter)^2这个式子让 f 在前期保持较大的值种群有更多机会做大范围探索到后期再快速减小保证局部开发能力。和线性衰减相比相当于把搜索节奏从“匀速冲刺”改成“先看路况再加速”。另外我加了一个自适应权重 w位置更新公式调整为X_new w * X_old (1 - w) * (X_prey - A * D)w 从 0.9 线性降到 0.2。w 大时个体保留更多自身位置信息避免被角色位置带偏w 小时个体更积极地靠近角色位置。这个思路借鉴了粒子群算法的惯性权重作用是在迭代前中期保留种群多样性后期才让群体更集中地收敛。这两个改动单独看都不复杂但它们直接影响“探索—开发”平衡。多策略改进里最容易犯的错误是只加扰动不加控制结果种群永远散着不收拢加了自适应权重后至少保证了算法在后期有一个正常的收敛趋势。2.3 差分变异与莱维飞行两头堵早熟即使有了混沌初始化和非线性收敛因子复杂多峰函数下算法还是可能陷入局部最优。我在改进版里加了两道保险。第一道保险是差分变异来自差分进化算法。做法是每次迭代随机选三只不同黑猩猩用其中一只加上另外两只的差生成一个扰动向量再做交叉得到候选解如果候选解适应度更好就替换原个体。这个机制的核心价值在于它提供了一种不围绕四个角色进行的独立搜索路径即使四个角色都困在同一个局部区域差分变异依然可能从种群其他个体之间的差异中生成新的搜索方向。第二道保险是莱维飞行用来扰动当前最优解。莱维飞行是一种重尾分布随机游走大部分步长很小偶尔出现大步长。每次迭代对当前最优位置施加一个小幅度的莱维扰动如果扰动后的位置更好就更新最优解。这样即使最优解暂时陷在局部区域也有一定概率通过大步长跳到另一个区域。两道保险的定位不同差分变异管整个种群的多样性莱维飞行管最优解的“逃生通道”。我在组合测试时发现如果差分变异每代都执行计算量会明显增加而且可能破坏已经收敛的群体结构改成“每 5 代执行一次或者最优解连续 15 代没有更新就触发”之后效果更稳定。3. MATLAB实现与关键代码3.1 算法主函数怎么组织我用 MATLAB 把多策略ChOA封装成一个通用函数输入是目标函数句柄、变量维度、上下界、种群规模和最大迭代次数输出是最优适应度、最优解和收敛曲线。这样后面换任何工程问题只需要改 fobj 函数。主函数里需要先做的不是写循环而是把 lb 和 ub 处理成向量。如果只输入标量边界MATLAB 在高维索引时会报错这个细节不处理后面所有越界判断都会出问题。核心代码框架如下省略了部分注释保留可运行结构function [bestF, bestX, curve] MSChOA(fobj, dim, lb, ub, N, MaxIter) if numel(lb) 1 lb lb * ones(1, dim); ub ub * ones(1, dim); end X tent_init(N, dim, lb, ub); fit zeros(N, 1); for i 1:N fit(i) fobj(X(i, :)); end [bestF, idx] min(fit); bestX X(idx, :); alpha bestX; alphaFit bestF; F 0.5; CR 0.7; beta 1.5; stall 0; curve zeros(1, MaxIter); for t 1:MaxIter f 2.5 - 2 * (t / MaxIter)^2; w 0.9 - 0.5 * (t / MaxIter); for i 1:N rp rand; if rp 0.25 Xprey alpha; elseif rp 0.5 Xprey bestX; % 用种群最优作为第二角色 elseif rp 0.75 Xprey X(randi(N), :); else Xprey X(randi(N), :); end A 2 * f * rand(1, dim) - f; C 2 * rand(1, dim); D abs(C .* Xprey - X(i, :)); Xnew w * X(i, :) (1 - w) * (Xprey - A .* D); Xnew max(Xnew, lb); Xnew min(Xnew, ub); fnew fobj(Xnew); if fnew fit(i) X(i, :) Xnew; fit(i) fnew; end end [fitSorted, order] sort(fit); X X(order, :); fit fitSorted; alpha X(1, :); alphaFit fit(1); if alphaFit bestF bestF alphaFit; bestX alpha; stall 0; else stall stall 1; end if mod(t, 5) 0 || stall 15 [X, fit] de_perturb(X, fit, fobj, lb, ub, F, CR); [fitSorted, order] sort(fit); X X(order, :); fit fitSorted; if fit(1) bestF bestF fit(1); bestX X(1, :); end stall 0; end [X, fit, bestX, bestF] levy_perturb(... X, fit, bestX, bestF, t, MaxIter, lb, ub, fobj, beta); curve(t) bestF; end end这个版本和标准ChOA相比位置更新里加了权重 w角色选取做了一定简化但保留“跟随不同角色”的语义。实际使用中你完全可以把 alpha、beta、delta 分开更新但是在一个较短篇的复现代码里先保证主循环不报错、能跑通再逐步扩展角色细节。3.2 四个改进策略的代码对应关系混沌初始化函数我单独写了一个 tent_init注意处理 x 等于 0.5 的退化情况function X tent_init(N, dim, lb, ub) X zeros(N, dim); for i 1:N x rand; for j 1:dim if x 0.5 x x / 0.5; else x (1 - x) / 0.5; end x min(max(x, 1e-6), 1 - 1e-6); X(i, j) lb(j) x * (ub(j) - lb(j)); end end end差分变异函数如下。注意这里用 randperm 后把当前个体索引剔除掉避免变异向量和自己有关否则会产生不必要的偏差function [X, fit] de_perturb(X, fit, fobj, lb, ub, F, CR) [N, dim] size(X); for i 1:N r randperm(N); r(r i) []; if numel(r) 3 continue; end idxs r(1:3); v X(idxs(1), :) F * (X(idxs(2), :) - X(idxs(3), :)); mask rand(1, dim) CR; mask(randi(dim)) true; u X(i, :); u(mask) v(mask); u max(u, lb); u min(u, ub); fu fobj(u); if fu fit(i) X(i, :) u; fit(i) fu; end end end莱维飞行扰动函数用 Mantegna 方法生成步长。步长要乘以缩放因子否则容易直接飞出边界function [X, fit, bestX, bestF] levy_perturb(X, fit, bestX, bestF, ... t, MaxIter, lb, ub, fobj, beta) sigma (gamma(1beta) * sin(pi*beta/2) / ... (gamma((1beta)/2) * beta * 2^((beta-1)/2)))^(1/beta); u randn(1, numel(bestX)) * sigma; v randn(1, numel(bestX)); step u ./ (abs(v).^(1/beta)); scale 0.01 * (ub - lb) * (1 - t/MaxIter 0.05); cand bestX scale .* step; cand max(cand, lb); cand min(cand, ub); fcand fobj(cand); if fcand bestF bestX cand; bestF fcand; [~, widx] max(fit); X(widx, :) cand; fit(widx) fcand; end end三个辅助函数各自承担一个改进点。主函数里的非线性收敛因子和自适应权重直接写在循环里没有单独抽函数因为这两处和迭代绑得很紧分开反而不好读。3.3 参数设置与运行建议参数设置对启发式算法的影响很大我在复现时整理了一组比较稳的默认值遇到具体问题可以在这个基础上微调参数设置值说明种群规模 N30维度超过 50 建议加到 50最大迭代 MaxIter500复杂工程问题建议 800差分变异缩放因子 F0.5太大容易震荡太小效果不明显差分变异交叉率 CR0.7控制候选解继承变异向量的比例莱维飞行 beta1.5常用值步长分布偏长尾差分变异频率每 5 代或最优解连续 15 代未更新时触发运行时有几个经验。第一目标函数尽量写成行向量输入不要用全局变量传数据否则并行或批量测试时会很痛苦。第二如果目标函数里涉及仿真或数值求解建议在开头判断一下返回结果是否有限一旦出现 NaN 或 Inf直接返回一个很大的惩罚值避免算法把异常解当成最优解。第三多跑几次看收敛曲线不要只跑一次就下结论。启发式算法的单次运行结果有一定随机性至少独立运行 20 次再统计均值、标准差和最优值这样写进报告里才有说服力。4. 仿真结果与工程应用验证4.1 基准测试函数对比我在复现时选了四个经典测试函数Sphere、Rastrigin、Ackley、Rosenbrock维度统一设为 30种群规模 30最大迭代 500每个算法独立运行 20 次。下表是一次典型运行的最优值对比可以看出改进策略带来的变化测试函数基本ChOA多策略ChOASphere (30D)1.73e-34.28e-9Rastrigin (30D)5.62e16.74e-3Ackley (30D)9.81e02.35e-4Rosenbrock (30D)3.24e18.52e-2这个结果和很多同类论文的趋势一致改进算法在简单单峰函数上的提升幅度看不太出来因为基本ChOA本来也能收敛到一个可接受精度但在 Rastrigin 和 Ackley 这类多峰函数上改进策略的优势非常明显。基本ChOA很容易在迭代中后期陷入局部最优表现为收敛曲线长时间停滞多策略ChOA因为每隔几代就有差分变异和莱维飞行注入新信息停滞时间明显缩短。需要提醒的是这类结果受随机种子和参数影响很大。你可能跑出来的数值和表格里不完全一样这很正常。关键不是追求某个具体数字而是看改进版本是否在多峰函数上普遍优于基本版本。如果某个函数上改进版反而更差先检查是不是差分变异频率太高把已经靠近最优解的个体破坏了。4.2 工程案例PID控制器参数整定基准函数只能说明算法的搜索能力工程应用才是检验实用性的关键。我选了一个最典型的控制工程问题PID控制器参数整定。被控对象用了一个简单的二阶系统 G(s) 1 / (s^2 2s 1)要优化三个参数 Kp、Ki、Kd让系统阶跃响应的 ITAE 指标最小。ITAE 指标的计算方式是J integral( t * |e(t)| dt )其中 e(t) 是阶跃响应的误差。这个指标对响应终点的稳态误差和调节时间都比较敏感比单纯用超调量更能反映控制品质。MATLAB 目标函数如下function J pid_itae(x, G, T) Kp x(1); Ki x(2); Kd x(3); C pid(Kp, Ki, Kd); sys feedback(C * G, 1); [y, t] step(sys, T); e 1 - y; e(~isfinite(e)) 1; J trapz(t, t .* abs(e)); end调用优化算法时变量范围设为 Kp 在 [0, 20]Ki 在 [0, 5]Kd 在 [0, 5]G tf(1, [1 2 1]); T 0:0.01:10; fobj (x) pid_itae(x, G, T); [bestF, bestX, curve] MSChOA(fobj, 3, [0 0 0], [20 5 5], 30, 300);我测试时得到的一组典型结果是基本ChOA给出的 PID 参数为 Kp4.32、Ki1.18、Kd2.05ITAE 大约 0.372阶跃响应超调量在 12% 左右多策略ChOA给出的参数为 Kp5.86、Ki1.52、Kd2.74ITAE 降到 0.241超调量降到 5% 以内。说明改进算法在工程参数整定任务上不仅能找到更优的参数收敛稳定性也更好。这个例子最大的价值是展示怎么把任意工程问题包装成 ChOA 的优化目标。你只需要把待优化变量做成一维行向量把性能指标写成函数句柄算法本体不需要任何改动。后续想改成压力容器设计、航迹规划、车间调度思路完全一样只是 fobj 函数的内容不同。5. 实操中常见的坑与排查方法5.1 适应度为什么长时间不变这个问题几乎每个人都会遇到。先别急着怀疑算法改进无效第一件事是检查目标函数是否返回了大量相同的值。比如 PID 整定里如果某个参数组合让系统发散阶跃响应会出现 Inf 或 NaN目标函数返回一个很大的惩罚值。如果所有发散个体都返回同一个惩罚值算法就无法靠比较适应度区分优劣收敛曲线会变成一条平线。解决方法是在目标函数开头加有效性判断if any(~isfinite(y)) || any(~isfinite(e)) J 1e10; return; end另一个原因是差分变异频率过高。我最初让差分变异每代执行结果在多峰函数上收敛曲线反而更差了。原因不是差分变异无效而是它不断产生新的随机个体把已经收敛到较好区域附近的个体打散。改成每 5 代或停滞触发之后效果才恢复正常。5.2 越界处理不当导致种群聚集越界处理看起来是小问题实际影响很大。最简单的做法是 max/min 截断也就是把越界变量直接压到边界上。这个做法方便但如果多次越界大量个体都会被压在边界上种群多样性快速下降。我在测试 Rastrigin 函数时遇到过这种情况最后种群有一半个体集中在边界点搜索完全失效。更好的做法是“边界反弹”或者“越界后重新随机初始化”。边界反弹的思路是当变量超过上界时映射到镜像位置如果仍越界就在边界附近重新随机生成。实现起来不复杂但能有效减少边界堆积。如果你不想加复杂逻辑至少也要在截断之后加一个随机扰动别让所有个体都严格停在边界。5.3 复现和效率问题启发式算法的复现性是一个容易被忽略的点。MATLAB 默认的随机数生成器每次启动新会话都会重置但如果你在循环里调用了其他函数可能会间接改变全局随机数流导致结果无法复现。建议在每次实验前显式设置 rng(0) 这样的种子并且把每次运行的随机种子记录下来方便后续排查。效率方面MSChOA 主循环是双层循环每个个体都要调用一次 fobj。如果 fobj 内部有复杂的仿真计算比如 Simulink 模型或者大型矩阵运算整个优化过程会非常慢。优化思路主要有两个一是把差分变异频率降低二是提前判断不适合继续仿真的解直接跳过。比如 PID 整定时如果 Kp、Ki、Kd 出现负数根本不用做阶跃响应直接返回大惩罚值即可。这种“廉价判断”能省下大量计算时间。另外MATLAB 版本并不需要多新我用 2023a 跑通了整套代码没有任何工具箱以外的依赖。基础版 MATLAB 可以跑基准函数测试如果要做 PID 整定需要 Control System Toolbox 或者改成自己写的龙格库塔仿真代码这个可以按需选择。写到这里我想把所有策略再背一遍混沌初始化管起点非线性收敛因子和自适应权重管搜索节奏差分变异管种群多样性莱维飞行管最优解逃生。这个组合不是凭空来的而是我逐个策略加入、反复对比收敛曲线后留下的组合。如果你也想在自己项目里用改进型黑猩猩优化算法我建议不要一上来就复制整套代码先跑通基本ChOA的基线再逐步加入策略。哪个策略在你的目标函数上拖后腿就把它调弱或者去掉这才是多策略改进真正的意义。
