多策略黑猩猩优化算法MATLAB实现与工程应用详解

黑猩猩优化算法（Chimp Optimization Algorithm，ChOA）是这几年低成本启发式算法里讨论度比较高的一个，核心是模拟黑猩猩群体狩猎时的分工协作。最近我在复现黄倩那篇《多策略黑猩猩优化算法研究及其工程应用》时，把基本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 - f，C 为 2 * r2，r1 和 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.5；x 大于等于 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(1+beta) * sin(pi*beta/2) / ... (gamma((1+beta)/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 参数设置与运行建议

参数设置对启发式算法的影响很大，我在复现时整理了一组比较稳的默认值，遇到具体问题可以在这个基础上微调：

运行时有几个经验。第一，目标函数尽量写成行向量输入，不要用全局变量传数据，否则并行或批量测试时会很痛苦。第二，如果目标函数里涉及仿真或数值求解，建议在开头判断一下返回结果是否有限，一旦出现 NaN 或 Inf，直接返回一个很大的惩罚值，避免算法把异常解当成最优解。第三，多跑几次看收敛曲线，不要只跑一次就下结论。启发式算法的单次运行结果有一定随机性，至少独立运行 20 次再统计均值、标准差和最优值，这样写进报告里才有说服力。

4. 仿真结果与工程应用验证

4.1 基准测试函数对比

我在复现时选了四个经典测试函数：Sphere、Rastrigin、Ackley、Rosenbrock，维度统一设为 30，种群规模 30，最大迭代 500，每个算法独立运行 20 次。下表是一次典型运行的最优值对比，可以看出改进策略带来的变化：

这个结果和很多同类论文的趋势一致：改进算法在简单单峰函数上的提升幅度看不太出来，因为基本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 参数为 Kp=4.32、Ki=1.18、Kd=2.05，ITAE 大约 0.372，阶跃响应超调量在 12% 左右；多策略ChOA给出的参数为 Kp=5.86、Ki=1.52、Kd=2.74，ITAE 降到 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的基线，再逐步加入策略。哪个策略在你的目标函数上拖后腿，就把它调弱或者去掉，这才是多策略改进真正的意义。