GENERAL ARTICLES

 

Approximations to performance measures in queuing systems

 

 

N.S. Kambo^I; A. Rangan^II; E. Moghimihadji^{III, *}

^IDepartment of Industrial Engineering, Eastern Mediterranean University, Famagusta, Turkey
^IIDepartment of Industrial Engineering, Eastern Mediterranean University, Famagusta, Turkey
^IIIDepartment of Industrial Engineering, Istanbul Aydin University, Istanbul, Turkey. ehsanhadji@aydin.edu.tr

 

 

---

ABSTRACT

Approximations to various performance measures in queuing systems have received considerable attention because these measures have wide applicability. In this paper we propose two methods to approximate the queuing characteristics of a GI/M/1 system. The first method is non-parametric in nature, using only the first three moments of the arrival distribution. The second method treads the known path of approximating the arrival distribution - by a mixture of two exponential distributions - by matching the first three moments. Numerical examples and optimal analysis of performance measures of GI/M/1 queues are provided to illustrate the efficacy of the methods, and are compared with benchmark approximations.

---

OPSOMMING

Benaderings tot verskeie prestasiemaatstawwe in toustaanstelsels ontvang beduidende aandag weens die wye toepasbaarheid daarvan. In hierdie artikel word twee metodes voorgestel om die toustaaneienskappe van die GI/M/1-stelsel te benader. Die eerste metode is nie-parametries van aard en gebruik slegs die eerste drie momente van die aankomsverdeling. Die tweede metode volg die bekende roete om die aankomsverdeling te benader deur die eerste drie momente te pas, deur middel van 'n kombinasie van twee eksponensiële verdelings. Die doelteffendheid van die metodes word aan die hand van numeriese voorbeelde en optimale ontleding van die prestasiemaatstawwe van GI/M/1-stelsels bewys, en word vergelyk met benaderings.

---

 

 

1. INTRODUCTION

The growth of queuing theoretic applications has been phenomenal, ranging from communication and multimedia systems to inventory and reliability theory. This has led to a sustained interest in methods to evaluate performance measures in queuing theory. In the case of non-Markovian queues, the computations of these measures involve the arrival and/or service distributions explicitly. However, in practical applications like management, optical, and communication networks, the specific forms of these distributions might not be known. At best, one might only be in possession of the moments of the underlying distribution. There are a number of cases where the moments of a distribution are easily obtained, but theoretical distributions are not available in closed forms [9]. Alternatively, from the sample data observed, efficient estimators for the various moments of the underlying distributions can be calculated. Thus, computation of performance measures based on the first two or three moments of the arrival and/or service distributions is very useful. Whitt [16], in a classic exposition, discussed approximations using extremal distributions giving the upper and lower bounds for the performance measures in a GI/M/1 system. Smith [13] proposed a two-moment approximation for the probability distribution of M/G/1/K systems, and extended it to the analysis of M/G/1/K queuing networks. Sohn & Lee [14] conducted a Monte Carlo simulation in order to study the relation between various performance measures in a G/G/1 queue. Recent work on such systems with working vacations for the server have immense applications in ATM machines and internet systems, such as optical nets, electric nets, and communication nets [8, 4, 2]. In these applications, the arrival epochs could be observed or, at worst, simulated. Our motivation in this paper has thus been to obtain approximations to the performance measures of a GI/M/1 system using only the first three moments of the arrival distribution, without explicit recourse to the arrival distribution

We will discuss our problem with specific reference to a GI/M/1 queuing system, even though our methods also work in a similar way for other non-Markovian queues. Consider a GI/M/1 queue whose traffic intensity is ρ= E (service time)/E (arrival time), L is the expected equilibrium queue length, and σ is the steady state probability that a customer will have to wait to begin his service. It is well known that

where σ is the unique root in the open interval (0, 1) of the equation

with µ=1/E(service time) and Φ(s), the Laplace-Stieltjes Transform of the inter-arrival distribution function, say F, given by:

We note that the evaluation of the performance measures σ and L require prior knowledge of the inter-arrival distribution function F, and not just the moments of F. As mentioned in the beginning of this section, many queuing applications are likely to produce only the moments of F and not the distribution itself. Thus the problem is to find σ and L on the basis of the first few moments of F only. Whitt [16] showed that there is a considerable reduction in the range of possible values of σ and L when the third moment is also used, compared with using just two moments of F. We propose simple and accurate methods to evaluate σ and L based on the first three moments of the distribution function F in the absence of any knowledge of the form of F. In section 2, we propose a non-parametric method based only on the first three moments of F - without recourse to approximate F - by another distribution function. Numerical illustrations are provided to compare the values of σ and L, using the present method, with their exact values. The method provides exact- results for certain important arrival distributions like Erlang of order 2, Coxian (K₂), a mixture of two exponentials, and exponential distribution. We also provide two optimisation illustrations to obtain economic performance measures in the application of GI/M/1 queuing systems.

Approximations of probability distributions by phase type distributions (by matching moments up to a certain order) have attracted the attention of researchers because of necessity, and because of their wide applicability. Pioneering work on phase type distributions and their various applications was done by Neuts [11]. Among the various members of the family of phase type distributions that have been studied, mixtures of two exponential distributions (known as H₂ distributions) play a key role in many approximations used in queuing theory. These distributions are log convex in nature, and can approximate highly-skewed distributions quite accurately. In section 3, we suggest a simple nonlinear programming method in which the first two moments are matched exactly, while the third moment is matched as closely as possible. This method works in the entire region of possible values (of m₁, m₂, m₃), and provides exact three-moment match wherever possible. The approximated H₂ distribution with the given three moments of the inter-arrival distribution is then used to calculate σ and L. Numerical illustrations are provided to validate the approximation and computation of the performance measures.

 

2. A NON-PARAMETRIC METHOD

We observe from equation (2) that the computation of the performance measures σ and L in the GI/M/1 system requires the use of Φ(s), the Laplace Transform of the density function f corresponding to the distribution function F. However, without prior knowledge of F, and armed only with the first three moments of F, an approximation to 0(s) is obtained from the proposition that follows.

Proposition

Suppose that the first three raw moments m₁(*0), m₂, and m₃ (3m₂²=2m₁m₃) of the distribution function F exist and are known. Then the following approximation to the Laplace Transform of the distribution function F holds.

Proof

In the classical renewal theory, the renewal density m(t) of a renewal process with interval density f(t) satisfies the integral equation:

Applying Laplace Transform to both sides of (6), we obtain the Laplace Transform of m(t) as:

where Φ (s) is the Laplace Transform of the density function f.

Now Φ (0) =1 and is non zero. This means that the denominator 1-Φ (s) of (7) has a simple zero at s=0. Thus we may approximate m (s) by a rational function of the form

where A, B, and s₀ are constants determined as follows. Assuming the existence of moments of the density function f(t), we can express Φ (s) as

where m₀=1 and m_n is the n^th order moment about the origin of f. Using (9) in (7) and (8), and comparing the coefficient of s⁰, s¹ and s² on both sides, we obtain (using some algebra) the values of A, B, and s₀ as given in (5). This completes the proof.

Note1: For the approximation to hold, it is necessary that s₀<0. Simple calculations show that this condition implies that Φ₂>2 and (see Figure 1), where C² is the squared coefficient of variation of inter-arrival time, Φ₂= C²+1, and

 

Figure 1: Feasible regions for the approximation

 

Note2: It is worth mentioning that the error in the approximation (4) is of o(1) as s -0.

Note3: The condition that s0 is non-positive, which is necessary for the approximation to hold, is satisfied by many standard arrival distributions: uniform, gamma, mixed exponential, lognormal, Coxian (K₂), mixture of Erlangian, and Weibull. Also s0 is non-positive for truncated normal and inverse Gaussian probability density functions under certain conditions.

Using (4) in equations (2) and (1) with the values of m1, m2, and m3, the performance measures σ and L were obtained immediately. To illustrate the efficiency of the proposed method, we present (in Table 1) the values of σ and L computed using (4) for certain choices of the set {m₁, m₂, m₃}. In order to compare the approximations with exact values, we have considered the values of the moments of commonly used arrival distributions -namely, gamma (Table 1) and PH4 distributions (Table 2). The values of σ and L are calculated for various values of traffic intensity p. Eckberg [5] specified upper and lower bound distributions that yield the maximum and minimum mean queue length in a steady state among all inter-arrival time distributions with two and three moments specified. Using these distributions, Whitt [16] calculated the maximum relative error for σ and L using the formula

 

Table 1: Approximations for σ and L with gamma inter-arrival distribution

 

where L_t and L_u are the minimum and maximum values of L using the lower and upper bound distributions. Table 1 also presents the upper and lower bounds for σ and L given two and three moments, and the corresponding maximum relative errors specified in (10). It can be seen that our method captures the values of σ and L with low relative errors.

 

3. INTER-ARRIVAL DISTRIBUTION APPROXIMATION BY MATCHING MOMENTS

Approximating general distributions by phase type distributions is important in queuing theory because their structure leads to Markovian state description and consequently analytical tractability. Although several phase type distributions have been used in the literature, two distributions that are prime candidates for such approximations (because of their simplicity and suitability) are mixtures of two exponentials (H₂) and Coxian (K₂) distributions. On this point, we use the former for analysis, as these distributions provide a fairly accurate match when the C² is large - which is true for arrival distribution in a queuing system. Furthermore, in some of the examples discussed by Whitt [17], the maximum relative error (MRE) when two moments are fitted was found to be 200 percent, while working with mixtures of exponentials reduced the MRE to 50 percent. Specifying the third moment reduced the MRE to 5 percent. Thus, a three-moment match using H₂ distributions for inter-arrival distributions seems to provide useful results.

Using an empirical study, Bere [3] showed that when the service time distribution is approximated using the first two of its moments, the third moment has a considerable effect on the probability distribution of the number of customers in an M/G/1 queue if C²>1 . He also showed that the probability distribution of the number of customers and the average number of customers in Á(n)/G/1/N system are highly sensitive to the third moment of the service time distribution if C²>1. Altiok [1], in justifying the inclusion of the third moment in matching, refers to Bere's empirical work. Whitt [15] empirically showed that the effect of the third moment on the average number in the system in a GI/G/1 queue becomes considerable as C² increases. If the service time distribution has C²<1, the impact of the third moment is not significant [3, 15]. Since the present work deals with matching an H₂ distribution with three moments, we confine our attention to the range C²>1 only.

It is well known [15] that three numbers m₁, m₂, and m₃ can be the first three raw moments of a distribution function F, provided that. Furthermore, if the first three moments exist for a distribution F, then an H₂ distribution exists with these three moments if, and only if, the first three moments of F satisfy the conditions. However, Karlin & Studden [7] have shown that m₁, m₂, and m₃ are the moments of some probability distributions on the positive real line if, and only if, m₁> 0, Φ₂> 1, and Φ₃> . Thus, in the region where exact three-moment matching is not possible, researchers have used adhoc methods to find the approximate H₂ distribution. These regions are clearly shown in Figure 2.

 

Figure 2: Region of three-moment matching (Region-I: Exact three-moment match possible; Region-II: Exact three-moment match not possible; Region-III: Infeasible region for three moments of a distribution function)

 

The parameters p, µ1, and µ2 of H2 (see (11) below) when (m₁, m2, m3) falls in Region-I, and where the three moments can be matched exactly, are well documented. However, when these moments fall in Region-II, such that the three moments cannot be matched exactly, methods suggested in the literature are recipes in nature. Lopez-Herrero [10], in the absence of information on service distribution, used the maximum entropy principle approach to estimate the true distribution of the number of customers served during the busy period in an M/G/1 retrial system. Whitt [15] suggests that "if m3 turns out to be too small when attempting an H₂-fit, one procedure is to replace m3 by something slightly larger than" . Altiok [1] suggests the use of but does not indicate how to calculate the perturbation parameter ε. In the following algorithm, we propose a goal programming procedure of matching the first two moments exactly, and matching the third moment as closely as possible in Region-II. This procedure subsumes Region-I as a particular case. We note that Region-III is an infeasible region for the existence of (m₁, m2, m3).

Algorithm

Given the first three moments of a distribution function F (say, m₁, m₂, and m₃, which are assumed to be finite), the parameters p (0<p<1), µ₁, and µ₂ of the H₂ distribution given by (11) - whose first three moments either match m₁, m_2¡ and m₃ exactly, or match the first two moments exactly, and match the third moment as closely as possible in the sense of squared differences - are given by the following algorithm:

Step 1: Find optimal p using the Golden section method [12, 6] to solve the following one dimension optimisation problem:

Step 2: Using the value of p found in step 1, compute:

The steps of the algorithm are justified using the following arguments. Consider the following probability density function of a mixture of two exponentials (H₂ distribution):

In the following, we do not consider the trivial cases of p=0, 1 as they lead us to exponential density.

The first three moments of the density function above are:

From (12) we have

Substituting (15) in (13), and after simple algebra

we have two cases. First we consider

We know that

, which implies that C²>1

Also, µ₁>0 implies that

Substituting (17) in (15) results in

The third parameter p is obtained by matching the third moment as closely as possible. Thus, we minimise

subject to:

In the second case, we set

Using algebra similar to the first case results in

and the third parameter p is determined by solving

It is easily seen that both cases lead to the same result, but with the roles of p and 1-p interchanged.

In Table 2 we continue with PH4 inter-arrival distribution. To illustrate this, we fit H₂ distribution for this distribution by matching the moments specified by the algorithm. The performance measures σ and L obtained by using the fitted H₂ distribution in steps 1 and 2 are given in Table 2. Attention is drawn to the MRE for these measures when three moments are matched [16] and to the actual relative error in using the approximated H₂ distributions. Also note that these two methods improve in accuracy for increasing p as expected.

 

Table 2: Approximations for o and L with PH4 inter-arrival distribution

 

4. TWO OPTIMISATION ILLUSTRATIONS

4.1 Illustration 1

In practice, queuing managers are generally interested in optimising the model parameters under their control by minimising the operating cost or maximising the business profit. In the first illustration, we will be interested in obtaining the optimal service rate in a cost minimisation problem for a GI/M/1 queuing system. The objective cost function consists of two components: the cost due to customers waiting in line (known as the delay cost), and the service cost rate. Thus the cost function to be minimised is given by:

where λ and µ are the arrival and service rates respectively, W is the expected waiting time of a customer in the system, c₁ is the expected cost per unit time of a customer's wait, and c₂ is the service cost rate. Using Little's formula, (24) reduces to

The optimal µ* of the objective function above was computed using our moments matching method (introduced in section 3) by assuming the first three moments of the arrival distribution only, and the cost rates. However, in order to compare our results with the exact values, the moments were chosen to correspond to Coxian (K₂) and Inverse Gaussian distributions commonly used in queuing theory. The results are presented in Tables 3 and 4. When the Coxian arrival distribution was used, our method provided the exact values of µ*. In the case of Inverse Gaussian distribution, the relative errors were significantly small.

 

Table 3. The optimal service rate µ* with Coxian inter-arrival distribution. (The optimal µ* computed, using our method and using the Coxian distribution, match exactly)

 

 

Table 4. The optimal service rate µ* with Inverse Gaussian inter-arrival distribution

 

4.2 Illustration 2

In the second illustration, we consider an optimisation problem in a GI/M/1 queue with working vacation for the server. These problems have wide application in Internet systems such as optical, electrical, and communication nets [8]. We consider a single server queuing system that has the general arrival process. The working vacation and vacation interruption are connected, and the server enters into vacation when there are no customers, and it can take service at the lower rate during the vacation period. If there are customers in the system at the instant of a service completion during the vacation period, the server will- return to the normal working level regardless of whether the vacation ends. Otherwise, it continues the vacation. The performance measure L, the mean queue length, and P(J=0) and P(J=1) - which are the state probabilities of a server in the steady state - have been derived by Li et al. [8]. We refer the reader to their paper for the relevant expressions. Li et al. [8] considered the problem of optimising the service rate η during the server's vacation period for a given cost structure. Let c_w represent the unit time cost of every waiting customer, and c₁ and c₂ the service costs per unit time during the normal working level and vacation period respectively. The expected net cost function to be optimised can be seen to be

where µ is the service rate during the service period. The optimal service rate η* was computed using our non-parametric method of section 2 for certain values of the model parameters and cost parameters. We have used the Coxian arrival distribution and its moments (as used in Illustration 1) to obtain the optimal service rate η*. However, the Inverse Gaussian distribution could not be used, as the objective function (26) loses its convexity and becomes monotonic. We have used Erlangian of order 2 (used by Li et al. [8]) in its place. In Figures 3 and 4, we present the values of η versus the associated cost. The optimal η* and the corresponding cost obtained using our method are very close to the exact values.

 

Figure 3: (c_w=4, c₁=15, c₂=10, Θ=1, ρ=0.65, and Coxian distribution parameters are p=0.8, λ₁=2, λ₂=0.2) Optimal service rate η* during servers vacation period with Coxian inter-arrival distribution

 

 

Figure 4: (c_w=4, c₁=15, c₂=10, Θ=1, ρ=0.65, and Erlangian distribution parameters are K=2, M=2.5) Optimal service rate η* during servers vacation period with Erlangian of order 2 inter-arrival distribution

 

5. CONCLUSION

This study introduces two methods for evaluating performance measures in a GI/M/1 queuing system in the absence of information on the arrival distribution, and when only the first three moments are known. The first method is non-parametric as it does not use the distribution function, whereas the second method uses an H2 distribution obtained by moment matching procedure. This procedure involves the computationally economical Golden section method. Note that a Coxian (K₂) distribution is also a good phase type distribution to consider. It is worth pursuing the regions (Φ2, Φ3) in which each of these approximations scores higher than the others in terms of relative errors. The usefulness of the methods in optimisation procedures has been illustrated with examples.

 

REFERENCES

[1] Altiok, T. 1985. On the phase-type approximations of general distributions. IIE Transactions, 17(2), pp 110-116.         [ Links ]

[2] Baba, Y. 2005. Analysis of a GI/M/1 queue with multiple working vacations. Operations Research Letters, 33(2), pp 201-209.         [ Links ]

[3] Bere, B. 1981. Influence du moment d'ordre 3 sur les files d'attente a lois de services generales. Thesis (MS). Universite de Rennes, France.         [ Links ]

[4] Chae, K.C., Lim, D.E. & Yang, W.S. 2009. The GI/M/1 queue and the GI/Geo/1 queue both with single working vacation. Performance Evaluation, 66(7), pp 356-367.         [ Links ]

[5] Eckberg, A.E. 1977. Sharp bounds on Laplace-Stieltjes transforms, with applications to various queuing problems. Mathematics of Operations Research, 2(2), pp 135-142.         [ Links ]

[6] Kambo, N.S. 1991. Mathematical programming techniques, revised edition, East-West Press Private Limited.         [ Links ]

[7] Karlin, S. & Studden, W.J. 1966. Tchebycheff systems: With applications in analysis and statistics. New York: John Wiley & Sons.         [ Links ]

[8] Li, J., Tian, N. & Ma, Z. 2008. Performance analysis of GI/M/1 queue with working vacations and vacation interruption. Applied Mathematical Modeling, 32(12), pp 2715-2730.         [ Links ]

[9] Lindsay, B.G., Pilla, R.S. & Basak, P. 2000. Moment-based approximations of distributions using mixtures: Theory and applications. Ann. Inst. Statist. Math, 52(2), pp 215-230.         [ Links ]

[10] Lopez-Herrero, M.J. 2002. On the number of customers served in the M/G/1 retrial queue: First moments and maximum entropy approach. Journal of Computers and Operations Research, 29(12), pp 1739-1757.         [ Links ]

[11] Neuts, M.F. 1981. Matrix geometric solutions in stochastic models: An algorithmic approach. John Hopkins University Press.         [ Links ]

[12] Rao, S.S. 2009. Engineering optimization: Theory and practice. 4^th edition, John Wiley & Sons.         [ Links ]

[13] Smith, J.M. 2011. Properties and performance modeling of finite buffer M/G/1/K networks. Journal of Computers and Operations Research, 38(4), pp 740-754.         [ Links ]

[14] Sohn, S.Y. & Lee, S.H. 2004. Sensitivity analysis for output performance measures in long-range dependent queuing system. Journal of Computers and Operations Research, 31(9), pp 1527-1536.         [ Links ]

[15] Whitt, W. 1982. Approximate a point process by a renewal process, I: Two basic methods. Operations Research, 30(1), pp 125-147.         [ Links ]

[16] Whitt, W. 1984. On approximations for queues, I: Extremal distributions. AT & T Bell Laboratory Technical Journal, 63(1), pp 115-138.         [ Links ]

[17] Whitt, W. 1984. On approximations for queues, III: Mixtures of exponential distributions. AT & T Bell Laboratory Technical Journal, 63(1), pp 163-175.         [ Links ]

 

 

* Corresponding author