基于排队网络模型的混流制造系统服务率优化*

张岩岩，白金花，武 福，李忠学

(兰州交通大学 机电工程学院，甘肃 兰州 730070)

1 引 言

随着市场竞争激烈化以及顾客需求个性化的增加，混流生产方式成为制造业普遍采用的一种生产组织方式，它可在基本不改变设施布局的前提下，在同一制造系统内加工出多种不同型号和数量的产品，具有很高的灵活性，在汽车、家电、玩具等产品制造企业中有着广阔的应用前景。拥有较好的服务率水平，是制造企业巩固和发展自身优势的一种持久能力，也是其最为关心和追求的效用指标，如何对制造系统的服务率水平进行优化也就成为人们普遍关心的问题。但在实际生产过程中，制造系统的服务率水平受到许多外部及内部复杂因素的制约，影响着生产计划、成本、供求关系等。在混流制造系统的实际运行中，不同类型的产品同时投入生产：对于同一类型的产品来说，不同工作单元的平均加工时间不同；而对于不同类型的产品来说，同一工作单元的平均加工时间又存在差异。因此对混流制造系统，服务率不合理造成资源浪费以及生产效率低下的问题较突出。

排队是制造系统中的一个常见现象。排队理论应用广泛，是研究一切服务系统的有效工具，这使得利用排队理论对混流制造系统优化问题进行研究成为可能。在已有的基于排队理论优化服务率的研究中， Grassmann，Chen和Kashyap[1]对具有状态无关到达率的M/G/1排队模型进行了研究，提出了确定其在稳定状态下最优服务率的方法；Jo[2]利用拉格朗日方法得到了开放式Jackson排队网络中的每一个节点在相邻两个稳定状态下最优服务率的相互关系，之后给出了计算每个节点最优服务率的有效算法；Song，Xing和Sun[3]研究了随机需求不可靠制造系统的最优服务率控制问题等。而基于排队理论对更为复杂的混流制造系统进行服务率优化的相关研究则显得不足。鉴于混流制造系统的效益性、现实性以及排队理论的实用性、有效性，本文应用排队理论对混流制造系统的服务率优化问题进行研究。

2 排队理论及排队网络

排队理论又称为随机服务系统理论，它通过对服务对象的到达模式及服务时间的统计分析，得到等待时间、排队长度、忙期长短等相关参数的统计规律，根据这些参数的统计规律来改进服务系统的结构或者重组服务系统，使系统的特定指标达到最优。

一个排队系统中存在多项服务，每项服务都有自己的服务机构。在这种情况下，系统中也就存在多个排队队列，整个排队系统便形成排队网络。排队网络是一个比较复杂的服务系统，但其应用广泛。排队网络研究的首个突破是20世纪五六十年代的Jackson网络模型[4]，使得计算机系统建模工作得到了很大的简化。随后，针对更加复杂的系统建模，如引入不同类型的顾客和新的服务规则，Baskett，Chandy，Muntz和Palacios[5]等学者提出了BCMP定理。而对于更加一般的排队网络，由于对其性能参数的评价难以获得准确解，许多学者研究并提出了一些如均值分析、分解等近似方法。本文则基于开放式排队网络理论对混流制造系统建立物理模型和数学模型。

3 系统建模

在混流制造系统中，各种类型的产品混合连续地到达各个工作单元接受服务，不同类型的产品在各个工作单元接受的服务时间不同。依据开放式排队网络模型的特征(顾客由外部进入网络，接受完所有的预定服务后离开网络)，一个混流制造系统可以等效为一个开放式排队网络模型。

对于混流排队网络模型，在计算过程中不同类型的产品之间相互干扰，使得单独针对每一类产品求解排队系统的相关指标之后再线性相加的策略成为不可能，因此在建立数学模型之前，聚合到达每一个工作单元的所有各类产品为一类混合产品，并将其看作单一类型的产品，如图1所示，则整个混流排队网络就简化成为了单一产品类型的排队网络，其中的每一个工作单元都等效为G/G/1排队模型。考虑工作单元i，进而可以根据其典型队列模型图和式(1)～(3)求得每一类产品和混合产品到达每一个工作单元的平均到达率λir和λi以及混合产品在制造系统中的转移概率pij。图2所示为i工作单元的典型队列模型图。

图1 不同类型产品的聚合

图2 i工作单元的典型队列模型图

(1)

(2)

(3)

其中:式(1)可由r类产品在排队网络中的转移概率矩阵确定。

产品在系统中逗留会产生一定的损失费用，提高工作单元的服务率虽然可以减少产品的逗留时间，但也意味着运行成本的增加，不利于企业实现高效益的目标，如图3所示。

图3 排队系统的期望费用图

综合考虑上述因素，如果取目标函数F(μi1，μi2，…，μiR)为工作单元i单位时间内的服务费用与产品逗留损失费用之和的期望值，则基于排队理论建立混流制造系统的期望费用模型如下：

minF(μi1，μi2，…，μiR)=Fsiμi+FwiE(Li)

(4)

约束条件：

μir>λir

(5)

其中:Fsi为工作单元i单位时间内提高一个单位的服务率时所增加的服务费用，Fwi为混合产品在工作单元i内逗留单位时间的损失费用，μi为工作单元i对混合产品的平均服务率，E(Li)为工作单元i单位时间内的平均产品数。

4 模型分析

4.1 E(μi)的计算

当混流制造系统处于稳态时，工作单元i加工r类产品的概率和平均加工时间分别为λir/λi和E(tir)，则工作单元i对混合产品的平均加工时间E(ti)为：

(6)

因此，得到工作单元i对混合产品的平均服务率为：

(7)

4.2 E(Li)的计算

Little定理 设产品按照泊松流到达系统，参数为λ，服务机构的服务时间服从参数为μ的负指数分布。排队系统中单位时间内的平均产品数为E(L)，平均逗留时间为E(T)，那么两者之间的关系为：

E(L)=λE(T)

(8)

对于一般服务时间分布或一般到达的排队系统，并非在任何时刻都有马尔可夫性，对其队列长度的平稳分布进行求解相当困难，也没有明显可利用的结果。由于G/G/1排队模型中仅知道均值和变异系数的平方两个参数，本文利用Kingman[6-7]近似式(9)计算混合产品在接受服务前的平均等待时间E(Tqi)：

(9)

(10)

由于工作单元对产品的服务规则是先到先服务，任何产品都不具有优先权，所有的产品在同一个排队队列中等待接受服务，因此混合产品在工作单元i内的平均逗留时间为：

(11)

根据Little定理，得到工作单元i单位时间内的平均产品数为：

(12)

1-qki}]

(13)

ωi=[1+4(1-ρi)2(υi-1)]-1

(14)

(15)

(16)

根据变异系数的定义得到工作单元i对r类产品加工时间的变异系数平方为：

(17)

(18)

因此，工作单元i对混合产品加工时间的二次矩为：

(19)

则工作单元i对混合产品加工时间的变异系数平方为：

(20)

因此，最终的目标函数可以记为：

(21)

约束条件：

μir>λir

(22)

4.5 E(Li)的简化处理

(23)

因此，目标函数式(21)可以修正记为：

(24)

约束条件：

μir>λir

(25)

选择合适的方法求解目标函数式(24)，就使得计算过程在很大程度上得到简化。

5 算 例

一混流制造系统同时生产2种类型的产品，r(r=1，2)类产品按照参数为λar=0.0014(件/s)的泊松过程从外部到达制造系统，加工路径如图4所示[11]，其中圆形节点表示制造系统中的工作单元，节点中的数字标注表示工作单元编号，工作单元之间转移路径上的数字标注表示产品在一个工作单元服务完成后进入下一个工作单元接受服务的概率。

图4 两类产品的加工路径

现假设两类产品在各个工作单元接受的服务时间服从负指数分布且所有的到达过程和服务过程相互独立。为了在减小产品逗留可能性的同时又不至于因其服务率过高而增加运行成本，对此混流制造系统各个工作单元的服务率进行控制优化。表1所示为每一个工作单元单位时间内每一类产品的逗留费用与服务费用。

表1 每个工作单元单位时间内每类产品 的逗留费用与服务费用

5.1 确定混合产品访问每个工作单元平均到达率

(1) 两类产品在排队网络中的转移概率矩阵P1和P2

本算例中，转移概率矩阵Pr(r=1，2)中的元素pijr(i=1，2，…，4；j=1，2，…，4)表示r类产品在工作单元i服务完成后进入工作单元j接受服务的概率，p0ir表示r类产品从网络外部到达工作单元i接受服务的概率，pi5r表示r类产品在工作单元i服务完成后离开网络的概率。

(2) 计算两类产品访问每一个工作单元的平均到达率

根据式(1)，结合转移概率矩阵Pr可以得到两类产品访问每一个工作单元的平均到达率为：

(λ11，λ21，λ31，λ41)=(0.0016，0.0012，0.0015，0.0000)

(λ12，λ22，λ32，λ42)=(0.0017，0.0019，0.0016，0.0035)

因此混合产品访问每一个工作单元的平均到达率为：(λ1，λ2，λ3，λ4)=

(0.0033，0.0031，0.0031，0.0035)

5.2 确定混合产品在制造系统中的转移概率矩阵P

由式(3)计算可得混合产品在制造系统中的转移概率矩阵P为：

P=

5.3 确定每一个工作单元对混合产品加工时间的变异系数平方

在满足约束μir>λir的前提下，为方便计算，选择μir的初始值为：

(2) 计算每一个工作单元对混合产品加工时间的变异系数平方

在本算例中，每一类产品在系统中各个工作单元接受的服务时间服从负指数分布，因此：

=(1.0000，1.0000，1.0000，1.0000)

5.4 确定混合产品到达每一个工作单元的平均间隔时间的变异系数平方

(ρ1，ρ2，ρ3，ρ4)=

(0.0033，0.0031，0.0031，0.0035)

由于每一类产品都是按照泊松流到达制造系统的，而相互独立的泊松流的合成流仍然为泊松流，因此：

(1.0000，1.0000，1.0000，1.0000)

由式(13)～(16)可计算出混合产品到达每一个工作单元的平均间隔时间的变异系数平方为：

(1.0056，0.9946，0.9999，0.9996)

5.5 建立目标函数

将上面所得参数代入式(24)就可得到各个工作单元最小化期望费用的目标函数。以工作单元1为例，其目标函数为：

minF(μ11，μ12)

(26)

约束条件：

μ11>0.0016

(27)

μ12>0.0017

(28)

5.6 目标函数求解算法设计

目标函数式(26)是一个不等式约束非线性规划问题，可利用内点惩罚函数法将原约束优化问题转化为无约束优化问题，然后选取坐标轮换法对目标函数进行求解，在坐标轮换法的每一轮迭代过程中利用黄金分割法确定各个步长，而黄金分割法实施的初始搜索区间[a，b]则通过进退法确定。

(1) 内点惩罚函数法的实现步骤

(29)

② 在可行域之内选取初始点x(0)，令k=1；

③ 选取适当大小的初始惩罚因子r(0)；

④ 从x(k-1)点出发，利用坐标轮换法求解minφ(x，r)的极小点xk1(r(k))；

⑤ 判断精度|xk1(r(k))-xk1(r(k-1)) |≤ε1，若满足，停止迭代计算，并以xk1(r(k))作为原问题的近似极小值，否则转向步骤⑥；

⑥ 取r(k+1)=cr(k)(c≈0.1～0.7)，x(k-1)=xk1(r(k))，k=k+1，转向步骤④。

(2) 坐标轮换法的实现步骤

(3) 黄金分割法的实现步骤

① 在进退法确定的初始搜索区间[a，b]内取两个试算点α1和α2，计算函数值φ(α1)和φ(α2)：

α1=a+0.382(b-a)

(30)

α2=a+0.618(b-a)

(31)

② 比较φ(α1)和φ(α2)

第一种情况，φ(α1)<φ(α2)：

置换符号：

α2→b，α1→α2，φ(α1)→φ(α2)

在新区间内取一个新点α1，计算函数值φ(α1)：

α1=a+0.382(b-a)

(32)

以缩小后的新区间作为原区间，重复步骤②。

第二种情况，φ(α1)>φ(α2)：

置换符号：

α1→a，α2→α1，φ(α2)→φ(α1)

在新区间内取一个新点α2，计算函数值φ(α2)：

α2=a+0.618(b-a)

(33)

以缩小后的新区间作为原区间，重复步骤②。

(4) 进退法的实现步骤

① 给定初始点a(0)和初始步长h0，令h0→h，a(0)→a(1)，0→k2；

② 令a(4)=a(1)+h，k2+1→k2；

③ 比较φ(a(4))和φ(a(1))，若φ(a(4))<φ(a(1))，转向步骤④，否则转向步骤⑤；

④ 令a(1)→a(2)，a(4)→a(1)，φ(a(1))→φ(a(2))，φ(a(4))→φ(a(1))，h=2h，转向步骤②；

⑤ 若k2=1，转向步骤⑥，否则转向步骤⑦；

⑥ 令h=-h，a(4)→a(2)，φ(a(4))→φ(a(2))，转向步骤②；

⑦ 令a(2)→a(3)，a(1)→a(2)，a(4)→a(1)，终止计算，极小点所在的两个可能初始搜索区间为[a，b]=[a(1)，a(3)]或[a，b]=[a(3)，a(1)]。

5.7 求解目标函数

根据内点惩罚函数法，引入惩罚因子，将目标函数式(26)转化为无约束优化问题。构造的惩罚函数为：

(34)

算法程序由MATLAB语言编写，求解过程中的所有程序在MATLAB7.5.0环境中运行。对工作单元1服务率优化的最终结果为：

(μ11，μ12)=(0.0 288，0.0 358)

同样地，可以求出其他工作单元加工每一类产品的最优服务率分别为：

(μ21，μ22)=(0.0 303，0.0 314)

(μ31，μ32)=(0.0 340，0.0 364)

(μ41,μ42)=(0，0.0 314)

即完成对算例中混流制造系统的服务率优化工作。

6 结 语

应用排队理论对混流制造系统的服务率优化问题建立了物理模型和数学模型，并利用最优化方法和MATLAB程序对所建数学模型进行了优化求解。在此过程中，将混流制造系统等效为了开放式排队网络模型，每一类产品到达网络的时间间隔以及在各个工作单元接受的服务时间都是服从负指数分布的。事实上，在每一类产品的到达时间间隔以及服务时间都服从一般分布的情形下，文中所建模型和相关计算过程同样适用。算例计算结果表明了应用排队理论对混流制造系统进行服务率优化的方法是可行有效的，从而为制造系统服务率优化问题的解决提供了一种新的思路。

参考文献：

[1] Grassmann W K，Chen X, Kashyap B R K. Optimal Service Rates for the State-Dependent M/G/1 Queues in Steady State[J]. Operations Research Letters，2001，29(2)：57-63.

[2] Jo K Y. A Lagrangian Algorithm for Computing the Optimal Service Rates in Jackson Queuing networks[J]. Computers and Operations Research，1989，16(5)：431-440.

[3] Song D P，Xing W, Sun Y X. Optimal Service Rate Allocation Policy of an Unreliable Manufacturing System with Random Demands[J].控制理论与应用，1998，15(4)：621-626.

[4] Jackson J R. Networks of Waiting Lines[J]. Operations Research，1957，5(4)：518-521.

[5] Baskett F，Chandy K M，Muntz R R, et al. Open，Closed，and Mixed Networks of Queues with Different Classes of Customers[J].Journal of the Association for Computing Machinery，1975，22(2)：248-260.

[6] Kingman J F C. On Queues in Heavy Traffic[J]. Journal of the Royal Statistical Society. Series B：Methodological，1962，24(2)：383-392.

[7] Kingman J F C. The Heavy Traffic Approximation in the Theory of Queues[C]. Proceedings of the Symposium on Congestion Theory. Chapel Hill：The University of North Carolina Press，1964：137-169.

[8] Whitt W. The Queueing Network Analyzer[J]. The Bell System Technical Journal，1983，62(9)：2779-2815.

[9] Bitran G R, Tirupati D. Tradeoff curves，Targeting and Balancing in Manufacturing Queueing Networks[J].Operations Research，1989，37(4)：547-564.

[10] Bitran G R, Morabito R. State-of-the-art survey：Open Queueing Networks：Optimization and Performance Evaluation Models for Discrete Manufacturing Systems[J]. Production and Operations Management，1996，5(2)：163-193.

[11] Curry G L, Feldman R M. Manufacturing Systems Modeling and Analysis[M]. USA: Springer，2010.