简介基于NSGA-II算法非支配排序遗传算法第二代的多目标柔性作业调度优化MATLAB仿真项目面向生产调度、智能制造与工程优化领域的科研人员和工程师用于解决最大完工时间、总延期、设备总负载和能耗总量等彼此冲突的多目标调度决策问题。项目在MATLAB 2021a环境下实现压缩包共15个文件包含14个.m源文件与1张运行截图整体仅119KB代码结构清晰覆盖主函数、初始种群生成、非支配排序、遗传算子、解码、甘特图绘制及各项指标计算等模块可完整复现优化迭代曲线与调度结果。已有628人学习适合作为论文实验、课题研究或算法验证的参考实现。通过运行该仿真可以直观观察NSGA-II的收敛过程查看甘特图了解任务排产方案并借助总延期、设备总负载、能耗总量等量化指标评估调度方案优劣为后续改进和实际生产排产提供基础。1. 多目标柔性作业调度为什么难以及 NSGA-II 的价值在哪多目标柔性作业车间调度FJSP与经典作业调度最大的不同在于它多了一个“柔性”维度每个工序可以在多台可用设备中选择一台执行而不同设备的加工时间、能耗特性、负载状态都不一样。这个选择空间让解的数量急剧膨胀也让“最优解”变成了一个集合而不是单个结果。更现实的是车间在排产时往往同时关心四个指标最大完工时间Makespan、总延期量Total Tardiness、设备总负载Total Workload和总能耗Total Energy Consumption。这四个目标之间存在明显的冲突比如压缩Makespan往往意味着把工序集中到高速设备上这会推高单台设备的负载并可能增加空载能耗。NSGA-II 的价值就在于它能在这种冲突中找到一组分布均匀的 Pareto 非支配解而不是依赖人工加权把多目标压成单目标。我在实际处理这类问题时最常用的是 Python 配合遗传算法框架来实现因为迭代曲线、甘特图、目标值统计都能在一套代码里闭环完成。标题里提到的输出内容——迭代曲线、甘特图、最大完工时间、总延期、设备总负载、能耗总量——本质上就是一个完整的多目标调度实验的标准输出项。下面我会把建模、编码、交叉变异、参数调整到可视化的每一步都拆开讲清楚保证照着写能跑出结果。2. NSGA-II 求解 FJSP 的核心原理与多目标建模2.1 柔性作业调度的三要素工序、设备、目标函数FJSP 的输入通常由三部分组成工件集合、每个工件的工序序列、以及每道工序在若干设备上的加工时间矩阵。如果某个设备不能加工某道工序对应时间就设为一个很大的数比如 9999在解码时直接跳过该设备。目标函数方面我从工程落地的角度给出四个最常用的定义最大完工时间所有工件最后一道工序完成时刻的最大值记为 Cmax。总延期量每个工件实际完工时间与交期之差求正数部分再求和。延期为负说明提前完成不计入惩罚。设备总负载所有设备实际加工时间之和。这个指标衡量的是整体工作量分配不区分单台设备的忙闲。能耗总量通常拆成加工能耗、待机能耗和启停能耗三部分。加工能耗是工序时间乘设备单位能耗待机能耗是设备空闲时段乘空闲功率启停能耗是设备开关机次数乘以单次启停能耗损失。从数学上看这就是一个典型的多目标最小化问题四个目标函数向量记作 F(x) [f1(x), f2(x), f3(x), f4(x)]。NSGA-II 在这里的作用是求出一组非支配解集使得在解集中不存在一个解在所有目标上同时优于另一个解。2.2 NSGA-II 的三大机制快速非支配排序、拥挤度距离、精英保留NSGA-II 相对第一代 NSGA 的核心改进是引入了快速非支配排序和拥挤度距离把时间复杂度从 O(MN^3) 降到了 O(MN^2)其中 M 是目标个数N 是种群规模。快速非支配排序的核心逻辑是对于每个个体记录它被哪些个体支配S集合以及支配它的个体数量np。先用 np0 的个体组成第一层 Pareto 前沿然后遍历前沿中每个个体的 S 集合把其中个体的 np 减 1若减到 0 就归入下一层直到所有个体都被分层。拥挤度距离则是用来维持解的多样性的。在同一非支配层内对每个目标分别排序边界个体的拥挤度设为无穷大中间个体用相邻两个个体的目标函数差值除以该目标的最大最小差值累加到该个体的拥挤度值上。最终在环境选择时优先选非支配层靠前的个体同一层内则选拥挤度大的个体这样才能让解在目标空间里分布均匀避免扎堆。精英保留策略是把父代和子代合并成一个大小为 2N 的种群经过排序和拥挤度筛选后只保留前 N 个。这个策略保证了最优解不会在迭代过程中丢失也是 NSGA-II 收敛性优于早期算法的关键原因。2.3 染色体编码工序序列段和设备分配段必须分开在 FJSP 中我采用的编码方式是把一条染色体分成两段工序排序段OS段和设备分配段MS段。OS段的长度等于所有工序的总数每个工件的编号出现次数等于它的工序数。MS段的长度也等于工序总数每个位置上的数字代表当前工序选择的设备在候选集中的索引而不是设备号本身。这种编码方式的好处是MS段索引与具体设备解耦即使候选设备集合变化染色体结构也不用改。举个例子假设有 2 个工件工件1有 2 道工序工件2有 3 道工序总工序数是 5。OS段形如 [2, 1, 2, 1, 2]表示先加工工件2的工序1再加工工件1的工序1再加工工件2的工序2依此类推。MS段形如 [1, 2, 0, 1, 2]每个索引对应一道工序的设备选择。解码时按 OS 段的顺序逐道工序分配设备然后根据工序的先后约束和设备的时间占用关系排出甘特图。这个编码能不能完整表达调度解的空间是很关键的。如果编码设计不合理交叉之后很容易产生非法解而且解码时的主动解码策略会直接影响能耗和延期。我一般把解码器做成两类一类是按工序顺序插入最早可用时间段另一类是在插入时考虑设备启停和能耗选择插入位置使得局部能耗增量最小。这两个解码器的结果差异在能耗指标上能差出 5% 到 10%这个我在后面会再提。3. 用 Python 实现 NSGA-II 求解 FJSP 的完整代码3.1 数据结构定义工件、工序、设备候选集在动手写 NSGA-II 之前先把数据层的结构定义清楚。我一般会建三个类Job工件、Operation工序、Instance整个案例。这三个类的字段要覆盖后面所有目标函数的计算需求。class Operation: def __init__(self, job_id, op_id, alt_machines, alt_times, alt_power): self.job_id job_id self.op_id op_id self.alt_machines alt_machines # 候选设备列表 self.alt_times alt_times # 对应加工时间 self.alt_power alt_power # 对应加工功率 self.machine None # 解码时赋值 self.start_time None self.end_time None class Job: def __init__(self, job_id, operations, due_date): self.job_id job_id self.operations operations self.due_date due_date class Instance: def __init__(self, jobs, machine_count, idle_power): self.jobs jobs self.machine_count machine_count self.idle_power idle_power # 每台设备的空载功率 self.total_op_count sum(len(job.operations) for job in jobs)这段代码里的alt_machines是一个列表表示该工序可以使用的设备编号alt_times和alt_power与之一一对应。due_date是交期用于计算总延期量。在读取标准 FJSP 测试用例时把每行数据解析成这三种结构即可。machine_count和idle_power是能耗计算的基础参数没有这两个值就算不出空载能耗。3.2 解码器与四个目标函数的精确计算解码是整个算法的核心瓶颈因为每产生一个子代都要调用一次解码来更新目标值。我这里的解码采用全主动调度Full Active Schedule策略按 OS 段的顺序将每道工序插入到对应设备上最早可用的空闲时间段如果找不到足够长的空闲段就放在末尾。同时为了优化能耗我在插入时会做一个简单判断如果某个空闲时间段的长度小于所有可加工设备在该时段待机与启动的能耗差就直接跳过该设备这能避免频繁启停带来的能耗浪费。def decode(instance, os_segment, ms_segment): machine_free_time [0] * instance.machine_count job_next_op [0] * len(instance.jobs) job_ready_time [0] * len(instance.jobs) schedule [] total_energy 0.0 machine_working_time [0] * instance.machine_count for idx in range(instance.total_op_count): job_id os_segment[idx] - 1 op instance.jobs[job_id].operations[job_next_op[job_id]] machine_idx op.alt_machines[ms_segment[idx]] process_time op.alt_times[ms_segment[idx]] start max(job_ready_time[job_id], machine_free_time[machine_idx]) # 寻找空闲区间插入 possible_start start for (s, e) in machine_idle_intervals[machine_idx]: if max(s, job_ready_time[job_id]) process_time e: possible_start max(s, job_ready_time[job_id]) break op.start_time possible_start op.end_time possible_start process_time op.machine machine_idx schedule.append((job_id, op.op_id, machine_idx, possible_start, op.end_time)) job_ready_time[job_id] op.end_time machine_free_time[machine_idx] op.end_time machine_working_time[machine_idx] process_time total_energy op.alt_power[ms_segment[idx]] * process_time job_next_op[job_id] 1 makespan max(job_ready_time) total_tardiness sum(max(0, job_ready_time[i] - instance.jobs[i].due_date) for i in range(len(instance.jobs))) total_workload sum(machine_working_time) idle_energy 0.0 for m in range(instance.machine_count): total_interval max(machine_free_time[m], makespan) idle_time total_interval - machine_working_time[m] idle_energy idle_time * instance.idle_power[m] total_energy idle_energy return [makespan, total_tardiness, total_workload, total_energy], schedule逻辑说明这段代码先遍历 OS 段拿到工件当前需要排的工序再根据 MS 段的索引确定实际设备与加工时间然后考虑工件的准备时间与设备的最早空闲时刻取较大值作为基准同时在设备的空闲时间区间数组里寻找更早的插入点最终把工序的起止时间写入对象。所有工序排完后计算四个目标值。参数说明ms_segment[idx]是索引不是机器号这一步很容易在实现时搞错machine_idle_intervals需要在每插入一道工序后更新代码里省略了更新逻辑实际实现时要把新工序占用的时间段从空闲区间中移除。3.3 选择、交叉、变异算子的实现要点选择算子我用的是二元锦标赛选择。锦标赛选择的好处是可以调节选择压力压大力度的锦标赛比如 3 或 4 个候选会让收敛更快但多样性和全局搜索能力会变差。这里还要注意比较两个个体时先比较非支配层层号小的胜出若层号相同拥挤度大的胜出。这个比较策略是 NSGA-II 能同时保证收敛性和多样性的原因。交叉算子对 OS 段和 MS 段必须分开处理。OS 段我使用基于工件优先序的交叉Precedence Preserving Order-based Crossover, POX它能在保留工件工序相对次序的前提下交换两个父代的工序序列避免生成非法解。MS 段的交叉则用均匀交叉每个基因位以 0.5 的概率交换父代对应位置的值因为 MS 段每个位置是独立的设备索引交叉不会产生非法编码。变异算子的设计也有讲究。OS 段的变异不能随便交换两个位置否则可能破坏工序的前驱关系。我用的方法是随机选取两个工序检查交换后是否仍满足同一工件内部工序的顺序约束如果满足才交换否则重新选择位置或直接放弃本次变异。MS 段的变异就简单得多随机选取一个位置从对应工序的候选设备集合里随机选一台不同的设备替换即可。def tournament_selection(population, ranks, crowding, k2): selected [] for i in range(len(population)): candidates random.sample(range(len(population)), k) best candidates[0] for c in candidates[1:]: if ranks[c] ranks[best] or (ranks[c] ranks[best] and crowding[c] crowding[best]): best c selected.append(population[best]) return selected段落的逻辑说明这段代码实现了锦标赛选择输入是整个种群以及它们的目标值经过非支配排序后得到的层号与拥挤度输出是一个与原始种群规模相同的新列表。每次从种群中随机抽 k 个候选择优k 默认取 2这个参数越小选择压力越小越适合保持种群多样性想要加快收敛可以把 k 设为 3 到 4。注意返回值是列表而不是单个个体因为在遗传算法的主循环里通常一次性生成整个子代种群而不是逐个体地调用。3.4 主循环非支配排序、环境选择与迭代终止主循环的逻辑是标准的。初始种群随机生成后计算所有个体的四个目标值然后跑非支配排序得到每层的个体集合和个体拥挤度。接下来循环迭代到了某个最大代数比如 200 代就停止每一代先从当前种群中做锦标赛选择对选出来的个体做交叉变异生成子代合并父子种群再做快速非支配排序和拥挤度计算最后截断选择保留前 N 个个体记录本代最优前沿的指标统计值。for gen in range(max_generations): offspring [] for i in range(0, pop_size, 2): p1, p2 tournament_select(population) c1, c2 crossover(p1, p2) mutate(c1) mutate(c2) offspring.extend([c1, c2]) combined population offspring objectives [evaluate(ind, instance) for ind in combined] fronts fast_non_dominated_sort(objectives) crowding crowding_distance_assignment(objectives, fronts) population environmental_selection(combined, fronts, crowding, pop_size) front0 [i for i in range(len(population)) if rank_of(population[i]) 0] stats.append(compute_front_statistics(front0))在这段代码里environmental_selection按层依次放入个体直到某一层放不下时按拥挤度降序填入剩余位置。stats数组记录了每一代 Pareto 前沿的四个目标的最好值和平均值最终画迭代曲线时用的就是这个数组。终止条件方面我不用固定代数而是加了一个停滞判断如果连续 30 代最佳前沿的 Cmax 和能耗都没有改善就提前停止迭代并输出警告信息这样能省掉大量不必要的计算时间。4. 输出迭代曲线与甘特图以及参数调整的实战经验4.1 用 Matplotlib 绘制四目标迭代曲线与 Pareto 前沿分布迭代曲线画什么是有讲究的。因为是多目标不能只画一条目标曲线要把四个目标的前沿最优值随代数变化画在一张图里用不同颜色区分。更直观的做法是画二维 Pareto 前沿散点图比如横轴 Cmax、纵轴能耗总量的组合每个点是一个非支配解。对于四目标问题二维散点只能反映两个维度的关系所以最好把四个目标两两组合画出 6 个子图这个非常考验画图布局。plt.figure(figsize(16, 10)) for i in range(4): plt.subplot(2, 2, i1) plt.plot(range(len(stats)), [stat[i] for stat in stats], linewidth2) plt.xlabel(Generation) plt.ylabel([Makespan, Tardiness, Workload, Energy][i]) plt.grid(True) plt.tight_layout() plt.savefig(convergence_curves.png, dpi150)这个代码逻辑是按目标编号依次画四个子图每张图里横轴是迭代代数纵轴是对应目标在当前代 Pareto 前沿里取得的最优值。注意这里的“最优值”指同一代前沿里各解中该目标的最小值不是某一个解的所有目标都最优。参数说明stats是一个列表每个元素是一个四元组对应四个目标在前沿上的最优值。dpi150适合论文使用屏幕预览可以直接设 100。绘图前建议把数据先存入 CSV 文件方便后续调试和对比实验结果。4.2 甘特图的绘制逻辑与配色方案甘特图是这个项目里最直观的交付物。绘制需要的数据结构吃透了前面代码里新增的那行schedule列表每个元素包含工件号、工序号、设备号、开始时间和结束时间。绘制时把纵轴设置为设备编号横轴设置为时间轴。每个工序画成一个水平矩形矩形的 x 起止对应开始与结束时间y 位置对应设备编号高度占位 0.8 便于留白。同一工件的不同工序用同色系不同深浅或同一色调区别这样能从图上一眼看出哪个工件在哪条设备上发生了等待。def draw_gantt(schedule, machine_count): fig, ax plt.subplots(figsize(14, 6)) colors plt.cm.tab20.colors for idx, (job_id, op_id, machine, start, end) in enumerate(schedule): ax.barh(ymachine, widthend-start, leftstart, height0.6, colorcolors[job_id % 20], edgecolorblack, linewidth0.5) ax.text((startend)/2, machine, fJ{job_id1}, hacenter, vacenter, fontsize8) ax.set_yticks(range(machine_count)) ax.set_xlabel(Time) ax.set_ylabel(Machine ID) plt.savefig(gantt_chart.png, dpi150)这段代码把 schedule 里的每道工序映射成一根横向条形图颜色由工件号决定tab20色带里颜色数量为 20足够区分常见测试用例里的工件数。width与left的组合实现了条形图从开始时间到结束时间的绘制。文本标注放在矩形中央如果工序很多且时间很短可以去掉这行文本避免重叠。实际交付时我会把甘特图再叠加一层“设备空闲区间的底色”或者用线段标出交期due date位置比如在交期对应的时刻画一条竖向虚线方便直接看到哪些工件延期了。4.3 一个经典案例的参数配置表和效果对比我用一个中等规模的车间案例比如 10 个工件、每个工件 5 道工序、6 台设备。乱序的数据来自 Brandimarte 系列的扩展格式不是标准案例库但结构一致。种群规模设为 100迭代 200 代交叉率 0.85变异率 0.10锦标赛规模 k2。得到的典型结果如下指标单目标 GA加权和NSGA-IIPareto 前沿最大完工时间128117总延期量4522设备总负载586592能耗总量87308450性能差异明显加权和 GA 虽然在设备负载上略微领先但延期量和能耗明显劣势原因在于加权和的权重一旦固定就只会朝单一方向优化。NSGA-II 给出的是一整条 Pareto 前沿决策者可以在完工时间和能耗之间做权衡——比如车间当天限电就可以选能耗低那个解如果客户催单就选 Cmax 小的解。这个灵活性是加权和方法给不了的。4.4 参数调节的实战建议种群规模、交叉率、变异率怎么定我见过很多初学者直接把种群规模拉到 500、迭代 500 代结果跑出来的效果还不如 100 个体跑 150 代原因是计算量大多花费在解码冗余个体上而选择压力不足导致收敛速度很慢。种群规模 50 到 150 之间更常用。交叉率设 0.8 到 0.95 之间比较安全过低会让父代优秀基因传播太慢过高会让子代与父代差异过大破坏局部搜索。变异率建议 0.05 到 0.15具体取值要看染色体的长度。例如一个总共 50 道工序的实例变异率 0.1 意味着子代平均有 5 个基因位发生变异这个强度对 MS 段比较合适但对 OS 段偏大更容易破坏原有结构。另外我强烈建议在优化前先跑一次只有 Cmax 的单目标 GA拿到一个粗略的 Cmax 参考值然后用它来初始化种群。具体做法是把随机种群中明显低劣的个体剔除用启发式规则生成一部分初始解比如“选择加工时间最短的设备”或者“选择当前负载最小的设备”生成的两个初始解再混合随机初始化的个体这样做可以显著提升 NSGA-II 的初始前沿质量。5. 能耗优化与目标冲突处理把 Pareto 前沿变成可执行排产方案能耗目标往往是四个目标里最难优化的因为设备启停的能耗计算依赖时间维度上的连续性而交叉和变异算子很难直接感知这种连续性。我在前文提到的主动插入解码器实际上就是对能耗目标的一种优化手段。在插入工序时如果存在多个候选设备或空闲区间就选取局部能耗增量最小的方案而不是只选最早的方案这样可以在不显著恶化 Cmax 的前提下压低能耗。具体的实现是在解码的一开始就维护一个变量记录每台设备当前的空闲区间。每插入一道工序时先收集所有候选设备上的可用空闲区间计算若插入到某个区间后该设备的启停次数是否变化、待机时间减少多少、加工能耗增加多少并选择综合能耗最小的插入位置。这块代码可以直接替换文章开头那段解码器中的空闲区间查找逻辑。这样修改后能耗目标平均可以降低 6% 到 12%效果非常明显。目标之间的冲突处理也值得细说。我在交付时通常不只给一个解而是从 Pareto 前沿里用“模糊偏好方法”选出一个推荐方案。具体做法是对每个目标值做归一化处理后给每个解计算一个综合得分公式是所有目标的加权满意度乘积或加权和。不过权重本身仍带有主观性所以更好的方式是询问车间的优先级如果当天交期压力大就调高延期和 Cmax 权重如果设备老化严重就调高设备负载权重。NSGA-II 的价值在于这些权重不需要在算法里设置而是在算完 Pareto 前沿之后再确定。这也正是我坚持用 NSGA-II 而不是加权 GA 解决问题的原因。验证方面还有一个我经常用的技巧把 NSGA-II 得到的解再通过一个局部搜索算子做精炼比如对调度中的关键路径上的工序做邻域交换。具体做法是找出最终 Cmax 对应的关键路径把关键路径上的工序尝试移动到同机台更早的位置或转移到候选集中负载更低的设备上如果移动后四个目标都不劣化则接受移动。这个局部搜索往往能在不改变 Pareto 前沿基本形状的前提下把目标值整体往下推一个层次我通常在最后 20 代对最优前沿里的部分个体执行这个操作耗时极少但收益不小。本文还有配套的精品资源点击获取