在半导体、钢铁、化工等实际生产系统中，在制品(Work In Process, WIP)在进行某些工序加工之前的等待时间往往不能超过设定的阈值，即存在等待时间约束(queue time constraints)。若超出该约束，则产品直接废弃或需要返工，这将降低生产效率。此外，设备劣化程度的加剧将会导致产品缺陷率的增大，从而进一步降低系统的有效产出。如何设置合理的生产控制策略以提高系统产出已成为制造企业面临的一项难题。

考虑实际生产中的等待时间约束正成为学术界和工业界的研究热点。Tu和Chen^[1]将设备的平均失效时间(Mean Time Between Failure, MTBF)和平均修复时间(Mean Time To Repair, MTTR)整合进GI/G/m排队模型，为带等待时间约束的串行工作站作了能力规划，并研究了不同控制策略下各工作站能力规划水平的变化。Ono等^[2]基于M/M/1模型和利特尔法则(Little’s law)，根据缓冲区WIP的阈值大小，建立了一种基于风险的评估方法。更多关于排队模型在这类问题上的应用可参见文献[3]。Wu等^[4]基于马尔可夫决策过程，通过最小化持有成本和废弃成本之和，为一个等待时间约束(含加工时间)呈指数分布的串行系统建立了最优控制策略；之后又将模型扩展至工作站内含多台并行机的情形^[5]。Shi和Gershwin^[6]假设上下游设备的MTBF和MTTR均服从几何分布，基于离散时间马尔可夫链推导出了WIP等待时间的概率分布。文献[7]同样考虑设备的随机故障，借鉴经济订货批量(Economic Order Quantity, EOQ)思想，通过权衡WIP的返工成本和下游瓶颈工作站的能力损失，推导出了缓冲区容量的最优控制限。上述文献从多种角度对带等待时间约束的生产系统进行了分析研究，如能力规划、风险、缓冲区控制等，并且大都考虑了设备的随机故障。但是，仍存在2个问题：一是仅选用MTBF和MTTR 2个参数对设备的劣化过程及修复进行刻画，与实际偏差较大；二是就其维护方式而言，属于事后维修，这一假设在自动化程度较高的行业里并不合理，如半导体制造。

关于生产设备的退化和可靠性问题一直是学术界研究的热点。Peng等^[8]综合考虑系统的连续劣化过程和受到的突发性外界冲击，建立了可靠性和维护模型。王严等^[9]分析保障设备故障对装备保障活动的影响, 建立了保障设备站点内的保障过程模型。张卓琦等^[10]则考虑了系统中部件间存在的故障相关性，提出了一种年龄预防性维修与机会维修相结合的维修策略。上述文献从设备故障的不同内在机理出发，建立故障模型，但并未将设备的退化过程与实际生产活动(如引起的产品质量问题)相联系起来。

综合以上2点，本文针对一个带等待时间约束的生产系统，用更为准确的劣化模型来刻画其劣化过程，并考虑其劣化过程引起的质量损失问题，通过引入预防性维护(Preventive Maintenance, PM)策略和M/G/1/K排队模型对预防性维护阈值和中间缓冲区容量作联合优化，以最大化系统的有效产出。

1 问题描述

考虑如图 1所示的一个由2个工作站串联的带缓冲区的生产系统。其中，上游工作站1为非能力约束资源(Non-Capacity Constraint Resources, NCCR)，下游工作站2为能力约束资源(Capacity Constraint Resources, CCR)，即瓶颈工作站。每一工作站前面都有一定数量的WIP，分别组成队列Q₁、Q₂，由工作站1加工后的WIP在Q₂中的等待时间不得超过某一设定的时间上限，即存在“等待时间约束T_C”。

定义1 有效产出：以图 1所示的生产系统为研究对象，将总期望产出去除质量损失、被阻塞未进入Q₂的潜在工件以及因超出T_C被废弃的WIP后剩余的产出，如图 2所示。

为进一步描述研究问题域，作如下假设：① 工作站1和2一次均只能加工一个工件。工作站1前面的缓冲区容量为无限大，为方便后文计算，工作站2前面的缓冲区容量设为K-1；② 工件到达工作站1的时间间隔服从指数分布，经工作站1加工后进入缓冲区Q₂, 当Q₂已满时，则关停工作站1以阻止新的WIP加入队列，直至Q₂中WIP数量小于K-1；③ 若Q₂中WIP的等待时间超过“等待时间约束T_C”，则直接废弃；④ 随着生产活动的进行，工作站1和2都不可避免地发生劣化。对非瓶颈工作站1采取事后维修(替换)策略，替换后恢复至全新状态，并且替换时间忽略不计；瓶颈工作站2服从连续的劣化过程，并对其状态进行连续监测，在其劣化量达到失效阈值D_F之前对其进行PM，并且每次PM均不能使其恢复至初始全新状态；⑤ 经工作站2加工出的产品并非全都是完美的，存在质量损失，并且质量损失随着其劣化量的增大而增大。

基于以上假设，本文的决策目标是：在瓶颈工作站2的预期寿命范围T_L内，综合决策工作站2的PM阈值D_P和缓冲区Q₂容量的控制上限K-1，以最大化系统的有效产出。

2 模型构建 2.1 生产系统的劣化过程及维护建模

考虑到工作站1为NCCR，仅对到达的工件进行初步的预处理，因此对其进行事后维修(替换)，该工作站内一般均为通用型设备，较易获得，忽略其替换时间。

基于约束理论(Theory of Constraints, TOC)，瓶颈决定系统产能，为此选用更为准确的劣化模型对工作站2的劣化过程进行刻画。伽马过程具有非负的独立增量特性，被广泛应用于模拟单调累积增长的渐变式退化过程，如磨损、疲劳、腐蚀、裂纹生长、蠕变等^[11]。更多关于伽马过程在设备劣化及维护领域的应用可参见文献[12]。文献[11]对金属表面点蚀生长过程的研究表明：点蚀坑深的增量与时间呈非线性关系。为此本文选取非定态伽马过程来对工作站2的连续劣化过程进行建模，假设其形状参数α(t)=αt^q(q>0)，尺度参数β>0。设工作站2在初始时刻为全新状态，即X(0)=0，则其t时刻劣化量X(t)的概率密度函数(PDF)为

式中：为α的伽马函数。设服役时间为T_L，则对于给定的失效阈值D_F，相应的失效时间的分布函数为

式中：Γ(t, x)为上不完全伽马函数，Γ(t, z)= 。

定理1 如果某设备(系统)的劣化过程服从伽马过程，并存在一个失效阈值(上界)，则其失效率函数是关于时间的增函数。

证明 不失一般性，考虑如下伽马过程X(t)：初始状态X(0)=0，形状参数α(t)>0，尺度参数β>0，失效阈值为D_F。根据定义，失效率函数可表示为

式中：f(t)为失效时间的概率密度函数；F(t)为可靠度函数。

令ΔX=X(t+Δt)-X(t)表示在区间[t, t+Δt]内的劣化增量，设其概率密度函数为g(y, Δt), y∈(0, D_F)。于是有

式中：为下不完全伽马函数，且。设

对t求偏导，可得

于是，G(x, y, t)是关于t的增函数，那么 , y∈(0, D_F)是关于t的减函数，而h(t)是关于t的增函数。    证毕

定义2 设X，Y是2个随机变量，若对任意正数a，有P(X>a)≥P(Y>a)，则称X随机的大于Y，记作X≥_stY^[13]。

进一步地，设{X_n, n=1, 2, …}表示某一随机过程，若对所有的n=1, 2, …，均满足X_n≥_st(≤_st)X_n+1，则称{X_n, n=1, 2, …}是随机递减(递增)的。

定理2 设X(M_i)表示第i次PM后，系统劣化量恢复至的中间某一水平。假设对，均有X(M_j)≥_stX(M_i)，并且h_k(t)|_x, k∈{1, 2, …, i, i+1, …, j, j+1, …}是关于x的增函数。那么，在任意时刻t，都有h_j(t)≥h_i(t)，并且期望长度E(T_j+1)≤E(T_i+1)。

证明 由文献[11]可知，若X≤_stY，φ(u)是一个单调增函数并且φ(u)≥0，则有E[φ(X)]≤E[φ(Y)]。根据定理2中的前提条件，易得h_j(t)=E[h_j(t)X(M_j)]≥E[h_i(t)X(M_i)]=h_i(t)其中：h_k(t)|_x显然恒为正数。设F_i+1(t)为系统第i+1个运行周期T_i+1的寿命分布函数，R_i+1(t)为其可靠度函数，则

工作站2从初始全新状态第1次运行达到PM阈值所经历时间的期望值为

第i次PM后(i=1, 2, …，N-1)，第i+1个运行周期的期望长度为

式中：f_X(Mi)为第i次PM后，工作站2的劣化量恢复至X(M_i)(0 < X(M_i) < D_P)的概率密度函数。

假设第i次PM时间M_i(i=1, 2，…, N-1) 是关于维护阈值和维护次数的增函数，M_i服从参数为λ_i的指数分布(其期望 ，γ₁、γ₂为比例系数)，其概率密度函数f_Mi(t)=λ_ie^－λ_it。

由以上分析可知，将工作站2的整个预期寿命T_L划分为N个维护周期。随着N的增大，运行周期的期望长度越来越短，而每个运行周期结束后进行PM所花费的时间越来越长，直至达到预期寿命对其进行更新(替换)。

2.2 质量损失建模

由第1节假设⑤，随着工作站2劣化程度的加剧，其所生产出产品的质量性能也随之下降，即存在质量损失。根据文献[14]，可以用田口质量损失函数(二次损失函数)衡量质量损失的相关成本：

式中：L(X)为质量损失函数；X₀为系统全新状态时的劣化量(本文X₀=0)；C为相关的成本系数。

在工作站2的整个寿命周期内，总质量损失为

为描述产品质量损失与产出损失L(O)间的关系，在文献[15]基础上，构造两者间函数关系如下：

式中：k₁、k₂为比例系数。于是，在工作站2的整个寿命周期内，由质量损失引起的总产出损失为

2.3 基于M/G/1/K排队模型的缓冲区建模

如图 1所示，工件经工作站1加工后进入队列Q₂，按照到达的先后顺序依次在工作站2进行加工，可以视为一排队系统。设α(n)表示在工作站2加工某工件期间，有n个新工件进入Q₂的概率，则

式中：λ为工件的到达率；g(t)为工作站2加工单位产品所需时间s₂的概率密度函数。由第1节假设②，工作站1和2之间的缓冲区最多可容纳K-1个WIP。设缓冲区的状态i表示在工作站2加工完某一工件的瞬间，队列中恰好存在i个WIP(i=0, 1，…, K-1)。p_{i, j}表示图 3所示马尔可夫链从状态到状态的转移概率，p_{i, j}定义为

设p_n(n=0, 1, …, K－1) 表示当系统处于稳态时，在工作站2加工完某工件的瞬间Q₂中恰好有n个工件的稳态概率。p_n可以通过式(10) 求解：

式(10) 中2个公式分别表示状态平衡方程和归一化等式。在工作站2加工完工件的瞬间，系统有0, 1，…, K-1共K种可能状态，即共有K个独立变量。但式(10) 中平衡方程只有K-1个是相互独立的，因此必须加入归一化等式才能求解。选取前K-1个平衡方程以及归一化等式，并改写成以下形式：

设向量p=[p₀, p₁, …, p_K－1]^T，b=[0, 0, …, 1]^T，A为一个K×K的矩阵：

根据式(10) 中K个相互独立的线性方程组，可将问题转化为求解线性方程组Ap=b，其解为p=A^－1b。注意到此时的p_n(n=0, 1，…, K)表示在工作站2加工完某工件的瞬间，Q₂中存在n个WIP的稳态概率，而非任意时刻系统中存在n个WIP的概率。根据泊松过程的PASTA性质，任意时刻系统中恰好有n个WIP的稳态概率p_n(n=0, 1，…, K)可通过式(11) 计算：

式中： 为工作站2加工单位产品所需时间的期望值。工件因超出缓冲区容量限制被阻塞的概率即为

缓冲区WIP超出T_C而废弃的概率为

式中：g_k(t)为工作站2加工完k个WIP所需时间的概率密度函数，可由其加工单个WIP时间的概率密度函数g(t)求得。

3 数值实例

本节通过一个数值实例来对所提模型进行演示。工作站2的预期寿命T_L=650 h，失效阈值D_F=150 μm，工件的到达率服从参数λ=5.5的指数分布，工作站2的加工单位产品的时间s₂服从[0.1, 0.2]的均匀分布，中间缓冲区WIP的等待时间约束T_C=1.2 h。参考文献[11]中对点蚀生长过程所服从非定态伽马过程的参数估计方法，用最大似然法对参数α、β和q进行估计。假设进行了n次不同时段的观测，似然函数为

对式(14) 取对数，得

最大化式(15)，利用MATLAB编程得到各参数的估计值：α=3.2，β=4.3，q=0.5。选取PM阈值D_P=120 μm，工作站2劣化量达到PM阈值所经历时间的累积概率函数(CDF)和概率密度函数如图 4所示。可以看到：在160 h以后，工作站2劣化量达到PM阈值的概率将稳定趋近于1。

根据第2节构建的劣化及维护模型，当D_P=120 μm时，工作站2在前5个维护周期(N=1~5) 内失效率函数的变化如图 5所示，可以得出以下结论：① 在任一维护周期内，其失效率函数是关于时间的增函数；② 随着维护周期的增大，同一时刻工作站2的失效率越来越大。图 5验证了定理2的正确性，并且可以得出推论：随着维护周期个数的增大，工作站2的期望运行时间越来越短。

图 6验证了上述推论的正确性，并且可以看到：当PM阈值D_P分别取为[110,115,120,125,130] μm时，对应工作站2在寿命范围内的维护周期个数分别为^{[12, 11, 10, 10, 9]}，这是由于D_P值设置的越大，则对应的运行时间和PM时间均越长，则在有限的寿命范围内总维护周期个数越少。图中，D_P=125 μm和130 μm对应的2条曲线在最后一个维护周期快速下降的原因是最后一个维护周期是不完整的，分别运行了5.4 h和22.2 h便达到了预期寿命T_L(更新点)。

图 7表明了在理想情况下(即不考虑质量损失、阻塞及WIP等待时间超出T_C)，不同的PM阈值对应的总期望产出及质量损失。可以看到：随着PM阈值的增大，系统在理想情况下的总期望产出也越大，原因在于总运行时间相对变长，总PM时间相对变短；但是引起的质量损失成本也随之增大，这是因为系统在进行PM前的劣化程度加剧，而2.2节的质量损失模型表明：质量损失成本与劣化量的二次方成正比关系。

工件的阻塞概率(P_B)及因超出T_C而废弃的概率(P_S)随缓冲区容量K-1的变化情况如图 8所示。随着K值的增大，工件的阻塞概率变小，但缓冲区内WIP超出T_C约束的概率变大。容易理解：将缓冲区容量设置的越大，则经工作站1加工后的工件有更大的概率加入队列Q₂；但同时意味着Q₂中WIP在被工作站2加工前的等待时间更长，则超出T_C的概率越大。

D_P及K的变化均会影响系统的有效产出，但是难以将两者整合进一个数学模型中，以寻求最大化目标函数对应的最优组合。为此采用数值迭代算法寻找两者的最优组合[D_P^*, K^*]，见表 1。

可以看到：当[D_P, K]取值为[130,7]时，对应的系统有效产出最大，为3412个。

4 结论

1) 定义“有效产出”这一概念，使制造企业关注实际生产过程中由于设备劣化引起的质量损失、缓冲区的阻塞以及由于WIP等待时间过长而废弃等问题，避免盲目追求生产率。

2) 对工作站2全寿命周期内劣化过程及预防性维护的建模，能较好地符合生产实际，为制造企业开展维护活动提供参考。

3) 在结论2) 的基础上考虑由于设备劣化引起的质量损失，进而提出了一种基于状态的预防性维护策略，减少系统因产品质量而损失的产出。

4) 引入M/G/1/K排队模型求解缓冲区WIP的阻塞及超出等待时间约束的概率，为存在等待时间约束的实际生产系统的缓冲区控制提供了一种新的优化思路，进而提高系统的有效产出。

数值实例表明：本文所提模型是切实、有效的，对于制造企业的设备维护、缓冲区控制以及产出提高有一定的指导意义。为使模型更符合实际生产情况，可进一步考虑实际生产成本、维护成本等。