Performance Evaluation and Dimensioning of 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 Systems through Kernel Estimation

Abstract

We extend the analysis of queueing systems for real-life situations, where the arrival pattern of customers is unknown. In real systems, we must understand how the choice of a method of estimation influences the configuration of the system. Using kernel smoothing, we evaluate algorithms to estimate performance measures of a 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 system, including the invariant probability distribution of the number of customers in the system, the blocking probability, the average queue size, and the average client queue time. We successfully apply the method to the arrivals to a call center to plan and improve the performance of these important queueing systems.

1. Introduction

There is a large practical interest in investigating the behavior of general-arrival queueing systems, namely, those of 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 type, because when managing real queueing systems, the behavior of the arrival process generally is not Markovian. In Kendall [1] notation, in these queueing systems, the interarrival times are independent and follow a general distribution (𝐺𝐼), the service times follow an exponential distribution (𝑀), we have 𝑐 identical servers working in parallel, and a maximum capacity of 𝑁 users that are simultaneously allowed in the system, including those in service. Finally, 𝑋 is a random variable representing the size of group (bulk) arrivals. Such queueing systems could be used in situations where we have relative control over how the servers work, but we do not know beforehand how the customers arrive at the system.

Naturally, the mathematical model depends on the type of queueing system considered, and there are several methods to obtain such models. The most widely used methods are those that attempt to explain the density functions of interarrival and service times by means of parametric statistical models. Nevertheless, real data sometimes do not fit well into parametric models; instead, they produce intractable models. Exact results for performance evaluation of Markovian and some simple general queueing systems are known [2], but such systems sometimes are not found in real life.

An alternative approach is using nonparametric methods to study queueing systems. Nonparametric methods that use kernel smoothing have received much attention lately [3]. Kernel estimators provide a simple way of finding structure in data sets without imposing a specific parametric model [4], which gives us flexibility to handle virtually any data set. There is extensive literature discussing queueing systems [5] and kernel smoothing [4] as separate concepts, but virtually no study has brought the two concepts together.

The contribution of this paper is twofold. First, we implement algorithms to calculate performance measures of queueing systems, where the density of the interarrival time is estimated by the kernel method. Second, we evaluate the performance of these algorithms as a function of the kernel estimator, the smoothing window, the intensity rate, and the system size. We also present a case study in a call center that illustrates the usefulness of the method.

First, we present a literature review about queueing systems and kernel estimators as well as the fundamental concepts required to understand the proposed queueing system model. We also discuss the use of kernel estimators, the chosen models, and the issue of selecting the smoothing parameter. Then, we describe how we estimated system performance and present the comparative results of simulations with the different methods discussed. We then apply our methods to a call center case and end with our main conclusions and some ideas for future work in the area.

2. Literature Review and Fundamental Concepts

2.1. Previous Works

There are many situations in real life where queues occur and queueing models may be helpful. Recently, queueing models have been used successfully in manufacturing processes [6, 7], transportation [8], airports, ports, and product distribution systems [9], computer and telecommunication systems [10], call center modeling [11], and the analysis of health systems [12]. Queues may cause the quality of the services or the prices of the goods to rise or fall, depending on their efficiency [8], which may be estimated by means of the mathematical tools developed in queueing theory.

As we mentioned earlier, there is not much literature that incorporates both general arrival queueing systems and kernel smoothing. Takács [13] analyzed a closed solution for various systems that have nonspecific distributions, including some multiserver queues such as 𝐺𝐼/𝑀/𝑐 and 𝑀/𝐺/𝑐. Hokstad [14] established some closed form results to the 𝐺𝐼/𝑀/𝑐/𝑁 system. Chaudhry and Templeton [15] analyzed various types of queue with bulk arrivals. Zhao [16] proposed a closed form solution for the 𝐺𝐼^𝑋/𝑀/𝑐 system, and Bareche and Aïssani [17] proposed a method to evaluate the proximity of 𝐺𝐼/𝑀/1 and 𝑀/𝑀/1 systems when the density of the interarrival time is estimated by kernel estimators.

Concerning kernel smoothing, Wand and Jones [4] introduced general fundamental concepts. Regarding the issues related to the asymmetry of the random variable under analysis (e.g., nonnegativity), Zhang et al. [18] proposed a boundary corrected kernel estimator based on pseudodata generation, transformation, and reflection around the 𝑌-axis. Chen [19] proposed the use of a gamma kernel to avoid boundary problems present in certain situations. Scaillet [20] studied the application of other asymmetric kernels. Bouezmarni and Scaillet [21] were concerned about the consistency of these asymmetrical estimators. For recent developments in the area of kernel smoothing and a thorough literature review, see Atuncar et al. [22] and de Lima and Atuncar [3].

2.2. The 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 Model with Partial Blocks

Vijaya Laxmi and Gupta [23] described a generalization of the system 𝐺𝐼/𝑀/𝑐/𝑁 when customers arrive in groups of size 𝑋 with 𝑃(𝑋=𝑖)=𝑔_𝑖 (𝑖≥1) and mean 𝐸(𝑋)=𝑔. The 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 is a finite capacity system such that a customer that arrives to the saturated system is refused with a probability that we will call 𝑃_BL, or partial blocks. This refers to the case in which an arrival group of a size greater than the remaining spots in the system is partially denied according to the number of the remaining vacancies until the system is complete. Let 𝐵_𝑛 be the number of clients who were served between the arrival of the 𝑛th customer and its successor. Therefore, the number of clients the 𝑛th customer finds in the system at the arrival, 𝑌_𝑛, would depend on 𝑋_𝑛 and 𝐵_𝑛 such that 𝑌_𝑛+1=[min(𝑌_𝑛+𝑋_𝑛,𝑁)−𝐵_𝑛]⁺. Because 𝑌_𝑛+1 depends only on 𝑌_𝑛, 𝐵_𝑛, and 𝑋_𝑛 but not on 𝑌_𝑛−1, 𝑌_𝑛−2, and so on, the stochastic process {𝑌_𝑛} is a first-order Markov Chain.

Vijaya Laxmi and Gupta [23] report that when the traffic intensity rate 𝜌=(𝜆𝑔/𝑐𝜇), where 𝜆 is the arrival rate and 𝜇 service rate, is smaller than 1, this Markov Chain has an invariant probability distribution of 𝜋_𝑘=lim_𝑛→∞𝑃(𝑌_𝑛=k), 𝑘=0,1,2,… associated with the number of clients an arbitrary customer finds in the system at arrival. The 𝜋_𝑘’s are often called prearrival probabilities.

Prearrival probabilities can be determined by the following system of linear equations:𝜋_𝑘=^𝑁_𝑗=0𝑝_𝑗𝑘𝜋_𝑗,(2.1) for 𝑘=0,1,2,…,𝑁−1, and^𝑁_𝑗=0𝜋_𝑗=1,(2.2) which gives⎡⎢⎢⎢⎢⎣𝑝^0,0−1𝑝_1,0⋯𝑝_𝑁,0𝑝𝑝⋮_0,1𝑝_1,1−1⋱⋮⋮_0,𝑁−11⋯⋯⋯⋯𝑝1_{𝑁,𝑁−1}−1⎤⎥⎥⎥⎥⎦⎡⎢⎢⎢⎣𝜋⁰𝜋₁𝜋⋮⋮_𝑁⎤⎥⎥⎥⎦=⎡⎢⎢⎢⎣⋮0001⎤⎥⎥⎥⎦,(2.3) where 𝑝_𝑗𝑘 are named transition probabilities such that𝑝_𝑗𝑘=𝑃𝑌_𝑛+1=𝑘∣𝑌_𝑛=𝑗=^∞_𝑖=1=𝑃𝑋_𝑛=𝑖,𝑌_𝑛+1=𝑘∣𝑌_𝑛=𝑗_∞_𝑖=1𝑔_𝑖𝑃𝑌_𝑛+1=𝑘∣𝑌_𝑛=𝑗,𝑋_𝑛=𝑖=^∞_𝑖=1=𝑔_𝑖𝑃𝐵_𝑛=min(𝑗+𝑖,𝑁)−𝑘_𝑁−𝑗_𝑖=ma_{x(1,𝑗−𝑘)}𝑔_𝑖𝑃𝐵_𝑛=𝑗+𝑖−𝑘+^∞_𝑖=𝑁_−𝑗+1=𝑔_𝑖𝑃𝐵_𝑛=𝑁−𝑘_𝑁−𝑗_𝑖=ma_{x(1,𝑗−𝑘)}𝑔_𝑖^∞0𝜋_{𝑗+𝑖,𝑘}(𝑧)𝑑𝜏(𝑧)+^∞_𝑖=𝑁_−𝑗+1𝑔_𝑖^∞0𝜋_𝑁,𝑘(𝑧)𝑑𝜏(𝑧),(2.4) where, 𝑗≤𝑁, 𝑘≤𝑁, 𝜋_{𝑗+𝑖,𝑘}(𝑧) is the serving probability of (𝑗−𝑘+𝑖) clients under the assumption of the interarrival time 𝜏_𝑛=𝑧 and 𝜏(𝑧) is the interarrival time distribution.

We will analyze how 𝜋_{𝑗+𝑖,𝑘}(𝑧) behaves for all possibilities of (𝑗+1) and (𝑘) in regard to the number of servers 𝑐.

(a) (𝑗+𝑖)≥𝑐 and 𝑘≥𝑐
When 𝑗+𝑖≥𝑐 and 𝑘≥𝑐, there will be more clients than the servers can handle in the entire interval. Because the service process is Markovian, we can treat the server group as a single unit that serves customers at a rate 𝑐𝜇𝑧 and a Poisson distributed transition probability 𝜋_{𝑗+𝑖,𝑘}(𝑧)=𝑒^{−𝑐𝜇𝑧}(𝑐𝜇𝑧)^{𝑗+𝑖−𝑘}(𝑗+𝑖−𝑘)!.(2.5)

(b) (𝑗+𝑖)<𝑐
When 𝑗+𝑖<𝑐, all clients within the system are being served and only 𝑘 customers will remain in the system to time 𝑧. Knowing that the probability of a service time greater than 𝑧 is 𝑒^−𝜇𝑧, we can describe this transition probability as a Binomial distribution 𝜋_{𝑗+𝑖,𝑘}(𝑧)=𝑗+𝑘𝑖𝑒^{−𝑘𝜇𝑧}(1−𝑒^−𝜇𝑧)^{𝑗+𝑖−𝑘}.(2.6)

(c) (𝑗+𝑖)≥𝑐 and 𝑘<𝑐
When 𝑗+𝑖≥𝑐 and 𝑘<𝑐, there will be (𝑗+𝑖−𝑐) customers waiting and 𝑐 customers being served at the beginning of the interval, but there will be (𝑐−𝑘) free servers at the end. Let 𝑦 be the interval ending immediately before (𝑗+𝑖−𝑐+1) clients are served. If each service time is exponentially distributed with a rate 𝑐𝜇, then 𝑦 is gamma distributed with a shape (𝑗+𝑖−𝑐+1) and rate a 𝑐𝜇. The other 𝑐 customers will be served in a time (𝑧−𝑦), and only 𝑘 will remain. The transition probability of this subinterval will follow a Binomial distribution with success probability 𝑒^{−𝜇(𝑧−𝑦)}. The values 𝜋_{𝑗+𝑖,𝑘}(𝑧) can be obtained by the convolution of these two variables 𝜋_{𝑗+𝑖,𝑘}(𝑧)=^𝑧0𝑒^{−𝑐𝜇𝑦}(𝑐𝜇)^{𝑗+𝑖−𝑐+1}𝑦^{𝑗+𝑖−𝑐}(𝑗+𝑖−𝑐)!𝑘𝑐𝑒^{−𝑘𝜇(𝑧−𝑦)}1−𝑒^{−𝜇(𝑧−𝑦)}^𝑐−𝑘=𝑘𝑐𝑒𝑑𝑦_{−𝑘𝜇𝑧}^𝑧0(𝑐𝜇𝑦)^{𝑗+𝑖−𝑐}(𝑗+𝑖−𝑐)!𝑒^𝑘𝜇𝑦=𝑘𝑐𝑒𝑒^{−𝑐𝜇𝑦}(1−𝑒^𝜇𝑦𝑒^−𝜇𝑧)^𝑐−𝑘𝑐𝜇𝑑𝑦_{−𝑘𝜇𝑧}^𝑧0(𝑐𝜇𝑦)^{𝑗+𝑖−𝑐}(𝑗+𝑖−𝑐)!𝑒^−𝜇𝑦=𝑘𝑐𝑒(1−𝑒^𝜇𝑦𝑒^−𝜇𝑧)^𝑐−𝑘𝑐𝜇𝑑𝑦_{−𝑘𝜇𝑧}^𝑧0(𝑐𝜇𝑦)^{𝑗+𝑖−𝑐}(𝑗+𝑖−𝑐)!(𝑒^−𝜇𝑦−𝑒^−𝜇𝑧)^𝑐−𝑘𝑐𝜇𝑑𝑦.(2.7) Therefore, from (2.4), we can get the transitions probabilities 𝑝_𝑗𝑘 given by (2.4) when we take each transition probability 𝜋_{𝑗+𝑖,𝑘}(𝑧) as 𝜋_𝑁,𝑘(𝑧) when 𝑖>𝑁−𝑗𝑝_𝑗𝑘=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩^𝑁−𝑗_𝑖=ma_{x{1,𝑘−𝑗}}𝛽_{𝑗+𝑖−𝑘}𝑔_𝑖+𝛽_{∞𝑁−𝑘𝑘=𝑁}_−𝑗+1𝑔_𝑖,𝑘≥𝑐,_{𝑁−𝑗𝑖=ma}_{x{1,𝑘−𝑗}}𝑉_{𝑗+𝑖,𝑘}𝑔_𝑖+𝑉_{∞𝑁,𝑘𝑘=𝑁}_−𝑗+11−_𝑁𝑔_𝑖,0<𝑘<𝑐,_𝑟=1𝑝_𝑗𝑟,𝑘=0,(2.8) where 𝑉_{𝑗+𝑖,𝑘} and 𝛽_{𝑗+𝑖−𝑘} are as follows: 𝑉_{𝑗+𝑖,𝑘}=⎧⎪⎪⎨⎪⎪⎩0,𝑗+𝑘𝑖𝑒^{−𝑘𝜇𝑧}𝑗<𝑘<𝑐,^∞0(1−𝑒^−𝜇𝑧)^{𝑗+𝑖−𝑘}𝑑𝜏(𝑧),𝑘≤𝑗≤𝑐,_∞0^𝑧0𝑘𝑐𝑒^{−𝑘𝜇𝑧}(𝑐𝜇𝑦)^𝑗−𝑐(𝑗−𝑐)!𝑐𝜇(𝑒^−𝜇𝑦𝛽−𝑒^−𝜇𝑧)^𝑐−𝑘𝑑𝑦𝑑𝜏(𝑧),𝑘<𝑐<𝑗,_{𝑗+𝑖−𝑘}=^∞0𝑒^{−𝑐𝜇𝑧}(𝑐𝜇𝑧)^{𝑗+𝑖−𝑘}(𝑗+𝑖−𝑘)!𝑑𝜏(𝑧),𝑗+𝑖−𝑘≥0.(2.9) There is a relationship between the vector of prearrival probabilities 𝝅 and the vector of arbitrary time probabilities 𝐏, related to the number of people that an outside observer finds in the system. Viajaya Laxmi and Gupta [23] established a method that proved the relationship 𝑃_𝑘=⎧⎪⎪⎨⎪⎪⎩𝜌𝑐min{𝑘,𝑐}𝑔^𝑘−1_𝑖=0𝜋_𝑖∞𝑗_=𝑘−𝑖1−_𝑁𝑔_𝑗,0<𝑘≤𝑁,_𝑖=1𝑃_𝑖,𝑘=0.(2.10) Performance measures (below) are used, including the average queue length 𝐿_𝑞, blocking probability of an arbitrary customer 𝑃_BL, and average waiting time in the queue 𝑊_𝑞, to analyze the efficiency of queueing systems 𝐿_𝑞=^𝑁_𝑖=0𝑃_BL(𝑖−𝑐)𝑃_𝑖,=^𝑁_𝑖=0𝜋_𝑖∞𝑗=_𝑁−𝑖𝑔^∞−1_𝑘_=𝑗+1𝑊_𝑞𝑔_𝑗,=𝐿^𝑞𝑔𝜆1−𝑃_𝐵𝐿).(2.11)

2.3. Kernel Estimators

Suppose that we have a sample of the interarrival times, 𝑋₁,…,𝑋_𝑛, with an unknown density 𝜏(𝑡). The kernel estimator is an analytical tool that provides an effective way of revealing the structure behind such a sample.

2.3.1. Gamma Kernel Estimator

Recently, Chen [19] suggested an asymmetric kernel with naturally varying shape as a way to avoid allocating weight for negative values. The gamma kernel estimators are always nonnegative, free of boundary bias, and achieve the optimal rate of convergence for the mean square error (MSE) in the nonnegative kernel estimator class. Bouezmarni and Scaillet [21] showed that this estimator is consistent and able to avoid boundary bias. Let 𝐾_𝐺(𝑝,𝑞) be the gamma density function with parameters 𝑝 (shape) and 𝑞 (rate). The gamma kernel considered is𝐾_𝐺𝑡𝑏+1,𝑏𝑋^𝑗=𝑋^{𝑗𝑡/𝑏}𝑒^{−𝑋𝑗/𝑏}𝑏^{(𝑡/𝑏)+1}Γ[(𝑡/𝑏)+1],(2.12) where 𝑏 is a smoothing parameter that satisfies the condition 𝑏→0, 𝑛𝑏→∞ as 𝑛→∞. The gamma kernel estimator is𝜏̂(𝑡;𝑏)=𝑛^𝑛−1_𝑗=1𝐾_𝐺𝑡𝑏+1,𝑏𝑋^𝑗.(2.13) The smoothing parameter 𝑏 is critical for the overall performance of the kernel estimator considered. A small 𝑏 leads to a relatively bumpy density, while a large one results in a smooth density. There are several methods to determine the best fit, from a minimization of the mean integrated squared error (MISE) of 𝜏̂(𝑡;𝑏) to the asymptotic behavior of the MISE (AMISE).

(a) Least Squares Cross-Validation (LSCV) Method
The least squares cross-validation (LSCV) method starts from the MISE expansion MISE{𝜏̂(𝑥;𝑏)}=𝐸𝜏̂(𝑥;𝑏)²𝑑𝑥−2𝐸𝜏̂(𝑥;𝑏)𝜏(𝑥)𝑑𝑥+𝐸𝜏(𝑥)²𝑑𝑥.(2.14) The minimization of the first term is equivalent to the minimization of MISE{𝜏̂(𝑥;ℎ)}−𝐸𝜏(𝑥)²𝑑𝑥=𝐸𝜏̂(𝑥;ℎ)²𝑑𝑥−2𝜏̂(𝑥;ℎ)𝜏(𝑥)𝑑𝑥.(2.15) The right-hand side is unknown, because it depends on 𝜏. However, an unbiased estimator for this quantity is LSCV(ℎ)=𝜏̂(𝑥;ℎ)²𝑑𝑥−2𝑛^𝑛−1_𝑖=1𝜏̂_−𝑖𝑋_𝑖;ℎ,(2.16) where 𝜏̂_−𝑖=(𝑋_𝑖;ℎ) is the density estimate based on the sample with 𝑋_𝑖 deleted; this is often called the “leave-one-out” density estimator. A disadvantage of this method is that it suffers from high variation.

(b) Asymptotic Behavior of the MISE (AMISE) Method
An alternative parameter selector is to consider the asymptotic behavior of the MISE of the gamma kernel estimator. Chen [19] uses some aspects of the gamma distribution and a Taylor expansion to determine the MISE as follows: MISE(𝜏̂)=𝑏²^∞0𝑥𝜏^(𝑥)+12𝑥𝜏^+2𝑛√(𝑥)²𝑑𝑥𝑏𝜋⁻¹^∞0𝑥^−1/2𝜏(𝑥)𝑑𝑥+𝑂𝑛⁻¹𝑏^−1/2+𝑏².(2.17) The asymptotic MISE disregards the last term; therefore, the optimal 𝑏  that minimizes the leading terms above is 𝑏_AMISE=⎡⎢⎢⎢⎣2√𝜋⁻¹∫^∞0𝑥^−1/2𝜏(𝑥)𝑑𝑥4𝑛∫^∞0𝑥𝜏^(𝑥)+2⁻¹𝑥𝜏^(𝑥)²𝑑𝑥⎤⎥⎥⎥⎦^2/5,(2.18) where the functions 𝜏, 𝜏^, and 𝜏^ are unknown. These quantities are obtained from the fitted gamma density with parameters adjusted from the sample. This solution still requires further study, but our paper shows promising results.

2.3.2. Zhang et al. [18] Estimator

Zhang et al. [18] submitted a model that works particularly well when 𝜏(0)>0 and combines pseudodata creation, its transformation, and its reflection around the 𝑌-axis in the following three steps.

Step 1. Transform the original data 𝑋₁,…,𝑋_𝑛 to 𝑔(𝑋₁),…,𝑔(𝑋_𝑛), while keeping the original data, where 𝑔 is a nonnegative, continuous, and monotonically increasing function from [0,∞) to [0,∞).
Based on extensive simulations, the transformation that best suits a broad variety of densities is 𝑔(𝑥)=𝑥+𝑑𝑥²+𝐴𝑑𝑥³,(2.19) where 𝐴>1/3 and 𝑑=𝑓^(0)/𝑓(0).

Step 2. Reflect the pseudodata, 𝑔(𝑋₁),…,𝑔(𝑋_𝑛), around the origin.

Step 3. Based on the enlarged data sample, −𝑔(𝑋₁),…,−𝑔(𝑋_𝑛),𝑋₁,…,𝑋_𝑛, define the new estimator as 𝜏̂_𝑛(𝑥,ℎ)=1𝑛ℎ^𝑛_𝑗=1𝐾𝑥−ℎ𝑋^𝑗+𝐾𝑥+𝑔𝑋^𝑗ℎ,𝑥≥0,(2.20) where ℎ is a smoothing parameter and 𝐾 is a symmetric probability function with support [−1,1] like the Epanechnikov kernel 𝐾(𝑡)=341−𝑡²𝐼_[−1,1].(2.21) Notice that the transformation 𝑔 defined above is not available in practice, because 𝑑 is unknown. A good estimator can be obtained when 𝑑 is written as (𝑑/𝑑𝑥)log𝑓(𝑥)∣_𝑥=0, 𝑑_𝑛=log𝑓^𝑛(ℎ)−log𝑓_𝑛(0)ℎ,(2.22) where 𝑓_𝑛(ℎ)=𝑓^∗𝑛(ℎ)+1/𝑛², 𝑓_𝑛(0)=max(𝑓^∗𝑛(0),1/𝑛²), 𝑓^∗𝑛(ℎ)=(1/𝑛ℎ)∑^𝑛𝑗=1𝐾((𝑥−𝑋_𝑗)/ℎ), 𝑓^∗𝑛(0)=(1/𝑛ℎ₀)∑^𝑛𝑗=1𝐾₀(𝑋_𝑗/ℎ₀), and 𝐾₀ is a so-called endpoint kernel, satisfying ⁰⁻¹𝐾₀(𝑡)𝑑ℎ𝑡=1,⁰⁻¹𝑡𝐾₀(𝑡)𝑑𝑡=0,⁰⁻¹𝑡²𝐾₀(𝑡)𝑑𝑡≠0,₀=⎧⎪⎨⎪⎩∫¹⁻¹𝐾(𝑡)𝑑𝑡²∫⁰⁻¹𝐾₀(𝑡)²𝑑𝑡∫⁰⁻¹𝑡²𝐾₀(𝑡)𝑑𝑡²∫¹⁻¹𝐾(𝑡)²𝑑𝑡⎫⎪⎬⎪⎭ℎ.(2.23) Zhang et al. [18] proved that for 𝑡≥ℎ, the effect of reflected pseudodata is insignificant, and the estimator can be reduced to the Parzen-Rosenblatt estimator 𝜏̂_𝑛(𝑡,ℎ)=1𝑛ℎ^𝑛_𝑗=1𝐾𝑡−𝑋^𝑗ℎ.(2.24) They also established that ^∞0𝜏̂_𝑛(𝑡)𝑑𝑡=1+1𝑛^𝑛_𝑖=1^{−𝑔(𝑋𝑖)/ℎ−𝑋}_𝑖/ℎ𝐾(𝑧)𝑑𝑧.(2.25) Thus, 𝜏̂_𝑛(𝑡) only integrates to 1 when 𝑑_𝑛=0, so 𝑔_𝑛(𝑋_𝑖)=𝑋_𝑖, or when 𝑋_𝑖=0 for all 𝑋_𝑖’s, because 𝑔_𝑛(0)=0. However, when 𝑛→∞, both limits of the second term will eventually converge to 0 and 𝜏̂_𝑛(𝑡) will integrate to 1 asymptotically.
Zhang and Karunamuni [24] used the endpoint kernel 𝐾₀(𝑡)=12(1+𝑡)12+𝑡𝐼^[−1,0](2.26) and showed that this kernel minimizes the MSE when estimating 𝜏(0). Therefore, ℎ₀=2ℎ is approximately the optimal smoothing parameter for estimating 𝜏(0) except when 𝜏(0)=0.
Chiu [25] described a parameter selecting method that considered the optimal ℎ that minimizes the asymptotic MISE when 𝐾 is a symmetric probability function with up to the fourth moment being finite ℎ_AMISE=∫𝐾(𝑥)²𝑑𝑥𝑛∫𝑥²𝐾(𝑥)𝑑𝑥²∫𝜏^(𝑥)²𝑑𝑥^1/5,(2.27) where the function ∫𝜏^(𝑥)²𝑑𝑥 is unknown.
Chiu’s [25] “plug-in” method consists of estimating this quantity through the characteristic function of the sample 𝜑(𝜆)=𝑛^𝑛−1_𝑗=1𝑒^{𝑖𝜆𝑋𝑗},(2.28) and calculating the optimal ℎ with the previous formula. The characteristic function of 𝜏 is 𝜑(𝜆)=𝑒^𝑖𝜆𝑥𝜏(𝑥)𝑑𝑥.(2.29) By the inversion formula, we have 𝜏(𝑥)=(2𝜋)⁻¹𝑒^𝑖𝜆𝑥𝜑(𝜆)𝑑𝜆,(2.30) thus, 𝜏^(𝑥)=(2𝜋)⁻¹𝜆²𝜏𝑒^𝑖𝜆𝑥𝜑(𝜆)𝑑𝜆,_(𝑥)²𝑑𝑥=(2𝜋)⁻¹𝜆²=𝜆𝑒^𝑖𝜆𝑥𝜑(𝜆)𝑑𝜆²𝑑𝑥₄(2𝜋)⁻¹𝑒^𝑖𝜆𝑥=𝜆𝜑(𝜆)𝑑𝜆²𝑑𝑥₄[𝜏(𝑥)]²𝑑𝑥.(2.31) Using Parseval’s identity, we can show that 𝜆⁴[𝜏(𝑥)]²𝑑𝑥=(2𝜋)⁻¹𝜆⁴||𝜑(𝜆)||²𝑑𝜆.(2.32) Chiu [25] introduced a cutoff value Λ for 𝜆 such that |𝜑(𝜆)|²<𝑐/𝑛. Extensive computational experiments from Bessegato et al. [26] show that 𝑐=3 is the value that minimizes the estimator variance. The final “plug in” estimator is then ℎ_AMISE=⎡⎢⎢⎣∫𝐾(𝑥)²𝑑𝑥𝑛∫𝑥²𝐾(𝑥)𝑑𝑥²𝜋⁻¹∫^Λ0𝜆⁴||𝜑(𝜆)||²−𝑛⁻¹𝑑𝜆⎤⎥⎥⎦^1/5.(2.33)

3. Experimental Results

This section presents some results of simulations for 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 systems with partial blocks, where the interarrival time is estimated through the following kernel methods:(i)gamma kernel estimator with LSCV method,(ii)gamma kennel estimator with optimal 𝑏_AMISE,(iii)Zhang et al. [18] estimator with Chiu’s [25] “plug in” method.

To evaluate the performance of the estimators above, we will compare the mean square error (MSE) of each estimated prearrival probability in the following way.

Step 1. Generate a sample of size 𝑛 of general interarrival distribution 𝜏.

Step 2. Calculate the mean service rate 𝜇=𝑔/𝜌𝑐E(𝜏).

Step 3. Estimate the optimal smoothing parameter ℎ or 𝑏.

Step 4. Use the kernel density method to estimate the theoretical density function 𝜏(𝑥).

Step 5. Find each estimated transition probability.

Step 6. Solve the linear system: ⎡⎢⎢⎢⎢⎣𝑝̂^0,0−1𝑝̂_1,0⋯𝑝̂_𝑁,0𝑝𝑝̂⋮̂_{0,10,𝑁−1}𝑝̂_1,1−1⋱⋮⋮1⋯⋯⋯⋯𝑝̂1_{𝑁,𝑁−1}−1⎤⎥⎥⎥⎥⎦⎡⎢⎢⎢⎣𝜋⁰𝜋₁𝜋⋮⋮_𝑁⎤⎥⎥⎥⎦=⎡⎢⎢⎢⎣⋮0001⎤⎥⎥⎥⎦.(3.1) The algorithm above was coded in R 2.8.0 (or earlier versions, see [27]) and is available upon request from the authors for educational and research purposes.
The experiments are based on random samples of interarrival times of size 𝑛=100. Slightly larger and smaller samples were also tested, but the results (not shown) are similar. The theoretical interarrival distributions considered in this experiment were(i)Weibull distribution with shape = 2 and rate = 20,(ii)gamma distribution with shape = 10 and rate = 2,(iii)gamma mixture distribution of 0.45 × gamma (5; 2) + 0.55 × gamma (30; 1).The Weibull density has 𝜏^(0)≠0, the pure gamma has 𝜏^(0)=0, and the gamma mixture is bimodal. We choose Weibull and gamma distributions because of their well-known flexibility for modeling real databases. For simplicity, the group size 𝑋 is constant and equal to 1. Two different numbers of servers were considered, 𝑐=5 and 10, with maximum capacities of 𝑁=20 and 25, respectively, (which result in buffers of fixed sizes equal to 15). Figures 1 and 2 show the results, where 𝜋_𝑘 is the theoretical distribution of the number of clients and MSE(𝜋_𝑘) is the mean square error of the estimates.

(a) Weibull (2; 20)

(b) Gamma (10; 2)

(c) 0.45 Gamma (5; 2) + 0.55 Gamma (30; 1)

Figure 1 

𝜋𝑘 and MSE(𝜋𝑘) for 𝑐=5.

(a) Weibull (2; 20)

(b) Gamma (10; 2)

(c) 0.45 Gamma (5; 2) + 0.55 Gamma (30; 1)

Figure 2 

𝜋𝑘 and MSE(𝜋𝑘) for 𝑐=10.

4. Discussion

In general, the errors decreased as the number of servers 𝑐 increased, and they were dependent on the theoretical distribution considered. We obtained the largest errors for the last distribution (gamma mixture). The LSCV method is better than 𝑏_AMISE on the first and last distributions. Zhang et al.’s [18] estimator has the worst performance on the first and a competitive performance on the second distribution (gamma distribution). All the estimators performed well with the second distribution. The last distribution (gamma mixture) has a particular behavior. Although the errors are large for most of the estimates, the error is small for the blocking probability (i.e., 𝜋₂₀, for the system with 𝑐=5, Figure 1, and 𝜋₂₅, for 𝑐=10, Figure 2). The blocking probability is an important performance measure, because it indicates the fraction of costumers lost by the system.

Finally, simulations showed that when 𝜏^(0)≠0 the gamma kernel method had the best performance. This suggests the Zhang et al. [18] method does not work well when 𝜏(0)>0. At the same time, its behavior with the bimodal density showed a very low MSE for probabilities near the maximum state. This implies a good estimation of the blocking probability and other performance measures.

The method used to select the smoothing parameter for the gamma kernel estimator had no effect. A better selector would have the function of its AMISE optimal parameter estimated like the “plug-in” method for symmetric kernels. The combined use of smoothing parameter selection and Bayesian techniques is promising (see de Lima and Atuncar [3]).

5. Application to a Real Call Center Data Set

We analyzed a real database of 7,761 phone calls made on April 7, 2006, from 8:00 am until 1:00 pm to a call center (data available from the authors upon request). We worked on a set of interarrival and service times to find the required minimum size of the server facility taking into account that the arrival rates are not homogeneous. In fact, larger interarrival times were detected toward the beginning of the day rather than throughout the day. This will lead to the adjustment of a different system for each hour of service that had an approximately homogeneous arrival rate. Using this data, we will model a 𝐺𝐼^𝑋/𝑀/𝑐/𝑁 system with a kernel interarrival density estimate for each hour of service.

The data set is presented in seconds, which is the precision given by the data acquisition system. This constraint will lead to ties in the arrival times (i.e., more than one call can arrive in the same second). Therefore, we treated same-time calls as part of a single arrival group and fit a discrete probability distribution. Table 1 shows the observed frequency of group sizes 𝑋 and the fitted distribution 𝐷 by hour. These arrivals were observed in each second.

Table 1 

Observed frequency of group sizes 𝑋 and fitted distributions 𝐷.

We used the gamma kennel-estimator with optimal 𝑏_AMISE as the kernel density estimation method. This method gave estimated densities bounded at 𝜏=0 and illustrates the information we miss when data is rounded. We also considered the gamma kernel estimator with LSCV method, but it did not behave as well with discrete data. Table 2 shows the smoothing parameter calculated, and Figure 3 shows the estimated densities.

Table 2 

AMISE smoothing parameter for gamma kernel.

(a) 8 am-9 am

(b) 9 am-10 am

(c) 10 am-11 am

(d) 11 am-12 am

(e) 12 am-1 pm

Figure 3 

Estimated interarrival times by hour.

Service time distribution fits an exponential distribution with parameter 𝜆=0.003339 and standard error of 3.8×10⁻⁵. The behavior of the system then becomes restricted to the choice of the number of servers 𝑐 and capacity of the system 𝑁. In Table 3, the minimum number of servers required is set to maintain the stability of the system at each hour considered.

Table 3 

Minimum number of servers.

To optimize system performance, it is necessary to establish some criterion function for the adopted design. If we only consider the effective arrival rate (𝜆_EF=𝜆[1−𝑃_BL]), we can describe this performance measure according to the number of servers and the maximum queue size. Figure 4 shows this relationship from 8 to 9 am. We saw similar results (not shown) for different periods. A good criterion would be the 𝐶=(𝐶₁/𝜆_EF+𝐶₂𝑐+𝐶₃𝑁), where 𝐶₁, 𝐶₂, and 𝐶₃ are costs related to each parameter and defined by the system environment. As an example, if we have 𝐶₁=500, 𝐶₂=1 and 𝐶₃=1, the state with minimal criterion would be 𝑐=19 and 𝑁=20. If cost 𝐶₁ is raised to 600, the best configuration would be 𝑐=21 and 𝑁=21.

Figure 4 

Effective arrival rate as a function of 𝑁 and 𝑐.

Figure 5 shows the invariant distributions estimated for other periods with different queue and server facility sizes. The method seems helpful for adjusting the maximum size 𝑁, and it is clear that the period plays a key role.

(a) 9 am-10 am

(b) 10 am-11 am

(c) 11 am-12 am

(d) 12 am-1 pm

Figure 5 

Estimated invariant distributions.

6. Conclusions

The main focus of this paper is on using kernels in the analysis of queuing systems. We studied the adequacy of the kernel estimator methods for calculating the invariant probability distribution and performance measures of queueing systems that have general interarrival distribution times with bulk arrivals. Discussions in simulations showed that the method is effective. Also, we successfully applied the method to the calls from a call center to plan and improve the performance of these important queueing systems.

Future research may take other directions—for instance, determining how the estimate of the prearrival invariant distribution moves away from its real value. One approach is from the variance of the bias and variance of each probability estimated. For example, the variance of an estimate of transition probability 𝑝_𝑖𝑗 could depend on the variance of term 𝛽̂₀, so we would need to findvar𝛽̂₀=var^∞0𝑒^{−𝑐𝜇𝑧}𝜏̂(𝑧,ℎ)𝑑𝑧.(6.1) This research can be applied to several areas of practical interest, including health and industry. For example, it could be critical in a medical emergency room that must optimize resource allocation.

Another possible future direction for research in this field involves the development of models to deal with dependency for the interarrival times. Such a case is particularly true in computer networks. Thus, in computer communications networks, people are interested in finding proper bounds of arrivals (see Li and Zhao [28]), with the assumption that arrivals are dependent (see Li and Zhao [29]), which are recently published results based on the so-called deterministic queuing theory (for details, see Le Boudec and Thiran [30]).

Funding

The Brazilian government funding agencies mentioned above had no role in the study.

Acknowledgments

This research is supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico; Grants no. 201046/1994-6, 301809/1996-8, 307702/2004-9, 472066/2004-8, 304944/2007-6, 561259/2008-9, 553019/2009-0, 550207/2010-4, 501532/2010-2, and 303388/2010-2), by CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior; Grant no. BEX-0522/07-4), by FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais, Grants nos. CEX-289/98, CEX-855/98, TEC-875/07, CEX-PPM-00401/08, and CEX-PPM-00390-10), and by PRPq-UFMG (Pró-Reitoria de Pesquisa da Universidade Federal de Minas Gerais).