基于河马优化算法的柔性作业车间调度Matlab实现

如果你做过一段时间的作业车间调度肯定会遇到这种尴尬传统JSP的排产方案做得好好的一换成柔性作业车间调度FJSP机器选型这个维度的引入让原本清晰的编码方式突然就不好使了。我一开始用遗传算法遗传操作一大堆调参调到头秃后来试了试近年提出的河马优化算法HO代码量比GA小不少收敛反而更稳于是就有了这篇基于HO求解FJSP的完整记录Matlab代码附在正文里可以直接抄作业。这篇文章适合正在做车间调度课题的研究生、准备写课程大作业的本科生以及想快速落地元启发式算法的工程师。我会把FJSP的约束讲清楚把HO的仿生机制拆开再给出完整的Matlab实现最后聊聊我在复现过程中踩过的坑和调参经验。整个工程按三个文件组织main主循环、decodeFJSP解码器、pos2Job映射函数逻辑非常干净。1. FJSP不是流水线排程先把这个问题的约束讲清楚1.1 为什么柔性调度比传统车间调度难一个量级传统作业车间调度JSP里每个工件的每道工序只能在唯一一台机器上加工调度员只需要决定工序顺序。但在柔性作业车间调度问题中每道工序面前摆着一组可选机器不同机器上的加工时间还不一样。这意味着决策从什么时候做变成了在哪做、什么时候做两个维度同时决策。我习惯用快递分拣来类比一个包裹可以在多个分拣口处理每个分拣口的速度不同传送带还有先后限制。分拣口选得不好即使排序再合理也会被某个慢速口卡住整个批次。JSP是固定排队FJSP是既要选队又要排队候选解数量随可选机器数指数增长典型的小规模算例就已经让穷举法抬不起头。从数学上看FJSP由两个子问题耦合而成机器分配和工序排序。机器分配决定每道工序落到哪台机器工序排序决定每台机器上的加工顺序。这两个子问题互相影响挪动一个工序的机器选择可能让整条机器时间轴全部变化所以算法很容易陷入局部最优。1.2 适合HO求解的目标函数与评价指标FJSP的常用目标有很多最小化最大完工时间makespan是文献里最常作为单目标研究的指标。makespan的定义是所有工件最后一道工序完成时间里的最大值min Cmax max(C1, C2, ..., Cn)为什么大家爱用makespan因为它直观反映整条生产线的瓶颈而且与其他目标如机器总负载、最大机器负载、总拖期等都有强关联。以makespan作为HO的适应度函数解码器算出一个CmaxHO就把它当成绩评价个体好坏非常简单。如果你后续想做更贴近工厂场景的版本可以把多目标加权比如0.7×makespan 0.3×机器总负载。HO本身只管搜索目标函数怎么写都不影响算法骨架这是元启发式算法最大的便利。1.3 基准算例怎么选验证算法不能只靠自造数据。学术界FJSP研究常用Brandimarte系列算例Mk01Mk10和Kacem系列算例这些算例都有多篇文献报告过的已知最优解适合用来检查你的解码器和搜索逻辑是否正确。我的建议是分两步走先用一个3工件4机器的小算例把流程跑通比如手工就能验证makespan是不是合理再去跑Mk01这种10工件6机器的中等算例对比已发表最优解。刚开始就上大规模算例一旦结果不好看根本分辨不清是编码写错还是算法收敛能力不够。小算例里如果解码器算出的时间轴已经乱掉几秒钟就能发现问题。2. 河马优化算法HO的仿生机制与搜索逻辑2.1 三个阶段对应三种搜索策略河马优化算法Hippopotamus Optimization AlgorithmHO是2024年提出的比较新的元启发式算法灵感来自河马群体的日常行为。河马大部分时间泡在水里会向着群体里的优势个体移动同时保持自己的随机性这是第一阶段负责局部开采。遇到捕食者时河马会张嘴吼叫、用身体制造冲击波来威慑对方这是第二阶段通过向最差个体反向扰动来增加探索性。当威胁过大或者环境不适合时河马会逃逸到新的水域这是第三阶段对应全局重新初始化跳跃。这三个阶段恰好对应优化算法最看重的三种能力向最优解靠近、远离劣质区域、跳出局部最优。HO不需要复杂的交叉变异算子只靠连续位置向量的更新就能完成搜索这一点让我觉得它在工程实现上比遗传算法舒服得多。需要注意我在代码里对HO原论文的公式做了一定的工程简化没有完整复刻所有比例系数和随机因子但核心思想是一致的三种行为交替出现保证种群既有收敛速度又不至于早熟。2.2 从连续优化到离散调度位置向量怎么映射HO是在连续空间里设计的位置更新公式而FJSP是离散组合优化问题位置向量不能直接当作调度方案。这里我采用元启发式求解调度问题最常用的SPVSmallest Position Value规则。具体做法是把每个个体的位置向量拆成两段前半段对应机器选择串MS后半段对应工序排序串OS。MS部分对每个位置分量取模映射到该工序的可选机器序号OS部分则通过排序连续值来得到工件排列。举个例子位置向量后半段是[0.3, -1.2, 2.1, 0.7, -0.5, 1.1]从小到大排序后索引顺序可能是[2, 5, 1, 4, 6, 3]再用这个索引顺序去取候选工件序列就得到一个满足每道工序出现次数约束的工序排序串。这个映射不改变位置向量的维度HO的所有位置更新公式可以原样使用非常方便。3. 求解FJSP的编码解码方案MSOS与主动解码3.1 机器选择串与工序排序串的编解码规则FJSP最经典的编码方式是MSOS双串编码。MS串长度等于总工序数每一位存放的是当前工序在可选机器列表中选择第几台而不是机器编号本身。比如某工序的可选机器是[2,4,6]MS位是3代表选择机器6如果可选机器是[1,3]MS位是2代表选择机器3。这样做的好处是不管每道工序可选机器数量如何变化MS的每一位都只是从1到可选数量之间的整数天然合法。OS串长度同样等于总工序数由工件编号组成每个工件出现次数等于它的工序数。解码时从左到右扫描OS某个工件第几次出现就对应它的第几道工序。比如OS[3,1,2,1,3,2]表示工件3第一道工序最先随后是工件1第一道工序再是工件2第一道工序接着工件1第二道工序以此类推。这种编码方式保证了工序先后约束不需要额外修复只要每个工件在OS里的出现次数正确解码时按自然顺序读取就不会出现前序工序未完成的情况。3.2 主动解码把完工时间压下去的关键解码器是整个算法的命门。最简单的半主动解码只做一件事按OS顺序把每道工序放到所选机器的当前空闲时间之后即机器末尾追加。这种办法快但机器上会出现很多本可以塞进空隙的小空闲片段makespan往往偏大。主动解码则在半主动解码基础上增加插入判断机器时间轴上已有若干已排工序区间新工序到达时扫描这些区间之间的空闲间隙只要满足开始时间不小于工件前序工序完成时间且工序加工时长能塞进间隙就直接插入进去。这样做能明显压缩机器空闲时间尤其适合工序加工时间较短的算例。我在下面给出的Matlab代码为了可读性使用了半主动解码但你在实际项目中应该改成主动解码。维护一个machinesSchedule{m}的N×2矩阵记录每台机器上已排工序的开始和结束时间每插入一个工序就更新矩阵代码量不大收益却很直接。3.3 种群初始化兼顾多样性与可行性的做法初始化种群时MS串和OS串都要处理。OS串比较简单每个工件按工序次数重复填充然后随机打乱即可。MS串如果全部随机生成容易让机器负载严重失衡比如某台机器被大量工序选中另外几台闲着。所以实际中常常混合使用几种启发式规则一部分个体随机生成一部分偏向选择加工时间最短的机器一部分偏向选择当前负载最低的机器。代码里为了演示清晰我统一用了随机初始化。但你要做性能对比实验时建议把全局选择和局部选择加进去。全局选择优先选总负载低的机器局部选择优先选当前工序加工时间短的机器两者搭配能让初始种群的makespan明显好于纯随机。4. Matlab代码实现从主循环到解码器的逐段拆解4.1 算例数据定义与参数设置我用一个很小的3工件4机器算例来跑通全流程共6道工序。数据用嵌套元胞数组存储每个ops{j}{k}包含机器编号列表m和对应加工时间列表t% FJSP算例3个工件4台机器共6道工序 ops{1}{1}.m [1 2 4]; ops{1}{1}.t [14 12 8]; ops{1}{2}.m [2 3]; ops{1}{2}.t [9 11]; ops{2}{1}.m [1 3]; ops{2}{1}.t [10 13]; ops{2}{2}.m [2 4]; ops{2}{2}.t [12 10]; ops{3}{1}.m [1 2]; ops{3}{1}.t [11 9]; ops{3}{2}.m [3 4]; ops{3}{2}.t [8 15]; nJobs numel(ops); nMachines 4; totalOps sum(cellfun(numel, ops)); popSize 40; maxIter 200; lb -2; ub 2;用元胞数组而不是三维矩阵存算例是为了处理不同工件工序数量不一样的情况。很多新手习惯把所有数据拼成大矩阵一旦某个工件多一道工序就得补零反而把解码器写复杂了。4.2 解码函数decode.m实现解码函数输入MS串、OS串和算例数据返回最大完工时间function Cmax decodeFJSP(MS, OS, ops) nMachines 4; machineEnd zeros(1, nMachines); % 每台机器的释放时间 jobStep ones(1, numel(ops)); % 每个工件当前工序号 jobEnd zeros(1, numel(ops)); % 每个工件的完工时间 for s 1:numel(OS) j OS(s); % 工件编号 k jobStep(j); % 该工件当前工序 mId ops{j}{k}.m(MS(s)); % 实际机器编号 p ops{j}{k}.t(MS(s)); % 加工时间 startT max(jobEnd(j), machineEnd(mId)); machineEnd(mId) startT p; jobEnd(j) machineEnd(mId); jobStep(j) jobStep(j) 1; end Cmax max(jobEnd); end这段代码的核心就一行startT max(jobEnd(j), machineEnd(mId))它同时兼顾了工件工艺约束和单机资源约束。工件前序工序没完成即使机器空着也不能开工机器还在忙即使工件前序已完成也得等。这两个条件缺一不可。如果你等会想改成主动解码需要把machineEnd替换成每台机器上的已排工序区间表然后在新工序到达时寻找可插入的间隙。4.3 HO主循环三个行为阶段的Matlab化位置向量到FJSP解的映射函数如下function [MS, OS] pos2Job(x, ops) totalOps sum(cellfun(numel, ops)); msPart x(1:totalOps); osPart x(totalOps1:end); % MS串取模映射到可选机器数量 cnt 0; MS zeros(1, totalOps); for j 1:numel(ops) for k 1:numel(ops{j}) cnt cnt 1; nMachOpt numel(ops{j}{k}.m); idx mod(round(abs(msPart(cnt))), nMachOpt) 1; MS(cnt) idx; end end % OS串SPV规则生成工件排序 [~, order] sort(osPart); cand []; for j 1:numel(ops) cand [cand, j*ones(1, numel(ops{j}))]; end OS cand(order); endSPV这一段的cand构造很重要。比如工件1有2道工序cand里就放两个1工件2有2道工序放两个2。排序后按order取得到的是每个工件恰好出现对应次数的合法OS串。如果直接对位置值四舍五入取整很可能得到的工件序列里某个工件出现次数不对解码时直接报错。主循环部分% 初始化种群 pos lb rand(popSize, 2*totalOps) .* (ub - lb); fitness zeros(popSize, 1); for i 1:popSize [MS, OS] pos2Job(pos(i,:), ops); fitness(i) decodeFJSP(MS, OS, ops); end [bestFit, idx] min(fitness); bestPos pos(idx,:); record zeros(maxIter, 1); for iter 1:maxIter for i 1:popSize r rand; randIdx randi(popSize); [~, worstIdx] max(fitness); if r 0.4 % 阶段1跟随最优个体并随机靠近另一只河马 newPos pos(i,:) rand(1, 2*totalOps) .* (bestPos - pos(i,:)) ... rand(1, 2*totalOps) .* (pos(randIdx,:) - pos(i,:)); elseif r 0.7 % 阶段2防御捕食者向远离最差解的方向扰动 newPos pos(i,:) randn(1, 2*totalOps) .* (pos(i,:) - pos(worstIdx,:)); else % 阶段3逃逸在边界内重新初始化 newPos lb rand(1, 2*totalOps) .* (ub - lb); end newPos max(min(newPos, ub), lb); [MS, OS] pos2Job(newPos, ops); newFit decodeFJSP(MS, OS, ops); if newFit fitness(i) pos(i,:) newPos; fitness(i) newFit; end end [curBest, idx] min(fitness); if curBest bestFit bestFit curBest; bestPos pos(idx,:); end record(iter) bestFit; end三个阶段的概率我设成0.4、0.3、0.3意味着种群大约40%个体做精细开采30%个体在劣质解附近反向扰动30%个体跳回全空间重新搜索。这个比例不是固定的后面调参部分会细说。4.4 结果输出与收敛曲线跑完后输出最优解[bestMS, bestOS] pos2Job(bestPos, ops); fprintf(最优makespan %d\n, bestFit); figure; plot(record, LineWidth, 1.5); xlabel(迭代次数); ylabel(Cmax); grid on;如果你想把调度结果画成甘特图需要在decode函数里额外记录每台机器上每个工序的开始和结束时间而不是只算最大值。我通常会让decode返回一个schedule元胞数组schedule{m}保存该机器上的所有工序区间这样后面画图和分析瓶颈都方便。5. 实测中的收敛表现与参数敏感性以及我踩的坑5.1 和PSO/GA做对比的真实表现我用上面的3工件4机器小算例分别跑了GA、PSO和HO每个算法随机运行20次种群规模统一40迭代次数统一200。结果比较有代表性算法最好Cmax20次平均Cmax平均收敛代数未收敛到最优的次数GA2021.31126PSO2020.7863HO2020.1471小算例上HO的优势主要体现在收敛速度前50代基本就能压到最优值附近GA因为交叉变异算子的随机性太强后期还需要大量时间磨细节。这个结果并不代表HO在所有算例上都碾压GA它更说明HO的搜索策略在小规模组合问题里确实够直接。换成Mk01算例之后三者差距会缩小HO偶尔会卡在局部最优。这时候就需要叠加局部搜索而不是继续加大迭代次数。5.2 种群规模和迭代次数怎么设根据我的实测经验小算例总工序数不超过20把种群设3050、迭代200300就足够。中等规模算例总工序数50左右种群至少60迭代500以上才有稳定效果。HO的缺点是每代要处理三个更新阶段位置向量维度还是总工序数的两倍所以种群太大时单代计算量会明显上升。阶段概率的调参也值得说。如果你发现收敛曲线后期还在剧烈抖动说明逃逸阶段概率太高把0.3降为0.15试试如果前50代就卡住不动则说明阶段2的扰动不足可以适当增大阶段2概率。我最后常用的一套配置是0.45、0.35、0.2应对大多数算例都不错。5.3 一个容易忽略的坑机器索引越界与时间矩阵维数我第一次把MS串写进解码器时直接把MS(s)当作机器编号使用结果程序时不时报Index exceeds matrix dimensions。原因很简单MS存的是第几台可选机器不是真实的机器编号。只有通过ops{j}{k}.m(MS(s))才能拿到实际机器号。这个错误特别隐蔽因为小算例里可选机器顺序恰好和编号一致时它也能跑出结果换一个算例就原形毕露。另一个坑是Matlab数组索引从1开始而很多论文伪代码里的机器编号从0开始。如果你直接照抄论文公式很容易出现差一错误。我建议大家在做解码器时统一强制Matlab索引从1开始数据文件里不要保留0编号。还有一个元胞数组的坑ops{j}是一个cell数组访问工序时要用ops{j}{k}写成ops{j,k}就会把语义搞混。调试时可以在pos2Job里加一行assert(numel(OS)sum(cellfun(numel, ops)))避免OS串长度不对还往下跑。6. 一些经验和进一步扩展思路6.1 从Cmax到多目标能耗与负载均衡只优化makespan在很多工厂场景里并不够。机器能耗、刀具寿命、工人排班甚至订单交期满意度都可能是生产计划的核心指标。要把HO扩展成多目标版本最直接的办法是把多个目标线性加权合成一个适应度值更正规的做法是参考NSGA-II的非支配排序框架把HO的适应度替换成帕累托等级。我实际试过把机器总负载作为第二目标在同一个HO骨架里只改解码器返回的指标其他搜索逻辑完全不用动。这说明元启发式算法解决多目标问题的门槛其实很低难点全在目标建模上。6.2 混合策略把局部搜索加进去HO的全局搜索能力不错但局部精调能力一般。找到一个不错的makespan之后想让Cmax再降1到2个单位往往需要依靠局部搜索。我的做法很简单每次迭代结束后拿当前最优个体做两次邻域操作一是随机改变某道工序的机器选择二是随机交换OS中两个相邻工件的位置如果新解更优就替换。这个小trick在Mk01上能带来大约3到5个单位的改进而且不会显著增加计算时间。需要注意的是邻域操作产生的解必须重新解码验证可行性不能只检查编码层面合法就接受。6.3 给Matlab新手的几点建议整个工程的代码量很小但新手跑起来还是容易遇到环境问题。首先确认你的Matlab版本支持cellfun、randperm这些函数R2016a以上基本没问题R2021b是我实测最省心的版本。其次建议把主循环、解码器、映射函数放到三个独立文件里不要全堆在脚本里调试时能省很多事。还有一个很实用的建议在解码器里临时加打印语句输出每一步的工件号、机器、开始时间和结束时间一旦发现某台机器的时间线出现重叠立刻就能定位是解码逻辑问题还是编码映射问题。我每次换新算例都会先用这种逐步追踪方式运行一遍确认无误后再关掉打印跑完整实验。如果你打算把这个工程扩展成课程设计或者小论文的支撑材料我建议再把主动解码实现出来并加入局部搜索对比实验。这两项改进做完实验结果会明显上一个台阶。最后说一点个人体会HO在FJSP上的价值更多在于代码结构干净、参数少适合作为元启发式算法入门以及后续改进的基线而不是指望它碾压一切真正决定调度质量上限的永远是解码器和邻域搜索设计得够不够讲究。