Blocking probability computation in reversible Markovian bufferless multi-server systems

Performance Evaluation 67 (2010) 121–140

Contents lists available at ScienceDirect

Performance Evaluation journal homepage: www.elsevier.com/locate/peva

Blocking probability computation in reversible Markovian bufferless multi-server systems Miguel de Vega Rodrigo a,∗ , Guy Latouche a , Marie-Ange Remiche b a

Université Libre de Bruxelles, SMG - CP 210/01, Blvd du Triomphe, 1050 Bruxelles, Belgium

b

FUNDP - The University of Namur, Faculty of Computer Science, Rue Grandgagnage 21, B-5000-Namur, Belgium

article

info

Article history: Received 14 February 2008 Received in revised form 9 March 2009 Accepted 25 September 2009 Available online 8 October 2009 Keywords: Multi-server bufferless system Blocking probability Long-range dependency Self-similarity Quasi-birth-and-death processes

abstract We study the blocking probability in a Markovian bufferless queuing system with a finite number of servers when the packet arrival process is a Markov-modulated Poisson process (MMPP), a superposition of a finite number N of independent MMPPs. We present an algorithm to compute the blocking probability in the system. The algorithm is exact if some related Quasi-birth-and-death (QBD) processes are reversible. We prove that the complexity of our algorithm scales linearly with N, whereas that of the standard solution does it exponentially. We illustrate the use of our algorithm in a bufferless multi-server system receiving longrange dependent (LRD) input traffic. This problem finds applications in the study of the blocking probability in bufferless optical packet switching (OPS) and optical burst switching (OBS) networks. We emulate LRD traffic with an MMPP, a superposition of N independent ON/OFF sources, which are also MMPPs. Our algorithm is in this case approximative since the reversibility condition is not fulfilled with the proposed MMPP. An extensive numerical analysis suggests that our algorithm can accurately approximate the blocking probability in the original queueing system with LRD input traffic. We provide an insight into the reason why our approximation is accurate. © 2009 Elsevier B.V. All rights reserved.

In this paper we study a queueing system with K < ∞ servers and no buffering capability. We focus on the case in which customers arrive at the system according to a Markov-modulated Poisson process (MMPP), a superposition of a finite number N of independent MMPPs, each one representing a source. If an arriving customer finds all K servers in the system busy, then it is blocked (i.e., lost). Otherwise, an available server is occupied by the customer during a phase-type (PH)distributed period of time, where the parameters of the PH-distribution may depend on the source the customer comes from. Our objective is to develop an analytical method for the computation of the main performance parameter in the bufferless multi-server system described above: the blocking probability. The standard solution to this problem presents serious scalability issues with the number N of traffic sources in the queueing system. More specifically, its complexity scales with N exponentially. We present an analytical method that does not exhibit such scalability problems, since its complexity increases with N only linearly. Let us denote MMPPi and PHi to the arrival process and service time distribution of the i-th traffic source in the queueing system, with 1 ≤ i ≤ N. The analytical method presented in this paper is exact if each QBD i describing the number of busy servers in the MMPPi /PHi /K /K system is reversible. This condition turns out to be rather restrictive, and so our method may be considered in the praxis to be approximative. In [1–5], long-range dependent (LRD) traffic is emulated with the superposition of a finite number of independent ON/OFF sources, which are also MMPPs. Such arrival processes naturally fit in the description of our queueing system. Thus, the

∗

Corresponding author. Tel.: +32 27352303. E-mail addresses: [email protected] (M. de Vega Rodrigo), [email protected] (M.-A. Remiche).

0166-5316/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.peva.2009.09.008

122

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

queueing system addressed in this paper can be used in order to approximate the behavior of a bufferless multi-server queueing system receiving LRD input traffic. For this reason, our analytical method finds applications in the approximative computation of the blocking probability in bufferless multi-server systems with LRD input traffic. Such queueing systems are of practical importance for two reasons. First, because of the ubiquitous presence of LRD in modern telecommunication networks. For instance, numerous empirical studies suggest the presence of LRD in future optical packet switching (OPS) and optical burst switching (OBS) networks [6–10]. Second, because of the significant negative impact of LRD traffic on network performance [11,12]. In order to illustrate the applicability of our method, we study the blocking probability in a bufferless OBS/OPS network. This problem has received much attention in the literature [13–20]. To our knowledge, our study is the first that approximates the blocking probability in the presence of LRD traffic. The paper is structured as follows. In Section 1 the basic definitions and notation are introduced. Section 2 presents the problem under study. In Section 2.1 we present the standard approach for the computation of the blocking probability in the bufferless multi-server system. The new method that we propose is based on the combined use of two concepts. The first concept is the superposition of QBD processes, which is formally defined in Section 3. The second concept is called the simplified Birth-Death (BD) process of a QBD process, and is addressed in Section 4. In Section 5 we present the new method for the computation of the blocking probability, together with a comparison of the complexities of the standard and new approaches in Section 5.1. Section 6 contains a numerical study which empirically evaluates the accuracy of our analytical method in a particular MMPP /M /K /K system that models an optical transmission link in a bufferless OPS/OBS network. In this system the MMPP results from the superposition of a finite number of independent MMPPs, defined according to [3]. As expected, the MMPPi /M /K /K systems for 1 ≤ i ≤ N are not reversible, and thus our analytic method provides an approximative solution. We evaluate the accuracy of such an approximation and compare also the accuracy of our analytical method when it is used to compute the blocking probability in a bufferless multi-server system receiving real LRD traffic. The conclusions of the paper are presented in Section 7. 1. Definitions and notation In this section we briefly define the main mathematical concepts used in the paper. We begin with the definition of Markov modulated Poisson processes (MMPP). We say that {(N (t ), J (t )), t ∈ R+ } is a MMPP with representation MMPP (D0 , D1 ) if [21]: – {(N (t ), J (t )), t ∈ R+ } is a bidimensional Markov process defined on N × {1, . . . , m}, with m ∈ N finite. – The corresponding generator is



D0

D1 D0

 

 D1

..

.

..

.

 ,

(1)

where D1 is a diagonal matrix. As a consequence of the MMPP being Markovian we have that the elements of D1 and the off diagonal elements of D0 are

≥0, the diagonal elements of D0 are <0, and (D0 + D1 )1 = 0, where the symbols 0 and 1 represent hereinafter a column vector of zeros and ones, respectively. The matrix D0 + D1 corresponds to the infinitesimal generator of the phase process {J (t ), t ∈ R+ }. Usually D0 + D1 is irreducible, in which case the stationary distribution ε associated to the MMPP is the unique nonnegative solution of ε(D0 + D1 ) = 0 and ε 1 = 1. The MMPP is stationary if α = ε , where α is the distribution of J (0). We assume hereinafter that all MMPPs are irreducible and stationary. We proceed now with the definition of phase-type (PH) distributions. Let {L(t ), t ∈ R+ } be an absorbing Markov process with generator



0 t

0T T



,

(2)

and initial phase distribution (0τ ), defined on the state space {0, 1, . . . , n} with n ∈ N finite, where 0T and τ are row vectors of size n, and t is a column vector of size n. The state 0 is the absorbing state. The time to absorption Z is a continuous phasetype random variable of order n represented as PH (τ , T ) with P [Z ≤ z ] = 1 − τ exp(Tz )1,

(3)

and z ∈ R+ [21]. A finite nonhomogeneous continuous-time Birth-death (BD) process {X (t ), t ∈ R+ } is a Markov process defined on the set of nonnegative integers {0, . . . , K }, with K < ∞ and an infinitesimal generator of the form

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140



(0)



(1)

a1

a 2 q=  

(0)

a0

(1)

a1

..

 .. ..

.

123

. .

(K )

a2

  ,  (K −1)  a 0

(4)

(K )

a1

where K + 1 is the number of states in the BD, also referred to as levels. The term finite denotes the fact that the number of (n) levels is finite. The term nonhomogeneous refers to the fact that the transition rates {ai , i ∈ {0, 1, 2}, 0 ≤ n ≤ K } in Eq. (4) may depend on the level n. According to Eq. (4), in a finite BD process with K + 1 levels only transitions of the form n → n + 1 or n + 1 → n are allowed, for any state 0 ≤ n ≤ K − 1. When the BD process is in state 0 ≤ n ≤ K − 1, the time to the next transition of the (n) form n → n + 1 is exponentially distributed with parameter a0 . When the BD process is in state 1 ≤ n ≤ K , the time to the (n)

next transition of the form n → n − 1 is exponentially distributed with parameter a2 . The holding time in state 0 ≤ n ≤ K (n)

is exponentially distributed with average −1/a1 . Assuming that the BD process is irreducible and positive recurrent, the stationary probability vector π = (π (0) , π (1) , . . . , π (K ) ) of the BD process is the unique solution of π q = 0 and π 1 = 1. Quasi-birth-and-death (QBD) processes constitute a generalization of BD processes in which each level is expanded into a series of states called phases. In particular, a nonhomogeneous continuous-time QBD is a Markov process on the two-dimensional state space {(n, p) : 0 ≤ n ≤ η, 1 ≤ p ≤ ϕ(n)} which we partition as ∪0≤n≤η l(n), where l(n) = {(n, 1), (n, 2), . . . , (n, ϕ(n))} for 0 ≤ n ≤ η [21]. The first coordinate n is called the level, and the second coordinate p is called the phase of the state (n, p). The number η of levels might be finite or infinite, and the number of phases ϕ(n) may depend on the level n and is finite. The infinitesimal generator of the QBD is block-tridiagonal and has the form:



(0)

A1

 (1) A 2 Q =  

(0)

A0

(1)

A1

..

.

 .. ..

. .

(η)

  .  (η−1)  A

A2

0

(5)

(η)

A1

We will assume throughout the paper that all QBDs are irreducible and positive recurrent. Then the system π Q = 0, π 1 = 1, has a unique solution π = (π (0) , π (1) , . . . , π (η) ) referred to as the stationary probability vector of the QBD. A continuous-time Markov process {X (t ), t ∈ R+ } with state space E, transition rates q(x, y), x, y ∈ E is reversible if there is a positive measure π on E that satisfies the detailed balance equations [22, Definition 1.4]

π (x)q(x, y) = π (y)q(y, x),

(6)

with x, y ∈ E. In an ergodic process the average number of transitions per unit of time of the process from state i to state j is equal to the proportion of time π (i) that the process spends on state i times the rate at which it makes a transition towards state j [22, Chapter 2]. Thus, for an ergodic process, the equation above can be interpreted as the equality between the average number of transitions per unit of time of the process from state x to state y and in the reverse direction from y to x. We use hereinafter uppercase and lowercase letters to distinguish between QBD and BD processes. We also denote by {Xi(n) , i ∈ {0, 1, 2}, 0 ≤ n ≤ ηX } to the transition rate matrices of a QBD X . The transition rates of a BD x are denoted by

{x(i n) , i = {0, 1, 2}, 0 ≤ n ≤ K }. In addition, ⊗ and ⊕ stand for Kronecker product and sum, respectively [21]. 2. The problem

In this section we present the problem under study. Consider a system of K servers and no buffering capability. Customers arrive at the system from N different sources, with N < ∞. The customer arrival process from source i is an independent MMPP (D0,i , D1,i ), with 1 ≤ i ≤ N and a finite number of phases per level. When a customer from source i arrives at the system one of the following two things may happen: – If there is an available server in the system, the server is occupied by the customer during a random period of time called the service time, distributed according to a PH (τ i , Ti ) distribution. – If there is no available server in the system the customer is lost, since there is no buffering capability in the system. We want to compute the stationary probability that an arriving customer from any of the N sources finds all K servers in the system busy. That is, the stationary probability that an arriving customer is lost. We shall refer to this as the customer blocking probability β .

124

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

The customer blocking probability refers to the blocking of a customer from any of the sources. Thus, we must first identify the customer arrival process at the system, regardless of the source that customers come from. The counting function {N (t ), t ∈ R+ } of the number of customers arriving during (0, t ] at the system is given by N X

N (t ) =

Ni (t ),

(7)

i=1

where {Ni (t ), t ∈ R+ , 1 ≤ i ≤ N } is the counting function of each MMPP (D0,i , D1,i ). From this it follows that the customer arrival process is a MMPP (D0 , D1 ) with [23]: D0 =

N M

D0,i ;

D1 =

N M

i=1

where

LN

i =1

N M

D1,i ,

(8)

i=1

Ai is defined as:

Ai , A1 ⊕ · · · ⊕ AN .

(9)

i =1

Let us consider first the degenerate case in which K is equal to ∞, that is, we are dealing with a queuing system with an infinite number of servers and no buffering capability. We refer hereinafter to this system as the infinite server system. In this case the process {X (t ), t ∈ R+ } counting the number of busy servers in the system corresponds to the level of a QBD process X with η = ∞ levels and state space ∪n≥0 l(n), where l(n) = {(n, 1), (n, 2), . . . , (n, ϕ(n))}, and ϕ(n) denotes the number of phases at level n, ϕ(n) < ∞. The following result is rather straightforward, but we present it as a proposition for later reference. Proposition 1. The QBD process {X (t ), t ∈ R+ } is the superposition X (t ) =

N X

Xi (t ),

(10)

i =1

of N independent QBD processes {Xi , t ∈ R+ }1≤i≤N . Each QBD process Xi has a level {Xi (t ), t ∈ R+ } which counts the number of busy servers in an independent MMPP (D0,i , D1,i )/PH (τ i , Ti )/∞/∞ system. Proof. On one hand, it is possible to create a finite disjoint partition ζ1 ∪ · · · ∪ ζN of the set of infinite servers such that each subset ζi , i ∈ {1, . . . , N } is also infinite. The queueing system receives customer arrivals from the superposition of N independent MMPPs. Let us redirect customer arrivals from the i-th MMPP to the subset ζi , for 1 ≤ i ≤ N. From the independence of the N arrival processes and their associated service times we have that the state of each group of servers ζi evolves independently from one another. Thus, the number of busy servers in group ζi is equal to the level {Xi (t ), t ∈ R+ } of a QBD in an independent MMPP (D0,i , D1,i )/PH (τ i , Ti )/∞/∞ system, with 1 ≤ i ≤ N. Summing up the number of busy servers in each group leads to {X (t ), t ∈ R+ }, which concludes the proof.  The infinitesimal generator Qi of QBD Xi is given by:



(0,i)

X1  (1,i) Qi = X2

(0,i)

X0 (1,i) X1

..

.

 (1,i)

X0

..

..

.

.

 ,

(11)

with (n,i)

= D1,i ⊗ IT⊗i n ⊗ τ i ;

(n,i)

= D0,i ⊗ IT⊗i n + ID0,i ⊗

X0 X1

(n,i)

X2

= ID0,i ⊗

n X

n ≥ 0, n X

fTi (z , n);

n ≥ 0,

z =1

ft i (z , n);

(12)

n ≥ 1,

z =1

where IA denotes the identity matrix with the size of matrix A. In (12), t i = 0 − Ti 1, IT⊗i n is given by IT⊗i n = ITi ⊗ · · · ⊗ ITi ,

|

{z

n times

(13)

}

and fA (z , n) is obtained from IT⊗i n in (13) after replacing the z-th matrix ITi by matrix or vector A. By convention, IT⊗i 0 = 1 and

P0

z =1

(·) = 0.

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

125

The infinitesimal generator QX of QBD X has the form



(0)

X1  (1) QX = X2

(0)



X0 (1) X1

..

(1)

X0

..

.

..

.

.

 .

(14)

From Eq. (10) it follows that the generator QX can be obtained from the superposition of the infinitesimal generators

{Qi }1≤i≤N in (11). The superposition of independent QBD processes is studied in Section 3. From this section (see Table 1) (n) and (11), the rate matrices {Xi , i ∈ {0, 1, 2}, n ≥ 0} in the generator QX are fully characterized. Let us consider now the queueing system under study under the assumption that the number of servers K is finite. We refer to this system as the finite server system, as opposed to the infinite server system described before. The only difference between both systems is in the value of K ; the rest of the system parameters remain unchanged. We denote by e X the QBD that counts the number of busy servers in the finite server system. The following proposition relates the infinitesimal generator e Qe X of QBD X to the infinitesimal generator QX of QBD X from the infinite server system, given by Eq. (14).

e Proposition 2. The infinitesimal generator Qe X of QBD X , which counts the number of busy servers in the finite server system is given by: 

(0)

X1

 (1) X

 Qe X = 

(0)



X0

2



.. ..

..

.

.

(K −1)

.

X1 (K ) X2

  ,  (K −1)  X0 (K ) e X

(15)

1

e(K )

(K )

(K )

where the off-diagonal elements in X1 are equal to those in X1 and the diagonal elements are set so that (X2 The rest of the blocks in (15) are equal to their corresponding counterparts in (14).

(K ) +e X1 )1 = 0.

Proof. The finite and infinite server systems exhibit the same behavior when there are less than K customers in the system. This explains why the rate matrices in Qe X are equal to those in QX for all levels between 0 and K − 1. Let us assume now that both, the finite and infinite server systems have K customers and are in the same phase. In this case three possible events may take place. First, a service from one of the K customers in the system is completed. This event takes place with the same rate and produces the departure of the served customer in both systems. For this reason, (K ) the matrix X2 is the same in QX and Qe X. Second, a level transition takes place in the MMPP (i.e., a customer arrives at the system). This event takes place with the same rate in both systems, but produces a different effect in each one. In the infinite server system, the new incoming customer enters an available server and QX increases its level. In the finite server system, the new incoming customer is (K ) (K ) blocked (i.e., lost), and Qe must be different from e X1 . Notice that we are X remains in the same level. This means that X1 dealing with an MMPP, where by definition a customer arrival does not imply a change of phase in the MMPP. This implies that when the above mentioned event takes place, the phases of QX and Qe X are still the same, although QX is at one level higher than Qe X . Therefore, the different effect that this event produces on each system translates only in different values for (K )

(K )

the diagonal elements of X1 and e X1 . Third, a phase transition is produced, either in the service time of one of the K customers in the system or in the MMPP arrival process, without triggering a service time completion or a customer arrival. This event takes place with the same rate in both systems and produces the same effect. That is, a change of phase of the corresponding service time or MMPP in both (K ) (K ) systems. For this reason, the off-diagonal elements of e X1 are equal to their counterparts in X1 . (K )

Finally, the condition (X2

(K ) +e X1 )1 = 0 in the proposition comes from the requirement that Qe X 1 = 0. (0) e X

(K ) ) and π X e X

(0)



(K )

In what follows denote by π e ···π = (π X · · · π X ) to the stationary probability vectors of QBDs X = (π e X and X , respectively. It is well-known that if QBD X is reversible (see Section 1), then its truncation e X at level K is also e reversible. Moreover, the stationary distribution π e of X is proportional to π . The following proposition formally presents X X this result for further reference.

Proposition 3 (Truncation of a Markovian Reversible Process [22, Proposition 2.12]).: If X is a Markovian reversible process with respect to πX (x), x ∈ E, then its restriction e X to e E ⊆ E is also reversible with respect to πX (x), x ∈ e E. If in addition, X is ergodic and e X is irreducible, then e X is ergodic and its stationary distribution is

πeX (x) = πX (x)

X

πX (y), x ∈ e E.

y∈e E

We proceed now to describe the standard method for computing β in the following section.

(16)

126

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

2.1. The standard solution The customer blocking probability β can be expressed as the ratio of the average rate of customers blocked by the system rB over the average rate of customers arriving at the system rI [24]. That is, β = rB /rI . The average rate rI of customers arriving at the system is equal to rI = ε D1 1, where ε is the stationary probability vector of the phase process associated to the MMPP (D0 , D1 ) [21]. Since the MMPP (D0 , D1 ) is the superposition of N independent MMPPs, this can be also expressed as rI =

N X

rI ,i ,

(17)

i=1

where rI ,i = ε i D1,i 1 represents the average arrival rate of MMPP (D0,i , D1,i ) (see [21]). The average rate rB of customers blocked by the system can be expressed as rB = rI − rO , where rO represents the average rate of customers not blocked by the system, or the average customer departure rate. Recall that the level {e X (t ), t ∈ R+ } of QBD e X is equal to the number of busy servers at time t in the system. Each time that QBD e X increases its level represents the arrival of a customer that is not blocked by the system. Thus, rO can be expressed as the average rate at which QBD e X (n) increases its level. When QBD e X is at level n and phase p it increases its level with an average rate (X0 1)p , where vp denotes the p-th element in vector v . Thus, conditioning on the level and phase of QBD e X we conclude that the average arrival rate rO is given by: rO =

K −1 X n =0

π e(Xn) X0(n) 1.

(18)

Summarizing, the customer blocking probability β is given by

β =1−

rO

, (19) rI where rO and rI are given by Eqs. (18) and (17), respectively. In the M /M /K /K system β simplifies to the well-known Erlang B formula. As stated in the introduction, the main goal of this work is to find an alternative method to compute the customer blocking probability β in the system under study. In particular, we focus on the computation of rO since according to Section 5.1 it is the term in Eq. (19) presenting the highest computational cost. In the following two sections we introduce the two mathematical tools on which the new method relies: the superposition of QBDs and the simplified BD of a QBD process. 3. The superposition of 2 QBD processes Let us consider two independent QBD processes A and B, and let {XA (t ), t ∈ R+ } and {XB (t ), t ∈ R+ } denote, respectively, their level as a function of t. The superposition of QBD processes A and B is another QBD process C such that its level {XC (t ), t ∈ R+ } will satisfy: XC (t ) = XA (t ) + XB (t ),

(20)

for all t ∈ R+ . That is, QBD C increases/decreases its level whenever either QBD A or QBD B increases/decreases its level, and QBD C is at level 0 iff QBDs A and B are both at level 0. In this section we characterize QBD C as a function of the parameters of QBDs A and B. In order to characterize QBD C we define its state space and describe all the feasible transitions in QBD C with the corresponding transition rates. The characterization of a QBD C fulfilling Eq. (20) is not unique. We chose a characterization that allows one to express the transition rates of QBD C by means of Kronecker products and sums. We write the state of QBD C as a 4-tuple (nA , nB , pA , pB ), where nA and nB are the levels of QBDs A and B, respectively, and pA and pB are the phases of QBDs A and B, respectively. The state space SC of QBD C is given by

SC = {(nA , nB , pA , pB )},

(21)

with 0 ≤ nA ≤ ηA , 0 ≤ nB ≤ ηB , 1 ≤ pA ≤ ϕA (nA ) and 1 ≤ pB ≤ ϕB (nB ). According to (20), when QBD C is at level nC , QBDs A and B can be at any levels nA and nB such that nA + nB = nC , 0 ≤ nA ≤ ηA , and 0 ≤ nB ≤ ηB . Equivalently, when QBD C is at level nC , nA = nC − nB and m(nC ) ≤ nB ≤ M (nC ) , with m(nC ) = max(0, nC − ηA ), M (nC ) = min(ηB , nC ).

(22)

That is, for a given nC we have that nB will be constrained to be at least m(nC ) and at most M (nC ) . We partition the state space of QBD C in levels according to:

SC = ∪0≤nC ≤ηC lC (nC ),

(23)

where ηC is the maximum level of QBD C and it is equal to ηC = ηA + ηB (see Eq. (20)). The set lC (nC ) contains all possible combinations of states (nA , pA ) and (nB , pB ) from QBDs A and B, respectively, such that nA + nB = nC . We enumerate these

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

127

Table 1 Transition rates of the QBD superposition of QBDs A and B. Here, (A)(x,y) represents the element in the x-th row and y-th column of matrix A. Moreover, IA (nA )

and IB represent the identity matrix of size A1

(nB )

and B1

, respectively.

Transition

Rate

Condition

( nA )

?

(nA , nB , pA , pB ) → (nA + 1, nB , pA , pB ) (nA , nB , pA , pB ) → (nA , nB + 1, pA , p?B ) (nA , nB , pA , pB ) → (nA , nB , p?A , pB ) (nA , nB , pA , pB ) → (nA , nB , pA , p?B ) (nA , nB , pA , pB ) → (nA − 1, nB , p?A , pB ) (nA , nB , pA , pB ) → (nA , nB − 1, pA , p?B )

nA < ηA , 1 ≤ p?A ≤ ϕA (nA + 1)

(A0 ⊗ ) (n ) (IA ⊗ B0 B )(pB ,p?B ) ( nA ) (A1 ⊗ IB )(pA ,p?A ) (n ) (IA ⊗ B1 B )(pB ,p?B ) ( nA ) (A2 ⊗ IB )(pA ,p?A ) (n ) (IA ⊗ B2 B )(pB ,p?B ) IB (pA ,p? ) A

nB < ηB , 1 ≤ p?B ≤ ϕB (nB + 1) 1 ≤ p?A ≤ ϕA (nA ) 1 ≤ p?B ≤ ϕB (nB )

nA > 0, 1 ≤ p?A ≤ ϕA (nA − 1) nB > 0, 1 ≤ p?B ≤ ϕB (nB − 1)

combinations by the smallest nB first; we start with nB = m(nC ) and increase nB until it is equal to M (nC ) , where m(nC ) and M (nC ) are given in Eq. (22). For each level combination of QBDs A and B we enumerate their phases in ascending order, increasing first the phases in QBD B and then in QBD A. That is, lC (nC ) = {(nC − nB , nB , 1, 1), . . . , (nC − nB , nB , 1, ϕB (nB )), (nC − nB , nB , 2, 1), . . . ,

(nC − nB , nB , 2, ϕB (nB )), . . . , (nC − nB , nB , ϕA (nC − nB ), 1), . . . , (nC − nB , nB , ϕA (nC − nB ), ϕB (nB )) : m(nC ) ≤ nB ≤ M (nC ) }. According to this ordering the stationary distribution π

(24) (η )

C

= (π (C0) , . . . , π C C ) of QBD C has the form:

π (Cn) = (π (An−i) ⊗ π B(i) : m(n) ≤ i ≤ M (n) ),

(25)

(ηA )

(0)

(ηB )

(0)

for 0 ≤ n ≤ ηC , where π A = (π A , . . . , π A ) and π B = (π B , . . . , π B ) are the stationary distributions of QBDs A and B, respectively. In order to complete the characterization of QBD C, we provide in Table 1 its feasible transitions and their corresponding rates. We denote by means of an operator
the infinitesimal generator QC of QBD C as a function of the infinitesimal generators QA and QB of QBDs A and B, respectively: QC = QA
QB .

(26)

The operator
is fully characterized from Table 1 and the ordering of states given in Eqs. (23) and (24). (n) For further reference, we now characterize the level-up rate matrix C0 of QC for all levels n ≥ 0 such that n − m(n) < ηA (n) and M < ηB . Let us remark that this last condition ensures that when QBD C is at level n, QBDs A and B can increase their levels nA and nB in at least one unit. That is, for any nA , nB such that nA + nB = n we have that nA < ηA and nB < ηB . Let hn = M (n) − m(n) + 1 be the number of combinations of levels nA and nB from QBDs A and B, such that nA + nB = n, for 0 ≤ n ≤ ηC . We introduce the following matrices: (i)

(n−m(n) −i)

U1 = A0

⊗I

(m(n) +i)

(i)

U2 = I ⊗ B 0

(27)

,

(28) (n)

(n)

with i ∈ {0, 1, . . . , hn − 1} and 0 ≤ n ≤ ηC , and where A0 and B0 are the level-up rate matrices in QA and QB , respectively. Then,

 (n)

C0

  = 

(0)

U1

(0)

U2 (1) U1

 (1)

U2

..

.

..

.

(hn −1)

U1

(hn −1)

  , 

(29)

U2

for all n ≥ 0, such that n − m(n) < ηA and M (n) < ηB . In other cases in which this condition is not fulfilled it is possible to (n) (i) (i) express C0 in a similar way as a function of U1 and U2 , i ∈ {0, 1, . . . , hn − 1}. These cases and the characterization of

{C1(n) }0≤n≤ηC and {C2(n) }1≤n≤ηC are not relevant for our computations and thus fall outside the scope of this paper. 4. The simplified BD at level K of a QBD

This section introduces the novel concept of the simplified BD at level K of a QBD. Recall from Section 2 that the level of QBD X represents the number of busy servers in the degenerate version of our bufferless system (i.e., with ∞ servers). Intuitively, the simplified BD at level K of QBD X is a BD process that captures all the information from QBD X needed to

128

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

compute the customer blocking probability β in the above mentioned system when the number of servers is finite and equal to K . In what follows, we implicitly refer to the ‘‘simplified BD at level K of a QBD’’ by the shorter version ‘‘simplified BD of a QBD’’. We proceed now to provide a precise definition for the simplified BD of a QBD, and introduce a basic property related to the superposition of reversible QBD processes. The connection between simplified BDs and the computation of the blocking probability will become clear later in Section 5. Consider a nonhomogeneous QBD X with infinitesimal generator QX and ηX + 1 levels (possibly ηX = ∞). Let π e X =

(π e(X0) · · · π e(XK ) ) be the stationary probability vector associated to the generator QeX representing the truncation of QX at level K.

Definition 4. We define the simplified BD of a QBD X with generator QX as a BD process y with stationary probability vector π y = (πy(0) · · · πy(K ) ) given by:

πy(n) = π e(Xn) 1,

(30)

for 0 ≤ n ≤ K , and a generator qy given by Eq. (4) with level-up transition rates (n)

y0 =

π e(Xn) X0(n) 1 πy(n)

,

(31)

for 0 ≤ n ≤ K − 1. (n)

Observe that this definition is complete since {y2 }1≤n≤W may be computed from the global balance condition: (n+1)

y2

= y(0n) πy(n) /πy(n+1) ,

(32)

(n)

with 1 ≤ n ≤ K − 1, and {y1 }0≤n≤W must be set so that qy 1 = 0. The following proposition immediately follows from the definition of the simplified BD and it is given here for further reference. Proposition 5. The simplified BD of a QBD with η ≥ K levels is equal to the simplified BD of a truncation of the original QBD at level η? , with K ≤ η? ≤ η. In what follows we assume that η ≥ K for every QBD involved, so that all simplified BDs are well-defined according to Definition 4. The simplified BD has two properties. First, it spends the same proportion of time in each level as the original QBD process regardless of the phase (Eq. (30)). Second, it performs the same average number of transitions from any level n ≤ K − 1 to level n + 1 per unit of time than the original QBD does from any phase at level n to any phase at level n + 1 (Eq. (31)). From the positive recurrence of the QBD this immediately implies that the average number of transitions from any level n + 1 to level n per unit of time are also the same in the QBD and BD processes. The following theorem presents a useful property for simplified BDs related to the superposition of reversible QBDs. Theorem 6. The simplified BD of the superposition of two independent reversible QBDs is equal to the simplified BD of the superposition of the simplified BDs of each one of the two QBDs. (0)

(K )

Proof. We denote by e π X = (e πX , . . . , e π X ) to the stationary distribution of the truncation e X of QBD X at level K . We refer to Fig. 1 for a diagram of the proof. We have two independent reversible QBDs A and B with generators QA and QB , stationary probability vectors π A and π B , and ηA + 1 and ηB + 1 levels, respectively. We denote by BD a and BD b to their simplified BDs, with generators qa and qb and stationary probability vectors π a , π b , respectively. Let QBD C with stationary probability vector π C and ηC + 1 = ηA + ηB + 1 levels (see Eq. (20)) be the superposition of QBDs A and B. That is, QC = QA
QB . Moreover, let QBD D with stationary probability vector π D be the superposition of BDs a and b. That is, QD = qa
qb . We prove that the generators qc and qd of the simplified BDs c and d of QBDs C and D, respectively, are equal. From Definition 4 it follows that both, BD c and d have K + 1 levels (i.e. from level 0 to level K ). We proceed by proving that π c = π d and that

{c0(n) = d(0n) , 0 ≤ n < K }. This univocally defines the rest of the parameters in each generator qc and qd and implies their

equality (see Definition 4). (n) (n) Let us begin with the proof that π c = π d . On one hand we have πc = π e 1, for 0 ≤ n ≤ K , which follows from C Eq. (30). Since QBD C is the superposition of two independent reversible processes, it is also reversible [22, Proposition PK (n) (n) −1 (n) 2.14]. Therefore, we have that πc = ξC π C 1, with ξC = n=0 π C 1 (see Proposition 3). Expanding the expression for π

πc(n) = ξC−1

n X i =0

(n) C

π (An−i) 1 · π B(i) 1.

according to Eq. (25) leads to: (33)

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

129

Fig. 1. Diagram of the proof of Theorem 6.

(n)

On the other hand, we have for 0 ≤ n ≤ K that πd

= π e(Dn) 1, which follows from Eq. (30). QBD D is the superposition (n)

of two BD processes, and therefore it is also reversible. Using this fact and (see Proposition 3) gives πd

ξD =

(n)

= ξD−1 π (Dn) 1, with

n=0 π D 1. (n) Expanding the expression for π D according to Eq. (25) leads to:

PK

πd(n) = ξD−1

n X

πa(n−i) πb(i) .

(34)

i=0

(i)

= π e(Ai) 1. Since QBD A is reversible, we have πa(i) = ξA−1 π (Ai) 1, with PK (n) (n) (i) = ξB−1 π (Bi) 1, with ξB = ξA = n=0 π A 1. A similar argument permits us to write πb n=0 π B 1. Substituting these (i) (i) expressions for πa and πb in Eq. (34) leads to According to Eq. (30) we have in addition that πa

PK

πd(n) = ξD−1 ξA−1 ξB−1

n X

π (An−i) 1 · π B(i) 1.

(35)

i =0

In order to establish the equality between Eqs. (35) and (33) it remains to prove that:

ξC = ξD ξA ξB . (n)

Expanding π

ξC =

C

(36) in ξC =

K X n X

PK

(n)

n=0

π C 1 according to Eq. (25) leads to

π A(n−i) 1 · π (Bi) 1.

n=0 i=0

PK Pn (n−i) (i) π (Dn) 1 according to Eq. (25) leads to ξD = πb . As n=0 i =0 π a (i) (i) (i) −1 −1 (i) mentioned before, since QBDs A and B are reversible, we have πa =ξA π A 1 and πb =ξB π B 1. Substituting this in the expanded expression for ξD proves the equality in Eq. (36). (n) (n) We proceed now to prove that {c0 = d0 , 0 ≤ n < K }. (n) On one hand, for 0 ≤ n < K , c0 can be expressed as: Analogously, expanding π

(n)

c0

=

π (Cn) C0(n) 1 ξC · πc(n)

(n) D

in ξD =

PW

n=0

,

(37)

which follows from Eq. (31) and from the fact that QBD C is reversible. The simplified BDs a and b from QBDs A and B are well-defined and so ηA ≥ K and ηB ≥ K (see Definition 4). Moreover, for 0 ≤ n ≤ K − 1 we have m(n) = 0 and M (n) = n (see Eq. (22)). Thus, the conditions n − m(n) < ηA and M (n) < ηB for Eq. (n) (n) (n) (29) to be valid are fulfilled. We therefore expand C0 in C0 1 in Eq. (37) from Eq. (29). We also expand π C according to Eq. (25), which leads to: n P

(n)

c0

=

(π A(n−i) ⊗ π B(i) )[A(0n−i) ⊗ IB(i) + IA(n−i) ⊗ B(0i) ]1 0

i =0

(n)

ξC · πc

0

.

(38)

From the use of the distributive property of the Kronecker product over the sum and of its matrix-multiplication property

(A ⊗ B)(C ⊗ D) = (AC ⊗ BD), we have:

130

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140 n P

(n)

c0

=

[π (An−i) A(0n−i) ⊗ π (Bi) + π (An−i) ⊗ π B(i) B0(i) ]1

i =0

.

ξC · πc(n)

(39)

Identifying Eqs. (30) and (31) for the simplified BDs a and b and using the fact that QBDs A and B are reversible processes leads to:

(n)

c0

ξA ξB

n P

πa(n−i) a(0n−i) πb(i) + πa(n−i) πb(i) b(0i)

i =0

=

.

ξC · πc(n)

(40)

(n)

On the other hand, for 0 ≤ n < K , d0 can be expressed as: (n)

d0 =

π (Dn) D(0n) 1 ξD · πd(n)

,

(41)

which derives from Eq. (31) and the fact that QBD D is reversible. Once again, the conditions n − m(n) < ηA and M (n) < ηB for Eq. (29) to be valid are fulfilled. We therefore identify Eq. (n) (n) (29) in the expression D0 1 in Eq. (37). We also expand π D according to Eq. (25), which leads to: n P

(n)

d0

=

ξD · πd(n) n P

=

πa(n−i) πb(i) [a(0n−i) + b0(i) ]1

i=0

πa(n−i) a(0n−i) πb(i) + πa(n−i) πb(i) b(0i)

i=0

ξD · πc(n)

.

The equality between Eqs. (40) and (42) follows immediately from Eq. (36).

(42) 

We introduce now new notation and represent by ΨK (Q ) = q to the computation of the generator q of the simplified BD of a QBD from its generator Q . Recall from Eq. (26) that the notation QA
QB represents the superposition of QBDs A and B. With this notation, we can rewrite Theorem 6 as:

ΨK (QA
QB ) = ΨK [ΨK (QA )
ΨK (QB )],

(43)

where QBDs A and B are independent and reversible. This representation suggests that it is possible to extend Theorem 6 to any number N of independent QBDs. This extension is built according to the following corollary. Corollary 7. The simplified BD with generator q = ΨK (Q1
· · ·
QN ), of the superposition of a number N of independent reversible QBDs with generators {Qi , i ∈ {1, . . . , N }} fulfills: q = ΨK [ΨK (Q1
· · ·
QN −1 )
ΨK (QN )]. Proof. The proof follows immediately from Theorem 6 after identifying in Eq. (43) QA = Q1
· · ·
QN −1 and QB = QN .



We now propose an algorithm to compute the generator q of the BD obtained by simplifying the superposition of a number N of QBD. Let us remark that Eq. (44) can be iteratively applied to ΨK (Q1
· · ·
QN −1 ), ΨK (Q1
· · ·
QN −2 ) and successive terms until we have the superposition ΨK (ΨK (Q1 )
ΨK (Q2 )). This leads to Algorithm 1. Algorithm 1 Computation of q = ΨK (Q1
· · ·
QN ) Require: The infinitesimal generators {Qi }1≤i≤N of the N independent QBDs. Ensure: The infinitesimal generator q of the simplified BD of the superposition Q1
· · ·
QN . for i = 2 to N do Q1 ← ΨK (Q1 )
ΨK (Qi ) end for q ← Ψ K (Q1 ) return q

The functioning of the algorithm is illustrated in Fig. 2. Starting from the infinitesimal generators of the different QBDs the algorithm proceeds by iterating three simple actions: First, compute simplified BD of two QBDs. Second, superpose the resulting pair of BD processes. Third, compute the simplified BD of the superposition. When the for loop is over, the BD with generator ΨK (Q1 ) is equal to q, by direct application of Corollary 7.

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

131

Fig. 2. Diagram of the algorithm to compute the simplified BD of the superposition of the generators {Q1 , . . . , Q5 } of 5 independent QBDs.

Remark 8. Observe that according to Eq. (30) the last statement q ← ΨK (Q1 ) in Algorithm 1 also provides the stationary distribution π of the simplified BD of the QBD with generator Q1
· · ·
QN . As it is shown in Section 5.1, computing the simplified BD using Algorithm 1 is much faster than obtaining it from the superposition of the N QBDs. The reason for this is that in the second case the state space grows as more and more QBDs are superposed. In the first case its size remains bounded since after each superposition we substitute the resulting QBD by its simplified BD, which has always K + 1 states. This fact is exploited in the next section to derive an algorithm for the computation of rO in our system showing only linear complexity with N. 5. The proposed solution We now present our approach for computing rO , and therefore for computing the customer blocking probability β in the system. Recall from Section 2 that the level {X (t ), t ∈ R+ } of QBD X represents the number of busy servers in our bufferless system PN when K = ∞. From Eq. (10) we have that {X (t ), t ∈ R+ } is the superposition X (t ) = i=1 Xi (t ) of N independent processes {Xi (t ), t ∈ R+ }1≤i≤N . Each process {Xi (t ), t ∈ R+ } is the QBD of an independent MMPP (D0,i , D1,i )/PH (τ i , Ti )/∞/∞ system, for 1 ≤ i ≤ N. e Let Qe Xi be the infinitesimal generator of a QBD Xi describing the number of simultaneous packet transmissions on a Xi corresponds to the truncation of QBD Xi MMPP (D0,i , D1,i )/PH (τ i , Ti )/K /K system, for every 1 ≤ i ≤ N. That is, each QBD e to the state space ∪0≤n≤K lXi (n) (see Proposition 2). We are now ready to introduce the main theorem of this work, which presents the proposed solution for the computation of rO . Theorem 9. The average customer departure rate rO in the MMPP (D0 , D1 )/PH /K /K system can be computed from the infinitesimal generator qy and stationary distribution π y of the simplified BD of QBD X according to: rO =

K −1 X

πy(n) y(0n) .

(44)

n =0

Moreover, if the QBDs {Xi }1≤i≤N with infinitesimal generators {QXi }1≤i≤N are reversible, qy and π Algorithm 1 by using the initial sequence of infinitesimal generators:

{QeXi }1≤i≤N .

y

can be obtained from (45)

Proof. Eq. (44) immediately follows after identifying Eqs. (30) and (31) from the definition of the simplified BD of QBD X . We proceed to prove the second part of the theorem. From Eq. (10) it follows that QX = QX1
· · ·
QXN ,

(46)

where QX is the generator of the QBD associated to QBD X and QXi is the generator of QBD Xi , 1 ≤ i ≤ N. From this it follows that we can compute qy = ΨK (QX1
· · ·
QXN ) and π y (see Remark 8) from Algorithm 1. The generators {QXi }1≤i≤N have an infinite number of levels, and thus are not appropriate for numerical computations. However, Algorithm 1 provides the same result if we use it to compute ΨK (Qe X1
· · ·
Qe XN ). This holds since Proposition 5 implies that ΨK (QXi ) = ΨK (Qe  Xi ), for every 1 ≤ i ≤ N. Once we have rO from Theorem 9, we can obtain the customer blocking probability β from Eq. (19). In the following subsection we estimate the complexity of this solution and compare it with that of the standard solution in Section 2.1. 5.1. Complexity evaluation Let us now proceed to estimate the complexity for the computation of the customer blocking probability β in the bufferless system with K servers presented in Section 2. The complexity is measured in terms of the number of scalar products needed. For the sake of simplicity we make the following assumptions. First, the parameters {D0,i }1≤i≤N and {D1,i }1≤i≤N of the N MMPPS are all ζ × ζ matrices. Second, customer transmission times have the same PH distribution with ϑ phases, independently of the MMPP that the customer comes from.

132

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

According to Eq. (19), in order to compute the customer blocking probability β in the system under study, one needs to compute the average rate of the arrival (i.e. rI ) and departure (i.e. rO ) processes. In the computation of rI , common to both solutions, we exploit the fact that the MMPP (D0 , D1 ) is a superposition of N MMPPs and express rI according to Eq. (17). This equation requires the computation of ε(i), the stationary distribution of the phase process associated to MMPP (D0,i , D1,i ), which scales as O(ζ 3 ). Therefore, the complexity of the computation of rI is O(N · ζ 3 ), showing only a linear increase of complexity with N. The computation of rO is the most troublesome part, and it is where the new method provides an improvement. According to the standard method, rO is computed from Eq. (18). This requires finding the stationary distribution π e X . That is, the unique non-negative solution to the system π e X = 0, normalized with π e X Qe X 1 = 1. This system can be solved by means of the Linear Level Reduction algorithm [21,25], which is a version of the well-known Gaussian elimination method that acknowledges P 3 the block-tridiagonal structure of Qe . The complexity of this method is O ( W X n ), where Wn is the size of the matrix 0≤n≤K

    (K ) n+ϑ−1 n+ϑ−1 N e X1 in the generator Qe . In our case it can be seen that W = · ζ . The term comes from the fact that n X ϑ−1 ϑ−1 ( K ) e X1 must keep track of the number of PH distributions from the n available which are in each one of the ϑ possible phases, (K ) regardless of which distribution is in which phase. The term ζ N comes from the fact that e X1 must keep track of the number of MMPPs processes from the N available which are in each one of the ζ possible phases. In this case it is also important to know which MMPP is in which phase, since the MMPPs are different. Let us define S according to: S=

X  n + ϑ − 1 3 . ϑ −1 0≤n≤K

(47)

Then, we can express the complexity of the direct solution for the computation of rO as O(S · ζ 3N ).

(48)

As it can be clearly appreciated, this complexity scales exponentially with the number N of sources. This is very likely to lead to scalability problems for moderate to large values of N if ζ > 1. And usually, ζ > 1 because if ζ = 1 we are dealing with the trivial case in which each MMPP represents a Poisson process. We now estimate the complexity of our solution in Theorem 9. It can be observed from Fig. 2 that Algorithm 1 performs the following operations, listed with no specific order. Operation 1 consists of the computation of N −1 superpositions of pairs of BD processes. Operation 2 consists of the computation of N − 1 simplified BDs of QBDs that result from the superposition of a pair of BD processes. Operation 3 consists of the computation of N simplified BDs, one for each original QBD. We begin with Operation 1. The superposition of a pair of BD processes is computationally very efficient. This superposition results in a QBD process, and its computation comprises two steps. The first step leads to the characterization of the level-up transition matrices of the QBD from levels 0 until K − 1 (see Eq. (29). Since the QBD is the superposition of two BDs, all Kronecker products required in this step are with identity matrices (see Section 3). This means that the coefficients of the generators of the two BDs are just rearranged according to the Kronecker operations, but that no scalar products take place. The second step leads to the characterization of the stationary distribution of the QBD between levels 0 and K . The QBD has Wn = n + 1 phases at level n, 0 ≤ n ≤ K (see Eq. (24)). The value of the stationary distribution at each phase is the result of a scalar product (see Eq. (25)). Thus, P the number G of scalar products needed to compute the stationary distribution of the superposition of two BDs is equal to 0≤n≤K n + 1. That is, 3 K 2 + K + 1. (49) 2 2 Thus, the overall complexity of Operation 1 is O((N − 1) · G). We address now the complexity of Operation 2; the computation of N − 1 simplified BDs of the superposition of pairs of BD processes. The stationary distribution of the QBD is known from Operation 1 (see Eq. (25)). Therefore, the computation of the simplified BDs with Eqs. (30) and (31) is not expensive in terms of scalar products. More specifically, Eq. (30) does not require the computation of any scalar product. Regarding Eq. (31), it requires the computation of Wn 2 scalar products, where Wn = n + 1 is the number of phases at level n in every QBD, 0 ≤ n ≤ K − 1 (see Eq. (24)). Thus, the number P J of scalar products needed to compute the level-up transition rates in Eq. (31) for each QBD in Operation 2 is equal to 0≤n≤K −1 (n + 1)2 . That is, G=

1



1

1

1



· K. (50) 3 2 3 Thus, the overall complexity of Operation 2 is O((N − 1) · J ). Finally, we proceed with Operation 3, which requires the computation of a number of N simplified BDs. This requires the computation of the stationary distribution of each one of the N QBDs with generators {Qe Xi }1≤i≤N (see Theorem 9). For each J =

K2 +

K+

QBD e Xi , 1 ≤ i ≤ N we can use the linear level reduction algorithm [21,25], which scales with O(

P

the size of the

(K ) matrix e X1,i

0≤n≤K

Wn3 ), where Wn is

e in the generator Qe Xi . Since QBD Xi is fed with a single MMPP we have that Wn =

Thus, the complexity of Operation 2 is O(N · S · ζ ), where S is defined in Eq. (47). 3



n+ϑ−1 ϑ−1



· ζ.

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

133

Therefore, the overall complexity for computing rO with the proposed solution in Theorem 9 is: O((N − 1) · (G + J ) + N · S · ζ 3 ),

(51)

where G and J are given by Eqs. (49) and (50), respectively. The conclusion is that Theorem 9 effectively reduces the exponential increase of complexity with N in Eq. (48) to a linear increase. Concerning the memory complexity of Algorithm 1, it is independent of N. Indeed, the statement inside the for can be implemented with three matrices, A, B and C in the following way. Assume for instance that Q1 , obtained from the previous iteration, is stored in matrix C . We then compute ϕK (C ) and store the result in matrix A. We read matrix Qi (e.g., from a file or from the input parameter list) and store its result in matrix C . We compute ϕK (C ) and store the result in matrix B. We then compute ϕK (A
B) and store the result in matrix C . We are then ready to begin the next iteration. 6. A case study: Long-range dependent input traffic In this paper we present a method to compute the blocking probability in a multiserver bufferless system when the arrival process consists of the superposition of a finite number of independent MMPPs, each MMPP having a finite number of phases. In our analysis of this system we have highlighted the role played by the QBDs {Xi }1≤i≤N with infinitesimal generators {QXi }1≤i≤N , defined in Proposition 1. Theorem 9 states that our method is exact when these QBDs are reversible, and approximative otherwise. In the praxis, the reversibility assumption from Theorem 9 turns out to be rather strong and many MMPPs of interest lead to nonreversible QBDs {Xi }1≤i≤N . The main objectives of this section are two. First, we want to show how our method may be applied in order to approximate the blocking probability in a queueing system of practical interest. In particular, we chose a bufferless multiserver queueing system with LRD traffic as the system under study. In what follows we refer to this system as the LRD queueing system. Second, we want to study the accuracy of our method by means of numerical examples. The analysis of the LRD queueing system is far from being straightforward. A common approach to tackle this problem is to approximate the LRD process with the superposition of a finite number of MMPPs, each MMPP having a finite number of phases. In this case we refer to the resulting queueing system as the pseudo LRD queueing system. The term pseudo denotes here the fact that we are approximating LRD traffic by using a finite number of Markovian traffic sources, as opposed to the original LRD queueing system, where true LRD traffic is used. We take this approach and apply our method to the analysis of the pseudo LRD queueing system. This section is structured as follows. Section 6.1 presents the queueing system under study. It provides the details concerning how our method is used to compute the blocking probability in the pseudo LRD queueing system. Section 6.2 presents the main results of the numerical experiments. It compares the blocking probability values computed from our method with those measured from simulations of the pseudo LRD queueing system and of the original LRD queueing system. 6.1. The pseudo LRD queueing system We are interested in the study of the blocking probability in a bufferless multiserver queuing system with LRD input traffic, that is, in the LRD queueing system described above. The practical relevance of this system can be found for instance in the study of the blocking probability in future optical networking solutions such as OBS and OPS, since in both solutions LRD traffic is expected to be present at the optical domain [6–10]. In particular, the bufferless system with K servers can model the output port of an optical switch with no fiber delay lines in an OBS/OPS network with K wavelength channels per fiber. From now on, customers are interpreted as packets and service times as packet transmission times over a link of capacity C . The problem of computing the blocking probability in a bufferless OPS/OBS network has received much attention in the literature [13–20]. The standard assumption in these studies is that input traffic is Poisson [13–18,20]. The work in [19] constitutes an exception, since the authors consider an ON/OFF traffic model with exponentially distributed ON and OFF periods. To our knowledge, no previous analytical studies exist that compute the blocking probability in the presence of LRD traffic in the OPS/OBS network. The analysis of a bufferless multiserver queueing system with true LRD traffic is very challenging. A common approach is to emulate LRD traffic with the superposition of a finite number N of MMPP sources, each MMPP having a finite number of phases. This is the approach taken in [1–5]. In this case, we are working with the pseudo LRD queueing system introduced above. According to Section 5.1, the direct analysis of the pseudo LRD queuing system easily becomes prohibitive with moderate to high values of the number of traffic sources N, since its complexity exponentially scales with N. Therefore, we may consider using our method to compute the blocking probability. Indeed, Theorem 9 constitutes an efficient mechanism to speed-up this analysis, since its complexity linearly scales with N. However, the resulting QBDs {Xi }1≤i≤N are not reversible when the LRD traffic is emulated with any of the approaches presented in [1–5]. Thus, the solution provided by Theorem 9 must be regarded in this case as approximative.

134

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

In order to study the approximation error of Theorem 9 we use a continuous-time version of the traffic model presented in [3] to emulate the LRD traffic in the pseudo LRD queueing system. The details concerning this traffic model are presented below. 6.1.1. Emulating LRD traffic In [3] it is proven that the superposition of an infinite number of Markovian ON/OFF processes (which are also MMPPs) may indeed be LRD. The authors provide also the transition matrices for two examples of ON/OFF processes fulfilling this property. We take the transition matrices of each source from Example 2 in [3] in order to illustrate the use of Theorem 9. According to this example an ON/OFF source i, with i ∈ {1, . . . , ∞} is defined to be a two-state discrete-time Markov chain {Zi [n], n ∈ N} with irreducible and aperiodic transition matrix Pi :

 Pi =

1 − (1/(i + 1))p (1/(i + 1))q

 (1/(i + 1))p , 1 − (1/(i + 1))q

(52)

where q, p ∈ R. In the OFF state the process generates 0 packets/slot, and in the ON state it generates 1 packet/slot. The duration of the ON state is geometrically distributed with parameter 1 − βi = (1/(i + 1))q . Similarly the duration of the OFF state is geometrically distributed with parameter 1 − αi = (1/(i + 1))p . The superposition of an infinite number of such ON/OFF sources is LRD iff p < 2q + 1 and p ≥ q + 1, in which case, the Hurst parameter H is such that 1/2 < H ≤ 1 and it is given by H =

(3q − p + 1) 2q

.

(53)

The evaluation of the resulting bufferless multi-server queueing system is not immediate in the discrete-time domain due to the presence of simultaneous events. Indeed, in this case the Markov chain counting the number of busy servers is not a QBD, but a M /G/1-type, or skip-free to the left chain [21]. This makes the direct application of our results, valid for QBDs, more difficult. For this reason we switch hereinafter to the continuous-time domain. In what follows we define and use a continuous-time version {Zi (t ), t ∈ R+ } of each discrete ON/OFF process {Zi [n], n ∈ N} in which the duration of a time slot is substituted by an exponential distribution of parameter λ. Then, the duration of the OFF period is equal to the sum of F exponential random variables with parameter λ, where F is geometrically distributed with parameter 1 − αi . This random sum of random variables is equal to an exponential random variable of parameter λ(1 − αi ). Using this fact, the process {Zi (t ), t ∈ R+ } is a stationary PH-renewal process with inter-renewal distribution PH (τ i , Ti ), where:

 Ti =

 −λ λ(1 − βi ) . λ(1 − αi ) −λ(1 − αi ) τ i = (1, 0).

(54)

Consider now a finite number N of such ON/OFF processes and denote {Z (t ), t ∈ R+ } to their superposition: Z (t ) =

N X

Zi (t ).

(55)

i=1

This superposition constitutes the pseudo LRD process. LRD processes typically are characterized by three parameters: the Hurst parameter H, the average arrival rate m, and the coefficient of variation a of the number of arrivals falling in a constant interval of time [26]. The Hurst parameter of the LRD process approximated by {Z (t )} can be directly computed from Eq. (53). We proceed now to provide approximate expressions for the computation of the other two parameters (i.e., m and a) of a LRD process from its approximation {Z (t )}. e of {Z (t )}, which is The average arrival rate m of the LRD process is approximated by means of the average arrival rate m simply equal to the sum of the arrival rate of each source {Zi (t )}1≤i≤N :

e=λ m

N X i=1

(i + 1)q − 1 . (i + 1)q + (i + 1)p

(56)

In order to obtain from {Z (t )} a simple approximation e a of the parameter a of the LRD process we introduce a slight modification in the definition of a. This modification consists of interpreting a as the coefficient of variation of the number of arrivals falling in an interval of time which is exponentially distributed with parameter υ (in the original definition the interval was constant). This alternative definition should capture in essence the same second-order statistics as the original one, especially when the size of the constant interval is chosen as sufficiently small. We proceed now to obtain an approximation e a from {Z (t )} for this coefficient of variation a. We define Mi as the number of packet arrivals from the stationary ON/OFF source {Zi (t ), t ∈ R+ } in an interval of length (t , t + Υ ), where Υ has an exponential distribution with parameter υ , independent of the PH-renewal process. The number

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

135

of arrivals Mi has a geometric distribution with parameter

ζi = π i [υ(υ I − Ti )−1 ]1,

(57)

where π i = (τ i Ti 1) τ i Ti is the stationary distribution of the PH-renewal process (see Example 3.2.5 in [21]). Consider now the superposition {Z (t )} of N ON/OFF sources. From the independence of all ON/OFF sources, Mi is PN independent of Mj , for i 6= j, 1 ≤ i, j ≤ N. From this we have that the coefficient of variation e a of the number M = i=1 Mi + of arrivals from {Z (t ), t ∈ R } in an interval with exponential distribution of parameter υ is given by: −1

N P

e a=

−1

−1

(1 − ζi )ζi−2

i=1 N P

,

(58)

ζi

−1

i =1

that is, the quotient of the sum of the variances over the sum of the averages of {Mi }1≤i≤N . 6.1.2. Analysis of the pseudo LRD queueing system Let us assume now that the service times are exponentially distributed with parameter µ. Notice that the i-th PH-renewal process defined in Eq. (54) is also a MMPP, for 1 ≤ i ≤ N. Indeed, the parameters of the i-th MMPP are D0,i = Ti and D1,i = t i τ i , where t i = 0 − Ti 1. This process is a MMPP since from Eq. (54) we have that D1,i is diagonal. Thus, the customer arrival process at the system is a MMPP, superposition of N independent MMPPs (see Section 2). From this and Proposition 2 we have that the number of busy servers in a bufferless system with K servers (K finite), arrival process {Zi (t ), t ∈ R+ }, and the above mentioned service times can be described by the truncated QBD e Xi from Eq. (15). In particular, QBD e Xi has an infinitesimal generator Qe Xi with the following rate matrices: (j,i)

X1

(j,i)

X0

(j,i)

X2

 • λ(1 − βi ) = , λ(1 − αi ) • λβi 0 = , 0 0   jµ 0 = , 0 jµ 

0 ≤ j ≤ K, 0 ≤ j ≤ K − 1,

(59)

1 ≤ j ≤ K,

where the symbol • represents a scalar, such that Qe Xi 1 = 0. As mentioned before, each QBD e Xi with generator Qe Xi is not reversible. For instance, the average number of direct jumps (j,i)

from an OFF state at level j to an OFF state at level j + 1 is zero, since element (2, 2) in matrix X0 is zero. However, the average number of direct jumps from an OFF state at level j + 1 to an OFF state at level j can be different from zero, since (j,i) element (2, 2) in matrix X2 is different from zero. Thus, the detailed balance equations in Eq. (6) do not hold. 6.1.3. Fitting the pseudo LRD queueing system We fit the QBD processes {e Xi }1≤i≤N presented in Section 6.1.1 to a well-known real IP traffic trace showing asymptotic second-order self-similar scaling behavior (i.e., LRD). In particular, we work with the August 1989 Bellcore trace pAug of 106 interarrival times, as measured by Leland et al. [27]. Although slightly dated, this data set provides a well-known benchmark useful for examining the LRD of network traffic. For the fitting process we use the Hurst parameter estimation b H ≈ 0.8 from [28] for this trace, and we estimate m and a according to Eqs. (56) and (58). We describe now in more detail how to fit the parameters from the QBDs {e Xi }1≤i≤N to the data from this trace. The unknown parameters from all QBD processes {e Xi }1≤i≤N are λ, µ, p and q. The parameter µ can be directly fitted to the average packet transmission time E [st ]. This is proportional to the average packet size from the trace E [ps] according to E [st ] = E [ps]/C , where C is the capacity of the link through which traffic is sent. The rest of the parameters (λ, p and q) b and b are fitted with estimations b H, m a from the trace of the parameters H, m and a characterizing a typical LRD process. We now describe each one in turn. The estimation b H of the Hurst parameter can be obtained from a variety of estimators. We refer to [29] for a performance b of the average arrival rate is obtained from the inverse comparison of some of the most important ones. The estimation m of the sample mean of the packet interarrival times in the trace. The estimation b a of the coefficient of variation is taken from the sample mean and variance of the number of arrivals from the trace falling in adjacent intervals of which length constitutes the samples from an exponential distribution of parameter υ . As in the standard procedure with constant sized intervals, the model allows for arbitrary choice of the size of the intervals. In our case we have chosen 1/υ to be equal to the mean interarrival time in the traffic trace. b and b Once we have the estimations b H, m a we fit the parameters λ, p and q as follows. The parameter p depends on q b and b through b H (Eq. (53)), so only two degrees of freedom q and λ remain available for fitting m a. Solving for λ in Eq. (56) and substituting its value in the generator given by Eq. (59) permits us to compute q from Eq. (58). The value of λ is then obtained form Eq. (56).

136

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

Blocking probability β

10-1 10-2 10-3 10-4 10-5 Proposed Method Simulation

10-6 40

60 80 100 120 Utilization factor ρ (Erlangs)

140

160

Fig. 3. Packet blocking probability as a function of the utilization factor for a bufferless system with K = 160 channels. Packets arrive at the system according to the superposition of N independent PH-renewal sources characterized in Eq. (54). From top to down the curves refer to N = 1000, N = 100, N = 50, and N = 10 sources superposed. The dash-dotted curves present simulation results with 95%-confidence intervals for 100 sample paths of length 222 corresponding to the superposition of N PH-renewal sources. The marked curves present the blocking probability as computed from Theorem 9 corresponding to the superposition of N PH-renewal sources.

6.2. A numerical example In the analysis of the LRD queueing system we use two levels of approximation. First, we approximate the blocking probability in the LRD queueing system by considering the blocking probability in a pseudo LRD queueing system. This approximation is a consequence of not using true but emulated LRD traffic in our queueing system (see Section 6.1.1). Second, we approximate the blocking probability in the pseudo LRD queueing system by using our method in Theorem 9. This approximation is a consequence of the fact that our method is in this case not exact since the corresponding QBDs {Xi }1≤i≤N in the pseudo LRD queueing system are not reversible. In this section we test the accuracy of each level of approximation. We begin by evaluating the accuracy of the method proposed in Theorem 9 compared to the blocking probability values obtained from simulations of the pseudo LRD queueing system. The corresponding numerical results are presented in Section 6.2.1. We then focus on the general problem of computing the blocking probability in the original LRD queueing system. That is, we evaluate the accuracy of the method proposed in Theorem 9 compared to the blocking probability values measured from simulations of the original LRD queueing system. The corresponding numerical results are presented in Section 6.2.2. 6.2.1. Numerical evaluation of the pseudo LRD queueing system In this section we fit the parameters of the QBD processes {e Xi }1≤i≤N in Section 6.1.2 to the well-known traffic trace introduced in Section 6.1.3. We then use Theorem 9 to compute the customer blocking probability in the bufferless multi-server system. This computation is an approximation of the real value since the QBD processes {e Xi }1≤i≤N defined in Section 6.1.2 are not reversible. In order to evaluate the efficiency of this approximation we measure the customer blocking probability from the simulation of the corresponding pseudo LRD queueing system with K servers and exponential service times. In particular, we synthesize packet arrivals from each of of the N independent MMPPs defined in Eq. (54), which in this case are PH-renewal processes. The parameters of each MMPPs are obtained after fitting the QBD processes {e Xi }1≤i≤N defined in Eq. (59) to the input traffic trace, according to the method presented in Section 6.1.3. We then superpose the synthesized traces from the N MMPPs in order to generate packet arrivals for the simulation of the pseudo LRD queueing system. In what follows we refer to the utilization factor ρ as the ratio of the average arrival rate rI of the input process to the service rate µ of the exponential service times. The Erlang B formula is used as a representative value for the case in which the arrival process is Poisson. Fig. 3 shows simulation results for the cases of N = 10, N = 50, N = 100 and N = 1000 sources. These simulation results are compared to the analytical results obtained from our method. Fig. 3 shows that the blocking probability curves for N = 10, 50, 100, 1000 computed from Theorem 9 fall within the 95%-confidence interval of the corresponding simulation values. This implies that in our example Theorem 9 provides a good approximation of the pseudo LRD queuing system. We give now some insight into the reason why Theorem 9 provides in Fig. 3 such a good approximation to the blocking probability in the pseudo LRD queueing system, especially for large values of N. When the reversibility condition does not hold, Theorem 6 is approximative since we cannot apply [22, Proposition 2.31] in its proof. Proposition 2.31 in [22] states that the stationary distribution of a reversible Markov process in a set A is proportional to the stationary distribution of the truncation of this process to the set A. We refer hereinafter to this as the property of proportionality. The property of

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

a

b

100

100

10-10

10-10

10-20 0

10-20 60 10 Tr

un

ca ti

20 on

30 Le ve 40 lK

40 30 20

50 60

10

Sou

b

um

N rce

er i

Ut 55 iliz ati 50 on Fa 45 cto r ρ 40 (E rla 35 ng s)

137

40 30 20 30

10

ber

i

um

N rce

Sou

Fig. 4. (a) Euclidean distance between the vectors π e Xi,K and the truncated and normalized version of π e Xi,M at level K for 1 ≤ K ≤ 60, 1 ≤ i ≤ 40 and M = 10000 levels. (b) Euclidean distance between the vectors π e Xi,K and the truncated and normalized version of π e Xi,M at level K = 40 for 30 ≤ ρ ≤ 60,

e 1 ≤ i ≤ 40 and M = 10000 levels. In both figures the size of π e Xi,K is equal to 3 ∗ K , since from Eq. (59) each level of QBD Xi,K has 3 phases.

proportionality is the only property from reversible processes that is needed in order for Theorem 6, and ultimately also Theorem 9, to be exact. We show now empirically that although the QBD processes {Xi }1≤i≤N we are using are not reversible, they almost fulfill the property of proportionality. We argue that this is the reason why Theorem 9 provides such accurate results. Let us proceed now with the description of the experimental setup. The main idea is to quantify how close a Markov process is of fulfilling the property of proportionality. We do so by computing the Euclidean distance between the normalized stationary distribution of the Markov process in a set A and the stationary distribution of the truncation of this process in A. More specifically, consider the QBD process Xi with ∞ levels constructed from the rate matrices presented in Eq. (59), 1 ≤ i ≤ N. Let us introduce new notation and refer by QBD e Xi,K to the truncation at level K of QBD Xi from source i. The QBD processes {Xi }1≤i≤N have an infinite number of levels, and thus are not easy to manipulate numerically. From now on we use QBD e Xi,M with large M as an approximation of QBD Xi . Let π e Xi,M and π e Xi,K denote the stationary probability vectors of QBDs

e Xi,M and e Xi,K , respectively. We show in Fig. 4(a) the Euclidean distance between π e Xi,M , normalized in the state space subset {0, . . . , K }, and π eXi,K , with K  M and M = 10000.

We observe from Fig. 4(a) that the Euclidian distance rapidly decreases with increasing truncation level K . This suggests that the probability mass of the stationary distribution of QBD Xi,M is mostly located in its lower levels, at least for 1 ≤ i ≤ 40. We also observe from Fig. 4(a) that the number i of the source has an impact on the Euclidean distance as well. Indeed, for moderate to large values of K the Euclidean distance decreases with decreasing source number i. The reason for this can be found in the rate matrices in Eq. (59). As the source number i decreases, so does βi = 1−(1−/(i+1))q , making the term λβi in (j,i) (j,i) X0 decrease at every level j, 0 ≤ j ≤ K − 1. Since the rate matrix X2 does not depend on i, when i decreases the probability mass of the stationary distribution of the QBD of source i is shifted to the lower levels of the QBD. This makes successive (normalized) truncations at level K of the stationary distribution of QBD Xi,M resemble the stationary distribution of QBD Xi,K more. In other words, it makes the Euclidean distance between both stationary distributions decrease with decreasing source number i. We observe from Fig. 4(b) that the Euclidian distance rapidly decreases with decreasing utilization factor ρ . This is because with decreasing ρ the probability mass of the stationary distribution of QBD Xi,M is shifted towards its lower levels. This produces the same effect as in Fig. 4(a). That is, it makes the normalized truncation at level K of the stationary distribution of QBD Xi,M resemble the stationary distribution of QBD Xi,K more, reducing their Euclidean distance. The parameter i produces the same effect in Fig. 4(b) as in Fig. 4(a). That is, the Euclidean distance between both stationary distributions decrease with decreasing source number i. The conclusion from Fig. 4 is that as K increases and ρ decreases, the Euclidean distance decreases and thus Theorem 9 provides a better approximation to the packet blocking probability in our experimental setup. In the praxis, the parameter K is the number of servers or channels in the bufferless system under study. Therefore, Fig. 4(a) suggests that the larger the number of servers in the bufferless system the more accurate is the approximative analysis for LRD traffic presented in Section 6. This result is extremely convenient since in most cases of interest the number K of servers in the system under study is rather large. For instance, in OBS/OPS networks K represents the number of wavelengths in the DWDM fiber and it may be equal to 40, 80 or even 160 [6,7].

138

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

b

Blocking probability β

a 100

Blocking probability β

10-1

Erlang B

10-2 Erlang B MWM Beta/Exp (Simulation) N= 5 N= 10 N= 50 N= 100 N= 150

10-3

MWM Beta/Exp (Simulation) N= 10 N= 50 N= 100

10-4

10-1 0

10

20

30

40

50

Numbers of Servers K

60

70

80

10

20

30

40

50

60

70

80

Numbers of Servers K

Fig. 5. Customer blocking probability as a function of the number of servers in the system for an utilization factor of (a) ρ = 80 Erlangs and (b) ρ = 50 Erlangs. The MWM Beta/Exp (Simulation) curve shows the 95%-confidence intervals for 10 sample paths of 219 packet arrivals each.

6.2.2. Numerical evaluation of the LRD queueing system As in Section 6.2.1 we fit the parameters of the QBD processes {e Xi }1≤i≤N defined in Eq. (59) to the well-known traffic trace introduced in Section 6.1.3. We then use Theorem 9 to compute the customer blocking probability in the bufferless multi-server system. This time the computation is treated as an approximation of the blocking probability in the original LRD queueing system. Recall that this approximation is double. First, our method approximates the blocking probability in the pseudo LRD queueing system. Second, the blocking probability in the pseudo LRD queueing system is itself an approximation of that in the original LRD queueing system. In order to evaluate the efficiency of this approximation we measure the customer blocking probability from the simulation of the corresponding LRD queueing system with K servers and exponential service times. The packet arrival times are synthesized from an asymptotic second-order self-similar (i.e., a LRD) traffic model fitted to the same pAug trace. This LRD traffic model is the Beta Multifractal Wavelet Model (MWM Beta) presented in [28]. The reason for using the MWM Beta model and not the original pAug trace in the simulation is to be able to synthesize packet arrival times for several simulations and to compute confidence intervals for the customer blocking probability. The particular choice of the MWM Beta model relies on the following two facts. First, the B-MWM produces a sequence of positive values. This permits us to easily interpret these values as the interarrival times of the arrival process. Other LRD models such as fractional Gaussian noise (fGn) produce also negative values, and therefore are not so easy to interpret in terms of interarrival times. Second, the accuracy of the MWM Beta model fitted to the same pAug trace we are using has already been shown to be excellent in a comparative study presented in [28]. Fig. 5 shows a performance comparison in terms of customer blocking probability in a bufferless multi-server system as a function of the number K of servers for a constant utilization factor. The range of values for K is inspired from OBS/OPS networks [6]. According to our results, the difference between the MWM Beta and Erlang B curves is strongly affected by the magnitude of K . For instance, for K > 20 in Fig. 5(a) the Erlang B formula does not fall within the 95%-confidence interval of the blocking probability measured from MWM Beta traces (detail not appreciated in the Figure). In addition, notice in Fig. 5 that increasing the number of sources N in the analytical model leads to results closer to the simulation of the traces generated by the LRD model (i.e., the MWM Beta/Exp curve in the Figure). In Fig. 5(a) and (b) the analytical results from Theorem 9 fall within the 95%-confidence interval for 1 ≤ K ≤ 80 with N = 100 and N = 150 sources, respectively. The impact of the utilization factor ρ is reflected in Fig. 6. This figure plots the blocking probability as a function of the utilization factor for two different configurations with K = 40 and K = 160 servers. It can be observed once more that increasing the number of sources N in the analytical model leads to results closer to the simulation of the traces generated by the LRD model. The analytical results from Theorem 9 fall within the 95%-confidence interval for N = 100 and N = 1000 sources in Fig. 6(a) and Fig. 6(b), respectively. From Figs. 5 and 6 we observe that the number of sources N needed in the analytical model in order to have results within the 95%-confidence interval increases as K increases and ρ decreases. 7. Conclusions We have presented a new method to compute the customer blocking probability in bufferless multi-server systems for which the arrival process results from the superposition of a number N of independent MMPPs. This method is based on the novel concept of the simplified BD of a QBD. Its complexity scales linearly with N, in contrast with that of the standard solution, which does it exponentially. The method relies on an assumption concerning the reversibility of some related

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

a

139

b Blocking probability β

Blocking probability β

10-1

10-2

10-3

Erlang B MWM Beta/Exp (Simul.) N= 5

10-4

10-1

Erlang B MWM Beta/Exp (Simul.) N= 5

10-2

N= 10

N= 10

N= 50

N= 50

N= 100

N= 100

N= 500 N= 1000

25

30 35 40 45 50 Utilization factor ρ (Erlang)

55

60

80

100

120 140 160 Utilization factor ρ (Erlang)

180

200

Fig. 6. Customer blocking probability as a function of the utilization factor for a bufferless system with (a) K = 40 servers and (b) K = 160 servers. The MWM Beta/Exp (Simulation) curve shows the 95%-confidence intervals for (a) 10 sample paths of 219 packet arrivals each and (b) 100 sample paths of 222 packet arrivals each.

QBD processes. In Section 6 we presented empirical evidence suggesting that the new method constitutes an excellent approximation to the computation of the customer blocking probability in a system with LRD input traffic, which for instance finds applications in modern optical network design. Acknowledgements We thank anonymous Reviewer 1 for the current version of Algorithm 1, which requires considerably less memory resources than the previous version proposed by the authors. References [1] A.T. Andersen, B.F. Nielsen, A markovian approach for modeling packet traffic with long-range dependence, IEEE Journal on Selected Areas in Communications 16 (1998) 719–732. [2] T. Yoshihara, S. Kasahara, Y. Takahashi, Practical time-scale fitting of self-similar traffic with markov-modulated poisson process, Telecommunication Systems 17 (2001) 185–211. [3] F. Geerts, C. Blondia, Superposition of markov sources and long range dependence, in: IFIP Conference, vol. 121, 1998. [4] A. Nogueira, P. Salvador, R. Valadas, A. Pacheco, Markovian approach for modeling IP traffic behavior on several time scales in: SPIE Conference on Performance and Control of Next Generation Communication Networks, ITCom, vol. 5244, Orlando, FL, USA, 2003. [5] P. Salvador, R. Valadas, A. Pacheco, Multiscale fitting procedure using markov modulated poisson processes, Telecommunications Systems 23 (2003) 123–148. [6] Y. Chen, C. Qiao, X. Yu, Optical burst switching: A new area in optical networking research, IEEE Network 18 (2004) 16–23. [7] C. Gauger, H. Buchta, E. Patzak, J. Saniter, OBS nodes: Technology meets performance, in: Proceedings of TransiNet-Workshop, HHI, Berlin, 2002. [8] M. de Vega Rodrigo, S. Spadaro, M.-A. Remiche, C.D., J. Barrantes, J. Goetz, On the statistical nature of highly-aggregated Internet traffic, in: Proceedings of 4th International Workshop on Internet Performance, Simulation, Monitoring and Measurement, Salzburg, Austria, 2006. [9] A. Ge, F. Callegati, L. Tamil, On optical burst switching and self-similar traffic, IEEE Communications Letters 4 (2000) 98–100. [10] Z. Sui, Q. Zeng, S. Xiao, The investigation of self-similar traffic for OBS assembly, Journal of Optical Communications 26 (2005) 58–62. [11] V. Paxson, S. Floyd, Wide area traffic: The failure of Poisson modeling, IEEE/ACM Transactions on Networking 3 (3) (1995) 226–244. [12] J. Gao, I. Rubin, Multiplicative multifractal modelling of long-range-dependent network traffic, International Journal on Communication Systems 14 (2001) 783–801. [13] K. Dolzer, C.M. Gauger, J. Späth, S. Bodamer, Evaluation of reservation mechanisms for optical burst switching, AEÜ International Journal of Electronics and Communications 55 (1). [14] K. Dolzer, C.M. Gauger, On burst assembly in optical burst switching networks — A performance evaluation of just-enough-time, in: Proceedings of the 17th International Teletraffic Congress, ITC, Salvador, Brazil, 2001. [15] L. Xu, H. Perros, G. Rouskas, A queuing network model of an edge optical burst switching node, in: IEEE INFOCOM, San Francisco, USA, 2003. [16] N. Akar, E. Karasan, Exact calculation of blocking probabilities for bufferless optical burst switched links, in: BROADNETS, San José, USA, 2004. [17] Z. Rosberg, H. Vu, M. Zukerman, J. White, Performance analyses of optical burst-switching networks, IEEE Journal on Selected Areas in Communications 21 (7) (2003) 1187–1197. [18] H.L. Vu, M. Zukerman, Blocking probability for priority classes in optical burst switching networks, IEEE Communications Letters 6 (5) (2002) 214–216. [19] R. Vallejos, A. Zapata, M. Aravena, Fast blocking probability evaluation of end-to-end optical burst switching networks with non-uniform ON–OFF input traffic model, Photonic Network Communications 13 (2007) 217–226. [20] V. Tintor, P. Matavulj, J. Radunovic, Analysis of blocking probability in optical burst switched networks, Photonic Network Communications. [21] G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM, 1999. [22] R. Serfozo, Introduction to Stochastic Networks, Springer Verlag, 1999. [23] W. Fischer, K. Meier-Hellstern, The markov-modulated poisson process (MMPP) cookbook, Performance Evaluation 18 (1992) 149–171. [24] G. Lin, T. Suda, F. Ishizaki, Loss probability for a finite buffer multiplexer with the m/g /∞ input process, Telecommunication Systems 29 (2005) 181–197. [25] D. Gaver, P. Jacobs, G. Latouche, Finite birth-and-death models in randomly changing environments, Advances in Applied Probability 16 (1984) 715–731.

140

M. de Vega Rodrigo et al. / Performance Evaluation 67 (2010) 121–140

[26] I. Norros, On the use of fractional brownian motion in the theory of connectionless networks, IEEE Journal of Selected Areas in Communications 13 (6) (1995) 953–962. [27] W.E. Leland, M.S. Taqqu, W. Willinger, D.V. Wilson, On the self-similar nature of Ethernet traffic, D.P. Sidhu (Ed.), ACM SIGCOMM, San Francisco, CA, 1993. [28] R.H. Riedi, M.S. Crouse, V.J. Ribeiro, R.G. Baraniuk, A multifractal wavelet model with application to network traffic, IEEE Transactions on Information Theory 45 (1999) 992–1018. [29] T. Karagiannis, M. Faloutsos, R. Riedi, Long-range dependence: Now you see it, now you don’t!, in: Proceedings of GLOBECOM, Taipei, Taiwan, 2002.

Miguel de Vega Rodrigo received his Master degree in Telecommunications in 1999 from the Universidad Politécnica de Madrid, Spain. He was a senior engineer at the Information Communication department at Siemens AG., Munich. Currently he works as an IT consultant and is pursuing his Ph.D. on traffic engineering in Optical Burst Switching networks at the Université Libre de Bruxelles (ULB), Belgium. His research interests are mainly on network planning and performance.

Guy Latouche received the Ph.D. degree in Mathematics from the Université Libre de Bruxelles (ULB) in 1976. He is professor in the Department of Informatics at the ULB, and has recently been Dean of the Faculty of Sciences. His research interests include various aspects of computational probability: matrix methods in Markov models, traffic models for telecommunication systems, stochastic processes in the plane and nearly completely decomposable systems.

Marie-Ange Remiche received her Ph.D. degree in mathematics in 1999 from the Université Libre de Bruxelles (ULB), Belgium. Currently she is a Chargée de Cours in the Faculty of Computer Science of the University of Namur, Belgium. Her current research interests concern the modelling and the performance evaluation of systems using probability theory, and in particular, queuing theory or matrix analytic techniques.