Performance Analysis of a Kitting Process as a Paired Queue

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.

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. Section 2 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 Figure 1. Each queue has a finite capacity—let 𝐶_ℓ 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 𝒞₁×𝒞₂×𝒦, whereby 𝒞_ℓ={0,…,𝐶_ℓ} for ℓ=1,2 and with 𝒦={1,2,…,𝐾} being the state space of the modulating process. At any time, the state of the kitting process is described by the triplet [𝑚,𝑛,𝑖], with 𝑚 and 𝑛 being the number of parts in the first and second buffer, respectively, and 𝑖 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 𝛼_𝑖𝑗 denote the transition rate from state 𝑖 to state 𝑗(𝑖,𝑗∈𝒦,𝑖≠𝑗) and let 𝐀 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 𝜆^(ℓ)_𝑖𝑗 denote the (marked) transition rate from state 𝑖 to state 𝑗 when there is an arrival at buffer ℓ, ℓ=1,2. Moreover, let 𝐋_ℓ denote the corresponding generator matrix. Note that such marked transitions from state 𝑖 to state 𝑖 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 𝜇_𝑖𝑗 and 𝐌 denote the corresponding transition rate and generator matrix, respectively. Summarising, arrivals at and departures from the buffers are described by the generator matrices 𝐀, 𝐋₁, 𝐋₂, and 𝐌. So far, no diagonal elements of 𝐀 have been defined. To simplify notation, it will be further assumed that the diagonal elements are chosen such that the row sums of 𝐀+𝐋₁+𝐋₂+𝐌 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, 𝐌^(𝑚,𝑛) denotes the generator matrix of state transitions with departure when there are 𝑚 parts in buffer 1 and 𝑛 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, 𝜆^(ℓ)_𝑖𝑗(𝑚,𝑛) denotes the arrival rate at buffer ℓ=1,2 from state 𝑖 to state 𝑗 when there are 𝑚 parts in buffer 1 and 𝑛 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 𝜆₁ and 𝜆₂ 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, (𝑚,𝑛). We have 𝐌=𝜇,𝐋₁=𝜆₁,𝐋₂=𝜆₂,𝐀=−𝜆₁−𝜆₂−𝜇. (1)

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 Λ₁ of transitions with arrivals at buffer 1, the generator matrix Λ₂ with arrivals at buffer 2, and the generator matrix Λ₀ without arrivals. As usual, the diagonal elements of Λ₀ are negative and ensure that the row sums of Λ₀+Λ₁+Λ₂ are zero. Retaining exponentially distributed kit assembly times, we have 𝐋₁=Λ₁,𝐋₂=Λ₂,𝐀=Λ₀−𝜇𝐈,𝐌=𝜇𝐈. (2) Here 𝐈 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 𝐓 which corresponds to non-absorbing transitions [14]. Let 𝐭^′=−𝐓𝐞^′ be the column vector with the rates to the absorbing state and let 𝐟 be a row vector with zero elements except the first one. Assuming Poisson arrivals in both buffers (with rates 𝜆₁ and 𝜆₂, resp.), we get the following matrices: 𝐋₁^(𝑚,𝑛)=𝜆₁𝐈1−𝟙_{{𝑚=0,𝑛>0}}+𝜆₁𝐞^′𝝉𝟙_{{𝑚=0,𝑛>0}},𝐋₂^(𝑚,𝑛)=𝜆₂𝐈1−𝟙_{{𝑚>0,𝑛=0}}+𝜆₂𝐞^′𝝉𝟙_{{𝑚>0,𝑛=0}},𝐀^(𝑚,𝑛)=𝐓𝟙_{{𝑚>0,𝑛>0}}−𝜆₁𝐈−𝜆₂𝐌𝐈,_(𝑚,𝑛)=𝐭^′𝝉𝟙_{{𝑚>1,𝑛>1}}+𝐭^′𝐟1−𝟙_{{𝑚>1,𝑛>1}}. (3) 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 [𝑄₁(𝑡),𝑄₂(𝑡),𝑆(𝑡)], 𝑄_ℓ(𝑡) being the number of parts in buffer ℓ at time 𝑡 and 𝑆(𝑡) being the state of the background process at time 𝑡.

Let 𝜋_𝑖(𝑚,𝑛)=lim_𝑡→∞Pr[𝑄₁(𝑡)=𝑚,𝑄₂(𝑡)=𝑛,𝑆(𝑡)=𝑖] be the steady-state probability to be in state [𝑚,𝑛,𝑖] and let 𝝅(𝑚,𝑛) be the vector with elements 𝜋_𝑖(𝑚,𝑛), for 𝑖∈𝒦. Figure 2 shows a fragment of the transition rate diagram of the kitting model in state [𝑚,𝑛,𝑖]. As mentioned above, two independent input streams arrive at the buffers with intensity 𝜆^(ℓ)_𝑖𝑗 and are processed into kits with intensity 𝜇_𝑖𝑗. 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 𝑖. Moreover, possible dependence of the transition rates on the buffer sizes is not indicated.

Figure 2 

Fragment of the transition rate diagram for state (𝑚,𝑛,𝑖).

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<𝑛<𝐶₁ and 0<𝑚<𝐶₂. We have 𝜋_𝑖(𝑚,𝑛)^𝐾_𝑗=1𝜆⁽¹⁾_𝑖𝑗(𝑚,𝑛)+^𝐾_𝑗=1𝜆⁽²⁾_𝑖𝑗(𝑚,𝑛)+𝜇_𝑖𝑗(𝑚,𝑛)+^𝐾_{𝑗=1,𝑗≠𝑖}𝛼_𝑖𝑗(=𝑚,𝑛)_𝐾𝑗=1𝜋_𝑗(𝑚−1,𝑛)𝜆⁽¹⁾_𝑗𝑖()+𝑚−1,𝑛_𝐾𝑗=1𝜋_𝑗(𝑚,𝑛−1)𝜆⁽²⁾_𝑗𝑖()+𝑚,𝑛−1_𝐾𝑗=1𝜋_𝑗(𝑚+1,𝑛+1)𝜇_𝑗𝑖()+𝑚+1,𝑛+1_{𝐾𝑗=1,𝑗≠𝑖}𝜋_𝑗(𝐦,𝐧)𝛼_𝑖𝑗(𝑚,𝑛), (4) for 𝑖,𝑗∈𝒦 or equivalently, 𝝅(𝑚−1,𝑛)𝐋₁^{(𝑚−1,𝑛)}+𝝅(𝑚,𝑛−1)𝐋₂^{(𝑚,𝑛−1)}+𝝅(𝑚+1,𝑛+1)𝐌^{(𝑚+1,𝑛+1)}+𝝅(𝑚,𝑛)𝐀^(𝑚,𝑛)=𝟎. (5)(ii) If buffer 1 is empty and buffer 2 is neither empty nor full (𝑚=0 and 0<𝑛<𝐶₂), we have 𝜋_𝑖(0,𝑛)^𝐾_𝑗=1𝜆⁽¹⁾_𝑖𝑗(0,𝑛)+^𝐾_𝑗=1𝜆⁽²⁾_𝑖𝑗()+0,𝑛_{𝐾𝑗=1,𝑗≠𝑖}𝛼_𝑖𝑗(=0,𝑛)_𝐾𝑗=1𝜋_𝑗(0,𝑛−1)𝜆⁽²⁾_𝑗𝑖()+0,𝑛−1_𝐾𝑗=1𝜋_𝑗(1,𝑛+1)𝜇_𝑗𝑖()+1,𝑛+1_{𝐾𝑗=1,𝑗≠𝑖}𝜋_𝑗(0,𝑛)𝛼_𝑗𝑖(0,𝑛), (6)?for 𝑖,𝑗∈𝒦 or equivalently, 𝝅(0,𝑛−1)𝐋₂^(0,𝑛−1)+𝝅(1,𝑛+1)𝐌^(1,𝑛+1)+𝝅(0,𝑛)𝐀^(𝟎,𝐧)+diag𝐌^(0,𝑛)𝐞^′=𝟎. (7)(iii) Similarly, if buffer 2 is empty and buffer 1 is neither empty nor full (𝑛=0 and 0<𝑚<𝐶₁), we have 𝝅(𝑚−1,0)𝐋₁^(𝑚−1,0)+𝝅(𝑚+1,1)𝐌^(𝑚+1,1)+𝝅(𝑚,0)𝐀^(𝐦,𝟎)+diag𝐌^(𝑚,0)𝐞^′=𝟎. (8)(iv) If both buffers are empty (𝑚=0 and 𝑛=0), we have 𝝅(1,1)𝐌^(1,1)+𝝅(0,0)𝐀^(0,0)+diag𝐌^(0,0)𝐞^′=𝟎. (9)(v) If buffer 1 is empty and buffer 2 is full (𝑚=0 and 𝑛=𝐶₂), we get 𝝅0,𝐶₂−1𝐋^(0,𝐶2₂⁻¹⁾+𝝅0,𝐶₂𝐀^(0,𝐶2)+diag𝐌^(0,𝐶2)𝐞^′+𝐋^(0,𝐶2)₂𝐞^′=𝟎. (10)(vi) Similarly, if buffer 1 is full and buffer 2 is empty (𝑚=𝐶₁ and 𝑛=0), we have 𝝅𝐶₁−1,0𝐋^(𝐶1₁^−1,0)+𝝅𝐶₁,0𝐀^(𝐶1,0)+diag𝐌^(𝐶1,0)𝐞^′+𝐋^(𝐶1₁^,0)𝐞^′=𝟎. (11)(vii) Finally, if both buffers are full (𝑚=𝐶₁ and 𝑛=𝐶₂), we find 𝝅𝐶₁−1,𝐶₂𝐋^{(𝐶1−1,𝐶2)}₁+𝝅𝐶₁,𝐶₂−1𝐋^(𝐶1,𝐶2₂⁻¹⁾+𝝅𝐶₁,𝐶₂𝐀^{(𝐶1,𝐶2)}+diag𝐋^{(𝐶1,𝐶2)}₁𝐞^′+𝐋^{(𝐶1,𝐶2)}₂𝐞^′=𝟎. (12)

3.2. Performance Measures

Given the steady-state vectors 𝝅(𝑚,𝑛), we can now obtain a number of interesting performance measures for the kitting system. For ease of notation, let 𝜋(𝑚,𝑛)=𝝅(𝑚,𝑛)𝐞^′ denote the probability to have 𝑚 parts in buffer 1 and 𝑛 parts in buffer 2. Moreover, let 𝝅⁽¹⁾(𝑚)=∑_𝑛𝝅(𝑚,𝑛) and 𝝅⁽²⁾(𝑛)=∑_𝑚𝝅(𝑚,𝑛) 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 𝜋⁽¹⁾(𝑚)=𝝅⁽¹⁾(𝑚)𝐞^′ and 𝜋⁽²⁾(𝑛)=𝝅⁽²⁾(𝑛)𝐞^′.

The following performance measures are of interest. (i) The mean buffer 1 and 2 content is 𝐸𝑄₁ and 𝐸𝑄₂, respectively, 𝐸𝑄₁=^𝐶1_𝑚𝜋⁽¹⁾(𝑚)𝑚,𝐸𝑄₂=^𝐶2_𝑛𝜋⁽²⁾(𝑛)𝑛. (13)(ii) The variance of the buffer 1 and 2 is Var𝑄₁ and Var𝑄₂, respectively, Var𝑄₁=^𝐶1_𝑚𝜋⁽¹⁾(𝑚)𝑚²−𝐸𝑄₁²,Var𝑄₂=^𝐶2_𝑛𝜋⁽²⁾(𝑛)𝑛²−𝐸𝑄₂². (14)(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 𝜌_eff=1−𝜋⁽¹⁾(0)−𝜋⁽²⁾(0)+𝜋(0,0). (15)(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 𝜂=^𝐶1𝐶_𝑚=1²_𝑛=1𝝅(𝑚,𝑛)𝐌^(𝑚,𝑛)𝐞^′. (16)(v) The shortage probability 𝐾 is the probability that one of the buffers is empty as 𝐾=𝜋⁽¹⁾(0)+𝜋⁽²⁾(0)−𝜋(0,0). (17)(vi) The loss probability 𝑏_𝑖 in buffer 𝑖 is the probability that an arriving part cannot be stored in buffer 𝑖, 𝑖=1,2. By noting that the accepted arrival load equals the departure rate, we find 𝑏_𝑖=𝜌_𝑖−𝜂𝜌_𝑖, (18) where 𝜌_𝑖 is the arrival load of part 𝑖={1,2}. If the arrival load in both buffers is the same, then the loss probability is also the same: 𝑏₁=𝑏₂. 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 𝐴𝑥=𝑏, the GMRES algorithm approximates 𝑥 by the vector 𝑥_𝑛∈𝐾_𝑛 in a Krylov subspace 𝐾_𝑛=span{𝑏,𝐴𝑏,…,𝐴^𝑛−1𝑏} which minimises the norm of the residual ‖𝐴𝑥_𝑛−𝑏‖. 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 𝑛 iterations grows by 𝑂(𝑛²) which is a major drawback of GMRES. This limitation is usually overcome by restarting the algorithm. After a chosen number of iterations 𝑚, the accumulated data is cleared and the intermediate results are used as the initial data for the next 𝑚 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 Example 2 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: 𝛽𝜎=1𝛼+𝛽,𝜅=𝛼+1𝛽,𝜌=𝜆𝜎. (19) 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.

Figure 3 shows the mean number of parts in buffer 1 versus the arrival load, for various values of the buffer capacities 𝐶₁ and 𝐶₂ 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 𝐶₁=𝐶₂. 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. Figure 4 depicts the probability that the buffer is full versus the buffer capacity 𝐶₁=𝐶₂. 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 𝜌₁ 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 𝜌₂/𝜎. As the loss probability in both buffers is the same, we have Pr[𝑄₁=𝐶₁]𝜌₁=Pr[𝑄₂=𝐶₂,𝑆=1]𝜌₂/𝜎, or equivalently Pr[𝑄₁=𝐶₁]=Pr[𝑄₂=𝐶₂∣𝑆=1]. As the second queue is more likely to be full when there are arrivals we have Pr[𝑄₂=𝐶₂∣𝑆=1]≥Pr[𝑄₂=𝐶₂] which shows Pr[𝑄₁=𝐶₁]≥Pr[𝑄₂=𝐶₂].

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 Example 3 of Section 2.

Figures 5 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 𝜎² 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 Sections 4.1 and 4.2 are analysed.

The proposed cost function is 𝑐₁𝐸𝑄₁+𝐸𝑄₂+𝑐₂𝐾+𝑐₃𝑏₁+𝑏₂, (20) where 𝑐₁ is the holding cost of a part in the buffer, 𝑐₂ is the shortage cost in one or in both of the buffers, and 𝑐₃ is the loss cost. Note that for all figures, the input parameters are symmetric for both parts such that 𝐸𝑄₁=𝐸𝑄₂ and 𝑏₁=𝑏₂.

4.3.1. Poisson Arrivals

In Figure 7, 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, 𝜌_𝑖=𝜆_𝑖=0.8 for part 𝑖={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 (𝑐₁+𝑐₂) 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 totalcost=2𝑐₁+𝑐₃𝜇+2𝜆+3𝜇𝑐₂3𝜇+2𝜆⁻¹. (21)

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 Figure 4. 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 𝑐₁ equal to 2, a shortage cost 𝑐₂ equal to 55, and a loss cost 𝑐₃ 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.

(a)

(b)

Figure 8 

Production inefficiency results in higher required storage space.

4.3.3. Phase-Type Kit Assembly Times

Finally, Figures 9(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 𝑐₁=2, 𝑐₂=55 and 𝑐₃=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.

(a)

(b)

Figure 9 

Shape of the service time distribution does not impact significantly the value of the optimal buffer capacity.

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 𝜆_𝑖=0.8 for part 𝑖={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 𝐶₁=𝐶₂=1 to 𝐶₁=𝐶₂=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.