Introduction

Fork-join queueing systems are widely used to model parallel and distributed systems. Some examples of these systems include RAID, distributed web applications, MapReduce, GRID, multipath routing networks, multiprocessor systems and etc. In these systems, arriving jobs are split for service by numerous servers and joined before departure.

In papers [1,2,3], fork-join queueing systems are applied to the analysis of MapReduce clusters and multipath routing. Results on modeling of RAID arrays and distributed databases are provided in [4,5]. Papers [6,7] describe an application of fork-join queueing systems to the analysis of multiprocessor systems and parallel algorithms. Fork-join queueing systems also arise in the performance analysis of automated manufacturing systems. In these systems, parallelism can be found in processes such as product assembly and logistic operations that involve multiple suppliers.

Paper [4] presents an overview of the main results for fork-join and related queueing models. Most of the papers consider parallel-server fork-join queueing systems [8,9,10,11,12].

The parallel-server fork-join queueing system consists of homogeneous servers [9,10,13] or heterogeneous servers [1,2,8,11,12,14]. An arriving job is split into a number of tasks, one for each server. Previous research in this area reports three basic types of parallel-server fork-join queueing systems [15], such as the centralized splitting model, the distributed splitting model, and the split-merge model.

In [11] upon arrival, a job of size S bringing tasks to S ≤ K of K servers is immediately split so that each of its S tasks is allocated to exactly one server. In [14] each job is split into a number of tasks, one for each free server. Paper [2] considers both deterministic and probabilistic scheduling strategies. Within probabilistic strategy, each task chooses a server with some probability. Instead, the present paper considers the allocation of tasks between the heterogeneous servers according to a threshold control policy. Under a threshold policy, a task is accepted at server i if and only if the number of tasks in queue is at or above a given threshold q i before acceptance.

The threshold control policy is used as optimal solution for the job routing problem named the slow server problem [16,17]. This problem focuses on routing jobs to heterogeneous servers so as to minimize the average response time of jobs or to minimize the average number of jobs in the system. As it shown in [16], the optimal policy which minimizes the average response time of jobs in the system with two heterogeneous servers is of threshold type, i.e., the fast server should always be used (as long as there are jobs in the queue), and that the slower server should only be used if the number of jobs in queue exceeds a certain threshold value. The generalization of this rule was obtained in [17,18] for the case of multiserver queueing system with heterogeneous servers. The optimal control policy has the monotonicity and the threshold property, i.e., there exists a queue length threshold such that a new server (with the minimal service cost and with the maximal service rate) should be used only after the queue reaches this threshold. Later, the results were extended to retrial queueing systems with heterogeneous servers [19], systems with unreliable servers [20] and with heterogeneous servers in speed and in quality of resolution [21].

In this paper, we consider a centralized splitting fork-join queueing model [15] with heterogeneous servers and a threshold control policy. We apply the matrixgeometric approach to analyze the considered system. A similar approach was used in [10,14].

The remainder of the paper is organized as follows. Section 2 describes the fork-join queueing system under consideration. In Sections 3 and 4 we obtain the stationary distribution and the main performance measures. Section 5 provides an illustration of the above results, we also discuss some control policies and compare them. Finally, a section of conclusions commenting the main research contributions of this paper is presented. We consider a system which consists of M heterogeneous servers S i , i = 1, . . . , M , and a FCFS infinite capacity queue as shown in Fig. 1. Jobs arrive at the system according to a Poisson process with rate λ. Each job is split into d homogeneous tasks, d ≤ M . An arriving task is placed in the queue. The service times at server S i have exponential distribution with parameter µ i , i = 1, . . . , M , and

µ 1 > µ 2 > • • • > µ M .

Let b denote the number of waiting tasks in the queue (queue length). At the arrival times of new tasks the control consists in sending tasks from the queue to the fastest idle server or leaving them in the queue. A task from the front of the queue is assigned to a fastest idle server S i , i ∈ {1, . . . , M }, if the queue length immediately after the arrival is not less than threshold level q i , i.e., When a task finishes its service at server S i , a task from the queue will be assigned to server S i if the queue length satisfies b ≥ q i . Otherwise server S i will be idle. Tasks of a job are served independently, when all the tasks associated with the job are finished, the job service is completed.

q i ≤ b + 1, (1 = q 1 ≤ q 2 ≤ • • • ≤ q M < q M +1 = ∞). µ 1 µ 2 . . . µ M . . . λ

The Stationary Distribution

The system state is described by a vector x = (r, n 0 , n 1 , . . . , n M ), where r = r(b) = b div d denotes the number of waiting jobs in the queue, n 0 = n 0 (b) = b mod d is the number of tasks for the served job in the queue,

n i = 0, if server S i is idle, 1, if server S i is busy.

The process {x(t), t ≥ 0} is a (M + 2)-dimensional continuous-time Markov chain on the state space X .

For each state x, we denote by F(x) and F(x) the sets of idle and busy server indices, respectively,

F(x) = {i > 0 : n i = 0}, F(x) = {i > 0 : n i = 1}.

Denote by γ(b) the mandatory (minimal) number of busy servers for queue length b. We can write

γ(b) =      0, b = 0, M, b ≥ q M , max{i : q i ≤ b}, otherwise. Define ϕ(x) by ϕ(x) = the minimal element of F(x), if F(x) = ∅, ∞, otherwise.

If there are idle servers in the system, a task will be assigned to server S ϕ(x) .

Let us define now the transition rate q(x, x ) from state x ∈ X to state x ∈ X . Each transition is associated with an event in the system, there are two types of events.

A job arrives and splits into d tasks(a) if b ≥ q M , q (x, (r + 1, n 0 , n 1 , . . . , n M )) = λ,(1)(b) if b < q M , q(x, x ) = λ,(2)

where x = x (d) can be obtained from the following recurrence relation

x (i+1) = r(b (i) + 1), n 0 (b (i) + 1), n (i) 1 , . . . , n (i) M , ϕ(x (i) ) > γ(b (i) + 1), x (i) + e 2+ϕ(x (i) ) ,

otherwise, with the initial value x (0) = x, i = 0, . . . , d − 1.

Here,

x (i) = r (i) , n (i) 0 , . . . , n (i) M and b (i) = dr (i) + n (i)

0 ; e j denotes the vector of the appropriate dimension with all components equal to 0, except the jth, which is 1.

A task finishes its service at serverS i (n i = 1) (a) if b < q i , q (x, (r, n 0 , n 1 , . . . , n i−1 , 0, n i+1 , . . . , n M )) = µ i ,(3)(b) if n 0 > 0, {i ∈ F(x) : b ≥ q i } = ∅, q (x, (r, n 0 − 1, n 1 , . . . , n M )) = i∈ F (x): b≥qi µ i ,(4)(c) if r > 0, n 0 = 0, {i ∈ F(x) : b ≥ q i } = ∅, q(x, (r − 1, d − 1, n 1 , . . . , n M )) = i∈ F (x): b≥qi µ i .(5)

Let state space X be lexicographically ordered. The subset X l is called level l, l = 0, 1, . . . , and defined as

X l = {(r, n 0 , n 1 , . . . , n M ) ∈ X : r = l}. The cardinality L of X 0 is L = 2 M + (q 2 − q 1 )2 M −1 + (q 3 − q 2 )2 M −2 + • • • + (d − q σ0 )2 M −σ0 , where σ 0 = γ(d − 1).

Let r be the minimal number r ∈ {1, 2, . . . } such that rd ≥ q M . The cardinalities of X l , l ≥ r are d.

The cardinalities of X l , 0 < l < r are

K l = 2 M −σ − l (q σ − l +1 − ld) + 2 M −σ − l −1 (q σ − l +2 − q σ − l +1 ) + • • • + + 2 M −σ + l +1 (q σ + l − q σ + l −1 ) + 2 M −σ + l (ld + d − q σ + l ), where σ − l = γ(ld), σ + l = γ(ld + d − 1)

. The system states are presented in Table 1.

M q1 ≤ n0 < q2 1, n2, . . . , nM (q2 − q1)2 M −1 . . . . . . . . . qσ 0 ≤ n0 ≤ d − 1 1, . . . , 1, nσ 0 +1, . . . , nM (d − qσ 0 )2 M −σ 0 0 < l < r 0 ≤ n0 < q σ − l +1 − ld 1, . . . , 1, n σ − l +1 , . . . , nM 2 M −σ − l (q σ − l +1 − ld) q σ − l +1 − ld ≤ n0, n0 < q σ − l +2 − ld 1, . . . , 1, n σ − l +2 , . . . , nM 2 M −σ − l −1 (q σ − l +2 − q σ − l +1 ) . . . . . . . . . q σ + l − ld ≤ n0 < d 1, . . . , 1, n σ + l +1 , . . . , nM 2 M −σ + l (ld + d − q σ + l ) l ≥ r 0 ≤ n0 ≤ d − 1 1, . . . , 1 d

It is easily shown that for X l , l > 0, there are transitions to X l−1 and X l+1 . Thus, the Markov chain {x(t), t ≥ 0} is a quasi-birth-and-death (QBD) process, the generator matrix Q of the process has the following form

Q =                   B 00 B 01 0 0 • • • 0 0 0 0 0 0 0 • • • B 10 B 11 B 12 0 • • • 0 0 0 0 0 0 0 • • • 0 B 21 B 22 B 23 • • • 0 0 0 0 0 0 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 • • • 0 B r−1,r−2 B r−1,r−1 B r−1,r 0 0 0 • • • 0 0 0 0 • • • 0 0 B r,r−1 B r,r A 0 0 0 • • • 0 0 0 0 • • • 0 0 0 A 2 A 1 A 0 0 • • • 0 0 0 0 • • • 0 0 0 0 A 2 A 1 A 0 • • • 0 0 0 0 • • • 0 0 0 0 0 A 2 A                 

, where blocks A 2 , A 1 , A 0 are square matrices of order d, B 00 is the square matrix of order L, blocks B 01 , B 10 are L × K 1 , K 1 × L matrices respectively, B l,l , l = 1, . . . , r, are square matrices of order K l (we set K r = d), blocks B l,l−1 , l = 2, . . . , r, are K l × K l−1 matrices, and blocks B l,l+1 , l = 1, . . . , r − 1, are K l × K l+1 matrices. The elements of matrices B 01 , B l,l+1 , l = 1, . . . , r −1, and A 0 are determined by ( 1) and ( 2). The elements of matrices B 10 , B l,l−1 , l = 2, . . . , r, and A 2 are determined by (5). Note that A 0 = λI, where I is an identity matrix of the appropriate order; A 2 is the square matrix of order d, where the (1, d)th element of the matrix is equal to M i=1 µ i , and other elements are zero. The non-diagonal elements of matrices B 00 , B l,l , l = 1, . . . , r, and A 1 are determined by ( 2)-(4); the diagonal elements of these matrices are determined from conditions (we set B r,r+1 = A 0 ):

B 00 1 + B 01 1 = 0, B l,l−1 1 + B l,l 1 + B l,l+1 1 = 0 , A 0 1 + A 1 1 + A 2 1 = 0,

where 1 is an ones column vector of the appropriate order. Note that A 1 is the bidiagonal square matrix of order d, where the (j, j)th element of the matrix is equal to −(λ + M i=1 µ i ), and the (j, j − 1)th element is equal to M i=1 µ i . It is known [22], that the QBD process is ergodic if and only if

π A A 0 1 < π A A 2 1, where π A A = 0, π A 1 = 1, A = A 0 + A 1 + A 2 . The condition implies λd < µ 1 + • • • + µ M .

The stationary distribution π = (π 0 , π 1 , . . . ) of the QBD process can be computed using the matrix-geometric method [22]. The row vector π l , l = 0, 1, . . . , defines probabilities of states for level l, according to the lexicographic order, π l = (π(x) : x ∈ X l ). We solve linear system πQ = 0 and π1 = 1 to find the stationary probabilities. Taking the advantage of the block structure in Q, we obtain

π 0 B 00 + π 1 B 10 = 0; π 0 B 01 + π 1 B 11 + π 2 B 21 = 0; π 1 B 12 + π 2 B 22 + π 3 B 32 = 0;

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

π i−1 B i−1,i + π i B i,i + π i+1 B i+1,i = 0; . . . . . . . . . . . . . . . . . . π r−2 B r−2,r−1 + π r−1 B r−1,r−1 + π r B r,r−1 = 0; π r−1 B r−1,r + π r B r,r + π r+1 A 2 = 0; π l−1 A 0 + π l A 1 + π l+1 A 2 = 0, l = r + 1, r + 2, . . . .(6)

We know [22], the solution of ( 6) has the matrix-geometric form given by

π r+l = π r R l , l = 1, 2, . . . ,

where R is the minimal non-negative solution to equation

A 0 +RA 1 +R 2 A 2 = 0, vectors π 0 , π 1 , . . . , π r satisfy (π 0 , . . . , π r )         B 00 B 01 0 • • • 0 0 B 10 B 11 B 12 • • • 0 0 0 B 21 B 22 • • • 0 0 • • • • • • • • • • • • • • • • • • 0 0 0 • • • B r−1,r−1 B r−1,r 0 0 0 • • • B r,r−1 B r,r + RA 2         = (0, . . . , 0), π 0 1 + π 1 1 + • • • + π r−1 1 + π r (I − R) −1 1 = 1.

Performance Measures

Once stationary distribution π is computed, a variety of other performance measures may be obtained.

The average number B of jobs in the queue

B = ∞ l=1 lπ l 1 = r l=1 lπ l 1 + ∞ l=1 (r + l)π r R l 1 = = r l=1 lπ l 1 + π r R r(I − R) −1 + (I − R) −2 1.

The average job waiting time W in the queue (the average time that the job waits in the queue before the first of its tasks is assigned to a server)

W = B/λ.

Denote by T the average response time, thenT = W + V,

where V is the average job service time, which is defined as the average total time required to process all of the tasks associated with the job once the first task is assigned to a server.

The average number N of jobs in the system is

N = λT.

To obtain the average job service time, we consider the service process of all d tasks of an arbitrarily chosen "tagged" job. There are two cases considered.

1. We assume that when the first task is assigned, the queueing system is in a state x such that r > r. Then all servers are busy, and tasks of the tagged job sequentially occupy available servers. The service process of the tagged job is described by an absorbing Markov chain ν = {ν(t), t ≥ 0}. Denote the state of chain ν by (k, S), where k denotes the number of tasks of the tagged job in the queue, S denotes the set of indices for busy servers that process tasks of the tagged job. State (0, ∅) is absorbing. The state space N of Markov chain ν is given by

N = d j=1 N j {(0, ∅)},

where N j = {(k, S) : k ∈ {0, 1, . . . , d − j}, S ∈ S j }, S j denotes the set of all j-element subsets of {1, . . . , M },

|N | = d j=1 j M d − j + 1 + 1.

The initial probability distribution defined on state space N is given by the vector (α 0 , α), where

α 0 = α 0 (0, ∅) = 0, α = (α(k, S) : (k, S) ∈ N \ {(0, ∅)}), α (d − 1, {i}) = µ i M i=1 µ i , i = 1, . . . , M,

and other initial probabilities are zero. Markov chain ν is determined by generator matrix Q ν with non-diagonal elements q ν ((k, S), (k , S ))

q ν ((k, S), (k , S )) =          µ i , if k = k = 0, S = S \ {i}, i ∈ S, µ i , if k = k − 1 ≥ 0, S = S {i}, i / ∈ S,

We denote by T the square matrix of order |N | − 1 which is a submatrix of matrix Q ν and contains transition rates among the transient states.

The absorption time ξ for Markov chain ν is a random variable which has a PH-distribution with PH-representation (α, T ) [22]. The expected absorption time is

E[ξ] = −αT −1 1.

2. We assume that when the first task is assigned, the queueing system is in a state x such that 0 ≤ r ≤ r. Then service of the tagged job starts right after the events:

the job arrives, -a task finishes and leaves the system, -a job arrives and this leads to a task of the job is assigned to an idle server. In this case, the service process of the tagged job is described by an absorbing Markov chain ζ = {ζ(t), t ≥ 0}. Denote the state of chain ζ by (k, S, x), where k denotes the number of tasks of the tagged job in the queue, S denotes the set of indices for busy servers that process tasks of the tagged job, x denotes the system state, x ∈ X , X = {x ∈ X : r ≤ r}. States (0, ∅, x), x ∈ X are absorbing. The state space Z of Markov chain ζ is given by The initial probability distribution defined on state space Z is given by the vector (β 0 , β), where β 0 = (β 0 (0, ∅, x) : x ∈ X ) = 0, β = (β(k, S, x ) : (k, S, x ) ∈ Z \ Z 0 ):

Z = d−1 k=0 Z k ,β(k, S, x ) =   λ x∈A(x ) π(x) + M i=1 µ i x∈Di(x ) π(x)   G −1 ,

x ∈ X presents states from which the system moves to state x as a result of a job arrival (x ∈ A(x )), or as a result of a task completion at server S i (x ∈ D i (x )), i = 1, . . . , M ; G denotes the normalizing constant. The Markov chain ζ is determined by generator matrix Q ζ with non-diagonal elements q ζ ((k, S, x), (k , S , x )).

Suppose there is a transition from x to x as a result of an event described in Section 3. So we have the corresponding transition from (k, S, x) to (k , S , x ), where k and S will be described below.

Assume the system moves from state x to state x as a result of a job arrival, we denote by H(x, x ) = F(x ) \ F(x) the index set of idle servers which will be busy right after the event. Let H

(x, x ) = {i 1 , i 2 , . . . , i m }, i 1 < i 2 < • • • < i m , then H j (x, x ) ⊆ H(x, x ) denotes the j-element set H j (x, x ) = {i 1 , . . . , i j }.

Thus, we can write (a) A job arrives

q ζ ((k, S, x), (k , S , x )) = λ, where k = max {k − |H(x, x )|, 0}, S = S H k−k (x, x ), (b) A task finishes its service at S i i. if k > 0, n i = 1, i / ∈ S, q ζ ((k, S, x), (k − 1, S ∪ {i}, x )) = µ i , ii. if k > 0, q ζ ((k, S, x), (k − 1, S, x )) = i∈S:n i =1 µ i , iii. if k > 0, n i = 0, q ζ ((k, S, x), (k, S \ {i}, x )) = µ i , iv. if k = 0, q ζ ((0, S, x), (0, S \ {i}, x )) = µ i .

We denote by η the absorption time for Markov chain ζ. The expected ab-

sorption time E[η] is computed similarly to E[ξ].

Finally, we have

V = E[ξ]c 1 + E[η]c 2 c 1 + c 2 ,

where

c 1 = ω 1 , c 2 = β1G, (ω 1 , . . . , ω d ) = π r (I − R) −1 − π r − π r+1 M i=1 µ i .

Consider a fork-join queueing system which consists of M = 3 heterogeneous servers, where µ 1 = 5, µ 2 = 2, µ 3 = 1.

The performance measures depend on threshold levels q i , i = 2, . . . , M . Let q * i be the optimal threshold levels with respect to the minimization of the average response time T . They can be obtained by the exhaustive search method. Assume that λ = 1 and d = 3, then the optimal threshold levels are q * 2 = 2, q * 3 = 8, and T = 0.838.

To illustrate the advantages of the threshold control policy, we consider the fastest free server policy with q i = 1, i = 1, . . . , M . Figure 2 presents the average response times for the optimal threshold control policy and the fastest free server policy for several values of the arrival rate λ. It may be seen that the optimal threshold control policy provides lower values of the average response time T . When the arrival rate λ is quite large, the difference between the policies can be neglected. As λ increases, the threshold levels decrease that leads to the fastest free server policy is near optimal one in this case.

Conclusions

We studied a Markovian fork-join queueing system with heterogeneous servers and threshold control policy. Applying a matrix-geometric approach, we obtained the stationary distribution of the system and performance measures. The results can be used for the performance analysis of multiprocessor systems and other modern distributed systems. An interesting topic for future research is to obtain the optimal threshold control policy which minimizes the long-run average number of jobs in the system or the average response time.