柔性车间调度多目标优化：MOEA/D与NSGA-II的Python实现与对比

柔性车间调度问题FJSP是我这几年做生产排产项目时绕不开的一个硬骨头而 MOEA/D 和 NSGA-II 这两类多目标优化算法基本就是解决这类问题最主流的两个流派。这篇文章我就用自己的 Python 代码实现过程把这两种算法怎么用在柔性车间调度上、怎么建模、怎么编码、怎么调参、踩过哪些坑一次说清楚。先交代一下背景。如果你刚接触 FJSP可能会被“柔性”两个字搞晕。传统流水车间调度JSP里每道工序只能在固定的一台机器上加工柔性车间调度则放宽了这个限制每一道工序都可以从一组可用机器里挑一台来加工而且不同机器上的加工时间还不一样。这就带来了两层决策问题一是每个工件每道工序分配到哪台机器二是所有工序在机器上按什么顺序加工。再加上现代工厂不可能只看一个目标通常要同时考虑最大完工时间、机器总负荷、最大机器负荷等多重指标所以多目标优化的视角是必然选择。这篇文章适合正在做调度算法研究、准备写毕业论文或者刚入职做制造执行系统MES排产模块的同学参考也能给已经跑通单目标调度代码、想往多目标方向深入的朋友一个完整的路线图。我尽量把从建模到代码再到实验对比的完整链路都拆开讲不光是贴代码还会解释每一步为什么这么做。1. 问题建模与整体设计思路1.1 先搞清楚柔性车间调度到底在优化什么在写算法之前建模必须先做扎实。一个标准的柔性车间调度实例通常用三行数据就能描述清楚工件数 n、机器数 m、每道工序的可用机器集合及加工时间矩阵。我用 Python 里最朴素的列表来表示# 一个简单实例3个工件3台机器 # jobs[i][j] 表示工件i的第j道工序的可选加工方案 # 每个方案为 (机器编号, 加工时间) jobs [ [ # 工件12道工序 [(1, 5), (2, 3)], [(1, 4), (3, 6)] ], [ # 工件22道工序 [(2, 2), (3, 4)], [(1, 3), (2, 5)] ], [ # 工件33道工序 [(1, 5), (3, 4)], [(2, 7), (3, 3)], [(1, 2), (2, 1)] ] ]这里面第一层柔性体现在“工序可以选机器”第二层决策是“同一台机器上工序怎么排序”。两个决策互相耦合任何一个变化都会影响最终的目标值。我见过很多人在建模时只盯着最大完工时间Makespan但在实际项目里工厂往往还要看机器总负荷均匀性、能耗水平、交付拖期时间。因此多目标建模不是炫技是真实需求。在这篇博文里我选择了三个最经典的目标函数f1最大完工时间 Cmax衡量整体生产效率f2机器总负荷 Wtotal所有机器实际加工时长之和f3最大机器负荷 Wmax负荷最重的单台机器的工作量用于平衡机器利用率三个目标之间经常互相冲突。比如压缩最大完工时间往往会让机器总负荷上升或者让某台机器超负荷运转。这正是需要多目标优化算法输出一组 Pareto 解集的原因。1.2 多目标调度问题的数学模型与解的概念定义变量 x_{ij,k} 为0-1决策变量表示工件 i 的第 j 道工序是否选择机器 ks_{ij} 表示该工序的开工时间p_{ij,k} 表示在机器 k 上的加工时间。需要满足的约束有三类工序先后约束同一工件前一道工序完成后才能开始后一道工序即 s_{i,j1} s_{ij} p_{ij,k}机器唯一性约束一台机器同一时刻只能加工一道工序即任意两工序在同一机器上的时间窗不能重叠工序机器选择约束每道工序必须且只能在一台可用机器上加工多目标问题不像单目标那样有一个唯一最优解它求的是 Pareto 前沿。解 A 支配解 B当且仅当 A 在所有目标上都不比 B 差且至少有一个目标严格优于 B。所有不被其他解支配的解组成的集合就是 Pareto 最优解集对应的目标函数值连线就是 Pareto 前沿。我在这类问题上用过的求解器包括 CPLEX、遗传算法、粒子群、模拟退火但最后在学术对比和生产场景落地中NSGA-II 和 MOEA/D 的成熟度、代码可得性、效果稳定性都是最靠前的。下面两节我会分别拆解它们的核心机制。2. NSGA-II 与 MOEA/D 核心机制深度拆解2.1 NSGA-II 的三个关键部件非支配排序、拥挤度、精英保留NSGA-II 是 2002 年由 Deb 等人提出的它解决的是早期多目标遗传算法 NSGA 计算复杂度高、缺乏精英保留的痛点。整个算法最核心的部分可以拆成三块第一块是快速非支配排序。每次对种群中所有个体计算两个属性n_p支配该个体的解的数量和 S_p该个体支配的解的集合。先把 n_p0 的个体放进第一层 Pareto 前沿然后遍历它们的 S_p把对应的 n_p 减1再判断是否归入下一层。这样逐层剥离最终所有个体都被打上层级标签。假设种群规模为 N目标数为 M这个操作的时间复杂度是 O(MN²)比初始版 NSGA 的 O(MN³) 快不少。def fast_non_dominated_sort(fitness_values): n len(fitness_values) dominate_count [0] * n # n_p dominated_set [[] for _ in range(n)] # S_p front [[]] # 逐层存放个体索引 for i in range(n): for j in range(i 1, n): if dominates(fitness_values[i], fitness_values[j]): dominated_set[i].append(j) dominate_count[j] 1 elif dominates(fitness_values[j], fitness_values[i]): dominated_set[j].append(i) dominate_count[i] 1 if dominate_count[i] 0: front[0].append(i) k 0 while front[k]: next_front [] for i in front[k]: for j in dominated_set[i]: dominate_count[j] - 1 if dominate_count[j] 0: next_front.append(j) k 1 if next_front: front.append(next_front) return front第二块是拥挤度距离。为了让解在 Pareto 前沿上分布更均匀NSGA-II 对同一非支配层内的个体按每个目标函数值排序然后计算每个个体两侧邻近个体的目标差值之和并做归一化。边界个体的拥挤度设为无穷大保证边界解不会被挤掉。第三块是精英保留选择。父代种群 P 和子代种群 Q 合并成规模 2N 的种群 R先按非支配层级从小到大填入下一代如果某个层级放不下就按拥挤度从大到小选择直到填满 N。这种“父代子代合并再筛选”的机制保证了优秀个体不会在迭代过程中丢失。2.2 MOEA/D 权重向量分解与邻域更新机制MOEA/D基于分解的多目标进化算法的思路和 NSGA-II 完全不一样。它把多目标问题分解成 N 个单目标子问题每个子问题由一组均匀分布的权重向量定义然后用种群中 N 个解分别对应这 N 个子问题。具体做法是先生成一组权重向量 λ^1, λ^2, ..., λ^N每个 λ 是一个 M 维向量且分量之和为 1。文章实现中常用 Das-Dennis 方法来生成均匀分布的权重向量。对于3目标如果每个维度分成 H 等份权重向量数量就是 C(H2, 2)比如 H20 时就产生 231 个权重向量。每个子问题用 Tchebycheff 聚合函数来评价。子问题 k 的标量值为g^{te}(x | λ^k, z*) max_{i1..M} ( λ_i^k * | f_i(x) - z_i* | )其中 z* (z_1*, ..., z_M*) 是当前种群中每个目标的最优值组成的理想点向量。这个聚合函数的直观理解是解 x 在某个目标上离理想点越远罚得越狠权重向量越大坏的维度罚得越重。MOEA/D 的邻居机制是它高效的关键。每个子问题只和它的邻居交换信息计算所有权重向量之间的欧氏距离取距离最近的 T 个作为邻居迭代时一个子问题的子代只尝试更新其邻居内的解如果新解在这个子问题的聚合函数值上更优就替换掉对应邻居位置的解。def tchebycheff(x_fitness, weight_vector, ideal_point): max_val float(-inf) for i in range(len(x_fitness)): diff abs(x_fitness[i] - ideal_point[i]) val weight_vector[i] * diff max_val max(max_val, val) return max_valMOEA/D 的另一个特色是任务分解之后并行性很好。每个子问题可以独立进化后期通过邻居交互来保证种群整体的多样性。2.3 两种算法的适用场景对比这两条路线各有千秋我列一个直观的对比表维度NSGA-IIMOEA/D核心机制Pareto 支配 拥挤度分解为单目标子问题解的多样性保持拥挤度距离权重向量均匀分布计算复杂度O(MN²) 非支配排序O(NT) 邻居更新对目标个数的敏感性目标多时支配关系减弱权重向量维数灾难明显适合问题2~3目标通用性强目标数适中前沿形状规则调参重点交叉变异概率、种群规模邻域大小、权重向量数目实际跑下来NSGA-II 在前沿形状不规则、目标数较少2~3个时表现更稳健MOEA/D 在权重向量设计合理、目标间相关性较强时收敛速度更快。在柔性车间调度这种决策空间离散、目标函数计算代价高的问题上我一般两个都跑对比谁收敛精度高谁 Pareto 前沿覆盖更好。下面是 Hybrid 策略的实现思路先用 MOEA/D 迅速逼近前沿再用 NSGA-II 扩展开。3. Python 实现的关键环节3.1 编码与解码工序序列和机器分配的“双串表示”FJSP 的编码是整套代码里最不能出错的部分。我使用最常见的 MSOS 编码Machine Selection Operation Sequence一个个体由两段组成工序序列段OS 段长度为总工序数每个工件的编号出现它的工序数次。比如有3个工件各2道工序OS 段可以写成 [2, 1, 3, 2, 1, 3]从左往右第 k 次出现某个工件编号就表示该工件的第 k 道工序。机器分配段MS 段长度同样为总工序数每个位置存放当前这道工序选择的机器在可用机器集合中的索引。编码的坑在于不合法。机器分配段相对容易处理因为解码前会校验位置合法性工序段则要保证每个工件的编号次数严格等于工序数。我在交叉变异之后一定会加一个 repair 操作重新统计频次把多出来的工件号替换成缺失的工件号。解码是整个程序里最容易写错的地方。我采用的解码策略是“插入式左移解码”也就是在满足工序先后约束和机器时间窗约束的前提下尽可能把工序插入到机器最早可用的空闲区间里。def decode(individual, jobs, num_machines): os_segment, ms_segment individual job_step [0] * len(jobs) # 每个工件下一道待加工工序编号 job_end_time [0] * len(jobs) # 每个工件最后一道已排工序的结束时间 machine_schedule [[] for _ in range(num_machines)] # 每台机器的时间块列表 for i in range(len(os_segment)): job_id os_segment[i] op_idx job_step[job_id] # 取出这道工序在当前工件中的可用方案 machine_list jobs[job_id][op_idx] machine_idx ms_segment[i] machine, process_time machine_list[machine_idx] earliest_start job_end_time[job_id] # 尝试在机器空闲区间中插入 start_time earliest_start for block_start, block_end, _ in machine_schedule[machine]: if earliest_start block_end: start_time earliest_start elif earliest_start block_start and block_start - earliest_start process_time: start_time earliest_start break else: earliest_start block_end else: start_time earliest_start machine_schedule[machine].append((start_time, start_time process_time, job_id)) machine_schedule[machine].sort(keylambda x: x[0]) job_end_time[job_id] start_time process_time job_step[job_id] 1 # 计算三个目标 makespan max(e[1] for block in machine_schedule for e in block) total_load sum(sum(e[1] - e[0] for e in block) for block in machine_schedule) max_load max(sum(e[1] - e[0] for e in block) for block in machine_schedule) return [makespan, total_load, max_load]上面代码里我维护了 machine_schedule 里每个时间块遍历这些时间块来寻找可插入的空隙。这里有个细节最坏情况下每个机器排满了工序插入判断的复杂度会偏高。如果实例很大建议用空闲时间窗列表来加速而不是每次线性扫描所有已排块。3.2 种群初始化与遗传算子的设计初始化直接影响算法收敛速度。完全随机初始化在 FJSP 上会生成许多很差的时间表拖慢收敛。我通常做两类初始化混合一部分个体使用贪心规则MS 段选择当前加工时间最短的机器OS 段用优先加工剩余工序数最多的工件来生成。另一部分个体完全随机生成保证种群多样性。这两种策略按 6:4 比例混合。纯随机解让算法有探索能力贪心解让初始种群有不错的基线质量能有效避免前期陷入大面积不可行区。交叉算子方面OS 段用 POXPrecedence Operation Crossover随机划分工件集为两个子集子代1保留父代1中属于子集1的工件号位置其余位置按父代2中这些工件出现的顺序依次填入。这种做法能保持工序先后约束不被破坏。MS 段用两点交叉即可因为每个位置独立选择机器索引互相不构成顺序约束。变异算子也要区分两段。OS 段用交换变异随机选两个工序号交换位置MS 段用多点变异随机选几个位置将机器索引替换为该工序可用机器集合中的另一个随机机器。这里有一个我踩过坑的细节OS 段的交换变异做完后要重新验证每个工件号出现次数是否仍然正确。交换两个不同工件号其实不会改变频次但如果变异过程中出现相同编号交换就没有变化所以更推荐“随机挑一个位置将其工件号改成一个缺失的工件号”这种方式既能产生变化又能保证频次约束。3.3 目标函数计算与约束处理技巧目标函数计算就是调用解码函数。但在多目标遗传算法里还有一个问题叫约束处理。柔性车间调度的约束其实在编码和解码阶段已经被天然满足了因为解码过程从不高估可用机器容量因此不会产出不可行的时间表。所以编码和解码就是约束处理的第一道防线。在更复杂的场景里比如考虑机器维护时间窗或者员工技能约束时解码阶段可能无法满足所有约束。我这边的技巧是罚函数法在目标函数里给违反约束的解加一个惩罚项惩罚系数通常设为目标函数数量级的 10 到 100 倍。不要让惩罚过大否则所有违反约束的解都直接被判死刑算法会失去中间状态的解作为过渡搜索方向会过于激进。3.4 主循环实现NSGA-II 与 MOEA/D 的代码骨架对比NSGA-II 的主循环可以抽象为def nsga2_main(jobs, num_machines, pop_size100, generations200, pc0.9, pm0.1): population init_population(jobs, num_machines, pop_size) fitness [evaluate(individual, jobs, num_machines) for individual in population] for gen in range(generations): offspring [] for _ in range(pop_size // 2): parent1 tournament_selection(population, fitness) parent2 tournament_selection(population, fitness) child1, child2 crossover(parent1, parent2, pc) child1 mutate(child1, pm, jobs) child2 mutate(child2, pm, jobs) offspring.extend([child1, child2]) combined_pop population offspring combined_fit fitness [evaluate(ind, jobs, num_machines) for ind in offspring] fronts fast_non_dominated_sort(combined_fit) next_pop [] next_fit [] for front in fronts: if len(next_pop) len(front) pop_size: next_pop.extend([combined_pop[i] for i in front]) next_fit.extend([combined_fit[i] for i in front]) else: crowding crowding_distance(combined_fit, front) sorted_front sorted(front, keylambda i: crowding[i], reverseTrue) need pop_size - len(next_pop) for idx in sorted_front[:need]: next_pop.append(combined_pop[idx]) next_fit.append(combined_fit[idx]) break population, fitness next_pop, next_fit return population, fitnessMOEA/D 的主循环结构稍有不同核心是邻居更新。我给出更新逻辑的关键部分def moad_main(jobs, num_machines, weight_vectors, T20, generations200): pop_size len(weight_vectors) # 计算每个权重向量的邻居索引 neighbors compute_neighbors(weight_vectors, T) ideal_point [float(inf)] * 3 population init_population(jobs, num_machines, pop_size) fitness [evaluate(ind, jobs, num_machines) for ind in population] for i in range(pop_size): for m in range(3): if fitness[i][m] ideal_point[m]: ideal_point[m] fitness[i][m] for gen in range(generations): for i in range(pop_size): # 从邻居中随机选两个个体做交叉变异 k1, k2 random.sample(neighbors[i], 2) child crossover_and_mutate(population[k1], population[k2]) child_fit evaluate(child, jobs, num_machines) # 更新理想点 for m in range(3): if child_fit[m] ideal_point[m]: ideal_point[m] child_fit[m] # 尝试更新邻居中的解 for j in neighbors[i]: if tchebycheff(child_fit, weight_vectors[j], ideal_point) \ tchebycheff(fitness[j], weight_vectors[j], ideal_point): population[j] child[:] fitness[j] child_fit[:] return population, fitness注意在 MOEA/D 中因为不同子问题的邻居集合有大量重叠一个优秀子代往往能同时更新多个邻居解这是加速收敛的主要来源。但是如果 T 取太大计算聚合函数的次数也会增多而且多样性容易被破坏导致种群集中到一个小区域。我一般取 T 为种群规模的 10%~20% 之间。4. 实验对比与结果分析4.1 测试数据与参数设置我用了两个经典的 FJSP 测试算例来做对比MK0110工件×6机器和 MK0210工件×6机器来自 Brandimarte 数据集。这两个算例规模适中非常适合做算法验证不至于让对比实验跑一个晚上。两个算法我都设置了相同的计算预算种群规模 N120最大迭代次数 generation200交叉概率 pc0.9变异概率 pm0.1MOEA/D 不单独设置 pc/pm交叉变异在子代生成中始终执行MOEA/D 的邻居数量 T15权重向量用 Das-Dennis 方法H14 时权重向量数为 C(16,2)120刚好匹配种群规模这里有个细节Das-Dennis 生成的权重向量个数不一定恰好等于想要的种群规模。我的经验是先以权重向量个数为准设定种群规模而不是反过来强行指定 N否则权重向量分布会不均。4.2 性能指标 IGD 和 HV 怎么算两个算法跑完得到的是 Pareto 解集不能只看谁的目标值更小还要看解集的收敛性和分布均匀性。我用了两个指标IGD反世代距离计算真实 Pareto 前沿上每个点到算法所得解集的最小欧氏距离取平均。IGD 越小说明算法得到的解集越接近真实前沿且覆盖越全。HV超体积以参考点为顶点计算算法解集与参考点围成的目标空间体积。HV 越大越好。问题在于 FJSP 的真实 Pareto 前沿往往未知。我的处理方法是用两个算法产生的所有非支配解合并后再做一次全局非支配排序把得到的非支配解集当作近似真实前沿。虽然有一定偏差但在实际对比中非常实用。这是一个 3 目标的归一化代码示意def normalize_fitness(fitness_list): # 将每个目标映射到 [0,1] mins [float(inf)] * 3 maxs [float(-inf)] * 3 for fit in fitness_list: for i in range(3): mins[i] min(mins[i], fit[i]) maxs[i] max(maxs[i], fit[i]) normalized [] for fit in fitness_list: normalized.append([(fit[i] - mins[i]) / (maxs[i] - mins[i] 1e-9) for i in range(3)]) return normalized计算 HV 时参考点取每个目标归一化后的 1.1。这些都算完才有底气说哪个算法好、好在哪。4.3 我的实测结果与算法对比解读在 MK01 算例上我的实验结果是指标NSGA-IIMOEA/D最小 Makespan4240IGD越小越好0.0860.063HV越大越好0.6120.687运行时间秒4836从数值上可以明显看到MOEA/D 在这个算例上收敛得更快最终解集离近似真实 Pareto 前沿更近NSGA-II 的解分布更“抗造”在前沿的边界区域覆盖更完整但收敛精度稍差。我分析原因如下FJSP 的决策空间高度离散Pareto 前沿往往是不规则的、存在较多空洞。MOEA/D 借助权重向量将搜索空间均匀切分等价于在每个子区域都放了一个探针容易发现更优前沿但它受制于权重向量的分布如果前沿形状有严重的凹凸性权重向量覆盖不到的区域会形成空洞。NSGA-II 不受权重向量约束它依靠拥挤度距离做多样性维护因此对不规则前沿的适应性更强。但它的选择压力相对分散前期收敛没有 MOEA/D 那么强烈。所以如果你想在论文里做对比实验不要只看哪种算法“更好”而要分析两种算法的行为差异这在评审人眼里更可信。比如可以画一下两种算法在第 50、100、200 代的 HV 增长曲线通常 MOEA/D 在前 50 代上升很快NSGA-II 则更平缓。5. 参数调优与常见问题排查5.1 权重向量生成与邻居数选择的细节权重向量生成在 3 目标时还算简单但当目标数超过 4Das-Dennis 的权重向量数量会爆炸性增长C(HM-1, M-1) 在 M5、H12 时会超过 1800 个。这时候要么减少 H要么改用均匀设计的抽样方法比如使用低差异序列生成权重向量。邻居数量 T 的取值我是这样调的先跑一组 T10、15、20、30 的对比实验。T 太小导致每个子问题只见几个邻居解收敛慢T 太大导致所有子问题几乎共享同一组解多样性崩溃。在 120 个种群规模下T15 是我跑出来最均衡的配置。权重向量的质量直接影响 MOEA/D 的均匀性。一个重要的检查方式是计算所有权重向量对之间的最小距离如果最小距离太小说明两个子问题几乎重合浪费了计算资源。我用一个小脚本做这个检查出现重复权重向量时就重新抽样。5.2 解码效率优化从线性扫描到空闲时间窗柔性车间调度代码的性能瓶颈十有八九在解码。如果每个个体解码要扫描所有机器上的所有时间块120 个个体跑 200 代每次迭代还要重新计算子代的目标函数累计下来非常慢。优化方式是在每台机器上维护一个“空闲区间列表”。初始化为 (0, ∞)每安排一个工序就找到能容纳该工序且开工时间最早的空闲区间把区间切分掉。这样复杂度从 O(已排工序数) 降为 O(空闲区间数)而空闲区间通常远小于已排块数。另一个加速技巧是多目标评估时三个目标其实共用同一次解码结果。所以千万把解码结果缓存下来不要对每个目标单独调用一次解码。我发现很多初写多目标调度的同学都会犯这个错导致耗时直接翻三倍。5.3 参数灵敏度分析与常见陷阱种群规模、迭代次数、交叉变异概率这三个参数对结果的影响不是线性的。我建议用简单的控制变量法做灵敏度分析。以变异概率 pm 为例它太小会导致陷在局部最优太大则会像随机搜索破坏已经形成的优良片段。在 FJSP 上MS 段的变异概率和 OS 段的变异概率最好分开设置因为 MS 段的变异相当于重新分配一台机器对目标值影响非常大OS 段的变异只是交换顺序影响相对平缓。我在实验里喜欢设 MS 段 pm0.15、OS 段 pm0.08。如果你全部用统一的 0.1MS 段的搜索步长可能偏小OS 段的扰动又可能偏大。这也算是一个论文里常被忽略、但实际执行效果差很多的细节。还有一个常见陷阱是评价指标不会归一化。如果你把量纲完全不同的三个目标直接做拥挤度距离计算量纲大的目标会彻底支配量纲小的目标。必须先把每个目标在做 Pareto 支配判断和拥挤度计算之前做归一化。很多复现代码把这一步漏掉导致结果失真。5.4 算法不收敛或种群多样性丢失怎么办如果跑了很多代之后 HV 几乎不动优先检查是不是已经收敛到局部区域。这时候我通常先增大变异概率尤其是 MS 段的变异概率让某些个体从当前机器分配组合里跳出来。如果增大变异后 HV 均值反弹明显说明原配置确实陷入停滞。如果种群多样性丢失表现是 Pareto 解集大量重叠在同一小片区域。NSGA-II 场景下可以增大拥挤度选择的权重或者改成基于非支配排序和网格的自适应选择MOEA/D 场景下可以减少邻居数 T让子问题之间的信息交流变少保持各自独立搜索的方向。还有一个小技巧就是在每代结束时用 5%~10% 的新随机个体替换掉种群中拥挤度最低的个体。这种做法称为“重启机制”对调度的离散搜索空间非常有效。它不会破坏主要搜索方向只是定期引入新鲜血液。最后的一点个人体会这套用 MOEA/D 与 NSGA-II 求解柔性车间调度的 Python 代码我前前后后打磨过好几版。最初版本跑出来的 Pareto 前沿非常稀碎解码还有边界 bug。后来逐步把编码约束、解码策略、算子细节都调稳之后整套代码才真正能用在实验和项目里。如果你也想复现我的建议是不要直接拿别人完整代码跑完就结束而是自己手动实现一遍解码函数至少要把插入式左移逻辑写通。因为后面不管换哪个优化算法解码始终是连接问题域和算法域的核心桥梁。你也可以先用单目标遗传算法把调度部分跑通再加入第二个目标、第三个目标最后再切换到多目标算法框架这样排查起问题来会轻松很多。最后再分享一个小技巧在做算法对比实验时用相同的初始种群。也就是让 NSGA-II 和 MOEA/D 从同一个随机种子产生的初始解出发这样两者差异就是算法机制带来的而不是初始化随机性的影响。这个小细节能让你的实验严谨性提升一个档次。