基于河马优化算法的柔性作业车间调度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的适应度函数，解码器算出一个Cmax，HO就把它当成绩评价个体好坏，非常简单。

如果你后续想做更贴近工厂场景的版本，可以把多目标加权，比如"0.7×makespan + 0.3×机器总负载"。HO本身只管搜索，目标函数怎么写都不影响算法骨架，这是元启发式算法最大的便利。

1.3 基准算例怎么选

验证算法不能只靠自造数据。学术界FJSP研究常用Brandimarte系列算例（Mk01～Mk10）和Kacem系列算例，这些算例都有多篇文献报告过的已知最优解，适合用来检查你的解码器和搜索逻辑是否正确。

我的建议是分两步走：先用一个3工件4机器的小算例把流程跑通，比如手工就能验证makespan是不是合理；再去跑Mk01这种10工件6机器的中等算例，对比已发表最优解。刚开始就上大规模算例，一旦结果不好看，根本分辨不清是编码写错还是算法收敛能力不够。小算例里如果解码器算出的时间轴已经乱掉，几秒钟就能发现问题。

2. 河马优化算法HO的仿生机制与搜索逻辑

2.1 三个阶段对应三种搜索策略

河马优化算法（Hippopotamus Optimization Algorithm，HO）是2024年提出的比较新的元启发式算法，灵感来自河马群体的日常行为。河马大部分时间泡在水里，会向着群体里的优势个体移动，同时保持自己的随机性，这是第一阶段，负责局部开采。遇到捕食者时，河马会张嘴吼叫、用身体制造冲击波来威慑对方，这是第二阶段，通过向最差个体反向扰动来增加探索性。当威胁过大或者环境不适合时，河马会逃逸到新的水域，这是第三阶段，对应全局重新初始化跳跃。

这三个阶段恰好对应优化算法最看重的三种能力：向最优解靠近、远离劣质区域、跳出局部最优。HO不需要复杂的交叉变异算子，只靠连续位置向量的更新就能完成搜索，这一点让我觉得它在工程实现上比遗传算法舒服得多。

需要注意，我在代码里对HO原论文的公式做了一定的工程简化，没有完整复刻所有比例系数和随机因子，但核心思想是一致的：三种行为交替出现，保证种群既有收敛速度又不至于早熟。

2.2 从连续优化到离散调度：位置向量怎么映射

HO是在连续空间里设计的位置更新公式，而FJSP是离散组合优化问题，位置向量不能直接当作调度方案。这里我采用元启发式求解调度问题最常用的SPV（Smallest 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(totalOps+1: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); end

SPV这一段的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。结果比较有代表性：

小算例上HO的优势主要体现在收敛速度，前50代基本就能压到最优值附近；GA因为交叉变异算子的随机性太强，后期还需要大量时间磨细节。这个结果并不代表HO在所有算例上都碾压GA，它更说明HO的搜索策略在小规模组合问题里确实够直接。

换成Mk01算例之后，三者差距会缩小，HO偶尔会卡在局部最优。这时候就需要叠加局部搜索，而不是继续加大迭代次数。

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

根据我的实测经验，小算例（总工序数不超过20）把种群设30～50、迭代200～300就足够。中等规模算例（总工序数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上的价值更多在于代码结构干净、参数少，适合作为元启发式算法入门以及后续改进的基线，而不是指望它碾压一切，真正决定调度质量上限的，永远是解码器和邻域搜索设计得够不够讲究。