随着300 mm晶圆加工技术的普及，单作业多订单(multiple orders per job,MOJ)调度问题^[1]作为半导体晶圆生产线上常见的NP-hard问题^[2]，近年来受到了业界和学术界越来越多的关注.

为有效解决单机MOJ调度问题，学者们已提出了多种调度方法.Qu等^[1]对单一品种、不考虑到达时间等因素的理想化单机MOJ调度模型进行了基于禁忌搜索算法的调度方法研究.Sobeyko等^[2]对理想化单机MOJ调度模型，构造了分组遗传算法和随机键遗传算法.Jampani等^[3]对带准备时间的单机MOJ调度问题，构造了列生成启发式算法.上述文献提及的调度方法，探讨的问题规模均不超过100个客户订单，而Mason等^[4] 提出的启发式算法对订单规模数为100以上甚至达到10 000的超大规模单机MOJ调度问题进行了有效求解.

在单机MOJ调度研究领域，现有Erramilli等^{[5, 6]}在批调度研究中考虑了产品的多样性.但批调度的研究重点是对由不同作业形成的批次的调度，与本文探讨的调度问题有较大区别.

Chen^[7]对考虑多品种、与次序相关的换模时间的流水线车间MOJ调度问题，建立了混合0-1整数规划模型，但未考虑设备的衰退效应.对带衰退效应的MOJ调度问题进行研究，具有普遍性和实用性，而目前鲜有文献探讨衰退效应背景下的MOJ调度方法.

本文研究的单机MOJ调度问题以最小化总加权提前/拖期惩罚成本为目标函数，同时考虑订单品种的多样性,p-s-d(past-sequence-dependent)换模时间以及设备的衰退效应.

单机MOJ调度问题包含订单成组形成作业和作业排序两个子问题，现举例说明，见图 1.图 1中有包含两种晶圆①和②的订单O₁,O₂,…,O₅和承载上限值为3片晶圆的作业J₁,J₂,J₃.首先进行订单分组形成作业J₁={O₁}，J₂={O₂,O₃,O₄}，J₃={O₅}；然后对分配好订单的作业J₁,J₂,J₃进行排序，图 1中假设进入机器加工的作业顺序为J₂,J₃,J₁.

图 1(Figure 1)

Fig. 1

图 1 单机MOJ调度问题实例Fig. 1 Instance for single machine MOJ scheduling problem

为有效表述单机MOJ调度问题，在文献[7]的基础上，现做如下假设：①调度过程中，订单需保持完整；②每个订单只含有一种产品，不同订单的产品种类可以不同；③每个作业只含有产品种类相同的订单；④每个作业中载有的晶圆数不能超过其承载上限K；⑤作业的数量能够满足调度要求且参与调度的作业能被有效利用；⑥作业加工时间随时间的恶化过程为非定态伽马过程；⑦机器在加工不同品种产品时考虑p-s-d换模时间；⑧机器一旦开始工作，除换模外，不会中断；⑨同一个作业中的订单有相同的完工时间.

对相关符号定义如下，以更清晰地描述调度问题：

M为一个较大的正数；N为订单数；K为晶圆的承载上限值；φ为产品种类数；Φ为作业(基因或节点)个数；λ_il为如果订单i的类型属于产品l，则其值为1，反之为0；f_i为订单i的产品种类标号；n_i为订单i中含有晶圆的数量；d_i为订单i的交货期；w_i为订单i的权重；ρ^f_i为批加工方式下订单i单位晶圆所需的加工时间；a为单位提前时间的惩罚金额 ；b为单位延迟时间的惩罚金额(b>a)；γ为p-s-d换模时间的规模系数；β为非定态伽马过程的尺度参数；α为非定态伽马过程的形状参数；q为非定态伽马过程的幂指数；p_j为作业j的实际加工时间；σ为在作业j之前加工的作业的序号；h_j为加工作业j时所需的换模时间；e_j为作业j的恶化加工时间；C_i为订单i的完工时间；s_j为作业j的开始加工时间；X_ij为如果订单i被分配到了作业j，则其值为1，反之为0；U_j为如果有订单分配给作业j，则其值为1，反之为0；v_jj′为如果作业j和作业j′产品种类相同，则其值为0，否则为1；E_i为订单i的提前完工时间；T_i为订单i的拖期完工时间.

由假设①可知，一个客户订单对应只能分配到一个作业中，即需满足式(1)：

由假设②和③可知，只有产品相同的订单才能分配到同一个作业中，即需同时满足式(2)~式(4)：

由假设④易知，每个作业中所有订单的晶圆总数不能超过K，即需满足式(5)：

由假设⑤可知，参与调度的每一个作业，其含有的晶圆个数不能为0，即需同时满足式(6),式(8)：

在假设⑥和文献[8]的基础上，得到作业恶化加工时间的计算需满足式(9)~式(11)：

在假设⑦和文献[9]的基础上可知，换模时间需满足式(12)，式(13)：

在假设⑧的基础上可以得到批加工方式下每个作业的实际加工时间，计算公式为

由假设⑨知，订单完工时间的计算需满足式(15)：

订单的提前期或者拖延期由其完工时间决定，计算公式分别为

调度目标为总加权提前/拖期惩罚成本最小，即

因此，研究的单机MOJ调度问题是以式(18)为目标，以式(1)~式(17)为约束的非线性规划问题.

将遗传算法融合到动态自适应调整信息素的蚁群算法的迭代过程中，为改善算法的寻优性能，构造ATCS (apparent tardiness cost with setups)准则进行检验和修正，从而建立改进型遗传蚁群算法(modified genetic algorithm-ant colony optimization algorithm,MGA-ACO).为方便描述，现做如下定义.

定义1 O_j=(X_1j,2X_2j,…,NX_Nj),j∈{1,2,…,Φ}表示由分配到作业j的订单序号组成的一维数组.

定义2 蚂蚁的搜索节点或者遗传操作中的基因π_j,j∈{1,2,…,Φ}指调度过程中形成的可行作业.

定义3 蚂蚁搜索得到的一条运动路径或者遗传操作中的一条染色体即为一个可行调度方案，用禁忌表Tabu_m,m∈{1,2,…,Φ}记录.

定义4 g(Tabu_m),m∈{1,2,…,Φ}为运动路径(染色体)Tabu_m的目标函数值，其中

算法的核心操作主要为编码、初始化算子、遗传操作和修正算子.具体步骤如下.

步骤1 编码：分离数学模型中的决策变量，采用双层嵌套编码机制，下层的节点(基因)编码由约束条件得到初始解，并由遗传算法的变异算子搜索，上层的作业加工序号主要由蚁群算法和ATCS修正准则编码.上、下两层编码通过作业序号形成映射关系.为满足订单、作业序号的标记需求，简便算法操作，采用自然数编码.

步骤1.1 节点(基因)编码：依据订单分配原则(见式(1)~式(8))可得蚂蚁搜索节点或遗传基因π_j=(j,f_j,O_j).节点(基因)集合为Γ={π_j|j=1,2,…,Φ}.图 2为对图 1实例的编码示例.

图 2(Figure 2)

Fig. 2

图 2 节点(基因)编码示例Fig. 2 Example for coding node (gene)

步骤1.2 路径(或染色体)编码:将Φ只蚂蚁随机放置在Φ个搜索节点上，蚂蚁m,m∈{1,2,…,Φ}在时刻t从节点j转移到节点j′(j′≠j)的概率为

在计算式(9)~式(17),式(20),式(21)的基础上，蚂蚁m遍历完所有节点后得到运动路径(染色体)Tabu_m={Tabu_m[j]|j=1,2,…,Φ}，其中Tabu_m[j]表示第j个被记录下的节点(基因)序号.图 3为与图 1、图 2对应的一个路径(染色体)编码示例Tabu_m.

图 3(Figure 3)

Fig. 3

图 3 路径(染色体)编码示例Fig. 3 Example for coding route (chromosome)

通过图 2和图 3的映射关系，得到蚂蚁m搜索得到的可行调度方案为{π₂,π₃,π₁}.

双层嵌套编码机制，相比近年由Sobeyko等^[2]和Qu等^[10] 提出的将作业与订单信息合并在同一层的繁冗编码方式更清晰、简便.

步骤2 外部循环终止判断:如果蚂蚁搜索的循环次数θ₁达到上限值θ，则输出结果，结束运算，否则，进入步骤3.

步骤3 初始化算子:Φ只蚂蚁通过步骤1.2遍历完Φ个节点后，得到初始种群population={Tabu_m|m=1,2,…,Φ}.

步骤4 适应度评估:基于式(18)和式(19)，得到适应度函数:

由式(9)~式(17)、式(19)和式(22)计算出初始种群中每一个个体的适应度.

步骤5 内部迭代终止判断:如果遗传操作的迭代次数ε₁达到规定上限ε，则各路径进行如下信息素更新:

步骤6 遗传操作.

步骤6.1 选择算子：依据个体适应度，采用轮盘赌法选择出父本.

步骤6.2 多点交叉：对每一对交叉父本Tabu_m和Tabu_m′(m′≠m)采用多点交叉，生成新的子代.

步骤6.3 变异算子

步骤6.3.1 随机选择变异基因π_j,j∈{1,2,…,Φ}，此时选择次数ζ记为1.

步骤6.3.2 搜索品种标号同为f_j的其他基因.

步骤6.3.3 若有符合步骤6.3.2的基因π_j′,j′≠j，则随机互换O_j和O_j′中的部分值.如果得到的新基因满足假设④或者循环次数ζ达到上限值δ，则进入步骤7，否则，进入步骤6.3.4.

步骤6.3.4 随机选择另一个变异基因π_j″,j″∈{1,2,…,Φ}且j″≠j，令ζ=ζ+1，回到步骤6.3.3.

步骤7 修正算子：在前人研究的ATCS准则^[11]的基础上，针对本文的问题域，构造了在时刻t，未确定加工位置的作业j,j∈{1,2,…,Φ}所具有的紧急度指数I_j(t)，批加工方式下的计算公式见式(25)，设k₁,k₂分别为2和1.

对种群中的个体以一定概率进行如下ATCS修正:按紧急度指数(I_max(t₁),I_max(t₂),…,I_max(t_Φ))排列对Tabu_m,m∈{1,2,…,Φ}进行修正，其中I_max(t_n),n∈{1,2,…,Φ}表示t_n时刻即Tabu_m中的作业Tabu_m[n]开始加工时刻，所有未确定加工位置的作业中的最大紧急度指数.

计算新种群的适应度，若子代的适应度更高，则接纳子代，否则用父代取代子代.经过ATCS修正形成新种群后，令ε₁=ε₁+1，再回到步骤5.

为有效评价算法的性能，以求解时间和解的优度 (performance ratio,PR)^[10]作为算法评价指标.

PR=V(H,T,Y)Best(Y).其中：V(H,T,Y)表示用算法H在迭代次数为T时对实例Y的求解结果；Best(Y)表示对实例Y求解所能得到的参照最优解.PR值越小，表明算法的收敛性能越好.

本文采用C++编程语言实现算法，在主频为2.2 GHz，内存为2 GB的PC机上进行仿真实验.

实验参数参照文献[2, 8, 9, 10, 11]中的参数设计规则设定.在批加工方式下随机生成了10个实例，并分别用MGA-ACO，GA-ACO和ACO三种算法对实例进行求解，得到迭代次数对PR值和求解时间(文中的PR值和求解时间均为平均值)的影响如图 4所示.

图 4(Figure 4)

Fig. 4

图 4 迭代次数对PR值和求解时间的影响Fig. 4 Effect of iterations on PR and computation time

由图 4可知，用MGA-ACO求解实例所得的PR值明显优于GA-ACO和ACO的结果，而三者的求解时间非常相近，因此由算法评价指标可得，MGA-ACO的求解性能明显优于GA-ACO和ACO.当迭代次数由600增加到1 000的过程中，MGA-ACO得到的PR值几乎没有改善，但求解时间成本几乎呈线性增加.考虑保证算法的优良求解性能，在后续实验中，将迭代次数选择为600次.

算法的求解时间主要由它的时间复杂度决定，算法的时间复杂度越高，其求解时间越长.通过计算得到MGA-ACO的时间复杂度为O(φ⁴+θ×ε×(N+φ))，由此易知，算法的运算时间与调度模型中的晶圆品种数呈正相关关系，而衰退效应和p-s-d换模时间的特征参数值的取值变化不会给算法的运算时间带来影响.

对产品种类数φ=5,6,7,8,9,10，而对应的订单数分别为中小规模的20，30和50及大规模的100，150和200的调度问题进行仿真实验.在批加工方式下的实验结果见图 5.

图 5(Figure 5)

Fig. 5

图 5 品种数对PR值的影响Fig. 5 Effect of wafers variety on PR value

由图 5易知，当订单数相同时，随着品种数的增加，PR值会随之减小，这与随着品种数增加，订单分配的限制也会增加，调度组合减少，算法更容易求得优秀调度结果相符；当晶圆品种一定时，随着订单数规模的变大，PR值会随之增大，这与随着订单规模的增大，问题复杂性增强，不易于算法的求解一致.

由实验易得，改变衰退效应和p-s-d换模时间的特征参数值后得到的实验结果与图 5的结果几乎一样.

1) 采用双层嵌套编码机制，以作业序号为纽带及上、下层次间编码信息映射统一的编码方式，为MOJ调度问题传统地将所有编码信息繁冗集中在同一层次提供了新的研究思路.

2) 将遗传算法融合到动态自适应调整信息素的蚁群算法的迭代过程中，同时构造ATCS修正准则，进一步改善算法的收敛性能，建立基于MGA-ACO的MOJ调度方法，丰富了解决此类调度问题的理论方法.

3) 仿真实验表明，MGA-ACO对于模型中特征参数的取值变化具有良好的适应性，算法求解结果稳定且令人满意，验证了MGA-ACO的可行性和有效性.