Abstract
We deal with finite-buffer queueing systems fed by a Markovian point process. This class includes the queues of type M/G/1/N, /G/1/N, PH/G/1/N, MMPP/G/1/N, MAP/G/1/N, and BMAP/G/1/N and is commonly used in the performance evaluation of network traffic buffering processes. Typically, such queueing systems are studied in the stationary regime using matrix-analytic methods connected with M/G/1-type Markov processes. Herein, another method for finding transient and stationary characteristics of these queues is presented. The approach is based on finding a closed-form formula for the Laplace transform of the time-dependent performance measure of interest. The method can be used for finding all basic characteristics like queue size distribution, workload distribution, loss ratio, time to buffer overflow, and so forth. To demonstrate this, several examples for different combinations of arrival processes and characteristics are presented. In addition, the most complex results are illustrated via numerical calculations based on an IP traffic parameterization.
1. Introduction
Since the beginning of the 1990s, when the strong auto-correlation of the Internet traffic was discovered, a variety of processes have been developed or adapted for proper teletraffic modeling. For instance, fractional Brownian motion [1], chaotic maps [2], FARIMA [3], and multifractal wavelets [4] have been applied in wide range of tasks connected with performance evaluation of buffering processes, traffic predictability, congestion and admission control, buffer sizing, and so forth.
However, none of the aforementioned processes suits as well for the teletraffic modeling as the famous class of Markovian point processes [5] or one of its well-known reparameterizations or subclasses (MMPP, MAP, BMAP, etc.). First of all, this is connected with the fact that processes are analytically tractable. They are also easy to simulate, and a variety of parameter fitting procedures have been developed for them [6β12].
The only disadvantage of Markovian models used for teletraffic modeling is that they are not truly self-similar or long-range dependent. However, for practical purposes, it is typically enough to mimic the self-similarity over a few time scales. This can be easily accomplished using Markovian processes and, as shown in [9], the resulting model can be reliable in terms of its marginal distribution, autocovariance function, and queueing behaviour.
One of the main reasons for developing traffic models is finding their queueing performance characteristics. In this paper, we deal with the finite-buffer queue whose arrival process is given by a Markovian point process from the class . So far such queueing systems have been solved typically in their stationary regime using matrix-analytic methods connected with M/G/1-type Markov chains [13β16] (this set of papers is not intended to be exhaustive, the literature devoted to the subject is vast).
Herein, a different, unified method for solving these queueing systems is described. Its main advantage is that it gives formulas for the characteristics of interest in a closed, easy-to-use form. It is devoted to computing transient characteristics, but the steady-state measures could also be obtained from the transient solutions. It can be used for all processes in the class (e.g., Poisson processes, batch Poisson processes, phase-type renewal processes, MMPPs, MAPs, BMAPs) and for many queueing performance characteristics, including queue size, virtual waiting time, loss ratio, time to buffer overflow, buffer overflow period. To demonstrate this, three queueing systems with different combinations of arrival processes and characteristics of interest are solved using the proposed method. In every next queueing system, an arrival process of growing complexity is used. Namely, we will start with the Poisson arrival process and the queue size distribution and finish with the MAP arrival process and the workload distribution.
The remaining part of the paper is structured in the following way. In Section 2, a description of the proposed method is presented. In Sections 3, 4, and 5, three detailed examples of its applications are shown. In particular, in Section 3 a formula for the transient queue size distribution in the M/G/1/N system is proven and illustrated via numerical examples. In Section 4, a formula for the time to buffer overflow in the /G/1/N queue is presented. In Section 5, a formula for the workload distribution in the MAP/G/1/N model is shown and illustrated via numerical example based on an IP traffic parameterization. Finally, remarks concluding the paper are gathered in Section 6.
2. Method
The proposed method can be sketched in the following three-step scheme.(I)In the beginning, we apply the total probability formula with respect to the first departure time. This allows us to utilize the Markovian structure of the arrival process and develop a system of integral equations for the characteristic of interest.(II)Then, by using the Laplace transform technique, we reduce the problem to a system of linear equations.(III)In the next step the solution of the resulting system of equations is presented in a closed-form formula using recurrent sequences.
By means of the resulting formula, we can compute the steady-state characteristic at once using basic properties of the Laplace transform or we can compute the transient characteristic applying an inversion algorithm.
The third step in this scheme is based on the following lemma (for proof, see [17, page 201]).
Lemma 2.1. Assume that , , is a sequence of matrices such that is nonsingular and , is a sequence of column vectors of size . Then every solution of the system of equations has the form: where is a column vector that does not depend on and the sequence is defined to be and denotes the matrix of zeroes.
It is easy to check that if the system (2.1) is indexed from 0, namely then its every solution has the form
Now it is time to show how this method works in practice.
3. Poisson Arrivals and Queue Size Distribution
In the first example, we will find a formula for the transient queue length distribution in the M/G/1/N model, that is, for the system with Poisson arrivals (with intensity ), general type of the service time distribution (given by distribution function ), and finite capacity (the total number of customers in the system, including service position, must not exceed ).
To the best of our knowledge, a closed-form formula for the transient queue size in the M/G/1/N system has not been reported in the English literature yet. A transient solution for the infinite-buffer M/G/1 system can be found in [18, Section 1.7] and [19, Chapter 3].
The transient behavior of queueing systems depends on the initial buffer content. Herein the initial buffer occupancy is not further specified and can be zero or nonzero. It is assumed that the time origin corresponds to a departure epoch. Thus, if the initial buffer content is non-zero, then the service begins at the time origin. Otherwise, the service begins at the first arrival time. The service time distribution can have any particular form, but for practical reasons we restrict this study to the class of service time distributions with explicit Laplace-Stieltjes transform.
Let denote probability, the queue size process, and
(I) We start from using the total probability formula with respect to the first departure time. For , we obtain The first sum in (3.2) describes the situation where the first departure time is before and there is no buffer overflow by the time , which means that the number of arrivals in must be less than . The second sum describes the situation where the first departure time is before and an overflow occurs by the time . Finally, describes the case where the first departure time is after .
Using the total probability formula with respect to the first arrival time for the initially empty system (), we have where is the Kronecker symbol, that is, if and 0 otherwise.
(II) In the second step, we apply the Laplace transform to both sides of (3.2) and (3.3). Therefore, for the transform we obtain where Substituting , we get from (3.5), (3.6) with
(III) Now, the system (3.8) has exactly the same form as (2.4). Thus, its solution has the same form as (2.5), namely, where does not depend on and is given in (2.3) for the sequence .
Now we only need to find unknown , , and . In order to find , we put into (3.11) and get . Putting into (3.8) and observing that yield . Then, substituting these results into (3.11) we have
To calculate , we set and in (3.12) and use the boundary condition (3.9). This gives Using (3.16) and (3.12) with , we obtain the final result.
Theorem 3.1. The transform of the queue size distribution in the M/G/1/N system has the form where and are given in (3.13) and (3.14), respectively.
Now we can obtain some numerical results.
3.1. Numerical Example
In this example, we will observe the transient queue size distributions and check how long it takes to stabilize the initially overflowed queue.
We assume that the system capacity is 20 (i.e., ), the arrival rate is 1 (i.e., ), and the service time is constant and equal to 0.9. Therefore, we have βthe traffic intensity is moderate. However, we assume that the system is initially full (i.e., ).
Using Theorem 3.1 and the Laplace transform inversion proposed in [20], we obtain the results depicted in Figure 1 and Table 1. In Figure 1 we can observe the queue size distribution after 10, 20, 50, 100, and 200 seconds of the system work and in the steady state (the thick curve). The distribution converges from shapes concentrated around 20 to the steady-state distribution. Shapes close to the steady-state are achieved after about 200βs of the system work.
As we can see, for high values of , the distribution of the queue size reaches its maximum at . To explain this, we note first that . This causes that for growing the probability mass moves towards the level 0. However, the maximum is not at the level 0, which is connected with the fact that the level 0 is a reflecting barrier for the queue size process. Roughly speaking, the level 1 can be reached either from the level 2 (job departure) or from the level 0 (arrival of job to an empty system), but the level 0 can be reached from the level 1 only (job departure). Therefore, the probability at is higher than at .
In Table 1 we can observe the convergence of the standard deviation to its steady-state value. We may notice that the standard deviation does not change monotonically and reaches a maximum for some . This can be explained in the following way. As we start from a queue size of 20, for a small the probability mass is concentrated around 20. Moreover, as we also have , the queue cannot get longer than 20. Therefore, for a small , the distribution has only the left tail, and its variance is small. Now, for a large , as explained before, the probability mass is concentrated around 1, the distribution has only the right tail and its variance is also relatively small. On the other hand, for moderate values of , the probability mass is distributed more uniformly between 0 and 20, which results in a higher variance. Thus, at least one maximum is to be expected for moderate values of .
4. Batch Poisson Arrivals and Time to Buffer Overflow
In the second example, we will find a formula for the distribution of the time to buffer overflow in the /G/1/N system. In this model, groups of customers arrive according to the Poisson process with rate . Sizes of consecutive groups are independent, identically distributed with discrete distribution , where . The partial rejection scheme is assumed. This means that, in the case of insufficient remaining buffer capacity for all the customers included in an arriving group, only a part of it is accepted and the rest is lost. We assume again that the service time of one customer is distributed according to distribution function , which is not further specified.
We are interested in the distribution of the time to buffer overflow in this system, namely, in the distribution of , where is defined as follows: and denotes the number of customers in the system at time .
(I) Using the total probability formula with respect to the first departure epoch for initially nonempty system, , we have where denotes the probability that customers arrive in interval .
The first term in (4.2) describes the situation where the first departure time is before and there is no buffer overflow by the time , which means that the number of arrivals in must not exceed . The second term describes the situation where the first departure time is after and there is no buffer overflow by the time . Naturally, the situation where an overflow occurs in interval is not taken into account now, as in this case we have .
If the system is initially empty, then conditioning on the first arrival epoch we get
(II) The Laplace transform applied to (4.2) and (4.3) reduces the problem to with Substituting we obtain where .
(III) Now, applying Lemma 2.1 the general solution of the system (4.6) has the form where does not depend on and is given in (2.3) for defined in (4.5).
Putting in (4.8), we can observe that . Then, using condition (4.7) together with (4.8), we have Finally, rewriting (4.8) as we arrive at the following theorem.
Theorem 4.1. The transform of the time to buffer overflow in the /G/1/N system has the form where and are given in (4.9) and (4.10), respectively.
To make this theorem useful, we have to be able to compute and . Computing these coefficients is not very demanding and may be carried out, for instance, using generating functions. It is easy to check that A very effective algorithm for generating function inversion can be found in [20]. Namely, if we have a generating function , then the original values of can be restored as where while and are used to control the roundoff error. (Typically, we use , .)
An alternative way of computing and is the uniformization technique [21]. Applying this technique to , we obtain where and can be computed as follows: Similarly, for , we get with
For the bibliography on other computational methods for , see [22].
4.1. Numerical Example
To demonstrate how (4.12) can be used in practice, let us assume that we have batch Poisson arrivals parameterized as follows: Therefore, the average batch size is equal to 4 and the total arrival rate is 2. We assume that the service time is constant and equal to (which gives ) and the buffer size is 100.
Suppose we want to compute the average time to buffer overflow, starting from an empty buffer, that is, It is easy to see that
The value of can be computed using the uniformization technique. From (4.17) and (4.20), we obtain, respectively, where denotes the incomplete gamma function. Now, using (4.16) and (4.19) we can compute and . Finally, applying (4.12), we get As we have a moderate traffic intensity and a quite big buffer, a large time to buffer overflow was to be expected.
5. MAP Arrivals and Workload Distribution
In the third example, we will find a formula for the workload distribution in the MAP/G/1/N model, that is, for the model with MAP arrivals, general type of the service time distribution (given by distribution function ), and finite capacity .
The Markovian arrival process (MAP) is one of the most flexible arrival processes from the class of Markovian processes. It enables a very precise fitting to network trace files in terms of not only the basic statistical parameters (mean, variance, higher moments) but also the shape of the marginal distribution and autocorrelation function. (For the newest, excellent parameter fitting procedures for MAP processes, see [11].)
The MAP is parametrized by two matrices, and , such that is nonnegative, has nonnegative off-diagonal elements, and negative diagonal elements and is an irreducible infinitesimal generator. We will use to denote the state of the underlying Markov chain, to denote the number of arrivals in , and to denote the counting function, that is, We will also use intensities and probabilities , , defined as
By the workload we mean the length of time a job (packet) which arrives at time waits before entering service. This is one of the most important characteristics from the practical point of view as it can be used to compute the queueing delay for packets or cells in network devices. The workload in a MAP queue has been studied so far either in the infinite-buffer model [23] or in the steady state [24]. We assume herein that the workload of a blocked cell is zero.
We will study the workload using its Laplace transform in the column vector form
As previously, denotes the number of customers in the system at time .
(I) As in the previous sections, we start from using the total probability formula with respect to the first departure moment. For , , we obtain where is the -fold convolution of the distribution function with itself.
The first double sum in (5.5) describes the case where the first departure time is before and there is no buffer overflow by the time . The second double sum describes the case where the first departure time is before and an overflow occurs by the time , which means that the number of arrivals is equal to or more. The third double sum describes the case where the first departure time is after and there is no overflow by the time . In this case, we have , where is the number of arrivals in .
If the system is initially empty, then for we get
(II) Applying transforms and matrix notation to (5.5) and (5.6) we obtain where Replacing and (5.7) gives
(III) We can see now that (5.10) has the same form as (2.4). Therefore, its solution is given in (2.5). Proceeding in the same way as in the previous sections, we arrive at the final result.
Theorem 5.1. The transform of the workload distribution in the MAP/G/1/N system has the form with
Note that matrices , , and can be computed effectively by means of the uniformization technique [21]. Using the elementary properties of the Laplace transform we can easily obtain the average workload in steady stateβsimply by calculating . Putting into and inverting the result with respect to only, we can compute the transient average workload. Finally, using a two-dimensional inversion algorithm, we may obtain the shape of the workload distribution for an arbitrary .
It is easy to check that the number of floating-point operations needed to compute (5.11) (time complexity) grows as . This estimate is a consequence of the form of (5.11) and (2.3) and the fact that matrix multiplication and inversion are of order. Thus, the approach proposed herein reduces the numerical complexity when comparing it to the brute-force solution of the system (5.7), which is of order.
5.1. Numerical Example for MAP Arrivals
For numerical purposes, we are going to utilize a parameterization of the MAP based on a recorded IP traffic sample. To accomplish that, the AMP-1138809025-1.tsh trace file, recorded at the AMP aggregation point run by the Passive Measurement and Analysis Project, has been used. Using an implementation of the EM algorithm [7] written for Mathematica environment, the following MAP parameterization was obtained: The average rate of the fitted MAP is where is the stationary distribution for the underlying Markov chain and .
It is assumed that the service time is constant and equal to . Manipulating we can easily obtain different traffic intensities .
In Figure 2, the mean workload as a function of time, , for the initial buffer occupancy of 0%, 25%, 50%, 75%, and 100% is depicted. The traffic intensity was set for , the initial phase for , and the buffer size for 100 packets.
As we can see, no matter what the initial buffer occupancy was, the steady-state value (1.549βms) was reached after about 0.5βs, which is equivalent to about 3800 packet arrivals.
In Figure 3, the stationary mean workload, , as a function of the buffer size, is shown for four traffic intensities, namely, 0.75, 0.90, 0.95, and 0.99. In each case, the curve becomes flat starting from some threshold value of the buffer size. For this border buffer size is about 20, for about 50, for about 100, and for about 500.
(a)
(b)
(c)
(d)
There is an obvious explanation of this behaviour of the workloadβfor large buffers the finite-buffer system is practically equivalent to the infinite-buffer one; thus, the constant workload observed for large buffers is equal to the infinite-buffer value. However, this behaviour is of some practical importance, especially when the border buffer size is known. Decreasing the buffer size below this border value, we can shorten the queueing delay of the system. The cost paid for this is a higher loss ratio, but it can be beneficial in some applications. In order to evaluate the tradeoff precisely, we have to know the loss ratio, which also can be computed using the method presented in this paper.
6. Conclusions
We presented a unified method for solving queues with Markovian arrivals. The most important features of this approach are the following(i)it can be applied for finding both steady-state and transient characteristics;(ii)it produces results in a closed, easy-to-use form;(iii)reduced numerical complexity comparing to the brute-force solution;(iv)it is suitable for computing many performance measures of finite-buffer queues, including queue size, workload, loss ratio, time to buffer overflow, buffer overflow period.
The main disadvantage of the method is that it cannot be used directly in solving infinite-buffer queues. This is connected with the necessity to invert the order of the system of equations (for instance, the substitution in (3.5)), which cannot be carried out in the infinite-buffer model.