Fluid Models of Parallel Service Systems under FCFS Hanqin Zhang - - PowerPoint PPT Presentation

Fluid Models of Parallel Service Systems under FCFS

Hanqin Zhang

Business School, National University of Singapore Joint work with Yuval Nov and Gideon Weiss from The University of Haifa, Israel Queueing and Networks, IMA, 2018

1 / 39

Parallel Service System

◮ Parallel servers S = {s1, . . . , sJ} of various skills; ◮ Various customer types C = {c1, . . . , cI}; ◮ A bipartite compatibility graph G where (sj, ci) ∈ G if

server sj can serve customers of type ci

λ c4 c3 c2 c1 s1 s2 s3

µ1,1 µ1,2 µ2,2 µ3,1 µ1,4 µ2,3 µ3,3 µ3,4 α1 α2 α3 α4

Figure: A parallel skilled based service system with 3-server and 4-customer-type

2 / 39

λ c4 c3 c2 c1 s1 s2 s3

µ1,1 µ1,2 µ2,2 µ3,1 µ1,4 µ2,3 µ3,3 µ3,4 α1 α2 α3 α4

Figure: A parallel skilled based service system with 3-server and 4-customer-type

• C(sj): the customer types compatible with sj;

For example, C(s3) = {c1, c3, c4}.

• S(ci): the servers compatible with customers of type ci.

For example, S(c2) = {s1, s2}. S(C) =

ci∈C S(ci), C(S) = sj∈S C(sj), U(S) = C(S),

αC =

ci∈C αci.

3 / 39

◮ a(ℓ) denotes the arrival time of the ℓth customer; ◮ u(ℓ) = a(ℓ) − a(ℓ − 1) is the interarrival times, where

ℓ = 0, ±1, ±2, . . ., and a(0) ≤ 0 < a(1);

◮ {u(ℓ) : ℓ = 0, ±1, ±2, . . .} is a sequence of iid with mean

1/λ;

◮ A(t) = max{ℓ : a(ℓ) ≤ t}; ◮ −A(t) is the number of customers arriving in (t, 0] for t < 0; ◮ A(t) gives the number of customers arriving in (0, t] for

t > 0;

◮ A(t) − A(s) counts the total number of arrivals in (s, t]

with s < t.

4 / 39

Type ci has probability αci, i = 1, . . . , I, let ξ(ℓ) be a unit vector of length I such that ξi(ℓ) = 1 if customer ℓ is of type ci, for ℓ = 0, ±1, ±2, . . .. The counts of arrivals of customers of each type are then given by Aci(t) =             

A(t)

• ℓ=1

ξi(ℓ), t ≥ 0, −

• ℓ=A(t)+1

ξi(ℓ), t < 0.

5 / 39

◮ vsj,ci(0) be the remaining service time of server sj if he is

serving a customer of type ci at time 0, and vsj,ci(0) = 0

• therwise;

◮ vsj,ci(k), k = 1, 2, . . . , be the processing time of the kth

customer of type ci that server sj is serving after time 0 with mean 1/µsj,ci;

◮ Xsj,ci(t) = max{k + 1 : k ℓ=0 vsj,ci(ℓ) ≤ t}.

6 / 39

1 1 1 3 1 2 2 2 1 1 3 3 3 1 3 3 2 1 2 3 1 1 1 2 2 3 3 2 2

time customers a(0) a(1) a(2) a(3) a(4) a(5) a(6) a(7) a(8) a(9) a(10) a(-1) a(-2) a(-3) a(-4) a(-5) a(-6)

t

Customer in service Arrival Service Start Departure Server Position at 0 Server Position at t Customer type Arrival time a( ) Customer in queue Key

S C

Bipartite Graph

Figure: Dynamics of a 3-server and 3-customer-type system under FCFS-ALIS

7 / 39

◮ Psj(t) is the position of server sj at time t, where we let

Psj(t) = ℓ if the server is serving at time t the ℓth customer in the sequence of arrivals. If servers sj1, . . . , sjk are idle at time t then their positions are defined as A(t) + 1, . . ., A(t) + k, ordered by duration of idleness, with A(t) + k the longest idle;

◮ Yj(t) is the current jth level, where we let

Y1(t) < . . . < YJ(t) be the ordered set of the positions of J servers at time t;

◮ Tsj,ci(t) is the cumulative time over (0, t) that server sj has

served customers of type ci. Let M1(t) = arg min{Psj(t) : sj ∈ S}, Mj(t) = arg min {Psℓ(t) : sℓ ∈ S \ {M1(t), · · · , Mj−1(t)}} . Clearly, PMj(t)(t) = Yj(t) for j = 1, · · · , J.

8 / 39

We let Qci,j(t) denote the number of customers of type ci which are waiting between servers Mj(t) and Mj+1(t) at time t. These are given by: Qci,j(t) =                   

Yj+1(t)−1

• ℓ=Yj(t)+1

ξi(ℓ) I{ci ∈ U(M1(t), . . . , Mj(t))}, j = 1, . . . , J − 1,

A(t)

• ℓ=YJ(t)+1

ξi(ℓ), j = J. (1) Let Qci(t) be the number of type-ci customers in the system at time t. Then, Qci(t) = Qci(0) + Aci(t) −

• sj∈S(ci)

Xsj,ci(Tsj,ci(t)). (2) Furthermore, the work-conserving principle gives that for t ≥ 0, t

• A(x)+1−Psj(x)

+ d

• x−
• ci∈C(sj)

Tsj,ci(x)

• = 0, sj ∈ S. (3)

9 / 39

◮ U(t): the remaining time at time t until next arrival; ◮ Vsj,ci(t): the remaining processing time of ci by server sj if

it is processing a type ci customer at time t, and Vsj,ci(t) = 0 otherwise;

◮ Rsj,ci(t): the number of type-ci customers served by server

sj by time t;

◮ Z(t) =

• A(t) − Psj(t), Qci,j(t), U(t), Vsj,ci(t)
• are Markov

processes. We say that the parallel service system is stable (ergodic) if {Z(t) : t ≥ 0} is positive Harris recurrent (ergodic).

• Determining condition for the stability of {Z(t) : t ≥ 0};
• Finding condition for

Pr

• limt→∞

Ps1(t) t

= · · · = limt→∞

PsJ (t) t

• = 1;
• Finding condition and determining constant rsj,ci such that

Pr

• limt→∞

Rsj,ci(t) t

= rsj,ci

• = 1 for (sj, ci) ∈ G.

10 / 39

Corresponding to (1)-(3), our fluid model given by

• ¯

Tsj,ci(t), ¯ Psj(t), ¯ Yj(t), ¯ Qci,j(t)

• : t ≥ 0
• is:

¯ P ¯

Mj(t)(t) = ¯

Yj(t), j = 1, . . . , J, (4) ¯ Qci,j(t) =    αci( ¯ Yj+1(t) − ¯ Yj(t)) I{ci ∈ U( ¯ M1(t), . . . , ¯ Mj(t))}, j = 1, . . . , J − 1, αci(λt − ¯ YJ(t)), j = J, (5) ¯ Qci(t) = ¯ Qci(0) + λαcit −

• sj∈S(ci)

µsj,ci ¯ Tsj,ci(t), (6) t

• λx − ¯

Psj(x) + d

• x −
• ci∈C(sj)

¯ Tsj,ci(x)

• = 0,

sj ∈ S. (7) In addition, we assume that ¯ Tsj,ci(x) are zero for (sj, ci) / ∈ G, and nondecreasing with

ci∈C(sj)( ¯

Tsj,ci(t2) − ¯ Tsj,ci(t1)) ≤ t2 − t1. (8)

11 / 39

Theorem 1: For almost all sample paths ω and any subsequence of {( ¯ T n

sj,ci(t), ¯

P n

sj(t), ¯

Y n

j (t), ¯

Qn

ci,j(t)) : n ≥ 1}

contains a subsequence along which {( ¯ T n

sj,ci(t), ¯

P n

sj(t), ¯

Y n

j (t), ¯

Qn

ci,j(t)) : n ≥ 1}

converges u.o.c. in (0, ∞) and any of their limits, ( ¯ Tsj,ci(t), ¯ Psj(t), ¯ Yj(t), ¯ Qci,j(t)), jointly satisfies relations (4)-(8). Moreover, ¯ Psj(t), ¯ Yj(t), and ¯ Qci,j(t) are Lipschitz continuous on (0, ∞), and ¯ Tsj,ci(t) are Lipschitz continuous on [0, ∞) for sj ∈ S, ci ∈ C, j = 1, . . . , J, and (sj, ci) ∈ G.

12 / 39

Definition 1: Consider the fluid model given by (( ¯ Psj(0), λαci, µsj,ci) : si ∈ S, ci ∈ C, (sj, ci) ∈ G). Let | ¯ P(0)| = − J

j=1 ¯

Psj(0). (i) We say that the fluid model is stable if starting from any initial state ¯ P(0) with | ¯ P(0)| = 1, there exists t0 such that for any ( ¯ Y1(t), · · · , ¯ YJ(t)) satisfying (4)-(8), λt − ¯ Y1(t) = 0 for all t > t0. (ii) We say that the fluid model has complete resource pooling if for all values of λ, starting from any initial state ¯ P(0) with | ¯ P(0)| = 1, there exists t0 such that for any ( ¯ Y1(t), · · · , ¯ YJ(t)) satisfying (4)-(8), ¯ YJ(t) − ¯ Y1(t) = 0 for all t > t0. (iii) We say that the fluid model has complete weak resource pooling if for all values of λ, starting from any initial state ¯ P(0) with | ¯ P(0)| = 1, there exists t0 such that for any ( ¯ Y1(t), · · · , ¯ YJ(t)) satisfying (4)-(8), ˙ ¯ Y1(t) = · · · = ˙ ¯ YJ(t) for all t > t0.

13 / 39

Time Customers Server 2 Server 3 Server 1 Arrivals

1 1 2 2 3 3 S C

Bipartite Graph

Figure: Conjectured Fluid Dynamics under FCFS-ALIS

14 / 39

Definition 2: (i) We say that a cumulative probability distribution H(·) of a non-negative random variable has unbounded support if H(t) < 1 for all t. (ii) We say that a cumulative probability distribution H(t) of a non-negative random variable is spread out if there exists a function q(t) > 0 such that ∞

0 q(t)dt > 0 and a constant k,

such that H(k)(b) − H(k)(a) ≥ b

a q(t)dt for all 0 ≤ a < b < ∞,

where H(k) is the k-fold convolution of H(·). Theorem 2 Assume that the distributions of the interarrival times and of the service times have unbounded support and are spread out. (i) If the fluid limit model is stable, then {Z(t) : t ≥ 0} is ergodic. (ii) If the fluid limit model has complete resource pooling, then for λ large enough, the joint distribution of (Qci,j(t), ci ∈ U(M1(t), . . . , Mj(t)), 1 ≤ j ≤ J − 1) converges to a stationary distribution as t → ∞.

15 / 39

Proposition 1: Consider a fluid limit ( ¯ Tsj,ci(t), ¯ Psj(t), ¯ Yj(t), ¯ Qci,j(t)). If fluid server positions at levels k, . . . , ℓ merge at time τ > 0, i.e., ¯ Yk−1(τ) < ¯ Yk(τ) = · · · = ¯ Yℓ(τ) < ¯ Yℓ+1(τ) (or if ℓ = J, ¯ Yℓ(τ) < λτ), for some k ≤ ℓ, then the fluid limit must satisfy

• ci∈C(Mj)\C(Mℓ+1,...,MJ)

˙ ¯ TMj,ci(τ) = 1, j = k, . . . , ℓ, ˙ ¯ Yk(τ) = · · · = ˙ ¯ Yℓ(τ) = 1 αci

ℓ

• j=k

µMj,ci ˙ ¯ TMj,ci(τ) for ci ∈ C(Mk, . . . , Mℓ) \ C(Mℓ+1, . . . , MJ).

16 / 39

Product-Form Service Rates:

λ c4 c3 c2 c1 s1 s2 s3

µ1,1 µ1,2 µ2,2 µ3,1 µ1,4 µ2,3 µ3,3 µ3,4 α1 α2 α3 α4

Assume µsj,ci = µsj × νci, where we can think of 1/νci as the average work required by a customer of type ci and µsj as the speed of server sj. This models the case where customer types that share the same server differ only in the amount of work that they require, and where servers that can serve the same type of customer differ only in the speed at which they perform the service.

17 / 39

λ c4 c3 c2 c1 s1 s2 s3

µ1,1 µ1,2 µ2,2 µ3,1 µ1,4 µ2,3 µ3,3 µ3,4 α1 α2 α3 α4

Service Rates Depend Only on Server: µsj = µsj,ci for all (sj, ci) ∈ G. Service Rates Depend Only on Customer Type: νci = µsj,ci for all (sj, ci) ∈ G.

18 / 39

Condition for complete resource pooling in the product-form service rate (µsj,ci = µsj × νci) case: Introduce the following notation: for any ci, θci = αci

νci ci∈C αci νci ; for any C ⊆ C, θC = ci∈C θci; for any sj,

βsj = µsj

sj∈S µsj; and for any S ⊆ S, βS = sj∈S βsj.

The following three conditions are equivalent (Condition A): (i) βS(C) > θC; (ii) θC(S) > βS; (iii) βS > θU(S).

19 / 39

Theorem 3: Consider a parallel service fluid model with product-form service rate. (i) Assume that Condition A holds, then complete resource pooling holds, that is, for any initial conditions there exists t0 such that for λ large enough, every fluid limit satisfies ¯ Y1(t) = · · · = ¯ YJ(t) and ˙ ¯ Y1(t) = · · · = ˙ ¯ YJ(t) = ¯ µ for all t > t0, where ¯ µ =

• sj∈S

µsj

ci∈C

αci νci . (ii) Assume that Condition A holds only as weak inequalities. Then complete weak resource pooling holds, and furthermore, for any initial conditions there exists t0 such that for λ large enough, every fluid limit satisfies ˙ ¯ Y1(t) = · · · = ˙ ¯ YJ(t) = ¯ µ.

20 / 39

Network with Complete Bipartite Compatibility Graph:

Every server can serve all types of customers, i.e., the compatibility graph is a complete bipartite graph. msj =

• ci∈C

αcimsj,ci, µsj = m−1

sj ,

µ =

J

• j=1

µsj. (9) Theorem 4: Consider the FCFS service discipline fluid model with a complete bipartite compatibility graph. We have (i) It is always complete resource pooling; (ii) The fluid server-level process is given by ¯ Y1(t) = · · · = ¯ YJ(t) = min ¯ YJ(0) + µt, λt

• ,

t > 0, where µ is given in (9); (iii) When ¯ YJ(t) < λt, the matching rates are given rsj,ci = µsjαci/(µ

• ck∈C

αck).

21 / 39

Network with Tree Bipartite Compatibility Graph:

A tree graph is a connected graph with no loops. With K nodes it will have exactly (K − 1) edges, and it will always have at least two leaves

• ci∈C(sk)

˙ ¯ Tsk,ci(τ) = 1, ˙ ¯ Psk(τ) =

• sj∈S(cm)

µsj,cm αcm ˙ ¯ Tsj,cm(τ), ˙ ¯ Ps1(τ) = · · · = ˙ ¯ PsJ(τ). This means that the set of linear equations

ci∈C(sk) ηsk,ci = 1

sk ∈ S,

• sj∈S(cm)

µsj,cm αcm ηsj,cm = µ,

cm ∈ C (10) with the (I + J − 1) unknowns ηsj,ci, (sj, ci) ∈ G and an additional unknown µ has a positive solution. Thus complete resource pooling will relate to the set of linear equations (10).

22 / 39

Theorem 5: Consider a fluid model with tree bipartite compatibility graph. (i) The system will have complete resource pooling if and only if (10) has a positive solution, and it will have complete weak resource pooling if the solution is non-negative. (ii) If complete resource pooling holds then µ given by (10) is the pooled service rate, and the matching rates are given by rsj,ci = µsj,ciηsj,ci/µ.

23 / 39

Foss and Chernova’s Example (stability depends on distributions rather than the first moments):

(On the stability of a partially accessible multistation queue with state-dependent routing. Queueing Systems, 29(1998), 55-73) 3-server, 3-customer-type and an almost complete bipartite compatibility graph as illustrated in the following Figure.

1 2 3 2 3 1 S C λ λ λ FL FR

αc1 = αc2 = αc3 = 1/3, and the service time distributions are vs1,c2 ∼ FL, vs1,c3 ∼ FR; vs2,c1 ∼ FR, vs2,c3 ∼ FL; vs3,c1 ∼ FL, vs3,c2 ∼ FR with means aL and aR respectively. fL(t) = 1 M e−t/M, fR(t) = 1 1 + A 1 Ae−t/A + A 1 + AAe−At

24 / 39

Let L : service time = exp( 1 M ); R1 : service time = exp(A); R2 : service time = exp( 1 A). Then the state-space can be written as for n = 0, 1, · · · , (n, L, L, L), (n, L, L, R1), (n, L, L, R2), (n, L, R1, L), (n, L, R1, R1), (n, L, R1, R2), (n, L, R2, L), (n, L, R2, R1), (n, L, R2, R2); (n, R1, L, L), (n, R1, L, R1), (n, R1, L, R2), (n, R1, R1, L), (n, R1, R1, R1), (n, R1, R1, R2), (n, R1, R2, L), (n, R1, R2, R1), (n, R1, R2, R2); (n, R2, L, L), (n, R2, L, R1), (n, R2, L, R2), (n, R2, R1, L), (n, R2, R1, R1), (n, R2, R1, R2), (n, R2, R2, L), (n, R2, R2, R1), (n, R2, R2, R2).

25 / 39

The generator is given by        D0 D1 D2 D3 · · · C0 C1 C2 C3 · · · C0 C1 C2 · · · C0 C1 · · · . . . . . . . . . . . . ...        where the first block-row corresponds to state n = 0, the second block-row corresponds to state n = 1, and so on.

26 / 39

Dk =   Dk

11

Dk

12

Dk

13

Dk

21

Dk

22

Dk

23

Dk

31

Dk

32

Dk

33

  , k = 0, 1, 2, · · · ; Ck =   Ck

11

Ck

12

Ck

13

Ck

21

Ck

22

Ck

23

Ck

31

Ck

32

Ck

33

  , k = 0, 1, 2, · · · . First consider the transition rate matrix D0. More specifically, D0

11 is the transition rates between (0, L, L, L), (0, L, L, R1),

(0, L, L, R2), (0, L, R1, L), (0, L, R1, R1), (0, L, R1, R2), (0, L, R2, L), (0, L, R2, R1), and (0, L, R2, R2).

27 / 39

D0

11 =

                       δ0

1 1 3a

δ0

2 1 3 1 a

δ0

3 1 3a 1 3 1 M

δ0

4 2 3a

δ0

5 1 3 1 a 1 3a

δ0

6 1 3 1 a 1 3 1 M

δ0

7 1 3 1 a 1 3a

δ0

8 1 3 2 a

δ0

9

                      

28 / 39

Now define G to be the minimal nonnegative solution to

∞

• k=0

CkGk = 0, ¯ Dk =

∞

• i=0

Dk+iGi, ¯ Ck =

∞

• i=0

Ck+iGi, k = 1, 2, · · · . Then πn =

• π0 × ¯

Dn +

n−1

• j=1

πj × ¯ Cn+1−j

• ×
• ¯

C1 −1 , n = 2, 3, · · · , =

• π0, π1
• ·

D0 ∞

j=0 D1+jGj

C0 ∞

j=0 C1+jGj

• .

29 / 39

10 20 30 40 50 0.495 0.505 0.515

M = 1

A alpha hat

Figure: The fraction of customers served with service time distribution FL, as a function of A.

30 / 39

10 20 30 40 50 0.446 0.450 0.454

M = 5

A alpha hat

Figure: The fraction of customers served with service time distribution FL, as a function of A.

31 / 39

Matching Rates Depends on Service Distribution:

1 2 3

C S

1 2 3 1 2 3 4

C S

1 2 3 4

C S

1 2 3 4 5 6 1 2 3 4 5 6

) 2 ( ) 1 ( (3)

Figure: Topologies of the systems for the simulation study.

For each system, we used four service time distributions:

◮ System 1: α = (.2, .6, .2) and µ = (.4, .2, .4); ◮ System 2: α = (.1, .4, .4, .1) and µ = (.4, .3, .2, .1); ◮ System 3: α = (.1, .2, .2, .1, .2, .2) and

µ = (.05, .1, .15, .2, .2, .3);

32 / 39

◮ Exponential; ◮ Pareto (denoted by the subscript ‘p’), density

f(x) = 3γ(γx + 1)−4, x ≥ 0 with mean 1/2γ;

◮ Uniform 1: U(0, 2/µ) with µ = 2γ; ◮ Uniform 2: U(.9/µ, 1.1/µ) with µ = 2γ.

System Exponential Pareto Uniform 1 Uniform 2 1 .285 .299 .535 .074 2 .528 * .0078 * 3 .636 * * *

Table: Resulting p-values of Hotelling’s T 2 test. * denote p-value < 10−15.

33 / 39

System 1

r =   .1 .1 .3 .3 .1 .1   ,

• rp =

  .09996 .10006 .29987 .30013 .1 .09997  

• ru1 =

  .1 .09996 .30002 .30003 .10001 .09998   ,

• ru2 =

  .10002 .10009 .29999 .29991 .10002 .09998   (see Adan and Weiss, Exact FCFS matching rates for two infinite multi-type sequences. Operations Research, 60(2011), 475–489)

34 / 39

System 2

r =     .06443 .03557 .3356 .06443 .2356 .1644 .03557 .06443     ,

• rp =

    .06477 .03524 .33519 .0648 .23523 .16478 .03531 .06468    

• ru1 =

    .06447 .03553 .33549 .06447 .23556 .16447 .03554 .06446     ,

• ru2 =

    .06465 .03537 .33535 .06461 .2354 .16465 .03535 .06463    

35 / 39

System 3

r =         .004584 .009497 .08592 .03825 .06356 .09819 .02694 .04084 .1322 .01097 .02815 .06087 .03963 .06184 .09854 .007164 .07729 .1155        

• rp =

        .00462 .01039 .08499 .03853 .06318 .0983 .02638 .04062 .13298 .01098 .02787 .0612 .03923 .06149 .09927 .00681 .07741 .11575        

36 / 39

• ru1 =

        .00465 .00894 .08645 .03781 .06386 .09834 .02718 .04078 .132 .0109 .02819 .06092 .03979 .06198 .09823 .00754 .07713 .11531        

• ru2 =

        .00476 .00864 .08661 .0374 .06409 .09857 .02727 .04055 .13213 .01088 .02802 .0611 .03985 .06206 .09811 .00784 .07684 .11529        

37 / 39

In summary, by Hotelling’s T 2-test, the matching rates of systems 2 and 3 under the Pareto or uniform service time assumption are different from the exponential service time assumption in a statistically significant manner.

38 / 39

Summary:

◮ Stability for the parallel service system; ◮ Completely resource pooling; ◮ Matching rate; ◮ The first moment is not sufficient to determine the stability,

the completely resource pooling, and matching rate.

39 / 39