# Analysis of the Symmetric Join the Shortest Orbit Queue

This work introduces the join the shortest queue policy in the retrial setting. We consider a Markovian single server retrial system with two infinite capacity orbits. An arriving job finding the server busy, it is forwarded to the least loaded orbit. Otherwise, it is forwarded to an orbit randomly. Orbiting jobs of either type retry to access the server independently. We investigate the stability condition, the stationary tail decay rate, and obtain the equilibrium distribution by using the compensation method.

There are no comments yet.

## Authors

• 7 publications
04/16/2021

### The generalized join the shortest orbit queue system: Stability, exact tail asymptotics and stationary approximations

We introduce the generalized join the shortest queue model with retrials...
10/17/2018

### Load balancing with heterogeneous schedulers

Load balancing is a common approach in web server farms or inventory rou...
10/20/2016

### Non-Asymptotic Delay Bounds for Multi-Server Systems with Synchronization Constraints

Multi-server systems have received increasing attention with important i...
03/01/2019

### Stabilizing a Queue Subject to Action-Dependent Server Performance

We consider a discrete-time system comprising an unbounded queue that lo...
12/07/2016

### An Infinite Dimensional Model for a Many Server Priority Queue

We consider a Markovian many server queueing system in which customers a...
02/20/2020

### Queueing Subject To Action-Dependent Server Performance: Utilization Rate Reduction

We consider a discrete-time system comprising a first-come-first-served ...
01/16/2021

### Sensitivity of Mean-Field Fluctuations in Erlang loss models with randomized routing

In this paper, we study a large system of N servers each with capacity t...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this work, we introduce the concept of the join the shortest queue

(JSQ) policy in the retrial setting. More precisely, we consider a single server retrial system with two infinite capacity orbit queues. The service station can handle at most one job, and thus, arriving jobs that find the server busy, they join the least loaded orbit queue. In case both orbits have the same occupancy the job is routed to an orbit with probability

. Orbiting jobs retry independently to connect with the service station after some random time period. Our primary aim is to investigate the stationary behaviour of such a Markov-modulated two-dimensional random walk by using the compensation method.

The (non-modulated) two-dimensional JSQ problem was initially studied in [16, 19, 14], A compact mathematical method using generating functions was provided in [10, 13]. However, it does not lead to an explicit characterization of the equilibrium probabilities, and cannot be easily used for numerical purposes. In [4, 5, 6, 3], the authors introduced the compensation method (CM), an elegant and direct method to obtain explicitly the equilibrium join queue-length distribution as infinite series of product form terms. We particularly mention the work in [2], where Erlang arrivals fed a two-queue system under the JSQ policy. The queueing system in [2] is described by a multilayer random walk in the quarter plane. For other related important works, see also [1, 9, 26, 25]. Alternatively, numerical/approximation methods were applied. The power series algorithm (PSA) was applied in JSQ systems in [7, 8]; see also [12, 17] (non-exhaustive list) for more complicated models. We also mention the matrix geometric method; see e.g., [15, 23], for which connections with CM was recently reported in [18].

### Fundamental contribution:

In this work, we provide an exact analysis that unifies two queueing models: The JSQ model and the two-class retrial model. Our primary aim is to extend the applicability of the compensation method to random walks in the quarter plane modulated by a two-state Markov process, and in particular to retrial queueing systems with two infinite capacity orbit queues; see Section 2. In such a case, the phase process represents the state of the server, and affects the evolution of the level process, i.e, the orbit queue lengths in two ways: i) The rates at which certain transitions in the level process occur depend on the state of the phase process. Thus, a change in the phase might not immediately trigger a transition of the level process, but changes its dynamics (indirect interaction). ii) A phase change does trigger an immediate transition of the level process (direct interaction). For this modulated two-dimensional random walk we investigate its stationary behaviour by using the compensation method. We further study the ergodicity conditions and investigate its stationary tail decay rate.

### Application oriented contribution:

This work is also application oriented, since introduces the concept of JSQ in the retrial setting. Applications of this model can be found in relay-assisted cooperative communication system: There is a source user that transmits packets to a common destination node (i.e., the single service station), and a finite number of relay nodes (i.e., the orbit queues) that assist the source user by retransmitting its blocked packets, e.g., [11, 22]. The JSQ protocol serves as the cooperation strategy among the source and the relays, under which, the user chooses to forward its blocked packet to the least loaded relay node. This works serves as a first step towards the analysis of even general retrial models operating under the JSQ policy.

The paper is organized as follows. In Section 2 we describe the model in detail, and investigate the stability condition. The stationary tail decay rate for the JSQ with retrials is investigated in Section 3. The stationary behaviour of this Markov-modulated random walk in the quarter plane using CM is presented in Section 4. A numerical example is given in Section 5, and a conclusion is presented in Section 6.

## 2 Model description and stability condition

Consider a single server retrial system with two infinite capacity orbit queues. Jobs arrive at the system according to a Poisson process with rate

. If an arriving job finds the server idle, it starts service immediately. Otherwise, it is routed to the least loaded orbit queue. In case both orbit queues have the same occupancy, the blocked job is routed to either orbit with probability 1/2. Orbiting jobs of either type retry independently to occupy the server after an exponentially distributed time period with rate

, i.e., we consider the constant retrial policy. Service times are independent and exponentially distributed with rate . Denote by the number of jobs stored in orbit at time , and by the state of the server, i.e., , when it is busy and when it is idle at time , respectively. Let ,

, the dynamics of our system is described by the following continuous time Markov chain

, with state space . Our aim is to determine the equilibrium distribution

 qm,n(k)=limt→∞P((X1(t),X2(t),C(t))=(m,n,k)),(m,n,k)∈S.

Consider the column vector

, where denotes the transpose of vector (or matrix) . The equilibrium equations are written in matrix notation as follows

 A0,0q(0,0)+A0,−1q(0,1)=0, (1) B0,0q(0,1)+A0,−1q(0,2)+2A−1,1q(1,0)+A0,1q(0,0)=0, (2) B0,0q(0,n)+A0,−1q(0,n+1)+A−1,1q(1,n−1)=0,n≥2 (3) C0,0q(m,0)+A0,−1q(m,1)+A1,−1q(m−1,1)=0,m≥1 (4) (5) C0,0q(m,n)+A0,−1q(m,n+1)+A1,−1q(m−1,n+1)+A−1,1q(m+1,n−1)=0,m≥1,n≥2, (6)

where

 A1,−1=A0,1=(000λ),A0,−1=A−1,1=(00α0),A0,0=(−λμλ−(λ+μ)),B0,0=A0,0−H,C0,0=A0,0−2H,H=(α000).

The transition rate matrix of is given by

 Q=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝¯T0T1T−1T0T1T−1T0T1⋱⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟⎠,

where,

 ¯T0=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝AT0,0AT0,1AT0,−1BT0,0AT0,−1BT0,0⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,T1=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝OOAT1,−1OAT1,−1O⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,T0=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝CT0,0AT0,1AT0,−1CT0,0AT0,−1CT0,0⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,T−1=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝O2AT−1,1OAT−1,1OAT−1,1⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

so that is a quasi birth death (QBD) process with repeated blocks , , .

Using truncation arguments we show that is ergodiciff . To proceed, we modify the original system with the shortest queue discipline, such that the shortest orbit queue will not retry when the difference of the two orbit queues attains a predefined constant . Namely, the state is truncated as by removing the state transitions from to .

For this modified model, the longer orbit queue is set for the level instead of the shortest queue. Let the queue length in orbit for the truncated model, and . Since the difference between the queue lengths at the orbits is restricted by , we can construct sample paths such that , with probability one in case . Therefore, is positive recurrent if is positive recurrent for sufficient large . Thus, we focus on the investigation of ergodicity for the truncated model. Then, we have the following result.

###### Lemma 1

The truncated model (and so as the original one) is stable iff .

See Appendix A.

## 3 Decay rate

We first focus on the single server retrial queue with a single orbit of infinite capacity, and a limited classical retrial policy, called the reference model. We show that the decay rate of the stationary orbit queue length distribution of this model equals . Then, we will prove that the tail decay rate for the shortest orbit queue model (i.e., the original model) is .

In the reference model, jobs arrive according to a Poisson process with rate , and service times are exponentially distributed with rate . Arriving jobs that find the server busy, are routed to the orbit, from which, they retry to access the server according to the limited classical retrial policy: If there is only one job in orbit, it retries after an an exponentially distributed time with rate . If there are at least two jobs in orbit, the retrial rate changes to , i.e., limited classical retrial policy.

Let be the equilibrium queue length in the this special single server, single orbit, retrial queue. The following lemma, provides results regarding the stationary distribution of the reference model

###### Lemma 2

For , the stationary distribution of the reference model is given by:

 πm(0):=P(Q=m,C=0)=π1(0)λ+αλ+2αρi,m≥2,πm(1):=P(Q=m,C=1)=π1(0)λ+αμρi,m≥2,π1(0):=P(Q=1,C=0)=λ22αμ(1−ρ)2αμ(1−ρ)(λ2+α(λ+μ))+λ2(2α−λ)(λ+α),π1(1):=P(Q=1,C=1)=π1(0)λ+αμ,π0(0):=P(Q=0,C=0)=π1(0)αμλ2,π0(1):=P(Q=0,C=1)=π1(0)αλ. (7)
###### Proof 2

See Appendix B.

The following theorem summarizes the main result of this section.

###### Theorem 1

For , , we have , , are positive constants independent of .

###### Proof 3

See Appendix C.

Theorem 1, can be understood by comparing the original model with the reference model described above. Let , be the equilibrium queue lengths of the original model. Since both models work at full capacity whenever the total number of customers grows, we expect that , will have the same decay rate. We also expect that for increasing values of ,

 P(min{Q1,Q2}=m)≃P(Q1+Q2=2m)≃P(Q=2m), (8)

since the JSQ policy constantly aims at balancing the lengths of the two orbit queues over time. Having in mind the results of Lemma 2, equation (8) leads to the following conjectured behaviour of the minimum orbit queue length:

 P(min{Q1,Q2}=m,C=k)≃ykρ2m,as m→∞,

for a positive constant , . Therefore, the decay rate of the tail probabilities for is conjectured to be equal to the square of the decay rate of the tail probabilities of . Proposition 1 states this conjecture for fixed values of the server state and of the difference of the queue lengths.

## 4 The compensation approach

We conjecture that the inner balance equations have a product form solution. To show this, we construct a closely related model that has the same behaviour in the interior as the original model. The modified model is constructed as follows. Starting from the state space of the original model we bend the vertical axis such that the modified model has the same equilibrium equations in the interior and on the horizontal boundary.

###### Lemma 3

For , the equilibrium equations (4)-(6) have a unique up to a constant solution of the form

 q(m,n)=ρ2mu(n),m,n≥0, (9)

with non-zero such that , .

###### Proof 4

We consider a modified model, which is closely related to the original one described by and that is expected to have the same asymptotic behavior. This modified model is considered on a slightly different grid, namely .

In the interior and on the horizontal boundary, the modified model has the same transition rates as the original model. A characteristic feature of the modified model is that its balance equations for are exactly the same as the ones in the interior (i.e., the modified model has no “vertical” boundary equations) and both models have the same stability region. Therefore, the balance equations for the modified model are given by (4)-(6) for all , with only the equation for state , being different due to the incoming rates from the states with , .

Observe that the modified model, restricted to an area of the form embarked by a line parallel to the axis, yields the exact same process. Hence, we can conclude that the equilibrium distribution of the modified model, say , satisfies , and therefore

 ^qm,n=γm^q0,n, m≥−n, n≥0. (10)

We further observe that . To determine the term we consider levels of the form , and let . The balance equations among the levels are:

 C0,0^qL+A1,−1^qL−1+2A0,−1^qL+1=0, L≥1, (11)

Moreover, equation (10) yields

 (12)

Substituting (12) into (11) yields . Combining (12) with (11) with yields

 det(ρC0,0+A1,−1+ρ22A0,−1)=0⇔2αμ(1−ρ)(ρ−λ(λ+2α)2αμ)=0,

which implies that indeed .

Thus, it is shown that the equilibrium distribution of the modified model has a product-form solution which is unique up to a positive multiplicative constant. Returning to the original process , we immediately assume that the solution of (4)-(6) is identical to the expression for the modified model as given in (9). Moreover, the above analysis implies that this product-form is unique, since the equilibrium distribution of the modified model is unique.

In the following, we describe the form of the solution satisfying the inner equilibrium equations, and specify the form of .

###### Lemma 4

The product form , , , , satisfies the inner equilibrium equation (6) if

 (13)
###### Proof 5

The desired result is obtained directly by substituting the product form in (6). For , the matrix has rank equal to 1. From the system of linear equations , and the form of we derive the second in (13).

###### Remark 1

Note that the value of eigenvector

is independent of the values of , that satisfy (13).

The next step is to determine s and s such that , for which there exists a non-zero solution of (13), i.e., . The next lemma gives information about the location of the zeros of .

###### Lemma 5

For , and for every fixed with , the equation takes the form

 γδ2(ρ+1)−2ρδ2−γ2−γδ2=0, (14)

and has exactly one root in the plane such that . For , and for every fixed with , the equation (14) has exactly one root in the plane such that .

###### Proof 6

Starting from , simple algebraic computations yield in (14). Surprisingly, the form (14) is exactly the same as the one of equation (8) in Lemma 1 of the seminal paper [4], and thus, the assertions , can be proven similarly, so further details are omitted.

Lemmas 4, 5 characterize the basic solutions satisfying the inner equation (6). It is easily seen that the pair satisfy the equilibrium equations (5), (6). In the following, we construct a linear combination of elements from the basis of these solutions, which is a formal solution to the balance equations. The proof of Lemma 6 follows by the substitution of (15) in (4)-(6).

###### Lemma 6

The solution

 q(m,n)= {h0γm0δn0θ,m,n≥1,γm0ξ,m≥1,n=0, (15) ξ= −h0γ0C−10,0[A1,−1+γ0A0,−1]δ0θ, (16)

satisfies the equilibrium equations (4)-(6).

It is easily seen that the solution in (15) does not satisfy the vertical boundary equation (3). To compensate for this error we add a new term such that the sum of two terms satisfies (3), (6). Thus, we assume that satisfies both (3), (6). Substituting this form in (3) yields

 [h0V(γ0,δ0+c1V(γ,δ)]θ=0,n≥2, (17)

where . Hence, , and for such , we obtain from (13) , such that , so that satisfies (3). Thus, the solution satisfies (3). The constant can be obtained by substituting it in (3), or equivalently, by using (17) and the fact that , are the roots of (14). This procedure yields after some algebra

 c1=−γ1−δ0(λ+μμ)γ0−δ0(λ+μμ)h0. (18)

Adding the new term, we violate the horizontal boundary equations (4), (5). Thus, we compensate for this error by adding a product form generated by the pair , such that . The new solution is now given by

 q(m,n)={h0γm0δn0θ+c1γm1δn0θ+h1γm1δn1θ,m,n≥1,γm0ξ+γm1ξ1,m≥1,n=0, (19)

where , are obtained such that to satisfy (4)-(6). In particular, by substituting (19) to (4) yields

 ξ1=−1γ1[γ0ξ+C−10,0(A1,−1+γ1A0,−1)(c1δ0+h1δ1)θ]. (20)

Now, note that (4) reads

 q(m,0)=−C−10,0[A1,−1q(m−1,1)+A0,−1q(m,1)],m≥1.

Substituting back to (5) yields for ,

 [C0,0−(A0,1C−10,0A0,−1+2A−1,1C−10,0A1,−1)]q(m,1)+A1,−1q(m−1,2)+A0,−1q(m,2)−A0,1C−10,0A1,−1q(m−1,1)−2A−1,1C−10,0A0,−1q(m+1,1)=0. (21)

Substituting (19) in (21) yields after tedious algebra that should satisfy

 [h1L(γ1,δ1)+c1L(γ1,δ0)]θ=0,L(γ,δ):=δ[γ[C0,0−(A0,1C−10,0A0,−1+2A−1,1C−10,0A1,−1)]+δA1,−1+γδA0,−1−δA0,1C−10,0A1,−1−γ22A−1,1C−10,0A0,−1]. (22)

Thus, (22) implies that , and having in mind (14), we obtain after some algebra that

 h1=−(ρ+γ1)/δ1−(1+ρ)(ρ+γ1)/δ0−(1+ρ)c1. (23)

We continue in this manner until we construct the entire formal series,

 q(m,n)={∑∞i=0(hiγmi+ci+1γmi+1)δniθ,m≥0,n≥1,γm0ξ+∑∞i=1γmiξi,m≥1,n=0,q(0,0)=−A−10,0A0,−1q(0,1). (24)

We now have to show that the solutions (24) converge in two steps: to show that the sequences , converge to zero exponentially fast, and the formal solution converges absolutely. The following theorem summarizes the main result:

###### Theorem 2

For ,

 q(m,n)∝∑∞i=0(hiγmi+ci+1γmi+1)δniθ,(% pairs with the same δ-term),m≥0,n≥1,∝(h0γm0δn0+∑∞i=0(hi+1δni+ci+1δni+1)γmi+1θ,(%pairswiththesame$γ$−term)m≥0,n≥1,q(m,0)∝(γm0ξ+∑∞i=1γmiξi),m≥1, (25)

and , where the symbol () means “directly proportional”. Moreover, , , and the sequences , , , , , are obtained recursively based on the analysis above.

The next task is to show that formal solution in (25) converges absolutely. To show that, we need some preliminary results. Since (14) has the same form as in [4, eq. (8)], the sequences , satisfy .

###### Proposition 1

The sequences , in (25) satisfy: , and .

###### Proof 7

We first show that for a fixed , with , , and then, for a fixed , with , we have .
For a fixed , set on . Under this transform, (14) reads . Set , . Note that

 |f(z)|=2(1+ρ)|z|=1+ρ,|g(z)|=|(2ρ+|γ|)z2+1|≤(2ρ+|γ|)|z|2+1=(2ρ+|γ|)14+1.

Moreover, . Note that . Thus, since on , Rouché’s theorem [27] completes the proof of .
For a fixed , we show that , by setting now in (14) on the domain . Then, (14) reads, . Set , . Note that has a single zero in the interior of . Then, and

 |m(w)|≥23|2(1+ρ)−|δ|−3ρ|=23|2−ρ−|δ||=23(2−ρ−|δ|).

Note that, . This completes the proof that on . Rouché’s theorem [27] completes the proof of . Then applying , iteratively yields,

 |γi|≤23|δi−1|≤2312|γi−1|≤…≤(2312)i|γ0|=(13)iρ2,|δi|≤12|γi|≤1223|δi−1|≤…≤(2312)i|δ0|=(13)iρ23=(13)i+1ρ2.

Proposition 1 states that , as .

In the following, we focus on the asymptotic behaviour of and . This result is important to investigate the convergence of the series in (25).

###### Lemma 7

Let fixed and the root in (14) such that . As , then , the smallest root of

 2ρw2−2(1+ρ)w+1=0. (26)

Let fixed and the root in (14) such that . A , then , with the larger root of (26).

###### Proof 8

See [4], since (14) that generates the sequences , has the same form as in [4, eq. (8)].

The final ingredient to check the convergence of the series (25) is to determine the values of the ratios , as .

###### Lemma 8
1. Let , , be roots of (14) such that . Then, as , .

2. Let , , be roots of (14) such that . Then, as , .

3. As , .

4. For , and , the vector , is such that , as .

###### Proof 9

Using the indexing of the compensation parameters (18),

since as , Lemma 7 implies that , , and where the last equality follows from .
Similarly, from (23), we have .

Since as , , , and