随着制造业的快速发展，展开智能化制造成为制造业发展的历史必然，具有自动装卸搬运机器人的制造单元已经成为高端智能设备不可缺少的一部分。该类制造单元具有如下特点：1) 只有在下一工序的加工中心处于空闲并且机器人没有其他任务时，工件才能由机器人搬到下一工序进行加工，即存在着加工中心和机器人2种资源的同时性约束；2) 该制造单元具有复杂的随机不确定性，如工件的到达时间，工件的加工时间，机器人的移动及装卸速率。

对于双重资源约束(dual resource constrained，DRC)系统，文献[1]从任务投放机制、任务指派规则、工作者柔性、工作者安排、运送成本5个方面对DRC系统的计划调度研究进行综述。在制造系统的范畴内考虑双重资源问题才刚开始，到2003年为止，在该领域内只有有限数量的研究报道^[2]。文献[3]建立了混合整数规划模型描述工人分配到生产单元的问题。文献[4]针对虚拟制造单元的零件、机器、工人分配问题，建立了具有双重资源约束的多目标规划数学模型。Li等^[5]对于DRC制造系统的生产调度问题进行了深入的研究。前人对于DRC制造系统的研究取得了不错的成果，但都是针对确定性的DRC制造系统，而对于具有随机性因素的DRC制造系统却未见有报道。因此，如何对随机性的DRC制造系统进行有效的建模是非常重要而且十分困难的。

对于随机性的制造系统的研究，一种有效的方法是建立排队网模型。文献[6]对排队网建模在制造系统中的早期应用进行了全面的综述，根据系统的在制品数量是否确定，排队网可分为开排队网、闭排队网和混合排队网。根据缓冲区的容量是否有限，排队网又可分为有限缓冲和无限缓冲的排队网。而由于系统的堵塞现象，有限缓冲的开排队网模型往往不具有乘积形式的解^[7]。对于有限缓冲的开排队网建模，国外学者Van Vuuren等^[8]针对具有堵塞影响的有限缓冲多服务台串联排队网模型，提出了基于两站子系统分解的近似求解方法。Manitz^[9]考虑有限缓冲和一般分布服务时间的多阶段装配线排队网模型，同样是基于两站子系统分解提出了一种近似求解方法。Cruz等^[10]提出了一种分析M/G/c/c状态相关排队网模型的近似方法。Smith^[11]利用M/G/c/c的排队网模型对传送带式储运系统进行建模和分析。文献[12]采用两矩法计算有限缓冲的M/G/c/K的排队网模型的稳态性能指标并进行分析。文献[13]针对含自动物料运输环节的生产系统，建立了带分段环形路径AGV的制造系统排队网模型。现阶段学者对于有限缓冲的开排队网建模在制造系统中的应用取得了很大进展，但是还没有见到考虑制造系统中机器人与加工中心2种资源同时性约束的排队网建模。

而建立仿真模型是研究随机性的制造系统的另一种有效方法，陈冠中等^[14]针对智能制造系统AGV小车及缓冲区容量配置优化问题建立仿真模型。谢洁明等^[15]建立机器人单元的仿真模型，分析机器人单元性能参数的随机性对系统性能的影响。

虽然对于随机性的DRC制造系统建立仿真模型是一种相对简单且广泛适用的方法，但是仿真往往需要消耗大量的时间和资源，因此需要对随机性的DRC制造系统的排队网建模进行探索。为此，本文针对具有自动装卸搬运机器人的制造单元建立有限缓冲的开排队网模型。而针对有限缓冲的开排队网模型，常用的求解方法主要有精确解法，状态空间分解法^[16-17]和广义扩展法^[18-19]。本文拟采用精确解法求解该问题。

1 问题描述与模型假设

该智能制造单元由待加工件缓冲区 ${B_0}$ 、已加工件缓冲区 ${B_3}$ 、加工中心 ${W_1}$ 、加工中心 ${W_2}$ 、机器人 $R$ 和单向环形的机器人轨道组成。待加工件缓冲区 ${B_0}$ 是一个有限容量的缓冲区，待加工件从外部经 ${B_0}$ 进入制造单元，然后工件在 ${B_0}$ 上等待机器人移动进来装载并搬运至加工中心 ${W_1}$ 进行加工。待工件在 ${W_1}$ 完成加工后，再由机器人把 ${W_1}$ 上的工件搬运到 ${W_2}$ 进行加工。待工件在 ${W_2}$ 完成加工后，由机器人把 ${W_2}$ 上的工件搬运到 ${B_3}$ 。工件进入 ${B_3}$ 后，将会马上离开制造单元。

建立该智能制造单元的排队网模型如图1所示。该系统满足以下假设。

1) 到来的工件服从泊松分布，工件到达速率为 $\lambda $ ，即工件到达的间隔时间服从均值为 $t_0=\dfrac{1}{\lambda }$ 的指数分布，其概率密度函数为 $f(x) = \lambda {{\rm e}^{ - \lambda x}}, x {\text{＞}} 0$ 。

2) 工件之间都是相互独立的。

3) 当一个工件到达系统后，如果 ${B_0}$ 上的工件数量还没到达最大容量 ${N_0}$ 时，则该工件进入 ${B_0}$ 的队列等待，如果 ${B_0}$ 上的工件数量到达了最大容量 ${N_0}$ 时，则该工件将被拒绝进入系统。

4) 2个加工中心均只有一台机器，只能同时加工1个工件，加工时间服从均值为 $t_i=\dfrac{1}{{{\mu _i}}}$ $(i = 1, 2)$ 的指数分布，概率密度函数为 $f(x) = {\mu _i}{{\rm e}^{ - {\mu _i}x}}, x {\text{＞}} 0$ ，其加工时间都是相互独立的。

5) 对于 ${B_0}$ 上的工件，服务原则是先到先服务，不存在插队，即先进入系统的工件会优先被机器人搬运至下一节点。

6) 系统服从后阻塞机制。

7) 机器人的容量为1，即每次只能搬运1个工件。机器人在轨道上相邻两节点移动及装卸的总时间均服从均值为 $\dfrac{1}{v}$ 的指数分布，其概率密度函数为 $f(x) = v{{\rm e}^{ - vx}}, x {\text{＞}} 0$ 。而由于机器人的装卸时间相较于移动时间是非常小的，所以对于机器人在轨道上不相邻的两节点移动及装卸的总时间，可以近似地认为其服从均值为 $\dfrac{2}{v}$ 的指数分布。机器人每次运行之间是相互独立的。

8) 机器人每次只能执行一个任务，且每次执行完任务后，若有新的任务可以执行，则继续执行新的任务，若无新的任务，则机器人停留在原地变成空闲状态。

2 排队网模型的分析及求解 2.1 系统的状态空间分析

利用精确解法进行排队网建模，可以得到该排队网模型的系统状态空间为

$\begin{split}&\quad\quad S = \Big\{ {({n_0}, {w_1}, {w_2}, {w_4}, {b_4})|0 {\text{≤}} {n_0} {\text{≤}} {N_0}, {w_1} = 0, 1, 2, {w_2} = }\\ &{0, 1, 2, {w_4} = 0, 1, 2, 3, {b_4} = - 1, 0, 1, 2, 3} \Big\}{\text{。}}\end{split}$

其中， ${n_0}$ 为 ${B_0}$ 中工件的数量； ${w_1}$ 为 ${W_1}$ 的状态， ${w_1} = 0$ 为 ${W_1}$ 中没有工件， ${w_1} = 1$ 为 ${W_1}$ 中有工件正在加工， ${w_1} = 2$ 为 ${W_1}$ 中有工件已经完成加工； ${w_2}$ 为 ${W_2}$ 的状态，其含义与 ${w_1}$ 类似； ${w_4}$ ， ${b_4}$ 为机器人正在执行从 ${w_4}$ 前往 ${b_4}$ 的任务，其中， ${w_4} = 0$ 为 ${B_0}$ ， ${w_4} = 1$ 为 ${W_1}$ ， ${w_4} = 2$ 为 ${W_2}$ ， ${w_4} = 3$ 为 ${B_3}$ ， ${b_4} = - 1$ 为机器人没有可执行的任务，即机器人在 ${w_4}$ 处于空闲状态， ${b_4} = 0$ 为 ${B_0}$ ， ${b_4} = 1$ 为 ${W_1}$ ， ${b_4} = 2$ 为 ${W_2}$ ， ${b_4} = 3$ 为 ${B_3}$ ，如 $({n_0}, {w_1}, {w_2}, 3, 0)$ 为机器人正在执行从 ${B_3}$ 到 ${B_0}$ 的任务。

对于整个系统来说，因为工件一旦进去 ${B_3}$ 后，将马上离开系统，不会对前面的节点造成堵塞，因此可以不考虑 ${B_3}$ 的状态，只考虑 ${B_0}$ ， ${W_1}$ ， ${W_2}$ 和机器人的状态。那么，整个系统中可容纳的工件的最大数量为 ${N_0} + 2$ 。因此，把系统的状态分为4类，共包含 $17{N_0} + 16$ 种状态，其状态转移分析如下。1) 系统中有0个、1个和2个工件，包含26种状态。2) 系统中有 $i$ 个工件，其中， $3 {\text{≤}} i {\text{≤}} {N_0} - 1$ ，包含 $17({N_0} - 3)$ 种状态。3) 系统中有 ${N_0}$ 个工件，包含17种状态。4) 系统中有 ${N_0} + 1$ 个和 ${N_0} + 2$ 个工件，包含16种状态。其中，图2表示系统中有 $i(3 {\text{≤}} i {\text{≤}}{N_0} -$ $ 1)$ 个工件时的状态空间转移图。

2.2 系统的稳态分析及状态平衡方程

从系统的状态转移分析可以看出，系统状态空间中的所有状态都是常返的，系统符合马尔可夫过程。系统中每个状态流出的速率与其他状态流入该状态的速率相等，因此可以推导出CTMC状态转移平衡方程，并且系统中每个状态都拥有自己独特的状态转移平衡方程，一共有 $17{N_0} + 16$ 个状态转移平衡方程。如针对系统的状态 $(i - 2, 1, 2, 1, 2)$ ，其中 $3 {\text{≤}} i {\text{≤}} {N_0} - 1$ ，一共有 ${N_0} - 3$ 种状态，以 ${P_{({n_0}, {w_1}, {w_2}, {w_4}, {b_4})}}$ 记为系统处于状态 $({n_0}, {w_1}, {w_2}, {w_4}, {b_4})$ 的概率，其状态转移平衡方程为

$\begin{split}&\quad\quad{P_{(i - 2, 1, 2, 1, 2)}} \times (\lambda + v + {\mu _1}) = {P_{(i - 2, 1, 1, 1, - 1)}} \times {\mu _2} + \\ &{P_{(i - 2, 0, 2, 0, 1)}} \times v + {P_{(i - 3, 1, 2, 1, 2)}} \times \lambda {\text{。}}\end{split}$

对于CTMC，所有状态都类似地列出其状态转移平衡方程，其系数矩阵 ${{Q}} = [{q_{ij}}]$ 可以转换成一个离散时间马尔可夫链的概率转移矩阵 ${{Q}}' = {{Q}}/q + {{I}}$ ，其中， $q$ 是一个常数且 $q {\text{＞}} \max \left| {{q_{ij}}} \right|$ ， $I$ 是单位矩阵。对于离散时间马尔可夫链的状态概率向量 $P$ ，对初始的状态概率向量采用反复的矩阵乘法进行迭代，即 $P(n) = P(n - 1){{Q}}'$ ， $P(n)$ 表示系统 $n$ 步转移后的状态概率。当系统的状态概率由 $n - 1$ 步和 $n$ 步转移后足够接近的时候，系统则达到了期望的精度水平的稳态，此时 $P(n)$ 就是所要求的系统稳态概率。因此，系统稳态概率 $P(n) = P(n){{Q}}'$ ，即 $P(n) = P(n)({{Q}}/q + {{I}})$ ，可以化简为 $0 = P(n){{Q}}$ 。因此把CTMC所有状态的状态转移平衡方程联立，组合成状态转移平衡齐次线性方程组，就可以求解出系统各个状态的稳态概率^[20]。

对于马尔可夫链的分析求解，一般其线性方程组的求解方法有直接求解法和迭代求解法。对于本文要求解的状态转移平衡方程组，因其系数矩阵是大型稀疏矩阵，因此本文将采用收敛速度更快的Gauss-Seidel迭代法来求解状态转移平衡方程组。对于CTMC状态转移平衡方程在迭代的过程中的收敛性证明，其证明过程可参考文献[20]，这里不阐述详细的证明过程。

2.3 系统性能指标的求解

在该模型中，系统的平均产出率和设备平均利用率可直接利用系统的状态概率求解出来，而平均生产周期则需要应用利特定律(Little’s Law)求解。利特定律的表达式为 $L = \lambda \times T$ 。其中， $L$ 表示系统的平均在制品数量， $\lambda $ 表示系统的有效到达率， $T$ 表示平均生产周期，即工件进入系统到离开系统的平均时间。而当系统处于稳态时，系统的有效到达率等于产出率。因此，通过求解出系统的状态概率以及系统性能指标的定义，可以推导出系统性能指标的计算公式。

平均产出率 $\eta = v \times \displaystyle\sum\limits_{{n_0} = 0}^{{N_0}} {\displaystyle\sum\limits_{{w_1} = 0}^2 {\displaystyle\sum\limits_{{w_2} = 0}^2 {{P_{({n_0}, {w_1}, {w_2}, 2, 3)}}} } } $ 。

平均在制品数量 $L = \displaystyle\sum\limits_{n = 0}^{{N_0} + 2} {n {P_n}} $ ，其中， $n$ 为系统中工件的数量， ${P_n}$ 为系统中有 $n$ 个工件的概率。

平均生产周期 $T = L/\eta $ 。

设备的平均利用率

$\begin{split}&\quad\quad\gamma =\Bigg(\displaystyle\sum\limits_{{n_0} = 0}^{{N_0}} {\displaystyle\sum\limits_{{w_2} = 0}^2 {\displaystyle\sum\limits_{{w_4} = 0}^3 {\displaystyle\sum\limits_{{b_4} = - 1}^3 {{P_{({n_0}, 1, {w_2}, {w_4}, {b_4})}}} } } } + \\ &\displaystyle\sum\limits_{{n_0} = 0}^{{N_0}} {\displaystyle\sum\limits_{{w_1} = 0}^2 {\displaystyle\sum\limits_{{w_4} = 0}^3 {\displaystyle\sum\limits_{{b_4} = - 1}^3 {{P_{({n_0}, {w_1}, 1, {w_4}, {b_4})}}} } } }\Bigg)/2{\text{。}}\end{split}$

机器人的平均利用率

$\begin{split}&\quad\quad\beta = 1 - \sum\limits_{{n_0} = 0}^{{N_0}} {\sum\limits_{{w_1} = 0}^2 {\sum\limits_{{w_2} = 0}^2 {\sum\limits_{{w_4} = 0}^3 {{P_{({n_0}, {w_1}, {w_2}, {w_4}, - 1)}}} } } }{\text{。}}\end{split}$

3 算例及系统性能分析

为验证排队网模型及其求解方法的有效性，采用Tecnomatix Plant Simulation仿真软件建立对应的仿真模型，并仿真统计得到系统的性能指标。本文通过设计一系列算例，分别进行排队网模型数值求解与仿真模型实验统计，并对2种方法的求解结果进行对比分析，说明排队网模型及其求解的有效性和精确度。所有的数值计算和仿真实验均在同一台计算机上进行，其操作系统环境为PC windows 7，硬件环境为Intel (R) CPU 1.6 GHz，4.00 GB RAM。排队网模型采用matlab软件编程求解，仿真模型在Tecnomatix Plant Simulation平台上进行实验统计。

3.1 仿真模型

建立的该制造单元排队网模型对应的仿真模型如图3所示，排队网模型中定义的系统状态空间在仿真模型有着相应的状态。比如，当 ${B_0}$ 中有2个工件， ${W_1}$ 和 ${W_2}$ 都在进行加工，机器人停在 ${W_1}$ 处并且空闲，此时仿真模型的状态则对应于排队网模型中的状态 $(2, 1, 1, 1, - 1)$ ，如果仿真模型的下一个事件为有一个工件到达了，则此时仿真模型的状态对应于排队网模型中的状态 $(3, 1, 1, 1, - 1)$ ，如果在这之后仿真模型的下一个事件为 ${W_2}$ 完成加工在等待，机器人需要从 ${W_1}$ 移动到 ${W_2}$ ，则此时仿真模型的状态对应于排队网模型中的状态 $(3, 1, 2, 1, 2)$ 。在仿真模型中，对于工件到达的间隔时间、加工时间、机器人的移动及装卸时间服从指数分布，是采用随机数发生器来进行仿真模拟的。

在仿真模型中，系统的平均产出率是工件的总产出量除以仿真运行时间，系统的平均生产周期是一个工件从进入系统到离开系统的平均时间，设备的平均利用率是设备的总加工时间除以仿真运行时间，机器人的平均利用率是机器人的总移动及装卸时间除以仿真运行时间。

在仿真实验中，为了防止随机性因素造成的影响，每组算例分别在95 %置信水平下进行5次重复的仿真实验(每次仿真实验分别釆用不同的随机数种子流)。为了防止系统初始状态对系统稳态数据造成的影响，每次实验仿真1 000 d后再进行统计。

3.2 排队网模型的有效性验证

本文设计一系列不同的算例进行计算，具体参数如表1所示。其中第0组为参照组。第1组在第0组的基础上保持其他参数不变，改变工件的到达时间。第2组在第0组的基础上保持其他参数不变，改变加工中心的加工时间。第3组在第0组的基础上保持其他参数不变，改变机器人的移动和装卸速率。第4组在第0组的基础上保持其他参数不变，改变缓冲区 ${B_0}$ 的容量。

对上面一系列算例采用排队网模型进行数值计算和仿真统计求解，其结果对比如表2所示。

从表2的数据对比中可以看出，本文针对具有自动装卸搬运机器人的制造单元建立的排队网模型及系统性能的计算程序是有效的，并且其2种方法计算结果相对误差非常小，所有的性能指标的相对误差都在0.25%以下。

3.3 基于排队网模型对系统性能的分析 3.3.1 到达时间对系统性能的影响

结合第1组和第0组的结果可以看出，随着工件到达时间增大，系统的平均产出率、平均生产周期、设备平均利用率、机器人平均利用率都会逐渐减小，如图4所示。

当工件的到达时间减小到一定程度时，系统的平均产出率增大的趋势越来越小。这是因为此时系统的资源几乎已经到了最大的利用程度，通过设备平均利用率、机器人平均利用率变化趋势也能说明这点。而缓冲区容量有限，当缓冲区容量达到最大时，系统会拒绝新来的工件，因此平均生产周期增大的趋势也越来越小。

当工件的到达时间增大到一定程度时，由于系统空闲的时间远大于系统加工的时间，因此系统各性能指标下降的趋势也越来越小。

3.3.2 加工时间对系统性能的影响

结合第2组和第0组的结果可以看出，随着设备的加工时间增大，系统的平均产出率减小，平均生产周期、设备平均利用率增大，机器人平均利用率先增大后减小，如图5所示。

当加工时间减小到一定程度时，系统的平均产出率增大的趋势越来越小，这是因为此时系统的其它资源几乎已经到了最大的利用程度。

当加工时间增大时，由于平均生产周期主要是由加工时间决定的，所以平均生产周期几乎呈线性增长。当加工时间很小时，加工时间增大，机器人利用率出现了小幅的增大然后减小。这是因为加工时间很小时系统空闲的时间很多，由于系统空闲的时间逐渐变少而机器人利用率逐渐增大，之后又随着加工时间增大，机器人因为要等待设备加工而利用率降低。随着加工时间增大，设备利用率增大的趋势越来越小，并且只趋近61%，这是由于加工时间服从指数分布造成的。

3.3.3 机器人速率对系统性能的影响

如图6所示，随着机器人速率增大，系统的平均产出率、设备平均利用率、机器人平均利用率都会逐渐增大，平均生产周期逐渐减小。

当机器人速率增大到一定程度时，系统的平均产出率、设备平均利用率、机器人平均利用率增大的趋势都越来越小，平均生产周期减小的趋势越来越小，这是因为受到了系统其它条件限制的原因。

3.3.4 缓冲区容量对系统性能的影响

如图7所示，随着缓冲区的容量增大，系统的平均产出率、平均生产周期、设备平均利用率、机器人平均利用率都会逐渐增大。

在缓冲区的容量大于3时，随着缓冲区的容量增大，系统的平均产出率增大的趋势越来越小，这是因为此时系统的资源几乎已经到了最大的利用程度，通过设备平均利用率、机器人平均利用率变化趋势也能说明这点。

随着缓冲区的容量增大，平均生产周期几乎呈线性增长，这是因为更多的工件进入系统等待了，因此平均每个产品在系统中等待的时间更加长了。

4 结论与展望

本文针对具有加工中心和机器人双重资源约束的制造单元，采用有限缓冲区的开排队网理论建立制造系统的模型，快速计算系统的性能指标，并且建立了对应的仿真模型。通过一系列的算例实验，验证排队网模型及其求解方法的有效性和精确度，最后基于排队网模型对机器人制造单元的系统性能进行分析。通过一系列的算例实验结果表明，在不同参数设定的情况下，具有双重资源约束的制造系统的排队网模型都是有效的，并且对于系统的性能计算精确度都非常高，相对误差都在0.25%以下。

本文建立的排队网模型虽然考虑系统的随机性和交互性，但没有考虑因紧急插单、工件返修、机器故障、员工休假及其他干扰因素等情形。因此，对于更加复杂的制造单元，以及对于更大规模的制造系统，需要提出新的建模思路和方法，并且计划将研究拓展至具有一般分布的输入输出的有限缓冲区开排队网模型，提出更加通用有效的求解方法。