Analysis of an Infinite-Server Queue with Markovian Arrival Streams

Analysis of an Infinite-Server Queue with Markovian Arrival Streams

Guidance

Professor Masao FUKUSHIMA Associate Professor Tetsuya TAKINE

Assistant Professor Nobuo YAMASHITA

Hiroyuki MASUYAMA

1999 Graduate Course in

Department of Applied Mathematics and Physics Graduate School of Informatics

Kyoto University

February 2001

Contents

1 Introduction 1

2 Model 1

3 Time-dependent joint distribution of the number of customers 3 3.1 Joint distribution of the number of customers . . . . 5 4 Numerically Feasible Formulas for Phase-Type Service Times 5 4.1 Time-dependent joint binomial moments . . . . 6 4.2 Time-dependent formula for phase-type services . . . . 8 4.3 Limiting formula for phase-type services . . . . 12

5 Numerical Examples 14

5.1 Impact of service time distribution on Var[N ] . . . . 15

5.2 Impact of arrival process . . . . 16

5.3 Impact of correlation in service time sequence . . . . 18

Abstract

This paper considers an infinite-server queue with multiple batch Markovian arrival streams.

The service time distribution of customers may be different for different arrival streams, and si-

multaneous batch arrivals from more than one stream are allowed. For this queue, we first derive

a system of ordinary differential equations which the time-dependent matrix joint generating func-

tion of the number of customers in the system satisfies. Next assuming phase-type service times,

we derive explicit and numerically feasible formulas for the time-dependent and limiting joint bi-

nomial moments. Further some numerical examples are provided to discuss the impact of system

parameters on the performance.

1 Introduction

Over the past few decades, several studies have been done on infinite-server queues with batch arrivals, e.g., [4, 5, 6, 7, 9]. In particular, Liu and Templeton [7] studied an infinite-server queue, where the arrival process is assumed to be a Markov renewal process, the batch size distribution and the service times of individual customers in the batch may depend on the state of the Markov renewal process, and service times of individual customers in a batch are independent and identically distributed (i.i.d.). To the best of our knowledge, this queue is the most general one among those studied in the past. They derived the time-dependent generating function and binomial moments of the number of customers in the system. However, most of those results are given in terms of integrals with respect to an Markov renewal function, and therefore numerical computations with those results are not easy to conduct, as stated in [7].

The main purpose of this paper is to develop numerically tractable formulas for infinite-server queues. For this purpose, this paper considers an infinite-server queue with multiple batch Markovian arrival streams. Roughly speaking, customer arrivals occur in the following way. There exists an irreducible finite-state Markov chain that governs the arrival process, and customers from respective arrival streams arrive with a predefined probability when a transition of the Markov chain happens.

Note that this arrival process includes a superposition of independent phase-type renewal batch arrivals and Markov modulated batch Poisson arrivals as special cases. We assume that the service time distribution of customers may be different for different arrival streams and simultaneous batch arrivals from more than one stream are allowed.

For this queue, we first derive a system of ordinary differential equations for the time-dependent matrix joint generating function of the number of customers in the system. This result is considered as a generalization of the known result for the P H/GI/∞ queue with single arrivals [8]. As shown in [8], applying a general-purpose numerical algorithm, we can compute the time-dependent joint distribution and the joint factorial moments of the number of customers in the system. Next we derive the time-dependent joint binomial moment matrices of customers in the system, and assuming phase-type service times, we obtain explicit and numerically feasible expressions of the time-dependent and the limiting joint binomial moment matrices. Furthermore, through numerical examples, we reveal the impact of system parameters on the time-dependent and the limiting performance.

The remainder of this paper is organized as follows. In Section 2, we describe the mathematical model. In Section 3, we derive a system of ordinary differential equations for the time-dependent matrix joint generating function of the number of customers in the system. In Section 4, we derive the time-dependent joint binomial moment matrices, and then assuming phase-type service times, we obtain numerically feasible expressions of the time-dependent and the limiting joint binomial moment matrices. Finally, in Section 5, we show some numerical examples and discuss the impact of system parameters on the performance.

2 Model

We consider an infinite-server queue fed by K arrival streams. Hereafter we call customers arriving from the νth (ν = 1, . . . , K ) stream as class ν customers. Let K denote a finite set of class indices, i.e., K = {1, . . . , K }.

Customer arrivals are governed by a time homogeneous, stationary Markov chain, which is called

the underlying Markov chain hereafter. The underlying Markov chain has a finite state space M =

{1, . . . , M } and it is assumed to be irreducible. The underlying Markov chain stays in state i ∈ M for

an exponential interval of time with mean µ ⁻¹ _i , and then changes its state to j ∈ M with probability

σ i,j . We define a nonnegative 1 × K vector n = (n ₁ , . . . , n K ) ∈ Z ⁺ , where Z = {n = (n ₁ , . . . , n _K ); n _ν = 0, 1, . . . for all ν ∈ K}, Z ⁺ = Z − {0}.

When a transition from state i to state j happens, n ν customers in class ν simultaneously arrive at the queue with probability σ i,j (n)/σ i,j , where σ i,i (0) is assumed to be zero for all i ∈ M and σ i,j (n) satisfies

σ i,j = ^X

∈Z

σ i,j (n).

This batch arrival process is considered as an extension of Markovian arrival streams (MASs) in [1].

We now introduce some notations which describe the above MAS. Let C denote an M × M matrix whose (i, j)th element C i,j (i, j ∈ M) is given by

C i,j =

( −µ i , i = j, σ i,j (0)µ i , i 6= j.

Further, for n ∈ Z ⁺ , we define D(n) as an M × M matrix whose (i, j)th element D _i,j (n) (i, j ∈ M) is given by

D i,j (n) = σ i,j (n)µ i , n ∈ Z ⁺ . Thus our marked batch MAS is characterized by (C, D(n)).

The infinitesimal generator of the underlying Markov chain is given by C + D, where D is defined as

D = ^X

∈Z ⁺

D(n).

We denote the (i, j )th element of D by D _i,j . We define π as the stationary probability vector of the underlying Markov chain. Note that π satisfies

π(C + D) = 0, πe = 1,

where e denotes a 1 × M column vector whose elements are all equal to one. We define λ ν (∈ K) as the arrival rate of class ν customers, i.e.,

λ _ν = π ^X

≥

ν

n _ν D(n)e, where e ν denotes the νth unit vector:

e _ν = (0,. . . ,0, 1, 0,. . . ,0).

νth

(2.1) Service times of class ν (∈ K) customers are assumed to be i.i.d. according to a distribution function H _ν (t) with finite mean b _ν . For simplicity, we assume that no customers are present in the system at time 0 throughout the paper.

Remark 2.1 Customer arrivals defined above can be viewed as follows. Let S t denote the state of the underlying Markov chain at time t. Then, given S t = i, the conditional joint probability that n ν

customers in class ν (ν ∈ K) simultaneously arrive at the queue during time interval (t, t + δt] and

S _t+δt = j is given by D i,j (n)δt + o(δt). Besides, given S t = i, event S _t+δt = j (j 6= i) with no arrivals

happens with probability C i,j δt + o(δt). The assumption σ i,i (n) = 0 for all i ∈ M implies that at

least one customer arrives whenever a transition from any state i to itself happens.

3 Time-dependent joint distribution of the number of customers

In this section, we consider the time-dependent joint generating function of the number of customers of each class in the system. Let T _m denote the mth (m = 1, 2, . . . , k) arrival epoch after time 0, where 0 < T ₁ < T ₂ < · · ·. Let X m,ν (ν ∈ K) denote the number of class ν customers arriving at time T m . We define N ν (t) as the number of class ν customers in the system at time t.

When no arrivals happen in time interval (0, t], we have N _ν (t) = 0 for all ν ∈ K. Thus E

"

Y

ν∈K

z _ν ^{N ν (t)} 1{S t = j, T ₁ > t} S ₀ = i

#

= ^h e ^{t i}

i,j , (3.1)

where 1{ξ} denotes an indicator function of event ξ.

Next we consider the case T _k ≤ t < T _k+1 for some k ≥ 1. Let ξ _k (k = 1, 2, . . .) denote the event T k ≤ t < T _k+1 . Given the event ξ k happens, we define x _m as (X _m,1 , . . . , X m,K ). Because customers are served independently, we have for |z ν | ≤ 1 (ν ∈ K)

E

"

Y

ν∈K

z ^N _ν ^{ν (t)}

T _m = t _m (m = 1, . . . , k), ξ _k , x _m = n _m (m = 1, . . . , k)

#

= Y k m=1

Y

ν m ∈K

h H ν m (t − t m ) + z ν m H ν m (t − t m ) ^{i n m,νm}

= Y k

m=1

Y

ν m ∈K

h 1 + (z _{ν m} − 1)H _{ν m} (t − t _m ) ^{i n m,νm} , (3.2)

where n m = (n _m,1 , . . . , n m,K ) ∈ Z ⁺ for all m = 1, 2, . . . , k and H ν (t) = 1 − H ν (t). Because events ξ k

(k = 1, 2, . . .) are exclusive, it follows from (3.2) that E

"

Y

ν∈K

z _ν ^{N ν (t)} 1{S t = j and T ₁ < t} S ₀ = i

#

= X ∞ k=1

X

1 ∈Z ⁺

X

2 ∈Z ⁺

· · · ^X

k ∈Z ⁺

Z t

0

dt ₁ Z t

t 1

dt ₂ · · · Z t

t _k−1

dt _k

· Y k m=1

Y

ν m ∈K

h 1 + (z ν m − 1)H ν m (t − t m ) ^{i n m,νm}

· ^h e ^{t 1} D(n ₁ )e ^{(t 2 −t 1 )} D(n ₂ ) · · · e ^{(t k −t k−1 )} D(n _k )e ^{(t−t k ) i}

i,j . (3.3)

We now define the time-dependent matrix joint generating function G ^∗ (t, z) of the number of customers of each class in the system at time t, where z = (z ₁ , . . . , z K ) with |z ν | ≤ 1 for all ν ∈ K.

Namely, G ^∗ (t, z) denotes an M × M matrix whose (i, j)th element represents E

"

Y

ν∈K

z _ν ^{N ν (t)} 1{S t = j} S ₀ = i

#

, i, j ∈ M.

It then follows from (3.1) and (3.3) that G ^∗ (t, z) = e ^t +

X ∞ k=1

X

1 ∈Z ⁺

X

2 ∈Z ⁺

· · · ^X

k ∈Z ⁺

Z t 0 dt ₁

Z t

t 1

dt ₂ · · · Z t

t k−1

dt k

· Y k

m=1

Y

ν m ∈K

h 1 + (z ν m − 1)H ν m (t − t m ) ^{i n m,νm}

· e ^{t 1} D(n ₁ )e ^{(t 2 −t 1 )} D(n ₂ ) · · · · · e ^{(t k −t k−1 )} D(n k )e ^{(t−t k )} . (3.4) Further letting u j = t − t _k+1−j (j = 1, 2, . . . , k) and rearranging terms, we have

G ^∗ (t, z) = e ^t + X ∞ k=1

X

k ∈Z ⁺

X

k−1 ∈Z ⁺

· · · ^X

1 ∈Z ⁺

Z t 0 du k

Z u _k

0 du _k−1 · · · Z u 2

0 du 1

· Y k

m=1

Y

ν m ∈K

h 1 + (z ν m − 1)H ν m (u m ) ^{i n m,νm}

· e ^{(t−u k )} D(n k )e ^{(u k −u k−1 )} D(n _k−1 ) · · · · · e ^{(u 2 −u 1 )} D(n ₁ )e ^{u 1}

= e ^t + X ∞ k=1

Z t 0 du k

Z u _k

0 du _k−1 · · · Z u 2

0 du ₁ e ^{(t−u k )} D ^∗ (u k , z)

· e ^{(u k −u k−1 )} D ^∗ (u _k−1 , z) · · · · · e ^{(u 2 −u 1 )} D ^∗ (u 1 , z)e ^{u 1} , (3.5) where

D ^∗ (t, z) = ^X

∈Z ⁺

Y

ν∈K

h 1 + (z ν − 1)H ν (t) ^{i n ν} D(n). (3.6)

Theorem 3.1 For any closed interval of t where all of H _ν (t) is continuous, G ^∗ (t, z) in (3.5) satisfies the following differential equation:

∂

∂t G ^∗ (t, z) = [C + D ^∗ (t, z)] G ^∗ (t, z), (3.7) Further (3.7) has a unique continuous solution in [0, ∞) with G ^∗ (0, z) = I.

Proof. Pre-multiplying both sides of (3.5) by exp(−Ct), we have e ^{− t} G ^∗ (t, z) = I +

X ∞ k=1

Z t

0

du k

Z u k

0

du _k−1 · · · Z u 2

0

du ₁ e ^{− u k} D ^∗ (u k , z)

· e ^{(u k −u k−1 )} D ^∗ (u _k−1 , z) · · · · · e ^{(u 2 −u 1 )} D ^∗ (u ₁ , z)e ^{u 1} .

(3.8) Differentiating both sides of (3.8) with respect to t yields

e ^{− t} ∂

∂t G ^∗ (t, z) − e ^{− t} CG ^∗ (t, z)

= e ^{− t} D ^∗ (t, z)e ^t + e ^{− t} D ^∗ (t, z)

· X ∞ k=2

Z t

0

du _k−1 Z u k−1

0

du _k−2 · · · Z u 2

0

du ₁ e ^{(t−u k−1 )} D ^∗ (u _k−1 , z)

· e ^{(u k−1 −u k−2 )} D ^∗ (u _k−2 , z) · · · · · e ^{(u 2 −u 1 )} D ^∗ (u ₁ , z)e ^{u 1}

= e ^{− t} D ^∗ (t, z)G ^∗ (t, z). (3.9)

Rearranging terms in (3.9) and pre-multiplying both sides by exp(C t), we obtain (3.7). The uniqueness of the solution of (3.7) with G ^∗ (0, z) = I follows from the well-known results of ordinary differential

equations (see [2], p.167, Theorem 1). 2

Remark 3.1 Ramaswami and Neuts [8] obtained a system of ordinary differential equations which the generating function of the number of customers in the P H/GI/∞ queue satisfies. Theorem 3.1 is considered as a generalization of the result in [8].

3.1 Joint distribution of the number of customers

We now consider the time-dependent joint distribution of the number of customers of each class in the system. Let D(t, n) (n ∈ Z ) denotes an M × M matrix which satisfies

D ^∗ (t, z) = X ∞ n 1 =0

· · · X ∞ n _K =0

z ^{n 1} · · · z ^{n K} D(t, n),

where D ^∗ (t, z) is given by (3.6). Also, let L(t, n) (n ∈ Z ⁺ ) denote an M × M matrix whose (i, j)th element represents

Pr[N ₁ (t) = n ₁ , . . . , N K (t) = n K , S t = j | S ₀ = i].

For a given n = (n ₁ , . . . , n K ) ∈ Z , comparing the coefficient matrices of z ^{n 1} · · · z ^{n K} on both sides of (3.7), we obtain

d

dt L(t, n) = (C + D(t, 0))L(t, n) + ^X

0 ≤ ≤

6= 0

D(t, m)L(t, n − m).

Therefore the time-dependent joint distributions L(t, n) (0 ≤ n ≤ m) are given by the solution of a system of ordinary differential equations which can be solved numerically by a general-purpose numerical algorithm (e.g., see [8]).

Also, from (3.7), we can derive a system of ordinary differential equations whose solution provides the joint factorial moments for the number of customers of each class in the system at time t. We define Ω(t, m) and Ψ(t, m) as

Ω(t, m) = lim

z 1 →1 · · · lim

z _K →1

∂ ^{m 1}

∂z ^m ₁ ¹ · · · ∂ ^{m K}

∂z _K ^{m K}

∂

∂t G ^∗ (t, z), Ψ(t, m) = lim

z 1 →1 · · · lim

z K →1

∂ ^{m 1}

∂z ^m ₁ ¹ · · · ∂ ^{m K}

∂z _K ^{m K} D ^∗ (t, z).

By differentiating both sides of (3.7) m ν (ν ∈ K) times with respective to z ν and setting z ν = 1 for all ν ∈ K, we obtain

∂

∂t Ω(t, m) = (C + D)Ω(t, m) + ^X

0 ≤ ≤ 6= 0

Y

ν∈K

m _ν n ν

!

Ψ(t, n)Ω(t, m − n).

Thus the time-dependent joint factorial moment matrices Ω(t, n) (0 ≤ n ≤ m) are given by the solution of a system of ordinary differential equations, too.

4 Numerically Feasible Formulas for Phase-Type Service Times

In this section, we consider the joint binomial moments of the number of customers of each class in

the system and we develop numerically feasible formulas for the joint binomial moments, assuming

phase-type service times. We define B(t, m) (m ∈ Z ⁺ ) as an M × M matrix whose (i, j)th element represents

[B(t, m)] _i,j = E

"

Y

ν∈K

N _ν (t) m ν

!

1{S _t = j}

S ₀ = i

#

. (4.1)

B(t, m) is called the mth joint binomial moment matrix of the number of customers of each class in the system at time t.

4.1 Time-dependent joint binomial moments

We define G ^∗ _B (t, ω) as the matrix joint binomial moment generating function of the number of cus- tomers of each class in the system at time t, i.e.,

G ^∗ _B (t, ω) = e ^{( + )t} + ^X

∈Z ⁺

Y

ν∈K

ω _ν ^{m ν} B(t, m), (4.2)

where ω = (ω ₁ , ω ₂ , . . . , ω K ) and |ω ν + 1| ≤ 1 for all ν ∈ K. Note here that G ^∗ _B (t, ω) is given in terms of G ^∗ (t, z):

G ^∗ _B (t, ω) = G ^∗ (t, ω ₁ + 1, . . . , ω K + 1), (4.3) Theorem 4.1 The time-dependent matrix joint binomial moment generating function G ^∗ _B (t, ω) is given by

G ^∗ _B (t, ω) = e ^{( + )t} +

X ∞ k=1

Z t 0 du k

Z u k

0 du _k−1 · · · Z u 2

0 du ₁ e ^{( + )(t−u k )} D ^∗ _B (u k , ω)

· e ^{( + )(u k −u k−1 )} D ^∗ _B (u _k−1 , ω)

· · · · · e ^{( + )(u 2 −u 1 )} D ^∗ _B (u 1 , ω)e ^{( + )u 1} , (4.4) where

D ^∗ _B (t, ω) = D ^∗ (t, ω 1 + 1, . . . , ω K + 1) − D. (4.5) Proof. It is easy to see from (3.7) that G ^∗ _B (t, ω) in (4.3) satisfies the following differential equation.

∂

∂t G ^∗ _B (t, ω) = [C + D + D ^∗ _B (t, ω)] G ^∗ _B (t, ω), (4.6)

G ^∗ _B (t, ω) = I , (4.7)

where D ^∗ _B (t, ω) is given by (4.5). In what follows, we shall show (4.4) is a solution of (4.6).

Pre-multiplying both sides of (4.4) by exp[−(C + D)t], we have e ^{−( + )t} G ^∗ _B (t, ω) = I +

X ∞ k=1

Z t 0 du k

Z u _k

0 du _k−1 · · · Z u 2

0 du 1 e ^{−( + )u k} D ^∗ _B (u k , ω)

· e ^{( + )(u k −u k−1 )} D ^∗ _B (u _k−1 , ω)

· · · · · e ^{( + )(u 2 −u 1 )} D ^∗ _B (u ₁ , ω)e ^{( + )u 1} . (4.8)

Differentiating both sides of (4.8) with respect to t yields e ^{−( + )t} ∂

∂t G ^∗ _B (t, ω) − e ^{−( + )t} (C + D)G ^∗ _B (t, ω)

= e ^{−( + )t} D ^∗ _B (t, ω)

·

"

e ^{( + )t} + X ∞ k=2

Z t

0

du _k−1 Z u k−1

0

du _k−2

· · · · · Z u 2

0

du ₁ e ^{( + )(t−u k−1 )} D ^∗ _B (u _k−1 , ω)e ^{( + )(u k−1 −u k−2 )}

· D ^∗ _B (u _k−1 , ω) · · · · · e ^{( + )(u 2 −u 1 )} D ^∗ _B (u ₁ , ω)e ^{( + )u 1 i}

= e ^{−( + )t} D ^∗ _B (t, ω)G ^∗ _B (t, ω). (4.9)

Thus, pre-multiplying both sides of (4.9) by exp[(C + D)t], we see that G ^∗ _B (t, ω) in (4.4) satisfies the

differential equation (4.6) and (4.7). 2

We now rewrite G ^∗ _B (t, ω) in (4.4) to be a more appealing form. Note first that D ^∗ _B (t, ω) = ^X

∈Z ⁺

Y

ν∈K

h 1 + ω ν H ν (t) ^{i n ν} D(n) − D

= ^X

∈Z ⁺

Y

ν∈K





n ν

X

m ν =0

n ν

m ν

!

ω ^m _ν ^ν H ^m _ν ^ν (t)



 D(n) − D

= ^X

∈Z

"

Y

ν∈K

ω ^m _ν ^ν H ^m _ν ^ν (t)

# X

≥ 6= 0

Y

ν∈K

n ν

m _ν

!

D(n) − D

= ^X

∈Z ⁺

D(n) + ^X

∈Z ⁺

"

Y

ν∈K

ω ^m _ν ^ν H ^m _ν ^ν (t)

# X

≥

Y

ν∈K

n ν

m ν

!

D(n) − D

= ^X

∈Z ⁺

"

Y

ν∈K

ω ^m _ν ^ν H ^m _ν ^ν (t)

# X

≥

Y

ν∈K

n ν

m _ν

! D(n)

= ^X

∈Z ⁺

ω ^{( )} H ^{( )} (t) ^c D(m), (4.10) where

ω ^{( )} = ^Y

ν∈K

ω ^m _ν ^ν , m ∈ Z ⁺ , H ^{( )} (t) = ^Y

ν∈K

H ^m _ν ^ν (t), m ∈ Z ⁺ , (4.11)

D(m) = c ^X

≥

Y

ν∈K

n ν

m ν

!

D(n), m ∈ Z ⁺ . (4.12) Thus, with (4.10), G ^∗ _B (t, ω) in (4.4) is rewritten to be

G ^∗ _B (t, ω) = e ^{( + )t} + X ∞ k=1

X

1 ∈Z ⁺

X

2 ∈Z ⁺

· · · ^X

k ∈Z ⁺

ω ^{( 1 +···+ k )}

· Z t

0 du k

Z u _k

0 du _k−1 · · · Z u 2

0 du 1 H ^{( k )} (u k )H ^{( k−1 )} (u _k−1 )

· · · · · H ^{( 1 )} (u ₁ )e ^{( + )(t−u k ) c} D(l _k )e ^{( + )(u k −u k−1 ) c} D(l _k−1 )

· · · · · e ^{( + )(u 2 −u 1 ) c} D(l ₁ )e ^{( + )u 1}

= e ^{( + )t} + ^X

∈Z ⁺

ω ^{( )}

| |

X

k=1

X

~ _k ∈L k ( )

· Z t

0

du k

Z u k

0

du _k−1 · · · Z u 2

0

du ₁ H ^{( k )} (u k )H ^{( k−1 )} (u _k−1 )

· · · · · H ^{( 1 )} (u ₁ )e ^{( + )(t−u k )} D(l ^c k )e ^{( + )(u k −u k−1 )} D(l ^c _k−1 )

· · · · · e ^{( + )(u 2 −u 1 )} D(l ^c 1 )e ^{( + )u 1} , (4.13) where

|m| = ^X

ν∈K

m ν , m ∈ Z, l _j = (l _j,1 , . . . , l _j,K ) ∈ Z ⁺ ,

~ l _k = ⁿ (l ₁ , . . . , l _k ) | l _j ∈ Z, j = 1, 2, . . . , k ^o , L k (m) = ⁿ ~ l _k = (l ₁ , . . . , l _k ) l ₁ + · · · + l _k = m,

c D(l _j ) 6= O, l _j ∈ Z ⁺ , j = 1, 2, . . . , k ^o . Therefore comparing (4.2) and (4.13), we obtain the following theorem.

Theorem 4.2 The mth joint binomial moment matrix B(t, m) is given by B(t, m) =

| |

X

k=1

X

~ k ∈L k ( )

Z t 0 du k

Z u k

0 du _k−1 · · · Z u 2

0 du ₁

· H ^{( k )} (u k )H ^{( k−1 )} (u _k−1 ) · · · · · H ^{( 1 )} (u ₁ )

· e ^{( + )(t−u k ) c} D(l k )e ^{( + )(u k −u k−1 ) c} D(l _k−1 )

· · · · · e ^{( + )(u 2 −u 1 ) c} D(l ₁ )e ^{( + )u 1} . (4.14) Remark 4.1 Note that H ^{( )} (t) has probability meanings. To see this, suppose m ν (ν ∈ K) customers of class ν simultaneously arrive. Let H ν,i (ν ∈ K, i = 1, . . . , m ν ) denote a random variable representing a service time of the ith customer of class ν. Then H ^{( )} (t) represents the complementary distribution function of min ν∈K,i=1,...,m ν H ν,i , because

H ^{( )} (t) = ^Y

ν∈K m ν

Y

i=1

Pr(H ν,i > t)

= Pr

ν∈K,i=1,...,m min ν

H ν,i > t

.

Remark 4.2 We see from (4.14) that B(t, m) is independent of all D(l) (l ^c ≥ m, l 6= m).

4.2 Time-dependent formula for phase-type services

We now assume that the service time distribution of class ν customers is a phase-type distribution with representation (β _ν , T _ν ):

H ν (x) = β _ν e ^{ν x} e, ν ∈ K. (4.15)

Then using the properties of Kronecker product ⊗ and Kronecker sum ⊕:

(U 1 U ₂ · · · · · U _n ) ⊗ (V 1 V ₂ · · · · · V _n )

= (U ₁ ⊗ V ₁ )(U ₂ ⊗ V ₂ ) · · · · · (U _n ⊗ V _n ), ∀n ≥ 1, exp(U ) ⊗ exp(V ) = exp(U ⊕ V ),

we rewrite H ^{( )} (x) to be

H ^{( )} (x) = ^Y

ν∈K

β _ν e ^{ν x} e ^{l ν}

= β ^{< >} exp(T ^{[ ]} x)e, where β ^{< >} and T ^{[ ]} (l ∈ Z ⁺ ) are given by

β ^{< >} = β ₁ ⊗ · · · ⊗ β ₁

| {z }

l 1

⊗ β ₂ ⊗ · · · ⊗ β ₂

| {z }

l 2

⊗ · · · ⊗ β _K ⊗ · · · ⊗ β _K

| {z }

l K

, T ^{[ ]} = T ₁ ⊕ · · · ⊕ T ₁

| {z }

l 1

⊕ T ₂ ⊕ · · · ⊕ T ₂

| {z }

l 2

⊕ · · · ⊕ T _K ⊕ · · · ⊕ T _K

| {z }

l K

.

Thus B(t, m) in (4.14) is rewritten to be B(t, m) =

| |

X

k=1

X

~ _k ∈ L k ( )

F k (t,~ l k ), (4.16)

where

F k (t,~ l k ) = Z t

0

du k

Z u k

0

du _k−1 · · · Z u 2

0

du ₁ Y k

j=1

h β ^{< j >} exp(T ^{[ j ]} u j )e ⁱ

· e ^{( + )(t−u k ) c} D(l _k )e ^{( + )(u k −u k−1 ) c} D(l _k−1 )

· · · · · e ^{( + )(u 2 −u 1 )} D(l ^c ₁ )e ^{( + )u 1} . (4.17) To obtain a numerically feasible formula for F _k (t,~ l _k ), we shall rewrite (4.17) by considering the time-reversed arrival process. Note that the time-reversed arrival process is a batch marked MAS with representation (C ⁻ , D ⁻ (n)), where

C ⁻ = diag (π) ⁻¹ C ^T diag (π) , D ⁻ (n) = diag (π) ⁻¹ D(n) ^T diag (π) ,

with an M × M diagonal matrix diag (π) whose (i, i)th element is equal to the ith element of π.

Therefore (4.17) is rewritten to be F k (t,~ l k ) =

Y k η=1

b(l η )d(l η ) · diag (π) ⁻¹

·

" Z t 0 du k

Z u _k

0 du _k−1 · · · Z u 2

0 du ₁ Y k

j=1

κ ^{( j )} exp(T ^{[ j ]} u j )(−T ^{[ j ]} )e

· exp Q ⁻ u ₁ ^c D ⁻ (l ₁ ) exp Q ⁻ (u ₂ − u ₁ ) ^c D ⁻ (l ₂ )

· · · · · exp Q ⁻ (u k − u _k−1 ) ^c D ⁻ (l K ) · exp Q ⁻ (t − u k )

# T

· diag (π) , (4.18)

where

b(l) = Z _∞

0

dx β ^{< >} exp(T ^{[ ]} x)e = β ^{< >} (−T ^{[ ]} ) ⁻¹ e, (4.19) d(l) = max

i∈M

h diag (π) ^{−1 c} D(l) ^T diag (π) e ⁱ

i , (4.20)

Q ⁻ = diag (π) ⁻¹ (C + D) ^T diag (π) ,

c D ⁻ (l) = d(l) ⁻¹ diag (π) ⁻¹ D(l) ^{c T} diag (π) , (ν ∈ K, l ∈ Z ⁺ ), κ ^{( )} = β ^{< >} (−T ^{[ ]} ) ⁻¹

β ^{< >} (−T ^{[ ]} ) ⁻¹ e ,

Note that Q ⁻ denotes the infinitesimal generator of the time-reversed process of the underlying Markov chain that governs an arrival process and satisfies πQ ⁻ = 0. Note also that ^c D ⁻ (l) is a non-negative matrix whose row sums are all equal to or less than one, i.e., D ^{c −} (l) is a sub-stochastic matrix. Further

κ ^{( j )} exp(T ^{[ j ]} u j )(−T ^{[ j ]} )e, j = 1, 2, . . .

is considered as the density function of a phase-type distribution with representation (κ ^{( j )} , T ^{[ j ]} ).

We now define F ⁻ _k (t,~ l _k ) (k ≥ 1) as F ⁻ _k (t,~ l _k ) =

Z t

0

du _k Z u _k

0

du _k−1 · · · Z u 2

0

du ₁ Y k

j=1

κ ^{( j )} exp(T ^{[ j ]} u _j )(−T ^{[ j ]} )e

· exp Q ⁻ u ₁ ^c D ⁻ (l ₁ ) exp Q ⁻ (u ₂ − u ₁ ) ^c D ⁻ (l ₂ )

· · · · · exp Q ⁻ (u _k − u _k−1 ) ^c D ⁻ (l _k ) exp Q ⁻ (t − u _k ) . (4.21) Then F _k (t,~ l _k ) is rewritten in terms of F ⁻ _k (t,~ l _k ).

F _k (t,~ l _k ) = Y k

η=1

b(l _η )d(l _η ) · diag (π) ^{−1 h} F ⁻ _k (t,~ l _k ) ^{i T} diag (π) . (4.22) In what follows, we derive a numerically feasible formula for F ⁻ _k (t, l _k ).

We first consider F ⁻ ₁ (t, l ₁ ):

F ⁻ ₁ (t, l ₁ ) = Z t

0 du ₁

κ ^{( 1 )} exp(T ^{[ 1 ]} u ₁ )(−T ^{[ 1 ]} )e

· exp(Q ⁻ u ₁ ) ^c D ⁻ (l ₁ ) exp Q ⁻ (t − u ₁ ) . (4.23) Using the properties of Kronecker product and Kronecker sum, we rewrite (4.23) to be

F ⁻ ₁ (t, l ₁ ) = Z t

0 du ₁

κ ^{( 1 )} · exp(T ^{[ 1 ]} u ₁ ) · (−T ^{[ 1 ]} )e

·

I _Q · exp(Q ⁻ u 1 ) · D ^{c −} (l 1 )

exp Q ⁻ (t − u 1 )

= (κ ^{( 1 )} ⊗ I _Q ) Z t

0 du ₁ exp(T ^{[ 1 ]} ⊕ Q ⁻ u ₁ )

· ^h (−T ^{[ 1 ]} )e ⊗ ^c D ⁻ (l ₁ ) ⁱ exp Q ⁻ (t − u ₁ ) .

where I _Q denotes unit matrix whose size is the same as that of Q ⁻ .

Similarly, F ⁻ ₂ (t,~ l ₂ ) is rewritten to be F ⁻ ₂ (t,~ l ₂ ) =

Z t

0

du ₂ Z u 2

0

du ₁

· ^h κ ^{( 1 )} · exp(T ^{[ 1 ]} u ₁ ) · (−T ^{[ 1 ]} )e · 1 · 1 · 1 ⁱ

· ^h κ ^{( 2 )} · exp(T ^{[ 2 ]} u ₁ ) · I ₁ (l ₂ ) · exp[T ^{[ 2 ]} (u ₂ − u ₁ )] · (−T ^{[ 2 ]} )e · 1 ⁱ

· ^h I Q · exp(Q ⁻ u ₁ ) · D ^{c −} (l ₁ ) · exp[Q ⁻ (u ₂ − u ₁ )] · D ^{c −} (l ₂ )

· exp[Q ⁻ (t − u ₂ )]

= (κ ^{( 1 )} ⊗ κ ^{( 2 )} ⊗ I _Q ) Z t

0

du ₂ Z u 2

0

du ₁ exp(T ^{[ 1 ]} ⊕ T ^{[ 2 ]} ⊕ Q ⁻ u ₁ )

· (−T ^{[ 1 ]} )e ⊗ I ₁ (l ₂ ) ⊗ D ^{c −} (l ₁ ) exp[T ^{[ 2 ]} ⊕ Q ⁻ (u ₂ − u ₁ )]

· (−T ^{[ 2 ]} )e ⊗ ^c D ⁻ (l ₂ ) exp[Q ⁻ (t − u ₂ )], where I ₁ ( ~ l ₂ ) denotes unit matrix whose size is the same as that of T ^{[ 2 ]} .

Following the same manipulation as in the cases k = 1 and 2, we obtain for k = 1, 2, . . ., F ⁻ _k (t,~ l _k ) = J _k ( ~ l _k )

Z t

0

du _k Z u k

0

du _k−1 · · · Z u 2

0

du ₁

· exp[U _k,1 ( ~ l _k )u ₁ ]V _k,1 ( ~ l _k ) exp[U _k,2 ( ~ l _k )(u ₂ − u ₁ )]V _k,2 ( ~ l _k )

· · · · · exp[U _k,k ( ~ l _k )(u _k − u _k−1 )]V _k,k ( ~ l _k ) exp[Q ⁻ (t − u _k )], (4.24) where, with I _k−j (l _j+1 , l _j+2 , . . . , l _k ) (k − j ≥ 1) being an unit matrix whose size is the same as that of T ^{[ j+1 ]} ⊕ · · · ⊕ T ^{[ k ]} ,

J k ( ~ l k ) = κ ^{( 1 )} ⊗ · · · ⊗ κ ^{( k )} ⊗ I Q , (4.25)

U _k,j ( ~ l _k ) = T ^{[ j ]} ⊕ · · · ⊕ T ^{[ k ]} ⊕ Q ⁻ , j = 1, 2, . . . , k, (4.26) V _k,j ( ~ l _k ) =

( (−T ^{[ j ]} )e ⊗ I _k−j (l _j+1 , . . . , l _k ) ⊗ ^c D ⁻ (l j ), j = 1, 2, . . . , k − 1,

(−T ^{[ k ]} )e ⊗ ^c D ⁻ (l k ), j = k. (4.27)

Lemma 4.1 For j = 1, . . . , k, let U _j denote an m i × m i matrix and V _j denote an m j × m _j+1 matrix.

If all matrices U j and V j are bounded, the following equation holds for all k (k = 1, 2 . . . , ).

V ₀ Z t

0

dx k

Z x k

0

dx _k−1 · · · Z x 2

0

dx ₁ e ^{1 x 1} V ₁ e ^{2 (x 2 −x 1 )} V ₂ · · · e ^{k (t−x k )} · V k

= ^h V ₀ O · · · O ⁱ exp



 

 

 

 



 

 

 

 

U ₁ V ₁ O · · · O

O U ₂ V ₂ .. .

.. . . .. . .. O

O · · · O U _k−1 V _k−1

O O · · · O U _k



 

 

 

  t



 

 

 

 



 

 

 O

.. . O V k



 

 

 . (4.28)

The proof of Lemma 4.1 is given in Appendix. Applying Lemma 4.1 to (4.24) and taking account

of (4.16) and (4.22), we have the following theorem.

Theorem 4.3 B(t, m) is given by

B(t, m) =

| |

X

k=1

X

~ _k ∈L k ( )

Y k η=1

b(l η )d(l η )diag (π) ^{−1 h} F ⁻ _k (t,~ l k ) ^{i T} diag (π) ,

where b(l _η ) and d(l _η ) are given by (4.19) and (4.20), respectively. Further F ⁻ _k (t,~ l _k ) (k ≥ 1) is given by

F ⁻ _k (t,~ l _k ) = ^h J _k ( ~ l _k ) O · · · O ⁱ exp ^h A _k ( ~ l _k )t ⁱ



 

 

 O

.. . O I Q



 

 

 , (4.29)

with

A _k ( ~ l _k ) =



 

 

 

 



U _k,1 ( ~ l _k ) V _k,1 ( ~ l _k ) O · · · O O U _k,2 ( ~ l k ) V _k,2 ( ~ l k ) .. .

.. . . .. . .. O

O · · · O U k,k ( ~ l k ) V k,k ( ~ l k )

O O · · · O Q ⁻



 

 

 

 



, (4.30)

where J _k ( ~ l _k ), U _k,j ( ~ l _k ) (j = 1, 2, . . . , k) and V _k,j ( ~ l _k ) (j = 0, 1, . . . , k) are defined in (4.25), (4.26) and (4.27), respectively.

Remark 4.3 A _k ( ~ l _k ) in (4.30) is considered as a defective infinitesimal generator of an absorbing Markov chain. Namely, A k ( ~ l k ) has negative diagonal elements and non-negative off-diagonal elements, all row sums of it are non-positive and at least one row sum is strictly negative. Therefore applying the uniformization technique [10], we can readily compute exp[A k ( ~ l _k )t]. Further since exp[A k ( ~ l _k )t], J k ( ~ l k ), I Q , diag (π), diag (π) ⁻¹ , b(l η ) and d(l η ) are all non-negative, the computation of B(t, m) is numerically stable.

4.3 Limiting formula for phase-type services

In this subsection, assuming phase-type service times, we derive an explicit and numerically feasible formula for the limit B(m) of B(t, m):

B(m) = lim

t→∞ B(t, m).

We define F ⁻ _k ( ~ l _k ) as

F ⁻ _k ( ~ l k ) = lim

t→∞ F ⁻ _k (t,~ l k ).

Theorem 4.4 B(m) is given by

B(m) =

| |

X

k=1

X

~ k ∈ L k ( )

Y k

η=1

b(l η )d(l η ) · diag (π) ⁻¹ F ⁻ _k ( ~ l _k ) ^T diag (π) , (4.31) where

F ⁻ _k ( ~ l _k ) = J _k ( ~ l _k )(−U _k,1 ( ~ l _k )) ⁻¹ V _k,1 ( ~ l _k ) · · · · · (−U k,k ( ~ l _k )) ⁻¹ V _k,k ( ~ l _k )eπ,

with J _k ( ~ l _k ), U _k,j ( ~ l _k ) (j = 1, 2, . . . , k) and V _k,j ( ~ l _k ) (j = 0, 1, . . . , k) in (4.25), (4.26) and (4.27),

respectively.

Proof. (4.31) follows from (4.16) and (4.22). Note that A _k ( ~ l _k ) in (4.30) is rewritten to be

A _k ( ~ l _k ) =



 

 

 

 



A ⁰ _k ( ~ l _k )



 

 

 O

.. . O V _k,k ( ~ l _k )



 

 



O Q ⁻



 

 

 

 



, (4.32)

where

A ⁰ _k ( ~ l _k ) =



 

 

 

 



U _k,1 ( ~ l _k ) V _k,1 ( ~ l _k ) O · · · O O U _k,2 ( ~ l _k ) V _k,2 ( ~ l _k ) .. .

.. . . .. . .. O

O · · · O U _k,k−1 ( ~ l _k ) V _k,k−1 ( ~ l _k )

O O · · · O U _k,k ( ~ l _k )



 

 

 

 



. (4.33)

Because (see (5.10) in Appendix) exp

"

K ₁₁ K ₁₂ O K ₂₂

! t

#

=



 e ^{11 t} Z t

0

due ^{11 (t−u)} K ₁₂ e ^{22 u} O e ^{22 t}



 ,

It follows from (4.32) that

exp ^h A _k ( ~ l _k )t ⁱ =



 

exp ^h A ⁰ _k ( ~ l k )t ⁱ Θ k (t,~ l k )

O exp[Q ⁻ t]



  , (4.34)

where

Θ _k (t,~ l _k ) = Z t

0

dx exp ^h A ⁰ _k ( ~ l _k )x ⁱ



 

 

 O

.. . O V _k,k ( ~ l _k )



 

 

 exp[Q ⁻ (t − x)]. (4.35)

Note that A ⁰ _k ( ~ l _k ) given by (4.33) is regarded as the defective infinitesimal generator of an absorbing Markov chain in the transient state, and therefore all row sums of the above matrix are strictly negative. Thus we have _Z

∞ 0

exp ^h A ⁰ _k ( ~ l _k )t ⁱ dt = ^h −A ⁰ _k ( ~ l _k ) ^{i −1} . (4.36) On the other hand, since Q ⁻ is the infinitesimal generator of the irreducible Markov chain with the stationary probability vector π, we have

t→∞ lim exp(Q ⁻ t) = eπ. (4.37)

Thus from (4.35), (4.36) and (4.37), we obtain Θ _k ( ~ l _k ) = lim

t→∞ Θ _k (t,~ l _k )

= Z ∞

0 dx exp ^h A ⁰ _k ( ~ l _k )x ⁱ



 

 

 O

.. . O V _k,k ( ~ l _k )



 

 

 exp(−Q ⁻ x) · lim

t→∞ exp(Q ⁻ t)

= Z ∞

0

dx exp ^h A ⁰ _k ( ~ l k )x ⁱ ·



 

 

 O

.. . O V _k,k ( ~ l _k )



 

 

 eπ

= ^h −A ⁰ _k ( ~ l _k ) ^{i −1}



 

 

 O

.. . O V k,k ( ~ l k )



 

 

 eπ, (4.38)

where we use exp(Q ⁻ t)e = e (t ≥ 0). Therefore

t→∞ lim exp ^h A k ( ~ l k )t ⁱ =

"

O Θ k ( ~ l _k ) O eπ

#

. (4.39)

Finally, from (4.29), (4.38) and (4.39), we obtain F ⁻ _k ( ~ l k ) = lim

t→∞ F ⁻ _k (t,~ l k )

= ^h J _k ( ~ l _k ) O · · · O O ⁱ

"

O Θ _k ( ~ l _k ) O eπ

# "

O I Q

#

= ^h J _k ( ~ l _k ) O · · · O ^{i h} −A ⁰ _k ( ~ l _k ) ^{i −1}



 

 

 O

.. . O V k,k ( ~ l k )



 

 

 eπ,

from which Theorem 4.4 follows. 2

5 Numerical Examples

In this section, we show some numerical examples based on Theorems 4.3 and 4.4, and discuss the impact of system parameters on the mean and variance of the number of customers in the system under the assumption that the arrival process is stationary. Note that when the arrival process is stationary, the mean number E[N ν (t)] (ν ∈ K) of class ν customers at time t is given by

E[N ν (t)] = πB(t, e ν )e = λ ν

Z t

0

dxH ν (x), where e _ν is given by (2.1). Also the limit E[N ν ] of E[N ν (t)] is given by

E[N _ν ] = πB(e _ν )e = λ _ν b _ν , where b ν denotes the mean service time of class ν customers.

b ν = Z ∞

0 du H ν (u).

Therefore the mean E[N (t)] of the total number of customers at time t and its limit E[N ] are given by

E[N (t)] = ^X

ν∈K

λ _ν Z t

0

dxH _ν (x), E[N ] = ^X

ν∈K

λ _ν b _ν .

On the other hand, the variance Var[N _ν (t)] (ν ∈ K) of the number of class ν customers at time t and its limit Var[N ν ] are given by

Var[N ν (t)] = 2πB(t, 2e ν )e + E[N ν (t)] − E[N ν (t)] ² , Var[N ν ] = 2πB(2e ν )e + E[N ν ] − E[N ν ] ² .

Also, the covariance Cov[N _ν (t), N _ν ⁰ (t)] (ν, ν ⁰ ∈ K) of the numbers of class ν and ν ⁰ customers at time t and its limit Cov[N ν , N ν ⁰ ] are given by

Cov[N ν (t), N _ν ⁰ (t)] = πB(t, e ν + e _ν ⁰ )e − E[N ν (t)]E[N _ν ⁰ (t)], Cov[N ν , N ν ⁰ ] = πB(e ν + e ν ⁰ )e − E[N ν ]E[N ν ⁰ ].

Therefore the variance Var[N (t)] of the total number of customers at time t and its limit Var[N ] are given by

Var[N (t)] = ^X

ν∈K

Var[N _ν (t)] + 2 ^X

ν,ν0∈K

ν6=ν ⁰

Cov[N _ν (t), N _ν ⁰ (t)],

Var[N ] = ^X

ν∈K

Var[N ν ] + 2 ^X

ν,ν0 ∈K

ν6=ν ⁰

Cov[N ν , N ν ⁰ ].

5.1 Impact of service time distribution on Var[ N ]

We first show the impact of the service time distribution on the limiting variance of the number of customers. To do so, we consider the following MMPP input. There is only one class, i.e., K = {1}, and

C =

"

−9 1

1 −3

#

, D =

"

8 0 0 2

# ,

D(n) = θ(n)D, θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise.

As for the service time distribution, we fix the mean service time to one (i.e., E[N ] = 10), and consider the followings. For 0 < C _v < 1, k-stage Erlang distributions:

H ₁ (t) = ^h 1 0 _{| {z }} . . . 0

k−1

i exp



 

 

 





 

 

 



−k k 0 · · · 0

0 −k k · · · 0

.. . . .. ... .. .

0 0 · · · −k k

0 0 · · · 0 −k



 

 

 





 

 

 





 

 

 

 1 1 .. . 1 1



 

 

 



k = 2, 3, . . . , (5.1)

for C v = 1, the exponential distribution:

H 1 (t) = e ^−t , (5.2)

and for C v > 1, a 2-state balanced hyper-exponential distribution:

H ₁ (t) = ^h p 1 − p ⁱ exp

"

−2p 0

0 −2(1 − p)

!# "

1 1

#

, 0 < p < 0.5. (5.3) Figure 1 plots the limiting variance Var[N] as a function of the coefficient of variation C v of the service time distribution. We observe that as C v increases, the variance of the number of customers in the system decreases.

15 20 25 30 35 40 45 50

0 5 10 15 20 25 30

Var[N ]

coefficeint of variation C v

Figure 1: Limiting variance of the number of customers.

Remark 5.1 For the M (t)/G/∞ queue with a sinusoidal arrival rate, Eick et al. [3] show that as C v

decreases, the deviation of the arrival process has greater impact on the mean number of customers in the system. Thus our observation coincides with that in [3].

5.2 Impact of arrival process

We consider the impact of the correlation in the arrival process on the variance of the number of customers. For this purpose, we assume that batches arrive according to the following 2-state Markov modulated Poisson process:

C =

"

−8 − c c

c −2 − c

#

, D =

"

8 0 0 2

#

, c > 0,

D(n) = θ(n)D, θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise.

Note that the time-correlation of the arrival process increases with the mean sojourn time 1/c in each state of the underlying Markov chain. As for the service time distribution, we consider three cases:

(1) the 2-stage Erlang distribution in (5.1) with k = 2, (2) the exponential distribution in (5.2), and (3) the 2-stage hyper-exponential distribution in (5.3) with p = 0.25.

Figure 2 plots the limiting variance Var[N ] as a function of 1/c. We observe that the limiting

variance increases with the increase of correlation in arrivals.

0 10 20 30 40 50 60

0 5 10 15 20 25 30 35 40 45 50

Var[N ]

1/c

Erlang Exp.

Hyper-exp.

Figure 2: Limiting variance of the number of customers.

Next we consider the time-dependent variance Var[N (t)] in systems with two stationary arrival streams. We fix the marginal characteristics of each arrival streams, which is represented by

C =

"

−6 1

1 −1

#

, D =

"

5 0 0 0

# ,

θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise,

D(n) = θ(n)D,

H(t) = ^h 0.5 0.5 ⁱ exp

"

−0.75 0

0 −1.5

! t

# "

1 1

# .

where H(t) denotes the complementary distribution of service times. Under this restriction, we con- sider the following three cases.

Case 1: Negatively correlated arrival streams.

The two arrival streams are negatively correlated. Namely, two arrival streams are governed by the single two-state underlying Markov chain and when it is in state ν (ν = 1, 2), only class ν customers can arrive. More precisely, the overall arrival process is characterized by

C =

"

−6 1

1 −6

#

, D ₁ =

"

5 0 0 0

#

, D ₂ =

"

0 0 0 5

# ,

θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise.

D(n ₁ , n ₂ ) =

 



 

θ(n 1 )D 1 , if n 1 ≥ 1 and n 2 = 0, θ(n ₂ )D ₂ , if n ₁ = 0 and n ₂ ≥ 1,

O otherwise,

H ν (t) = ^h 0.5 0.5 ⁱ exp

"

−0.75 0

0 −1.5

! t

# "

1 1

#

, ν = 1, 2.

Case 2: Superposition of two independent arrival streams.

The two arrival streams are independent of each other. Thus the overall arrival process is charac- terized by

C =

"

−6 1

1 −1

#

⊕

"

−6 1

1 −1

#

=



 

 

−12 1 1 0

1 −7 0 1

1 0 −7 1

0 1 1 −2



 

  ,

D ₁ =

"

5 0 0 0

#

⊗

"

1 0 0 1

#

=



 

 

5 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0



 

  ,

D ₂ =

"

1 0 0 1

#

⊗

"

5 0 0 0

#

=



 

 

5 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0



 

  ,

θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise,

D(n ₁ , n ₂ ) =

 



 

θ(n ₁ )D ₁ , if n ₁ ≥ 1 and n ₂ = 0, θ(n ₂ )D ₂ , if n ₁ = 0 and n ₂ ≥ 1,

O, otherwise,

H ν (t) = ^h 0.5 0.5 ⁱ exp

"

−0.75 0

0 −1.5

! t

# "

1 1

#

, ν = 1, 2.

Case 3: Positively correlated arrival streams.

The two arrival streams are positively correlated. Namely, arrivals from both streams occur simul- taneously with probability one. More precisely, the overall arrival process is characterized by

C =

"

−6 1

1 −1

#

, D =

"

5 0 0 0

# ,

θ(n) =

 



 

1/2, if n = 1, 1/2, if n = 3, 0, otherwise,

D(n ₁ , n ₂ ) = θ(n ₁ )θ(n ₂ )D,

H ν (t) = ^h 0.5 0.5 ⁱ exp

"

−0.75 0

0 −1.5

! t

# "

1 1

#

, ν = 1, 2.

Because the marginal characteristics of each stream in these three cases are the same, the time- dependent and the limiting marginal distributions of the numbers of customers of respective class are identical among these three cases. However, the covariance of the numbers of customers in respective classes are different. Figure 3 show the time-dependent covariance Cov[N 1 (t), n 2 (t)] of the numbers of customers in respective class as a function of t. We observe that positive (resp. negative) correlation between two arrival streams leads to positive (resp. negative) covariance in any time t. Thus the positive (resp. negative) correlation in arrivals leads to large (resp. small) variance, compared with the superposition of independent streams, i.e., Case 2, as shown in Figure 4.

5.3 Impact of correlation in service time sequence

Finally we consider the M ^X /SM/∞ queue to investigate the correlation in the service time sequence.

We assume that batches arrive in a Poisson process with rate five, and batch sizes are i.i.d. according

-15 -10 -5 0 5 10 15 20 25 30

0 5 10 15 20

Cov[N (t)]

t

Case 3 Case 2 Case 1

Figure 3: Time-dependent covariance of the number of customers.

0 10 20 30 40 50 60 70 80 90

0 5 10 15 20

Var[N (t)]

t

Case 3 Case 2 Case 1

Figure 4: Time-dependent variance of the number of customers