# Performance Analysis of a Kitting Process as a Paired Queue

Authors: E. De Cuypere; K. De Turck; D. Fiems
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Nowadays, customers request more variation in a company’s product assortment leading to an increased amount of parts moving around on the shop floor. To cope with this tendency, a kitting process can be
implemented. Kitting is the operation of collecting the necessary parts for
a given end product in a specific container, called a kit, prior to arriving at
an assembly unit. As kitting performance is critical to the overall cost and
performance of the manufacturing system, this paper analyses a two-part
kitting process as a Markovian model. In particular, kitting is studied as a paired queue, thereby accounting for stochastic part arrivals, and kit assembly times. Using sparse matrix techniques, we assess the impact of kitting interruptions, bursty part arrivals and phase-type distributed kit assembly times on the behaviour of the part buffers. Finally, a cost-profit analysis of kitting processes is conducted and an approximation for a two-part kitting process is established.


---

## Body

## 1. Introduction

Nowadays manufacturing systems are often composed of multiple in-house fabrication units [1]. The semifinished products stemming from these units are the input materials for other fabrication units or for assembly lines. Hence, efficient transport of materials between the different stages of the production process is key for overall production cost minimisation. Kitting is a particular strategy for supplying materials to an assembly line. Instead of delivering in containers, each holding a single type part and all holding the same number of parts, kitting collects the necessary set of parts for an individual end product in a specific container, referred to as kit, prior to arriving at an assembly unit [1–6].Kitting mitigates storage space requirements at the assembly station since no part inventories need to be kept there. Moreover, parts are placed in proper positions in the container such that assembly time reductions can be realised. Additional benefits include reduced learning time of the workers at the assembly stations and increased quality of the product. Although kitting is a nonvalue adding activity, its application can reduce the overall materials handling time [6]. Indeed activities such as selecting and gripping parts are performed more efficiently. Furthermore, the whole operator walking time is drastically reduced or even eliminated since kits, each containing a complete set of components, are brought to the assembly station [7]. The advantages mentioned above do not come for free since the kitting operation itself incurs additional costs such as the time and effort for planning the allocation of the parts into kits and the kit preparation itself. Moreover, the introduction of a kitting operation in a production process involves a major investment and the effect on efficiency is uncertain. Therefore, it is important to analyse the performance of kitting in a production environment prior to its actual introduction. This is the subject of the present paper.In the literature, most authors consider a kitting process as a queueing system with stochastic part arrivals and kit assembly times. Hopp and Simon [8] developed a model for a kitting process with exponentially distributed processing times for kits and Poisson arrivals. They found accurate bounds for the required buffer capacity of kitting processes with two parts. Explicitly accounting for finite buffer capacities, Som et al. [3] further refined the results of Hopp and Simon.Of course real buffers always have a finite capacity, the capacity being constrained by the storage room. However, if the capacity is large enough, we can have a good approximation of a process with a finite capacity on the basis of a model with unlimited capacity. This means that there is always enough space for upcoming parts, which simplifies the analysis. Unfortunately, the assumption of an infinite buffer is not valid for kitting processes. If the capacity is assumed to be infinite, then the model will degrade to an unstable stochastic system. Harrison [9] showed for a multiple input generalisation of the GI/G/1 queue that it is necessary to impose a restriction on the size of the buffer to ensure stability in the operations of a kitting process. Under this assumption, the probability to have a certain long-term stock position is equal and independent of the current stock position. This was also demonstrated by Latouche [10] who studied waiting lines with paired customers. We can consider this analysis as an abstraction of a kitting process with two types of parts.In this work, we focus on a kitting process modulated by a Markovian environment. The introduction of this environment allows us to study kitting under more realistic stochastic assumptions: kitting interruptions, bursty part arrivals, phase-type distributed kit assembly times, and so forth. Our paper extends the results on kitting in a Markovian environment [11].The remainder of this paper is organised as follows. Section2 describes the kitting process at hand. In Section 3, Chapman-Kolmogorov equations are derived and their numerical solution is discussed. To illustrate our approach, Section 4 considers a number of numerical examples. In particular, we assess the impact of kitting interruptions, bursty part arrivals, and phase-type distributed kit assembly times on the behaviour of the part buffers. Then, a cost-profit analysis of kitting processes is conducted and an approximation for a two-part kitting process is established. Finally, conclusions are drawn in Section 5.


## 2. Model Description

In this paper, we study a two-queue kitting process, as depicted in Figure1. Each queue has a finite capacity—let Cℓ denote the capacity of buffer ℓ, ℓ=1,2—and models the inventory of parts of a single type. New parts arrive at the buffers and if both buffers are nonempty, a kit is assembled by collecting a part from each buffer. Hence, departures from the buffers are synchronised, and the buffers are paired. Operation of part buffers therefore considerably differs from other queueing systems.Figure 1
Kitting process: the buffers are on the left, and the triangularly shaped kitting process is on the right.Arrivals at both buffers are modelled by a Markovian arrival process and kit assembly is not instantaneous. For ease of modelling, it is assumed that there is a modulating Markov process, arrival and service rates depending on the state of this process. To be more precise, the kitting process is modelled as a continuous-time Markov process with state space𝒞1×𝒞2×𝒦, whereby 𝒞ℓ={0,…,Cℓ} for ℓ=1,2 and with 𝒦={1,2,…,K} being the state space of the modulating process. At any time, the state of the kitting process is described by the triplet [m,n,i], with m and n being the number of parts in the first and second buffer, respectively, and i being the state of the modulating process. We now describe the state transitions. (i)
The state of the modulating process can change when there are neither arrivals nor departures. Letαij denote the transition rate from state i to state j(i,j∈𝒦,i≠j) and let A denote the corresponding generator matrix.(ii)
The state of the modulating process may remain the same or may change when there is an arrival. Letλij(ℓ) denote the (marked) transition rate from state i to state j when there is an arrival at buffer ℓ, ℓ=1,2. Moreover, let Lℓ denote the corresponding generator matrix. Note that such marked transitions from state i to state i are allowed.(iii)
Analogously, the state of the modulating process may remain the same or may change when there is a departure (in each buffer). Letμij and M denote the corresponding transition rate and generator matrix, respectively.Summarising, arrivals at and departures from the buffers are described by the generator matrices A, L1, L2, and M. So far, no diagonal elements of A have been defined. To simplify notation, it will be further assumed that the diagonal elements are chosen such that the row sums of A+L1+L2+M are zero.The computational method employed here does not require any homogeneity of the generator matrices. When required by the applications at hand, intensities may depend on the buffer content. In this case, we introduce superscripts to make this dependence explicit. For example,M(m,n) denotes the generator matrix of state transitions with departure when there are m parts in buffer 1 and n parts in buffer 2. In addition, we use arguments for rates as we already used superscripts and subscripts to distinguish the arrival rates at the different queues. For example, λij(ℓ)(m,n) denotes the arrival rate at buffer ℓ=1,2 from state i to state j when there are m parts in buffer 1 and n parts in buffer 2.Example 1.
In the most basic setting, parts arrive at the buffers in accordance to an independent Poisson processes with rateλ1 and λ2 and kit assembly times are exponentially distributed with parameter μ. In this case, there is no need to have a modulating Markov process, and the state is completely described by the number of parts in each buffer, (m,n). We have
(1)M=[μ],L1=[λ1],L2=[λ2],A=[-λ1-λ2-μ].Example 2.
To account for burstiness in the arrival process of the parts at the different buffers, the modulating process allows the mitigation of the Poissonian arrival assumptions. We can replace the Poisson processes by a two-class Markovian arrival processes. Multiclass Markovian arrival processes allow for intricate correlation and can be efficiently characterised from trace data [12, 13]. As we have two types of arrivals, the Markovian arrival process is described by the generator matrix Λ1 of transitions with arrivals at buffer 1, the generator matrix Λ2 with arrivals at buffer 2, and the generator matrix Λ0 without arrivals. As usual, the diagonal elements of Λ0 are negative and ensure that the row sums of Λ0+Λ1+Λ2 are zero. Retaining exponentially distributed kit assembly times, we have
(2)L1=Λ1,L2=Λ2,A=Λ0-μI,M=μI.
Here I denotes the identity matrix.Example 3.
As for the arrival processes, the model at hand is sufficiently flexible to include phase-type distributed kit assembly times. The phase-type distribution is completely characterised by an initial probability vectorτ and the matrix T which corresponds to non-absorbing transitions [14]. Let t′=-Te′ be the column vector with the rates to the absorbing state and let f be a row vector with zero elements except the first one. Assuming Poisson arrivals in both buffers (with rates λ1 and λ2, resp.), we get the following matrices:
(3)L1(m,n)=λ1I(1-1{m=0,n>0})+λ1e′τ1{m=0,n>0},L2(m,n)=λ2I(1-1{m>0,n=0})+λ2e′τ1{m>0,n=0},A(m,n)=T1{m>0,n>0}-λ1I-λ2I,M(m,n)=t′τ1{m>1,n>1}+t′f(1-1{m>1,n>1}).
Here, it is assumed that the background state equals 1 if one of the buffers is empty. When service starts again, the background state is chosen in accordance with the probability vector τ.


## 3. Analysis

Having established the modelling assumptions and settled our notation, we now focus on the analysis of the kitting process.


### 3.1. Balance Equations

We aim to define a set of equations for the steady-state probability vector for the Markov process[Q1(t),Q2(t),S(t)], Qℓ(t) being the number of parts in buffer ℓ at time t and S(t) being the state of the background process at time t.Letπi(m,n)=limt→∞Pr[Q1(t)=m,Q2(t)=n,S(t)=i] be the steady-state probability to be in state [m,n,i] and let π(m,n) be the vector with elements πi(m,n), for i∈𝒦. Figure 2 shows a fragment of the transition rate diagram of the kitting model in state [m,n,i]. As mentioned above, two independent input streams arrive at the buffers with intensity λij(ℓ) and are processed into kits with intensity μij. Upon completion of a kit, the content of both buffers is decreased by 1. Note that we only show the transitions whereby the modulating Markov process remains in state i. Moreover, possible dependence of the transition rates on the buffer sizes is not indicated.Figure 2
Fragment of the transition rate diagram for state(m,n,i).Based on the transition rate diagram, we now derive the balance equations of the kitting process at hand.(i)
First, consider the case where both buffers are neither empty nor full,0<n<C1 and 0<m<C2. We have
(4)πi(m,n)(∑j=1K‍λij(1)(m,n)+∑j=1K‍λij(2)(m,n)+μij(m,n)+∑j=1,j≠iK‍αij(m,n))=∑j=1K‍πj(m-1,n)λji(1)(m-1,n)+∑j=1K‍πj(m,n-1)λji(2)(m,n-1)+∑j=1K‍πj(m+1,n+1)μji(m+1,n+1)+∑j=1,j≠iK‍πj(m,n)αij(m,n),
for i,j∈𝒦 or equivalently,
(5)π(m-1,n)L1(m-1,n)+π(m,n-1)L2(m,n-1)+π(m+1,n+1)M(m+1,n+1)+π(m,n)A(m,n)=0.(ii)
If buffer 1 is empty and buffer 2 is neither empty nor full (m=0 and 0<n<C2), we have
(6)πi(0,n)(∑j=1K‍λij(1)(0,n)+∑j=1K‍λij(2)(0,n)+∑j=1,j≠iK‍αij(0,n))=∑j=1K‍πj(0,n-1)λji(2)(0,n-1)+∑j=1K‍πj(1,n+1)μji(1,n+1)+∑j=1,j≠iK‍πj(0,n)αji(0,n),fori,j∈𝒦 or equivalently,
(7)π(0,n-1)L2(0,n-1)+π(1,n+1)M(1,n+1)+π(0,n)(A(0,n)+diag(M(0,n)e′))=0.(iii)
Similarly, if buffer 2 is empty and buffer 1 is neither empty nor full (n=0 and 0<m<C1), we have
(8)π(m-1,0)L1(m-1,0)+π(m+1,1)M(m+1,1)+π(m,0)(A(m,0)+diag(M(m,0)e′))=0.(iv)
If both buffers are empty (m=0 and n=0), we have
(9)π(1,1)M(1,1)+π(0,0)(A(0,0)+diag(M(0,0)e′))=0.(v)
If buffer 1 is empty and buffer 2 is full (m=0 and n=C2), we get
(10)π(0,C2-1)L2(0,C2-1)+π(0,C2)(A(0,C2)+diag(M(0,C2)e′+L2(0,C2)e′))=0.(vi)
Similarly, if buffer 1 is full and buffer 2 is empty (m=C1 and n=0), we have
(11)π(C1-1,0)L1(C1-1,0)+π(C1,0)(A(C1,0)+diag(M(C1,0)e′+L1(C1,0)e′))=0.(vii)
Finally, if both buffers are full (m=C1 and n=C2), we find
(12)π(C1-1,C2)L1(C1-1,C2)+π(C1,C2-1)L2(C1,C2-1)+π(C1,C2)(A(C1,C2)+diag(L1(C1,C2)e′+L2(C1,C2)e′))=0.


### 3.2. Performance Measures

Given the steady-state vectorsπ(m,n), we can now obtain a number of interesting performance measures for the kitting system. For ease of notation, let π(m,n)=π(m,n)e′ denote the probability to have m parts in buffer 1 and n parts in buffer 2. Moreover, let π(1)(m)=∑n‍π(m,n) and π(2)(n)=∑m‍π(m,n) denote the marginal probability vectors of the buffer 1 and 2 content, respectively, for the different background states. Finally, the probability mass functions of the buffer contents equal π(1)(m)=π(1)(m)e′ and π(2)(n)=π(2)(n)e′.The following performance measures are of interest.(i)
The mean buffer1 and 2 content is EQ1 and EQ2, respectively,
(13)EQ1=∑mC1‍π(1)(m)m,EQ2=∑nC2‍π(2)(n)n.(ii)
The variance of the buffer1 and 2 is VarQ1 and VarQ2, respectively,
(14)VarQ1=∑mC1‍π(1)(m)m2-(EQ1)2,VarQ2=∑nC2‍π(2)(n)n2-(EQ2)2.(iii)
The effective load of the systemρeff is the amount of time that kitting is ongoing. As kitting is only ongoing when none of the buffers is empty, we have
(15)ρeff=1-π(1)(0)-π(2)(0)+π(0,0).(iv)
Let the throughputη be defined as the number of kits departing from the system per time unit. Taking into account all possible states from which we can have a departure, we find
(16)η=∑m=1C1‍∑n=1C2‍π(m,n)M(m,n)e′.(v)
The shortage probabilityK is the probability that one of the buffers is empty as
(17)K=π(1)(0)+π(2)(0)-π(0,0).(vi)
The loss probabilitybi in buffer i is the probability that an arriving part cannot be stored in buffer i, i=1,2. By noting that the accepted arrival load equals the departure rate, we find
(18)bi=ρi-ηρi,
where ρi is the arrival load of part i={1,2}. If the arrival load in both buffers is the same, then the loss probability is also the same: b1=b2. If the arrival load is not the same, then the excess load in the most loaded buffer is lost as well.


### 3.3. Methodology: Sparse Matrix Techniques

Queueing models for kitting processes are rather complicated. Indeed, the modelled kitting process has a multidimensional state space. Even for relative moderate buffer capacity, the multidimensionality leads to huge state spaces; this is the so-called state-space explosion problem. For many queueing systems, infinite-buffer assumptions may mitigate this problem. Given some buffer system with finite capacity, more efficient numerical routines can be constructed for the corresponding queueing system with infinite capacity. Unfortunately, as mentioned above, the infinite-buffer capacity assumption is not applicable for kitting processes and therefore cannot simplify the analysis. Recall that the infinite-capacity model is always unstable. For all input parameters except trivial ones (no arrivals), some or all of the queues grow unbounded with positive probability.Consequently, the multidimensionality of the state space and the inapplicability of the infinite-buffer assumption yield Markov processes with a finite but very large state space. However, the number of possible state transitions from any specific state is limited. This means that most of the entries in the generator matrix are zero; the matrix is sparse. As illustrated by the numerical examples, using sparse matrices and their associated specialised algorithms result in manageable memory consumption and processing times, compared to standard algorithms.The method used here to solve the sparse matrix equation is the projection method GMRES (generalized minimum residual). Details can be found in Stewart [15] (p. 194–197); Philippe et al. [16], Buchholz [17], and Saad and Schultz [18]. To solve the linear equation Ax=b, the GMRES algorithm approximates x by the vector xn∈Kn in a Krylov subspace Kn=span{b,Ab,…,An-1b} which minimises the norm of the residual ∥Axn-b∥. Since the residual norm is minimised at every step of the method, it is clear that it is nonincreasing. However, the work and storage requirement per iteration increases linearly with the iteration count. Hence, the cost of n iterations grows by O(n2) which is a major drawback of GMRES. This limitation is usually overcome by restarting the algorithm. After a chosen number of iterations m, the accumulated data is cleared and the intermediate results are used as the initial data for the next m iterations [18].To ensure fast convergence, it is also key to properly choose the initial vector that is passed on to the algorithm. We rely on MATLAB’s build in GMRES algorithm and assume a uniform initial vector if no additional information about the solution is known. However, this is often the case. Calculation speed can be further improved as, in practice, performance measures are not calculated for an isolated set of parameters. For example, when a plot is created, a parameter is varied over a range of values. In this case, a previously calculated steady-state vector for some set of parameters can be used as a first estimate of the steady-state vector for a new “perturbed” set of parameters. Using previously calculated steady-state vectors is trivial if the state spaces corresponding to the parameter sets are equal. In this case, the previously calculated steady-state vector can be passed on unmodified. If the state space changes, the steady-state vector must be rescaled to the new state space. In general, adding zero-probability states if the state space increases or removing states if the state space decreases turns out to be ineffective. This is easily explained by a simple example. Assume that we increase the queue capacity of one of the part buffers. Typically, even for moderate load, a considerable amount of probability mass can be found for queue size equal to capacity. Increasing the queue size and assigning zero probability to the new states is not a good estimate for the new steady-state vector. Also for the system with higher capacity, a considerable amount of probability mass can be found when the queue size equals the capacity (while zero probabilities were assigned).


## 3.1. Balance Equations

We aim to define a set of equations for the steady-state probability vector for the Markov process[Q1(t),Q2(t),S(t)], Qℓ(t) being the number of parts in buffer ℓ at time t and S(t) being the state of the background process at time t.Letπi(m,n)=limt→∞Pr[Q1(t)=m,Q2(t)=n,S(t)=i] be the steady-state probability to be in state [m,n,i] and let π(m,n) be the vector with elements πi(m,n), for i∈𝒦. Figure 2 shows a fragment of the transition rate diagram of the kitting model in state [m,n,i]. As mentioned above, two independent input streams arrive at the buffers with intensity λij(ℓ) and are processed into kits with intensity μij. Upon completion of a kit, the content of both buffers is decreased by 1. Note that we only show the transitions whereby the modulating Markov process remains in state i. Moreover, possible dependence of the transition rates on the buffer sizes is not indicated.Figure 2
Fragment of the transition rate diagram for state(m,n,i).Based on the transition rate diagram, we now derive the balance equations of the kitting process at hand.(i)
First, consider the case where both buffers are neither empty nor full,0<n<C1 and 0<m<C2. We have
(4)πi(m,n)(∑j=1K‍λij(1)(m,n)+∑j=1K‍λij(2)(m,n)+μij(m,n)+∑j=1,j≠iK‍αij(m,n))=∑j=1K‍πj(m-1,n)λji(1)(m-1,n)+∑j=1K‍πj(m,n-1)λji(2)(m,n-1)+∑j=1K‍πj(m+1,n+1)μji(m+1,n+1)+∑j=1,j≠iK‍πj(m,n)αij(m,n),
for i,j∈𝒦 or equivalently,
(5)π(m-1,n)L1(m-1,n)+π(m,n-1)L2(m,n-1)+π(m+1,n+1)M(m+1,n+1)+π(m,n)A(m,n)=0.(ii)
If buffer 1 is empty and buffer 2 is neither empty nor full (m=0 and 0<n<C2), we have
(6)πi(0,n)(∑j=1K‍λij(1)(0,n)+∑j=1K‍λij(2)(0,n)+∑j=1,j≠iK‍αij(0,n))=∑j=1K‍πj(0,n-1)λji(2)(0,n-1)+∑j=1K‍πj(1,n+1)μji(1,n+1)+∑j=1,j≠iK‍πj(0,n)αji(0,n),fori,j∈𝒦 or equivalently,
(7)π(0,n-1)L2(0,n-1)+π(1,n+1)M(1,n+1)+π(0,n)(A(0,n)+diag(M(0,n)e′))=0.(iii)
Similarly, if buffer 2 is empty and buffer 1 is neither empty nor full (n=0 and 0<m<C1), we have
(8)π(m-1,0)L1(m-1,0)+π(m+1,1)M(m+1,1)+π(m,0)(A(m,0)+diag(M(m,0)e′))=0.(iv)
If both buffers are empty (m=0 and n=0), we have
(9)π(1,1)M(1,1)+π(0,0)(A(0,0)+diag(M(0,0)e′))=0.(v)
If buffer 1 is empty and buffer 2 is full (m=0 and n=C2), we get
(10)π(0,C2-1)L2(0,C2-1)+π(0,C2)(A(0,C2)+diag(M(0,C2)e′+L2(0,C2)e′))=0.(vi)
Similarly, if buffer 1 is full and buffer 2 is empty (m=C1 and n=0), we have
(11)π(C1-1,0)L1(C1-1,0)+π(C1,0)(A(C1,0)+diag(M(C1,0)e′+L1(C1,0)e′))=0.(vii)
Finally, if both buffers are full (m=C1 and n=C2), we find
(12)π(C1-1,C2)L1(C1-1,C2)+π(C1,C2-1)L2(C1,C2-1)+π(C1,C2)(A(C1,C2)+diag(L1(C1,C2)e′+L2(C1,C2)e′))=0.


## 3.2. Performance Measures

Given the steady-state vectorsπ(m,n), we can now obtain a number of interesting performance measures for the kitting system. For ease of notation, let π(m,n)=π(m,n)e′ denote the probability to have m parts in buffer 1 and n parts in buffer 2. Moreover, let π(1)(m)=∑n‍π(m,n) and π(2)(n)=∑m‍π(m,n) denote the marginal probability vectors of the buffer 1 and 2 content, respectively, for the different background states. Finally, the probability mass functions of the buffer contents equal π(1)(m)=π(1)(m)e′ and π(2)(n)=π(2)(n)e′.The following performance measures are of interest.(i)
The mean buffer1 and 2 content is EQ1 and EQ2, respectively,
(13)EQ1=∑mC1‍π(1)(m)m,EQ2=∑nC2‍π(2)(n)n.(ii)
The variance of the buffer1 and 2 is VarQ1 and VarQ2, respectively,
(14)VarQ1=∑mC1‍π(1)(m)m2-(EQ1)2,VarQ2=∑nC2‍π(2)(n)n2-(EQ2)2.(iii)
The effective load of the systemρeff is the amount of time that kitting is ongoing. As kitting is only ongoing when none of the buffers is empty, we have
(15)ρeff=1-π(1)(0)-π(2)(0)+π(0,0).(iv)
Let the throughputη be defined as the number of kits departing from the system per time unit. Taking into account all possible states from which we can have a departure, we find
(16)η=∑m=1C1‍∑n=1C2‍π(m,n)M(m,n)e′.(v)
The shortage probabilityK is the probability that one of the buffers is empty as
(17)K=π(1)(0)+π(2)(0)-π(0,0).(vi)
The loss probabilitybi in buffer i is the probability that an arriving part cannot be stored in buffer i, i=1,2. By noting that the accepted arrival load equals the departure rate, we find
(18)bi=ρi-ηρi,
where ρi is the arrival load of part i={1,2}. If the arrival load in both buffers is the same, then the loss probability is also the same: b1=b2. If the arrival load is not the same, then the excess load in the most loaded buffer is lost as well.


## 3.3. Methodology: Sparse Matrix Techniques

Queueing models for kitting processes are rather complicated. Indeed, the modelled kitting process has a multidimensional state space. Even for relative moderate buffer capacity, the multidimensionality leads to huge state spaces; this is the so-called state-space explosion problem. For many queueing systems, infinite-buffer assumptions may mitigate this problem. Given some buffer system with finite capacity, more efficient numerical routines can be constructed for the corresponding queueing system with infinite capacity. Unfortunately, as mentioned above, the infinite-buffer capacity assumption is not applicable for kitting processes and therefore cannot simplify the analysis. Recall that the infinite-capacity model is always unstable. For all input parameters except trivial ones (no arrivals), some or all of the queues grow unbounded with positive probability.Consequently, the multidimensionality of the state space and the inapplicability of the infinite-buffer assumption yield Markov processes with a finite but very large state space. However, the number of possible state transitions from any specific state is limited. This means that most of the entries in the generator matrix are zero; the matrix is sparse. As illustrated by the numerical examples, using sparse matrices and their associated specialised algorithms result in manageable memory consumption and processing times, compared to standard algorithms.The method used here to solve the sparse matrix equation is the projection method GMRES (generalized minimum residual). Details can be found in Stewart [15] (p. 194–197); Philippe et al. [16], Buchholz [17], and Saad and Schultz [18]. To solve the linear equation Ax=b, the GMRES algorithm approximates x by the vector xn∈Kn in a Krylov subspace Kn=span{b,Ab,…,An-1b} which minimises the norm of the residual ∥Axn-b∥. Since the residual norm is minimised at every step of the method, it is clear that it is nonincreasing. However, the work and storage requirement per iteration increases linearly with the iteration count. Hence, the cost of n iterations grows by O(n2) which is a major drawback of GMRES. This limitation is usually overcome by restarting the algorithm. After a chosen number of iterations m, the accumulated data is cleared and the intermediate results are used as the initial data for the next m iterations [18].To ensure fast convergence, it is also key to properly choose the initial vector that is passed on to the algorithm. We rely on MATLAB’s build in GMRES algorithm and assume a uniform initial vector if no additional information about the solution is known. However, this is often the case. Calculation speed can be further improved as, in practice, performance measures are not calculated for an isolated set of parameters. For example, when a plot is created, a parameter is varied over a range of values. In this case, a previously calculated steady-state vector for some set of parameters can be used as a first estimate of the steady-state vector for a new “perturbed” set of parameters. Using previously calculated steady-state vectors is trivial if the state spaces corresponding to the parameter sets are equal. In this case, the previously calculated steady-state vector can be passed on unmodified. If the state space changes, the steady-state vector must be rescaled to the new state space. In general, adding zero-probability states if the state space increases or removing states if the state space decreases turns out to be ineffective. This is easily explained by a simple example. Assume that we increase the queue capacity of one of the part buffers. Typically, even for moderate load, a considerable amount of probability mass can be found for queue size equal to capacity. Increasing the queue size and assigning zero probability to the new states is not a good estimate for the new steady-state vector. Also for the system with higher capacity, a considerable amount of probability mass can be found when the queue size equals the capacity (while zero probabilities were assigned).


## 4. Numerical Results

With the balance equations at hand, we now illustrate our numerical approach by means of some examples.


### 4.1. Bursty Part Arrivals

As a first example, we quantify the impact of production inefficiency on the performance of a kitting process. To this end, we compare part buffers with Poisson arrivals to corresponding kitting systems with interrupted Poisson arrivals. The arrival interruptions account for inefficiency in the production process. Kit assembly times are assumed to be exponentially distributed with service rate equal to one, this value being independent of the number of parts in the different buffers. This is a kitting process with Markovian arrivals as described in Example2 in Section 2.The interrupted Poisson process considered here is a two-state Markovian process. In the active state, new parts arrive in accordance with a Poisson process with rateλ whereas no new parts arrive in the inactive state. Let α and β denote the rate from the active to the inactive state and vice versa, respectively. We then use the following parameters to characterise the interrupted Poisson process:
(19)σ=βα+β,κ=1α+1β,ρ=λσ.
Note that σ is the fraction of time that the interrupted Poisson process is active, the absolute time parameter κ is the average duration of an active and an inactive period, and ρ is the arrival load of the parts.Figure3 shows the mean number of parts in buffer 1 versus the arrival load, for various values of the buffer capacities C1 and C2 for Poisson arrivals (for both buffers) as well as for interrupted Poisson arrivals (again for both buffers). The arrival load is set to ρ=0.8 for all curves. In addition, we set σ=0.4 and κ=10 for the interrupted Poisson processes. Clearly, the mean buffer content increases as the arrival load increases as expected. Moreover, if more buffer capacity is available, it will also be used: the mean buffer content increases for increasing values of C1=C2. Comparing interrupted Poisson and Poisson processes, burstiness in the production process has a negative impact on performance—more buffering is required—if the queues are not fully loaded (ρ<1). As for ordinary queues, the opposite can be observed for overloaded buffers. In this case, the effect of a large burst is mainly reflected in loss and not in additional queue content. Burstiness also yields larger periods without arrivals during which the buffer size decreases.Figure 3
When the queues are not fully loaded, burstiness in the production process has a negative impact on the performance.By numerical examples, we can quantify expected buffer behaviour—for example, more production yields higher buffer content, higher buffer capacity mitigates the loss probability, and so forth. However, less trivial behaviour can be observed as well. Figure4 depicts the probability that the buffer is full versus the buffer capacity C1=C2. We compare performance of kitting with Poisson arrivals to kitting with interrupted Poisson arrivals at one buffer and at both buffers. As in the preceding figure, the interrupted Poisson processes are characterised by σ=0.4 and κ=10. As expected, the probability that the buffer is full decreases for increasing values of the buffer capacities. Moreover, to reduce this probability, more buffer capacity is required for the case of two interrupted Poisson processes than for the case of two Poisson processes. For the kitting process with one Poisson and one interrupted Poisson process, non-trivial performance results can be observed. Namely, interruptions in the production of a part more negatively affect buffer performance of the other part. Indeed, buffer 1 is full with higher probability if the arrivals at buffer 2 are interrupted than if the arrivals at buffer 1 are interrupted. First note that the loss probabilities in both buffers are the same. For the Poisson buffer, parts are rejected at the arrival rate ρ1 whenever the buffer is full. For the IPP buffer, parts are rejected when the buffer is full and the arrival process is in the on-state at rate ρ2/σ. As the loss probability in both buffers is the same, we have Pr[Q1=C1]ρ1=Pr[Q2=C2,S=1]ρ2/σ, or equivalently Pr[Q1=C1]=Pr[Q2=C2∣S=1]. As the second queue is more likely to be full when there are arrivals we have Pr[Q2=C2∣S=1]≥Pr[Q2=C2] which shows Pr[Q1=C1]≥Pr[Q2=C2].Figure 4
Irregularity in the production of a part leads to a higher probability to have a full buffer for the other part.


### 4.2. Phase-Type Distributed Kit Assembly Times

The second numerical example quantifies the impact of the distribution of the kit assembly times on kitting performance. In particular, we here study Erlang-distributed kit assembly times. Limiting ourselves to Poisson arrivals to both buffers, this numerical example fits Example3 of Section 2.Figures5 and 6 depict the mean number of parts in buffer 1 and the loss probability in buffer 1 for the kitting process and, as a reference point, for the M/E/1/N queue as well. In both figures, the arrival load is varied and different values of the variance of the kitting time distribution are assumed as indicated. The mean kitting time is equal to 1 for all curves and the capacity of both buffers is equal to 20. In underload (ρ<1), kitting performs worse than the M/E/1/N queue: the mean buffer content and the loss probability have a higher value. This follows from the fact that kitting stops when one of the buffers is empty. By increasing the load, it is obvious that the buffer content converges to the capacity and the loss probability to one. It is most interesting to observe that the shape of the service time distribution only has a small effect on these performance measures. Indeed, there is no significant performance difference when σ2 equals 1/4 and when it equals 1/8.Figure 5
Variation of the kitting time distribution does not have a significant impact on the mean number of parts.Figure 6
Variation of the kitting time distribution does not have a significant impact on the loss probability.


### 4.3. Cost and Profit Analysis

We now add a cost structure to the kitting process under study. In particular, cost and profit for the kitting systems of Sections4.1 and 4.2 are analysed.The proposed cost function is(20)c1(EQ1+EQ2)+c2K+c3(b1+b2),
where c1 is the holding cost of a part in the buffer, c2 is the shortage cost in one or in both of the buffers, and c3 is the loss cost. Note that for all figures, the input parameters are symmetric for both parts such that EQ1=EQ2 and b1=b2.


#### 4.3.1. Poisson Arrivals

In Figure7, the total cost of applying kitting systems versus the buffer capacity (varying from 1 to 30) is depicted. Limiting ourselves to Poisson arrivals for both buffers, ρi=λi=0.8 for part i={1,2}. We compare kitting performance for different costs as depicted. As expected, higher cost values lead to higher total cost. The cost structure also affects the optimal buffer capacity. When looking at the different costs separately, we observe that a higher holding cost decreases the optimal buffer capacity (from 15 to 12) and a higher loss cost increases the optimal capacity (from 12 to 14), whereas the optimal buffer capacity remains the same (12) for a higher shortage cost. Obviously, buffering is interesting when the storage cost is low and when the cost of rejecting parts (because of a full buffer) is high. It is most interesting to observe that, as the capacity increases, the cost models converge when the sum of the holding and the shortage cost (c1+c2) is equal. Indeed, as the loss tends to zero, the loss cost tends to zero as well. Furthermore, when the capacity equals one, the state space of the kitting model has size 4 and hence the cost is easily calculated explicitly as
(21)totalcost=(2(c1+c3)(μ+2λ)+3μc2)(3μ+2λ)-1.Figure 7
Each type of cost has a different impact on the value of the optimal buffer capacity.


#### 4.3.2. Bursty Part Arrivals

Next, we analyse the cost and the profit of kitting systems with different Markovian arrivals. We use the same values to model the arrivals as in Figure4. In Figures 8(a) and 8(b) the total cost (left) and the profit (right) versus the buffer capacity are depicted. We consider a holding cost c1 equal to 2, a shortage cost c2 equal to 55, and a loss cost c3 equal to 40. On the left figure, the optimal buffer capacity for Poisson arrivals is 12, for an interrupted Poisson arrival at one buffer it equals 22, and for an interrupted Poisson process at both buffers the optimal capacity is 28. As expected, the optimal buffer capacity increases as the burstiness in the production process increases. Concerning the profit analysis, we assume that the yield equals the throughput multiplied by a sale unit equal to 100. Assuming a maximum storage room of 30 parts, we observe that the optimal buffer capacity is 5, 11, and 12 for the depicted arrival processes, respectively. Consequently, kitting systems under production inefficiency require much higher storage space, especially when profit is applied as the parameter to determine the optimal buffer capacity. Moreover, the optimal capacity is very sensitive to the burstiness parameters σ and κ. In the example at hand, the plots suggest that the cost function is a convex function of the buffer capacity. However, one can also easily choose the cost parameters to obtain nonconvex cost functions.Production inefficiency results in higher required storage space.
(a)
(b)


#### 4.3.3. Phase-Type Kit Assembly Times

Finally, Figures9(a) and 9(b) depict the cost (left) and the profit (right) for a kitting system with Erlang-distributed kit assembly times versus the buffer capacity. As in the preceding figure, the profit equals the difference between the yield (equal to the throughput multiplied by 100) and the cost (defined by the parameters c1=2, c2=55 and c3=40). Out of the numerical examples, we observe that the shape of the service time distribution has no significant impact on profit and cost. These results confirm those found in Section 4.2.Shape of the service time distribution does not impact significantly the value of the optimal buffer capacity.
(a)
(b)


### 4.4. Performance Analysis of Solution Methods

We compare the performance of (sparse) GMRES and solving the Markov chain by standard LU decomposition Strang [19] on a kitting process with Poisson arrivals with rate λi=0.8 for part i={1,…,2} and with exponentially distributed kit assembly times with rate μ. Figure 10 depicts both methods in terms of speed versus the state space for a kitting process starting from a symmetric buffer capacity C1=C2=1 to C1=C2=60. While LU performs better than GMRES when the capacity (and state space) is small, the figure clearly shows that GMRES performs considerably is better than LU decomposition for capacity.Figure 10
The GMRES method performs well in terms of speed.


## 4.1. Bursty Part Arrivals

As a first example, we quantify the impact of production inefficiency on the performance of a kitting process. To this end, we compare part buffers with Poisson arrivals to corresponding kitting systems with interrupted Poisson arrivals. The arrival interruptions account for inefficiency in the production process. Kit assembly times are assumed to be exponentially distributed with service rate equal to one, this value being independent of the number of parts in the different buffers. This is a kitting process with Markovian arrivals as described in Example2 in Section 2.The interrupted Poisson process considered here is a two-state Markovian process. In the active state, new parts arrive in accordance with a Poisson process with rateλ whereas no new parts arrive in the inactive state. Let α and β denote the rate from the active to the inactive state and vice versa, respectively. We then use the following parameters to characterise the interrupted Poisson process:
(19)σ=βα+β,κ=1α+1β,ρ=λσ.
Note that σ is the fraction of time that the interrupted Poisson process is active, the absolute time parameter κ is the average duration of an active and an inactive period, and ρ is the arrival load of the parts.Figure3 shows the mean number of parts in buffer 1 versus the arrival load, for various values of the buffer capacities C1 and C2 for Poisson arrivals (for both buffers) as well as for interrupted Poisson arrivals (again for both buffers). The arrival load is set to ρ=0.8 for all curves. In addition, we set σ=0.4 and κ=10 for the interrupted Poisson processes. Clearly, the mean buffer content increases as the arrival load increases as expected. Moreover, if more buffer capacity is available, it will also be used: the mean buffer content increases for increasing values of C1=C2. Comparing interrupted Poisson and Poisson processes, burstiness in the production process has a negative impact on performance—more buffering is required—if the queues are not fully loaded (ρ<1). As for ordinary queues, the opposite can be observed for overloaded buffers. In this case, the effect of a large burst is mainly reflected in loss and not in additional queue content. Burstiness also yields larger periods without arrivals during which the buffer size decreases.Figure 3
When the queues are not fully loaded, burstiness in the production process has a negative impact on the performance.By numerical examples, we can quantify expected buffer behaviour—for example, more production yields higher buffer content, higher buffer capacity mitigates the loss probability, and so forth. However, less trivial behaviour can be observed as well. Figure4 depicts the probability that the buffer is full versus the buffer capacity C1=C2. We compare performance of kitting with Poisson arrivals to kitting with interrupted Poisson arrivals at one buffer and at both buffers. As in the preceding figure, the interrupted Poisson processes are characterised by σ=0.4 and κ=10. As expected, the probability that the buffer is full decreases for increasing values of the buffer capacities. Moreover, to reduce this probability, more buffer capacity is required for the case of two interrupted Poisson processes than for the case of two Poisson processes. For the kitting process with one Poisson and one interrupted Poisson process, non-trivial performance results can be observed. Namely, interruptions in the production of a part more negatively affect buffer performance of the other part. Indeed, buffer 1 is full with higher probability if the arrivals at buffer 2 are interrupted than if the arrivals at buffer 1 are interrupted. First note that the loss probabilities in both buffers are the same. For the Poisson buffer, parts are rejected at the arrival rate ρ1 whenever the buffer is full. For the IPP buffer, parts are rejected when the buffer is full and the arrival process is in the on-state at rate ρ2/σ. As the loss probability in both buffers is the same, we have Pr[Q1=C1]ρ1=Pr[Q2=C2,S=1]ρ2/σ, or equivalently Pr[Q1=C1]=Pr[Q2=C2∣S=1]. As the second queue is more likely to be full when there are arrivals we have Pr[Q2=C2∣S=1]≥Pr[Q2=C2] which shows Pr[Q1=C1]≥Pr[Q2=C2].Figure 4
Irregularity in the production of a part leads to a higher probability to have a full buffer for the other part.


## 4.2. Phase-Type Distributed Kit Assembly Times

The second numerical example quantifies the impact of the distribution of the kit assembly times on kitting performance. In particular, we here study Erlang-distributed kit assembly times. Limiting ourselves to Poisson arrivals to both buffers, this numerical example fits Example3 of Section 2.Figures5 and 6 depict the mean number of parts in buffer 1 and the loss probability in buffer 1 for the kitting process and, as a reference point, for the M/E/1/N queue as well. In both figures, the arrival load is varied and different values of the variance of the kitting time distribution are assumed as indicated. The mean kitting time is equal to 1 for all curves and the capacity of both buffers is equal to 20. In underload (ρ<1), kitting performs worse than the M/E/1/N queue: the mean buffer content and the loss probability have a higher value. This follows from the fact that kitting stops when one of the buffers is empty. By increasing the load, it is obvious that the buffer content converges to the capacity and the loss probability to one. It is most interesting to observe that the shape of the service time distribution only has a small effect on these performance measures. Indeed, there is no significant performance difference when σ2 equals 1/4 and when it equals 1/8.Figure 5
Variation of the kitting time distribution does not have a significant impact on the mean number of parts.Figure 6
Variation of the kitting time distribution does not have a significant impact on the loss probability.


## 4.3. Cost and Profit Analysis

We now add a cost structure to the kitting process under study. In particular, cost and profit for the kitting systems of Sections4.1 and 4.2 are analysed.The proposed cost function is(20)c1(EQ1+EQ2)+c2K+c3(b1+b2),
where c1 is the holding cost of a part in the buffer, c2 is the shortage cost in one or in both of the buffers, and c3 is the loss cost. Note that for all figures, the input parameters are symmetric for both parts such that EQ1=EQ2 and b1=b2.


### 4.3.1. Poisson Arrivals

In Figure7, the total cost of applying kitting systems versus the buffer capacity (varying from 1 to 30) is depicted. Limiting ourselves to Poisson arrivals for both buffers, ρi=λi=0.8 for part i={1,2}. We compare kitting performance for different costs as depicted. As expected, higher cost values lead to higher total cost. The cost structure also affects the optimal buffer capacity. When looking at the different costs separately, we observe that a higher holding cost decreases the optimal buffer capacity (from 15 to 12) and a higher loss cost increases the optimal capacity (from 12 to 14), whereas the optimal buffer capacity remains the same (12) for a higher shortage cost. Obviously, buffering is interesting when the storage cost is low and when the cost of rejecting parts (because of a full buffer) is high. It is most interesting to observe that, as the capacity increases, the cost models converge when the sum of the holding and the shortage cost (c1+c2) is equal. Indeed, as the loss tends to zero, the loss cost tends to zero as well. Furthermore, when the capacity equals one, the state space of the kitting model has size 4 and hence the cost is easily calculated explicitly as
(21)totalcost=(2(c1+c3)(μ+2λ)+3μc2)(3μ+2λ)-1.Figure 7
Each type of cost has a different impact on the value of the optimal buffer capacity.


### 4.3.2. Bursty Part Arrivals

Next, we analyse the cost and the profit of kitting systems with different Markovian arrivals. We use the same values to model the arrivals as in Figure4. In Figures 8(a) and 8(b) the total cost (left) and the profit (right) versus the buffer capacity are depicted. We consider a holding cost c1 equal to 2, a shortage cost c2 equal to 55, and a loss cost c3 equal to 40. On the left figure, the optimal buffer capacity for Poisson arrivals is 12, for an interrupted Poisson arrival at one buffer it equals 22, and for an interrupted Poisson process at both buffers the optimal capacity is 28. As expected, the optimal buffer capacity increases as the burstiness in the production process increases. Concerning the profit analysis, we assume that the yield equals the throughput multiplied by a sale unit equal to 100. Assuming a maximum storage room of 30 parts, we observe that the optimal buffer capacity is 5, 11, and 12 for the depicted arrival processes, respectively. Consequently, kitting systems under production inefficiency require much higher storage space, especially when profit is applied as the parameter to determine the optimal buffer capacity. Moreover, the optimal capacity is very sensitive to the burstiness parameters σ and κ. In the example at hand, the plots suggest that the cost function is a convex function of the buffer capacity. However, one can also easily choose the cost parameters to obtain nonconvex cost functions.Production inefficiency results in higher required storage space.
(a)
(b)


### 4.3.3. Phase-Type Kit Assembly Times

Finally, Figures9(a) and 9(b) depict the cost (left) and the profit (right) for a kitting system with Erlang-distributed kit assembly times versus the buffer capacity. As in the preceding figure, the profit equals the difference between the yield (equal to the throughput multiplied by 100) and the cost (defined by the parameters c1=2, c2=55 and c3=40). Out of the numerical examples, we observe that the shape of the service time distribution has no significant impact on profit and cost. These results confirm those found in Section 4.2.Shape of the service time distribution does not impact significantly the value of the optimal buffer capacity.
(a)
(b)


## 4.3.1. Poisson Arrivals

In Figure7, the total cost of applying kitting systems versus the buffer capacity (varying from 1 to 30) is depicted. Limiting ourselves to Poisson arrivals for both buffers, ρi=λi=0.8 for part i={1,2}. We compare kitting performance for different costs as depicted. As expected, higher cost values lead to higher total cost. The cost structure also affects the optimal buffer capacity. When looking at the different costs separately, we observe that a higher holding cost decreases the optimal buffer capacity (from 15 to 12) and a higher loss cost increases the optimal capacity (from 12 to 14), whereas the optimal buffer capacity remains the same (12) for a higher shortage cost. Obviously, buffering is interesting when the storage cost is low and when the cost of rejecting parts (because of a full buffer) is high. It is most interesting to observe that, as the capacity increases, the cost models converge when the sum of the holding and the shortage cost (c1+c2) is equal. Indeed, as the loss tends to zero, the loss cost tends to zero as well. Furthermore, when the capacity equals one, the state space of the kitting model has size 4 and hence the cost is easily calculated explicitly as
(21)totalcost=(2(c1+c3)(μ+2λ)+3μc2)(3μ+2λ)-1.Figure 7
Each type of cost has a different impact on the value of the optimal buffer capacity.


## 4.3.2. Bursty Part Arrivals

Next, we analyse the cost and the profit of kitting systems with different Markovian arrivals. We use the same values to model the arrivals as in Figure4. In Figures 8(a) and 8(b) the total cost (left) and the profit (right) versus the buffer capacity are depicted. We consider a holding cost c1 equal to 2, a shortage cost c2 equal to 55, and a loss cost c3 equal to 40. On the left figure, the optimal buffer capacity for Poisson arrivals is 12, for an interrupted Poisson arrival at one buffer it equals 22, and for an interrupted Poisson process at both buffers the optimal capacity is 28. As expected, the optimal buffer capacity increases as the burstiness in the production process increases. Concerning the profit analysis, we assume that the yield equals the throughput multiplied by a sale unit equal to 100. Assuming a maximum storage room of 30 parts, we observe that the optimal buffer capacity is 5, 11, and 12 for the depicted arrival processes, respectively. Consequently, kitting systems under production inefficiency require much higher storage space, especially when profit is applied as the parameter to determine the optimal buffer capacity. Moreover, the optimal capacity is very sensitive to the burstiness parameters σ and κ. In the example at hand, the plots suggest that the cost function is a convex function of the buffer capacity. However, one can also easily choose the cost parameters to obtain nonconvex cost functions.Production inefficiency results in higher required storage space.
(a)
(b)


## 4.3.3. Phase-Type Kit Assembly Times

Finally, Figures9(a) and 9(b) depict the cost (left) and the profit (right) for a kitting system with Erlang-distributed kit assembly times versus the buffer capacity. As in the preceding figure, the profit equals the difference between the yield (equal to the throughput multiplied by 100) and the cost (defined by the parameters c1=2, c2=55 and c3=40). Out of the numerical examples, we observe that the shape of the service time distribution has no significant impact on profit and cost. These results confirm those found in Section 4.2.Shape of the service time distribution does not impact significantly the value of the optimal buffer capacity.
(a)
(b)


## 4.4. Performance Analysis of Solution Methods

We compare the performance of (sparse) GMRES and solving the Markov chain by standard LU decomposition Strang [19] on a kitting process with Poisson arrivals with rate λi=0.8 for part i={1,…,2} and with exponentially distributed kit assembly times with rate μ. Figure 10 depicts both methods in terms of speed versus the state space for a kitting process starting from a symmetric buffer capacity C1=C2=1 to C1=C2=60. While LU performs better than GMRES when the capacity (and state space) is small, the figure clearly shows that GMRES performs considerably is better than LU decomposition for capacity.Figure 10
The GMRES method performs well in terms of speed.


## 5. Conclusion

In this paper, we evaluate the performance of two-part kitting processes in a Markovian setting. Furthermore, a cost-profit analysis is conducted and an approximation for our model is given. Note that the particularity of the studied kitting systems is that the part buffers are paired. This means that each demand requires both parts such that a kit can only be assembled if both inventories are nonempty. Methodologically, we have applied sparse matrix techniques (e.g., GMRES) as most of the entries in the generator matrix have a value equal to zero. The solution is not exact but performs well in terms of solution speed and accuracy.As our numerical results show, the interplay between the different queues leads to complex performance behaviour. For example, interruptions in the production of a part more negatively affect buffer performance of the other part. Overall, we observe extreme sensitivity of kitting performance on arrival process parameters while performance is reasonable insensitive to variation of the kitting time distribution. Furthermore, we determine the optimal buffer capacity based on a cost-profit analysis and establish a quite good approximation for our two-part kitting process.



---

*Source: 843184-2013-11-21.xml*