**GENERAL ARTICLES**

**Approximations to performance measures in queuing systems**

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

^{I}Department of Industrial Engineering, Eastern Mediterranean University, Famagusta, Turkey

^{II}Department of Industrial Engineering, Eastern Mediterranean University, Famagusta, Turkey

^{III}Department 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 *G*I*/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 _{2}),* 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 _{2}* 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*and provides exact three-moment match wherever possible. The approximated

_{1}, m_{2}, m_{3}),*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}

**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 _{1}(*0), m_{2},* and

*m*of the distribution function

_{3}(3m_{2}^{2}=2m_{1}m_{3})*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 _{0}* are constants determined as follows. Assuming the existence of moments of the density function

*f(t),*we can express

*Φ*

*(s)*as

where *m _{0}=1* and

*m*is the

_{n}*n*order moment about the origin of f. Using (9) in (7) and (8), and comparing the coefficient of

^{th}*s*and

^{0}, s^{1}*s*on both sides, we obtain (using some algebra) the values of

^{2}*A, B,*and

*s*as given in (5). This completes the proof.

_{0}** Note1:** For the approximation to hold, it is necessary that

*s*Simple calculations show that this condition implies that

_{0}<0.*Φ*

*and (see Figure 1), where*

_{2}>2*C*is the squared coefficient of variation of inter-arrival time,

^{2}*Φ*

*and*

_{2}= C^{2}+1,

**Figure 1: Feasible regions for the approximation**

]]>

**It is worth mentioning that the error in the approximation (4) is of**

*Note2:**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

_{2}), 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_{1}, *m _{2},* m

_{3}}. 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*are the minimum and maximum values of

_{u}*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 _{2})* and Coxian

*(K*distributions. On this point, we use the former for analysis, as these distributions provide a fairly accurate match when the

_{2})*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

^{2}*H*distributions for inter-arrival distributions seems to provide useful results.

_{2}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 ^{2}>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*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

^{2}>1.*GI/G/1*queue becomes considerable as

*C*increases. If the service time distribution has

^{2}*C*the impact of the third moment is not significant [3, 15]. Since the present work deals with matching an

^{2}<1,*H*distribution with three moments, we confine our attention to the range

_{2}*C*only.

^{2}>1It is well known [15] that three numbers *m _{1}, m_{2},* 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

_{3}*F,*then an

*H*distribution exists with these three moments if, and only if, the first three moments of

_{2}*F*satisfy the conditions. However, Karlin & Studden [7] have shown that

*m*and

_{1}, m_{2},*m*are the moments of some probability distributions on the positive real line if, and only if,

_{3}*m*

_{1}> 0,*Φ*

*and*

_{2}> 1,*Φ*

*. Thus, in the region where exact three-moment matching is not possible, researchers have used adhoc methods to find the approximate*

_{3}>*H*distribution. These regions are clearly shown in Figure 2.

_{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

_{1},

*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*

_{2}-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

_{1},

*m2, m3).*

__Algorithm__

Given the first three moments of a distribution function *F* (say, *m _{1}, m_{2},* and

*m*which are assumed to be finite), the parameters

_{3},*p (0<p<1),*

*µ*

*and*

_{1},*µ*

*of the*

_{2}*H*distribution given by (11) - whose first three moments either match

_{2}*m*and

_{1}, m_{2¡}*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:

_{3}** 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 _{2}* distribution):

*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*

^{2}__>__1Also, *µ** _{1}>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 inand 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 _{2}* 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

_{2}*H*distributions. Also note that these two methods improve in accuracy for increasing

_{2}*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 _{1}* 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

_{2}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 _{2})* 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

_{1}*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

_{2} 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_{1}=15, c_{2}=10,**

*Θ=1,*

*ρ=0.65,*and Coxian distribution parameters are*p=0.8,*

*λ*_{1}=2,

*λ*Optimal service rate_{2}=0.2)

*η**during servers vacation period with Coxian inter-arrival distribution

]]>

**Figure 4:**

*(c*_{w}=4, c_{1}=15, c_{2}=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_{2}) 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

]]>