灰色狼群优化GWO的MATLAB完整实现与源码详解

简介：面向机器学习与智能优化方向的开发者，这份资源提供灰色狼群优化（GWO）算法的MATLAB源码及详细中文注解，可用于支持向量机（SVM）与支持向量回归（SVR）参数寻优，帮助读者从算法原理过渡到工程实现。压缩包共9个文件，包含6个m脚本、2个mat数据集和1个PDF理论文档，整体大小约1.49MB，轻量易用。m脚本分别演示SVM与SVR的完整优化流程，涵盖狼群初始化、领导狼更新、位置迭代和适应度计算，并配套MSE、MAE、MAPE等误差评估函数；mat文件为可直接加载的测试数据集，PDF文档则对算法原理做了系统梳理。代码注释详细、模块划分清晰，既适合初学者快速理解GWO机制，也可作为论文实验或项目验证的参考工具，省去从零编写优化算法的成本。已有6680人学习，是实践元启发式算法与机器学习结合的实用资源。 第一次把灰色狼群优化（GWO）写进MATLAB，是给一个参数辨识项目找替代方案的时候。当时遗传算法和粒子群算法我都调过，GA参数多、调优窗口长得让人头疼，PSO又时不时早熟。GWO给我的第一印象是：结构简单到不像进化算法，但实际跑起来收敛速度和稳定性都相当能打。所以这次我把带中文注解的MATLAB源码完整整理出来，从狼群初始化、三头领导狼的位置更新到边界处理，每一行关键代码都拆开讲清楚。这篇文章适合三类人：正在写智能算法课程作业的学生、想给项目换一个免梯度优化器的工程师、以及准备复现GWO论文但又不想从零抠公式的研究者。

1. 灰色狼群优化到底在做什么：算法思路与适用场景

1.1 三头领导狼和“社会等级”在优化里是什么含义

GWO模仿灰狼群体的社会等级和捕食行为。算法里每一只狼代表候选解，狼群会按照适应度值排名：当前最优解叫Alpha，第二叫Beta，第三叫Delta，剩下的全部叫Omega。搜索过程中，Omega狼不直接向真正的猎物位置靠近，而是跟随Alpha、Beta、Delta三头领导狼的方向更新。

这里有个很直觉的设计逻辑：如果只跟着最优解Alpha走，群体容易一头撞进局部最优；如果把前三名都纳入参考，相当于班里其他同学不是只抄第一名作业，而是把前三名的思路各抄一份再取平均。这样就算第一名思路暂时卡住，第二第三名也能把你往别的方向带。这个“多头指导”机制就是GWO在收敛速度和多样性之间取得平衡的核心原因。

实际编写代码时，Alpha、Beta、Delta三个位置变量会随着每一次迭代不断刷新。要注意的是，这三头狼只代表当前种群里的前三名，而不是历史全局最优；历史最优信息其实是通过“保留Alpha_score”这种形式间接存在的。明白这一点，后续写排序和更新逻辑才不会乱。

1.2 包围、猎捕与攻击三个公式到底改了什么

GWO的原始数学模型分三个阶段：包围猎物、猎捕猎物、攻击猎物。包围阶段的核心是距离公式：

D = |C * X_p(j) - X(i,j)|

位置更新公式是：

X_new = X_p(j) - A * D

其中A = 2a*r1 - a，C = 2*r2，r1和r2是[0,1]之间的随机数，a从2线性降到0。这两个公式看起来简单，但含义很深。

A的绝对值决定狼是“靠近”还是“走远”。当|A| < 1时，狼会向领导狼方向收缩，对应攻击猎物，也就是局部开发；当|A| > 1时，狼会偏离目标，对应探索新区域。C是一个随机权重，用来模拟捕食过程中遇到的障碍和风向扰动，它可以避免算法因为过于规律而陷入死循环。早期迭代a接近2时，A取值范围大，狼群会到处探索；后期a接近0时，A收敛到小范围，狼群集中攻击。这就是为什么a的衰减方式直接影响算法表现，后面我会专门展开。

2. 搭建可复用的GWO函数：MATLAB源码逐行注解

2.1 主函数GWO.m完整源码

以下是我在实际项目中整理过的版本，函数签名、输入输出都做了规范化，可以直接复制到你的工程目录里使用。

function [Best_score, Best_pos, Convergence_curve] = GWO(func, dim, LB, UB, SearchAgents_no, Max_iteration) % ============================================================= % GWO 灰色狼群优化算法主函数 % 输入： % func : 适应度函数句柄，例如 @(x) sum(x.^2) % dim : 决策变量维度 % LB : 下界，标量或1*dim向量 % UB : 上界，标量或1*dim向量 % SearchAgents_no : 灰狼个数（种群规模） % Max_iteration : 最大迭代次数 % 输出： % Best_score : 历史最优适应度值 % Best_pos : 最优解 % Convergence_curve : 每次迭代的最优适应度，画收敛曲线用 % ============================================================= if numel(LB) == 1 LB = repmat(LB, 1, dim); % 把标量下界扩展成行向量 end if numel(UB) == 1 UB = repmat(UB, 1, dim); end % ---------- 1. 初始化灰狼种群 ---------- Positions = rand(SearchAgents_no, dim) .* (UB - LB) + LB; Alpha_pos = zeros(1, dim); Beta_pos = zeros(1, dim); Delta_pos = zeros(1, dim); Alpha_score = inf; Beta_score = inf; Delta_score = inf; Convergence_curve = zeros(1, Max_iteration); % ---------- 2. 初始评估 ---------- for i = 1:SearchAgents_no fitness = feval(func, Positions(i, :)); % 评估每只狼 if fitness < Alpha_score Alpha_score = fitness; Alpha_pos = Positions(i, :); elseif fitness < Beta_score Beta_score = fitness; Beta_pos = Positions(i, :); elseif fitness < Delta_score Delta_score = fitness; Delta_pos = Positions(i, :); end end % ---------- 3. 主循环 ---------- for iter = 1:Max_iteration a = 2 - 2 * iter / Max_iteration; % 收敛因子a从2线性衰减到0 for i = 1:SearchAgents_no for j = 1:dim % 对 Alpha 狼的位置更新 r1 = rand(); r2 = rand(); A1 = 2 * a * r1 - a; C1 = 2 * r2; D_alpha = abs(C1 * Alpha_pos(j) - Positions(i, j)); X1 = Alpha_pos(j) - A1 * D_alpha; % 对 Beta 狼的位置更新 r1 = rand(); r2 = rand(); A2 = 2 * a * r1 - a; C2 = 2 * r2; D_beta = abs(C2 * Beta_pos(j) - Positions(i, j)); X2 = Beta_pos(j) - A2 * D_beta; % 对 Delta 狼的位置更新 r1 = rand(); r2 = rand(); A3 = 2 * a * r1 - a; C3 = 2 * r2; D_delta = abs(C3 * Delta_pos(j) - Positions(i, j)); X3 = Delta_pos(j) - A3 * D_delta; % 三头狼共同指导，取平均得到新位置 Positions(i, j) = (X1 + X2 + X3) / 3; end end % ---------- 4. 边界处理 + 重新评估 ---------- for i = 1:SearchAgents_no Positions(i, :) = min(max(Positions(i, :), LB), UB); fitness = feval(func, Positions(i, :)); if fitness < Alpha_score Alpha_score = fitness; Alpha_pos = Positions(i, :); elseif fitness < Beta_score Beta_score = fitness; Beta_pos = Positions(i, :); elseif fitness < Delta_score Delta_score = fitness; Delta_pos = Positions(i, :); end end Convergence_curve(iter) = Alpha_score; end Best_pos = Alpha_pos; Best_score = Alpha_score; end

这个版本没有做排序，因为GWO只需要前三个最优位置，不需要把整个种群按适应度排一遍。用if-elseif结构维护Alpha、Beta、Delta的rank，每次评估一只狼就更新一次，省下不少时间，尤其在高维问题上差别很明显。

2.2 为什么用feval、怎么处理LB/UB

代码里我用的是feval(func, Positions(i,:))，而不是直接写func(Positions(i,:))。原因有两个：第一，feval兼容函数句柄和字符串形式，万一别人传进来的是一串函数名，代码也不会立刻报错；第二，当你需要把GWO封装成一个工具箱接口时，feval的调用方式更统一，方便做函数参数检查。

当然如果你更习惯现代MATLAB风格，把feval(func, x)改成func(x)也完全没问题，效果一样。我保留feval纯粹是兼容性考虑，尤其是在接手师兄师姐老代码的时候，你永远不知道调用方会用哪种姿势传函数。

上下界的处理也是一个小细节。很多测试函数里LB和UB是标量，比如区间[-100,100]，但高维优化时需要把上下界广播到每个维度。我用numel(LB) == 1判断后配合repmat扩展成行向量，这样后面的Positions = rand(...).*(UB-LB)+LB就会自动做逐维运算，避免出现在维度上“串位”的怪问题。如果你用的是R2016b之后的MATLAB，隐式扩展其实也能处理标量和向量的混合运算，但写成repmat之后代码意图更明确，调试时也更清楚。

2.3 调用示例：Sphere与Rastrigin两步跑通

拿到函数后，先用两个经典测试函数验证环境通不通。Sphere函数是最简单的凸函数，Rastrigin则带大量局部极小值，适合观察GWO的跳出能力。

clear; clc; % 第一步：Sphere测试函数，最优解为全0，最优值0 sphere = @(x) sum(x.^2); [Best_score, Best_pos, Convergence_curve] = GWO(sphere, 30, -100, 100, 30, 500); figure; semilogy(Convergence_curve, 'LineWidth', 2); xlabel('迭代次数'); ylabel('最优适应度值（对数坐标）'); title('GWO收敛曲线——Sphere函数'); grid on; fprintf('Sphere最优值: %.6e\n', Best_score); % 第二步：Rastrigin测试函数，同样的最优位置和最优值 rastrigin = @(x) sum(x.^2 - 10*cos(2*pi*x) + 10); [Best_score, Best_pos, Convergence_curve] = GWO(rastrigin, 30, -5.12, 5.12, 40, 1000); figure; semilogy(Convergence_curve, 'LineWidth', 2); xlabel('迭代次数'); ylabel('最优适应度值（对数坐标）'); title('GWO收敛曲线——Rastrigin函数'); grid on; fprintf('Rastrigin最优值: %.6e\n', Best_score);

跑完这两个测试函数，你基本就能判断程序是否正常。Sphere的收敛曲线应该是一条平滑下滑的曲线，Rastrigin的曲线可能会有一段“平台期”，这是算法在跳出局部最优，属于正常现象。如果画semilogy时因为最优值变成0报了警告，可以在curve上加上1e-300这种极小值再画。

3. 从Demo到实际问题：参数调整与工程化改造

3.1 四个关键参数的取舍经验

GWO的参数已经比GA、DE少很多，但“少”不代表可以乱填。我实际使用中重点关注四个数：维度dim、边界LB/UB、狼群规模SearchAgents_no、最大迭代次数Max_iteration。

维度是问题本身决定的，不能改。边界要尽量贴合真实物理约束，太窄会导致最优解落在边界上，太宽会浪费大量迭代去探索无效区域。狼群规模我一般设30到50，少于15匹狼时多样性明显不足，多于100匹狼时计算量增长快但精度提升不明显。最大迭代次数则要看目标函数的计算成本：像Sphere这种秒出结果的函数，跑1000次迭代无所谓；但如果适应度函数是一次仿真，单次评估要好几秒，那就需要控制总评估次数SearchAgents_no * Max_iteration，必要时改成“达到精度就提前终止”的策略。

还有一个容易忽略的参数是收敛因子a。线性衰减是最经典做法，但很多工程问题中前期探索不够，后期又收敛太急。我在项目里经常改成非线性衰减，比如a = 2 * (1 - iter/Max_iteration)^2，或者用余弦曲线，实测对多峰函数有一定改善。想快速出结果时用线性版本，想冲精度时优先试非线性a。

3.2 把测试函数换成自己的目标函数：三个注意点

第一步是把适应度函数抽象成一个“输入一行向量、输出一个标量”的MATLAB函数。最常见的坑有三个：

第一个坑是优化方向。GWO默认做最小化，如果你的业务需求是最大化，不要改算法内部逻辑，直接在适应度函数里返回-f(x)就行，这样改动最小也最安全。第二个坑是目标函数返回维度。它必须返回标量，如果忘了对向量做sum或norm，MATLAB会在赋值时直接报维度不匹配，这个问题后面我会放到报错表里。第三个坑是额外参数传递。如果目标函数需要传入数据、权重矩阵或仿真参数，别改GWO主函数的签名，直接用匿名函数包一层，例如：

result = GWO(@(x) myCostFunction(x, data, weights), dim, lb, ub, 40, 500);

这样GWO函数保持通用，你的业务逻辑全部留在回调函数里，后续换数据集、换模型都只需要重新生成匿名函数，不需要碰算法本体。

3.3 真实工程里GWO常用来干嘛

GWO最舒服的应用场景是那些“求导困难、变量中等规模、评估成本可接受”的优化问题。我之前见过别人用它调PID控制器参数，决策变量就是Kp、Ki、Kd三个参数，适应度函数是阶跃响应的ITAE指标，效果比手动试凑稳定很多。它也被经常用在神经网络训练上，把网络权重和偏置当成一个长向量来优化，不过这个场景计算量较大，一般只用来替代BP的梯度下降做预训练或混合训练。

图像分割阈值选取、特征选择、路径规划、系统辨识与自适应控制参数估计，这些方向GWO都有论文验证过。实际工程里如果已经有降级方案，建议先用GWO做离线寻优，找到一组不错的初始解后再放到在线环境里做局部微调，这样既避免在线优化超时，又能发挥全局搜索的优势。

4. 常见报错与调试技巧实录

4.1 报错速查表

下面这张表是从我这几年使用MATLAB做优化踩过的坑里整理出来的，直接对照现象排查就行。

4.2 不收敛或早熟怎么排查

如果收敛曲线只降了一点就停在原地，先不要急着怀疑GWO本身，通常是边界范围给得不好，或者目标函数存在大片平坦区域。我习惯的做法是先用rng(1)固定随机种子复现问题，然后把初始种群的分布画出来看是否覆盖了解空间。另一种常见情况是a衰减过快，导致狼群还没探索充分就开始收拢，这时候把最大迭代次数适当增加，或者用我前面提过的非线性a方案。

还有一种特别隐蔽的问题来自适应度函数的“非平滑”特性。如果目标函数里有round、if分支或者查表操作，适应度会出现大量相同值，GWO的排名更新就不容易区分优劣。此时可以在适应度函数里加一个很小的正则项，例如x' * x * 1e-8，把相同适应度打散，优化器就能重新获得梯度方向感。

5. 顺手做的几个扩展：混合改进与后续学习路线

5.1 三种低成本改进思路

当你已经把基础GWO跑通，可以试试三种改动成本很低的改进方案。

第一种是上面提到的非线性a衰减，代码只需改一行，但对多峰函数往往有明显提升。第二种是边界反射处理，目前代码用的是截断式min(max(...))，简单但不能保留搜索惯性；改成反射式if x < lb, x = 2*lb - x; end可以防止大量狼堆积在边界上。第三种是和局部搜索混合，每隔固定代数，把当前Alpha位置作为初始点扔给fmincon或patternsearch跑几步，适合连续可导或中等规模的工程问题。

如果你想做二值版本的GWO，比如特征选择这类0/1优化问题，可以在位置更新后加一个Sigmoid转换：x_binary = 1 ./(1 + exp(-x))，再以0.5为阈值转成0或1。这个方法我在两个特征选择项目里用过，效果中规中矩，关键是代码复用率很高。

5.2 踩坑总结：给后来者的几条实在建议

最后再说几句掏心窝子的经验。在把GWO接到自己的工程之前，建议先把目标函数写成独立脚本，并用一组“已知最优解”验证返回值和手算结果一致。这一步看起来多余，实际能省掉后面所有Debug时间——我见过有人把目标函数写成返回向量，结果一路查到收敛曲线的绘制才发现问题，浪费了整整一个下午。

另外，跑对比实验时别只运行一次就下结论。GWO和大多数元启发式算法一样，单次结果受随机数影响很大，写论文或做技术报告时至少跑20到30次独立试验，记录均值、标准差和最好值，画Box-plot比只贴一条收敛曲线有说服力得多。如果时间有限，至少也要运行5次看趋势是否稳定，否则你很难区分“算法牛”还是“运气好”。按照这个流程走一遍，我猜你很快就能把GWO从Demo阶段搬到自己的任务里，踩过的坑和调好的参数就是你自己的方法论。

本文还有配套的精品资源，点击获取