An algorithm for calculating steady state probabilities of queueing systems
Abstract
This paper presents a method for calculating steady state probabilities of queueing systems. The infinitesimal generator matrix is used to define all possible states in the system and their transition probabilities. While this matrix can be written down immediately for many other queueing systems with phasetype service times (e.g. Coxian, Hypoexponential, …), it requires a more careful analysis for systems with Erlangian service times. The constructed matrix may then be used to calculate steady state probabilities using an iterative algorithm. The resulting steady state probabilities can be used to calculate various performance measures, e.g. the average queue length. Additionally, computational issues of the implementation are discussed and an example from the field of telecommunication callcenter queue length will be outlined to substantiate the applicability of these efforts. In the appendix, tables of the average queueing length given a specific number of service channels,
traffic density, and system size are presented.
Keywords: Phasetype distributions, Erlang queueing systems, steady state probabilities, generator matrix–performance measures
1 Introduction
Multiserver queueing systems with Poisson input and phasetype distributed service times ( queueing systems) are an important extension of simple queueing systems
The class of Erlang() distributed service times () are important, because for constant mean service time the family of Erlangian distributions interpolates infinitely many distributions between the negative exponential and constant service time . Many solution approaches for different types of Erlangian distributed queueing systems have been proposed in the literature. Shapiro [16] studied the system and proved that the probability of jobs in the system can be expressed as a linear combination of powers. Poyntz and Jackson [14] analyzed the and the system by applying generating functions. Tables of such results have been published by Sakasegawa [15]. The system was analyzed by Mayhugh and McCormick [11] and Heffer [5]. They determined the generating function of the stationary state probabilities. However, their results are from a computational point of view only useful for small values of and only. Yu [18] solved the queueing system with heterogeneous servers by generalizing the approach of Mayhugh and McCormick [11]. Hillier and Lo [6] presented some numerical results based on the above procedure for the special case of homogeneous servers. Approximation formulas for the average waiting time and average system size in equilibrium were given by Page [12] and Smith [17]. Hillier and Lo also presented tables and graphs of performance measures for different queueing systems. However, steady state probabilities for these queueing systems with limited waiting room are not given yet.
This paper presents an easy to implement algorithm to calculate the steady state probabilities of queueing systems. For this, we use a similar approach like Mayhugh and McCormick [11]. Firstly, all possible states in the queueing system are determined and a generator matrix is built up. This generator matrix determines  via a well known theorem  the transition matrix. This transition matrix is used to find the steady state probabilities for being in a particular stage in the queueing system. These can be summed up to get the steady state probabilities for a specific number of customers in the queueing system. These probabilities can then be used to calculate different performance measures like the average queueing length and, with the help of Little’s theorem [9], the average waiting time in the queue. In contrast to many other queueing systems with phasetype distributed service times, where the generator matrix can be written down immediately, the calculation of the matrix for systems needs a more detailed description, and is the main topic of this paper.
This paper is organized as follows. Section 2 provides an overview of queueing systems and summarizes the algorithms used to calculate the generator matrix and derive steady state probabilities. Section 3 gives an example how this method can be used in telecommunication callcenters to calculate performance measures and discusses computational issues of the algorithm presented above. Additionally, Appendix A contains tables for the average queueing length with different traffic densities for the system.
2 The queueing system
2.1 Assumptions
The queueing system is defined by the following assumptions:

The arrival process of customers follows a Poisson distribution with intensity rate .

The population of customers is infinite.

The service system consists of service channels.

Only one customer can be served in one service channel.

The maximum queue size is .

If an arriving customer finds all service channels busy and the maximum queue size is not reached, she joins the single waiting line which is served in the order of arrival.

An arriving customer is forced to balk if she arrives at a time when the queue size is at its limit .

If an arriving customer finds more than one service channel vacant, she randomly selects a free service channel.

The customer at the head of the waiting line is forwarded to the first vacant channel without delay, i.e. the queue discipline is FCFS.

The service process of customers follows an Erlang probability distribution of degree . The distribution of service time is equal in every channel.
2.2 Building the generator matrix
The probability density of an Erlang() distributed service time is
(1) 
Each service channel may be regarded as consisting of ordered stages such that the conditional probability of transition of a customer from any stage to the succeeding stage in the time interval is . Equivalently, the distribution (1) can be seen as the distribution of the sum of independent negative exponentially distributed random variables each with the same parameter . The reason for using Erlang distributed random variables in Queuing theory lies in the fact, that the phase method can be used to describe the process as a function of a Markov process. The cost is an enlarged state space for the model, and consequently an increased complexity of the numerical solution. The benefit, of course, is in enlarging the class of service time distributions for which the model is solvable.
If the service time is Erlang() distributed (1), the mean service time is and the mean traffic density per service channel in an system is thus
(2) 
can be regarded as the mean traffic density per stage. It is well known in queueing theory [1] that a steady state solution exists if and only if . This means e.g. for queueing systems with limited waiting room and that the average queueing length is near or at its maximum.
The transition matrix of a continuoustime Markov chain can be expressed [13] as
(3) 
where is the intensity matrix .
For an ergodic process, the limit of the transition matrix
(4) 
exists, is of rank 1 and its identical rows coincide with the stationary distribution of the Markov process. Instead of solving the linear equations , , numerical solutions also can be obtained by using the generator matrix, formula 3 and taking the power of until no changes in the rows occur anymore.
We use this method to calculate the stationary probability distributions of the queueing system with the assumptions made above. First, we define the different states with the vector where is the queueing length and is the total number of customers being in phase , this is a similar approach as in Mayhugh and McCormick [11].
The total number of possible states can be calculated through the equation:
(5) 
or in closed form, calculated with the Zeilberger algorithm
where the first part in (5) is the number of states where no customers are waiting and the second part is the number of states where customers are waiting to get served . The third part refers to the empty state for . A possible algorithm to calculate these states as matrix is given in Algorithm 1.
Algorithm 1. Initial Matrix Setup 
\li \li \li\While \li \If \kwthen \li \li \li \kwend if \li \li \For \To \li \If \kwthen \li \li \li \kwend if \li \kwnext \li\kwend while
It should be clear that other algorithms, like the lexicographically ordering in Mayhugh and McCormick [11], could also be used. However, the states created here are the same, but in different orders. We note, that we only have to find all possible states for a specific queuing system and it is not necessary to order the states in a specific way to build up our generating matrix.
After that, we have to create all possible states where customers are waiting. The number of states in which customers are waiting is independent from and therefore a constant equal to . Furthermore, the set of possible states for different numbers of waiting customers are equal except in . E.g. the set of states with are identical to the set of states with except for the number of customers in the queue. We can use the created states from the algorithm above, to find all possible states for the queuing system with waiting customers. Therefore, we are looking only at the created states where the , because in this situation arriving customers could not be served and have to wait in the queue. To create all possible states with we only have to take the states created by the above algorithm with and and set instead of , the same procedure can be used for until .
Because of (5) the generator matrix must have size . Below the possible transitions from one state to another are summarized:

If the queue is empty () and at least one service channel is free () an arriving customer starts phase 1 (from () to () with intensity ).

If the queue is not at its maximum () and all service channels are busy () the arriving customer joins the queue (from () to () with intensity )..

A customer being in phase transits to phase with intensity .

If a customer is in phase she leaves the system with intensity . If there is a customer waiting in the queue she starts in phase , otherwise the server is idle until the next customer arrives.
The rows of the matrix Q can be seen as the starting point of each state and the columns are filled in the way described above. After that, the () can be calculated by the formula:
(6) 
Matrix Setup Example
As an example we consider the queue. First we give a table of all possible states for this system, it was generated with the algorithm given above:
0  0  0  0  0  0  1  1  1  

0  1  2  0  1  0  2  1  0  
0  0  0  1  1  2  0  1  2 
The generator matrix build with the algorithm above is
The first row represents the state (0,0,0), as we can see the only transition to another state is given by (0,1,0) with intensity . The sixth row represents the state (0,0,2), this state can go to (0,0,1) with intensity and to (1,0,2) with intensity .
2.3 Calculating steady state probabilities
The algorithm for calculating steady state probabilities of the queueing system needs five steps and an optional sixth step:

Calculate all possible states

Generate the matrix with the steps explained above

Use formula (5) to get the transition matrix

Take the power of the transition matrix iteratively until the columns don’t change anymore (see Algorithm 2)

Take the sum of the different steady state probabilities for the stages to get the steady state probabilities for customers in the system

(Optional) Calculate performance measures
Step 4 is needed to calculate the steady state probabilities for the different stages. Let be the transition matrix, a possible algorithm could be a simple iteration like shown in algorithm 2 where is some critical level defined by the user.
Algorithm 2. Iterative Transition Matrix Calculation 
\li \li \li\While( \kwdo \li \li \li\kwend while
Step 5 in the algorithm above can be calculated in various ways. The chosen method depends on the algorithm which is used to create the possible states. If an algorithm is used which creates states like in Mayhugh and McCormick [11] formulas to calculate the probability of customers in the system can be devised easily. However, to calculate the possible states with the algorithm presented here, a different formula has to be used. Denote the state probabilities with . We first notice that , and if then can be calculated with the formula given in Mayhugh and McCormick [11], pp. 710. For with we have to sum up the states for which .
Step 6 is optional, if some performance measures are requested, they can be easily calculated with the steady state probabilities , see for example Gross and Harris [4]. E.g. the average system size can be calculated, using the well known formula:
2.4 Empirical computational issues
The current algorithm is easy to implement, but a large waiting queue with more than phases leads to high computational times. The following runtime experiments have been conducted on a Pentium IV computer with 2.6 Ghz and 1 GB RAM using Microsoft Windows XP Professional and MatLab 6.5. Figure 1 shows the calculation time in seconds for average queueing length of a system with , for and .
3 Telecommunication CallCenters
A call center is a service network in which agents provide telephonebased services. The typical call center setup is shown in figure 2. When all agents are busy customers seeking these services are delayed in queues, hence it is convenient to model call centers as a queueing system. Processwise these queues can be compared to inventories in manufacturing (justintime, timebasedcompetition, …). But human queues include personal preferences, complaints, abandonments and the like. Thus, customers are likely to base judgments about the serviceproviding company on their queueingexperience. Therefore the goal of a company providing such a service is to minimize the average waiting time of their costumers. This can always be accomplished by extending the number of agents, which in turn raises the cost of the call center significantly. The decision problem is schematically shown in figure 3. To take a good decision it is important to calculate the average queue length as correct as possible.
A survey of queueing system theory for call centers was published by Koole and Mandelbaum [10] and more detailed by Koole et al. [3]. Recently Ishay [7] applied phase type distributions to fit call center data and found out that phase type distributions of order can be used to fit the service durations of callcenter data for different priorities and servicetypes. The general structure of order already provides a reasonable fit to the overall service time. She used the program EMpht (see S. Asmussen et al. [2]) to fit phasetype distributions to available data from a call center of a large Israelian bank.
4 Conclusion
In this paper we presented an easytoimplement algorithm for calculating steady state probabilities of queueing systems. In contrast to other types of queueing systems with phasetype distributed service times, the main problem with the type of queueing system considered throughout this paper is how to generate the generator matrix. Hence, an algorithm to generate this matrix was presented. The resulting methodology can be used for several practical applications, including telecommunication callcenters. Phasetype distributions and associated average system size calculations are well suited to redesign and redimension existing call centers.
Future research includes an extension of this methodology to (especially ) queueing systems.
Appendix A Tables
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.400  1.198  1.958  2.606  2.876  3.110  3.214  3.274  3.293  
3  0.400  1.211  2.090  3.051  3.528  3.977  4.184  4.303  4.342  
K  5  0.400  1.212  2.132  3.384  4.190  5.071  5.511  5.768  5.852 
7  0.400  1.212  2.136  3.480  4.472  5.683  6.325  6.707  6.832  
10  0.400  1.212  2.136  3.527  4.664  6.212  7.090  7.622  7.797 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.600  1.798  2.952  3.941  4.346  4.689  4.840  4.924  4.951  
3  0.600  1.803  3.046  4.348  4.979  5.559  5.824  5.974  6.022  
K  5  0.600  1.803  3.076  4.646  5.612  6.648  7.159  7.455  7.552 
7  0.600  1.804  3.078  4.731  5.879  7.253  7.975  8.401  8.540  
10  0.600  1.804  3.079  4.773  6.060  7.775  8.739  9.320  9.511 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.800  2.399  3.957  5.309  5.857  6.315  6.512  6.621  6.656  
3  0.800  2.401  4.025  5.681  6.470  7.184  7.505  7.685  7.743  
K  5  0.800  2.401  4.045  5.950  7.078  8.268  8.848  9.181  9.290 
7  0.800  2.401  4.047  6.026  7.332  8.867  9.665  10.13  10.28  
10  0.800  2.401  4.048  6.062  7.503  9.382  10.42  11.05  11.26 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  1.000  2.999  4.965  6.696  7.395  7.971  8.216  8.350  8.393  
3  1.000  3.000  5.013  7.035  7.987  8.837  9.215  9.425  9.492  
K  5  1.000  3.000  5.028  7.279  8.571  9.915  10.56  10.93  11.05 
7  1.000  3.000  5.029  7.347  8.815  10.50  11.38  11.88  12.05  
10  1.000  3.000  5.029  7.380  8.977  11.01  12.14  12.81  13.03 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  1.500  4.499  7.482  10.20  11.30  12.19  12.56  12.76  12.82  
3  1.500  4.500  7.502  10.47  11.84  13.05  13.56  13.85  13.94  
K  5  1.500  4.500  7.508  10.66  12.38  14.11  14.92  15.38  15.53 
7  1.500  4.500  7.509  10.71  12.59  14.69  15.74  16.35  16.54  
10  1.500  4.500  7.509  10.74  12.74  15.19  16.50  17.28  17.53 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.400  1.198  1.959  2.610  2.881  3.117  3.223  3.282  3.301  
3  0.400  1.211  2.086  3.049  3.533  3.990  4.202  4.323  4.363  
K  5  0.400  1.211  2.121  3.353  4.163  5.068  5.527  5.795  5.883 
7  0.400  1.211  2.124  3.431  4.414  5.655  6.329  6.732  6.864  
10  0.400  1.211  2.124  3.466  4.576  6.151  7.074  7.640  7.827 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.600  1.798  2.953  3.945  4.352  4.697  4.848  4.933  4.960  
3  0.600  1.803  3.045  4.350  4.987  5.576  5.846  5.998  6.047  
K  5  0.600  1.803  3.070  4.626  5.596  6.655  7.184  7.491  7.591 
7  0.600  1.803  3.072  4.696  5.836  7.239  7.990  8.436  8.582  
10  0.600  1.803  3.073  4.727  5.989  7.730  8.737  9.350  9.552 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.800  2.399  3.958  5.313  5.863  6.322  6.521  6.630  6.665  
3  0.800  2.401  4.024  5.684  6.480  7.204  7.529  7.712  7.771  
K  5  0.800  2.401  4.042  5.936  7.068  8.281  8.878  9.222  9.334 
7  0.800  2.401  4.044  6.000  7.299  8.862  9.688  10.17  10.33  
10  0.800  2.401  4.044  6.028  7.446  9.349  10.43  11.09  11.31 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  1.000  2.999  4.965  6.700  7.400  7.978  8.224  8.359  8.402  
3  1.000  3.000  5.013  7.040  7.998  8.858  9.240  9.453  9.521  
K  5  1.000  3.000  5.026  7.269  8.567  9.933  10.59  10.97  11.10 
7  1.000  3.000  5.027  7.327  8.789  10.51  11.41  11.93  12.11  
10  1.000  3.000  5.027  7.353  8.929  10.99  12.16  12.86  13.09 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.400  1.198  1.959  2.612  2.885  3.121  3.227  3.286  3.306  
3  0.400  1.210  2.084  3.047  3.534  3.996  4.211  4.334  4.373  
K  5  0.400  1.211  2.116  3.336  4.147  5.065  5.534  5.809  5.899 
7  0.400  1.211  2.118  3.405  4.381  5.637  6.328  6.743  6.880  
10  0.400  1.211  2.118  3.434  4.525  6.114  7.062  7.647  7.841 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.600  1.798  2.953  3.947  4.355  4.701  4.853  4.937  4.965  
3  0.600  1.803  3.044  4.350  4.991  5.585  5.857  6.010  6.060  
K  5  0.600  1.803  3.068  4.614  5.585  6.657  7.195  7.509  7.611 
7  0.600  1.803  3.069  4.677  5.811  7.227  7.996  8.454  8.604  
10  0.600  1.803  3.069  4.703  5.949  7.701  8.732  9.364  9.573 
0.1  0.3  0.5  0.7  0.8  0.9  0.95  0.98  0.99  

1  0.800  2.399  3.958  5.315  5.866  6.326  6.525  6.635  6.670  
3  0.800  2.401  4.024  5.686  6.485  7.213  7.542  7.726  7.785  
K  5  0.800  2.401  4.041  5.928  7.061  8.286  8.893  9.243  9.357 
7  0.800  2.401  4.042  5.986  7.279  8.856  9.698  10.19  10.36  
10  0.800  2.401  4.042  6.009  7.413  9.327  10.43  11.11  11.33 
Footnotes
 using the notation suggested by Kendall [8]
References
 S. Asmussen. Applied probability and queues, volume 51 of Applications of Mathematics. Springer, 2nd edition, 2003.
 S. Asmussen, O. Nerman, and M. Olsson. Fitting phasetype distribution via the EM algorithm. Scandinavian Journal of Statistics, 23:419–441, 1996.
 N. Gans, G. Koole, and A. Mandelbaum. Telephone call centers: Tutorial, review, and research prospects. Manufacturing & Service Operations Management, 5:79–141, 2003.
 D. Gross and C. M. Harris. Applied probability and queues. Wiley, New York, 1974.
 J. C. Heffer. Steady state solution of the queueing system. INFOR J. Canadian O.R.S., 17:16–30, 1969.
 F.S. Hillier and F.D. Lo. Tables for multiserver queueing systems involving erlang distributions. Technical Report 14, Operations Research Department, Stanford University, 1971.
 E. Ishay. Fitting phasetype distributions to data from a telephone call center. Master’s thesis, Technion  Israel Institute Of Technology, 2003.
 D. Kendall. Some problems in the theory of queue. Journal of the Royal Statistical Society, 13:151–153, 1951.
 L. Kleinrock. Queueing systems, Volume 1: Theory. Wiley, New York, 1975.
 G. Koole and A. Mandelbaum. Queueing models of call centers: An introduction. Annals of Operations Research, 113:41–59, 2002.
 J. O. Mayhugh and R. E. McCormick. Steadystate solution of the . Management Science, 14:692–712, 1968.
 E. Page. Tables of waiting time for and and their use to give approxiate waiting times in more general queues. J. Op. Res. Soc., 33:453–473, 1982.
 G. Ch. Pflug. Stochastische Modelle in der Informatik. Stuttgart, 1986.
 C.D. Poyntz and R.R.P. Jackson. The steadystate solution for the queuing process . Operations Research Quarterly, 24:615–625, 1973.
 H. Sakasegawa. Numerical tables of the queueing systems 1:. Institute of Statistical Mathematics, Computer Science Monograph, 1978.
 S. Shapiro. The Mserver queue with poisson input and gamma distributed service of order two. Operations Research, 14:685–694, 1966.
 V. L. Smith. Approximating the distribution of customers in queues. J. Op. Res. Soc., 36:327–332, 1985.
 O.S. Yu. On the steadystate solution of an queue with heterogeneous servers. Technical Report 38, Operations Research Department, Stanford University, 1971.