有限小批量定制生产伯努利三机器装配系统实时性能分析

刘畅 贾之阳

装配系统是生产实践中最基本的生产系统结构之一.系统中最终的产品通常由两个或两个以上组件装配构成(例如:汽车、家电、消费电子).相比在稳态分析研究方面取得的大量成果,装配系统的暂态过程仍然未被深入研究.这主要是由于不同零件生产线的相互作用导致了对装配系统的分析要比传统的串行线复杂得多.装配系统的早期研究通常只考虑多队列单服务器的情况,即几种类型的零件到达一个装配机器,从而被执行组装操作[1].后来,Lipper等[2]和Kuo等[3]研究了有限队列容量的三机器系统.在这些研究中,两个服务器代表组件的生产,而另一个服务器代表装配操作.此外,Manitz[4]对基于排队模型的装配系统进行了研究.具有不可靠机器和有限缓冲区容量的装配系统的稳态性能评估在一些论文中进行了讨论[5−10].

需要注意的是,虽然在过去几十年里存在大量的串行线和装配生产线的研究工作,但是其中的绝大部分都是假设系统运行在稳定状态之下的.与此同时,具有有限缓冲区容量和不可靠机器的生产系统的暂态性能只在最近的一些论文中得到初步的研究,其中大部分集中在基于伯努利机器的串行线系统研究[11−14].此外,Meerkov等[15]研究了使用几何可靠性机器模型的串行线的暂态性能.当一个生产系统根据客户订单或需求预测,每次只生产一个批次(或批量)的同类产品时,我们称之为有限小批量生产运行,基于此,Jia等[16−18]近两年研究了串行生产线和闭环生产系统的实时系统性能.

基于暂态的装配系统的分析研究仅在以下几篇论文中出现,Alexander等[19]研究了一类具有无限队列容量的单机器马尔科夫类装配排队系统的暂态吞吐量.Jia等[20−21]研究了复杂装配系统的暂态性能,但是假设了系统具有无限的原材料供应量.与此同时,值得关注的是,近年来智能制造技术的发展对生产系统的暂态和动态特性研究提出了更高的要求,这对于研究相应的实时生产控制算法也至关重要.因此,本文的目的是研究有限小批量定制化生产运行下,具有有限缓冲区容量的三机装配系统的基于暂态的性能评价.

本文余下的部分组织如下:第1节介绍了本文所研究的系统的前提假设并定义了所关注的系统性能指标.第2节详细讨论了系统数学模型的建立以及相关性能指标的推导过程.然后,提出的一种基于分解的用于近似原始系统性能的方法在第3节中被详细介绍.所提出方法的准确性通过数值实验进行了验证并分析展示.最后,第4节给出了总结和未来工作展望.

1 数学模型和性能指标

1.1 数学模型

考虑如图1所示的一个三机装配系统,其中圆形表示机器,矩形表示缓冲区.

图1 三机装配系统Fig.1 Assembly production system with three Bernoulli machines

系统根据以下假设来进行定义:

1)系统的最终产品(F0)需要两个组件.一个组件(R1)由机器m1处理.我们称系统的这一部分(从机器m1到b1)为零件生产线1.类似地,另一个组件(R2)由机器m2处理.系统的这一部分(从机器m2到b2)称为零件生产线2.

2)机器m0从零件生产线1和零件生产线2各取一个完成的零件装配组成一个成品.

3)机器mi,i=0,1,2,拥有恒定且相同的周期时间τ.以一个加工周期τ为一段,将时间轴分段.所有机器在一个新的生产批次开始时运行.小批量定制生产下的每个批次具有有限的产量,每个生产批次的规模为B.每台机器在加工完规定数量的工件后立即停止工作.

4)机器遵循伯努利可靠性模型,即,机器mi,i=0,1,2,如果既没有被阻塞也没有饥饿,在一个时间间隙(即加工周期)里加工处理一个工件的概率是pi,未能加工处理一个工件的概率是1−pi.参数pi∈(0,1)称为机器mi的效率.

5)每一个在制品缓冲区,bi,i=1,2,可以用其容量Ni来表征,0

6)如果机器m0在时间间隙n内处于工作状态,缓冲区b1或者b2在时间间隙开始时为空,则机器m0在时隙n内会饥饿.机器m1和m2在一个批次生产结束前不会出现饥饿的情况.

7)如果机器mi,i=1,2,在时间间隙n内处于工作状态,缓冲区bi在时间间隙开始时有Ni个在制品工件,并且装配机器m0没能从其中取走一个工件进行处理(由于故障或源自另一条零件生产线的饥饿情况),则机器mi,i=1,2,在时隙n内被阻塞.即,加工前阻塞机制.同时假设机器m0任何时候都不会被阻塞.

注1.值得注意的是,在许多生产系统中,机器周期时间几乎是恒定或接近恒定的.这样的情况大多见于汽车、电子、电器等行业的生产系统.还需注意到,伯努利可靠性机器模型是适用于平均故障时间接近机器的加工周期的情况(参见使用伯努利模型为实例的文献[22−24]).具有其他可靠性机器模型(例如:几何型、指数型、威布尔型、对数正态型等)的生产系统将在今后的工作中进一步研究.

注2.基于批次的生产广泛用于各种制造系统(小规模,中等规模,甚至大规模生产,单型或多类型产品生产等).一个批次有时被称为一个分组或一个订单.

注3.由于通常定制化生产下每个批次生产数量是有限的,整个生产过程部分或完全是在暂态下进行的.因此,严格来说,稳定状态分析不再适用,而基于暂态的系统分析是必要的.

注4.上述的模型仅仅包括两条零件生产线和一个装配操作机器.每条零件生产线仅包含了一台机器和一个缓冲区.每条零件生产线拥有多台机器和缓冲区,以及拥有多条零件生产线的复杂装配系统具有类似的假设,并且这样的装配系统会在未来工作中进一步研究.

1.2 性能指标

在上述定义的模型框架下,我们感兴趣的性能指标包括:

1)生产率PR(n):在时间间隙n+1里,机器m0生产工件个数的期望;

2)消耗率CRi(n):在时间间隙n+1里,机器mi,i=1,2,消耗原材料个数的期望;

3)在制品库存水平,WIPi(n):在时间间隙n里缓冲区bi,i=1,2,中的在制品个数的期望;

4)阻塞率BLi(n):机器mi,i=1,2,在时间间隙n+1里被阻塞的概率.

由于机器m0可能由于任意一条零件生产线而饥饿,我们定义机器饥饿率为:

一种通过递归聚合来估计这些稳态性能值的方法在文献[22]中被提出.在本文中,我们提出了在有限量定制生产运行下评估这些暂态性能指标的方法.

此外,使ct表示机器m0完成生产B个产品的时间.将其均值表示为:

2 系统性能精确分析

2.1 性能分析

用fi(n)表示机器mi在时间间隙n结束时已经生产的工件总数量,用hi(n)表示在时间间隙n结束时缓冲区内的在制品工件数量.显而易见,

那么,不失一般性,系统可以用一个状态为(h1(n),h2(n),f0(n))的马尔科夫链来表征,其中,

显然,此马尔科夫链的最大系统状态数为

需要注意,有一些系统状态是不可达到的,比如,(1,1,B),因为机器m1和m2在加工好B个工件后立刻停止了运作.换句话说,在任意一个时间间隙里,h1+f0≤B,并且h2+f0≤B.

为了计算这一马尔科夫链中的状态间转移概率,我们首先如表1排列系统的状态.

表1 系统状态排序Table1 Arrangement of the system states

其中,

同时也需要注意,在每个时间间隙中,系统状态的样本空间是由机器23种的工作状态所组成的.那么,

因此,在每一个时间间隔开始时,对系统的每一个可达状态i,i∈{1,···,Q},如果,并且,可以枚举所有的23种机器状态的组合,根据系统动态性式(4)来确定相应的在这一时间间隔结束时的结果状态j,j∈{1,···,Q}.然后,对于得到相同结果状态的机器状态组合情况,使用式(5)来计算相应的转移概率,并将这些概率相加,最终得到一个时间间隔里,从起始的系统状态i到结果状态j的转移概率.对于所有符合条件的系统状态重复这一步骤.

系统的实时性能可以通过下式计算:

其中

其中,01,k和J1,k分别代表1×k的零矩阵和元素全为1的矩阵.与此同时,i×j维矩阵Ci×j=[Ii···Ii]表示由j/i个单位矩阵Ii组成的矩阵.

3 基于分解的性能评估

上面描述的精确分析可以扩展到更大的系统,即每个零件生产线中有多台机器的系统.然而,随着机器数量M,缓冲区容量,和生产规模B的增长,马尔科夫链状态的数量呈指数型增长,这将导致对大型的复杂装配系统的分析变得不可能.因此,本节提出了一种基于分解的算法,并将其应用于三台伯努利机器的小型装配系统.相应的研究结果将在未来的工作中扩展到更通用的大型系统中.

3.1 基于分解的概念

文献[8]提出一种分解方法,将原系统分解为一对串行线:上线和下线,研究了基于无限原材料供应的装配系统的稳态性能.此外,我们以前的工作[20−21]解决了这类系统的暂态性能研究的问题.与此同时,当考虑到小批量有限量生产运行下的串行线,基于暂态的系统性能近似评估也在我们以前的工作[16−17]中进行了讨论.在这一节中,我们将基于有限量生产运行下系统的暂态性能分析扩展到三台机器的装配系统性能分析研究中.对由多台机器组成的零件生产线或多条零件生产线以及多个装配操作的复杂装配系统,将在今后的研究中进行分析.

图2 辅助装配系统Fig.2 Auxiliary assembly system

具体而言,引入三种辅助系统/生产线来分析此类系统.辅助装配系统(图2所示)首先被引入,这一辅助装配系统具有所有原始的机器和缓冲区,但假设具有无限的原材料供应.只有原来的装配操作中装配机器m0处于工作状态并且缓冲区b1非空的情况下,才可能处于工作状态.因此,让hi(n)表示在时间间隙n结束时缓冲区bi中的在制品零件数,可以通过式(9)估算:

图3 辅助双机串行线Fig.3 Auxiliary two-machine lines

用于分析这种辅助两机串行线的方法在文献 [11]中被提出.具体来说,用和,分别来表示辅助两机生产线中b1和b2在时隙n结束时,有i个工件的概率.令. 根据文献 [11],的演化可表示为:

为了研究这一辅助装配系统的暂态性能,分析方法可以参考我们以前的工作(参阅文献[20]).具体而言,使用辅助两机线(图3)来近似分析.

上生产线通过移除辅助装配系统中的机器m2和缓冲区b2来构造.考虑到这种修改,组装机器m0由效率随时间变化的虚拟机器(图3(a))来代替.同样,下生产线可以通过移除机器m1和缓冲区b1,同时使用效率随时间变化的虚拟机器来构造.

最后,引入有限量生产运行下的辅助单机生产线(见图4).

图4 辅助单机生产线Fig.4 Auxiliary one-machine lines

由于无论是m1还是m2,能够生产一个工件的前提条件都是当且仅当它处于工作状态且不被阻塞,同时无论或者是,能够生产的条件是当且仅当其处于工作状态且不会饥饿,我们定义随时间变化的辅助单机机器效率如下:

为了分析辅助单机生产线,注意,它们每个都可由一个马尔科夫链来表征,其中,系统状态为已被这台机器加工过的工件数量(参阅文献[17]).让,其中表示在时隙n结束时已经加工了j个工件的概率.的演化可以通过以下线性时变方程给出:

其中初始状态是

其中初始状态是

综上,为了分析图1中的有限量运行下的三机装配线的暂态性能,我们将原始系统的动态特性进行分解和简化,通过分析一系列分解后相互影响的动态特性更加简单的系统,来近似评估原始系统的实时性能.具体来说,对于原始系统(图1),其动态特性包括两方面:缓冲区中在制品数量的演化和在每台机器上已完成加工处理的工件数量.首先引入使用原始系统机器和缓冲区参数的辅助装配系统,同时假设无限原材料(图2).在这个系统中,我们只关注系统中的缓冲区在制品数量的演化.为了分析图2所示系统,进一步引入辅助双机串行线(图3),其中,为了考虑移除相应机器和缓冲区所带来的影响,上生产线和下生产线中装配机器所在位置分别使用相应的参数时变的虚拟机器来替代.因此,通过分析辅助双机串行线(图3),事实上可以得到辅助装配系统(图2)中系统状态(缓冲区在制品数量)的实时分布情况.最后,引入辅助单机生产线(图4)来分析在相应机器上完成加工处理工件数量的动态特性.而每一台单机生产线的时变参数都是在考虑了辅助双机串行线中的系统状态的影响下,近似推导得出的.

3.2 性能评估近似公式

基于上述构造的辅助生产线或生产系统,我们提出了近似原系统性能指标的计算公式.首先,有限量生产运行下一个批次的生产完成时间通过使用辅助虚拟单机线中的任意一个来近似估算.不失一般性,使用,同时令表示原始系统中机器m0在时隙n结束时处理加工完整个批次所有工件的近似概率.那么,

其次,原系统的生产率和各个零件生产线的消耗率可由辅助单机生产线的生产率来近似:

为了估算WIPi(n),BLi(n)和ST0,i(n),两种辅助生产线需要结合起来.具体来说,这些性能评估使用相应的辅助两机生产线来近似估算,同时考虑相对应的机器在辅助单机生产线上还没有完成加工整个批次所有产品的概率:

最后,一个批次的完成时间期望可以被近似为:

其中,T满足以下条件:

综上,基于分解的计算方法流程图如图5所示.

图5 分解算法流程图Fig.5 Flow chart of the calculation procedure

对于所提出的性能近似的方法的精确程度,我们通过对10000条参数随机而均匀地从式(21)所示的集合或者区间中选取的三机伯努利装配系统,进行基于精确解析和基于分解的近似性能评估分析,来验证所提近似方法的精确性.

对于每一条参数随机产生的装配系统,我们分别通过精确分析式(8)和基于分解的近似分析式(17)～(20)来计算其各项性能指标.结果显示,对于这10000条装配系统的各项性能指标的平均相对误差,它们的中值都在1%以下.

作为一个例子,考虑图6所显示的装配系统.每个机器(圆形表示)上的数字表示其效率,而每个缓冲区(矩形)中的数字表示其容量.这些参数是随机生成的.在本例中,所有缓冲区都被假设在起始状态时是空的.首先需要注意的是,使用精确分析方法,根据式(2),系统的状态数量为1620;而经过分解后,我们只需要分析六个相对较小但相互影响的系统:一条双机上生产线,一条双机下生产线,两条上单机生产线,两条下单机生产线.六个较小的马尔科夫链的总状态数为333.在保证精确度的基础上,相较精确分析,基于分解的近似分析使系统状态数量有了极大的降低.与此同时,从计算时间的角度来看,使用MATLAB软件在同一台电脑配置为因特尔酷睿i7-6700的CPU和16GB的RAM上,基于精确分析和基于分解的近似分析,所需要的运算时间分别为13.35秒和0.11秒,近似算法在计算高效性上也显示出了极大的优势.系统的暂态性能如图7所示,从图中可以看出,整个生产运行过程分为三个阶段.在第一阶段,产品开始进入空系统.在此期间,生产率和在制品数量都从0上升到稳态值.同时,由于更多的工件进入系统,零件生产线1(或者零件生产线2)的消耗率从p1(或者p2)开始逐渐减小.在第二阶段中,系统运行接近稳定状态,所有暂态性能指标都或多或少地处于平稳状态.最后,当生产运行接近完成时,所有性能指标开始下降,最终达到0.基于该分解算法的高精度也可以从图中清晰地看到.需要注意的是,虽然精确的分析在这种小型装配系统中仍然可以被推导出来,然而随着系统参数的增长,精确分析也变得越来越不可能实现.基于分解思想的性能近似评估方法的计算高效性将在这样的大型装配系统中体现出来.深入的相关研究将在未来的工作中被进一步讨论.

图6 三机伯努利装配系统的数值实例Fig.6 Example of an assembly system with three Bernoulli machines

4 结束语

本文研究了具有三台伯努利机器,有限缓冲区容量和有限量生产运行下的装配系统的暂态性能评估问题.具体地,首先推导了系统性能评价的精确数学模型和解析公式.然后,提出了一种基于分解的性能评估算法,通过将系统转换成一系列相互作用的辅助串行线来近似评估原始系统的暂态性能.论文推导了基于分解的三机装配系统实时性能估计公式,并通过数值实验验证了算法的准确性和计算高效性.

图7 分解近似与精确分析的三机伯努利装配系统暂态实时性能评估对比Fig.7 Comparison of decomposition-based approxiamtion and exact analysis for transientperformance evaluation in assembly system with three Bernoulli machines

今后在这方面的工作包括将算法扩展到每个零件生产线具有多台机器和多个缓冲区的系统,或多条零件生产线和多装配操作的复杂装配系统.此外,还会将研究结果推广到具有其他机器可靠性模型(几何型、指数型、威布尔型等)的装配系统中.