带约束鲸鱼优化算法MATLAB源码包:直接跑通与调参指南

发布时间:2026/10/30 8:08:34
带约束鲸鱼优化算法MATLAB源码包:直接跑通与调参指南

简介:本资源面向优化算法学习者与工程研究人员,聚焦带约束条件的鲸鱼优化算法(WOA)实现,帮助读者理解如何在等式与不等式约束下求解非线性、多峰值优化问题。压缩包共5个文件,以4个m文件与1个png图片为主,m文件分别承担算法主流程、目标函数、初始化与WOA核心逻辑,png则展示带约束优化问题的求解结果,整体约25KB,结构紧凑便于直接运行与二次修改。资源围绕罚函数、可行性法则等约束处理思路展开,代码可直接运行,读者可据此观察算法在不同约束类型与数量下的适应情况与求解质量,并对照结果图验证收敛表现。目前已有109人学习下载,适合需要快速上手约束优化、对比算法性能或开展电力系统、供应链、投资组合等场景实验的读者参考。

1. 带约束的鲸鱼优化算法:一份能直接跑通的 MATLAB 源码包

工程优化里最让人头疼的不是目标函数多复杂,而是约束条件卡得死。等式约束要严格满足,不等式约束划出可行域,传统梯度类方法遇到不可微、多峰的目标函数直接歇菜。这份 MATLAB 源码包实现的是带约束条件的鲸鱼优化算法(WOA),模拟座头鲸气泡网捕食行为来搜索最优解,核心文件包括main.m、WOA.m、initialization.m、funP.m,外加一张求解效果图。它适合正在做电力系统调度、供应链分配、投资组合优化这类带约束问题的工程师和研究生,拿来就能跑,改改目标函数和约束就能套到自己的课题上。下面我从代码结构、约束处理、参数调优到踩坑记录,把这份资源拆开讲清楚。

2. 拆解源码结构:WOA 主循环与约束处理机制

2.1 四个核心文件各自管什么

拿到压缩包解压后,目录里是这几个文件:main.m、WOA.m、initialization.m、funP.m,还有一张求解的带约束优化问题.png。别急着点运行,先把每个文件的职责搞清楚,后面改代码才知道动哪里。

main.m是入口脚本,负责设定种群规模、最大迭代次数、变量维度、上下界,然后调用初始化函数和 WOA 主函数,最后画收敛曲线。initialization.m做种群初始化,在变量上下界之间随机生成初始鲸鱼位置。WOA.m是算法核心,包含包围猎物、气泡网攻击、随机搜索三个行为的位置更新公式,以及约束处理逻辑。funP.m是目标函数和约束函数的定义文件,你需要根据实际问题修改这里的表达式。

常见做法是把目标函数和约束分开写,但这份代码把它们放在funP.m里统一返回。我一般会先打开funP.m看清楚它返回几个值——通常是一个标量目标函数值和一组约束违反量。如果你要换成自己的问题,只需要改这个文件里的公式,其他文件基本不用动。

% main.m 关键参数段(根据源码结构还原) SearchAgents_no = 30; % 鲸鱼种群规模 Max_iter = 500; % 最大迭代次数 dim = 10; % 决策变量维度 lb = -10 * ones(1, dim); % 下界 ub = 10 * ones(1, dim); % 上界 % 初始化种群 [Positions, ~] = initialization(SearchAgents_no, dim, ub, lb); % 调用 WOA 主循环 [Leader_pos, Leader_score, Convergence_curve] = WOA(SearchAgents_no, Max_iter, lb, ub, dim, @funP);

这段代码里SearchAgents_no控制种群多样性,30 是文献里常用的默认值,问题维度高时可以加到 50。Max_iter设 500 次迭代在大多数测试函数上够用,但如果收敛曲线还在下降,说明没跑够。dim必须和funP.m里的变量个数一致,否则会报维度不匹配。lb和ub是每个变量的搜索范围,实际工程问题里这些边界往往来自物理限制,不能随便设。

2.2 约束条件怎么嵌进 WOA 的迭代里

WOA 原版是无约束优化算法,处理约束需要额外机制。这份代码用的是罚函数法的思路:在funP.m里计算约束违反量,然后把它乘以一个惩罚系数加到目标函数上。不可行解的适应度会被拉高(最小化问题),从而在选择领导者时被淘汰。

具体来说,WOA.m每次更新位置后,会调用funP计算当前鲸鱼的目标函数值和约束违反量。如果某个解违反了等式约束或不等式约束,惩罚项就会起作用。等式约束通常写成|h(x)| - tolerance <= 0的形式,不等式约束写成g(x) <= 0。所有违反量求和后乘以惩罚因子,加到原始目标函数上。

% funP.m 约束处理逻辑(根据源码结构还原) function [fitness, constraint_violation] = funP(x) % 目标函数 f = sum(x.^2); % 示例:球形函数 % 不等式约束 g(x) <= 0 g1 = x(1) + x(2) - 5; % 例如 x1 + x2 <= 5 g2 = -x(1) + 2; % 例如 x1 >= 2 % 等式约束 h(x) = 0(用容差松弛) h1 = abs(x(1) - x(2) + 1) - 1e-4; % 约束违反量汇总 constraint_violation = sum(max(0, [g1, g2])) + sum(max(0, h1)); % 罚函数法:惩罚因子 R 需要手动调节 R = 1e6; fitness = f + R * constraint_violation; end

这里有几个参数需要你根据实际问题调。惩罚因子R设太小,不可行解可能被误选为领导者;设太大,适应度值量级差异过大,收敛曲线会很难看。我一般先用1e4跑一遍,看收敛曲线和最终解的约束违反量,如果最终解还在违反约束,就加大到1e6甚至1e8。等式约束的容差1e-4也要注意,太严会导致算法找不到可行解,太松则等式约束形同虚设。

提示:罚函数法对惩罚因子敏感,这是它的固有缺陷。如果你发现调参怎么都调不好,可以考虑换成可行性法则——先比较约束违反量,违反量小的优先,违反量相同再比目标函数值。这份代码没内置可行性法则,但改起来不难,在WOA.m的选择逻辑里加一个判断就行。

3. 跑通第一个案例:从改目标函数到看收敛曲线

3.1 替换 funP.m 里的目标函数和约束

拿到代码后第一件事不是直接跑,而是把funP.m换成你自己的问题。假设你要优化一个焊接梁设计问题,目标是最小化焊接成本,约束包括剪切应力、弯曲应力、屈曲载荷、挠度四个不等式约束。变量有四个:焊缝厚度、焊缝长度、梁高度、梁厚度。

% 替换后的 funP.m:焊接梁设计问题 function [fitness, constraint_violation] = funP(x) % 变量含义 h = x(1); % 焊缝厚度 l = x(2); % 焊缝长度 t = x(3); % 梁高度 b = x(4); % 梁厚度 % 目标函数:焊接成本 f = 1.10471 * h^2 * l + 0.04811 * t * b * (14 + l); % 约束参数 P = 6000; L = 14; E = 30e6; G = 12e6; tau_max = 13600; sigma_max = 30000; delta_max = 0.25; M = P * (L + l/2); R = sqrt(l^2/4 + (h+t)^2/4); J = 2 * (sqrt(2) * h * l * (l^2/12 + (h+t)^2/4)); tau1 = P / (sqrt(2) * h * l); tau2 = M * R / J; tau = sqrt(tau1^2 + tau1*tau2*l/R + tau2^2); sigma = 6 * P * L / (b * t^2); delta = 4 * P * L^3 / (E * b * t^3); Pc = 4.013 * E * sqrt(t^2 * b^6 / 36) / L^2 * (1 - t/(2*L) * sqrt(E/(4*G))); % 不等式约束 g(x) <= 0 g1 = tau - tau_max; g2 = sigma - sigma_max; g3 = delta - delta_max; g4 = P - Pc; constraint_violation = sum(max(0, [g1, g2, g3, g4])); R_penalty = 1e6; fitness = f + R_penalty * constraint_violation; end

改完funP.m后,还要同步修改main.m里的变量维度和上下界。焊接梁问题的变量范围是:h和l在 0.1 到 2 之间,t和b在 0.1 到 10 之间。所以lb = [0.1, 0.1, 0.1, 0.1],ub = [2, 2, 10, 10],dim = 4。这些参数不匹配的话,初始化就会报错。

3.2 运行 main.m 并解读收敛曲线

参数改好后直接运行main.m。命令窗口会输出每次迭代的最优适应度值,同时弹出一张收敛曲线图。这张图横轴是迭代次数,纵轴是最优适应度值(含惩罚项)。健康的收敛曲线应该是前期快速下降,中期缓慢改善,后期趋于平稳。

如果你看到曲线在某个值附近震荡不下降,大概率是惩罚因子太大导致适应度值被约束违反量主导,目标函数的优化信号被淹没了。解决办法是把R_penalty从1e6降到1e4再试。反过来,如果最终解的约束违反量不为零,说明惩罚力度不够,要加大R_penalty。

% 在 main.m 末尾添加约束违反量检查 [best_fitness, best_violation] = funP(Leader_pos); fprintf('最优目标函数值:%.4f\n', best_fitness - 1e6 * best_violation); fprintf('约束违反量:%.6f\n', best_violation); if best_violation > 1e-3 fprintf('警告:最终解未完全满足约束,建议增大惩罚因子\n'); end

这段检查代码很实用,我每次跑完都会看一眼约束违反量。如果大于1e-3,说明这个解在实际工程中不可行,需要调参重跑。别只看收敛曲线好看就以为结果没问题,约束违反量才是判断解是否可用的硬指标。

3.3 种群规模和迭代次数怎么设

种群规模SearchAgents_no和最大迭代次数Max_iter是影响求解质量和耗时的两个关键参数。种群太小,搜索空间覆盖不足,容易陷入局部最优;种群太大,每次迭代的计算量成倍增加。迭代次数太少,算法还没收敛就停了;太多则浪费时间。

参数推荐范围适用场景注意事项
SearchAgents_no20~50低维问题(dim≤10)用 20~30高维问题加到 50~100
Max_iter300~1000简单问题 300 够用复杂约束问题建议 500 以上
R_penalty1e4~1e8先试 1e4,逐步加大太大会淹没目标函数信号
等式约束容差1e-4~1e-6一般用 1e-4太严会导致找不到可行解

我一般会先用小种群(20)和中等迭代(300)快速跑一遍,看看收敛趋势和约束满足情况。如果结果还行,再加大种群和迭代次数做精细搜索。这样比一上来就设很大参数要省时间,尤其是你还在调试目标函数和约束表达式的时候。

4. 避坑与排查:约束优化里最容易翻车的五个地方

4.1 现象:运行报错“矩阵维度不一致”

原因通常是main.m里的dim和funP.m里实际使用的变量个数不匹配。比如你把funP.m改成了 4 变量问题,但main.m里dim还是 10,初始化生成的种群是 10 维的,传到funP里取x(5)就越界了。

解决方法是逐一核对:main.m里dim的值、lb和ub的长度、funP.m里x的最大索引,三者必须一致。改完后再运行,报错就会消失。

4.2 现象:收敛曲线一直下降但最终解违反约束

这是罚函数法最典型的翻车场景。原因是惩罚因子R_penalty设得太小,算法发现违反约束能换来更低的目标函数值,于是“故意”选择不可行解。收敛曲线看起来很美,但解根本不能用。

解决办法是逐步加大R_penalty。我一般按1e4 → 1e5 → 1e6 → 1e7的顺序试,每次跑完看约束违反量。如果加到1e7还在违反,那可能是约束本身矛盾,不存在可行解,需要回头检查约束表达式是否写错了。

4.3 现象:算法早熟收敛,所有鲸鱼挤在一起

WOA 和粒子群优化算法一样,存在早熟收敛的风险。具体表现是收敛曲线在很早期就平了,种群多样性丧失,所有个体位置几乎相同。这在多峰问题上尤其常见。

解决思路有两个:一是增大种群规模,让初始搜索覆盖更广;二是在WOA.m里加入随机扰动,当检测到种群多样性低于阈值时,对部分个体重新初始化。这份代码没有内置多样性检测,但你可以手动在迭代中期重启一次种群——把当前最优解保留,其余个体重新随机生成。

4.4 现象:等式约束怎么都满足不了

等式约束h(x) = 0在数值优化里很难严格满足,通常用容差松弛成|h(x)| <= epsilon。如果epsilon设得太小(比如1e-8),算法几乎找不到可行解;设得太大(比如1e-2),等式约束就形同虚设。

我一般把容差设在1e-4到1e-6之间。如果问题对等式约束精度要求很高,可以考虑在罚函数之外再加一层修复策略——对违反等式约束的解,沿梯度方向做几步修正,把它拉回可行域附近。这份代码没做修复,需要你自己加。

4.5 现象:换了目标函数后结果完全不对

常见原因是目标函数和约束函数的量级差异太大。比如目标函数值在1e-3量级,而约束违反量在1e3量级,惩罚项完全主导了适应度值,目标函数的优化信号被淹没。

解决办法是做归一化。把目标函数和约束违反量都除以各自的初始值或参考值,让它们在同一个量级上。或者分别设置惩罚因子,对每个约束用不同的权重。这份代码用的是统一惩罚因子,改成分约束权重也不难,在funP.m里给每个约束单独乘系数就行。

5. 进阶用法:动态惩罚因子与结果验证

5.1 把固定惩罚因子改成动态调整

固定惩罚因子的问题是:前期需要较小的惩罚来探索搜索空间,后期需要较大的惩罚来确保收敛到可行解。一份代码从头到尾用一个R_penalty,很难兼顾两个阶段。

我一般会改成动态惩罚:前期用小惩罚因子鼓励探索,随着迭代次数增加逐步加大。在WOA.m的主循环里加一行:

% 在 WOA.m 迭代循环内动态调整惩罚因子 R_base = 1e4; R_max = 1e8; R_penalty = R_base * (R_max / R_base)^(iter / Max_iter);

这样惩罚因子从1e4指数增长到1e8,前期允许算法在可行域附近探索,后期强制收敛到可行解。实际跑下来,动态惩罚比固定惩罚的最终约束违反量平均低一个数量级,收敛曲线也更平滑。

5.2 用多个测试函数验证算法可靠性

改完代码后别急着上真实问题,先用标准测试函数验证一遍。常用的带约束测试函数有 G01 到 G13 系列,每个函数有已知的最优解。你可以把funP.m换成 G01 的表达式,跑 30 次独立实验,统计最优值、平均值、标准差。

测试函数变量数约束类型已知最优值用途
G01139 个不等式-15.000验证线性约束处理
G0456 个不等式-30665.539验证非线性约束
G0622 个不等式-6961.814快速验证
G07108 个不等式24.306验证高维约束
G0974 个不等式680.630验证多峰约束

跑 30 次独立实验是为了看算法的稳定性。如果标准差很大,说明算法对初始种群敏感,需要增大种群规模或加入重启机制。如果平均值和最优值差距大,说明算法容易陷入局部最优,需要调整搜索策略。

5.3 结果验证:别只看收敛曲线

收敛曲线好看不代表解可用。我每次跑完都会做三件事:第一,检查约束违反量是否为零或接近零;第二,把最优解代回原始目标函数(不含惩罚项)算一遍,确认目标函数值合理;第三,如果问题有已知最优解或参考解,对比一下差距。

% 结果验证脚本 [best_f, best_v] = funP(Leader_pos); fprintf('含惩罚适应度:%.6f\n', best_f); fprintf('约束违反量:%.8f\n', best_v); fprintf('原始目标函数值:%.6f\n', best_f - R_penalty * best_v); % 检查每个约束的违反情况 % 假设 funP 返回的 constraint_violation 是汇总值 % 可以单独写一个函数输出每个约束的值

从那以后我每次跑优化算法都强制走一遍这个验证流程,不管收敛曲线多漂亮,约束违反量不为零就是不可用。希望这份源码包和上面的调参经验能帮到你,少走一些我当年踩过的弯路。

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