- About this Journal ·
- Abstracting and Indexing ·
- Advance Access ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Mathematical Problems in Engineering
Volume 2013 (2013), Article ID 843184, 10 pages
Performance Analysis of a Kitting Process as a Paired Queue
Department of Telecommunications and Information Processing, Ghent University, St. Pietersnieuwstraat 41, 9000 Gent, Belgium
Received 19 February 2013; Revised 14 May 2013; Accepted 14 May 2013
Academic Editor: Zhengguang Wu
Copyright © 2013 E. De Cuypere et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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.
Nowadays manufacturing systems are often composed of multiple in-house fabrication units . 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 . 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 . 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  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.  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  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  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 .
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 , —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.
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 for and with 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 , . 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 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 from state to state when there are parts in buffer and parts in buffer .
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
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 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 . 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: 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 .
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 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.
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, and . We have for or equivalently, (ii) If buffer 1 is empty and buffer 2 is neither empty nor full ( and ), we have ?for or equivalently, (iii) Similarly, if buffer 2 is empty and buffer 1 is neither empty nor full ( and ), we have (iv) If both buffers are empty ( and ), we have (v) If buffer 1 is empty and buffer 2 is full ( and ), we get (vi) Similarly, if buffer 1 is full and buffer 2 is empty ( and ), we have (vii) Finally, if both buffers are full ( and ), we find
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 and 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 and content is and , respectively, (ii) The variance of the buffer and is and , respectively, (iii) The effective load of the system is the amount of time that kitting is ongoing. As kitting is only ongoing when none of the buffers is empty, we have (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 (v) The shortage probability is the probability that one of the buffers is empty as (vi) The loss probability in buffer is the probability that an arriving part cannot be stored in buffer , . By noting that the accepted arrival load equals the departure rate, we find where is the arrival load of part . 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  (p. 194–197); Philippe et al. , Buchholz , and Saad and Schultz . To solve the linear equation , the GMRES algorithm approximates by the vector in a Krylov subspace 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 .
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: 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 for all curves. In addition, we set and 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 (). 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.
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 and . 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 , or equivalently . As the second queue is more likely to be full when there are arrivals we have which shows .
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 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 . In underload (), 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 and when it equals .
4.3. Cost and Profit Analysis
The proposed cost function is 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 to ) is depicted. Limiting ourselves to Poisson arrivals for both buffers, for part . 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 to ) and a higher loss cost increases the optimal capacity (from to ), whereas the optimal buffer capacity remains the same () 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 and hence the cost is easily calculated explicitly as
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 , a shortage cost equal to , and a loss cost equal to . On the left figure, the optimal buffer capacity for Poisson arrivals is , for an interrupted Poisson arrival at one buffer it equals , and for an interrupted Poisson process at both buffers the optimal capacity is . 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 . Assuming a maximum storage room of parts, we observe that the optimal buffer capacity is , , and 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.
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 , and ). 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.
4.4. Performance Analysis of Solution Methods
We compare the performance of (sparse) GMRES and solving the Markov chain by standard LU decomposition Strang  on a kitting process with Poisson arrivals with rate for part 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 to . 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.
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.
- L. Medbo, “Assembly work execution and materials kit functionality in parallel flow assembly systems,” International Journal of Industrial Ergonomics, vol. 31, no. 4, pp. 263–281, 2003.
- Y. A. Bozer and L. F. McGinnis, “Kitting versus line stocking: a conceptual framework and a descriptive model,” International Journal of Production Economics, vol. 28, no. 1, pp. 1–19, 1992.
- P. Som, W. Wilhelm, and R. Disney, “Kitting process in a stochastic assembly system,” Queueing Systems, vol. 17, no. 3-4, pp. 471–490, 1994.
- H. Bryzner and M. I. Johansson, “Design and performance of kitting and order picking systems,” International Journal of Production Economics, vol. 41, no. 1–3, pp. 115–125, 1995.
- S. Ramachandran and D. Delen, “Performance analysis of a Kitting process in stochastic assembly systems,” Computers & Operations Research, vol. 32, no. 3, pp. 449–463, 2005.
- R. Ramakrishnan and A. Krishnamurthy, “Analytical approximations for kitting systems with multiple inputs,” Asia-Pacific Journal of Operational Research, vol. 25, no. 2, pp. 187–216, 2008.
- B. Johansson and M. I. Johansson, “High automated Kitting system for small parts: a case study from the Volvo uddevalla plant,” in Proceedings of the 23rd International Symposium on Automotive Technology and Automation, pp. 75–82, Vienna, Austria, 1990.
- W. J. Hopp and J. T. Simon, “Bounds and heuristics for assembly-like queues,” Queueing Systems, vol. 4, no. 2, pp. 137–155, 1989.
- J. M. Harrison, “Assembly-like queues,” Journal of Applied Probability, vol. 10, pp. 354–367, 1973.
- G. Latouche, “Queues with paired customers,” Journal of Applied Probability, vol. 18, no. 3, pp. 684–696, 1981.
- E. De Cuypere and D. Fiems, “Performance evaluation of a kitting process,” in Proceedings of the 17th International Conference on analytical and stochastic modelling techniques and applications, vol. 6751, Lecture Notes in Computer Science, Venice, Italy, 2011.
- P. Buchholz, P. Kemper, and J. Kriege, “Multi-class Markovian arrival processes and their parameter fitting,” Performance Evaluation, vol. 67, no. 11, pp. 1092–1106, 2010.
- D. Fiems, B. Steyaert, and H. Bruneel, A Genetic Approach To Markovian Characterisation of H. 264 Scalable Videopages, Multimedia Tools and Applications, 2011.
- G. Latouche and V. Ramaswami, Introduction To Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM Series on Statistics and Applied Probability, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999.
- W. J. Stewart, Introduction To the Numerical Solution of Markov Chains, Princeton University Press, Princeton, NJ, USA, 1994.
- B. Philippe, Y. Saad, and W. Stewart, “Numerical methods in markov chain modeling,” Operations Research, vol. 40, no. 6, pp. 1156–1179, 1992.
- P. Buchholz, “Structured analysis approaches for large Markov chains,” Applied Numerical Mathematics, vol. 31, no. 4, pp. 375–404, 1999.
- Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986.
- G. Strang, Linear Algebra and Its Applications, Academic Press, Orlando, FL,, 2nd edition, 1980.