# TheN×D-BMAP/G/1 Queueing Model: Queue Contents and Delay Analysis

Authors: Bart Steyaert; Joris Walraevens; Dieter Fiems; Herwig Bruneel
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We consider a single-server discrete-time queueing system with N sources, where each source is modelled as a correlated Markovian customer arrival process, and the customer service times are generally distributed. We focus on the analysis of the number of customers in the queue, the amount of work in the queue, and the customer delay. For each of these quantities, we will derive an expression for their steady-state probability generating function, and from these results, we derive closed-form expressions for key performance measures such as their mean value, variance, and tail distribution. A lot of emphasis is put on finding closed-form expressions for these quantities that reduce all numerical calculations to an absolute minimum.


---

## Body

## 1. Introduction

In this paper, we analyze a queueing model with generally distributed customer service times and a correlated arrival process.  This kind of model is, for instance, useful in assessing the performance of a packet-switched communication infrastructure, where messages carried by the network can have a variable transmission time, such as the current Internet, or a dedicated packet-based network carrying video-on-demand (VoD).  Multiplexer queues and/or switches and routers in such a network can in general be modelled by means of a discrete-time queueing system, where new packets (i.e., “customers”) are generated by a superposition of individual traffic sources.  The size of the packets is proportional to their transmission times, whose distribution depends on the specific application under consideration.  Analyzing such a system is mandatory in the design and evaluation of these networks.  However, this can be a difficult task because of the time-correlated behaviour that each of the individual sources may exhibit.To facilitate the queueing analysis, a source is often modelled as a discrete-time Markov modulated arrival process, such as the discrete-time batch Markovian arrival process (D-BMAP; [1, 2]) with the Markov modulated Bernoulli process (MMBP; [3–5]) and the switched batch bernoulli process (SBBP; [6, 7]) as frequently used special cases.  As a consequence, the analysis of a queueing model with a specific SBBP [8] or MMBP [4, 9–11] source model or with a general D-BMAP (e.g., [1, 12–21], and the vast amount of references therein) has become the focus of many research papers.  In addition to the characteristics of the arrival process, such queueing models can also be classified by means of the service process.  For customer service times equal to one slot, single-server queueing models are considered in [4, 9, 12, 15], while [18, 21] focus on the multiserver case.  For the single-server case, geometrically [1, 10] and phase-type [13] distributed customer service times have been considered as well, while these quantities are generally distributed in [8, 11, 14, 17, 19, 20].  In this paper, we investigate the scenario where arrivals are generated by an aggregation of N sources, where each source is characterized by means of a (two-state) D-BMAP, and the sources are assumed to be homogeneous, that is, they are all described by the same stochastic process.  The service times on the other hand are assumed to be generally distributed.Although the queueing model that is studied is not entirely new, the novelty of this contribution lies in the way by which the problem is tackled and solved.  A widely used method to deal with correlated arrival processes is thematrix-analytic solution technique [22, 23], which provides a broad and generic framework for the analysis of a wide set of queuing systems, both for finite- [4, 9–12, 14, 15, 19] and infinite-capacity [8, 24] queues.  This is based on an intelligent matrix representation of the quantities that are needed in these analyses and making full use of the properties that these matrices exhibit.  One impediment it may nevertheless sometimes suffer from is the size of the matrices that must be manipulated, which is (roughly speaking) determined by the state spaceof the system. An alternative solution technique, adopted in [1, 16–18], is the so-called spectral decompositionsolution method; however, the main results that are generated still contain operations that can be computationally demanding, such as Krönecker sums and products, or inversions, of matrices that can have large dimensions.  The approach that we adopt—based on the work presented in [21] for single-slot service times—is related to the latter one,  albeit that we follow a different path by invoking as much as possible a representation by means of (joint) probability generating functions(pgfs) in our derivations presented in the remainder of this paper.  The main advantage lies in the observation that the formulae for the performance indices that we finally obtain, such as the moments of the queue contents, unfinished work, and customer delay, are—to a large extent—expressed directly in terms of the system parameters and therefore often lend themselves to an interpretation that is helpful for understanding the mechanisms that determine the queue behaviour. In addition, in order to establish (semi)analytic expressions for the relevant performance measures that pertain to the queue behaviour, a (possibly large) set of boundary probabilities has to be calculated numerically. To reduce the numerical work as much as possible, we introduce approximations for these quantities that allow us to compute the aforementioned performance measures without great difficulty. By means of some numerical examples we will demonstrate the efficiency and accurateness of these approximations, as well as the influence of the system parameters on the queue behaviour.


## 2. System Description

We consider a discrete-time queueing model, that is, a system where the time axis is divided into fixed-length contiguous intervals, referred to asslots.  The system consists of one single server and an infinite waiting room for customers awaiting service.  Customers arrive in the queue according to a correlated Markovian arrival process, described further on.  It is assumed that the service of a customer requires a positive integer number of slots, described by a general service time distribution, and can start (and end) at slot marks only, that is, at the time points between consecutive slots, meaning that the customer service is synchronized with respect to the slot marks.  This is not necessarily the case as far as customer arrivals are concerned: customers may enter the queue at any continuous time instant.  Nevertheless, due to the synchronous nature of a customer’s service, it can commence at the earliest beginning of the slot following its arrival (if it arrives in an empty queue).  Therefore, the precise details of the position of the customer arrival instants within a slot are irrelevant, and it suffices to characterize the arrival process by a random variable describing the total number of customer arrivals during a slot.  We will now discuss the customer arrival and service processes in more detail.Consider a buffer model fed byN identical and independent sources generating customers with variable service times.  Each source is modelled as a 2-state Markov modulated arrival process, where the state of a source during a slot is represented by Si, 1≤i≤2.  The customer arrival process during a slot is completely characterized by the 2×2 matrix
(2.1)Q(z)=[q11(z)q12(z)q21(z)q22(z)].
State transitions occur at slot boundaries, and let us denote by pij, 1≤i,j≤2, the one-step transition  probability Pr[Si→Sj] at the end of slot k-1, where of course pj1+pj2=1,1≤j≤2. The elements qij(z) of Q(z) are then given by (2.2)qij(z)=̂Gij(z)pij,
where Gij(z), 1≤i,j≤2, is the probability generating function (pgf) describing the number of customers generated during a slot by a source, given that the source is in state Sj  during the tagged slot and was in state Si during the preceding slot, that is,(2.3)Gij(z)=̂E[zek(n)∣Si⟶Sjattheendofslotk-1],
where ek(n), 1≤n≤N, represents the number of customer arrivals generated by source n during slot k and E[·] denotes the expected value of the tagged quantity. For the Markov modulated Bernoulli process (MMBP), the number of customers generated by a source during any slot is either zero or equal to one, meaning that each of the generating functions Gij(z) is a linear function of z, that is, (2.4)Gij(z)=1-gij+gijz,
for some parameters gij satisfying 0≤gij≤1.  Although attention is focused on this specific arrival model in the numerical examples, the analysis is general and can also be applied when the Gij(z)’s have a more complex form.  Moreover, the notational conventions have been chosen such that the derivations and results that are presented in the subsequent sections can be readily extended to an L-state D-BMAP as well (see also [20]), L≥2. We confine ourselves to the case of L=2 in this paper in order to enhance the readability of the derivations that follow hereafter, while their level of complexity is that of arbitrary L.The aggregate customer arrival process is fully determined, once the probability generating matrixQ(z) has been specified.  Let us define ak,j,1≤j≤2, as the number of sources in states Sj during slot k. Note that, due to the fact that the total number of sources equals N, these random variables satisfy(2.5)∑j=12ak,j=N.
In the rest of the paper, we will denote by x the 2×1 column vector (x1,x2)T (where (·)T  represents the matrix transposition operation),  and similarly we write ak=(ak,1,ak,2).  Let us also denote by the 2×1 column vector B(x,z) the matrix product
(2.6)B(x,z)=̂(B1(x,z),B2(x,z))T=Q(z)x,
and, for a set of random variables r=(r1,r2), define the joint generating function as(2.7)E[xr]=̂E[∏j=12xjrj].
Then, with the previous definitions, it is not difficult to show that the joint probability generating function of the random variables denoting the number of customer arrivals and the state of the arrival process during a slot can be written in terms of the state of the arrival process during the previous slot, leading to(2.8)E[zek+1xak+1]=E[zek+1∏j=12xjak+1,j]=E[∏j=12Bj(x,z)ak,j]=E[(Q(z)x)ak],
where ek represents the total number of customer arrivals during slot k, that is,(2.9)ek=̂∑n=1Nek(n).
Equation (2.8) fully describes the number of customers generated during consecutive slots by the N customer sources.Starting from some initial state, the customer arrival process will evolve to a stochastic equilibrium after a sufficiently long period of time, and we defineA(x) as the joint probability generating function of the number of sources in states S1 and S2 during an arbitrary slot. If we define σ1 and σ2 as the steady-state probability that the Markov state of a source during an arbitrary slot is S1 and S2, respectively, and σ as the 2×1 vector (σ1,σ2)T, which is the solution of the matrix equation(2.10)σT=σTQ(1),σT1=1,
leading to(2.11)σ1=p21p12+p21,σ2=p12p12+p21,
where 1 is the 2×1 column vector with all elements equal to 1, then A(x) equals(2.12)A(x)=(σTx)N.
It is further assumed that the service times of consecutive customers that arrive in the queue form a set of independent and identically distributed random variables, denoted by s, with common probability mass function(2.13)s(l)=Pr[customer service time=lslots],l≥1,
and corresponding probability generating function(2.14)S(z)=E[zs]=∑l=1∞s(l)zl.
Obviously, we assume that the service time of a customer is at least one slot.Due to the infinite-queue capacity, the system will reach a stochastic equilibrium only if the equilibrium condition is satisfied.  Definingρ as the load of the system, this implies that(2.15)ρ=̂NpLp<1
must hold, where p denotes the mean number of customer arrivals per slot and per source, and Lp=S′(1) represents the mean customer service time.  Taking into account the customer arrival process described in this section, it follows that p can be calculated from(2.16)p=ddz[σTQ(z)1]z=1=∑i=12σi∑j=12qij′(1),
where primes denote derivatives with respect to the argument.


## 3. The Joint pgf of the State Vector

### 3.1. System Equations

We first establish the system equations that control the evolution of the number of customers in the queue during consecutive slots. This is inspired by the work presented in [25], where the case of an i.i.d. (i.e., uncorrelated) customer arrival process was considered and analysed by means of the so-called supplementary variable method.  Let us therefore define the random variable vk as the system contents at the beginning of slot k+1, which is the number of customers in the queue at the beginning of slot k+1, including the one being served (in case of a nonempty queue).  In addition to vk, we also need information about the amount of time the customer in service still requires at the beginning of slot k+1 before being completely served.  We therefore define the random variable hk as the residual service time, which is the number of slots required at the beginning of slot k+1 to complete the service of the customer being served in case of a nonempty queue; when vk=0, we automatically have hk=0.Together with the definitions of the previous section, we can now establish the system equations.  We must distinguish between the three following cases.(1)
hk=0(⇔vk=0). When the queue is empty at the beginning of a slot, the number of customers in the queue at the beginning of the next slot equals the number of new arrivals, and we find that
(3.1)vk+1=ek+1.
If no customers have entered the queue during slot k+1, the residual service time remains zero, otherwise it is equal to the service time of a new customer, which leads to
(3.2)hk+1={0if(ek+1=0),sif(ek+1>0).(2)
hk=1. This implies that the customer in service is completely served at the end of slot k+1:
(3.3)vk+1=vk+ek+1-1.
Similar to the previous case, if no customers have entered the queue during slot k+1 and the customer in service was the only one in the queue, then the queue becomes empty and the residual service time equals zero at the beginning of the next slot, otherwise the service of a new customer commences, that is,
(3.4)hk+1={0,if(ek+1=0),(vk=1),s,if(ek+1>0)or(vk>1).(3)
hk>1. In this case, the customer being served receives an extra slot of service without the service being completed, and we obtain
(3.5)vk+1=vk+ek+1,hk+1=hk-1.


### 3.2. Derivation of a Functional Equation

Clearly, in view of (3.1)–(3.5), we need to keep track of the random variables vk and hk in our state description.  Consequently, due to the correlated first-order Markov nature of the process controlling the number of customer arrivals during a slot, it becomes clear that the set of random variables (vk,hk,ak) constitutes a four-dimensional Markov state description of the system at the beginning of consecutive slots; this set of four random variables will be referred to as the state vector.  Note that, in view of (2.5), one of the random variables (ak,1,ak,2) can be omitted from the state description.  However, due to reasons of symmetry, we prefer to maintain both random variables in the state description.Let us now define the joint probability generating function of(vk,hk,ak) as(3.6)Pk(z,y,x)=̂E[zvkyhkxak]=E[zvkyhk∏j=12xjak,j].
Let us, for the sake of notation, define a partial generating function of a random variable X if an event A occurs as(3.7)E[zX{A}]=̂E[zX∣A]Pr[A].
From system equations (3.1) and (3.2), taking into account that hk=0⇔vk=0, we then first of all derive that
(3.8)E[zvk+1yhk+1xak+1{vk=0,hk=0}]=E[zek+1ysxak+1{vk=0,hk=0}]+E[(y0-ys)xak+1{hk=0,vk=0,ek+1=0}].
In a similar way, we obtain from (3.3) and (3.4)
(3.9)E[zvk+1yhk+1xak+1{hk=1}]=E[zvk+ek+1-1ysxak+1{hk=1}]+E[(y0-ys)xak+1{hk=1,vk=1,ek+1=0}].
Finally, (3.5) can also be translated into a relation between z-transforms, leading to
(3.10)E[zvk+1yhk+1xak+1{hk>1}]=E[zvk+ek+1yhk-1xak+1]-y-1E[zek+1xak+1{hk=0,vk=0}]-E[zvk+ek+1xak+1{hk=1}].
Summation of (3.8)–(3.10) yields(3.11)yPk+1(z,y,x)=E[zvk+ek+1yhkxak+1]+(yS(y)-1)E[zek+1xak+1{hk=0,vk=0}]+y∑i=01(1-S(y))E[zek+1xak+1{hk=i,vk=i,ek+1=0}]+y(S(y)-z)E[zvk+ek+1-1xak+1{hk=1}].
Note that, given that ak is known, the value of ak+1(and the number of customer arrivals ek+1) is independent of (vk,hk). Therefore, applying (2.8), the previous relation can be transformed into(3.12)yPk+1(z,y,x)=Pk(z,y,Q(z)x)+y(1-S(y))φk(Q(0)x)+(yS(y)-1)φk(Q(z)x)+y(1-S(y))ψk(0,Q(0)x)+y(S(y)-z)ψk(z,Q(z)x),
where
(3.13)ϕk(x)=̂E[xak{uk=0,hk=0}]=E[∏j=12xjak,j{uk=0,hk=0}],ψk(z,x)=̂E[zvk-1xak{hk=1}]=E[zvk-1∏j=12xjak,j{hk=1}].
First of all, due to the foregoing definitions, we have that Pk(0,y,x)≡Pk(0,0,x)=φk(x).  If we insert z=0 into (3.12), we derive that
(3.14)φk+1(x)=φk(Q(0)x)+ψk(0,Q(0)x).
Let us now assume that the equilibrium condition (2.15) is satisfied, implying that all the functions that occur in (3.12) evolve to a steady-state limit, which will be reflected in the remainder by suppressing the subscript k (which expressed the time dependence in the previous derivations).  Together with the previous relation, we finally obtain the functional equation
(3.15)yP(z,y,x)=P(z,y,Q(z)x)+y(1-S(y))φ(x)+(yS(y)-1)φ(Q(z)x)+y(S(y)-z)ψ(z,Q(z)x).
Equation (3.15), together with definition (3.13), defines a functional equation for P(z,y,x), the steady-state joint generating function of the random variables (vk,hk,ak), that is, the system contents and residual service time at the beginning of an arbitrary slot, and the state of the Markovian customer arrival process during the preceding slot.


### 3.3. Solution of the Functional Equation

In order to solve the functional equation, we follow a similar approach to that in [21].  First, if we substitute consecutive arguments x=Q(z)ix′,0≤i≤n, into (3.15), we obtain
(3.16)P(z,y,x)=y-nP(z,y,Q(z)nx)+(1-S(y))∑i=0n-1y-iϕ(Q(z)ix)+(yS(y)-1)∑i=1ny-iϕ(Q(z)ix)+y(S(y)-z)∑i=1ny-iψ(z,Q(z)ix).
We now let n approach infinity. Let us therefore express Q(z)n in terms of its eigenvalues and eigenvectors as W(z)Λ(z)nU(z). For further details concerning the eigenvalues and eigenvectors of Q(z), we refer to the appendix.  Define the functions Flm(z),  1≤l,  m≤N,  as(3.17)(∑j=12u1j(z)xj)l(∑j=12u2j(z)xj)N-l=̂∑m=0NFlm(z)(x1mx2N-m).
In the above expression, the functions uij(z) are the components of the column eigenvectors of Q(z), given by (A.5). One can easily derive that this relation implies that these functions can be written as
(3.18)Flm(z)=∑j=(l+m-N)+min⁡{l,m}(lj)(N-lm-j)u11(z)ju12(z)l-ju21(z)m-ju22(z)N-l-m+j,
where (xy)=x!/(y!(x-y)!) and (·)+=max{0,.}.  In view of definition (3.17) and defining
(3.19)pl=̂lim⁡k→∞Pr[vk=0,hk=0,ak,1=l,ak,2=N-l],
for n→∞, the first sum in the right-hand side of (3.16) can be written as
(3.20)∑i=0∞y-iφ(Q(z)ix)=∑i=0∞y-i∑l=0N(∑j=12u1j(z)λj(z)iwj(z)x)l(∑j=12u2j(z)λj(z)iwj(z)x)N-lpl,
where λj(z) and wj(z) constitute the jth eigenvalue and corresponding row eigenvector of Q(z). The Perron-Fröbenius (P-F) eigenvalue is denoted by λ1(z) (see the appendix for further details). Due to definition (3.17), this can be transformed into(3.21)∑i=0∞y-iφ(Q(z)ix)=∑i=0∞y-i∑l=0N∑m=0NFlm(z)(λ1(z)iw1(z)x)m(λ2(z)iw2(z)x)N-mpl.
Let us now consider values of y and z for which |λj(z)|<|y|≤1,1≤j≤2.  Such values of y and z exist; this is for instance the case if |y|=1, and |z|=1and z≠1 (e.g., see [26]) since(3.22)|λ1(z)|<1,|λ2(z)|<1,for|z|=1,z≠1,
unless some special cases for the generating functions Gij(z) defined in (2.3) are considered—for example, the case that the number of customer arrivals per slot is always a multiple of n,n>1—a situation that falls beyond the scope of this paper.  The sum for i in the above expression then converges, and we obtain
(3.23)∑i=0∞y-iφ(Q(z)ix)=∑m=0Ny(w1(z)x)m(w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)pl.
In a similar way we can derive that
(3.24)∑i=1∞y-iφ(Q(z)ix)=∑m=0N(λ1(z)w1(z)x)m(λ2(z)w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)pl,∑i=1∞y-iψ(z,Q(z)ix)=∑m=0N(λ1(z)w1(z)x)m(λ2(z)w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)ψl(z),
where
(3.25)ψl(z)=̂lim⁡k→∞E[zvk-1{hk=1,ak,1=l,ak,2=N-l}].
In addition, note that
(3.26)lim⁡n→∞y-nP(z,y,Q(z)nx)=∑j=0∞zj∑k=0∞yk∑l=0N∑m=0NFlm(z)Pr[v=j,h=k,a1=l,a2=N-l]⋅lim⁡n→∞y-n(λ1(z)nw1(z)x)m(λ2(z)nw2(z)x)N-m,
which, for |λj(z)|<|y|≤1,1≤j≤2, becomes equal to 0.  Summarizing, in view of the previous results, (3.16) becomes
(3.27)P(z,y,x)=∑m=0N(w1(z)x)m(w2(z)x)N-my-Am(z)×∑l=0NFlm(z){[y(S(y)-z)ψl(z)+(yS(y)-1)pl]Am(z)+y(1-S(y))pl},
where Am(z) has been defined as
(3.28)Am(z)=̂λ1(z)mλ2(z)N-m.
Expression (3.27) for the joint probability generating function P(z,y,x) still contains the unknown constants pl, as well as the unknown functions ψl(z).  Let us therefore, for each value of m, consider values of y for which y=Am(z), which will be represented by ym(z).  Note that if |z|≤1, due to (3.22) this solution for y also satisfies |ym(z)|≤1, and, for these values of z and y, P(z,y,x) is finite.  This means that, if we multiply both hand sides of (3.27) by y-Am(z) and consider y=ym(z), then the left-hand side becomes equal to zero, implying that the same must hold for the right-hand side.  For each value of 0≤m≤N, this leads to the following equation:
(3.29)∑l=0NFlm(z){Am(z)(S(Am(z))-z)ψl(z)+S(Am(z))(Am(z)-1)pl}=0.
If we insert this equation into (3.27), we finally obtain
(3.30)P(z,y,x)=∑m=0N(w1(z)x)m(w2(z)x)N-m⋅{1+zy(Am(z)-1)(S(y)-S(Am(z)))(y-Am(z))(z-S(Am(z)))}∑l=0NFlm(z)pl,
This final result for P(z,y,x) expresses the joint generating function of the state of the arrival process during an arbitrary slot, and the queue contents and residual service time at the beginning of the next slot, in terms of the unknown constants pl.  In the next section, it will become clear how these can be calculated or approximated.  Also, starting from this result, we will establish expressions for the generating functions of quantities such as the queue contents, unfinished work, and message delay, and their corresponding moments and tail distribution.


### 3.4. Calculation of the Unknownspl

We can calculate these probabilities, defined in (3.19), by expressing that P(z,y,x) is analytic when the complex variables y and z both lie inside the unit circle, that is, |z|≤1 and |y|≤1.  We have already expressed in the previous section that, for |y|≤1 and y=Am(z),P(z,y,x) is finite.  In view of expression (3.30), the only remaining potential singularities inside the complex unit disk are those values of z (and y) for which z=S(Am(z)).  Indeed, taking into account the definition for Am(z), then from (3.22) and Rouché’s theorem it follows that this equation has exactly one solution inside the complex unit circle.  Let us, for each value of 0≤m≤N, denote this solution by z{m}, that is,
(3.31)z{m}=S(Am(z{m})),|z{m}|≤1,0≤m≤N.
Then, by expressing that z{m}must also be a zero of the corresponding numerator in expression (3.30) for P(z,y,x), we obtain
(3.32)ξm(z{m})=0,0≤m≤N-1,
where(3.33)ξm(z)=̂11-ρ∑l=0NFlm(z)pl.
Note that z{N}=1 (since λ1(1)=1, see (A.6), implying AN(1)=1), which leads to no additional equation for the unknowns.  Together with the normalization condition, which yields ξN(1)=1 (as we will see further on), (3.32) forms a set of N+1 linear equations for the N+1 unknowns pl that can easily be solved.However, in deriving numerical results for the performance measures considered further in this paper, we want to avoid all numerical calculations as much as possible.  This makes sense in the particular case thatN is large, that is, when a high number of sources are multiplexed, since the number of equations to be solved equals N+1.  Therefore we also present an approximation for the boundary probabilities pl.  As before, denote by the random variables e and v the number of customer arrivals during a tagged arbitrary slot and the queue contents at the beginning of the following slot; also, a represents the state of the Markovian arrival process during the tagged slot.  We define the steady-state joint probability
(3.34)ql=̂Pr[a1=l,a2=N-l,e=0].
Obviously, v=0 implies that there have been no customer arrivals during the tagged slot, that is, v=0⇒e=0.  It is therefore clear that the inequality ql≥pl holds.  We will show by some numerical examples that approximating the boundary probabilities by
(3.35)pl≅1-ρPr[e=0]ql
yields excellent approximations for the performance measures of interest.  Here, we have taken into account that Pr[v=0]=1-ρ so that substituting pl by this formula in (3.30) still yields a normalized generating function.An additional advantage of the approximation proposed here is that it avoids the explicit calculation of theql’s, as well as the functions Flm(z) from (3.18), that occur in expression (3.30) for P(z,y,x).  From definition (3.17) and expression (2.8) it is not difficult to show that
(3.36)∑l=0NFlm(z)ql=(Nm)(σTQ(0)u1(z))m(σTQ(0)u2(z))N-m,
where uj(z) represents the jth column eigenvector of Q(z).  In view of the previous, this implies that the following approximation may be applied:(3.37)ξm(z)≅1Pr[e=0](Nm)(σTQ(0)u1(z))m(σTQ(0)u2(z))N-m.
To conclude, note that it immediately follows from (2.8) that Pr[e=0] satisfies
(3.38)Pr[e=0]=(σTQ(0)1)N.


## 3.1. System Equations

We first establish the system equations that control the evolution of the number of customers in the queue during consecutive slots. This is inspired by the work presented in [25], where the case of an i.i.d. (i.e., uncorrelated) customer arrival process was considered and analysed by means of the so-called supplementary variable method.  Let us therefore define the random variable vk as the system contents at the beginning of slot k+1, which is the number of customers in the queue at the beginning of slot k+1, including the one being served (in case of a nonempty queue).  In addition to vk, we also need information about the amount of time the customer in service still requires at the beginning of slot k+1 before being completely served.  We therefore define the random variable hk as the residual service time, which is the number of slots required at the beginning of slot k+1 to complete the service of the customer being served in case of a nonempty queue; when vk=0, we automatically have hk=0.Together with the definitions of the previous section, we can now establish the system equations.  We must distinguish between the three following cases.(1)
hk=0(⇔vk=0). When the queue is empty at the beginning of a slot, the number of customers in the queue at the beginning of the next slot equals the number of new arrivals, and we find that
(3.1)vk+1=ek+1.
If no customers have entered the queue during slot k+1, the residual service time remains zero, otherwise it is equal to the service time of a new customer, which leads to
(3.2)hk+1={0if(ek+1=0),sif(ek+1>0).(2)
hk=1. This implies that the customer in service is completely served at the end of slot k+1:
(3.3)vk+1=vk+ek+1-1.
Similar to the previous case, if no customers have entered the queue during slot k+1 and the customer in service was the only one in the queue, then the queue becomes empty and the residual service time equals zero at the beginning of the next slot, otherwise the service of a new customer commences, that is,
(3.4)hk+1={0,if(ek+1=0),(vk=1),s,if(ek+1>0)or(vk>1).(3)
hk>1. In this case, the customer being served receives an extra slot of service without the service being completed, and we obtain
(3.5)vk+1=vk+ek+1,hk+1=hk-1.


## 3.2. Derivation of a Functional Equation

Clearly, in view of (3.1)–(3.5), we need to keep track of the random variables vk and hk in our state description.  Consequently, due to the correlated first-order Markov nature of the process controlling the number of customer arrivals during a slot, it becomes clear that the set of random variables (vk,hk,ak) constitutes a four-dimensional Markov state description of the system at the beginning of consecutive slots; this set of four random variables will be referred to as the state vector.  Note that, in view of (2.5), one of the random variables (ak,1,ak,2) can be omitted from the state description.  However, due to reasons of symmetry, we prefer to maintain both random variables in the state description.Let us now define the joint probability generating function of(vk,hk,ak) as(3.6)Pk(z,y,x)=̂E[zvkyhkxak]=E[zvkyhk∏j=12xjak,j].
Let us, for the sake of notation, define a partial generating function of a random variable X if an event A occurs as(3.7)E[zX{A}]=̂E[zX∣A]Pr[A].
From system equations (3.1) and (3.2), taking into account that hk=0⇔vk=0, we then first of all derive that
(3.8)E[zvk+1yhk+1xak+1{vk=0,hk=0}]=E[zek+1ysxak+1{vk=0,hk=0}]+E[(y0-ys)xak+1{hk=0,vk=0,ek+1=0}].
In a similar way, we obtain from (3.3) and (3.4)
(3.9)E[zvk+1yhk+1xak+1{hk=1}]=E[zvk+ek+1-1ysxak+1{hk=1}]+E[(y0-ys)xak+1{hk=1,vk=1,ek+1=0}].
Finally, (3.5) can also be translated into a relation between z-transforms, leading to
(3.10)E[zvk+1yhk+1xak+1{hk>1}]=E[zvk+ek+1yhk-1xak+1]-y-1E[zek+1xak+1{hk=0,vk=0}]-E[zvk+ek+1xak+1{hk=1}].
Summation of (3.8)–(3.10) yields(3.11)yPk+1(z,y,x)=E[zvk+ek+1yhkxak+1]+(yS(y)-1)E[zek+1xak+1{hk=0,vk=0}]+y∑i=01(1-S(y))E[zek+1xak+1{hk=i,vk=i,ek+1=0}]+y(S(y)-z)E[zvk+ek+1-1xak+1{hk=1}].
Note that, given that ak is known, the value of ak+1(and the number of customer arrivals ek+1) is independent of (vk,hk). Therefore, applying (2.8), the previous relation can be transformed into(3.12)yPk+1(z,y,x)=Pk(z,y,Q(z)x)+y(1-S(y))φk(Q(0)x)+(yS(y)-1)φk(Q(z)x)+y(1-S(y))ψk(0,Q(0)x)+y(S(y)-z)ψk(z,Q(z)x),
where
(3.13)ϕk(x)=̂E[xak{uk=0,hk=0}]=E[∏j=12xjak,j{uk=0,hk=0}],ψk(z,x)=̂E[zvk-1xak{hk=1}]=E[zvk-1∏j=12xjak,j{hk=1}].
First of all, due to the foregoing definitions, we have that Pk(0,y,x)≡Pk(0,0,x)=φk(x).  If we insert z=0 into (3.12), we derive that
(3.14)φk+1(x)=φk(Q(0)x)+ψk(0,Q(0)x).
Let us now assume that the equilibrium condition (2.15) is satisfied, implying that all the functions that occur in (3.12) evolve to a steady-state limit, which will be reflected in the remainder by suppressing the subscript k (which expressed the time dependence in the previous derivations).  Together with the previous relation, we finally obtain the functional equation
(3.15)yP(z,y,x)=P(z,y,Q(z)x)+y(1-S(y))φ(x)+(yS(y)-1)φ(Q(z)x)+y(S(y)-z)ψ(z,Q(z)x).
Equation (3.15), together with definition (3.13), defines a functional equation for P(z,y,x), the steady-state joint generating function of the random variables (vk,hk,ak), that is, the system contents and residual service time at the beginning of an arbitrary slot, and the state of the Markovian customer arrival process during the preceding slot.


## 3.3. Solution of the Functional Equation

In order to solve the functional equation, we follow a similar approach to that in [21].  First, if we substitute consecutive arguments x=Q(z)ix′,0≤i≤n, into (3.15), we obtain
(3.16)P(z,y,x)=y-nP(z,y,Q(z)nx)+(1-S(y))∑i=0n-1y-iϕ(Q(z)ix)+(yS(y)-1)∑i=1ny-iϕ(Q(z)ix)+y(S(y)-z)∑i=1ny-iψ(z,Q(z)ix).
We now let n approach infinity. Let us therefore express Q(z)n in terms of its eigenvalues and eigenvectors as W(z)Λ(z)nU(z). For further details concerning the eigenvalues and eigenvectors of Q(z), we refer to the appendix.  Define the functions Flm(z),  1≤l,  m≤N,  as(3.17)(∑j=12u1j(z)xj)l(∑j=12u2j(z)xj)N-l=̂∑m=0NFlm(z)(x1mx2N-m).
In the above expression, the functions uij(z) are the components of the column eigenvectors of Q(z), given by (A.5). One can easily derive that this relation implies that these functions can be written as
(3.18)Flm(z)=∑j=(l+m-N)+min⁡{l,m}(lj)(N-lm-j)u11(z)ju12(z)l-ju21(z)m-ju22(z)N-l-m+j,
where (xy)=x!/(y!(x-y)!) and (·)+=max{0,.}.  In view of definition (3.17) and defining
(3.19)pl=̂lim⁡k→∞Pr[vk=0,hk=0,ak,1=l,ak,2=N-l],
for n→∞, the first sum in the right-hand side of (3.16) can be written as
(3.20)∑i=0∞y-iφ(Q(z)ix)=∑i=0∞y-i∑l=0N(∑j=12u1j(z)λj(z)iwj(z)x)l(∑j=12u2j(z)λj(z)iwj(z)x)N-lpl,
where λj(z) and wj(z) constitute the jth eigenvalue and corresponding row eigenvector of Q(z). The Perron-Fröbenius (P-F) eigenvalue is denoted by λ1(z) (see the appendix for further details). Due to definition (3.17), this can be transformed into(3.21)∑i=0∞y-iφ(Q(z)ix)=∑i=0∞y-i∑l=0N∑m=0NFlm(z)(λ1(z)iw1(z)x)m(λ2(z)iw2(z)x)N-mpl.
Let us now consider values of y and z for which |λj(z)|<|y|≤1,1≤j≤2.  Such values of y and z exist; this is for instance the case if |y|=1, and |z|=1and z≠1 (e.g., see [26]) since(3.22)|λ1(z)|<1,|λ2(z)|<1,for|z|=1,z≠1,
unless some special cases for the generating functions Gij(z) defined in (2.3) are considered—for example, the case that the number of customer arrivals per slot is always a multiple of n,n>1—a situation that falls beyond the scope of this paper.  The sum for i in the above expression then converges, and we obtain
(3.23)∑i=0∞y-iφ(Q(z)ix)=∑m=0Ny(w1(z)x)m(w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)pl.
In a similar way we can derive that
(3.24)∑i=1∞y-iφ(Q(z)ix)=∑m=0N(λ1(z)w1(z)x)m(λ2(z)w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)pl,∑i=1∞y-iψ(z,Q(z)ix)=∑m=0N(λ1(z)w1(z)x)m(λ2(z)w2(z)x)N-my-λ1(z)mλ2(z)N-m∑l=0NFlm(z)ψl(z),
where
(3.25)ψl(z)=̂lim⁡k→∞E[zvk-1{hk=1,ak,1=l,ak,2=N-l}].
In addition, note that
(3.26)lim⁡n→∞y-nP(z,y,Q(z)nx)=∑j=0∞zj∑k=0∞yk∑l=0N∑m=0NFlm(z)Pr[v=j,h=k,a1=l,a2=N-l]⋅lim⁡n→∞y-n(λ1(z)nw1(z)x)m(λ2(z)nw2(z)x)N-m,
which, for |λj(z)|<|y|≤1,1≤j≤2, becomes equal to 0.  Summarizing, in view of the previous results, (3.16) becomes
(3.27)P(z,y,x)=∑m=0N(w1(z)x)m(w2(z)x)N-my-Am(z)×∑l=0NFlm(z){[y(S(y)-z)ψl(z)+(yS(y)-1)pl]Am(z)+y(1-S(y))pl},
where Am(z) has been defined as
(3.28)Am(z)=̂λ1(z)mλ2(z)N-m.
Expression (3.27) for the joint probability generating function P(z,y,x) still contains the unknown constants pl, as well as the unknown functions ψl(z).  Let us therefore, for each value of m, consider values of y for which y=Am(z), which will be represented by ym(z).  Note that if |z|≤1, due to (3.22) this solution for y also satisfies |ym(z)|≤1, and, for these values of z and y, P(z,y,x) is finite.  This means that, if we multiply both hand sides of (3.27) by y-Am(z) and consider y=ym(z), then the left-hand side becomes equal to zero, implying that the same must hold for the right-hand side.  For each value of 0≤m≤N, this leads to the following equation:
(3.29)∑l=0NFlm(z){Am(z)(S(Am(z))-z)ψl(z)+S(Am(z))(Am(z)-1)pl}=0.
If we insert this equation into (3.27), we finally obtain
(3.30)P(z,y,x)=∑m=0N(w1(z)x)m(w2(z)x)N-m⋅{1+zy(Am(z)-1)(S(y)-S(Am(z)))(y-Am(z))(z-S(Am(z)))}∑l=0NFlm(z)pl,
This final result for P(z,y,x) expresses the joint generating function of the state of the arrival process during an arbitrary slot, and the queue contents and residual service time at the beginning of the next slot, in terms of the unknown constants pl.  In the next section, it will become clear how these can be calculated or approximated.  Also, starting from this result, we will establish expressions for the generating functions of quantities such as the queue contents, unfinished work, and message delay, and their corresponding moments and tail distribution.


## 3.4. Calculation of the Unknownspl

We can calculate these probabilities, defined in (3.19), by expressing that P(z,y,x) is analytic when the complex variables y and z both lie inside the unit circle, that is, |z|≤1 and |y|≤1.  We have already expressed in the previous section that, for |y|≤1 and y=Am(z),P(z,y,x) is finite.  In view of expression (3.30), the only remaining potential singularities inside the complex unit disk are those values of z (and y) for which z=S(Am(z)).  Indeed, taking into account the definition for Am(z), then from (3.22) and Rouché’s theorem it follows that this equation has exactly one solution inside the complex unit circle.  Let us, for each value of 0≤m≤N, denote this solution by z{m}, that is,
(3.31)z{m}=S(Am(z{m})),|z{m}|≤1,0≤m≤N.
Then, by expressing that z{m}must also be a zero of the corresponding numerator in expression (3.30) for P(z,y,x), we obtain
(3.32)ξm(z{m})=0,0≤m≤N-1,
where(3.33)ξm(z)=̂11-ρ∑l=0NFlm(z)pl.
Note that z{N}=1 (since λ1(1)=1, see (A.6), implying AN(1)=1), which leads to no additional equation for the unknowns.  Together with the normalization condition, which yields ξN(1)=1 (as we will see further on), (3.32) forms a set of N+1 linear equations for the N+1 unknowns pl that can easily be solved.However, in deriving numerical results for the performance measures considered further in this paper, we want to avoid all numerical calculations as much as possible.  This makes sense in the particular case thatN is large, that is, when a high number of sources are multiplexed, since the number of equations to be solved equals N+1.  Therefore we also present an approximation for the boundary probabilities pl.  As before, denote by the random variables e and v the number of customer arrivals during a tagged arbitrary slot and the queue contents at the beginning of the following slot; also, a represents the state of the Markovian arrival process during the tagged slot.  We define the steady-state joint probability
(3.34)ql=̂Pr[a1=l,a2=N-l,e=0].
Obviously, v=0 implies that there have been no customer arrivals during the tagged slot, that is, v=0⇒e=0.  It is therefore clear that the inequality ql≥pl holds.  We will show by some numerical examples that approximating the boundary probabilities by
(3.35)pl≅1-ρPr[e=0]ql
yields excellent approximations for the performance measures of interest.  Here, we have taken into account that Pr[v=0]=1-ρ so that substituting pl by this formula in (3.30) still yields a normalized generating function.An additional advantage of the approximation proposed here is that it avoids the explicit calculation of theql’s, as well as the functions Flm(z) from (3.18), that occur in expression (3.30) for P(z,y,x).  From definition (3.17) and expression (2.8) it is not difficult to show that
(3.36)∑l=0NFlm(z)ql=(Nm)(σTQ(0)u1(z))m(σTQ(0)u2(z))N-m,
where uj(z) represents the jth column eigenvector of Q(z).  In view of the previous, this implies that the following approximation may be applied:(3.37)ξm(z)≅1Pr[e=0](Nm)(σTQ(0)u1(z))m(σTQ(0)u2(z))N-m.
To conclude, note that it immediately follows from (2.8) that Pr[e=0] satisfies
(3.38)Pr[e=0]=(σTQ(0)1)N.


## 4. The Queue Contents

### 4.1. Derivation of the Probability Generating Function

First of all, due to property (A.2), we have that wj(z)1=1,1≤j≤2.  Therefore, inserting y=1 and x=1 into (3.30), we obtain the following expression for the probability generating function V(z) describing the number of customers in the queue at the beginning of an arbitrary slot:(4.1)V(z)=(1-ρ)∑m=0N(z-1)S(Am(z))z-S(Am(z))ξm(z).
From (A.7) (see the appendix) and (2.15), it follows that S’(1)AN’(1)=NS’(1)λ1’(1)=ρ, implying that the normalization condition for V(z) indeed requires that ξN(1)=1, that is, Pr[v=0]=1-ρ.  From (4.1), one can easily derive the performance measures of interest related to the queue contents, as shown next.


### 4.2. Moments of the Queue Contents

The mean and variance of the queue contentsv can be calculated by taking the appropriate derivatives of expression (4.1) for V(z).  First, note that (A.9) implies that Flm(1)=δ(N-m) (where δ(n) represents the discrete Krönecker-delta function, which equals 1 if n is 0 and 0 otherwise). Since this, in turn, implies that ξm(1)=δ(N-m), and defining ξm(z)=̂S(Am(z)), we obtain the following expression for mean and variance of the queue contents:
(4.2)E[v]=ξN′(1)+ξN′′(1)2(1-ξN′(1))+ξN′(1),Var[v]=ξN′′(1)+ξN′(1)-ξN′(1)2+ξN′′′(1)3(1-ξN′(1))+ξN′′(1)2(1-ξN′(1))(1+ξN′′(1)2(1-ξN′(1)))+ξN′′(1)+ξN′(1)-ξN′(1)2+T′(1),
where, with (A.7) and (2.15),(4.3)ξN′(1)=NpLp=ρ,ξN′′(1)=N2p2S′′(1)+(N-1)ρp+NLpλ1′′(1),ξN′′′(1)=N3p3S′′′(1)+3Np((N-1)Lpλ1′′(1)+NS′′(1)λ1′′(1)+N(N-1)S′′(1)p2)+ρ(N-1)(N-2)p2+NLpλ1′′′(1).
In these expressions, primes denote derivatives with respect to the argument, as before. From definition (3.33) for ξm(z) and the observation that FlN(z)=u11(z)lu12(z)N-l, we also find that
(4.4)ξN′(1)=11-ρ∑l=0N(lu11′(1)+(N-l)u21′(1))pl,ξN′′=11-ρ∑l=0N(lu11′′(1)+2l(N-l)u11′(1)u21′(1)+(N-l)u21′′(1)+l(l-1)u11′(1)2+(N-l)(N-l-1)u21′(1)2)pl.
If we choose to use approximation (3.37) for ξm(z) in order to avoid the numerical calculation of the pl’s, this becomes
(4.5)ξN′(1)≅N(σTQ(0)u1′(1))(σTQ(0)1)-1,ξN′′(1)≅N(σTQ(0)u1′′(1))(σTQ(0)1)-1+N(N-1)(σTQ(0)u1′)2(σTQ(0)1)-2.
In addition, the term T′(1) in the expression for Var[v] is given by
(4.6)T′(1)=(1-ρ)∑m=0N-1Am(1)1-Am(1)ξm′(1),
which can be considered as the first derivative with respect to z of T(z) for z=1, where, using (3.22) and the definition for ξm(z) and Am(z), T(z) satisfies
(4.7)T(z)=11-ρ∑m=0N-1Am(1)1-Am(1)ξm(z)=∑l=0N∑k=1∞{-FlN(z)+∑m=0NAm(1)kFlm(z)}pl=∑l=0N∑k=1∞{-u11(z)lu21(z)N-l+(∑j=12u1j(z)λj(1)k)l(∑j=12u2j(z)λj(1)k)N-l}pl,
where we have also taken into account that λ1(1)=AN(1)=1. Taking the first-order derivative, we thus obtain, keeping in mind that uj2(1)=0 and uj2′(1)=-uj1′(1),1≤j≤2 (see (A.9)-(A.10))
(4.8)T′(1)=-(1-ρ)λ2(1)1-λ2(1)ξN′(1).
By using expressions (4.1)–(4.8), we have now expressed the moments of the queue contents exclusively in terms of the system parameters and the derivatives of the P-F eigenvalue and eigenvector given in (A.7)–(A.11).


### 4.3. Tail Distribution of the Queue Contents

It has been observed in many cases (see, e.g., [16, 20, 26–29]) that approximating the tail distribution of the queue contents by a geometric form may be quite accurate if the dominant singularity (i.e., the singularity with the smallest modulus) of V(z) is a pole; as shown in [20], this will be the case if (4.9)lim⁡z→<RS(λ1(z)N)z>1,
where ℛ is the minimum of the radii of convergence of the pgf’s S(Gij(z)),1≤i,j≤2. If each of  the pgf’s S(·) and Gij(·) either is analytic in the entire complex plane or has a dominant singularity that is a pole, then ℛ corresponds to a pole as well and the above inequality is always satisfied, since S(λ1(z)N) will then tend to +∞ for z→ℛ.  Only if ℛ corresponds to a branch point, it may be so that this inequality is not satisfied; this calls for a case-by-case investigation of the corresponding tail behaviour of the queue contents and falls beyond the scope of the current paper.Approximating the distribution of the queue contents by a geometric form corresponds to approximatingV(z) by(4.10)V(z)≅Cvz-z0,v=-Cvz0,v∑n=0∞(zz0,v)n,
where we are particularly interested in sufficiently large values of n.  In this expression, z0,v is the pole with the smallest modulus of V(z), which is a zero of the denominator for m=N in expression (4.1) for V(z), which yields
(4.11)z0,v=S(AN(z0,v))=S(λ1(z0,v)N),
and z0,v is a real and positive number larger than 1. We assume that this is a pole with multiplicity 1; note that the above formula can be readily extended to the case where z0,v has multiplicity larger than 1 should such a scenario occur. Furthermore, the constant Cv can be calculated from the residue theorem, leading to
(4.12)Cv=lim⁡z→z0,v(z-z0,v)V(z)=(1-ρ)z0,v(z0,v-1)1-Nλ1′(z0,v)λ1(z0,v)N-1S′(λ1(z0,v)N)ξN(z0,v).
Combining these relations, we find the following approximate expression for the probability that the queue contents exceed a certain threshold n:
(4.13)Pr[v>n]≅(1-ρ)ξN(z0,v)1-Nλ1′(z0,v)λ1(z0,v)N-1S′(λ1(z0,v)N)z0,v-n.
Again, we can either use the exact value for ξN(z0,v) calculated from (3.33) or the approximation calculated from (3.37).


## 4.1. Derivation of the Probability Generating Function

First of all, due to property (A.2), we have that wj(z)1=1,1≤j≤2.  Therefore, inserting y=1 and x=1 into (3.30), we obtain the following expression for the probability generating function V(z) describing the number of customers in the queue at the beginning of an arbitrary slot:(4.1)V(z)=(1-ρ)∑m=0N(z-1)S(Am(z))z-S(Am(z))ξm(z).
From (A.7) (see the appendix) and (2.15), it follows that S’(1)AN’(1)=NS’(1)λ1’(1)=ρ, implying that the normalization condition for V(z) indeed requires that ξN(1)=1, that is, Pr[v=0]=1-ρ.  From (4.1), one can easily derive the performance measures of interest related to the queue contents, as shown next.


## 4.2. Moments of the Queue Contents

The mean and variance of the queue contentsv can be calculated by taking the appropriate derivatives of expression (4.1) for V(z).  First, note that (A.9) implies that Flm(1)=δ(N-m) (where δ(n) represents the discrete Krönecker-delta function, which equals 1 if n is 0 and 0 otherwise). Since this, in turn, implies that ξm(1)=δ(N-m), and defining ξm(z)=̂S(Am(z)), we obtain the following expression for mean and variance of the queue contents:
(4.2)E[v]=ξN′(1)+ξN′′(1)2(1-ξN′(1))+ξN′(1),Var[v]=ξN′′(1)+ξN′(1)-ξN′(1)2+ξN′′′(1)3(1-ξN′(1))+ξN′′(1)2(1-ξN′(1))(1+ξN′′(1)2(1-ξN′(1)))+ξN′′(1)+ξN′(1)-ξN′(1)2+T′(1),
where, with (A.7) and (2.15),(4.3)ξN′(1)=NpLp=ρ,ξN′′(1)=N2p2S′′(1)+(N-1)ρp+NLpλ1′′(1),ξN′′′(1)=N3p3S′′′(1)+3Np((N-1)Lpλ1′′(1)+NS′′(1)λ1′′(1)+N(N-1)S′′(1)p2)+ρ(N-1)(N-2)p2+NLpλ1′′′(1).
In these expressions, primes denote derivatives with respect to the argument, as before. From definition (3.33) for ξm(z) and the observation that FlN(z)=u11(z)lu12(z)N-l, we also find that
(4.4)ξN′(1)=11-ρ∑l=0N(lu11′(1)+(N-l)u21′(1))pl,ξN′′=11-ρ∑l=0N(lu11′′(1)+2l(N-l)u11′(1)u21′(1)+(N-l)u21′′(1)+l(l-1)u11′(1)2+(N-l)(N-l-1)u21′(1)2)pl.
If we choose to use approximation (3.37) for ξm(z) in order to avoid the numerical calculation of the pl’s, this becomes
(4.5)ξN′(1)≅N(σTQ(0)u1′(1))(σTQ(0)1)-1,ξN′′(1)≅N(σTQ(0)u1′′(1))(σTQ(0)1)-1+N(N-1)(σTQ(0)u1′)2(σTQ(0)1)-2.
In addition, the term T′(1) in the expression for Var[v] is given by
(4.6)T′(1)=(1-ρ)∑m=0N-1Am(1)1-Am(1)ξm′(1),
which can be considered as the first derivative with respect to z of T(z) for z=1, where, using (3.22) and the definition for ξm(z) and Am(z), T(z) satisfies
(4.7)T(z)=11-ρ∑m=0N-1Am(1)1-Am(1)ξm(z)=∑l=0N∑k=1∞{-FlN(z)+∑m=0NAm(1)kFlm(z)}pl=∑l=0N∑k=1∞{-u11(z)lu21(z)N-l+(∑j=12u1j(z)λj(1)k)l(∑j=12u2j(z)λj(1)k)N-l}pl,
where we have also taken into account that λ1(1)=AN(1)=1. Taking the first-order derivative, we thus obtain, keeping in mind that uj2(1)=0 and uj2′(1)=-uj1′(1),1≤j≤2 (see (A.9)-(A.10))
(4.8)T′(1)=-(1-ρ)λ2(1)1-λ2(1)ξN′(1).
By using expressions (4.1)–(4.8), we have now expressed the moments of the queue contents exclusively in terms of the system parameters and the derivatives of the P-F eigenvalue and eigenvector given in (A.7)–(A.11).


## 4.3. Tail Distribution of the Queue Contents

It has been observed in many cases (see, e.g., [16, 20, 26–29]) that approximating the tail distribution of the queue contents by a geometric form may be quite accurate if the dominant singularity (i.e., the singularity with the smallest modulus) of V(z) is a pole; as shown in [20], this will be the case if (4.9)lim⁡z→<RS(λ1(z)N)z>1,
where ℛ is the minimum of the radii of convergence of the pgf’s S(Gij(z)),1≤i,j≤2. If each of  the pgf’s S(·) and Gij(·) either is analytic in the entire complex plane or has a dominant singularity that is a pole, then ℛ corresponds to a pole as well and the above inequality is always satisfied, since S(λ1(z)N) will then tend to +∞ for z→ℛ.  Only if ℛ corresponds to a branch point, it may be so that this inequality is not satisfied; this calls for a case-by-case investigation of the corresponding tail behaviour of the queue contents and falls beyond the scope of the current paper.Approximating the distribution of the queue contents by a geometric form corresponds to approximatingV(z) by(4.10)V(z)≅Cvz-z0,v=-Cvz0,v∑n=0∞(zz0,v)n,
where we are particularly interested in sufficiently large values of n.  In this expression, z0,v is the pole with the smallest modulus of V(z), which is a zero of the denominator for m=N in expression (4.1) for V(z), which yields
(4.11)z0,v=S(AN(z0,v))=S(λ1(z0,v)N),
and z0,v is a real and positive number larger than 1. We assume that this is a pole with multiplicity 1; note that the above formula can be readily extended to the case where z0,v has multiplicity larger than 1 should such a scenario occur. Furthermore, the constant Cv can be calculated from the residue theorem, leading to
(4.12)Cv=lim⁡z→z0,v(z-z0,v)V(z)=(1-ρ)z0,v(z0,v-1)1-Nλ1′(z0,v)λ1(z0,v)N-1S′(λ1(z0,v)N)ξN(z0,v).
Combining these relations, we find the following approximate expression for the probability that the queue contents exceed a certain threshold n:
(4.13)Pr[v>n]≅(1-ρ)ξN(z0,v)1-Nλ1′(z0,v)λ1(z0,v)N-1S′(λ1(z0,v)N)z0,v-n.
Again, we can either use the exact value for ξN(z0,v) calculated from (3.33) or the approximation calculated from (3.37).


## 5. The Unfinished Work

### 5.1. Derivation of the Probability Generating Function

The unfinished work in the queue at the beginning of a slot is defined as the amount of work in the queue at the beginning of the slot, which is the number of slots required to empty the queue if no more customers were to arrive during the subsequent slots.  The amount of work required to process a customer whose service has not started yet is given by its service time described byS(z).  On the other hand, at the beginning of a slot, the customer currently being served (if any) still requires h slots before it is entirely sent, where h equals the residual service time.  The unfinished work uk at the beginning of slot k+1 is therefore given by
(5.1)uk={0,ifvk=0,hk=0,hk+∑i=1vk-1siifvk>0,
where each of the si’s represents a customer service time, and is therefore described by the pgf S(z).   Consequently, the corresponding steady-state pgf U(z) is given by
(5.2)U(z)=S(z)-1{P(S(z),z,1)+(S(z)-1)P(0,0,1)}.
The observation that P(0,0,1)=Pr[v=0]=1-ρcan be written as
(5.3)P(0,0,1)=∑l=0Npl=∑m=0N∑l=0NFlm(S(z))pl,
together with the previous relations, leads to the following expression for the steady-state probability generating function describing the unfinished work at the beginning of an arbitrary slot:
(5.4)U(z)=(1-ρ)∑m=0N(z-1)Am(S(z))z-Am(S(z))ξm(S(z)).
Again, due to A′(1)S′(1)=NSp=ρ and ξN(1)=δ(N-m), note that U(1)=1, as expected.  We would like to point out that each of the denominators in the above expression also has a zero inside the unit disk |z|≤1, which leads to  an additional linear equation for the boundary probabilities pl (since U(z) is bounded inside the complex unit circle |z|<1), as explained in Section 3.4.  However, note that if z̃{m} is a solution of z̃{m}=Am(S(z̃{m})) with |z̃{m}|≤1, we then have that S(z̃{m})=S(Am(S(z̃{m}))), implying that(5.5)z{m}=S(z̃{m}),z̃{m}=Am(z{m}),
where z{m} satisfies (3.35), and we obtain exactly the same set of linear equations as in (3.37) for the boundary probabilities.


### 5.2. Moments of the Unfinished Work

In a similar way to that in Section4.2, if we define χm(z)=̂Am(S(z)) and take the appropriate derivatives with respect to z of expression (5.4) for U(z), we find an expression for the mean and variance of the unfinished work:
(5.6)E[u]=χN′(1)+χN′(1)2(1-χN′(1))+LpξN′(1),Var[u]=χN′′(1)+χN′(1)-χN′(1)2+χN′′′(1)3(1-χN′(1))+χN′′(1)2(1-χN′(1))(1+χN′′(1)2(1-χN′(1)))+Lp2ξN′′(1)+S′′(1)ξN′(1)+LpξN′(1)-Lp2ξN′(1)2+2LpT′(1),
where
(5.7)χN′(1)=NpLp=ρ,χN′′(1)=ρ2(1-1N)+ρ(λ1′′(1)Lpp+S′′(1)Lp),χN′′′(1)=ρ3(1-1N)(1-2N)+3ρ2(1-1N)(λ1′′(1)Lpp+S′′(1)Lp)+ρ(λ1′′′(1)Lp2p+3λ1′′(1)S′′(1)p+S′′′(1)Lp).
The derivatives of ξN(z) and T′(1) in the previous expressions are given by (4.4) (or (4.5) in case the approximations for the boundary probabilities are applied) and (4.8), respectively.


### 5.3. Tail Distribution of the Unfinished Work

Under similar assumptions to those outlined in Section4.3, we can adopt a geometric tail approximation for the distribution of the unfinished work.  We then obtain an approximation for the probability that the unfinished work exceeds a threshold n:
(5.8)Pr[u>n]≅(1-ρ)ξN(S(z0,u))1-NS′(z0,u)λ1′(S(z0,u))λ1(S(z0,u))N-1z0,u-n,
where the smallest pole z0,u of U(z) is the solution of z0,u=λ1(S(z0,u))N, which is a real and positive quantity larger than 1, and, similar to (5.5), can also be calculated from z0,u=λ1(z0,v)N.  Also note that S(z0,u)=z0,v, implying that the above formula can be expressed in terms of the smallest pole z0,v of V(z) as well.


## 5.1. Derivation of the Probability Generating Function

The unfinished work in the queue at the beginning of a slot is defined as the amount of work in the queue at the beginning of the slot, which is the number of slots required to empty the queue if no more customers were to arrive during the subsequent slots.  The amount of work required to process a customer whose service has not started yet is given by its service time described byS(z).  On the other hand, at the beginning of a slot, the customer currently being served (if any) still requires h slots before it is entirely sent, where h equals the residual service time.  The unfinished work uk at the beginning of slot k+1 is therefore given by
(5.1)uk={0,ifvk=0,hk=0,hk+∑i=1vk-1siifvk>0,
where each of the si’s represents a customer service time, and is therefore described by the pgf S(z).   Consequently, the corresponding steady-state pgf U(z) is given by
(5.2)U(z)=S(z)-1{P(S(z),z,1)+(S(z)-1)P(0,0,1)}.
The observation that P(0,0,1)=Pr[v=0]=1-ρcan be written as
(5.3)P(0,0,1)=∑l=0Npl=∑m=0N∑l=0NFlm(S(z))pl,
together with the previous relations, leads to the following expression for the steady-state probability generating function describing the unfinished work at the beginning of an arbitrary slot:
(5.4)U(z)=(1-ρ)∑m=0N(z-1)Am(S(z))z-Am(S(z))ξm(S(z)).
Again, due to A′(1)S′(1)=NSp=ρ and ξN(1)=δ(N-m), note that U(1)=1, as expected.  We would like to point out that each of the denominators in the above expression also has a zero inside the unit disk |z|≤1, which leads to  an additional linear equation for the boundary probabilities pl (since U(z) is bounded inside the complex unit circle |z|<1), as explained in Section 3.4.  However, note that if z̃{m} is a solution of z̃{m}=Am(S(z̃{m})) with |z̃{m}|≤1, we then have that S(z̃{m})=S(Am(S(z̃{m}))), implying that(5.5)z{m}=S(z̃{m}),z̃{m}=Am(z{m}),
where z{m} satisfies (3.35), and we obtain exactly the same set of linear equations as in (3.37) for the boundary probabilities.


## 5.2. Moments of the Unfinished Work

In a similar way to that in Section4.2, if we define χm(z)=̂Am(S(z)) and take the appropriate derivatives with respect to z of expression (5.4) for U(z), we find an expression for the mean and variance of the unfinished work:
(5.6)E[u]=χN′(1)+χN′(1)2(1-χN′(1))+LpξN′(1),Var[u]=χN′′(1)+χN′(1)-χN′(1)2+χN′′′(1)3(1-χN′(1))+χN′′(1)2(1-χN′(1))(1+χN′′(1)2(1-χN′(1)))+Lp2ξN′′(1)+S′′(1)ξN′(1)+LpξN′(1)-Lp2ξN′(1)2+2LpT′(1),
where
(5.7)χN′(1)=NpLp=ρ,χN′′(1)=ρ2(1-1N)+ρ(λ1′′(1)Lpp+S′′(1)Lp),χN′′′(1)=ρ3(1-1N)(1-2N)+3ρ2(1-1N)(λ1′′(1)Lpp+S′′(1)Lp)+ρ(λ1′′′(1)Lp2p+3λ1′′(1)S′′(1)p+S′′′(1)Lp).
The derivatives of ξN(z) and T′(1) in the previous expressions are given by (4.4) (or (4.5) in case the approximations for the boundary probabilities are applied) and (4.8), respectively.


## 5.3. Tail Distribution of the Unfinished Work

Under similar assumptions to those outlined in Section4.3, we can adopt a geometric tail approximation for the distribution of the unfinished work.  We then obtain an approximation for the probability that the unfinished work exceeds a threshold n:
(5.8)Pr[u>n]≅(1-ρ)ξN(S(z0,u))1-NS′(z0,u)λ1′(S(z0,u))λ1(S(z0,u))N-1z0,u-n,
where the smallest pole z0,u of U(z) is the solution of z0,u=λ1(S(z0,u))N, which is a real and positive quantity larger than 1, and, similar to (5.5), can also be calculated from z0,u=λ1(z0,v)N.  Also note that S(z0,u)=z0,v, implying that the above formula can be expressed in terms of the smallest pole z0,v of V(z) as well.


## 6. The Customer Delay

### 6.1. Derivation of the Probability Generating Function

We define the customer delayd as the number of slots an arbitrary customer remains in the queue, including its service time, that is, the number of slots between the end of its arrival and departure slot, as depicted in Figure 1.  Consider an arbitrary tagged customer that enters the queue during a slot, say slot I. Note that this slot is not an arbitrary slot but a slot where an arbitrary customer arrives (which for instance implies that at least one customer arrives).  Due to the FCFS service discipline, the delay of a customer is determined by the amount of work in the queue upon the tagged customer’s arrival.  In particular, if we denote by ûI the unfinished work at the beginning of slot I and by fI the number of customers that have arrived during the same slot as the tagged one and will be served before it, then we find that the delay of the tagged customer is given by(6.1)d=(ûI-1)++∑i=1fI+1si,
where the random variables {si,1≤i≤fI+1} represent the service times of the tagged customer and the fI customers that have arrived during its arrival slot and will be served before it, which are all described by the pgf S(z). Let us denote by G(z,x) the steady-state joint generating function of the random variables (u,e) describing the unfinished work at the beginning and the number of customer arrivals during an arbitrary slot k, that is,
(6.2)G(z,x)=E[zuxe],
and by H(z,x) the steady state joint pgf of (ûI,fk).  Then, similar to [30], we can establish the following relation:
(6.3)H(z,x)=G(z,1)-G(z,x)Np(1-x).
System equation (6.1) can be transformed into a relation between z-transforms, which together with (6.3), leads to an expression for the steady-state pgf D(z) describing the delay of an arbitrary customer
(6.4)D(z)=S(z)z{H(z,S(z))+(z-1)H(0,S(z))}=S(z)Npz{G(z,1)-G(z,S(z))1-S(z)+(z-1)G(0,1)-G(0,S(z))1-S(z)}.
On the other hand, since the unfinished work satisfies the system equation
(6.5)uk=(uk-1-1)++∑i=1eksi
and assuming that a stochastic equilibrium is reached, this implies that the pgf U(z) (given by (5.4)) can also be written as (again, all si’s are described by S(z))
(6.6)U(z)=z-1{G(z,S(z))+(z-1)G(0,S(z))}.
In combination with the previous expression for D(z) (and making use of U(z)=G(z,1)), this finally leads to
(6.7)D(z)=S(z)zS′(1)(z-1)S(z)-1U0(z),
where in view of U(0)=Pr[u=0]=Pr[v=0]=1-ρ(6.8)U0(z)=E[zu∣u>0]=(U(z)-U(0))ρ.
This formula for D(z), which expresses the delay of an arbitrary customer in terms of the unfinished work at the beginning of an arbitrary slot and which apparently is independent of the details of the customer arrival process (correlation, etc.) was also derived in [31] via an alternative analysis.  From this result, it is not difficult to express the performance measures concerning the customer delay in terms of the corresponding performance measures for the unfinished work, as shown next.Figure 1
The delay of an arbitrary customer, with related quantities.


### 6.2. Moments of the Customer Delay

By taking the first two order derivatives of (6.7)-(6.8) with respect to z for z=1, we can establish the following relations between the mean and variance of the customer delay and the mean and variance of the unfinished work:
(6.9)E[d]=Lp-1-S′′(1)2Lp+E[u]ρ,Var[d]=Var[s]-S′′′(1)3Lp-S′′(1)2Lp(1-S′′(1)2Lp)+Var[u]ρ-(1-ρ)(E[u]ρ)2.
Through some simple calculations, one can verify that Little’s theorem indeed holds, that is, E[d]=E[v]/(Np).


### 6.3. Tail Distribution of the Customer Delay

Again, applying an analogous technique as the one used for determining the tail distribution of the queue contents and unfinished work, the asymptotic behavior of the customer delay tail distribution can be calculated from(6.10)D(z)≅Cdz-z0,d=-Cdz0,d∑n=0∞(zz0,d)n,Cd=̂lim⁡z→z0,d(z-z0,d)D(z),
where z0,d is the smallest pole of D(z) which, in view of (6.7)-(6.8), leads to the observation that z0,d=z0,u, and the probability that the customer delay exceeds a threshold n is approximated by(6.11)Pr[d>n]≅1-ρρLp(z0,u-1)S(z0,u)-1S(z0,u)ξN(S(z0,u))1-Nλ1′(S(z0,u))S′(z0,u)λ1(S(z0,u))N-1z0,u-n-1.


## 6.1. Derivation of the Probability Generating Function

We define the customer delayd as the number of slots an arbitrary customer remains in the queue, including its service time, that is, the number of slots between the end of its arrival and departure slot, as depicted in Figure 1.  Consider an arbitrary tagged customer that enters the queue during a slot, say slot I. Note that this slot is not an arbitrary slot but a slot where an arbitrary customer arrives (which for instance implies that at least one customer arrives).  Due to the FCFS service discipline, the delay of a customer is determined by the amount of work in the queue upon the tagged customer’s arrival.  In particular, if we denote by ûI the unfinished work at the beginning of slot I and by fI the number of customers that have arrived during the same slot as the tagged one and will be served before it, then we find that the delay of the tagged customer is given by(6.1)d=(ûI-1)++∑i=1fI+1si,
where the random variables {si,1≤i≤fI+1} represent the service times of the tagged customer and the fI customers that have arrived during its arrival slot and will be served before it, which are all described by the pgf S(z). Let us denote by G(z,x) the steady-state joint generating function of the random variables (u,e) describing the unfinished work at the beginning and the number of customer arrivals during an arbitrary slot k, that is,
(6.2)G(z,x)=E[zuxe],
and by H(z,x) the steady state joint pgf of (ûI,fk).  Then, similar to [30], we can establish the following relation:
(6.3)H(z,x)=G(z,1)-G(z,x)Np(1-x).
System equation (6.1) can be transformed into a relation between z-transforms, which together with (6.3), leads to an expression for the steady-state pgf D(z) describing the delay of an arbitrary customer
(6.4)D(z)=S(z)z{H(z,S(z))+(z-1)H(0,S(z))}=S(z)Npz{G(z,1)-G(z,S(z))1-S(z)+(z-1)G(0,1)-G(0,S(z))1-S(z)}.
On the other hand, since the unfinished work satisfies the system equation
(6.5)uk=(uk-1-1)++∑i=1eksi
and assuming that a stochastic equilibrium is reached, this implies that the pgf U(z) (given by (5.4)) can also be written as (again, all si’s are described by S(z))
(6.6)U(z)=z-1{G(z,S(z))+(z-1)G(0,S(z))}.
In combination with the previous expression for D(z) (and making use of U(z)=G(z,1)), this finally leads to
(6.7)D(z)=S(z)zS′(1)(z-1)S(z)-1U0(z),
where in view of U(0)=Pr[u=0]=Pr[v=0]=1-ρ(6.8)U0(z)=E[zu∣u>0]=(U(z)-U(0))ρ.
This formula for D(z), which expresses the delay of an arbitrary customer in terms of the unfinished work at the beginning of an arbitrary slot and which apparently is independent of the details of the customer arrival process (correlation, etc.) was also derived in [31] via an alternative analysis.  From this result, it is not difficult to express the performance measures concerning the customer delay in terms of the corresponding performance measures for the unfinished work, as shown next.Figure 1
The delay of an arbitrary customer, with related quantities.


## 6.2. Moments of the Customer Delay

By taking the first two order derivatives of (6.7)-(6.8) with respect to z for z=1, we can establish the following relations between the mean and variance of the customer delay and the mean and variance of the unfinished work:
(6.9)E[d]=Lp-1-S′′(1)2Lp+E[u]ρ,Var[d]=Var[s]-S′′′(1)3Lp-S′′(1)2Lp(1-S′′(1)2Lp)+Var[u]ρ-(1-ρ)(E[u]ρ)2.
Through some simple calculations, one can verify that Little’s theorem indeed holds, that is, E[d]=E[v]/(Np).


## 6.3. Tail Distribution of the Customer Delay

Again, applying an analogous technique as the one used for determining the tail distribution of the queue contents and unfinished work, the asymptotic behavior of the customer delay tail distribution can be calculated from(6.10)D(z)≅Cdz-z0,d=-Cdz0,d∑n=0∞(zz0,d)n,Cd=̂lim⁡z→z0,d(z-z0,d)D(z),
where z0,d is the smallest pole of D(z) which, in view of (6.7)-(6.8), leads to the observation that z0,d=z0,u, and the probability that the customer delay exceeds a threshold n is approximated by(6.11)Pr[d>n]≅1-ρρLp(z0,u-1)S(z0,u)-1S(z0,u)ξN(S(z0,u))1-Nλ1′(S(z0,u))S′(z0,u)λ1(S(z0,u))N-1z0,u-n-1.


## 7. Some Numerical Examples

We will illustrate the results derived throughout the previous sections by a few numerical examples.  Let us therefore consider a two-state ON/OFF MMBP arrival process for each source, which means that(7.1)Gj1(z)=1-γ+γz,Gj2(z)=1,1≤j≤2,
that is, customers are generated at a rate γ during an ON period (corresponding with state S1 in the modulating Markov chain) and no customers are generated during an OFF period (corresponding with state S2).  Let us also denote by L1 and L2 the mean length of an ON and OFF period, respectively, that is, L1=1/p12 and L2=1/p21. Furthermore, we will assume that customers have a constant service time equal to Lp slots, implying that S(z)=zLp. The load of the system is then given by
(7.2)ρ=NLpγL1L1+L2,
and any subset containing five parameters of the set (N,ρ,γ,L1,L2,Lp) fully describes the customer arrival process in the queue.  Such a model has for instance been used to effectively model a voice-over-IP (VoIP) packet arrival stream in a network node [3].In the remainder we consider a multiplexer queue fed byN=16 identical sources and focus on the customer delay behaviour, which is the primary performance measure for VoIP applications.  Similar results can be readily obtained for the queue contents and unfinished work. For the arrival model described above, we have plotted the mean and variance of the the moments of the customer waiting time w (which is defined as the number of slots a customer waits in the queue before its service commences and is therefore equal to the customer delay minus the customer service time) versus ρ in case of L1=L2=8 and for increasing customer service time Lp=1,2,4,8,16 in Figures 2 and 4, respectively.  Varying the value of the load in this case means that the arrival rate γ is adjusted.  From these figures we can conclude that the moments of the customer waiting time strongly depend on the customer service time and increase for increasing values of Lp.  In these figures, we have also included the curves when approximation (3.35) is used for the boundary probabilities (dashed lines), which leads to approximation (4.5) for the derivatives of ξN(z); the absolute error between the exact and approximate results is so small that the latter are hardly visible in Figures 2 and 4.  Therefore, in Figures 3 and 5, we have also plotted the relative error (RE) that is made when using the approximate formula for the mean and variance of the waiting time.  Apparently, except for Lp=1, the RE of mean and variance appears to be a decreasing function of ρ that approaches zero for ρ=1.  In addition, the larger the customer service times are, the smaller the relative error generally tends to be.  We also observed that, while for the mean value the approximation always forms an upper bound compared to the exact value, this is not the case as far as the variance is concerned.Figure 2
E[w] versus ρ,N=16,L1=L2=8.Figure 3
RE of E[w] versusρ,N=16,L1=L2=8.Figure 4
Var[w] versus ρ,N=16,L1=L2=8.Figure 5
RE of Var[w] versus ρ,N=16,L1=L2=8.In Figure6, we consider constant values for γ=0.05 and Lp=10 and increasing values of the mean length L1 of an active period.  Hence, the probability L1/(L1+L2) that a source is in the ON state during an arbitrary slot tends to 0 for ρ→0.  Clearly, the mean value of the customer waiting time, which is plotted versus ρ in these figures, strongly depends on the mean length of an active period.  In addition, we observe from these curves that the error that is made by using the approximation for the boundary probabilities (dashed line) can be significant when the load is low and the arrival rate γ and L1 are high.  Another feature is that the mean waiting time not necessarily tends to 0 for ρ→0.  Indeed, even under these conditions, we sporadically encounter the situation where a source becomes active and generates customers at a rate γ during a geometrically distributed ON period with mean length L1, and the intersection with the ordinate ρ=0 designates the mean waiting time for customers generated by such a single source.Figure 6
E[w] versus ρ,N=16,γ=0.05,Lp=10.In Figures7, 8 and 9, we have plotted the tail distribution of the customer waiting time for a wide range of system parameters.  Each curve for the exact value of the parameters of the geometric tail behaviour (full line) is accompanied by the corresponding approximation (dashed line) obtained by using approximation (3.37) for ξN(z).  We may conclude that, for the wide range of system parameters considered in these figures, the approximate results are extremely accurate.  Furthermore, it is observed that the approximation forms an upper bound, and, consequently, is on the safe side.  The only numerical calculation required for this upper bound is that of the dominant pole z0,u of the pgf describing the steady-state waiting time distribution.  Insertion of the system parameters in the corresponding formulas then immediately yields these results.Figure 7
Pr[w>n] versus n,N=16,ρ=0.8,γ=0.1,L1=50.Figure 8
Pr[w>n] versus n,N=16,γ=0.1,Lp=10,L1=50.Figure 9
Pr[w>n] versus n,N=16,ρ=0.8,γ=0.1,Lp=10.In Figure7, we have plotted the probability that the customer waiting time exceeds a certain threshold n versus n, for a multiplexer queue with N=16 sources, a system load ρ=0.8, an arrival rate γ=0.1 during an ON period, and a mean length of an ON period L1=50, for varying values of the customer service time Lp.  Again, as was also observed when discussing the moments, for these parameter values the results for the quantile of the waiting time (defined as the smallest value of n(x) for which Pr[w>n(x)]<10-x) significantly depend on the customer service time. From this figure, it also becomes apparent that the quantile of the waiting time is roughly proportional to Lp, and a wide range of case studies has revealed that this observation is not limited to the case of constant customer service times.In Figure8, we show the results for the customer waiting time in case of N=16, γ=0.1, L1=50 and Lp=10 and increasing values of the load (which means that the arrival rate γ during the ON periods increases), while the load is fixed to ρ=0.8 and L1=20,50,100 in Figure 9.   As can be expected, the slope of the customer waiting time tail distribution increases as both ρ and L1 increase. We also should point out that, whereas the accuracy of the approximate results for the moments of the customer waiting time—calculated by invoking the approximation for the boundary probabilities—deteriorates for increasing values of L1 (and L2) (e.g., Figure 6), this is not so as far as the tail distribution of the customer waiting time is concerned, and the approximate data that we thus compute remain highly accurate.


## 8. Conclusion

In this paper we have investigated and analysed the discrete-time (2-state) D-BMAP/G/1 queueing model. By combining a solution method that was previously developed for the D-BMAP/D/c queueing system, with the supplementary variable technique that allows us to cope with generally distributed customer service times, we were able to derive an expression for the steady-state pgfs of the queue contents, unfinished work, and customer delay.  This has led to (semi)closed-form expressions for a number of interesting performance measures such as the mean value, variance, and tail distribution of these quantities. In addition to the analysis of the system, a lot of emphasis was put on finding closed-form expressions for these quantities that reduce all numerical calculations to an absolute minimum, at the expense of introducing approximations for a number of boundary probabilities that must normally be calculated by solving a—potentially large—set of linear equations.  The only numerical computation that needs to be carried out is that of a dominant pole, situated on the positive real axis, when one calculates tail probabilities, and the complexity does not depend on the size of the state space of the model under consideration. To conclude, the accuracy of this approach, as well as the impact of the system parameters on the customer delay performance, was evaluated through a considerable set of numerical examples.



---

*Source: 401365-2011-09-26.xml*