Queues with Dropping Functions and Autocorrelated Arrivals

We present an analysis of the queueing system in which arriving jobs are dropped with probability depending on the queue size. The arrivals are assumed to be autocorrelated and they are modeled by the Markov-modulated Poisson process. Both transient and stationary distributions of the queue size, as well as the system loss ratio and throughput are obtained. The analytical results are accompanied with numerical examples based on the autocorrelated traffic recorded in an IP computer network.


Introduction
The queue with the dropping function is a simple FIFO queue with an additional mechanism. Namely, an arriving job (customer, packet etc.) can be dropped (rejected) with probability d(n), where n is the queue size observed upon arrival of this job (see Fig. 1). The function d(n) is called a dropping function.
The significance of the queueing system with the dropping function can be presented from two perspectives: (a) universal sense of the system and its general applicability, (b) direct applicability of the system in networking. As regards (a), it is well known that in the classic FIFO queueing model we cannot control the performance. Given the arrival and service processes, we can calculate characteristics like the queue length, system throughput, number of losses, etc., but we cannot control them. On the other hand, in some applications of queueing systems there is a need to control the performance of the queue, namely to set the mean queue size or the throughput etc. to an arbitrary value.
There are at least three ways to achieve that. Firstly, we may try to alter dynamically arrival rate, depending on the current or past system state. Secondly, we may try to alter dynamically the service rate. Thirdly, we may try to reject arriving jobs. The latter is similar to alternating the arrival rate, but the significant difference is that dropping jobs causes losses, i.e. jobs that were not served and never return to the queue.
There are several queueing models of the first and the second type analyzed in the literature. For example, threshold-based queueing models are studied in (Chydzinski 2002;Pacheco and Ribeiro 2008;Bekker 2009). In such models, the arrival or service rates alternate when the queue length reaches a threshold level.
Intuitively, the third approach is the simplest one -it is usually a simple matter to drop an arriving job. Furthermore, application of the dropping function enables a powerful control on the performance of the queueing system. For instance, in Chydzinski and Chrost (2011) it was shown that using dropping functions allows setting the average queue size. The same applies to other parameters, e.g. the queue size variance, the system throughput, etc.
Queueing systems with dropping functions have at least one direct application, namely the active queue management in Internet routers. It has been shown that simple FIFO tail-drop queues, commonly used in Internet routers by device vendors, have some important disadvantages. In particular, they cause large queueing delays, flow synchronization and unfairness between flows. To overcome these problems, the active queue management (AQM) for Internet routers was proposed. The idea was that the router can drop incoming packets even if the buffer is not full yet, thus preventing the queue from building up.
The router can drop incoming packets depending on several factors, but the simplest approach is that the incoming packet is dropped with the probability that is a function of the queue size. Now, it must be stressed that the analytical results obtained so far for queues with dropping functions fed by classic queueing traffic models (e.g. Poisson) are of little use when the real arrival process is autocorrelated. This is the case of networking -it is well known that Internet traffic possesses strongly autocorrelated structure. In Fig. 2  interarrival times autocorrelation function based on recorded traffic is presented 1 . As we can see in the figure, the autocorrelation decays slowly and is significant on several time scales.
Not taking this autocorrelation into account may lead to extreme underestimation (in a negative way) of characteristics of the queueing system. An example can be found in Fig. 2 of Chydzinski (2006a), where a comparison of computed full buffer probabilities in FIFO queues is presented. Both queues have the same buffers, arrival rates and service rates, but one of them uses uncorrelated arrivals (Poisson), whilst the other uses autocorrelated arrivals with the autocorrelation function as in Fig. 2. The obtained full buffer probability for the Poisson traffic is 8.97 × 10 −33 , while the value obtained for the autocorrelated traffic is 6.55×10 −3 . Therefore, not taking into account the autocorrelation resulted in the optimistic underestimation by 30 orders of magnitude.
For these reasons, in this paper we carry out an analysis of the queueing system with the dropping function and autocorrelated arrivals. We do not impose any assumptions on the dropping function nor the service time distribution -they both can have any form.
In order to make the results useful in practice, the following two requirements on the arrival process must be met: it has to be able to mimic arbitrary shape of the autocorrelation function and there should exist an algorithm, which allows fitting the model parameters to this particular shape.
Both these requirements are fulfilled by the Markov-modulated Poisson process (MMPP), Fischer and Meier-Hellstern (1992). In particular, the MMPP is analytically tractable, allows for precise fitting of complicated shapes of the autocorrelation function (see e.g. Salvador et al. (2003)) and several MMPP parameter fitting procedures have been proposed to date, (Salvador et al. 2003;Deng and Mark 1993;Ryden 1996;Yoshihara et al. 2001;Singh and Dattatreya 2004). Therefore, it will be used as the arrival traffic model.
The paper is organized as follows. In Section 2, references to papers on queues with dropping functions, the MMPP process, the MMPP queue and the active queue management are given. In Section 3, the definition and basic formulas on the MMPP arrival process are recalled. In Section 4, the queueing model is formally defined. Then, in Section 5, the results on transient and stationary queue size distributions, as well as the loss ratio, are presented. Section 6 is devoted to calculation of the counting function of the arrival process filtered by the dropping function -this is needed to use the results of Section 5 in practice. In Section 7, some numerical examples are presented. In particular, the examples demonstrating the abilities of the dropping function to maintain a given average queue size and a given throughput are shown. The final conclusions are gathered in Section 8.

Related Work
Queueing systems with dropping functions and uncorrelated arrivals were studied in the following papers. In Bonald et al. (2000), an approximate analysis of the model with batch Poisson arrivals, linear dropping function and exponential service time distribution was presented. In Hao and Wei (2005), an approximate analysis of the model with general batch arrivals and exponential service times was carried out. In Chydzinski and Chrost (2011) an exact, steady-state analysis of the model with arbitrary dropping function, Poisson arrivals and arbitrary service time was performed. In Kempa (2013a) an exact analysis of the model with arbitrary dropping function, general uncorrelated arrivals and exponential service time was presented. In Kempa (2013b), the transient analysis of the model with arbitrary dropping function and Poisson arrivals was carried out. Finally, in Tikhonenko and Kempa (2013) the system with the Poisson arrival stream and a general distribution of the job size has been solved.
As for the arrival process, we refer the reader to the excellent MMPP cookbook, Fischer and Meier-Hellstern (1992), and the references given there. In this cookbook, the main results on the infinite-buffer MMPP queue (without the dropping function), are presented as well. The finite-buffer MMPP queue in the steady state was analyzed in Baiocchi and Blefari-Melazzi (1992). The same system in the transient state was studied in Chydzinski (2006a).
Regarding the active queue management in routers, the famous RED algorithm, Floyd and Jacobson (1993), was the first exploiting the dropping function. In the case of RED, the dropping function is the simplest possible, i.e. the linear function. Besides the linear function, some other dropping functions were used in literature: the exponential dropping function, Athuraliya et al. (2001), or the doubly linear dropping function, Rosolen et al. (1999). The active queue management is a widely studied subject. In addition to such well-recognized algorithms as REM (Athuraliya et al. 2001), PI (Hollot et al. 2002;Wu et al. 2001), BLUE (Feng et al. 2002), GREEN (Wydrowski and Zukerman 2002) and AVQ (Kunniyur and Srikant 2004), several other propositions emerged recently, e.g. (Na et al. 2012;Suzer et al. 2012;Farzaneh et al. 2013;Kahe et al. 2013). The new algorithms are usually evaluated either by means of simulators (ns2, ns3), or by means of the control theory. We prefer a different approach, based on tools and results of queueing theory.
The methodology used herein is an extension of the methodology used in solving classic queues with Markovian arrival processes (see e.g. Chydzinski 2006aChydzinski , 2006bChydzinski , 2006c, based on formulating and solving a set of Volterra integral equations in the convolution form. The main difference and difficulty herein is a complicated characterization of the counting function of the arrival process filtered by the dropping function, which requires a different approach.

Arrival Process
The Markov-modulated Poisson process, Fischer and Meier-Hellstern (1992), is constructed by varying the rate of Poisson arrivals according to a continuous-time Markov chain, called the modulating chain. Assuming that the state space of the modulating chain is {1, . . . , m}, to parameterize an MMPP we have to provide the intensity matrix of the modulating chain, Q, and the vector of corresponding intensities of Poisson arrivals, [λ 1 , . . . , λ m ]. The latter will also be used in the form of a diagonal m × m matrix: The average rate of the MMPP can be obtained as where π is the stationary vector of the modulating chain, i.e.: − → 1 denotes the column vector of 1's. The MMPP counting function is defined as: where P denotes probability, N(t) denotes the number of events (jobs, customers, packets) in time interval (0, t], while J (t) denotes the state of the modulating chain at time t.
In this paper we pay a special attention to the autocorrelated structure of the MMPP. The k−lag autocorrelation of the MMPP equals:

Queueing Model
The system of interest is a single-server queue, in which the customers form the Markovmodulated Poisson process and they are served in the arrival order. A general type of service time distribution, given by a distribution function F (t), is assumed. The queue size is limited by a finite buffer. Namely, the total number of jobs present in the system is always smaller or equal to b. A job that arrives when the buffer is full (i.e. there are b jobs present in the system) is dropped and never returns.
In addition to that, any arriving job can be dropped. This happens with probability d(n), where n is the queue size at the time of arrival of this job (including the service position).
The function d(n) is called the dropping function. It may take any values in [0, 1] for n = 0, . . . , b − 1. The finite-buffer assumption forces that d(n) = 1 for n ≥ b.
The queue size at time t will be denoted by X(t). We adopt the convention that X(t) includes the service position, if occupied. If X(0) > 0, then it is assumed that the time origin corresponds to a service completion.
The load offered to the queue is then:

Queue Size and Loss Ratio
Let us denote by n,i (t, l) the probability that the queue size at time t equals l, provided that the initial queue size was n and the initial state of the modulating chain was i, i.e.: In order to compute n,i (t, l), we have to use function A n,k,i,j (u), which is the counting function of the arrival process filtered by the dropping mechanism. Namely, A n,k,i,j (u) is defined as the probability that in a system without service (arrivals only), exactly k jobs are accepted to the queue in time interval (0, u] and at the end of this interval the state of the modulating chain is j , provided that the queue size at t = 0 was n and the state of the modulating chain at t = 0 was i. Assuming that the system is not empty at the time origin and using the total probability formula with respect to the first departure time, u, we obtain for 1 where The first summand of Eq. 1 corresponds to the situation, where the first service completion time, u, occurs before t, while the second summand to the situation where there is no service completion by the time t.
Assuming that the system is empty at the time origin and using the total probability formula with respect to the first event in the arrival process, which may be a change of the modulating state or a job arrival, we get for 1 ≤ i ≤ m: where ij , Q ij denote the entries in the i-th row and j -th column of , Q, respectively, while δ ij = 1 if i = j and 0 otherwise.
The first summand of Eq. 2 corresponds to the case, where the first event happens before t and it is a change of the modulating state. The second summand corresponds to the case, where the first event happens before t, it is a job arrival, but the job is dropped. The third summand corresponds to the case, where the first event happens before t, it is a job arrival and it is accepted. Finally, the fourth summand corresponds to the case, where there are no events in the MMPP by the time t.
Introducing the following notation: and applying the Laplace transform to Eqs. 1 and 2 we get for 1 ≤ n ≤ b, 1 ≤ i ≤ m: and for 1 ≤ i ≤ m: respectively. Using the following m × m matrices: we can rewrite Eqs. 3 and 4 as φ n (s, l) = b−n k=0Ã n,k (s)φ n+k−1 (s, l) + r n (s, l), respectively. Then we may group all known coefficients of Eqs. 13 and 14 into one in the form: where I is the unit matrix of size m × m. Now the system Eqs. 13 and 14 can be rewritten as where R(s, l), φ(s, l) are the following column vectors of size (b + 1)m: In order to compute φ(s, l) from Eq. 16, matrix G(s) must have an inverse. As it can be hard to prove that G(s) is non-singular in general, we assume from now on that this is the case. This is not a problem in practice -we never came accross a singular G(s), despite performing a large number of numerical computations for different system parameters. Now Eq. 16 can be solved and we obtain the following result.

Theorem 1 The Laplace transform of the queue size distribution at time t of a finite-buffer queue with the dropping function has the form:
where G(s) and R(s, l) are given in formulas Eqs. 15 and 18 respectively.
The practical applicability of Theorem 1 depends on our ability to compute function A n,k,i,j (u) -all other components are simple functions of the parameters of the queueing system.
Therefore the whole next section will be devoted to the computation of A n,k,i,j (u). Now let us observe only, that having A n,k,i,j (u), we can compute several important characteristics of the queueing system using Theorem 1. Firstly, the stationary distribution of the queue size can be computed. Namely, denoting P l = lim t→∞ P(X(t) = l), and using the well-known properties of the Laplace transform we have where [·] 1 denotes the first element of a vector. In fact, any other element could be used in Eq. 20. Secondly, we can obtain the transient distribution of the queue size, i.e. the distribution at an arbitrary chosen time t, which requires numerical inversion of the Laplace transform (19).
Thirdly, using Eq. 20 we can compute the loss ratio and throughput. The loss ratio, L, is defined as the long-run fraction of jobs which were dropped upon arrival. In a time interval of length T , where T is large, the number of finished jobs is approximately equal to: In the same interval, there are approximately λT new arrivals. Therefore, the fraction of accepted jobs (the throughput) must be: while the loss ratio must be:

Calculating A
In this section we deal with finding function A n,k,i,j (u), which is a crucial task for practical applicability of Theorem 1. Firstly, let us consider the case, where in time interval (0, u] at least one job is accepted to the queue, i.e. where k ≥ 1 in A n,k,i,j (u). The first event in time interval (0, u] may be a change of the modulating state, or an arrival of a job which is immediately dropped, or an arrival of a job which gets accepted. Therefore, for any k ≥ 1, 0 ≤ n ≤ b, 1 ≤ i ≤ m, 1 ≤ j ≤ m we have: Secondly, assuming that there are no jobs accepted in (0, u], the first event in (0, u] may only be a change of the modulating state, or an arrival of a job which is immediately dropped. Alternatively, there may be no events in (0, u] at all. Summarizing, for any 0 ≤ n ≤ b, 1 ≤ i ≤ m, 1 ≤ j ≤ m we obtain: Applying the Laplace transform to Eqs. 23 and 24 we get for k ≥ 1: (25) and a n, where a n,k,i, Using the following m × m matrices: A n,k (s) = a n,k,i,j (s) i,j , for k ≥ 1 (25) yields: while Eq. 26 yields: A n,0 (s) = U n (s)A n,0 (s) + Z(s).
Hereafter we assume that matrix I − U n (s) is non-singular. As in the case of G(s), we verified this assumption in a large number of numerical calculations for different system parameterizations. Denoting and E n (s) = (I − U n (s)) −1 Z(s), we have for k ≥ 1: and A n,0 (s) = E n (s).
Finally, exploiting Eqs. 31, 32 and the mathematical induction we arrive at the following theorem.

Theorem 2 The Laplace transform of the counting function of the MMPP filtered by the dropping function is
where C n (s) and E n (s) are given in Eqs. 29 and 30, respectively.
It can be easily verified, that this is a generalization of the result for the ordinary Poisson process presented in Theorem 1 of Chydzinski and Chrost (2011). As for the numerical calculations of A n,k,i,j (u) values, we use the Spinelli inversion method (Spinelli 1966).
This parameterization was obtained in Chydzinski (2006a) as a result of fitting the MMPP to the recorded network traffic. The autocorrelation function of this MMPP is shown in Fig. 3, which is to be compared with Fig. 2, depicting the autocorrelation of the original traffic. The average rate of this MMPP is λ = 71729.36.
In all examples the buffer size is b = 200. The service time is constant and chosen in such a way that the queueing system is mildly overloaded, i.e. ρ = 1.1.
As we can see, the distributions are quite different and the average queue size can vary by an order of magnitude, depending on the shape of the dropping function, even though the throughput is kept at the level of 80 %. This indicate that we may search for dropping functions that provide both the required throughput and the required average queue size, e.g. γ = 80 % and the average queue = 30.
Finding dropping functions that provide only a predefined average queue size is very easy. For instance, the following three dropping functions assure the average queue size of 75:  The shapes of functions d 5 , d 6 and d 7 are presented in Fig. 6, while their main performance characteristics in Table 2. Now the values of the performance characteristics are close. This is connected with the fact that functions d 5 -d 7 are all monotonic and positive for n > 50. If we have dropping function d that provides the average queue size a, it is easy to obtain any average queue size smaller than a by scaling the function d, i.e. using: where c > 1. For instance, using function d 6 we can obtain functions d 8 and d 9 which give the average queue of 50 and 25, respectively. Namely, we have: Fig. 7. The performance of these functions is summarized in Table 3.
By means of the dropping function we may also force particular values of performance characteristics for different load values. For instance, assume that we want the average queue size to be 60 for an overloaded system (ρ = 1.1) and 20 for an underloaded system (ρ = 0.9).  Again, manipulating the shape of the dropping function we may achieve this goal. For instance, the following function meets the presented requirements: Function d 10 is depicted in Fig. 8, while its performance is summarized in Table 4. It can be interesting to see, in one figure, the trade-off between the average queue size and throughput for different dropping functions. For this purpose, a scatter plot is presented in Fig. 9 for functions d 0 -d 10 , where d 0 is the zero dropping function, i.e.: Of course, for d 0 we obtain the maximum value of the average queue size, 119.4, and the maximum value of the throughput, 89.5 %. As we can see, among the considered functions, d 2 provides the smallest average queue size, while maintaining a moderate throughput (only 9.5 % worse than the best possible). On the other hand, d 5 provides a very high throughput (only 1.4 % worse than the best possible) while maintaining a moderate average queue size. It is also interesting to draw a scatter plot for several dropping functions belonging to one family. For instance, we studied a family of linear functions in the form: where a is a parameter. The following values of a were used in computations: 0.0005, 0.001, 0.002, 0.003, 0.004, 0.005, 0.007, 0.01, for k = 11, 12, . . . , 18, respectively. Therefore 8 new dropping functions, d 11 -d 18 , were obtained. They are depicted in Fig. 10. As we can see in Fig. 11, the steeper the function, the smaller the queue size, but the smaller the throughput at the same time. Therefore, we cannot say that one of the considered functions is "better" than other.
In this way, a very interesting question arises: what is the border of the set of all possible points in such scatter plots? In other words, what is the maximal possible throughput, given the average queue size of x? What is the minimal possible average queue size, given the throughput of y? So far, such questions seem to be hard to answer. Of course, using the results for d 0 we know that we can obtain any average queue size in interval [0, 119.4] and any throughput in interval [0, 89.5 %]. But it is hard to say for instance, if the average queue of 20.0 and the throughput of 89.0 % can be obtained at the same time.

Conclusions
We presented an analysis of the queueing system with the dropping function and autocorrelated interarrival times. We argued that taking into account the autocorrelated structure of the arrival process is crucial for applicability of the model in some important areas of applications, e.g. in networking. The main results of the analysis are the transient and stationary distributions of the queue size and the stationary loss ratio and throughput in the model with the MMPP arrivals. We also presented several numerical results exploiting 18 different dropping functions, the buffer for 200 packets and a model of autocorrelated traffic recorded in an IP network. The main purpose of these examples was to demonstrate the practical usability of the model for solving systems with realistic parameterizations (properly parameterized traffic model, reasonable buffer size and load etc.).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.