# Simple and Yet Efficient Estimators for Markovian Multiserver Queues

Authors: Emilio Suyama; Roberto C. Quinino; Frederico R. B. Cruz
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Estimators for the parameters of the Markovian multiserver queues are presented, from samples that are the number of clients in the system at arbitrary points and their sojourn times. As estimation in queues is a recognizably difficult inferential problem, this study focuses on the estimators for the arrival rate, the service rate, and the ratio of these two rates, which is known as the traffic intensity. Simulations are performed to verify the quality of the estimations for sample sizes up to 400. This research also relates notable new insights, for example, that the maximum likelihood estimator for the traffic intensity is equivalent to its moment estimator. Some limitations of the results are presented along with a detailed numerical example and topics for future developments in this research area.


---

## Body

## 1. Introduction

One of the major research interests in queuing theory is to study the performance characteristics of the queues as a function of their parameters, namely, the arrival rateλ, the service rateμ, or else the traffic intensityρ. The mathematical problem dealt with here is how these parameters can be estimated using some statistical inferential method. This article focuses on theM/M/s queuing system, which in Kendal notation represents the Markovian number of arrivals, exponential (Markovian) service times,s identical servers working in parallel, and an infinite maximum capacity of users simultaneously allowed in the system, which is one of the classical basic queuing models [1]. From the above parameters and reminding that forM/M/s queues the traffic intensity is defined as the ratio ρ=λ/sμ, important performance measures can be derived, such as the empty system probabilityP0, the expected number of customers in the systemL, the expected number of customers in the queueLq, the expected time in the systemW, and the expected time in the queueWq.One of the classical basic queuing models [1],M/M/s queues have many applications and consequently great necessity of inference. In fact, many real-life systems are approximately described by queuing models, including computer and telecommunication networks [2–4], manufacturing and service systems [5–7], and heathcare systems [8–10].Classical inference results for queues (based on the maximum likelihood) are unknown before the results reported in parallel by Clark [11], with the development of the maximum likelihood estimates forλ andμ for Markovian single-served queues, that is,M/M/1 queues, in Kendal notation, and by Benes [12], for Markovian infinite queues, that is,M/M/∞ queues, in Kendal notation. More recently, Schruben and Kulkarni [13] proved that popular performance measures unfortunately do not have expected values and standard errors, Basawa and Prabhu [14] discovered some asymptotic properties of the estimators for the performance evaluation ofM/M/1 queues, Zheng and Seila [15] showed that under a certain upper limit for the traffic intensity (ρ<1), the estimators do exist and have the necessary properties, and Almeida et al. [16] proposed bias correction methods for classical estimators forM/M/1 queues.However, probably because of the inherent computational difficulties, Bayesian inference has not gained enough interest until the pioneering results reported by Armero and Bayarri [17–19], which were followed by the improvements from Armero and Conesa [20–23], McGrath et al. [24], McGrath and Singpurwalla [25], Choudhury and Borthakur[26], and, more recently, Almeida and Cruz [27], Cruz et al. [28], Quinino and Cruz [29], Cruz et al.[30], and Choudhury and Basak, [31], among many others.Almost all the literature that addresses estimation issues describes methods that require the continuous observation of the process over a fixed interval of time. Among the exceptions is the article by Basawa et al. [32], which considers the estimation forM/M/1 queues from the sojourn time data, and the articles by Choudhury and Borthakur [26] and Chowdhury and Mukherjee [33], in which the number of customers present in theM/M/1 queue at successive departure epochs is used to estimate the traffic intensity. In Ross et al. [34], a method that requires the numbers in the queue at successive time points, which are not necessarily equally spaced, has been developed. The specific results presented by Ross et al. [34] apply toM/M/s queues with large number of servers (in particular when s>40, a limit that generates the smallest estimation errors).Thus, this study aims to derive a maximum likelihood estimator forρ inM/M/s queues, from samples of the number of clients in the system at arbitrary points, which has been shown to be conveniently equivalent to the moment estimator. Additionally, with the use the additional information from samples of the sojourn times, a maximum likelihood estimator for the arrival and the service rates is derived. To the best of the authors’ knowledge, these have not been done before.Summarizing, the main contributions of this paper are as follows. The traffic intensity moment estimator is proved to be equal to the maximum likelihood estimator. Also, it is shown that if additional information of the sojourn times is used then estimators for the arrival and service rates can be deduced. Finally, it is demonstrated by extensive computational experiments that the method of moments is competitive to estimate arrival and service rates inM/M/s queues.The rest of this paper is organized as follows. In the next section, the probabilistic model of theM/M/s queue is detailed, the respective maximum likelihood estimator forρ is deduced, demonstrating its equivalence to the moment estimator, and the estimators forλ andμ are derived. In the following section, results from Monte Carlo simulations are presented to attest the efficiency of the estimators, along with an illustrative numerical example. In the last section, the paper is concluded, with discussions, final remarks, and the presentation of some topics for future research in this area.


## 2. Materials and Methods

### 2.1. Probabilistic Model

InM/M/s queues, the arrivals of the customers are assumed to follow a Poisson process with rateλ, which means that the times between successive arrivals are independent exponentially distributed random variables with mean 1/λ. Upon arrival, a customer goes immediately into service, if any of the servers are free, or else joins the queue if alls servers are busy. After service, a customer leaves the system and if there are any customers waiting, the next one in line enters the service. The successive service times are assumed to be independent exponential random variables with mean 1/μ.Denoting byX(t) the number in the system at timet, then {X(t),t≥0} is a birth-and-death process with μn=nμ, if 1 ≤n≤s; μn=sμ, ifn > s; and λn=λ, if n≥0. Defining ρ=λ/sμ as the traffic intensity (also known as the server utilization), it is required that ρ < 1, for the queue to be stable, which is a common assumption in queuing theory [1]. The traffic intensityρ may also be seen as the average proportion of time that each of the servers is occupied (assuming that the jobs finding more than one vacant server choose their servers randomly). If the traffic intensity is such that ρ<1, then the system is under equilibrium and it has a stationary distribution. From the theory developed for birth-and-death processes, the stationary distribution of the number of customersN in the system at the departure epochs is given by (for instance, see Gross et al. [1])(1)Pn≡PN=n=sρnn!P0,0<n≤s,ssρns!P0,n>s,in which P0≡PN=0 is given from the usual boundary condition that the probabilities must sum to 1; that is,(2)P0=∑j=0s-1sρjj!+sρss!11-ρ-1.


### 2.2. Maximum Likelihood Estimator for the Traffic Intensity

In this section, the traffic intensityρ is the focus; that is, the ratio ρ=λ/(sμ). To generate data, it is necessary to observe the system at arbitrary epochs. Such a scheme ensures that the data-generating process is consistent with the distribution given by (1). Assume that the number of customers in the queue is xi. If x={x1,x2,...,xl} constitutes our sample of size l, then the corresponding likelihood function is(3)Lx∣ρ=∏i=1lsρxixi!P0I0≤xi≤s+ssρxis!P0Ixi>s,in which I{•} is the indicator function. It is pertinent to note here that this data-generating scheme must ensure the independence of the sample observations, provided that they are sufficiently spaced in time.The maximum likelihood estimator for the traffic intensityρ is the value that maximizes the likelihood given by (3) for 0<ρ<1. Observe that (3) may be written as Lx∣ρ=α∏i=1lρxiP0, in which α is a constant, and that the value that maximizes (3) also maximizes(4)Ux∣ρ=∏i=1lρxiP0=P0nρ∑j=1lxj.The value that maximizes (4) is the maximum likelihood estimator, ρ^MLE. Applying the natural logarithm to (4), we obtain(5)ln⁡Ux∣ρ=nln⁡∑j=0s-1sρjj!+ssρss!1-ρ+∑i=0lxiln⁡ρ.Making the derivative of (5) equal to zero will produce(6)ux∣ρ=-nP0∑j=0s-1sjρj-1j-1!+sss!s1-ρρs-1+ρs1-ρ2+∑i=0lxiρ=0.The value that solves (6) is the maximum likelihood estimator, ρ^MLE. However, because (6) is not algebraically simple, it is necessary to implement a numerical search to find ρ^MLE, as we will show shortly. For our convenience, the software selected was MATLAB® [35], but many others could be used as well.To better understand the behavior of (6), it is interesting to closely analyze some special cases ofs. Table 1 illustrates some of these special cases fors = 1, 2, 3, and 4, in which x-=∑i=1lxi/l. Observe that each polynomial is of a degree equal to the number of servers,s. Figure 1 illustrates the likelihood function, L(x∣ρ), and the derivative of its natural logarithm with respect to ρ, u(x∣ρ), fors = 4 and sample value x-=4.03. The maximum likelihood estimated by the software is ρ^MLE=0.7180. It is noticeable that the functions were well behaved and seem to be unimodal. Other examples for different values of s were tested, and the results (not shown) were similar.Table 1
Special cases for (6).
s Eq. (6) 1 ρ x - + 1 - x - = 0 2 ρ 2 x - + 2 ρ - x - = 0 3 ρ 3 3 x - - 3 + ρ 2 6 + x - + ρ 6 + 2 x - - x - = 0 4 ρ 4 x - 2 - 1 + ρ 3 x - 4 + 3 4 + ρ 2 - 3 x - 16 + 3 2 + ρ 3 4 - 3 x - 8 - 3 x - 16 = 0Figure 1
Behavior of the likelihood function,L(x∣ρ), and the derivative of its natural logarithm with respect to ρ, u(x∣ρ), fors = 4, and sample value x-=4.03.


### 2.3. Traffic Intensity Moment Estimator versus the Maximum Likelihood Estimator

The moment estimator [36] for ρ can be obtained by setting the analytical expression of the expected number of clients in the queuing system,L, equal to the average of the sample of the number of clients at the arbitrary epochs, x-. In accordance with Gross et al. [1], the expected number of clients in anM/M/s queuing system is given by(7)L=ρρsss!1-ρ2P0+sρ.To obtain the moment estimator, ρ^ME, we set (7) equal to the average of the sample, x-. In other words, the moment estimator is the value ρ^ME that satisfies(8)L=ρρsss!1-ρ2P0+sρ=x-,which should be solved preferably by numerical methods since finding an analytical solution does not seem to be trivial.After some computational experiments (not shown), the authors observed that the results from the maximum likelihood and moment estimators were actually equal, which led to the conjecture that, inM/M/s queuing systems, the maximum likelihood estimator and the moment estimator for the traffic intensity ρ are equivalent. Indeed, it is easy to verify that both estimators are equal by showing that an algebraic manipulation of (6) will produce (7):(9)P0ρs∑j=0s-1sj-1ρj-1j-1!+ss-1ρs-1s-1!1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-2syρyy!+ss-1ρs-1s-1!1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-2syρyy!+ss-1ρs-1s-1!1+ρ1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-1syρyy!+ss-1ρs-1s-1!ρ1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-1sρyy!+sρss!1-ρ+ρsρss!1-ρ2=x-⇒L=sρ+ρsρss!1-ρ2P0=x-.Thus, the moment estimator is in fact equal to the maximum likelihood estimator, which is very convenient because, in general, obtaining the moment estimator is simpler than obtaining the likelihood estimator and this facilitates the explanation of the estimation process. In addition, in this particular case, the moment estimator for the traffic intensity inherits the good asymptotic properties shared by all maximum likelihood estimators.


### 2.4. Arrival and Service Rate Estimators

If the samplex={x1,x2,…,xl} related to the number of clients in the system at arbitrary epochs is collected and also the sojourn times of the clients, it is possible to estimate the arrival rateλ and the service rateμ as well. Define the sojourn time of a client in the system, which includes both the waiting time in line and the service time. In accordance with Gross et al. [1], the expected sojourn timeW may be expressed as(10)W=ρsss!sμ1-ρ2P0+1μ.With a sample of the sojourn times in the system, it is possible to estimateW as x-W. Also with the estimate ρ^, which was obtained in the previous section, it is possible to estimateμ as(11)μ^=ρ^ssP^0+s!s1-ρ^2x-Ws!s1-ρ2.Notice that in (11) P^0 may be calculated from (2) using the estimated traffic intensity, ρ^. Then, from the traffic intensity, ρ=λ/sμ, it is possible to estimate λ as(12)λ^=sμ^ρ^.


## 2.1. Probabilistic Model

InM/M/s queues, the arrivals of the customers are assumed to follow a Poisson process with rateλ, which means that the times between successive arrivals are independent exponentially distributed random variables with mean 1/λ. Upon arrival, a customer goes immediately into service, if any of the servers are free, or else joins the queue if alls servers are busy. After service, a customer leaves the system and if there are any customers waiting, the next one in line enters the service. The successive service times are assumed to be independent exponential random variables with mean 1/μ.Denoting byX(t) the number in the system at timet, then {X(t),t≥0} is a birth-and-death process with μn=nμ, if 1 ≤n≤s; μn=sμ, ifn > s; and λn=λ, if n≥0. Defining ρ=λ/sμ as the traffic intensity (also known as the server utilization), it is required that ρ < 1, for the queue to be stable, which is a common assumption in queuing theory [1]. The traffic intensityρ may also be seen as the average proportion of time that each of the servers is occupied (assuming that the jobs finding more than one vacant server choose their servers randomly). If the traffic intensity is such that ρ<1, then the system is under equilibrium and it has a stationary distribution. From the theory developed for birth-and-death processes, the stationary distribution of the number of customersN in the system at the departure epochs is given by (for instance, see Gross et al. [1])(1)Pn≡PN=n=sρnn!P0,0<n≤s,ssρns!P0,n>s,in which P0≡PN=0 is given from the usual boundary condition that the probabilities must sum to 1; that is,(2)P0=∑j=0s-1sρjj!+sρss!11-ρ-1.


## 2.2. Maximum Likelihood Estimator for the Traffic Intensity

In this section, the traffic intensityρ is the focus; that is, the ratio ρ=λ/(sμ). To generate data, it is necessary to observe the system at arbitrary epochs. Such a scheme ensures that the data-generating process is consistent with the distribution given by (1). Assume that the number of customers in the queue is xi. If x={x1,x2,...,xl} constitutes our sample of size l, then the corresponding likelihood function is(3)Lx∣ρ=∏i=1lsρxixi!P0I0≤xi≤s+ssρxis!P0Ixi>s,in which I{•} is the indicator function. It is pertinent to note here that this data-generating scheme must ensure the independence of the sample observations, provided that they are sufficiently spaced in time.The maximum likelihood estimator for the traffic intensityρ is the value that maximizes the likelihood given by (3) for 0<ρ<1. Observe that (3) may be written as Lx∣ρ=α∏i=1lρxiP0, in which α is a constant, and that the value that maximizes (3) also maximizes(4)Ux∣ρ=∏i=1lρxiP0=P0nρ∑j=1lxj.The value that maximizes (4) is the maximum likelihood estimator, ρ^MLE. Applying the natural logarithm to (4), we obtain(5)ln⁡Ux∣ρ=nln⁡∑j=0s-1sρjj!+ssρss!1-ρ+∑i=0lxiln⁡ρ.Making the derivative of (5) equal to zero will produce(6)ux∣ρ=-nP0∑j=0s-1sjρj-1j-1!+sss!s1-ρρs-1+ρs1-ρ2+∑i=0lxiρ=0.The value that solves (6) is the maximum likelihood estimator, ρ^MLE. However, because (6) is not algebraically simple, it is necessary to implement a numerical search to find ρ^MLE, as we will show shortly. For our convenience, the software selected was MATLAB® [35], but many others could be used as well.To better understand the behavior of (6), it is interesting to closely analyze some special cases ofs. Table 1 illustrates some of these special cases fors = 1, 2, 3, and 4, in which x-=∑i=1lxi/l. Observe that each polynomial is of a degree equal to the number of servers,s. Figure 1 illustrates the likelihood function, L(x∣ρ), and the derivative of its natural logarithm with respect to ρ, u(x∣ρ), fors = 4 and sample value x-=4.03. The maximum likelihood estimated by the software is ρ^MLE=0.7180. It is noticeable that the functions were well behaved and seem to be unimodal. Other examples for different values of s were tested, and the results (not shown) were similar.Table 1
Special cases for (6).
s Eq. (6) 1 ρ x - + 1 - x - = 0 2 ρ 2 x - + 2 ρ - x - = 0 3 ρ 3 3 x - - 3 + ρ 2 6 + x - + ρ 6 + 2 x - - x - = 0 4 ρ 4 x - 2 - 1 + ρ 3 x - 4 + 3 4 + ρ 2 - 3 x - 16 + 3 2 + ρ 3 4 - 3 x - 8 - 3 x - 16 = 0Figure 1
Behavior of the likelihood function,L(x∣ρ), and the derivative of its natural logarithm with respect to ρ, u(x∣ρ), fors = 4, and sample value x-=4.03.


## 2.3. Traffic Intensity Moment Estimator versus the Maximum Likelihood Estimator

The moment estimator [36] for ρ can be obtained by setting the analytical expression of the expected number of clients in the queuing system,L, equal to the average of the sample of the number of clients at the arbitrary epochs, x-. In accordance with Gross et al. [1], the expected number of clients in anM/M/s queuing system is given by(7)L=ρρsss!1-ρ2P0+sρ.To obtain the moment estimator, ρ^ME, we set (7) equal to the average of the sample, x-. In other words, the moment estimator is the value ρ^ME that satisfies(8)L=ρρsss!1-ρ2P0+sρ=x-,which should be solved preferably by numerical methods since finding an analytical solution does not seem to be trivial.After some computational experiments (not shown), the authors observed that the results from the maximum likelihood and moment estimators were actually equal, which led to the conjecture that, inM/M/s queuing systems, the maximum likelihood estimator and the moment estimator for the traffic intensity ρ are equivalent. Indeed, it is easy to verify that both estimators are equal by showing that an algebraic manipulation of (6) will produce (7):(9)P0ρs∑j=0s-1sj-1ρj-1j-1!+ss-1ρs-1s-1!1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-2syρyy!+ss-1ρs-1s-1!1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-2syρyy!+ss-1ρs-1s-1!1+ρ1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-1syρyy!+ss-1ρs-1s-1!ρ1-ρ+ρsρss!1-ρ2=x-⇒P0ρs∑y=0s-1sρyy!+sρss!1-ρ+ρsρss!1-ρ2=x-⇒L=sρ+ρsρss!1-ρ2P0=x-.Thus, the moment estimator is in fact equal to the maximum likelihood estimator, which is very convenient because, in general, obtaining the moment estimator is simpler than obtaining the likelihood estimator and this facilitates the explanation of the estimation process. In addition, in this particular case, the moment estimator for the traffic intensity inherits the good asymptotic properties shared by all maximum likelihood estimators.


## 2.4. Arrival and Service Rate Estimators

If the samplex={x1,x2,…,xl} related to the number of clients in the system at arbitrary epochs is collected and also the sojourn times of the clients, it is possible to estimate the arrival rateλ and the service rateμ as well. Define the sojourn time of a client in the system, which includes both the waiting time in line and the service time. In accordance with Gross et al. [1], the expected sojourn timeW may be expressed as(10)W=ρsss!sμ1-ρ2P0+1μ.With a sample of the sojourn times in the system, it is possible to estimateW as x-W. Also with the estimate ρ^, which was obtained in the previous section, it is possible to estimateμ as(11)μ^=ρ^ssP^0+s!s1-ρ^2x-Ws!s1-ρ2.Notice that in (11) P^0 may be calculated from (2) using the estimated traffic intensity, ρ^. Then, from the traffic intensity, ρ=λ/sμ, it is possible to estimate λ as(12)λ^=sμ^ρ^.


## 3. Performance Analysis of the Estimation Methods

### 3.1. Simulations

In this section, the performances of the estimation methods described earlier are evaluated. The algorithms were implemented in MATLAB [35], and the codes are available directly from the authors upon request for research and teaching purposes. A Monte Carlo simulation study was conducted as follows to evaluate the methodology described earlier in this paper.Random samples of the number of customers in the arbitrary epochs were obtained from simulations of anM/M/s queuing system, with sample sizes l={50,200,400}, traffic intensities ρ={0.1,0.5,0.9}, and number of servers s={3,9,27}. Without loss of generality, the arrival rate was fixed to λ = 10 for all experiments, where the service rate is given by μ=λ/sρ. With these values, it is possible to evaluate a broad range of different situations that are likely to occur in practice.Table2 shows the average values (means) and the standard deviations (SDs) of estimates for ρ, μ, and λ, from (8), (11), and (12), respectively, obtained from a Monte Carlo simulation procedure, repeated 1,000 times. From Table 2, it is observed that, with the increase of the sample size, the difference between the estimates and the true values reduces and the standard deviation decreases.Table 2
Averages (means) and standard deviations (SDs) of estimates forρ, λ, and μ.
λ=10 Sampling distribution ℓ ρ s μ ρ estimated μ estimated λ estimated Mean SD Mean SD Mean SD 50 0.1 3 33.3333 0.1006 0.0253 36.0326 4.7752 10.4347 1.4326 9 11.1111 0.0998 0.0150 11.9889 1.6397 10.3399 1.4758 27 3.7037 0.0996 0.0090 4.0059 0.5537 10.3394 1.4949 0.5 3 6.6667 0.4979 0.0452 7.1426 1.0592 10.2850 1.5888 9 2.2222 0.4992 0.0320 2.4031 0.3244 10.3540 1.4599 27 0.7407 0.4990 0.0186 0.8015 0.1112 10.3596 1.5006 0.9 3 3.7037 0.8744 0.0138 3.3016 0.5743 9.1311 1.5507 9 1.2346 0.8643 0.0137 1.0493 0.1328 8.7706 1.0757 27 0.4115 0.8661 0.0126 0.4256 0.0630 9.7771 1.5296  200 0.1 3 33.3333 0.1002 0.0130 35.9311 2.3430 10.3781 0.7029 9 11.1111 0.0998 0.0073 11.9786 0.8625 10.3374 0.7762 27 3.7037 0.1001 0.0042 3.9896 0.2456 10.3651 0.6631 0.5 3 6.6667 0.4990 0.0228 7.1038 0.4836 10.2856 0.7254 9 2.2222 0.4992 0.0167 2.3884 0.1575 10.3282 0.7088 27 0.7407 0.4994 0.0093 0.7973 0.0516 10.3430 0.6965 0.9 3 3.7037 0.8755 0.0067 3.2548 0.2719 9.0676 0.7341 9 1.2346 0.8649 0.0068 1.0335 0.0616 8.6964 0.4987 27 0.4115 0.8668 0.0064 0.4234 0.0310 9.7694 0.7543  400 0.1 3 33.3333 0.1001 0.0090 35.9298 1.6050 10.3695 0.4815 9 11.1111 0.1001 0.0051 11.9613 0.5730 10.3618 0.5157 27 3.7037 0.1002 0.0029 3.9871 0.1878 10.3697 0.5071 0.5 3 6.6667 0.4986 0.0158 7.0955 0.3618 10.2760 0.5427 9 2.2222 0.5001 0.0120 2.3877 0.1108 10.3486 0.4985 27 0.7407 0.4991 0.0064 0.7985 0.0377 10.3397 0.5085 0.9 3 3.7037 0.8758 0.0047 3.2426 0.1902 9.0476 0.5136 9 1.2346 0.8652 0.0045 1.0294 0.0423 8.6808 0.3424 27 0.4115 0.8671 0.0042 0.4229 0.0209 9.7613 0.5083


### 3.2. Numerical Example

Now an application based on data collected in a large supermarket network in a region of interest is discussed for a better understanding of the presented methodology. The cashier's desks in this network are separated into normal desks and fast desks. The focus of our example is the latter, which serve exclusively consumers with up to 15 items. The primary interest of those responsible for the supermarket is to evaluate the time between 15:00 and 22:00 on Saturdays, mainly because there is a large volume of buyers at that time, and the volume has increased because the supermarket chain has an advertising agency conducting regular campaigns. The fear is that users will give up shopping or not return to this supermarket because of dissatisfaction with the sojourn time. The past experiences of those responsible for supermarkets indicate that the traffic intensityρ should not exceed 85% to minimize the risk of compromising current and future sales. Based on several previous studies, it is reasonable to assume that consumers arrive during a given period of time approximately according to a Poisson process and that the service times are exponentially distributed. The objective is to evaluate the traffic intensity, arrival rates, and the service rates for anM/M/6 queuing system, when the current number of servers (that is, 6) is insufficient.In this scenario, 200 random observations of the number of customers in the system were collected, 28 at each of 8 consecutive Saturdays, at times sufficiently spaced and previously defined by the person responsible for data collection. The observed values (O) and frequencies (F) are presented in Table3. For example, from the 200 observations, at 2 times, only 1 customer was found in the system and at 6 times, 2 customers were found.Table 3
Observed values and frequencies.
O F O F O F O F 1 2 7 14 13 4 19 4 2 6 8 14 14 11 20 6 3 14 9 11 15 8 21 4 4 23 10 8 16 7 22 8 5 18 11 10 17 3 23 2 6 12 12 8 18 2 25 1Along with the observation of the number of clients in the system, the sojourn times for 200 randomly selected customers were also collected (which was performed simultaneously with the observation of the number of customers in the system). Figure2 depicts the histogram for the observed results.Figure 2
Histogram of the observed sojourn times for 200 clients.Using the collected data, the estimates from (8), (11), and (12) are ρ^ =0.9013, μ^ =2.0046 clients/min, and λ^ =10.8247 clients/min, respectively. The data suggests thats = 6 servers are insufficient to handle the clients and to keep the traffic intensity below 85%. Using the previous estimates for μ and λ, we can verify that increasing the number of servers tos = 7 would reduce the traffic intensity to ρ^=λ^/sμ^ =10.8247/(7×2.0046) = 0.7714, which would result in an acceptable performance for the queuing system (that is, below 85%).


## 3.1. Simulations

In this section, the performances of the estimation methods described earlier are evaluated. The algorithms were implemented in MATLAB [35], and the codes are available directly from the authors upon request for research and teaching purposes. A Monte Carlo simulation study was conducted as follows to evaluate the methodology described earlier in this paper.Random samples of the number of customers in the arbitrary epochs were obtained from simulations of anM/M/s queuing system, with sample sizes l={50,200,400}, traffic intensities ρ={0.1,0.5,0.9}, and number of servers s={3,9,27}. Without loss of generality, the arrival rate was fixed to λ = 10 for all experiments, where the service rate is given by μ=λ/sρ. With these values, it is possible to evaluate a broad range of different situations that are likely to occur in practice.Table2 shows the average values (means) and the standard deviations (SDs) of estimates for ρ, μ, and λ, from (8), (11), and (12), respectively, obtained from a Monte Carlo simulation procedure, repeated 1,000 times. From Table 2, it is observed that, with the increase of the sample size, the difference between the estimates and the true values reduces and the standard deviation decreases.Table 2
Averages (means) and standard deviations (SDs) of estimates forρ, λ, and μ.
λ=10 Sampling distribution ℓ ρ s μ ρ estimated μ estimated λ estimated Mean SD Mean SD Mean SD 50 0.1 3 33.3333 0.1006 0.0253 36.0326 4.7752 10.4347 1.4326 9 11.1111 0.0998 0.0150 11.9889 1.6397 10.3399 1.4758 27 3.7037 0.0996 0.0090 4.0059 0.5537 10.3394 1.4949 0.5 3 6.6667 0.4979 0.0452 7.1426 1.0592 10.2850 1.5888 9 2.2222 0.4992 0.0320 2.4031 0.3244 10.3540 1.4599 27 0.7407 0.4990 0.0186 0.8015 0.1112 10.3596 1.5006 0.9 3 3.7037 0.8744 0.0138 3.3016 0.5743 9.1311 1.5507 9 1.2346 0.8643 0.0137 1.0493 0.1328 8.7706 1.0757 27 0.4115 0.8661 0.0126 0.4256 0.0630 9.7771 1.5296  200 0.1 3 33.3333 0.1002 0.0130 35.9311 2.3430 10.3781 0.7029 9 11.1111 0.0998 0.0073 11.9786 0.8625 10.3374 0.7762 27 3.7037 0.1001 0.0042 3.9896 0.2456 10.3651 0.6631 0.5 3 6.6667 0.4990 0.0228 7.1038 0.4836 10.2856 0.7254 9 2.2222 0.4992 0.0167 2.3884 0.1575 10.3282 0.7088 27 0.7407 0.4994 0.0093 0.7973 0.0516 10.3430 0.6965 0.9 3 3.7037 0.8755 0.0067 3.2548 0.2719 9.0676 0.7341 9 1.2346 0.8649 0.0068 1.0335 0.0616 8.6964 0.4987 27 0.4115 0.8668 0.0064 0.4234 0.0310 9.7694 0.7543  400 0.1 3 33.3333 0.1001 0.0090 35.9298 1.6050 10.3695 0.4815 9 11.1111 0.1001 0.0051 11.9613 0.5730 10.3618 0.5157 27 3.7037 0.1002 0.0029 3.9871 0.1878 10.3697 0.5071 0.5 3 6.6667 0.4986 0.0158 7.0955 0.3618 10.2760 0.5427 9 2.2222 0.5001 0.0120 2.3877 0.1108 10.3486 0.4985 27 0.7407 0.4991 0.0064 0.7985 0.0377 10.3397 0.5085 0.9 3 3.7037 0.8758 0.0047 3.2426 0.1902 9.0476 0.5136 9 1.2346 0.8652 0.0045 1.0294 0.0423 8.6808 0.3424 27 0.4115 0.8671 0.0042 0.4229 0.0209 9.7613 0.5083


## 3.2. Numerical Example

Now an application based on data collected in a large supermarket network in a region of interest is discussed for a better understanding of the presented methodology. The cashier's desks in this network are separated into normal desks and fast desks. The focus of our example is the latter, which serve exclusively consumers with up to 15 items. The primary interest of those responsible for the supermarket is to evaluate the time between 15:00 and 22:00 on Saturdays, mainly because there is a large volume of buyers at that time, and the volume has increased because the supermarket chain has an advertising agency conducting regular campaigns. The fear is that users will give up shopping or not return to this supermarket because of dissatisfaction with the sojourn time. The past experiences of those responsible for supermarkets indicate that the traffic intensityρ should not exceed 85% to minimize the risk of compromising current and future sales. Based on several previous studies, it is reasonable to assume that consumers arrive during a given period of time approximately according to a Poisson process and that the service times are exponentially distributed. The objective is to evaluate the traffic intensity, arrival rates, and the service rates for anM/M/6 queuing system, when the current number of servers (that is, 6) is insufficient.In this scenario, 200 random observations of the number of customers in the system were collected, 28 at each of 8 consecutive Saturdays, at times sufficiently spaced and previously defined by the person responsible for data collection. The observed values (O) and frequencies (F) are presented in Table3. For example, from the 200 observations, at 2 times, only 1 customer was found in the system and at 6 times, 2 customers were found.Table 3
Observed values and frequencies.
O F O F O F O F 1 2 7 14 13 4 19 4 2 6 8 14 14 11 20 6 3 14 9 11 15 8 21 4 4 23 10 8 16 7 22 8 5 18 11 10 17 3 23 2 6 12 12 8 18 2 25 1Along with the observation of the number of clients in the system, the sojourn times for 200 randomly selected customers were also collected (which was performed simultaneously with the observation of the number of customers in the system). Figure2 depicts the histogram for the observed results.Figure 2
Histogram of the observed sojourn times for 200 clients.Using the collected data, the estimates from (8), (11), and (12) are ρ^ =0.9013, μ^ =2.0046 clients/min, and λ^ =10.8247 clients/min, respectively. The data suggests thats = 6 servers are insufficient to handle the clients and to keep the traffic intensity below 85%. Using the previous estimates for μ and λ, we can verify that increasing the number of servers tos = 7 would reduce the traffic intensity to ρ^=λ^/sμ^ =10.8247/(7×2.0046) = 0.7714, which would result in an acceptable performance for the queuing system (that is, below 85%).


## 4. Conclusions and Final Remarks

The estimation of the traffic intensityρ, the arrival rateλ, and the service rateμ, for Markovian multiservice queues, known in Kendall notation asM/M/s queues, is a challenging mathematical problem that was solved here by statistical inferential methods. For theM/M/s queues, even without knowledge of the arrival rate and queue length of service, it is possible to make inferences about the traffic intensity by collecting a random sample of the number of customers in the system at the arbitrary epochs. If, in addition to the number of customers, the sojourn times in the queue are also collected, then it is possible to accurately estimate the arrival and service rates by means of simple expressions, as demonstrated here. In fact, a comprehensive set of Monte Carlo simulations were performed to show that estimates forρ,μ, andλ may be obtained that are very close to the true values, especially when the sample sizes are close to 400.Finally, it is worthwhile to mention that future studies in this area could include extensions to more general single-server queuing systems, such asM/G/1 andG/M/1, and general multiserver queues, such asM/G/s andG/M/s. These are only a few of the many topics for future research in this area.



---

*Source: 3280846-2018-12-25.xml*