CN105956295A - 批处理流水线性能评估方法 - Google Patents

Abstract

本发明提出一种批处理流水线性能评估方法，用于评估批处理流水线平均产出和各缓冲区水平等指标。批处理机流水线性能评估问题，主要存在的问题是：“单机+批处理机”的两机器精确模型的求解和批处理机流水线分解方程的建立。本发明具体过程为：1、求出“单机+批处理机”以及“单机+单机”的两机器模型的精确解；2、将含多机器的流水线分解为多个含两机器单缓冲的构件块，并针对这些构件块，建立一些分解方程；3、通过一种迭代方法(PDDX法)对这些方程进行求解，并最终求得该流水线的性能指标：生产率和在制品数量；4、通过和仿真实验结果进行数值对比，验证本发明的合理性。

Description

批处理流水线性能评估方法

技术领域

本发明涉及生产系统性能评估领域，具体为一种批处理流水线性能评估方法，用于评估批处理流水线平均产出和各缓冲区水平指标。

背景技术

对生产系统进行严谨的性能分析、持续改进和精益设计可以给企业带来很大经济效益，研究意义重大。关于不可靠机器和有限缓冲容量缓冲区的流水线系统性能评估问题，已取得很多研究成果，这些研究通常都假设一台机器一次只能加工一个工件。而在实际情况中，也存在一台机器一次需要加工一批，即多个工件的情况。如齿轮制造的过程中，经过铣齿，磨削后，要进行热处理(淬火和退火)增大强度，淬火和退火这两个过程都是批处理的。

解析和仿真方法都可以对生产系统指标进行评估，但仿真方法只能面对某个具体问题，不具有一般通用性，且成本高，花费时间长；而解析方法可以面向同一类问题，且能够提供较为一致的解决方案。解析法包括精确分析方法和近似分析方法：1)精确分析方法适合简单的两工作站流水线性能评估，通过构建马尔科夫过程，求解得到系统状态稳态概率分布的精确解，并进一步得到性能评估指标；2)近似分析方法适合更复杂系统的性能分析，它是在简单系统的基础上进行递推迭代，主要有分解(Decomposition)和集结(Aggregation)两种方法。分解方法将原始生产系统分解为多个可以采用精确解析方法求解的子系统，根据各子系统间需满足的流失效方程、流修复方程、流加工方程等，构建迭代算法求解平均产出(Average Throughput)、平均缓冲水平(Average Buffer Level)等系统性能指标。分解子系统形式一般采用两机器流水线；集结方法采用了与分解方法相反的建模思路，将两机器生产单元近似为单个等效机器，并沿着串行生产线前向和后向递归进行，当聚合方法收敛时，能够获得系统的生产率等性能指标。

发明内容

要解决的技术问题

针对批处理机流水线性能评估问题，主要存在的问题是：“单机+批处理机”的两机器精确模型的求解和批处理机流水线分解方程的建立。

(1)“单机+批处理机”两机器精确模型求解

流水线分解模型是沿着物料流的流向，将初始m台机器m-1个缓冲的流水线，依次分解成m-1个两机器单缓冲构件块。针对含批处理机的特殊流水线，这些构件块由两类构成：一种是“单机+单机”的两机器模型，另一种是“单机+批处理机”的两机器模型。如何对“单机+批处理机”的两机器模型进行精确求解，是本发明技术方案要解决的技术问题之一。

(2)批处理机水线分解方程的建立

考虑到“单机+批处理机”构件块的特殊性，如何在上述两类器构件块之间，构造合理的分解方程是本发明技术方案要解决的另外一个重要技术问题。

技术方案

本发明提出了批处理流水线性能评估方法，并通过仿真对比实验对其合理性进行了说明。具体研究过程为：1、求出“单机+批处理机”以及“单机+单机”的两机器模型的精确解；2、将含多机器的流水线分解为多个含两机器单缓冲的构件块，并针对这些构件块，建立一些分解方程；3、通过一种迭代方法(PDDX法)对这些方程进行求解，并最终求得该流水线的性能指标：生产率和在制品数量；4、通过和仿真实验结果进行数值对比，验证本发明的合理性。

批处理流水线描述及假设：

考虑含m-1台单机，一台批处理机的流水线，如图1所示。其中M_i(i＝1,2…m-1)为单机，M_m为批处理机，B_i,(i∈[1,m-1])为M_i和M_i+1之间的缓冲。

我们使用以下约束来定义该系统的运行：

·物料流为离散的，且在加工和储存过程不会损失；

·所有机器具有相同的加工周期τ，且以一个加工周期为一段，将时间轴分段；

·所有机器为伯努利机器，机器M_i的独立效率为p_i，即在每个时间段开端，M_i处于工作状态的概率为p_i，处于故障状态的概率为1-p_i，与该机器之前任何时刻状态无关；

·每个缓冲B_i,(i∈[1,m-2])都具有一个容量N_i∈[1,+∞)，批处理机M_m的上游缓冲区B_m-1的容量为N_m-1∈[2,+∞)，在每个时间段的末端确定缓冲区的状态h_i∈[0,N_i]，令初始状态h_i＝0；

·上游缓冲区为空的机器称其处于饥饿(Starve)状态，下游缓冲区饱和的机器称其处于阻塞(Block)状态。第一台机器不会饥饿(原料充足)，最后一台机器不会被下游阻塞；

·批处理机M_m一次固定加工的工件个数为k∈[2,N_m-1]，当上游缓冲区工件个数不足k，即h_m-1＜k时，批处理机M_m会一直处于饥饿状态；

·缓冲区具有有限的缓冲区容量，阻塞机制采用服务前阻塞(Blocking BeforeService,BBS)的方式，即工件在被机器提取加工时已被考虑置于下游缓冲区内，若阻塞则停止加工；

·所有机器相互独立；

·机器的故障是与时间相关的故障(Time Dependent Failures,TDFs)。机器在阻塞或饥饿时仍然可能发生故障。

在时刻n，该流水线性能指标的计算公式为：

·缓冲B_i的在制品库存水平：

${\mathrm{WIP}}_i\left(n\right)=\underset{j=1}{\overset{N_i}{\Σ}}jP\[h_i\left(n\right)=j\].---\left(1\right)$

P[h_i(n)＝j]表示n时刻第i个缓冲区状态h_i(n)＝j的概率。

·缓冲B_i的饥饿率：

$S_i\left(n\right)=\left\{\begin{array}{c}P\[h_i\left(n\right)=0\],1\≤i\≤m-2\\ P\[h_{m-1}\left(n\right)\≤k-1\].\end{array}\right.---\left(2\right)$

·缓冲B_i的阻塞率：

B_i(n)＝P[h_i(n)＝N_i](1-p_i+1),1≤i≤m-1 (3)

·机器M_i的生产率：

${\mathrm{PR}}_i\left(n\right)=\left\{\begin{array}{c}{p}_i(1-S_{i-1}(n-1\left)\right),2\≤i\≤m-1,\\ {\mathrm{kp}}_m(1-S_{m-1}(n-1\left)\right).\end{array}\right.---\left(4\right)$

·机器M_i的效率：

$\begin{array}{c}{E}_i\left(n\right)=p_i(1-S_{i-1}(n-1)-B_i(n-1\left)\right),1\≤i\≤m\\ S_0\left(n\right)=B_n\left(n\right)=0\end{array}---\left(5\right)$

1.两机器精确模型求解

1.1“单机+批处理机”两机器精确模型求解

考虑由一台单机M_m-1和一台批处理机M_m组成的两机器构件块如图2所示，二者之间的缓冲容量为N_m-1∈[2,+∞)，批处理机M_m一次固定加工的工件个数为k∈[2,N_m-1]。该系统的运行满足上述假设条件。系统的状态空间共由N_m-1+1个状态组成：0,1,…N_m-1。系统的状态转移图如3所示，其中实线表示缓冲区的占用量在每个时间段的最大变化量不超过1个工件，虚线表示缓冲区的占用量在每个时间段的最大变化量不超过k个工件。

1)具有不变效率的“单机+批处理机”模型求解

假设机器M_i,(i＝m-1,m)在时刻n的效率为p_i,(i＝m-1,m)，如图2(a)所示，令x_η(n),η∈[0,N_m-1]表示系统在时间段n时处于状态η的概率，在这里x_η(n)＝P[h_m-1(n)＝η]，表示马尔科夫链的概率分布。x(n)随时间的演化过程可以通过下面带约束的线性动态系统来描述：

x(n+1)＝Ax(n)

$\underset{i=0}{\overset{N_{m-1}}{\Σ}}{x}_i(n+1)=1---\left(6\right)$

转移概率矩阵A定义如下：

2)具有可变效率的“单机+批处理机”模型求解

假设机器M_i,(i＝m-1,m)在时刻n的效率为p_i(n),(i＝m-1,m)，如图2(b)所示，则x(n)随时间的演化过程可以描述为：

$x\left(n\right)=A\left(p_{m-1}\right(n),p_m(n),N_{m-1})x(n-1),\underset{i=0}{\overset{N_{m-1}}{\Σ}}{x}_i\left(n\right)=1.---\left(8\right)$

其中，A(p_m-1(n),p_m(n),N_m-1)可以通过将公式(7)中的p_i替换为p_i(n)得到。

“单机+批处理机”在时刻n的性能评估计算公式为：

$\begin{array}{c}PR\left(n\right)=C_{1}x\left(n\right)=k\[Z_0,p_m\left(n\right)Z_1^{\′}\]x\left(n\right),\\ WIP\left(n\right)=C_{2}x\left(n\right)=\[0,1,2,\mathrm{...},N_{m-1}\]x\left(n\right),\\ S\left(n\right)=C_{3}x\left(n\right)=\[Z_0^{\′},Z_1\]x\left(n\right),\\ B\left(n\right)=C_{4}x\left(n\right)=\[0,\mathrm{...0},1\]x\left(n\right).\end{array}---\left(9\right)$

其中，Z₀表示1×k阶零矩阵；Z′₀表示1×k阶单位阵；Z₁表示1×(N_m-1+1-k)零矩阵；Z′₁表示1×(N_m-1+1-k)单位阵。

1.2“单机+单机”两机器精确模型求解

对于“单机+单机”模型如图4所示，可通过类似方法求得具有不变效率的精确解，及具有变动效率的精确解。

1)具有不变效率的“单机+单机”模型求解

假设机器M_i,(i＝1,2)在时刻n的效率为p_i,(i＝1,2)，如图4(a)所示，系统状态x(n)随时间的演化过程可以通过下面带约束的线性动态系统来描述：

x(n+1)＝A₁x(n),

$\underset{i=0}{\overset{N_1}{\Σ}}{x}_i(n+1)=1.---\left(10\right)$

其中，

2)具有可变效率的“单机+单机”模型求解

假设机器M_i,(i＝1,2)在时刻n的效率为p_i(n),(i＝1,2)，如图4(b)所示，则x(n)随时间的演化过程可以描述为：

$x\left(n\right)=A_1\left(p_1\right(n),p_2(n),N_1)x(n-1),\underset{i=0}{\overset{N_i}{\Σ}}{x}_i\left(n\right)=1.---\left(12\right)$

其中，A₁(p₁(n),p₂(n)_,N₁)可以通过将式(11)中的p_i替换为p_i(n)得到。

“单机+单机”在时刻n的性能评估计算公式为：

$\begin{array}{c}{\mathrm{PR}}_1\left(n\right)=C_1^{\′}x\left(n\right)=\[0,p_2\left(n\right)D_1\]x\left(n\right),\\ {\mathrm{WIP}}_1\left(n\right)=C_2^{\′}x\left(n\right)=\[0,1,2,\mathrm{...}{N}_1\]x\left(n\right),\\ S_1\left(n\right)=C_3^{\′}x\left(n\right)=\[1,D_0\]x\left(n\right),\\ B_1\left(n\right)=C_4^{\′}x\left(n\right)=\[D_0,1\]x\left(n\right).\end{array}---\left(13\right)$

其中，D₀表示1×N₁零矩阵；D₁表示1×N₁单位阵。

2.批处理流水线分解模型

将如图1所示含m-1台单机、一台批处理机的流水线l，分解为m-1个两机器构件块l(i)，其中l(i),(i∈[1,m-2])表示具有可变效率的“单机+单机”构件块，l(m-1)表示具有可变效率的“单机+批处理机”构件块，如图5所示。以下对构件块中参数和原始流水线中参数以括号加以区分。

每个构件块l(i)由缓冲B(i)、上游虚拟机器M^u(i)和下游虚拟机器M^d(i)构成。B(i)为B_i的复制品，即它们具有相同的缓冲容量、饥饿率和阻塞率。M^d(i-1)和M^u(i)由同一台真实机器M_i分解得到，分别表示M_i的流入和流出。上游虚拟机器的可变独立效率参数为p^u(i,n)，下游虚拟机器的可变独立效率参数为p^d(i,n)。构件块l(i)在时刻n的状态定义如下：

·上游虚拟机器效率E^u(i,n)；

·下游虚拟机器效率E^d(i,n)；

·缓冲B(i)饥饿率S_i(n)；

·缓冲B(i)阻塞率B_i(n)。

由于M^d(i-1)和M^u(i)由同一台真实机器M_i分解得到，所以在时刻n，三者具有相同的效率。即：

E_i(n)＝E^u(i,n)＝E^d(i-1,n),1≤i≤m (14)

其中，迭代公式为：

E_i(n)＝p_i(1-S_i-1(n-1)-B_i(n-1)),1≤i≤m

E^u(i,n)＝p^u(i,n)(1-B_i(n-1)),1≤i≤m-1

E^d(i-1,n)＝p^d(i-1,n)(1-S_i-1(n-1)),2≤i≤m

S₀(n-1)＝B_m(n-1)＝0

从而可得

$p^u(i,n)=\frac{E^u(i,n)}{1-B_i(n-1)}=\frac{E_i\left(n\right)}{\frac{E_i\left(n\right)}{p_i}+S_{i-1}(n-1)}=\frac{1}{\frac{1}{p_i}+\frac{S_{i-1}(n-1)}{E^d(i-1,n)}},i=2,\mathrm{3...}m-1---\left(15\right)$

$p^d(i-1,n)=\frac{E^d(i-1,n)}{1-S_{i-1}(n-1)}=\frac{E_i\left(n\right)}{\frac{E_i\left(n\right)}{p_i}+B_i(n-1)}=\frac{1}{\frac{1}{p_i}+\frac{B_i(n-1)}{E^u(i,n)}},i=m-1,\mathrm{...2}---\left(16\right)$

边界条件为

$\begin{array}{c}{p}^u(1,n)\≡p_1,n=0,1,\mathrm{...}\\ p^d(m-1,n)\≡p_m,n=0,1,\mathrm{...}\end{array}---\left(17\right)$

3.批处理流水线分解方程的求解(PDDX法)

1)初始化

考虑在时刻n＝0，将构件块l(i),(i＝1…m-1)中虚拟机器的参数初始化为原始流水线中真实机器相应的参数；

p^u(i,n)＝p_i,i＝1,2…m-1,

p^d(i,n)＝p_i+1,i＝1,2…m-1. (18)

边界条件：

$\begin{array}{c}{p}^u(1,n)\≡p_1,n=0,1,\mathrm{...}\\ p^d(m-1,n)\≡p_m,n=0,1,\mathrm{...}\end{array}---\left(19\right)$

然后，根据公式(9)和(13)，计算各个构件块在时刻n＝0的阻塞率B_i(0)、饥饿率S_i(0)、在制品库存水平WIP_i(0)和生产率PR_i(0)等指标。

2)令n＝n+1，依次更新构件块l(i),(i＝2…m-1)上游机器参数

将构件块l(i-1)的最新评估结果代入下述等式，更新构件块l(i)上游机器M^u(i)的参数，并利用公式(9)和(13)，对构件块l(i)重新进行评估。

$p^u(i,n)=\frac{1}{\frac{1}{p_i}+\frac{S_{i-1}(n-1)}{E^d(i-1,n)}},i=2...m-1.---\left(20\right)$

3)依次更新构件块l(i-1),(i＝m-1,…2）下游机器参数

将构件块l(i)的最新评估结果代入下述等式，更新构件块l(i-1)下游机器M^d(i-1)的参数，并利用公式(9)和(13)，对构件块l(i-1)重新进行评估。

$p^d(i-1,n)=\frac{1}{\frac{1}{p_i}+\frac{B_i(n-1)}{E^u(i,n)}},i=m-1,...2.---\left(21\right)$

4)收敛条件

令当Δ(n)＜10^-8迭代算法终止；否则，回到步骤2)。

5)输出结果

当迭代结果收敛时，输出该流水线的平均产出和缓冲区平均缓冲水平性能指标值。

综上所述，对如图1所示含m-1台单机、一台批处理机的流水线l进行性能评估步骤如下：

Step1 流水线分解

将批处理流水线l，依次分解为m-1个两机器构件块l(i),(i∈[1,m-1])，如图5所示；其中l(i),(i∈[1,m-2])表示具有可变效率的“单机+单机”构件块，l(m-1)表示具有可变效率的“单机+批处理机”构件块；每个构件块l(i)由缓冲B(i)、上游虚拟机器M^u(i)和下游虚拟机器M^d(i)构成；

Step2 参数初始化

Step2.1 根据公式

p^u(i,n)＝p_i,i＝1,2…m-1,n＝0

p^d(i,n)＝p_i+1,i＝1,2…m-1,n＝0

p^u(1,n)≡p₁,n＝0,1,…

p^d(m-1,n)≡p_m,n＝0,1,…

将构件块l(i),(i∈[1,m-1])中上下游机器参数初始化，p^u(i,n)为上游虚拟机器的可变独立效率参数，p^d(i,n)为下游虚拟机器的可变独立效率参数；其中，p_i(n),i＝1,2…m表示真实第i个机器的效率；

Step2.2 根据公式

PR(n)＝C₁x(n)＝k[Z₀,p_m(n)Z′₁]x(n),

WIP(n)＝C₂x(n)＝[0,1,2,…,N_m-1]x(n),

S(n)＝C₃x(n)＝[Z′₀,Z₁]x(n),

B(n)＝C₄x(n)＝[0,…0，1]x(n).

计算“单机+批处理机”构件块在时刻n＝0的阻塞率B_i(0)、饥饿率S_i(0)、在制品库存水平WIP_i(0)和生产率PR_i(0)，其中x(n)表示马尔科夫链的概率分布，Z₀表示1×k阶零矩阵；Z′₀表示1×k阶单位阵；Z₁表示1×(N_m-1+1-k)零矩阵；Z′₁表示1×(N_m-1+1-k)单位阵；

根据公式

PR₁(n)＝C′₁x(n)＝[0,p_i(n)D₁]x(n),i＝2,3,…,m-1

WIP₁(n)＝C'₂x(n)＝[0，1，2…N_i-1]x(n),

S₁(n)＝C′₃x(n)＝[1,D₀]x(n),

B₁(n)＝C'₄x(n)＝[D₀,1]x(n).

计算第i-1个“单机+单机”构件块在时刻n＝0的阻塞率B_i(0)、饥饿率S_i(0)、在制品库存水平WIP_i(0)和生产率PR_i(0)，i＝2,3,…,m-1，其中，D₀表示1×N₁零矩阵；D₁表示1×N₁单位阵。

Step3 递归迭代

Step3.1 向上迭代，更新上游机器参数：

令n＝n+1，利用公式

$p^u(i,n)=\frac{1}{\frac{1}{p_i}+\frac{S_{i-1}(n-1)}{E^d(i-1,n)}},i=2...m-1.$

依次更新构件块l(i),(i＝2…m-1)上游机器参数，其中

E^d(i-1,n)＝E_i(n)＝p_i(1-S_i-1(n-1)-B_i(n-1))

并利用公式

PR(n)＝C₁x(n)＝k[Z₀,p_m(n)Z′₁]x(n),

WIP(n)＝C₂x(n)＝[0,1,2,…,N_m-1]x(n),

S(n)＝C₃x(n)＝[Z′₀,Z₁]x(n),

B(n)＝C₄x(n)＝[0,…0，1]x(n).

计算“单机+批处理机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

利用公式

PR₁(n)＝C′₁x(n)＝[0,p_i(n)D₁]x(n),i＝2,3,…,m-1

WIP₁(n)＝C'₂x(n)＝[0，1，2…N_i-1]x(n),

S₁(n)＝C′₃x(n)＝[1,D₀]x(n),

B₁(n)＝C'₄x(n)＝[D₀,1]x(n).

计算第i个“单机+单机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

Step3.2向下迭代，更新下游机器参数

利用公式

$p^d(i-1,n)=\frac{1}{\frac{1}{p_i}+\frac{B_i(n-1)}{E^u(i,n)}},i=m-1,...2.$

依次更新构件块l(i-1),(i＝m-1,…2)下游机器参数，其中

E^u(i,n)＝E_i(n)＝p_i(1-S_i-1(n-1)-B_i(n-1))

并利用公式

PR₁(n)＝C′₁x(n)＝[0,p_i(n)D₁]x(n),i＝2,3,…,m-1

WIP₁(n)＝C'₂x(n)＝[0，1，2…N_i-1]x(n),

S₁(n)＝C′₃x(n)＝[1,D₀]x(n),

B₁(n)＝C'₄x(n)＝[D₀,1]x(n).

计算第i-1个“单机+单机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

Step4 收敛性

令当Δ(n)＜10-⁸迭代算法终止；否则，返回Step3

Step5 输出结果

当迭代结果收敛时，输出批处理流水线的平均生产率和缓冲区平均在制品库存水平。

有益效果

针对含批处理机流水线性能评估问题，本发明采用将PDDX结果与PlantSimulation仿真结果对比的方式验证本发明的有效性。

具体地，针对同一流水线分别采用本发明提出的PDDX、Plant Simulation仿真模型两种方法获得流水线线系统平均生产率、各缓冲区平均在制品库存水平等性能指标值，并以仿真结果为基准，计算PDDX方法的系统性能指标偏差百分比，作为流水线分解模型有效性的评价指标，系统性能指标偏差百分比计算公式为

$\begin{array}{c}{\∈}_{PR}=\left(\frac{{\mathrm{PR}}_{PDDX}-{\mathrm{PR}}_{Simulation}}{{\mathrm{PR}}_{Simulation}}\right)\×100\%\\ {\∈}_{B_i}=\left(\frac{B_{i_{PDDX}}-B_{i_{Simulation}}}{N_i}\right)\×100\%\end{array}---\left(22\right)$

为了保证有效性分析的客观和全面，申请人分别对含一台批处理机的两机器流水线、三机器流水线在不同参数下共16组实验的结果和用Plant Simulation8.2仿真结果进行对比，对所提算法进行性能分析。

两机器流水线线和三机器流水线实验参数如表1～表3所示，针对每组实验，批处理流水线分解模型根据PDDX方法中的收敛条件，用matlab软件运行至收敛，计算各缓冲区平均缓冲水平和系统平均产出；批处理流水线Plant Simulation仿真模型运行10days，统计各缓冲区平均缓冲水平和系统平均产出，所求得的结果如表4～表6所示。

表1两机器流水线各机器参数(一)

表3三机器流水线各机器参数

表4两机器批处理流水线的PDDX算法与仿真算法对比结果(一)

表5两机器批处理机流水线的PDDX算法与仿真算法对比结果(二)

表6三机器批处理流水线的PDDX算法与仿真算法对比结果

将仿真结果作为真实结果，代入公式(22)对比发现，由PDDX方法获得的缓冲区平均缓冲水平和系统平均产出满足：

1)由本专利的方法所得平均产出接近仿真结果，在16组试验中最大误差不超过3％左右，且仅有两组大于1％，其余平均产出与仿真结果偏差均小于1％；

2)对缓冲区在制品数量的估计精确度稍微低一些，在16组对比实验中与仿真结果偏差最大误差不超过15％；

虽然解析和仿真方法都可以对这16组生产系统指标进行评估，但仿真方法一次只能面对一个具体问题，不具有一般通用性，且成本高，花费时间长；而解析方法可以面向这一类问题，且能够提供较为一致的解决方案。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1：批处理流水线；

图2：“单机+批处理机”模型

a)具有不变效率的“单机+批处理机”模型；

b)具有可变效率的“单机+批处理机”模型；

图3：“单机+批处理机”状态转移图；

图4：“单机+单机”模型

a)具有不变效率的“单机+单机”模型；

b)具有可变效率的“单机+单机”模型；

图5：批处理流水线分解模型；

图6：三机器流水线Matlab实验结果。

具体实施方式

下面详细描述本发明的实施例，描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

本实施例以表6实验1的算例对本专利提出的批处理流水线的性能评估方法(PDDX)进行说明：

步骤1：流水线分解：

将含两台单机、一台批处理机的三机器批处理流水线l，分解为一个“单机+单机”构件块l(1)和一个“单机+批处理机”构件块l(2)。

步骤2：参数初始化：

构件块中虚拟机器的参数初始化定义为原始流水线中真实机器相应的参数(表3)。

l(1)构件块参数：

p^u(1,0)＝p₁＝0.9

p^d(1,0)＝p₂＝0.85

l(2)构件块参数：

p^u(2,0)＝p₂＝0.85

p^d(2,0)＝p₃＝0.8

根据公式

PR₁(n)＝C′₁x(n)＝[0,p_i(n)D₁]x(n),i＝2,3,…,m-1

WIP₁(n)＝C'₂x(n)＝[0，1，2…N_i-1]x(n),

S₁(n)＝C′₃x(n)＝[1,D₀]x(n),

B₁(n)＝C'₄x(n)＝[D₀,1]x(n).

计算“单机+单机”构件块l(1)在n＝0的评估结果

PR₁(0)＝0.8346

B₁(0)＝0.0727

S₁(0)＝0.0181

WIP₁(0)＝2.2562

根据公式

PR(n)＝C₁x(n)＝k[Z₀,p_m(n)Z′₁]x(n),

WIP(n)＝C₂x(n)＝[0,1,2,…,N_m-1]x(n),

S(n)＝C₃x(n)＝[Z′₀,Z₁]x(n),

B(n)＝C₄x(n)＝[0,…0，1]x(n).

计算“单机+批处理机”构件块l(2)在n＝0的评估结果

PR₂(0)＝0.8499

B₂(0)＝1.57×10^-4

S₂(0)＝0.4688

WIP₂(0)＝1.6086

利用“单机+单机”构件块l(1)在n＝0的评估结果，更新构件块l(2)上游机器的参数，并对构件块l(2)重新进行评估，得到l(2)在n＝1的评估结果：

$p^u(2,1)=\frac{1}{\frac{1}{p_2}+\frac{S_1\left(0\right)}{E^d(1,1)}}=\frac{1}{\frac{1}{0.85}+\frac{0.0181}{0.8346}}=0.8346$

p^d(2,1)＝p^d(2,0)＝0.8

PR₂(1)＝0.8345

B₂(1)＝1.4287×10-⁴

S₂(1)＝0.4784

WIP₂(1)＝1.5867

利用将最新得到的“单机+批处理机”构件块l(2)的评估结果，更新构件块l(1)下游机器的参数，并对构件块l(1)重新进行评估，得到l(1)在n＝1的评估结果：

$p^d(1,1)=\frac{1}{\frac{1}{p_2}+\frac{B_2\left(0\right)}{E^u(2,1)}}=\frac{1}{\frac{1}{0.85}+\frac{1.57\×{10}^{-4}}{0.8346}}=0.8499,$

p^u(1,1)＝p^u(1,0)＝0.9.

PR₁(1)＝0.8345

B₁(1)＝0.0728

S₁(1)＝0.0181

WIP₁(1)＝2.2567

因为所以继续迭代计算，直至满足收敛性条件，采用MATLAB编程实现，最终得到结果如图6所示。

仿真求解过程如下：

对三机器装配模型，用Plant Simulation8.2仿真软件建模，各机器参数定义如下效率：p，平均修复时间(MTTR)：1/p，处理时间：t＝1。

表7三机器装配线各机器仿真参数

设置仿真模型运行10天，最终得到结果如表6所示。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (1)

1.一种批处理流水线性能评估方法，其特征在于：包括以下步骤：

步骤1：流水线分解：

将含m-1台单机、一台批处理机的批处理流水线l，依次分解为m-1个两机器构件块l(i),i∈[1,m-1]，其中l(i),i∈[1,m-2]表示具有可变效率的“单机+单机”构件块，l(m-1)表示具有可变效率的“单机+批处理机”构件块；每个构件块l(i)由缓冲B(i)、上游虚拟机器M^u(i)和下游虚拟机器M^d(i)构成；

步骤2：参数初始化：

步骤2.1：根据公式

p^u(i,n)＝p_i,i＝1,2…m-1,n＝0

p^d(i,n)＝p_i+1,i＝1,2…m-1,n＝0

p^u(1,n)≡p₁,n＝0,1,…

p^d(m-1,n)≡p_m,n＝0,1,…

在时刻n＝0，将构件块l(i),i∈[1,m-1]中上下游机器参数初始化，p^u(i,n)为上游虚拟机器的可变独立效率参数，p^d(i,n)为下游虚拟机器的可变独立效率参数；p_i(n),i＝1,2…m表示真实第i个机器的效率；

步骤2.2：根据公式

PR(n)＝C₁x(n)＝k[Z₀,p_m(n)Z′₁]x(n),

WIP(n)＝C₂x(n)＝[0,1,2,…,N_m-1]x(n),

$S\left(n\right)=C_{3}x\left(n\right)=\[Z_0^{\′},Z_1\]x\left(n\right),$

B(n)＝C₄x(n)＝[0,…0,1]x(n).

计算“单机+批处理机”构件块在时刻n＝0的阻塞率B_i(0)、饥饿率S_i(0)、在制品库存水平WIP_i(0)和生产率PR_i(0)，其中x(n)表示马尔科夫链的概率分布，Z₀表示1×k阶零矩阵；Z′₀表示1×k阶单位阵；Z₁表示1×(N_m-1+1-k)零矩阵；Z′₁表示1×(N_m-1+1-k)单位阵；单机M_m-1和批处理机M_m之间的缓冲容量为N_m-1∈[2,+∞)，批处理机M_m一次固定加工的工件个数为k∈[2,N_m-1]；

根据公式

${\mathrm{PR}}_1\left(n\right)=C_1^{\′}x\left(n\right)=\[0,p_i\left(n\right)D_1\]x\left(n\right),i=2,3,...,m-1$

${\mathrm{WIP}}_1\left(n\right)=C_2^{\′}x\left(n\right)=\[0,1,2...N_{i-1}\]x\left(n\right),$

$S_1\left(n\right)=C_3^{\′}x\left(n\right)=\[1,D_0\]x\left(n\right),$

$B_1\left(n\right)=C_4^{\′}x\left(n\right)=\[D_0,1\]x\left(n\right).$

计算第i-1个“单机+单机”构件块在时刻n＝0的阻塞率B_i(0)、饥饿率S_i(0)、在制品库存水平WIP_i(0)和生产率PR_i(0)，i＝2,3,…,m-1，其中，D₀表示1×N_i-1零矩阵；D₁表示1×N_i-1单位阵，所计算的两个单机之间的缓冲容量为N_i-1；

步骤3：递归迭代：

步骤3.1：向上迭代，更新上游机器参数：

令n＝n+1，利用公式

$p^u(i,n)=\frac{1}{\frac{1}{p_i}+\frac{S_{i-1}(n-1)}{E^d(i-1,n)}},i=2...m-1.$

依次更新构件块l(i),i＝2…m-1上游机器参数，其中

E^d(i-1,n)＝E_i(n)＝p_i(1-S_i-1(n-1)-B_i(n-1))

并利用公式

$PR\left(n\right)=C_{1}x\left(n\right)=k\[Z_0,p_m\left(n\right)Z_1^{\′}\]x\left(n\right),$

WIP(n)＝C₂x(n)＝[0,1,2,…,N_m-1]x(n),

$S\left(n\right)=C_{3}x\left(n\right)=\[Z_0^{\′},Z_1\]x\left(n\right),$

B(n)＝C₄x(n)＝[0,…0,1]x(n).

计算“单机+批处理机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

利用公式

${\mathrm{PR}}_1\left(n\right)=C_1^{\′}x\left(n\right)=\[0,p_i\left(n\right)D_1\]x\left(n\right),i=2,3,...,m-1$

${\mathrm{WIP}}_1\left(n\right)=C_2^{\′}x\left(n\right)=\[0,1,2...N_{i-1}\]x\left(n\right),$

$S_1\left(n\right)=C_3^{\′}x\left(n\right)=\[1,D_0\]x\left(n\right),$

$B_1\left(n\right)=C_4^{\′}x\left(n\right)=\[D_0,1\]x\left(n\right).$

计算第i个“单机+单机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

步骤3.2：向下迭代，更新下游机器参数：

利用公式

$p^d(i-1,n)=\frac{1}{\frac{1}{p_i}+\frac{B_i(n-1)}{E^u(i,n)}},i=m-1,...2.$

依次更新构件块l(i-1),i＝m-1,…2下游机器参数，其中

E^u(i,n)＝E_i(n)＝p_i(1-S_i-1(n-1)-B_i(n-1))

并利用公式

${\mathrm{PR}}_1\left(n\right)=C_1^{\′}x\left(n\right)=\[0,p_i\left(n\right)D_1\]x\left(n\right),i=2,3,...,m-1$

${\mathrm{WIP}}_1\left(n\right)=C_2^{\′}x\left(n\right)=\[0,1,2...N_{i-1}\]x\left(n\right),$

$S_1\left(n\right)=C_3^{\′}x\left(n\right)=\[1,D_0\]x\left(n\right),$

$B_1\left(n\right)=C_4^{\′}x\left(n\right)=\[D_0,1\]x\left(n\right).$

计算第i-1个“单机+单机”构件块在新的n时刻下的阻塞率B_i(n)、饥饿率S_i(n)、在制品库存水平WIP_i(n)和生产率PR_i(n)；

步骤4：判断收敛性：

取当Δ(n)＜10^-8迭代算法终止，输出批处理流水线的平均生产率和缓冲区平均在制品库存水平；否则，返回步骤3。