摘要
第3章晶圆制造并行机调度问题建模
3.1引言
晶圆制造系统是典型的离散时间动态系统，存在各种各样的不确定性。设备情况变化、客户需求变化、加工瓶颈变化、新产品试制、临时工艺更改、返工等都会影响生产线正常运行，造成性能指标波动。晶圆制造过程集成混合了各种生产调度模型，大大增加了生产线复杂性和调度难度。为了设计有效的晶圆制造系统并行机调度方法，首先需要建立合理的调度模型对不同设备调度的本质特点进行分析。
具体来说，该模型应该具有以下功能： ①对并行设备的内在运行变化过程和机理进行分析，确定系统的静态和动态因素对系统运行性能的影响，进而为建立有效的调度方法提供基础； ②确定评价调度性能的评价指标，进而对调度方法的性能进行快速有效的评价，以验证调度方法的有效性。
目前，晶圆制造系统中加工系统并行机调度问题的建模方法主要有数学规划模型和排队论模型； 物流运输系统并行机调度问题的建模方法主要有复杂网络模型、马尔可夫模型及仿真模型等。
3.2晶圆制造加工系统建模
3.2.1基于数学规划模型的晶圆制造加工系统建模
3.2.1.1基本概念
数学规划是运筹学的一个重要分支，也是现代数学的一门重要学科。其基本思想出现在19世纪初，并由美国哈佛大学的Robert Dorfman于20世纪40年代末提出。数学规划的研究对象是数值很优化问题，这是一类古老的数学问题。古典的微分法已可以用来解决某些简单的非线性很优化问题。直到20世纪40年代以后，由于大量实际问题的需要和电子计算机的高速发展，数学规划才得以迅速发展起来，并成为一门十分活跃的新兴学科。今天，数学规划的应用极为普遍，它的理论和方法已经渗透到自然科学、社会科学和工程技术中。
数学规划问题实际上都是在满足某些约束条件的情况下求函数极值的问题。一个数学规划问题通常包含以下要素： 一组决策变量，变量的值就是要求解的对象； 一个目标函数，就是我们希望通过改变决策变量而优选化或最小化的对象； 一组约束条件，决策变量的取值必须满足这些约束条件。通常满足所有约束条件的解被称作可行解，所有可行解的集合叫作可行域。
数学规划模型一般可表示为
min(max)f(x,α,β)
s.t.g(x,α,β)≤0(3？1) 其中，f表示数学规划问题的目标函数，g为数学规划问题的约束函数，x为数学规划问题中的变量，α表示模型中的已知参数，β为模型中的随机参数。
数学规划主要包括线性规划、整数规划、混合整数规划以及动态规划等，按照很优化理论可将数学规划问题划分为无约束优化问题和约束优化问题两大类，可通过很优化理论和算法对数学规划模型进行化简和求解。
3.2.1.2数学规划模型在晶圆制造并行机调度中的应用
晶圆制造系统中的加工设备划分为： 单件加工型设备、批处理型设备等。Kumar［43］等学者曾利用线性规划的方法解决晶圆制造并行机调度问题，给出了生产过程中保证生产稳定进行的WIP水平界限和开环带重入生产线最小延迟的界限，对实际生产起到一定的指导作用。但总体而言，由于调度问题的复杂度为NP？难，所以数学规划方法的应用受到了。为此，提出了分支定界法、拉格朗日松弛法以及过滤束搜索法等方法以降低问题的难度。Sung和Choung［55］曾使用分支定界法解决半导体生产工厂中的一类批处理设备――炉管的动态调度问题，其优化目标是优选化所有工件的完成时间(Makespan)。Kaskavelies和Caramanis［56,57］使用拉格朗日松弛算法对超过10000个资源约束的半导体测试工厂进行调度。他们介绍了两类拉格朗日乘子更新过程的新特征，以帮助运算的进行。De和Lee［58］使用过滤束搜索方法，构建一个基于知识的调度系统，以调度半导体测试生产线。其调度目标是最小化总生产周期。在调度过程的搜索树上，所需要搜索的节点的数量取决于过滤宽度和束的宽度。
批处理机一般为瓶颈机，它制约着晶圆制造系统的绩效。这里介绍数学规划模型在具有重入特性的批处理机重组批调度问题中的应用，具体描述如下。
图3？1具有重入特性的批处理组批加工工艺流程简图
加工工艺流程简图(图3？1)主要包括两个设备组，设备组1和设备组2，且这两个设备组均是批处理机，设备组2的优选批处理能力大于设备组1； 设备组的产品组批包括两种类型： 产品相容和不相容： 设备组1具有产品不相容性组批、重入特性； 设备组2具有产品相容性组批，但由设备组1提供的批不拆分，具有优选批处理能力，需要对进入缓冲器中产品批进行重组批、排序［59］。主要假设如下：
(1) 设备组2由于加工时间长，故此工位是生产瓶颈，假设设备组2不会出现饥饿。 (2) 抢占不允许，即任意设备一旦开始加工，加工过程就不允许中断。
(3) 缓冲器里产品批到达是动态的。
(4) 设备组2的缓冲器里产品批存在等待时间，即为了避免由于工件表面接触空气时间过长而对产品的质量造成影响，产品批在设备组1中加工完成后必须在规定的时间内进入设备组2中加工。
(5) 三层混合调度算法中第一层完成缓冲器里产品批实时状态信息和设备组2实时状态信息的记录、保存和传送； 第二层完成进入缓冲器中产品重组批、排序； 第三层完成重组优先批的装载分配及其分配完毕后缓冲器里产品批实时状态信息和设备组2的实时状态信息更新、保存。
(6) 不计工件从设备1组到缓冲器，从缓冲器到设备组2和在缓冲器重组批的过程时间。
为了使术语具有通用性，对后面需要的术语进行统一如下定义。
K： 当设备组2出现设备空闲可用时，缓冲器里来自设备组1的不同批的数量集；
K′： 当设备组2出现设备空闲可用时，重组批的数量集；
S： 设备组2的优选加工能力；
M： 一个足够大的正整数；
t： 设备组2出现设备空闲可用时的当前时间；
qk： 批k在设备组1和设备组2之间等待时间；
P： 设备组2的加工时间；
sk： 批k包含的工件数量；
rk： 批k到达缓冲器的时间；
dk： 批k的交货时间；
ZSk： 批k总加工工序数量集；
DSk： 批k在设备组2上所属的加工工序数；
Pwk： 批k在第w个加工工序中的加工时间；
PRk： 批k在设备组2后纯剩余的加工时间；
drk′： 重组批k′的交货时间；
Crk′： 重组批k′的完工时间；
PRrk′： 重组批k′在设备组2后纯剩余加工时间；
∑T： 总拖延时间和。
注： k=1,2,…,K； k′=1,2，…,K′； w=1,2,…,ZSk。
以拖延时间和最小为目标，且具有等待时间和工件动态到达的批处理机重组批的NP？难调度问题，在分解规则下，将总的调度时间窗分解为许多时间窗，每个时间窗的长度等于相邻两个触发事件之间的时间间隔，其值是变化的，即变时间窗； 在滚动时域策略下，每个时间窗之间不断更新优化调度结果。
每个时间窗内的子问题，由于规模被减小，准确算法可以适用。在整数规划模型求解重组批及其排序优化调度问题中，重组批状态的0？1变量和重组批的排序位置序号是决策变量； 满足设备加工时间，同一设备不能同时加工两个重组批，同一重组批只能占有一个排序位置序号等是约束条件； 总拖延时间和最小为目标函数，具体整数规划模型如下：
目标函数
∑T=∑K′k′=1max［(Crk′-drk′),0］(3？2)
约束条件
Crk′=t+Position(k′)×P+PRrk′,？k′∈{1,2, …,K′}(3？3)
Position(k′)∈{1,2,…,K′},？k′∈{1,2,…,K′}(3？4)
如果k′≠l′，那么Position(k′)≠Position(l′),
？k′和l′∈{1,2,…,K′}(3？5)
xkk′=1，当缓冲器的批k被分配到重组批k′中
0,其他，
？k′∈{1,2,…,K′},？k∈{1,2,…,K}(3？6)
∑K′k′=1xkk′=1,？k∈{1,2,…,K}(3？7)
PRk=μw∑ZSkw=DSk+1Pwk,？k∈{1,2,…,K}(3？8)
μw=sk，如果第w的工序数不是批处理设备
skBw，如果第w的工序数是批处理设备，且其加工能力是Bw(3？9)
K′=
∑Kk=1skS
(3？10)
∑Kk=1xkk′？sk≤S,？k′∈{1,2,…,K′}(3？11)
PRrk′=min1≤k≤K(xkk′×PRk),？k′∈{1,2,…,K′}(3？12)
drk′=min1≤k≤K(xkk′×dk),？k′∈{12,…,K′}(3？13)
t-rk+xkk′×Position(k′)×P≤qk,
？k′∈{1,2,…,K′}，？k∈{1,2,…,K}(3？14)
式(3？2)是每个时间窗内的总拖延时间和最小的目标函数。约束条件式(3？3)用于计算重组批k′的完工时间。约束条件式(3？4)和式(3？5)用于确保每个重组批只占有一个排序位置号，且每个排序位置号只有一个重组批。约束条件式(3？6)引入重组批状态的0？1变量。约束条件式(3？7)确保缓冲器中每个批都被进行重组批，且只被进行1次重组批。约束条件式(3？8)用于计算批k在设备组2后剩余的加工时间，其中μw是调整系数。约束条件式(3？9)用于确定调整系数。约束条件式(3？10)用于计算当设备组2出现设备空闲可用时，重组批的数量集K′。约束条件式(3？11)用于确保任一重组批包含工件数量不大于设备组2的优选加工能力S。约束条件式(3？12)用于计算重组批k′在设备组2后剩余的加工时间，其值等于重组批中所包含产品族中在设备组2后剩余的加工时间的最小值。约束条件式(3？13)用于计算重组批k′的交货时间，其等于重组批中所包含产品族中交货期的最小值。约束条件式(3？14)用于确保被调度批k在缓冲器里的等待时间不超过设备组1和设备组2之间等待时间qk。
注： 约束条件式中
a
代表返回到不小于a的最小整数，例如，
(∑Kk=1sk)S
表示返回不小于∑Kk=1sk除以S的最小整数，当∑Kk=1sk=51和S=12，那么
∑Kk=1skS
=
51/12
=
3.5
=4。(Position(1),Position(2),…,Position(K′))是优化排序后位置序号的一维数组。
在CPLEX和开发程序联合平台中，线性的约束适合CPLEX求解。约束条件式(3？12)、式(3？13)和式(3？14)中表达关系式是非线性的，需要进行线性化。
Montagne介绍批调度时考虑用平均加工时间p-=∑ni=1pi
计算比率系数pin×p--di处理总拖延时间之和最小问题的Montagne比率方法［60］。这里基于Montagne比率方法，计算每个重组批内所有待调度批的平均加工时间，分别建立约束条件式(3？12)和式(3？13)的松弛线性约束条件式(3？15)和式(3？16)。
PRrk′=∑Kk=1(xkk′×PRk)∑Kk=1xkk′,？k′∈{1,2,…,K′}(3？15)
drk′=∑Kk=1(xkk′×dk)∑Kk=1xkk′,？k′∈{1,2,…,K′}(3？16)
采用大M法将约束条件式(3？14)线性化：
t-rk+Position(k′)×P-M(1-xkk′)≤qk,
？k′∈{1,2,…,K′},？k∈{1,2,…,K}(3？17)
综上，等式(3？2)、约束条件式(3？3)~式(3？11)和式(3？15)~式(3？17)构建基于松弛的混合整数线性规划模型，适合CPLEX和开发程序联合平台优化求解。实验结果表明，该模型能在CPLEX环境下以较短CPU运行时间取得令人满意的目标函数值。
3.2.2基于排队论模型的晶圆加工系统建模
3.2.2.1基本概念
排队论，或称随机服务系统理论，是通过对服务对象到来及服务时间的统计研究，得出这些数量指标(等待时间、排队长度、忙期长短等)的统计规律，然后根据这些规律来改进服务系统的结构或重新组织被服务对象，使得服务系统既能满足服务对象的需要，又能使机构的费用最经济或某些指标很优。它是数学运筹学的分支学科，也是研究服务系统中排队现象随机规律的学科。所谓排队，指在服务系统中要求服务的“顾客”所形成的等待线。其主要研究方向是如何合理地设计与控制服务系统，使之以最经济的服务机构满足顾客的需要。由于服务系统中顾客的到达及服务时间和次序都是随机的，因此，排队论又被称为随机服务系统理论。
D.G.Kendall首先用3个字母组成的符号A/B/C表示排队系统，其中： A表示顾客到达时间分布； B表示服务时间的分布； C表示服务机构中的服务台的个数。在此基础上，L.Takacs等人又将组合方法引进排队论，使它更能适应各种类型的排队问题。
排队论在晶圆制造中的研究在1950年也已经开始，用于解决在制品排队问题的方法层出不穷，从大的方向上来看有两个显著的研究趋势： ①基于运筹学的排队方法，强调数学规划的方法来建立排队模式，或利用启发式搜索法来求得很好或近似很好的排队方法，这里面包括派工法、线性规划、动态规划、分支界限法、启发式搜索等方法； ②基于知识的方法，即通过模拟人类解决问题的思考方法或人类某部分器官的功能来解决复杂度高的排队问题，这里面包括智能调度专家系统、神经网络、基因演算法、条件满足技巧、基于约束的方法以及一些基于智能搜索的算法。
3.2.2.2排队论模型在晶圆制造并行机调度中的应用
排队论模型能直观描述生产线的加工过程，因此，很多生产系统都用此方法建模。Kumar［61］将晶圆制造系统抽象成可重入的排队网络，Meyn等［62］用马尔可夫模型对可重入排队网络进行建模，对模型研究了相应的动态规划算法。黄敏等［63］提出了排队论模型和非线性整数规划模型相结合的方法，并通过遗传算法对模型求解。朱连章等［64］在基本排队论模型的基础上加入Monitor监视器监视网络状态，在仿真过程中收集状态数据进行性能评价，该方法在分析过程中不会对系统结构进行，具有很大的灵活性。
现代晶圆制造业正在朝着更高密度的集约化生产和更大规模的量产方向发展。在现代的半导体晶圆制造厂中，经常能看到晶圆在进行某些工艺流程的加工制造中，通常会有多台晶圆制造设备来对应同一工艺和工序。于是，这就给晶圆制造中的排队策略带来了新的问题： 如何达到晶圆制造并行设备的很优控制――是采用多制造设备、多队列的排队系统方案，还是应该采用多制造设备、单队列的排队系统方案。为了解决这个问题，通过建立符合实际情况的数学模型来求解［65］。
首先考虑当某个批次的晶圆加工到某一工艺和工序时，有M台晶圆制造设备正在独立的并行制造，这里只考虑晶圆制造设备状态处于正常和待机两种情况下，即有M台晶圆制造设备处于正常或待机。当晶圆到达时，若有处于待机状态的晶圆制造设备便立刻开始加工制造，若当前所有的晶圆制造设备都处于正常状态，则这个批次的晶圆进入排队等待系统，直到有处于待机状态的晶圆制造设备时再开始加工制造。
其次，为了使建立的数学模型更具有随机性和典型性，这里假设晶圆批次按照参数为λ(λ>0)的泊松流到达，每个批次晶圆需要待加工的制造时间是独立的，且服从参数为μ(μ>0)的负指数分布，不考虑排队等待系统的容量，而且到达与加工制造是彼此独立的。假设M(t)表示某时刻t，M台晶圆制造设备中的晶圆批次数，令
pk1k2(Δt)=P{M(t+Δt)=k2|M(t)=k1}(3？18)
易得
pk1k2(Δt)=PλΔt+o(Δt),k2=k1+1,k1≥0
k1μΔt+o(Δt),k2=k1-1,k1=1,2,…,M-1
k1μΔt+o(Δt),k2=k1-1,k1=M,M+1,…
o(Δt),|k1-k2|≥2
(3？19)
于是{M(t),t≥0}为E={0,1,2，…}上的生灭过程，其中：
λk2=λ,k2≥0;μk2=k2μ,1 Mμ,k2≥M(3？20)
再令ρ=λμ,ρM=λMμ，pk2=limt→∞P{M(t)=k2},k2=0,1,2,…，则当ρM<1，有{pk2,k2≥0}存在，与初始条件无关，且
pk2=1k2!ρk2p0,1≤k2≤M-1
1Mk2-M？M!ρk2p0,k2≥M，p0=∑M-1k2=0ρk2k2!+MρMM!(M-ρ)-1
(3？21)
因为有M台晶圆制造设备，所有晶圆批次到达时需要等待的概率为
p=∑∞k2=MPk2=11-ρMpM(3？22)
其中，ρM=λMμ<1， pM=ρMM!p0。
为了求得平均队长与平均等待队长，在统计平衡下，等待队长Mq呈现分布状态，为
P{Mq=0}=∑Mk2=0Pk2,P{Mq=μ}=pM+μ,μ=1,2,…
故当ρ1<1时，有以下关系存在：
Mq=EMq=∑∞k2=M(k2-M)pk2=∑∞k2=M(k2-M)1Mk2-M？M!ρk2p0
=ρMp0M!∑∞k2=M(k2-M)ρMk2-M(3？23)
经过方程式的进一步化解，可得
Mc=ρ-ρM(M-1)!p0-ρM+1M!(1-ρM)p0+ρM(M-1)!(1-ρM)p0
(3？24)
不难发现，上面的表达式即为ρ的泰勒二次展开式，所以Mc=ρ。
从而得出，当排队等待系统平衡时，正在加工的晶圆批次数与晶圆制造设备的个数(M)无关。
由式(3？23)和式（3？24）可得：
在某时刻t，平均队长M=Mq+Mc=ρ+ρM(1-ρM)2ρM。
其次，假设晶圆批次数是按照FIFO(优选先出)的规则进行加工制造，令p-k2表示晶圆批次数到达排队等待系统后发现已有k2个晶圆批次的平稳概率，则可以推出：
p-k2=pk2,k2=0,1,2,…，当ρ=λMμ<1，在统计平衡下，晶圆批次的排队等待时间可以用分布函数Wq(t)=1-pM1-pMe-μ(M-ρ)t来表示，则其平均等待时间为
Wq(t)=ρMλ(1-ρM)2pM(3？25)
在到达的晶圆批次看到已有k2(k2>M)个晶圆批次的条件下，由于M台晶圆制造设备均处于正常状态，所以该晶圆批次必须等待(k2-M+1)个晶圆批次加工完成后才能被加工。在M台晶圆制造设备均处于正常状态的条件下，由于每个晶圆批次的离去均是参数为μ的泊松流，这样相继离去晶圆批次的间隔时间服从参数为Mμ的负指数分布，所以该晶圆批次的排队等待时间等于(k2-M+1)个晶圆批次相继离去的间隔时间总和，其分布为参数Mμ的(k2-M+1)阶埃尔朗分布，即P{0 Wq(t)=W0(t)+∑∞k2=Mp0ρk2M!Mk2-M∫t0Mμ(Mμx)k2-M(k2-M)!e-Mμxdx
=1-pM1-ρM+pM1-ρM1-e-(1-ρx)Mμt(3？26)
经过简化后可得排队等待时间为
Wq(t)=1-pM1-ρMe-μ(M-ρ)t(3？27)
从而平均排队等待时间为： Wq(t)=∫∞0tdWq(t)=ρMλ(1-ρM)2pM。
由于逗留时间W=Wq+ χ，且Wq与 χ独立，因此根据分布函数的卷积公式可以得出逗留时间的分布函数为
W(t)=1-e-μt1+pM1-ρMμt,ρ=M-1
1-e-μt-pM(1-ρM)(M-ρM-1)［e-μt-e-μ(M-ρ)t］,ρ≠M-1
(3？28)
？= Wq+1μ(3？29)
根据上面公式推导的定性分析，可以通过一个实际的举例来进行定量判断。
在某个批次的晶圆加工到某一工艺和工序时，有两台晶圆制造设备正在独立地并行制造，且批次的晶圆按照参数λ=8批次/min的泊松流到达，而每台晶圆制造设备的加工制造时间都服从μ=5批次/min的负指数分布。可以比较两种排队方案的运行指标，并以此判断晶圆制造并行设备的很优控制。
（1） 多制造设备、多队列的排队系统方案。晶圆批次到达后，以50%的概率分成两个队列进行排队等待，如图3？2所示。
图3？2多制造设备、多队列的排队系统方案
在这种情况下，可以将排队系统方案看成两个独立的M/M/1/∞排队系统，每个晶圆批次的到达均为λ1=λ2=4(批次/min)的泊松流，且
ρ=λμ=45？p0=1-ρ=0.2
