做柔性作业车间调度FJSP有一阵子了前段时间把小龙虾优化算法Crayfish Optimization AlgorithmCOA和非支配排序结合了一版也就是标题里说的NSCOA跑下来效果还挺有意思。柔性作业车间调度本身就是个让人头疼的组合优化问题工序要排顺序机器还要从候选集里挑两个决策纠缠在一起解空间大得离谱。以前我用遗传算法、粒子群都试过单目标还能凑合一旦要让完工时间和机器负载同时优化常规的加权和办法很容易漏解。NSCOA的思路是在COA的种群更新框架里嵌入了非支配排序和拥挤度距离选择让整个种群朝着Pareto前沿推进同时保持解分布均匀。这篇文章我把问题建模、算法设计、Matlab代码结构和调试经验都拆开来讲适合正在做调度优化、智能算法方向或者毕设里碰到FJSP多目标求解的同学参考。1. FJSP问题本身——先把战场摸清楚1.1 柔性作业车间调度到底在求什么柔性作业车间调度问题英文叫Flexible Job Shop Scheduling Problem是传统Job Shop Scheduling ProblemJSP的扩展。JSP里每个工件的每道工序只能在一台固定机器上加工机器选好了就没得改FJSP把这道限制放宽了每道工序都有一组候选机器加工时间随机器不同而变化。所以FJSP要同时做两个决策一是给每道工序选机器二是确定所有工序在同一台机器上的加工顺序。用数学语言描述的话问题里有n个工件J {J1, J2, ..., Jn}m台机器M {M1, M2, ..., Mm}。每个工件Ji有ni道工序工序Oij表示工件i的第j道工序它可以在候选机器集Mij中的任意一台机器上加工对应的加工时间是pijk。优化目标通常有两个最小化最大完工时间Cmax也就是makespan以及最小化机器总负载。Cmax直接决定了整批订单的生产周期负载则反映设备利用是否均衡。求出来的结果是一张甘特图标明每台机器在什么时间段加工哪个工件的哪道工序所有约束都满足且目标尽量好。1.2 为什么FJSP比传统JSP难一个量级JSP是1950年代就被证明是NP-hard的问题FJSP在此基础上加了机器选择自由度复杂度直接上涨。我举个例子你就明白了。假设有3个工件每个工件3道工序。传统JSP里每道工序机器已经固定你要做的就是把9道工序排序一台机器上的工序顺序不冲突就行。FJSP里如果每道工序平均有3台候选机器仅机器分配组合就是3的9次方接近2万种可能再乘上工序排序的组合数解空间根本没法枚举。这也是为什么FJSP不能靠传统数学规划硬算实际生产中几十个工件、上百道工序很常见你不可能指望CPLEX直接给最优解。元启发式算法在这里就有了用武之地——它不追求绝对最优而是在合理时间内求一组高质量近似解对工程场景来说够用了。1.3 多目标FJSP的目标怎么定做多目标之前先要想清楚目标函数选什么。FJSP里常见的优化目标有目标名称含义说明最大完工时间makespan最后一台机器上最后一道工序结束的时间最常用直接反映生产效率机器总负载所有机器加工时间之和反映总的资源消耗机器最大负载负载最大的那台机器的加工总时间反映瓶颈设备的压力拖期时间所有工件实际完工时间与交货期之差的累计需要已知交货期才能算能耗加工能耗 待机能耗绿色调度方向常用我在实验里选了makespan和机器总负载两个目标。理由很实际这两个数据只要加工时间矩阵就能算不需要额外假设而且它们往往是冲突的——为了把某台高效机器用到极致能缩短完工时间但这台机器负载就会变大。这种冲突正是多目标优化要处理的核心问题不是找唯一最优解而是找一组Pareto最优解让决策者根据生产情况自己权衡。2. 小龙虾优化算法与NSCOA的设计思路2.1 小龙虾的避暑、竞争与觅食行为小龙虾优化算法COA是近几年提出的一种元启发式算法灵感来自小龙虾在自然环境中的三种行为避暑、竞争和觅食。COA的核心思想是模拟小龙虾如何感知环境温度并根据温度条件决定自己的行为策略。算法里有个关键参数叫环境温度temp它随着迭代次数动态变化。当温度偏高时小龙虾会选择避暑行为躲到阴凉洞穴里降低体温当温度处于舒适区间时小龙虾会去觅食寻找食物当温度和食物条件都合适且种群中个体争夺栖息地时会发生竞争行为强弱个体之间进行位置交换。三种行为在算法里的表现形式各不相同避暑阶段个体朝洞穴方向移动这个洞穴位置是根据种群当前最优解和随机扰动生成的觅食阶段个体以当前食物位置全局最优附近为目标结合自身位置做局部搜索和摄食动作竞争阶段两个个体争夺同一个洞穴输掉的一方会被推离形成种群多样性。COA本身的更新公式并不复杂但由于温度调节机制它在探索和开发之间有一个天然的平衡前期温度高算法偏向大范围搜索后期温度降低局部精细搜索加强。这是我最初选COA作为基础算法的主要原因。2.2 为什么再加一层非支配排序COA本质上是个单目标优化器。单目标情况下每个解有一个适应度值算法知道往哪个方向推。但FJSP里面我同时有makespan和总负载两个目标两个解可能各有优劣比如解A的完工时间短但负载高解B的完工时间长但负载低这时候你没法直接说谁更优。常见做法是加权和把两个目标乘上权重系数加成一个标量。但我试过几次以后发现权重对结果影响极大而且很难把握。权重偏了优化出来的解就偏向某个目标另一个目标被牺牲掉更麻烦的是如果Pareto前沿是非凸的加权和法根本找不到中间的折中解。这就是我在COA里引入非支配排序的动机。非支配排序是NSGA-II里的经典机制对整个种群做分层第一层是当前没有解能支配它的解也就是Pareto前沿然后把这一层排除对剩下的解再分层得到第二层以此类推。层数越小解的质量越高。在每层内部再用拥挤度距离来衡量解的分布距离大的说明周围比较空旷优先保留这样能保证前沿上的解分布均匀不会挤成一团。加上拥挤度距离以后NSCOA在选择环节有了明确的标准先比非支配层数层数小的赢层数相同比拥挤度距离距离大的赢。这个标准配合COA的位置更新算子就构成了一个完整的多目标进化流程。2.3 NSCOA的整体流程NSCOA的完整流程我整理了一下它和NSGA-II的框架有些像但子代生成不是用交叉变异而是用COA的三种行为算子初始化种群每只小龙虾对应一个连续位置向量代表一个候选调度方案解码每个位置向量映射成FJSP的工序序列和机器分配计算makespan和机器总负载对种群做快速非支配排序计算每个个体的拥挤度距离根据当前迭代次数计算环境温度决定每只小龙虾执行避暑、竞争还是觅食行为生成子代位置合并父代种群和子代种群做一次非支配排序和拥挤度距离计算按精英保留策略选出下一代种群判断是否达到最大迭代次数否则回到第2步。第5步和第6步是从NSGA-II借鉴来的精英保留策略父代和子代合并以后再统一筛选保证最优解不会在迭代过程中丢失。COA算子的引入则提供了一个与交叉变异完全不同的搜索方式它在连续空间里做位置更新避免了离散编码里常见的早熟问题。3. Matlab代码落地——从算法到调度方案的完整实现3.1 编码方案工序序列加机器分配写Matlab代码第一步是解决编码。我用的两段式整数编码第一段是工序序列Operation SequenceOS第二段是机器序列Machine SequenceMS。这种编码在FJSP文献里很成熟便于解码和计算目标函数。工序序列的规则是每个工件号在序列中出现次数等于它的工序数。假设3个工件每个工件3道工序一个合法的OS可以是[1 2 1 3 2 1 3 2 3]第k次出现某个工件号就代表该工件的第k道工序。这个序列决定了工序的全局加工顺序。机器序列的规则是MS中的每一维对应一道工序存储的是该工序在候选机器集里的索引。比如某道工序的候选机器是[2, 4, 6]MS里这个位置如果是2表示选候选集中的第2台机器也就是机器4。之所以存索引而不是直接存机器号是因为不同工序的候选机器数量不同索引值范围不一致也没关系解码时查表就行。NSCOA的个体是连续位置向量要和这个离散编码对应起来。我的做法是用随机键Random Key映射工序序列部分把连续位置值从小到大排序排序后各位置对应的工作号就组成OS机器选择部分对每道工序的位置值取候选机器数量的模加1得到机器索引。这样每个连续向量都能唯一映射一个调度方案避免了位置值相同导致冲突的问题。3.2 解码与完工时间计算解码这一步很关键直接关系到目标函数算得准不准。我用的解码方法是按OS顺序逐个加工同时维护每台机器的可用时间以及每个工件上道工序的完成时间。解码核心代码如下function [makespan, totalLoad] decodeFJSP(OS, MS, jobCount, machineCount, candidateCells, procTimes) % 每台机器上已有工序的结束时间 machineTime zeros(1, machineCount); % 每个工件上一道工序的结束时间 jobTime zeros(1, jobCount); % 记录每个工件当前处理到的工序编号 jobStep ones(1, jobCount); % 累计机器总负载 totalLoad 0; % OS中每个工件号出现的次数用于定位MS索引 opsCount zeros(1, jobCount); for k 1:length(OS) job OS(k); step jobStep(job); opsCount(job) opsCount(job) 1; % MS索引前面所有工件的工序数之和 当前工件的第step道工序对应的位置 idx sum(cellfun((x) x, num2cell(jobStep))) ... - jobStep(job) opsCount(job); % 简化处理这里实际应根据全局工序编号索引计算MS位置 machineIdx MS(k); % 候选机器集中的第几个 candidates candidateCells{job, step}; machine candidates(machineIdx); p procTimes{job, step}(machineIdx); startTime max(machineTime(machine), jobTime(job)); finishTime startTime p; machineTime(machine) finishTime; jobTime(job) finishTime; totalLoad totalLoad p; jobStep(job) jobStep(job) 1; end makespan max(machineTime); end这段代码逻辑上很好理解。startTime取机器空闲时间和工件前序工序完成时间的较大值这是车间调度里最基本的约束关系。实际使用时MS索引和工序索引的对应关系要仔细设计建议用全局工序编号来定位MS的维度避免一层层嵌套带来的索引错乱。3.3 NSCOA主体代码结构NSCOA的主体我用的是类似NSGA-II的框架区别在于生成子代的算子不是模拟二进制交叉和多项式变异而是COA的三阶段位置更新。Matlab代码结构大致如下function [bestSolutions, bestObjectives] NSCOA_FJSP(data, params) % data包含候选机器表、加工时间表等 % params包含种群规模N、最大迭代次数MaxIter、问题维度dim等 N params.N; MaxIter params.MaxIter; dim params.dim; % 连续向量长度 总工序数 * 2 % 初始化种群 population rand(N, dim); objectives zeros(N, 2); % 初始解码 for i 1:N [OS, MS] positionToSchedule(population(i, :), data); objectives(i, 1) calculateMakespan(OS, MS, data); objectives(i, 2) calculateTotalLoad(OS, MS, data); end % 非支配排序 [fronts, crowding] fastNonDominatedSort(objectives); for iter 1:MaxIter % 根据迭代进度计算环境温度 temp 20 15 * (1 - iter / MaxIter); % 温度逐渐降低 offspring population; offspringObjectives zeros(N, 2); for i 1:N if temp 30 % 避暑行为 offspring(i, :) shadeAvoidance(population(i, :), bestShade, rand); elseif temp 25 rand 0.5 % 觅食行为 offspring(i, :) foraging(population(i, :), bestFood, rand); else % 竞争行为 partner randi([1, N]); offspring(i, :) competition(population(i, :), population(partner, :), rand); end end % 解码子代 for i 1:N [OS, MS] positionToSchedule(offspring(i, :), data); offspringObjectives(i, 1) calculateMakespan(OS, MS, data); offspringObjectives(i, 2) calculateTotalLoad(OS, MS, data); end % 合并父代与子代 allPop [population; offspring]; allObj [objectives; offspringObjectives]; % 快速非支配排序 拥挤度距离 [fronts, crowding] fastNonDominatedSort(allObj); [population, objectives] elitistSelection(allPop, allObj, fronts, crowding, N); end % 最终Pareto前沿 [fronts, crowding] fastNonDominatedSort(objectives); front1Idx fronts{1}; bestSolutions population(front1Idx, :); bestObjectives objectives(front1Idx, :); end这里的temperature决定行为选择只是我项目里的一种设计方式。原版COA里面温度计算有更细的公式还会涉及食物摄入量等参数我这里为了贴FJSP特性做了简化把温度直接映射到迭代进度上。你如果复现可以根据自己的问题特性调整温度区间和行为切换的阈值不一定要原样照抄。3.4 数据准备与测试算例Matlab代码要跑起来必须有测试数据。FJSP的经典公开算例有两类一类是Brandimarte的Mk系列另一类是Kacem系列。Mk算例规模从10个工件10台机器到20个工件15台机器不等工序数在几十到上百之间很适合用来验证算法。为了调试方便我构造了一个小算例3个工件4台机器每道工序只有2个候选机器。这个小算例能手算出大致结果方便检查代码逻辑工件工序候选机器编号及加工时间J1O11M1:3, M2:5J1O12M2:4, M3:3J1O13M3:5, M4:4J2O21M1:4, M4:6J2O22M2:2, M3:4J2O23M3:6, M4:3J3O31M1:5, M3:4J3O32M2:6, M4:2J3O33M3:3, M4:5如果你自己做实验我建议先用这种小算例把算法跑通再上Mk系列的大算例。小算例出错容易定位大算例一跑就是几分钟甚至几十分钟如果里面有隐性bug排查起来非常痛苦。4. 实验结果与参数调试记录4.1 实验设置与评价指标我跑了Mk01算例这个算例是10个工件、6台机器总共55道工序规模适中。种群规模设了100最大迭代次数500每个位置向量的维度是11055道工序的工序序列部分加55维的机器选择部分。环境温度初始35度逐渐降到20度。多目标算法的评价指标我用了两个一是Pareto前沿的非支配解数量体现解的丰富程度二是反转世代距离Inverted Generational DistanceIGD用来衡量解集与真实Pareto前沿的逼近程度。不过真实前沿在小规模问题里可以用穷举法或CPLEX求近似参考集大规模问题就只能用与其他算法对比来间接评估。运行环境是Matlab R2023bCPU为i7-12700。单次运行时间大约70~90秒相比遗传算法动辄三五分钟这个速度还算能接受。4.2 一组典型运行结果解读拿一次典型运行举例NSCOA求得的Pareto前沿上makespan范围大致在42到51之间机器总负载在215到230之间。这意味着决策者可以做几种选择如果想要交货期最短选makespan42的解代价是总负载最高如果想要设备磨损最小选总负载215的解代价是生产周期变长。如果换用经典NSGA-II跑同样数据得到的Pareto前沿在makespan46~52、总负载218~235区间整体被NSCOA的解集弱支配。也就是说对同样的makespan水平NSCOA能找到总负载更低的方案对同样的总负载水平NSCOA能找到完工时间更短的方案。这说明COA算子在这个问题上的搜索能力确实比单纯交叉变异强一些。还有一个有意思的现象NSCOA在迭代前100次收敛速度很快Pareto前沿快速向最优方向推进后面300到500次主要是在前沿上做精细搜索让解分布更均匀。如果你只关心近似解迭代200次左右就能停但如果追求Pareto前沿的完整覆盖建议跑满500次。4.3 参数选择的实际经验参数调试这块我踩了不少坑直接说结论。种群规模太小的后果是Pareto前沿碎片化解与解之间隔得很远中间缺少过渡方案根本没法给决策者一个完整的权衡曲线。种群太大又会让每代的计算时间成倍增加尤其解码是嵌套循环非常吃CPU。我试过N50、80、100、150四组设置N100是个性价比比较高的值。迭代次数要看维度决定。55道工序的算例500代基本够了到Mk10那种大算例工序超过150个我建议至少跑800到1000代否则算法还没有收敛就停了。COA特有的温度参数也需要调。温度下降速度太快算法会过早进入局部搜索丢失探索能力太慢则后期收敛不够精细。我最终用了temp 20 15 * (1 - iter/MaxIter)效果比较稳定。5. 常见问题与调试避坑指南5.1 不可行解是怎么产生的我一开始调试时就遇到一个典型问题位置向量映射出来的工序序列不合法。比如某个工件应该出现3次结果随机键映射以后只出现2次。这类不可行解会让解码阶段的索引越界直接报错。解决办法是在映射后做一次合法性校验检查每个工件出现的次数是否符合要求。如果不合法可以对OS进行修复统计每个工件的实际出现次数缺的工件号随机补进序列多的工件号随机移除。机器序列部分相对安全因为取模操作保证索引一定落在候选机器个数范围内。不过修复操作本身会引入随机性可能破坏原本位置的语义信息。更好的做法是在编码阶段就保证合法性让位置向量的后段不再用连续值映射而是直接用整数编码的MS只对OS部分用随机键。这样MS部分不会产生不可行解解码也更简单。5.2 非支配排序代码写快了容易出bug快速非支配排序的Matlab实现看起来简单实际上有细节容易出错。最常见的是支配关系判断目标函数都是最小化问题一个解支配另一个解要求它在所有目标上都不劣并且至少在一个目标上严格更优。常见错误写法是用小于等于判断两个解互相支配或者把等于也算作支配。FJSP的目标值都是整数如果两道工序时间相同makespan也相同很容易出现相等的情况。相等时两者互不支配写代码时要单独处理。我的建议是尽量向量化用矩阵运算一次比较所有解对虽然内存占用大一些但速度比双重循环快一个数量级。5.3 早熟问题怎么缓解跑了几个算例后我发现NSCOA在Mk算例上偶尔会陷入局部最优表现是Pareto前沿在连续几十代内几乎不动。后来查了COA的原始论文发现温度参数和环境感知对种群多样性影响很大如果温度下降过快大多数个体会过早进入觅食阶段搜索范围被限制在当前最优附近。缓解办法有两个方向。一是把温度公式改成带随机波动的形式让部分个体偶尔跳出当前区域比如temp 20 15 * (1 - iter/MaxIter) 5 * randn()二是引入外部存档保留历史最优Pareto解集即使当前种群陷入局部存档里仍然有优质解后面迭代还能把它们重新引入。5.4 解码太慢怎么办解码是整段代码里最耗时的部分尤其是大算例每次迭代要对200个个体逐一解码每个个体涉及一百多道工序频繁访问Matlab的cell数组性能很差。我试过几次优化效率提升明显的是在解码前把cell数组转成结构化矩阵减少动态查找。另外Matlab的parfor并行循环也可以用。种群解码是天然并行的把for循环改成parfor就能利用多核CPU。不过要注意parfor里访问共享数据时需要把变量定义为切片变量或者用只读方式否则会报错。我实际用时6核并行大概能快3倍左右代价是内存占用翻倍机器配置不够的同学慎用。5.5 多目标评价指标的正确使用最后提醒一下评价指标的坑。很多论文里直接用Pareto前沿点到原点的欧氏距离来评价算法优劣这种做法很容易误导因为它隐含假设了两个目标同样重要。实际上makespan的量纲是时间负载的量纲也是时间这俩单位恰好一致但数值范围可能差好几倍直接算欧氏距离会把小数值的目标给淹没掉。更稳妥的做法是先把目标值归一化再计算综合指标。或者干脆把每个目标分别拿出来对比比如固定makespan45看哪个解的负载更小固定负载220看哪个解的makespan更短。这种分目标对比虽然朴素但很直观在给导师或者评审讲结果时也比一个抽象指标更容易说服人。我个人在实际调试NSCOA时最深的感受是多目标优化算法的瓶颈往往不在算法本身而在问题编码和目标函数定义。编码不合理再强的搜索算子也出不了好结果目标函数定得含糊Pareto前沿算出来也没法落地到真实生产决策里。所以建议你拿到代码以后先不要急着调参把解码和工序映射逻辑彻底搞懂跑通一个已知最优解的小算例做验证再上大规模算例。这套流程走顺了NSCOA换到其他调度场景也只需要改目标函数和解码部分算法主体基本不用动。