# A Unified Method of Analysis for Queues with Markovian Arrivals

Authors: Andrzej Chydzinski
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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,MX/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.


---

## Body

## 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 classN 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 N 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 classN. 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 classN (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 Section2, 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 MX/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 thatA0,  A1,  A2,… is a sequence of  m×m matrices such that  A0 is nonsingular and  ψ1,  ψ2,… is a sequence of column vectors of size  m. Then every solution of the system of equations
(2.1)∑k=-1n-1Ak+1xn-k-xn=ψn,n≥1,
has the form:
(2.2)xn=Rnc+∑k=1nRn-kψk,n≥1,
where c is a column vector that does not depend on n and the sequence Rk is defined to be
(2.3)R0=0,R1=A0-1,Rk+1=R1(Rk-∑i=0kAi+1Rk-i),k≥1,
and 0 denotes the m×m matrix of zeroes.It is easy to check that if the system (2.1) is indexed from 0, namely
(2.4)∑k=-1nAk+1xn-k-xn=ψn,n≥0,
then its every solution has the form
(2.5)xn=Rn+1c+∑k=0nRn-kψk,n≥0.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 F(t)), and finite capacity (the total number of customers in the system, including service position, must not exceed N).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.LetP(·) denote probability, X(t) the queue size process, and
(3.1)ϕ̂nl(t)=P(X(t)=l∣X(0)=n).(I) We start from using the total probability formula with respect to the first departure time. For0<n≤N, we obtain
(3.2)ϕ̂nl(t)=∑k=0N-n-1∫0tϕ̂n+k-1,l(t-u)e-λu(λu)kk!dF(u)+∑k=N-n∞∫0tϕ̂N-1,l(t-u)e-λu(λu)kk!dF(u)+ρnl(t),ρnl(t)=(1-F(t))⋅{0ifl<n,e-λt(λt)l-n(l-n)!ifn≤l<N,∑k=N-n∞e-λt(λt)kk!ifl=N.
The first sum in (3.2) describes the situation where the first departure time u is before t and there is no buffer overflow by the time u, which means that the number of arrivals in (0,u] must be less than N-n. The second sum describes the situation where the first departure time u is before t and an overflow occurs by the time u. Finally, ρnl(t) describes the case where the first departure time u is after t.Using the total probability formula with respect to the first arrival time for the initially empty system (X(0)=0), we have
(3.3)ϕ̂0l(t)=∫0tϕ̂1l(t-u)λe-λudu+δ0le-λt,
where δij is the Kronecker symbol, that is, δij=1 if i=j 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
(3.4)ϕnl(s)=∫0∞e-stϕ̂nl(t)dt
we obtain
(3.5)ϕnl(s)=∑k=0N-n-1ak(s)ϕn+k-1,l(s)+∑k=N-n∞ak(s)ϕN-1,l(s)+rnl(s),0<n≤N,(3.6)ϕ0l(s)=λs+λϕ1l(s)+δ0ls+λ,
where
(3.7)ak(s)=∫0∞e-(λ+s)t(λt)kk!dF(t),dk(s)=∫0∞e-(λ+s)t(λt)kk!(1-F(t))dt,rkl(s)={0ifl<k,dl-k(s)ifk≤l<N,(1-f(s))s-∑i=0N-n-1di(s)ifl=N,f(s)=∫0∞e-stdF(t).
Substituting φn(s)=ϕN-n,l(s), we get from (3.5), (3.6)
(3.8)∑k=0n+1ak(s)φn-k+1(s)-φn(s)=ψn(s),0≤n<N,(3.9)φN(s)=λs+λφN-1(s)+δ0ls+λ,
with
(3.10)ψn(s)=an+1(s)φ0(s)-∑k=n+1∞ak(s)φ1(s)-rN-n,l(s).(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,
(3.11)φn(s)=Rn+1(s)c(s)+∑k=0nRn-k(s)ψk(s),
where c(s) does not depend on n and Rk(s) is given in (2.3) for the sequence ak(s).Now we only need to find unknownc(s), φ0(s), and φ1(s). In order to find c(s), we put n=0 into (3.11) and get c(s)=φ0(s)/R1(s). Putting n=0 into (3.8) and observing that ∑k=0∞ak(s)=f(s) yield φ1(s)=(φ0(s)-rNl(s))/f(s). Then, substituting these results into (3.11) we have
(3.12)φn(s)=φ0(s)cn(s)+hnl(s),(3.13)hkl(s)=∑i=0kRk-i(s)[rNl(s)(1-1f(s)∑j=0iaj(s))-rN-i,l(s)],(3.14)ck(s)=Rk+1(s)a0(s)+∑i=0kRk-i(s)bi(s),(3.15)bk(s)=ak+1(s)+1f(s)∑i=0kai(s)-1.To calculateφ0(s), we set n=N and n=N-1 in (3.12) and use the boundary condition (3.9). This gives
(3.16)φ0(s)=λhN-1,l(s)-(s+λ)hNl(s)+δ0l(s+λ)cN(s)-λcN-1(s).
Using (3.16) and (3.12) with n′=N-n, 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(3.17)ϕnl(s)=cN-n(s)λhN-1,l(s)-(s+λ)hN,l(s)+δ0l(s+λ)cN(s)-λcN-1(s)+hN-n,l(s),0≤n≤N,
where hkl(s) and ck(s) 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.,N=20), the arrival rate is 1 (i.e., λ=1), and the service time is constant and equal to 0.9. Therefore, we have ρ=0.9—the traffic intensity is moderate. However, we assume that the system is initially full (i.e., X(0)=20).Using Theorem3.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.Table 1
The mean queue size and the standard deviation in the M/D/1/20 system at different moments in time.
Mean queue sizeStandard deviationt=119.080.29t=218.780.63t=518.231.23t=1016.601.88t=2014.982.92t=5011.304.95t=1007.195.33t=2004.974.40t=5004.624.13t=∞4.624.13Figure 1
Queue size distributions in the M/D/1/20 system at different moments in time.As we can see, for high values oft, the distribution of the queue size reaches its maximum at l=1. To explain this, we note first that ρ<1. This causes that for growing t 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 l=1 is higher than at l=0.In Table1 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 t∈(50,200). This can be explained in the following way. As we start from a queue size of 20, for a small t the probability mass is concentrated around 20. Moreover, as we also have N=20, the queue cannot get longer than 20. Therefore, for a small t, the distribution has only the left tail, and its variance is small. Now, for a large t, 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 t, 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 t.


## 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.,N=20), the arrival rate is 1 (i.e., λ=1), and the service time is constant and equal to 0.9. Therefore, we have ρ=0.9—the traffic intensity is moderate. However, we assume that the system is initially full (i.e., X(0)=20).Using Theorem3.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.Table 1
The mean queue size and the standard deviation in the M/D/1/20 system at different moments in time.
Mean queue sizeStandard deviationt=119.080.29t=218.780.63t=518.231.23t=1016.601.88t=2014.982.92t=5011.304.95t=1007.195.33t=2004.974.40t=5004.624.13t=∞4.624.13Figure 1
Queue size distributions in the M/D/1/20 system at different moments in time.As we can see, for high values oft, the distribution of the queue size reaches its maximum at l=1. To explain this, we note first that ρ<1. This causes that for growing t 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 l=1 is higher than at l=0.In Table1 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 t∈(50,200). This can be explained in the following way. As we start from a queue size of 20, for a small t the probability mass is concentrated around 20. Moreover, as we also have N=20, the queue cannot get longer than 20. Therefore, for a small t, the distribution has only the left tail, and its variance is small. Now, for a large t, 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 t, 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 t.


## 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 theMX/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 {p0,p1,p2,…}, where ∑i=0∞pi=1. 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 F(t), 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τn, where τn is defined as follows:
(4.1)τn=inf⁡{t>0:X(t)=N∣X(0)=n},
and X(t) denotes the number of customers in the system at time t+.(I) Using the total probability formula with respect to the first departure epoch for initially nonempty system,0<X(0)<N, we have
(4.2)P(τn>t)=∑k=0N-n-1∫0tP(τn+k-1>t-u)mk(u)dF(u)+(1-F(t))∑k=0N-n-1mk(t),0<n<N,
where mk(t) denotes the probability that k customers arrive in interval (0,t].The first term in (4.2) describes the situation where the first departure time u is before t and there is no buffer overflow by the time u, which means that the number of arrivals in (0,u] must not exceed N-n-1. The second term describes the situation where the first departure time is after t and there is no buffer overflow by the time t. Naturally, the situation where an overflow occurs in interval (0,t] is not taken into account now, as in this case we have P(τn>t)=0.If the system is initially empty, then conditioning on the first arrival epoch we get(4.3)P(τ0>t)=∑k=0N-1pk∫0tP(τk>t-u)λe-λudu+e-λt.(II) The Laplace transform applied to (4.2) and (4.3) reduces the problem to
(4.4)ln(s)=∑k=0N-n-1ln+k-1(s)ak(s)+d̃N-n(s),0<n<N,l0(s)=λs+λ∑k=0N-1pklk(s)+1s+λ,
with
(4.5)ln(s)=∫0∞e-stP(τn>t)dt,ak(s)=∫0∞e-stmk(t)dF(t),d̃k(s)=∑i=0k-1∫0∞e-stmi(t)(1-F(t))dt.
Substituting ln(s)=uN-n(s) we obtain
(4.6)∑k=-1n-1un-k(s)ak+1(s)-un(s)=ψn(s),0<n<N,(4.7)uN(s)=λs+λ∑k=0N-1pkuN-k(s)+1s+λ,
where ψn(s)=u1(s)an(s)-d̃n(s).(III) Now, applying Lemma2.1 the general solution of the system (4.6) has the form
(4.8)un(s)=c(s)Rn(s)+∑k=1nψk(s)Rn-k(s),n>0,
where c(s) does not depend on n and Rk(s) is given in (2.3) for ak(s) defined in (4.5).Puttingn=1 in (4.8), we can observe that c(s)=u1(s)/R1(s). Then, using condition (4.7) together with (4.8), we have
(4.9)u1(s)=vN(s)wN(s),vN(s)=λs+λ∑k=1NpN-k∑i=1kd̃i(s)Rk-i(s)-∑k=1Nd̃k(s)RN-k(s)-1s+λ,(4.10)wN(s)=λs+λ∑k=0NpN-k∑i=0kai(s)Rk-i(s)-∑k=0Nak(s)RN-k(s).
Finally, rewriting (4.8) as
(4.11)un(s)=u1(s)∑k=0nak(s)Rn-k(s)-∑k=1nd̃k(s)Rn-k(s),
we arrive at the following theorem.Theorem 4.1.
The transform of the time to buffer overflow in theMX/G/1/N system has the form
(4.12)ln(s)=vN(s)wN(s)∑k=0N-nak(s)RN-n-k(s)-∑k=1N-nd̃k(s)RN-n-k(s),0≤n<N,
where vN(s) and wN(s) are given in (4.9) and (4.10), respectively.To make this theorem useful, we have to be able to computeak(s) and d̃k(s). Computing these coefficients is not very demanding and may be carried out, for instance, using generating functions. It is easy to check that
(4.13)a(z,s)=∑k=0∞zkak(s)=f(s+λ(1-p(z))),p(z)=∑k=0∞zkpk,d̃(z,s)=∑k=1∞zkd̃k(s)=z[1-f(s+λ(1-p(z)))](1-z)(s+λ(1-p(z))).
A very effective algorithm for generating function inversion can be found in [20]. Namely, if we have a generating function q(z)=∑k=0∞qkzk, then the original values of qk can be restored as
(4.14)qk≈12klrk(a0(k,l,r)+(-1)kak(k,l,r)+2∑j=1k-1(-1)jRe(aj(k,l,r))),
where
(4.15)aj(k,l,r)=∑n=0l-1e-πin/lq(reπi(n+lj)/lk),
while l and r are used to control the roundoff error. (Typically, we use l=1, r=10-4/k.)An alternative way of computingak(s) and d̃k(s) is the uniformization technique [21]. Applying this technique to ak(s), we obtain
(4.16)ak(s)=∑j=0∞γj(s)Kk,j,
where
(4.17)γj(s)=1j!∫0∞e-t(s+1)tjdF(t),
and Kk,j can be computed as follows:
(4)K0,0=1,Kk,0=0,k≥1,K0,j=0,j≥1,Kk,j=∑i=0k-1Ki,j-1pk-i.
Similarly, for d̃k(s), we get
(4.19)d̃k(s)=∑i=0k-1∑j=0∞δj(s)Ki,j,
with
(4.20)δj(s)=1j!∫0∞e-t(s+1)tj(1-F(t))dt.For the bibliography on other computational methods forτn, 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:
(4.21)p1=16,p3=12,p7=13,λ=12.
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 2/5 (which gives ρ=4/5) and the buffer size is 100.Suppose we want to compute the average time to buffer overflow, starting from an empty buffer, that is,(4.22)Eτ0=∫0∞tdP(τ0<t).
It is easy to see that
(4.23)Eτ0=l0(0).The value ofl0(0) can be computed using the uniformization technique. From (4.17) and (4.20), we obtain, respectively,
(4.24)γj(0)=e-2/5(2/5)jj!,δj(0)=Γ(j+1,0)-Γ(j+1,2/5)j!,
where Γ(j,x) denotes the incomplete gamma function. Now, using (4.16) and (4.19) we can compute ak(0) and d̃k(0). Finally, applying (4.12), we get
(4.25)Eτ0=82596.64.
As we have a moderate traffic intensity and a quite big buffer, a large time to buffer overflow was to be expected.


## 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:
(4.21)p1=16,p3=12,p7=13,λ=12.
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 2/5 (which gives ρ=4/5) and the buffer size is 100.Suppose we want to compute the average time to buffer overflow, starting from an empty buffer, that is,(4.22)Eτ0=∫0∞tdP(τ0<t).
It is easy to see that
(4.23)Eτ0=l0(0).The value ofl0(0) can be computed using the uniformization technique. From (4.17) and (4.20), we obtain, respectively,
(4.24)γj(0)=e-2/5(2/5)jj!,δj(0)=Γ(j+1,0)-Γ(j+1,2/5)j!,
where Γ(j,x) denotes the incomplete gamma function. Now, using (4.16) and (4.19) we can compute ak(0) and d̃k(0). Finally, applying (4.12), we get
(4.25)Eτ0=82596.64.
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 functionF(t)), and finite capacity N.The Markovian arrival process (MAP) is one of the most flexible arrival processes from the classN 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 twom×m matrices, D0 and D1, such that D1 is nonnegative, D0 has nonnegative off-diagonal elements, and negative diagonal elements and D=D0+D1≠D0 is an irreducible infinitesimal generator. We will use J(t) to denote the state of the underlying Markov chain, N(t) to denote the number of arrivals in (0,t], and Pij(n,t) to denote the counting function, that is,
(5.1)Pij(n,t)=P(N(t)=n,J(t)=j∣N(0)=0,J(0)=i).
We will also use intensities λi and probabilities pi(k,j), k=0,1, defined as
(5.2)λi=-(D0)ii,pi(1,j)=1λi(D1)ij,1≤i,j≤m,pi(0,j)=1λi(D0)ij,1≤i,j≤m,j≠i.By the workloadV(t) we mean the length of time a job (packet) which arrives at time t 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(5.3)wn,i(s1,s2)=∫0∞e-s2tdt∫0∞e-s1xw̃n,i(x,t)dx,w̃n,i(x,t)=P(V(t)>x∣X(0)=n,J(0)=i),
in the column vector form
(5.4)wn(s1,s2)=(wn,1(s1,s2),…,wn,m(s1,s2))T.As previously,X(t) denotes the number of customers in the system at time t+.(I) As in the previous sections, we start from using the total probability formula with respect to the first departure moment. For0<n≤N, 1≤i≤m, we obtain
(5.5)w̃n,i(x,t)=∑j=1m∑k=0N-n-1∫0tw̃n+k-1,j(x,t-u)Pi,j(k,u)dF(u)+∑j=1m∑k=N-n∞∫0tw̃N-1,j(x,t-u)Pi,j(k,u)dF(u)+∑j=1m∑k=0N-n-1Pi,j(k,t)∫t∞(1-F(n+k-1)*(x-u+t))dF(u),
where F(k)* is the k-fold convolution of the distribution function F with itself.The first double sum in (5.5) describes the case where the first departure time u is before t and there is no buffer overflow by the time u. The second double sum describes the case where the first departure time u is before t and an overflow occurs by the time u, which means that the number of arrivals is equal to N-n or more. The third double sum describes the case where the first departure time u is after t and there is no overflow by the time t. In this case, we have P(V(t)>x)=1-F(n+k-1)*(x-u+t), where k is the number of arrivals in (0,t].If the system is initially empty, then for1≤i≤m we get
(5.6)w̃0,i(x,t)=∑j=1m∑k=01∫0tw̃k,j(x,t-u)pi(k,j)λie-λiudu.(II) Applying transforms and matrix notation to (5.5) and (5.6) we obtain
(5.7)wn(s1,s2)=∑k=0N-n-1Ak(s2)wn+k-1(s1,s2)+∑k=N-n∞Ak(s2)wN-1(s1,s2)+qn(s1,s2),0<n≤N,(5.8)w0(s1,s2)=∑k=01Yk(s2)wk(s1,s2),
where
(5.9)Ak(s)=[∫0∞e-stPi,j(k,t)dF(t)]i,j,Yk(s)=[λipi(k,j)s+λi]i,j,D¯k(s)=[∫0∞e-stPi,j(k,t)(1-F(t))dt]i,j,Ck(s1,s2)=[∫0∞e-s2tPi,j(k,t)dt∫0∞e-s1xdxF(x+t)]i,j,qn(s1,s2)=1s1∑k=0N-n-1[D¯k(s2)-fn+k-1(s1)Ck(s1,s2)]⋅1,1=(1,…,1)T.
Replacing νn(s1,s2)=wN-n(s1,s2) and (5.7) gives
(5.10)∑k=0n+1Ak(s2)νn-k+1(s1,s2)-νn(s1,s2)=Ψn(s1,s2),0≤n<N,Ψn(s1,s2)=An+1(s2)ν0(s1,s2)-∑k=n+1∞Ak(s2)ν1(s1,s2)-qN-n(s1,s2).(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(5.11)wn(s1,s2)=∑k=0N-nRN-n-k(s2)hk(s1,s2)+(∑k=-1N-nRN-n-k(s2)Bk(s2)+RN-n+1(s2))MN-1(s2)uN(s1,s2),0≤n≤N,
with
(5.12)hn(s1,s2)=A¯n+1(s2)(A¯0(s2))-1qN(s1,s2)-qN-n(s1,s2),A¯k(s)=∑i=k∞Ai(s),Bk(s)=Ak+1(s)-A¯k+1(s)(A¯0(s))-1,MN(s)=∑k=0N+1RN-k+1(s)Bk-1(s)+RN+1(s)-∑k=N-1NYN-k(s)[∑i=0k+1Rk-i+1(s)Bi-1(s)+Rk+1(s)],uN(s1,s2)=∑k=N-1NYN-k(s2)∑i=0kRk-i(s2)hi(s1,s2)-∑k=0NRN-k(s2)hk(s1,s2).Note that matricesAk, D¯k, and Ck 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 lims1,s2→0+s2wn(s1,s2). Putting s1=0+ into wn(s1,s2) and inverting the result with respect to s2 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 t.It is easy to check that the number of floating-point operations needed to compute (5.11) (time complexity) grows as O(m3N2). This estimate is a consequence of the form of (5.11) and (2.3) and the fact that matrix multiplication and inversion are of O(m3) 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 O(m3N3) 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:
(5.13)D0=[-11188.00145.53845.35816.76173.27-4786.97364.04202.79729.86739.55-11958.43236.94191.44791.87105.75-11481.10],D1=[2467.302907.271055.812949.981229.92609.921178.331028.702273.492442.683166.512369.401932.711522.553761.293175.49].
The average rate of the fitted MAP is
(5.14)Λ=π⋅D1⋅1=7608.46pkts/s,
where π=(0.18537,0.40467,0.20238,0.20758) is the stationary distribution for the underlying Markov chain J(·) and 1=(1,…,1)T.It is assumed that the service time is constant and equal tod. Manipulating d we can easily obtain different traffic intensities ρ=Λd.In Figure2, the mean workload as a function of time, EV(t), for the initial buffer occupancy of 0%, 25%, 50%, 75%, and 100% is depicted. The traffic intensity was set for ρ=0.95, the initial phase for J(0)=1, and the buffer size for 100 packets.Figure 2
Mean workload versus time for the initial buffer occupancy of 0%, 25%, 50%, 75%, and 100%, counting from the bottom.ρ=0.95, J(0)=1, N=100 pkts.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 Figure3, the stationary mean workload, limt→∞EV(t), 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 ρ=0.75 this border buffer size is about 20, for ρ=0.90 about 50, for ρ=0.95 about 100, and for ρ=0.99 about 500.Stationary mean workload versus the buffer size for different traffic intensities. (a):ρ=0.75. (b): ρ=0.90. (c): ρ=0.95. (d): ρ=0.99.
(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.


## 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:
(5.13)D0=[-11188.00145.53845.35816.76173.27-4786.97364.04202.79729.86739.55-11958.43236.94191.44791.87105.75-11481.10],D1=[2467.302907.271055.812949.981229.92609.921178.331028.702273.492442.683166.512369.401932.711522.553761.293175.49].
The average rate of the fitted MAP is
(5.14)Λ=π⋅D1⋅1=7608.46pkts/s,
where π=(0.18537,0.40467,0.20238,0.20758) is the stationary distribution for the underlying Markov chain J(·) and 1=(1,…,1)T.It is assumed that the service time is constant and equal tod. Manipulating d we can easily obtain different traffic intensities ρ=Λd.In Figure2, the mean workload as a function of time, EV(t), for the initial buffer occupancy of 0%, 25%, 50%, 75%, and 100% is depicted. The traffic intensity was set for ρ=0.95, the initial phase for J(0)=1, and the buffer size for 100 packets.Figure 2
Mean workload versus time for the initial buffer occupancy of 0%, 25%, 50%, 75%, and 100%, counting from the bottom.ρ=0.95, J(0)=1, N=100 pkts.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 Figure3, the stationary mean workload, limt→∞EV(t), 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 ρ=0.75 this border buffer size is about 20, for ρ=0.90 about 50, for ρ=0.95 about 100, and for ρ=0.99 about 500.Stationary mean workload versus the buffer size for different traffic intensities. (a):ρ=0.75. (b): ρ=0.90. (c): ρ=0.95. (d): ρ=0.99.
(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φn(s)=ϕN-n(s) in (3.5)), which cannot be carried out in the infinite-buffer model.



---

*Source: 831956-2012-02-01.xml*