Structured Markov Chains Ivo Adan and Johan van Leeuwaarden Where - - PowerPoint PPT Presentation

Where innovation starts

Structured Markov Chains

Ivo Adan and Johan van Leeuwaarden

2

Book on Analysis of structured Markov processes (arXiv:1709.09060)

I Basic methods – Basic Markov processes – Advanced Markov processes – Queues and transforms II Basic processes – Birth-and-death processes – Queueing networks – Quasi-birth-and-death-processes – Processes that are quasi-skip-free in one direction III Advanced processes – Single-server priority systems – Gated single-server systems – Production systems – Join the shortest queue

3

Markov chains in discrete time

A stochastic process {Xn, n = 0, 1, 2, ...} with state space {0, 1, 2, ...} is a Markov chain if P(Xn+1 = in+1|X0 = i0, X1 = i1, ... , Xn = in) = P(Xn+1 = in+1|Xn = in) for all i0, ... , in+1 and n ≥ 0

• Assumptions:

– Time homogeneous: P(Xn+1 = j|Xn = i) = pij – Irreducible: all states are reachable – Aperiodic

• Equilibrium distribution p0, p1, ... satisfy the equilibrium (or balance) equations

pi =

∞

• j=0

pjpji, i = 0, 1, ...

• If equilibrium distribution exists then it is equal to the limiting distribution limn→∞ P(Xn = i) = pi
• Vector-matrix notation: p = pP where p = (p0, p1, ...) and transition probability matrix P = (pij)

4

Markov chains in continuous time

A stochastic process {X(t), t ≥ 0} with state space {0, 1, 2, ...} is a Markov chain if P(X(t + s) = j|X(u) = x(u), 0 ≤ u < s, X(s) = i) = P(X(t + s) = j|X(s) = i) for all i, j, x(u) and s, t ≥ 0

• Assumptions:

– Time homogeneous: P(X(t + s) = j|X(s) = i) = pij(t) – Irreducible: all states are reachable

• Transition rates qij = limh→0

pij(h) h

where we let qii = 0 and qi = ∞

j=0 qij

• Equilibrium distribution p0, p1, ... satisfy the equilibrium (or balance) equations

piqi =

∞

• j=0

pjqji, i = 0, 1, ...

• If equilibrium distribution exists then it is equal to the limiting distribution limt→∞ P(X(t) = i) = pi
• Vector-matrix notation: 0 = pQ where p = (p0, p1, ...) and generator Q with Qij = qij, i = j, Qii = −qi

5

Balance equations

• Global balance equations:

piqi =

∞

• j=0

pjqji

• r for subset A and its complement Ac
• i∈A

piqij =

• j∈Ac

pjqji

• Partial balance equations:

pi

• j∈A

qij =

• j∈B

pjqji for subsets A and B

• Detailed (or local) balance equations:

piqij = pjqji

6

Balance equations

• Global balance equations:
• Partial balance equations:
• Detailed (or local) balance equations:

7

Transforms

• Generating function of non-negative discrete random variable X with probability distribution pi = P(X = i)

PX(z) = E(zX) =

∞

• i=0

pizi, |z| ≤ 1

• Laplace-Stieltjes transform a non-negative (continuous) random variable X with distribution function FX(x)

X(s) = E(e−sX) = ∞ e−sxdFX(x), Re(s) ≥ 0

• Properties:

– Probability distribution pi is uniquely determined by generating function PX(z) – Distribution function FX(x) is uniquely determined by Laplace-Stieltjes X(s) – E(X(X − 1) · · · (X − n + 1)) = P(n)

X (1) and E(X n) = (−1)nX (n)(0) where (n) is nth derivative

8

Transforms: Example

• Markov chain {Xn, n = 0, 1, 2, ...} with state space {0, 1, 2, ... , n + 1} and transition probabilities pij
• State n + 1 is single absorbing state

1 2 3 p11 = 1

2

p12 = 1

2

p23 = 1

3

p21 = 2

3

p33 = 1

• Random variable Ni is number of visits to state i before absorption into n + 1
• Time spent in state i during the kth visit is equal to αk−1

i

where 0 < αi < 1 (learning curve)

• Random variable T is total time spent in states {1, ... , n} before absorption in n + 1
• What is E(T|X0 = 1)?

9

Transforms: Example

• Note

T =

n

• i=1
• 1 + αi + · · · + αNi−1

i

• =

n

• i=1

1 − αNi

i

1 − αi so E(T|X0 = 1) =

n

• i=1

1 − E(αNi

i |X0 = 1)

1 − αi

• Introduce for i = 1, ... , n the generating function

Pi(z) = E(zN1

1 · · · zNn n |X0 = i),

|z1| ≤ 1, ... , |zn| ≤ 1, where z = (z1, ... , zn) so P1(1, ... , 1, αi, 1, ... , 1) = E(αNi

i |X0 = 1)

10

Transforms: Example

• One-step analysis:

Pi(z) = E(zN1

1 · · · zNn n |X0 = i)

=

n+1

• j=1

pijE(zN1

1 · · · z1+Ni i

· · · zNn

n |X0 = j)

=

n+1

• j=1

zipijE(zN1

1 · · · zNi i

· · · zNn

n |X0 = j)

=

n+1

• j=1

zipijPj(z), i = 1, ... , n where Pn+1(z) = 1 (since N1 = · · · = Nn = 0 if X0 = n + 1)

• Hence P1(1, ... , 1, αi, 1, ... , 1) can be solved from above equations with z = (1, ... , 1, αi, 1, ... , 1)

11

Finite Markov chains (no structure)

• Markov process {X(t), t ≥ 0} with state space {0, 1, ... , N} and transition rates qij
• Global balance equations

pi

N

• j=0

qij =

N

• j=0

pjqji, i = 0, 1, ... , N

12

Finite Markov chains (no structure)

i j N qiN qNj ¯ qij

• Numerical solution through embedding:

– Remove state N by embedding Markov chain on {0, 1, ... , N − 1} with transition rates ¯ qij = qij + qiN qNj qN0 + · · · + qNN−1 , i, j = 0, ..., N − 1 – Equilibrium probability pN follows from global balance pN

N−1

• j=0

qNj =

N−1

• j=0

pjqjN – Repeat this for N − 1, N − 2, ... , 1 and finally normalize p0, ... , pN

• This is Gaussian elimination with multiplication, division and addition of positive numbers

13

Infinite Markov chains (with structure): Example of impatient customers

µ λ θ

• Single server
• Poisson arrivals with rate λ
• Exponential service times with rate µ
• Exponential patience with rate θ: customer leaves if service does not start before patience expires
• Markov chain with states i where i is number in system
• Transition rates qi,i−1 = µ + (i − 1)θ and qi−1,i = λ

1 i − 1 i λ λ µ + (i − 1)θ µ · · · · · ·

14

Infinite Markov chains (with structure): Birth-and-Death processes

1 i − 1 i λ0 λi−1 µi µ1 · · · · · ·

• Markov chain on {0, 1, ...} with transitions to neighboring states with rates qi,i−1 = µi and qi−1,i = λi−1
• Global balance

p0λ0 = p1µ1 pi(λi + µi) = pi−1λi−1 + pi+1µi+1, i = 1, 2, ...

• Global balance for {0, 1, ... , i − 1}

pi−1λi−1 = piµi, i = 1, 2, ...

15

Infinite Markov chains (with structure): Birth-and-Death processes

1 i − 1 i λ0 λi−1 µi µ1 · · · · · ·

• Iterating yields

pi = pi−1 λi−1 µi = pi−2 λi−2λi−1 µi−1µi = · · · = p0 λ0 · · · λi−1 µ1 ... µi

• Normalization

p0 = 1 ∞

i=0

i

j=1 λj−1 µj

• Assumption

∞

• i=0

i

• j=1

λj−1 µj < ∞

• Example of impatient customers

pi = p0 λi µ(µ + θ) · · · (µ + (i − 1)θ)

16

Infinite Markov chains (with more structure): Constant birth and death rates

1 i − 1 i λ λ µ µ · · · · · ·

• Global balance equations

p0λ = p1µ pi(λ + µ) = pi−1λ + pi+1µ, i = 1, 2, ...

• Assume λ − µ < 0 and (without loss of generality) λ + µ = 1
• Generator Q of the form

Q =       −λ λ ... µ −(λ + µ) λ µ −(λ + µ) λ µ −(λ + µ) . . . ...      

17

Infinite Markov chains (with more structure): Constant birth and death rates

1 i − 1 i λ λ µ µ · · · · · ·

• Solution: Generating function P(z) = ∞

i=1 pizi

– Multiply balance equations by zi and add all equations P(z)(λ + µ) − p0µ = P(z)λz + (P(z) − p0)µz−1 – Solving this equation P(z) = p0µ(1 − z−1) (λ + µ) − λz − µz−1 = p0µ(1 − z−1) (1 − z−1)(µ − λz) = p0 1 − ρz where ρ = λ/µ – Normalization P(1) = 1 so p0 = 1 − ρ and P(z) = 1 − ρ 1 − ρz =

∞

• i=0

(1 − ρ)ρizi

18

Infinite Markov chains (with more structure): Constant birth and death rates

1 i − 1 i λ λ µ µ · · · · · ·

• Solution: Geometric

– Ratio pi+1/pi is expected time spent in i + 1 before first return to i, given initial state is i – Transition structure implies that ratio pi+1/pi does not depend on i pi+1 pi = α – Hence pi = p0αi and substitution in global balance equation yields αi(λ + µ) = αi−1λ + αi+1µ – Dividing by common power αi−1 α(λ + µ) = λ + α2µ – This leads to α = λ

µ = ρ

19

Infinite Markov chains (with more structure): Constant birth and death rates

1 i − 1 i λ λ µ µ · · · · · ·

• Solution: Difference equation

– Balance equations are second-order difference equation with constant coefficients – Solution is linear combination of powers pi = αi – Substitution of pi = αi in balance equation yields α(λ + µ) = λ + α2µ with roots α1 = ρ and α2 = 1 – General solution pi = c1αi

1 + c2αi 2 = c1ρi + c2

where c1 and c2 are constants – Convergence of pi implies c2 = 0 and normalization yields c1 = 1 − ρ

20

Infinite Markov chains: M/G/1 structure

1 i − 1 i · · · i + 1 i + 2 a0 a1 a2 a3 b0 b1

• Transition probabilities P = (pij)

P =       b0 b1 b2 b3 · · · a0 a1 a2 a3 a0 a1 a2 a0 a1 . . . ...      

• Example: Number in system just after departure in M/G/1 queue
• Let A(z) = ∞

i=0 aizi and B(z) = ∞ i=0 bizi and assume mean step size = A(1)(1) − 1 < 0

• Global balance equations

pi = pi+1a0 + pia1 + · · · + p1ai + p0bi, i = 0, 1, 2, ...

21

Infinite Markov chains: M/G/1 structure

1 i − 1 i · · · i + 1 i + 2 a0 a1 a2 a3 b0 b1

• Solution: Recursion

– Global balance for {0, 1, ... , i − 1} p0¯ bi + p1¯ ai−1 + · · · + pi−1¯ a1 = pia0, i = 1, 2, ... where ¯ bi =

∞

• j=i

bi, ¯ ai =

∞

• j=i

ai – Recursively calculate p1, p2, ... starting with p0 – p0 follows from (mean displacement is 0 in equilibrium) p0B(1)(1) + (1 − p0)(A(1)(1) − 1) = 0

22

Infinite Markov chains: M/G/1 structure

1 i − 1 i · · · i + 1 i + 2 a0 a1 a2 a3 b0 b1

• Solution: Generating function

– Let P(z) = ∞

i=0 pizi, A(z) = ∞ i=0 aizi, B(z) = ∞ i=0 bizi

– Multiply balance equations by zi and add all equations P(z) = z−1(P(z) − p0)A(z) + p0B(z) – Solving this equation P(z) = p0(B(z) − z−1A(z)) 1 − z−1A(z) = p0(zB(z) − A(z)) z − A(z) – Normalization P(1) = 1 so p0 = 1 − A(1)(1) B(1)(1) + 1 − A(1)(1)

23

Infinite Markov chains: G/M/1 structure

1 i − 1 i · · · i + 1 a2 a1 a0 i − 2 a3 ai bi · · ·

• Transition probabilities P = (pij)

P =       b0 a0 · · · b1 a1 a0 b2 a2 a1 a0 b3 a3 a2 a1 . . . ...       , bi = 1 − (a0 + · · · + ai) =

∞

• j=i+1

ai

• Example: Number in system just before arrival in G/M/1 queue
• Global balance equations

pi = pi−1a0 + pia1 + pi+1a2 + · · · =

∞

• j=0

pi−1+jaj, i = 1, 2, ...

24

Infinite Markov chains: G/M/1 structure

1 i − 1 i · · · i + 1 a2 a1 a0 i − 2 a3 ai bi · · ·

• Solution: Geometric

– Ratio pi+1/pi is expected time spent in i + 1 before first return to i, given initial state is i – Transition structure implies that ratio pi+1/pi does not depend on i pi+1 pi = α – Hence pi = p0αi and substitution in global balance equation yields α =

∞

• j=0

αjaj = A(α) – α is the unique root on (0, 1) of α = A(α) provided mean step size 1 − A(1)(1) < 0

25

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Single server
• Poisson arrivals with rate λ wich need exponential service times with rate µ
• Server switched off when system is empty
• Server switched on when first customer arrives which requires exponential setup time with rate θ

26

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Markov chain with states (i, j) with i number in system and j status of server: 0 means off, 1 means on
• Transition rate diagram

0, 0 1, 0 1, 1 2, 1 2, 0 i, 0 i, 1 λ λ λ λ λ µ µ θ θ θ · · · · · · 0, 1 λ µ

• Global balance equations

p(0, 0)λ = p(1, 1)µ p(0, 1)λ = 0 p(i, 0)(λ + θ) = p(i − 1, 0)λ p(i, 1)(λ + µ) = p(i − 1, 1)λ + p(i, 0)θ + p(i + 1, 1)µ, i = 1, 2, ...

27

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Order states as (0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1), ... and partition into levels i = {(i, 0), (i, 1)}
• Generator Q of the form

Q =       B1 A0 ... B2 A1 A0 A2 A1 A0 A2 A1 . . . ...       where B1 = −λ −λ

• ,

B2 = 0 0 µ 0

• ,

A0 = λ 0 0 λ

• ,

A1 = −(λ + θ) θ −(λ + µ)

• ,

A2 = 0 0 0 µ

• Let p = (p0, p1, ...) and pi = (p(i, 0), p(i, 1))

p0B1 + p1B2 = 0, pi−1A0 + piA1 + pi+1A2 = 0, i = 1, 2, ...

28

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Matrix-geometric

– Global balance for level 0, 1, ... , i p(i + 1, 1)µ = (p(i, 0) + p(i, 1))λ

• r in matrix form

pi+1A2 = piA3 with A3 = 0 λ 0 λ

• – Substituting in balance equation

0 = pi−1A0 + piA1 + pi+1A2 = pi−1A0 + pi(A1 + A3) – Hence pi = −pi−1A0(A1 + A3)−1 = pi−1R with R = −A0(A1 + A3)−1 = λ/(λ + θ) λ/µ λ/µ

29

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Matrix-geometric

– Iterating yields pi = p0Ri, i = 0, 1, 2, ... – Vector p0 follows from boundary equation p0B1 + p1B2 = p0(B1 + RB2) = 0 and normalization 1 =

• i,j

p(i, j) =

∞

• i=0

pie = p0

∞

• i=0

Rie = p0(I − R)−1e where I is identity matrix and e all-one vector

30

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Matrix-analytic

– Let gjk be probability of first passage to level i − 1 in (i − 1, k) for initial state (i, j) at level i G = (gjk) = 0 1 0 1

• – Then

pi+1A2 = piA0G – Substituting in balance equation 0 = pi−1A0 + piA1 + pi+1A2 = pi−1A0 + pi(A1 + A0G) – Hence pi = −pi−1A0(A1 + A0G)−1 = pi−1R

31

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Spectral expansion

– Seek solutions of balance equations pi−1A0 + piA1 + pi+1A2 = 0 of the form pi = y · xi, i = 0, 1, 2, ... , where y = (y(0), y(1)) is non-null vector and |x| < 1 – Substitution and dividing by common powers of x yields y

• A0 + xA1 + x2A2
• = 0

– Desired values of x are roots inside the unit circle of det(A0 + xA1 + x2A2) = (λ − (λ + θ)x)(µx − λ)(x − 1) = 0

32

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Spectral expansion

– Roots x1 = λ λ + θ , x2 = λ µ with corresponding non-null vectors y1 = (1, −θx1 λ − (λ + µ)x1 + µx2

1

), y2 = (0, 1) – Set pi = c1y1xi

1 + c2y2xi 2

where coefficients c1 and c2 follow from boundary equation p0B1 + p1B2 = 0 and normalization 1 =

• i,j

p(i, j) = c1y1e 1 − x1 + c2y2e 1 − x2

33

Examples of Quasi-Birth-and-Death processes: Single server with setup

• Solution: Generating function

– Let P(z) = ∞

i=0 pizi

– Multiply balance equations by zi and add all equations zP(z)A0 + (P(z) − p0)A1 + z−1(P(z) − p0 − p1z)A2 = 0 – Solving this equation P(z)A(z) = p0(A1 + z−1A2) + p1A2 where A(z) = zA0 + A1 + z−1A2 = λz − (λ + θ) θ λz − (λ + µ) + µz−1

• – Vectors p0, p1 and P(1) follow from

p0B1 + p1B2 = 0 P(1)A(1) = p0(A1 + A2) + p1A2, P(1)e = 1 P(1)A(1)(1)e = −p0A2e

34

Examples of Quasi-Birth-and-Death processes: Shortest queue with jockeying

• Two parallel exponential servers with rate µ, each with its own queue
• Customers arrive according to Poisson stream with rate λ and join shortest queue on arrival
• Jockeying: If difference in queues exceeds 1, then a customer jumps from longest to shortest queue

35

Examples of Quasi-Birth-and-Death processes: Shortest queue with jockeying

• Markov chain with states (i, j) with i number in shortest queue and j difference between queues
• Transition rate diagram

0, 0 1, 0 1, 1 · · · · · · 0, 1 i, 1 i, 0 λ 2µ λ λ µ λ 2µ i+1, 1 i+1, 0

36

Examples of Quasi-Birth-and-Death processes: Shortest queue with jockeying

• Partition into levels i = {(i, 0), (i, 1)} for i = 0, 1, ...
• Generator Q of the form

Q =       B1 A0 ... A2 A1 A0 A2 A1 A0 A2 A1 . . . ...       where B1 = −λ λ µ −(λ + µ)

• ,

A0 = 0 0 λ 0

• ,

A1 = −(λ + 2µ) λ 2µ −(λ + 2µ)

• ,

A2 = 0 2µ

• Let p = (p0, p1, ...) and pi = (p(i, 0), p(i, 1))

p0B1 + p1A2 = 0 pi−1A0 + piA1 + pi+1A2 = 0, i = 1, 2, ...

37

Examples of Quasi-Birth-and-Death processes: Single server serving n customer types

• Single server serving n types of customers, numbered 1, ... , n
• Type i customers arrive according to Poisson stream with rate λi
• Type i customers require exponential service with rate µi

38

Examples of Quasi-Birth-and-Death processes: Single server serving n customer types

• Markov chain with states (i, j) where i number waiting in the queue and j is customer type in service
• State (0, 0) is empty state
• Transition rate diagram where λ = λ1 + · · · + λn and pi = λi/λ

0, 0 λ λ · · · · · · . . . 0, 1 0, n 1, 1 1, n i, n i, 1 λ µnp1 µnp1 µnpn µnpn µ1p1 µ1p1 µ1 λ1

39

Examples of Quasi-Birth-and-Death processes: Single server serving n customer types

• Partition into level 0 = {(0, 0), ... , (0, n)} and levels i = {(i, 1), ... , (i, n)} for i = 1, 2, ...
• Generator Q of the form

Q =       B1 B0 ... B2 A1 A0 A2 A1 A0 A2 A1 . . . ...       where A0 =   λ ... λ   , A1 =   −(λ + µ1) ... −(λ + µn)   , A2 =   µ1p1 · · · µ1pn . . . . . . µnp1 · · · µnpn  

• Let p = (p0, p1, ...), p0 = (p(0, 0), ... , p(0, n)) and pi = (p(i, 1), ... , p(i, n))

p0B1 + p1B2 = 0 p0B0 + p1A1 + p2A2 = 0 pi−1A0 + piA1 + pi+1A2 = 0, i = 2, 3, ...

40

Quasi-Birth-and-Death processes

• Markov chain with states (0, 0), (0, 1), ... , (0, n), (1, 0), (1, 1), ... , (1, n), (2, 0), ...
• Partition into levels i = {(i, 0), ... , (i, n)} for i = 0, 1, ...
• Generator Q of the form

Q =       B1 A0 ... A2 A1 A0 A2 A1 A0 A2 A1 . . . ...       where B1, A0, A1 and A2 are n + 1 × n + 1 matrices

• Let p = (p0, p1, ...) and pi = (p(i, 0), p(i, 1), ... , p(i, n))

p0B0 + p1A2 = 0 pi−1A0 + piA1 + pi+1A2 = 0, i = 1, 2, ...

41

Quasi-Birth-and-Death processes

• Solution: Matrix-geometric

– Probability vector p is given by pi = p0Ri, i = 0, 1, 2, ... – Rate matrix R is minimal nonnegative solution of A0 + RA1 + R2A2 = 0 – Rjk is expected time spent in (i + 1, k) before first return to level i, given initial state is (i, j) – Vector p0 follows from boundary equation and normalization p0(B1 + RA2) = 0 p0(I − R)−1e = 1 – Rate matrix R can be determined by successive substitution from fixed point equation R = −(A0 + R2A2)A−1

1

starting with R = 0

42

Quasi-Birth-and-Death processes

• Solution: Matrix-analytic

– Let gjk be probability of first passage to level i − 1 in (i − 1, k) for initial state (i, j) at level i – Matrix G = (gjk) is is minimal nonnegative solution of A2 + A1G + A0G2 = 0 – Then pi+1A2 = piA0G – Substituting in balance equation 0 = pi−1A0 + piA1 + pi+1A2 = pi−1A0 + pi(A1 + A0G) – Probability vector p is given by pi = −pi−1A0(A1 + A0G)−1 = pi−1R with R = −A0(A1 + A0G)−1

43

Quasi-Birth-and-Death processes

• Solution: Matrix-analytic

– Vector p0 follows from boundary equation and normalization p0(B0 + RA2) = p0(B0 + A0G) = 0 p0(I − R)−1e = 1 – Note RA2 = A0G – Matrix G can be determined by successive substitution from fixed point equation G = −(A2 + A0G2)A−1

1

starting with G = 0

44

Quasi-Birth-and-Death processes

• Solution: Spectral expansion

– Seek solutions of balance equations pi−1A0 + piA1 + pi+1A2 = 0 of the form pi = y · xi, i = 0, 1, 2, ... , where y = (y(0), ... , y(n)) is non-null vector and |x| < 1 – Substitution and dividing by common powers of x yields y

• A0 + xA1 + x2A2
• = 0

– Suppose determinant equation det(A0 + xA1 + x2A2) = 0 has exactly n + 1 roots x0, ... , xn with |xk| < 1 – Let y0, ... , yn be the corresponding non-null vectors yk

• A0 + xkA1 + x2

kA2

• = 0,

k = 0, ... , n

45

Quasi-Birth-and-Death processes

• Solution: Spectral expansion

– Set pi =

n

• k=0

ckykxi

k

– Coefficients c0, ... , ck follow from boundary equation p0B1 + p1A2 =

n

• k=0

ck (ykB1 + ykA2xk) = 0 and normalization 1 =

• i,j

p(i, j) =

n

• k=0

ckyke 1 − xk

46

Quasi-Birth-and-Death processes

• Solution: Generating function

– Let P(z) = ∞

i=0 pizi

– Multiply balance equations by zi and add all equations zP(z)A0 + (P(z) − p0)A1 + z−1(P(z) − p0 − p1z)A2 = 0 – Solving this equation P(z)A(z) = p0(A1 + z−1A2) + p1A2 = p0(A1 − B1 + z−1A2) where A(z) = zA0 + A1 + z−1A2

47

Quasi-Birth-and-Death processes

• Solution: Generating function

– Suppose determinant equation det(A(z)) = 0 has exactly n + 1 roots z0 = 1, z1, ... , zn with |zk| < 1 for k = 1, ... , n – Let v0 = e, v1, ... , vn be non-null vectors satisfying A(zk)vk = 0 – Vector p0 follows from p0(A1 − B1 + z−1

k A2)vk = 0,

k = 1, ... , n P(1)A(1) = p0(A1 − B1 + A2), P(1)e = 1 P(1)A(1)(1)e = −p0A2e

48

Quasi-Birth-and-Death processes: More structure

• A0, A1 and A2 are upper triangular: A0(i, j) = 0 if j < i
• Implication: R is also upper triangular
• Example:

A0 = a0 a1 a2

• ,

A1 = b0 b1 b2

• ,

A2 = c0 c1 c2

• ,
• Then

R = r0 r1 0 r2

• where

a0 + r0b0 + r2

0c0 = 0

a2 + r2b2 + r2

2c2 = 0

a1 + r0b1 + r1b2 + r2

0c1 + (r0r1 + r1r2)c2 = 0

49

Quasi-Birth-and-Death processes: More structure

• Transition rates A2 from level i to level i − 1 have the form

A2 = rp =   r0 . . . rn   p0 · · · pn

• =

  r0p0 · · · r0pn . . . . . . rnp0 · · · rnpn  

• Example: Single server with setup

A2 = 0 0 0 µ

• =

µ 0 1

• Example: Shortest queue with jockeying

A2 = 0 2µ

• =

2µ 0 1

• Example: Single server serving n customer types for which

A2 =   µ1p1 · · · µ1pn . . . . . . µnp1 · · · µnpn   =   µ1 . . . µn   p1 · · · pn

50

Quasi-Birth-and-Death processes: More structure

• Transition rates A2 from level i to level i − 1 have the form

A2 = rp =   r0 . . . rn   p0 · · · pn

• =

  r0p0 · · · r0pn . . . . . . rnp0 · · · rnpn  

• Interpretation: Rate of jumping from (i, j) to level i − 1 is rj, this jump is to (i − 1, k) with probability pk
• Important: Jump probability does not depend on j
• This implies

G =   p0 · · · pn . . . . . . p0 · · · pn   = ep and R = −A0(A1 + A0G)−1 = −A0(A1 + A0ep)−1

51

Quasi-Birth-and-Death processes: More structure

• Transition rates A0 from level i to level i + 1 have the form

A0 = rp =   r0 . . . rn   p0 · · · pn

• =

  r0p0 · · · r0pn . . . . . . rnp0 · · · rnpn  

• Example: Shortest queue with jockeying

A0 = 0 0 λ 0

• =

λ 1 0

52

Quasi-Birth-and-Death processes: More structure

• Transition rates A0 from level i to level i + 1 have the form

A0 = rp =   r0 . . . rn   p0 · · · pn

• =

  r0p0 · · · r0pn . . . . . . rnp0 · · · rnpn  

• Interpretation: Rate of jumping from (i, j) to level i + 1 is rj, this jump is to (i + 1, k) with probability pk
• This implies

R = rq =   r0 . . . rn   q0 · · · qn

• =

  r0q0 · · · r0qn . . . . . . rnq0 · · · rnqn   for some vector q

• Hence

Ri = (qr)i−1R, pi = p0Ri = (qr)i−1p1, i = 1, 2, ...

• Number x = qr is spectral radius of R which is unique root on (0, 1) of

det(A0 + xA1 + x2A2) = 0

53

Quasi-Birth-and-Death processes: Linear rates

• Transition rates have linear form to neighboring phases

– Rate from (i, j) to (i + k, j + 1) is ak(n − j) – Rate from (i, j) to (i + k, j) is bk(n − j) + ckj – Rate from (i, j) to (i + k, j − 1) with rate dkj

• So A0 is of tri-diagonal form

A0 =               b1n a1n d1 b1(n − 1) + c1 a1(n − 1) d12 b1(n − 2) + c12 ... ... ... b12 + c1(n − 2) a12 d1(n − 1) b1 + c1(n − 1) a1 d1n c1n              

54

Quasi-Birth-and-Death processes: Linear rates

• Transition rate diagram

i, j

i + 1, j i, j + 1 i, j − 1

a1(n − j) a−1(n − j) b1j + c1(n − j) d0j

55

Quasi-Birth-and-Death processes: Example with linear rates

• n parallel servers serving two customer types 1 and 2
• Type 1 customers require exponential service with rate µ1, type 2 with rate µ2
• Customers arrive as Poisson stream with rate λ, fraction p1 is of type 1 and p2 of type 2
• States (i, j) with i number of customers in the queue and j number of type 1 customers in service
• Transition rate diagram

i, j

i + 1, j i, j + 1 i, j − 1

jµ1p2 (n − j)µ2p1 jµ1p1 + (n − j)µ2p2

i − 1, j

λ

• Non-zero coefficients a−1 = µ2p1, b1 = c1 = λ/n, b−1 = µ2p2, c−1 = µ1p1, d−1 = µ1p2

56

Quasi-Birth-and-Death processes: Example with linear rates

• n parallel servers and Poisson arrivals with rate λ which need exponential service times with rate µ
• Server switched off when queue is empty
• Idle server switched on when customer arrives which requires exponential setup time with rate θ
• States (i, j) with i number of customers in the queue and j number of active servers
• Transition rate diagram

i, j

i + 1, j i, j + 1 i, j − 1 i − 1, j

λ (n − j)θ jµ

• Non-zero coefficients a0 = θ, b1 = c1 = λ/n, c−1 = µ

57

Quasi-Birth-and-Death processes: Example with linear rates

• Solution: Spectral expansion

– Seek solutions of balance equations of the form pi = y · xi, i = 0, 1, 2, ... , where y = (y(0), ... , y(n)) is non-null vector and |x| < 1 – Substitution and dividing by common powers of x yields following equations for j = 0, · · · , n 0 = (j + 1)D(x)y(j + 1) + ((n − j)B(x) + jC(x))y(j) −((n − j)(A(1) + B(1)) + j(C(1) + D(1)))xy(j) + (n − j + 1)A(x)y(j − 1) where A(x) = a1+a0x +a−1x2, B(x) = b1+b0x +b−1x2, C(x) = c1+c0x +c−1x2, D(x) = d1+d0x +d−1x2 – Solve equations for y by generating function Y (z) =

n

• i=0

y(j)zj

58

Quasi-Birth-and-Death processes: Example with linear rates

• Solution: Spectral expansion

– Multiply equations by zj and add all equations D(x)Y ′(z) + nB(x)Y (z) + (C(x) − B(x))zY ′(z) − (A(1) + B(1))nxY (z) +(A(1) + B(1) − C(1) − D(1))xzY ′(z) + nA(x)zY (z) − A(x)z2Y ′(z) = 0 – This equation can be rewritten as Y ′(z) Y (z) = n[A(x)z + B(x) − (A(1) + B(1))x] A(x)z2 − ((A(1) + B(1) − C(1) − D(1))x + C(x) − B(x))z − D(x) = E(x) z − z1(x) + n − E(x) z − z2(x) where z1(x) and z2(x) are the roots of the denominator and E(x) satisfies 2E(x) − n = n · B(x) + C(x) − x (A(1) + B(1) + C(1) + D(1))

• F(x)2 + 4A(x)D(x)

, where F(x) = (A(1) + B(1) − C(1) − D(1))x + C(x) − B(x)

59

Quasi-Birth-and-Death processes: Example with linear rates

• Solution: Spectral expansion

– Differential equation Y ′(z) Y (z) = E(x) z − z1(x) + n − E(x) z − z2(x) – General solution Y (z) = C(z − z1(x))E(x)(z − z2(x))n−E(x) – Key observation: Y (z) is a polynomial in z, so exponent E(x) = k for k = 0, ... , n – For each k = 0, ... , n the equation 2k − n = n · B(x) + C(x) − x (A(1) + B(1) + C(1) + D(1))

• F(x)2 + 4A(x)D(x)

, has a unique solution x = xk in the interval (0, 1)

60

Quasi-Birth-and-Death processes: Example with linear rates

• Solution: Spectral expansion

– For each k = 0, ... , n the equation 2k − n = n · B(x) + C(x) − x (A(1) + B(1) + C(1) + D(1))

• F(x)2 + 4A(x)D(x)

, has a unique solution x = xk in the interval (0, 1)

−6 −4 −2 2 4 6 1 x0 x1 x2 x3 x4

61

Two-dimensional Markov chain: Re-entrant line (1)

∞ µ1 µ2 µ3

• Two stations, station 1 with two servers, labeled 1 and 3 and station 2 with single server, labeled server 2
• Queue i is served at exponential rate µi
• Queue 1 has infinite supply and customers subsequently go through queue 1, 2 and 3

62

Two-dimensional Markov chain: Re-entrant line (1)

µ3 µ3 µ1 µ1 µ1 µ2 µ2 i j

• Markov chain with states (i, j) where i number in queue 2 and j number in queue 3
• Global balance equations

p(i, j)(µ1 + µ2 + µ3) = p(i − 1, j)µ1 + p(i, j + 1)µ3 + p(i + 1, j − 1)µ2, i > 0, j > 0 p(i, 0)(µ1 + µ2) = p(i − 1, 0)µ1 + p(i, 1)µ3, i > 0 p(0, j)(µ1 + µ3) = p(0, j + 1)µ3 + p(1, j − 1)µ2, j > 0 p(0, 0)µ1 = p0,1µ3

63

Two-dimensional Markov chain: Re-entrant line (1)

µ3 µ3 µ1 µ1 µ1 µ2 µ2 i j

• Substitution of trial solution p(i, j) = αiβj in global balance equations

αβ(µ1 + µ2 + µ3) = βµ1 + αβ2µ3 + α2µ2 α(µ1 + µ2) = µ1 + αβµ3 β(µ1 + µ3) = β2µ3 + αµ2

• Equations solved by α = µ1

µ2 and β = µ1 µ3

• Normalization yields

p(i, j) = (1 − α)αi(1 − β)βj

64

Two-dimensional Markov chain: Re-entrant line (1)

• Global balance equations

p(i, j)(µ1 + µ2 + µ3) = p(i − 1, j)µ1 + p(i, j + 1)µ3 + p(i + 1, j − 1)µ2, i > 0, j > 0

• Partial balance

– Rate out of (i, j) due to arrival in queue 2 is equal to rate into (i, j) due to departure from queue 3 p(i, j)µ1 = p(i + 1, j)µ3 – Rate out of (i, j) due to departure from queue 2 is equal to rate into (i, j) due to arrival in queue 2 p(i, j)µ2 = p(i − 1, j)µ1 – Rate out of (i, j) due to departure from queue 3 is equal to rate into (i, j) due to arrival in queue p(i, j)µ3 = p(i + 1, j − 1)µ2

• Partial balance yields candidate solution

p(i, j) = µ2 µ3 j p(i + j, 0) = µ2 µ3 j µ1 µ2 i+j p(0, 0) = µ1 µ2 i µ1 µ3 j p(0, 0)

65

Two-dimensional Markov chain: Re-entrant line (2)

∞ µ1 µ2 µ3

• Two stations, both stations with single server
• Queue i is served at exponential rate µi
• Queue 1 has infinite supply and customers subsequently go through queue 1, 2 and 3
• Customers in queue 3 have preemptive priority over customers in queue 1 (Last Buffer First Served)
• Assumption 1

µ1 + 1 µ3 > 1 µ2

66

Two-dimensional Markov chain: Re-entrant line (2)

µ3 µ3 µ1 µ2 µ2 i j

• Markov chain with states (i, j) where i number in queue 2 and j number in queue 3
• Global balance equations

p(i, j)(µ2 + µ3) = p(i, j + 1)µ3 + p(i + 1, j − 1)µ2, i > 0, j > 0 p(i, 0)(µ1 + µ2) = p(i − 1, 0)µ1 + p(i, 1)µ3, i > 0 p(0, j)µ3 = p(0, j + 1)µ3 + p(1, j − 1)µ2, j > 0 p(0, 0)µ1 = p(0, 1)µ3

67

Two-dimensional Markov chain: Re-entrant line (2)

µ3 µ3 µ1 µ2 µ2 i j

• Substitution of trial solution p(i, j) = αiβj in global balance equations for i > 0

β(µ2 + µ3) = β2µ3 + αµ2 α(µ1 + µ2) = µ1 + αβµ3

• Equations solved by α = µ1

µ2(1 − β) and β is the root on (0, 1) of β2µ3 − β(µ1 + µ2 + µ3) + µ1 = 0

β = µ1 + µ2 + µ3 −

• (µ1 + µ2 + µ3)2 − 4µ1µ3

2µ3

68

Two-dimensional Markov chain: Re-entrant line (2)

µ3 µ3 µ1 µ2 µ2 i j

• Substitution of p(i, j) = αiβj for i > 0 in global balance equations for i = 0

p(0, j)µ3 = p(0, j + 1)µ3 + αβj−1µ2, j > 0

• Equations are solved by p(0, j) = cβj−1 where c = µ1

µ3

• Normalization yields

p(i, j) =     

µ3 µ1+µ3(1 − α),

i = 0, j = 0

µ1 µ1+µ3(1 − α)βj−1, i = 0, j > 0 µ3 µ1+µ3(1 − α)αiβj, i > 0, j ≥ 0

69

Two-dimensional Markov chain: Single server priority queue

• Single server serving type 1 and type 2 customers
• Type i customers arrive according to Poisson stream with rate λi
• Type i customers require exponential service with rate µi
• Type 1 customers have preemptive priority over type 2 customers

70

Two-dimensional Markov chain: Single server priority queue

i j λ2 λ2 λ2 λ1 λ1 λ1 µ1 µ1 µ2

• Markov chain with states (i, j) where i number of type 1 and j number of type 2
• Global balance equations

p(i, j)(λ1 + λ2 + µ1) = p(i − 1, j)λ1 + p(i + 1, j)µ1 + p(i, j − 1)λ2, i > 0, j > 0 p(i, 0)(λ1 + λ2 + µ1) = p(i − 1, 0)λ1 + p(i + 1, 0)µ1, i > 0 p(0, j)(λ1 + λ2 + µ2) = p(1, j)µ1 + p(0, j − 1)λ2 + p(0, j + 1)µ2, j > 0 p(0, 0)(λ1 + λ2) = p(0, 1)µ2 + p(1, 0)µ1

• Note that p(0, 0) = 1 − ρ1 − ρ2 where ρi = λi

µi

71

Two-dimensional Markov chain: Single server priority queue

• Solution: Difference equation

– Balance equations for j = 0 are homogeneous second-order difference equation p(i, 0)(λ1 + λ2 + µ1) = p(i − 1, 0)λ1 + p(i + 1, 0)µ1, i > 0 – General solution p(i, 0) = c0,0xi

0 + c0,1xi 1,

i ≥ 0 where 0 < x0 < 1 < x1 are roots of quadratic equation x(λ1 + λ2 + µ1) = λ1 + x2µ1 – Convergence of p(i, 0) implies c0,1 = 0 so p(i, 0) = c0,0xi = (1 − ρ1 − ρ2)xi

0,

i ≥ 0

72

Two-dimensional Markov chain: Single server priority queue

• Solution: Difference equation

– Balance equations for j = 1 are inhomogeneous second-order difference equation p(i, 1)(λ1 + λ2 + µ1) = p(i − 1, 1)λ1 + p(i + 1, 1)µ1 + p(i, 0)λ2 = p(i − 1, 1)λ1 + p(i + 1, 1)µ1 + c0,0xi

0λ2,

i > 0 – General solution p(i, 1) = c1,0xi

0 + c1,1(i + 1)xi

= c1,0 i + 1

• xi

0 + c1,1

i + 1 1

• xi

0,

i ≥ 0 where c1,1 = λ2c0,0 λ1 + λ2 − 2µ1x0 – Coefficient c1,0 follows from balance equation in (0, 0) c1,0 = λ1 + λ2 − µ1x0 µ2 c0,0 − c1,1

73

Two-dimensional Markov chain: Single server priority queue

• Solution: Difference equation

– Repeating procedure for balance equations for j = 2, 3, ... leads to p(i, j) =

j

• k=0

cj,k i + j k

• xi

0,

i, j ≥ 0 where coefficients cj,k can be determined recursively starting with c0,0 = 1 − ρ1 − ρ2

74

Two-dimensional Markov chain: Single server priority queue

• Solution: Generating function

– Let P(x, y) = ∞

i=0

∞

j=0 p(i, j)xiyj for |x|, |y| ≤ 1

– Multiply balance equations by xiyj and add all equations h1(x, y)P(x, y) = h2(x, y)P(0, y) + h3(x, y)P(0, 0) with h1(x, y) = λ1xy(1 − x) + λ2xy(1 − y) − µ1y(1 − x) h2(x, y) = −µ1y(1 − x) + µ2x(1 − y) h3(x, y) = −µ2x(1 − y) and P(0, 0) = p(0, 0) = 1 − ρ1 − ρ2

75

Two-dimensional Markov chain: Single server priority queue

• Solution: Generating function

– For every |y| ≤ 1 let x = ξ(y) be the unique root with |x| ≤ 1 of the quadratic equation 0 = h1(x, y) = λ1xy(1 − x) + λ2xy(1 − y) − µ1y(1 − x) – Then x = ξ(y) = λ1 + µ1 + λ2(1 − y) −

• (λ1 + µ1 + λ2(1 − y))2 − 4λ1µ1

2λ1 – Substituting x = ξ(y) and P(0, 0) = 1 − ρ1 − ρ2 in h1(x, y)P(x, y) = h2(x, y)P(0, y) + h3(x, y)P(0, 0) gives P(0, y) = −h3(ξ(y), y)(1 − ρ1 − ρ2) h2(ξ(y), y)

76

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix geometric

– Partition states space into levels i = {(i, 0), (i, 1), ...} – Generator Q of the form Q =       B1 A0 ... A2 A1 A0 A2 A1 A0 A2 A1 . . . ...       where (λ = λ1 + λ2) A0 =   λ1 λ1 ...   , A1 =   −λ − µ1 λ2 −λ − µ1 λ2 ... ...   , A2 =   µ1 µ1 ...   B1 =   −λ λ2 µ2 −λ − µ2 λ2 ... ... ...  

77

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix geometric

– A0, A1 and A2 are upper triangular and have repeating structure – Implication: R is also upper triangular and of the form R =     r0 r1 r2 ... r0 r1 ... r0 ... ...     – Elements can be determined recursively from matrix equation A0 + RA1 + R2A2 = 0 µ1r2

0 − (λ + µ1)r0 + λ1 = 0

µ1

k

• l=0

rk−lrl − (λ + µ1)rk + λ2rk−1 = 0, k ≥ 1

78

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix geometric

– Vector p0 = (p(0, 0), p(0, 1), ...) follows from boundary equations 0 = p0B1 + p1A2 = p0(B1 + RA2) yielding the recursion p(0, 1)µ2 = p(0, 0)(λ − µ1r0) p(0, j + 1)µ2 = (λ + µ2)p(0, j) − λ2p(0, j − 1) − µ1

j

• k=0

p(0, j − k)rk, j ≥ 1 starting with p(0, 0) = 1 − ρ1 − ρ2 – Vector pi+1 = (p(0, 0), p(0, 1), ...) can be obtained from pi+1 = piR yielding p(i + 1, j) =

j

• k=0

p(i, j − k)rk

79

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix analytic

– A0, A1 and A2 are upper triangular and have repeating structure – Implication: G is also upper triangular and of the form G =     g0 g1 g2 ... g0 g1 ... g0 ... ...     – gk is probability of first passage to level i − 1 in (i − 1, j + k) for initial state (i, j) at level i – One-step analysis g0 = µ1 λ + µ1 + λ1 λ + µ1 g2 gk = λ2 λ + µ1 gk−1 + λ1 λ + µ1

k

• l=0

glgk−l, k ≥ 1

80

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix analytic

– From pi+1A2 = piA0G follows p(i + 1, j)µ1 =

j

• k=0

p(i, j − k)λ1gk, i, j, ≥ 0

81

Two-dimensional Markov chain: Single server priority queue

• Solution: Matrix analytic

– Vector p0 follows from boundary equations 0 = p0B1 + p1A2 = p0(B1 + A0G) – B1 + A0G is generator of Markov chain embedded on level 0

λ2 λ1gk µ2 (0, j) (0, j + k)

– Global balance for {(0, 0), (0, 1), ... , (0, j)} p(0, j + 1)µ2 = p(0, j)λ2 +

j

• k=0

p(0, j − k)λ1

• 1 −

k

• l=0

gl

82

Two-dimensional Markov chain: Shortest queue

µ µ λ

• Two parallel exponential servers with rate µ, each with its own queue
• Customers arrive according to Poisson stream with rate λ and join shortest queue on arrival

83

Two-dimensional Markov chain: Shortest queue

λ µ µ µ λ λ 2µ

• Markov chain with states (i, j) where i number in shortest queue and j difference between queues
• Global balance equations

p(i, j)(λ + 2µ) = p(i − 1, j + 1)λ + p(i, j + 1)µ + p(i + 1, j − 1)µ, i > 0, j > 1 p(i, 1)(λ + 2µ) = p(i − 1, 2)λ + p(i, 2)µ + p(i, 0)λ + p(i + 1, 0)2µ, i > 0 p(i, 0)(λ + 2µ) = p(i − 1, 1)λ + p(i, 1)µ, i > 0 p(0, j)(λ + µ) = p(0, j + 1)µ + p(1, j − 1)µ, j > 1

84

Two-dimensional Markov chain: Shortest queue

λ µ µ µ λ λ 2µ

• Markov chain with states (i, j) where i number in shortest queue and j difference between queues
• Global balance equations

p(i, j)(λ + 2µ) = p(i − 1, j + 1)λ + p(i, j + 1)µ + p(i + 1, j − 1)µ, i > 0, j > 1 p(i, 1)(λ + 2µ) = p(i − 1, 2)λ + p(i, 2)µ +(p(i − 1, 1)λ + p(i, 1)µ) λ λ + 2µ + (p(i, 1)λ + p(i + 1, 1)µ) 2µ λ + 2µ, i > 0 p(0, j)(λ + µ) = p(0, j + 1)µ + p(1, j − 1)µ, j > 1

85

Intermezzo: Method of images

p ∆φ=0 φ=0 q

• What is the potential Φ in point p?

86

Intermezzo: Method of images

−q ∆φ=0 φ=0 q p

• What is the potential Φ in point p?
• Add image charge −q:
• Potential Φ in point p is

Φ(p) = q

r1 − q r2 where r1 and r2 are distances from point p to charge q and −q

87

Intermezzo: Method of images

p ∆φ=0 φ=0 q −q −q q

• For two planes at right angles

Φ(p) = q

r1 − q r2 + q r3 − q r4

88

Intermezzo: Method of images

a a c b A B Φ=Φ Φ=0

• What is the potential Φ in any point p outside the spheres A and B?

89

Intermezzo: Method of images

a c b A B Φ=Φ Φ=0 a α

• What is the potential Φ in any point p outside the spheres A and B?
• Add image charge α0 in center of A with α0 = aΦa: Then

Φ(p) ≈ α0

r0 where r0 is distance from p to charge α0 (correct if c = ∞)

• Potential of α0 does not vanish on B

90

Intermezzo: Method of images

β a c b A B Φ=Φ Φ=0 a α0

• Place new image charge β0 in B to compensate for this error

Φ(p) ≈ α0

r0 + β0 r1 where r0 and r1 are distances from p to charges α0 and β0

• Charge β0 alters potential on A

91

Intermezzo: Method of images

1 a c b A B Φ=Φ Φ=0 a α0 β α

• Place new image charge β0 in B to compensate for this error

Φ(p) ≈ α0

r0 + β0 r1 where r0 and r1 are distances from p to charges α0 and β0

• Charge β0 alters potential on A
• Place new image charge α1 in A to compensate for this error

Φ(p) ≈ α0

r0 + β0 r1 + α1 r2

92

Intermezzo: Method of images

1 a c b A B Φ=Φ Φ=0 a α0 β α

• Potential problem of two spheres is solved by infinite sequence of image charges in A and B

Φ(p) = α0

r0 + β0 r1 + α1 r2 + β1 r3 + · · ·

93

Two-dimensional Markov chain: Shortest queue

λ µ µ µ λ λ 2µ

• Image charges are products αiβj satisfying the global balance equations

p(i, j)(λ + 2µ) = p(i − 1, j + 1)λ + p(i, j + 1)µ + p(i + 1, j − 1)µ, i > 0, j > 1

• Substituting αiβj in this equation and dividing by common powers gives the curve

αβ(λ + 2µ) = β2λ + αβ2µ + α2µ

94

Two-dimensional Markov chain: Shortest queue

• Quadratic curve

αβ(λ + 2µ) = β2λ + αβ2µ + α2µ

95

Two-dimensional Markov chain: Shortest queue

• Unique product form satisfying global balance equations for i > 0, j > 1 and i > 0, j = 1:

p(i, j) ≈ c0αi

0βj

where α0 = ρ2, β0 = ρ2/(2 + ρ) with ρ = λ/2µ

• Product form c0αi

0βj 0 violates global balance equations for i = 0

96

Two-dimensional Markov chain: Shortest queue

• Add new product form c1αi

1βj 1 such that sum satisfies balance equations for i = 0

p(i, j) ≈ c0αi

0βj 0 + c1αi 1βj 1

• Then β1 = β0 and α1 is companion root of α0 and

c1 = −α1 − β0 α0 − β0 c0

• Product form c1αi

1βj 1 violates global balance equations for j = 1

97

Two-dimensional Markov chain: Shortest queue

• Add new product form c2αi

2βj 2 such that sum satisfies balance equations for j = 1

p(i, j) ≈ c0αi

0βj 0 + c1αi 1βj 1 + c2αi 2βj 2

• Then α2 = α1 and β2 is companion root of β1 and (with ρ = λ/(2µ))

c2 = −(ρ + α1)/β2 − (1 + ρ) (ρ + α1)/β1 − (1 + ρ)c1

• Product form c2αi

2βj 2 violates global balance equations for i = 0

98

Two-dimensional Markov chain: Shortest queue

• Iterating yields infinite sequence

p(i, j) = c0αi

0βj 0 + c1αi 1βj 1 + c2αi 2βj 2 + c3αi 3βj 3 + · · ·

99

Two-dimensional Markov chain: Random walk in quarter plane

i,j

q

• Equilibrium probabilities can be expressed as infinite sequence of product forms

p(i, j) = c0αi

0βj 0 + c1αi 1βj 1 + c2αi 2βj 2 + · · ·

100

Two-dimensional Markov chain: Random walk in quarter plane

i,j

q

• Equilibrium probabilities can be expressed as infinite sequence of product forms

p(i, j) = c0αi

0βj 0 + c1αi 1βj 1 + c2αi 2βj 2 + · · ·

• Provided there are no transitions to North, North-East and East

q0,1 = q1,1 = q1,0 = 0

101

Two-dimensional Markov chain: Random walk in quarter plane

• Condition q0,1 = q1,1 = q1,0 = 0 implies that curve for α and β is of the above form
• This form guarantees that αk → 0 and βk → 0 as k → ∞

102

Multi-skilled system

• Customer types C = {a, b, ...}
• Independent Poisson arrivals with rates λc, c ∈ C
• Exponential service requirements with mean 1
• J servers, S = {m1, ... , mJ}
• Server mj works at rate µmj
• Skill-based service:

– Server mi can serve customer types C(mj) – Type c customer can be served by servers S(c)

• Service discipline combination of FCFS and ALIS:

– Server picks longest waiting compatible customer – Customer assigned to longest idle compatible server

103

Multi-skilled system: Example

µm2 µm1 λa,λb,λc a µm3 a,c m1 m2 m3 a b c a,b m1 m2 m3

• C = {a, b, c} and S = {m1, m2, m3}
• Skill-based service:

– C(m1) = {a, b}, C(m2) = {a, c}, C(m3) = {a} – S(a) = {m1, m2, m3}, S(b) = {m1}, S(c) = {m2}

104

Multi-skilled system: Markov chain

• Markov chain with states (M1, n1, M2, ... , Mi, ni, Mi+1, ... , MJ)

– i busy servers M1, ... , Mi – J − i idle servers Mi+1, ... , MJ with increasing idle times – M1, ... , MJ is a permutation of m1, ... , mJ – nj number of waiting jobs between server Mj and Mj+1

105

Multi-skilled system: Examples of states

m3

* *

a/b/c

m1

a/b/c b/c b/c b/c b/c b b b

m2 m3

a/b a/c a

m1

b b b a/b

m2

*

m1

b b b a/b

m3

a

(i) All servers busy, servers move left to right, arrivals come from the right (ii) Two servers busy, no possible customers between m3 and m

1

(iii) Two servers are idle,m3 has been idle longest time

m2

(i) (m1, 3, m2, 4, m3, 2) (ii) (m3, 0, m1, 3, m2) (iii) (m1, 3, m2, m3)

106

Multi-skilled system: N-system

b a m1 m2 λa, λb m1 m2 µ1 µ2 a, b a

• Markov chain with states

– (m1, n1, m2, n2) with n1, n2 ≥ 0 – (m2, 0, m1, n2) with n2 ≥ 0 – (m1, n1, m2) with n1 ≥ 0 – (m2, 0, m1) – (m1, m2) and (m2, m1) (empty system)

• Assume λb < µ1 and λa + λb < µ1 + µ2

107

Multi-skilled system: N-system

λa + λb µ2p1 µ2p2 µ1 n1 (m1, n1, m2, n2) n2

• State transition from (m1, n1, m2, n2) to (m1, n1 + j − 1, m2, n2 − j) with rate µ2pj where pj = (1 − γ )γ j−1

γ = λb λa + λb

108

Multi-skilled system: N-system

λa + λb µ2p1 µ2p2 µ1 n1 (m1, n1, m2, n2) n2

• Global balance equation in (m1, n1, m2, n2) with n2 > 0

(λa + λb + µ1 + µ2)p(m1, n1, m2, n2) = (λa + λb)p(m1, n1, m2, n2 − 1) + µ1p(m1, n1 + 1, m2, n2) +µ2(1 − γ )

n1

• j=0

γ jp(m1, n1 − j, m2, n2 + j + 1) +µ2(1 − γ )γ n1p(m2, 0, m1, n1 + 1 + n2)

109

Multi-skilled system: Partial balance

(i) Rate out of state (M1, n1, M2, ... , Mi, ni, Mi+1, ... , MJ) due to an arrival activating a server equals rate into that state due to a departure deactivating Mi+1 (ii) Rate out of state (M1, n1, M2, ... , Mi, ni, Mi+1, ... , MJ) due to an arrival that is queueing equals rate into that state due to a departure from a server staying active (iii) Rate out of state (M1, n1, M2, ... , Mi, 0, Mi+1, ... , MJ) with ni = 0 due to a departure equals rate into that state due to an arrival activating Mi (iv) Rate out of state (M1, n1, M2, ... , Mi, ni, Mi+1, ... , MJ) with ni > 0 due to a departure equals rate into that state due to an arrival that is queueing

110

Multi-skilled system: N-system

(iv) Rate out of state (m1, n1, m2, n2) with ni > 0 due to departure equals rate into that state due to an arrival that is queueing (λa + λb + µ1 + µ2)p(m1, n1, m2, n2) = (λa + λb)p(m1, n1, m2, n2 − 1) + µ1p(m1, n1 + 1, m2, n2) +µ2(1 − γ )

n1

• j=0

γ jp(m1, n1 − j, m2, n2 + j + 1) +µ2(1 − γ )γ n1p(m2, 0, m1, n1 + 1 + n2)

111

Multi-skilled system: N-system

(ii) Rate out of state (m1, n1, m2, n2) due to an arrival that is queueing equals rate into that state due to a departure from a server staying active (λa + λb + µ1 + µ2)p(m1, n1, m2, n2) = (λa + λb)p(m1, n1, m2, n2 − 1) + µ1p(m1, n1 + 1, m2, n2) +µ2(1 − γ )

n1

• j=0

γ jp(m1, n1 − j, m2, n2 + j + 1) +µ2(1 − γ )γ n1p(m2, 0, m1, n1 + 1 + n2)

112

Multi-skilled system: N-system

• Partial balance dictates solution

p(m1, n1, m2, n2)

(iv)

= p(m1, n1, m2, n2 − 1) λa + λb µ1 + µ2

(iv)

= p(m1, n1, m2, 0) λa + λb µ1 + µ2 n2

(iii)

= p(m1, n1, m2) λa µ1 + µ2 λa + λb µ1 + µ2 n2

(iv)

= p(m1, 0, m2) λb µ1 n1 λa µ1 + µ2 λa + λb µ1 + µ2 n2

(i)

= p(m1, m2)λa + λb µ1 λb µ1 n1 λa µ1 + µ2 λa + λb µ1 + µ2 n2

113

Multi-skilled system: N-system

• Candidate solution

p(m1, n1, m2, n2) = p(m1, m2)λa + λb µ1 λb µ1 n1 λa µ1 + µ2 λa + λb µ1 + µ2 n2 p(m2, 0, m1, n2) = p(m1, m2)λa µ2 λa + λb µ1 + µ2 n2+1 p(m1, n1, m2) = p(m1, m2)λa + λb µ1 λb µ1 n1 p(m2, 0, m1) = p(m1, m2)λa µ2 p(m2, m1) = p(m1, m2) λa λa + λb

114

Multi-skilled system: N-system

• Set C = p(m1, m2)(λa + λb)λa then candidate solution

p(m1, n1, m2, n2) = C λn1

b

µn1+1

1

(λa + λb)n2 (µ1 + µ2)n2+1 p(m2, 0, m1, n2) = C 1 µ2 (λa + λb)n2 (µ1 + µ2)n2+1 p(m1, n1, m2) = C λn1

b

µn1+1

1

1 λa p(m2, 0, m1) = C 1 µ2 1 λa + λb p(m2, m1) = C 1 λa + λb 1 λa + λb p(m1, m2) = C 1 λa + λb 1 λa

• Checking partial balance equations proves that this is solution to the global balance equations

115

Multi-skilled system

• C is set of customer types {a, b, ...}
• S is set of servers {m1, m1, ... , mJ}
• λX =

c∈X λc where X ⊂ C

• µY =

M∈Y µM where Y ⊂ S

• S(X) is set of server types that can handle customers in X
• C(Y) is set of customer types that can be handled by servers in Y
• U(Y) is set of customer types that can be uniquely handled by servers in Y
• Assume (stability) λU(Y) < µY for every subset of servers Y ⊂ S

116

Multi-skilled system

• Solution to the global balance equations

p(M1, n1, ... , Mi, ni, Mi+1, ... , MJ) = C

i

• j=1

λnj

U({M1,...,Mj})

µnj+1

{M1,...,Mj} J

• j=i+1

λ−1

C({Mj,...,MJ})

where C is a normalizing constant

• Proof by checking partial balance equations

117

Multi-skilled system: Alternative Markov chain

• Markov chain with states (C1, ... , CL, M1, ... , MK)
• L customers waiting for service C1, ... , CL with decreasing waiting times (L ≥ 0)
• K idle servers M1, ... , MK with decreasing idle times (K ≤ J)
• Note that L = n1 + ... ni and M1, ... , MK correspond to MJ, ... , Mi+1
• Solution to the global balance equations

p(C1, ... , CL, M1, ... , MK) = C

L

• l=1

λCl µS({C1,...,Cl})

K

• k=1

µMk λC({M1,...,Mk}) where C is a normalizing constant

• Proof by checking partial balance equations
• Note that

p(C1, ... , CL|all servers busy) ∝

L

• l=1

λCl µS({C1,...,Cl})

118

FCFS matching queue

• Independent Poisson arrivals of customers with rates λc, c ∈ C
• Independent Poisson arrivals of servers with rates µmj, mj ∈ S
• Arriving server mj scans queue of customers

– Matches with the longest waiting compatible customer and the two leave the system immediately – If no match is possible he leaves immediately without customer

• Note that this system is equivalent with a redundancy service system

119

FCFS matching queue: Markov chain

• Markov chain with states (C1, ... , CL)
• L customers waiting for service C1, ... , CL with decreasing waiting times (L ≥ 0)
• Solution to the global balance equations

p(C1, ... , CL) = C

L

• l=1

λCl µS({C1,...,Cl}) where C is a normalizing constant