# A 3-Queue Polling System with Join the Shortest - Serve the Longest Policy

Efrat Perel1, Nir Perel1111Corresponding author, e-mail<EMAIL_ADDRESS>and
Uri Yechiali2

###### Abstract

In 1987, J.W. Cohen analyzed the so-called Serve the Longest Queue (SLQ)
queueing system, where a single server attends two non-symmetric $M/G/1$-type
queues, exercising a non-preemptive priority switching policy. Cohen further
analyzed in 1998 a non-symmetric 2-queue Markovian system, where newly
arriving customers follow the Join the Shortest Queue (JSQ) discipline. The
current paper generalizes and extends Cohen’s works by studying a combined
JSQ-SLQ model, and by broadening the scope of analysis to a non-symmetric
3-queue system, where arriving customers follow the JSQ strategy and a single
server exercises the preemptive priority SLQ discipline. The system states’
multi-dimensional probability distribution function is derived while applying
a non-conventional representation of the underlying process’s state-space. The
analysis combines both Probability Generating Functions and Matrix Geometric
methodologies. It is shown that the joint JSQ-SLQ operating policy achieves
extremely well the goal of balancing between queue sizes. This is emphasized
when calculating the Gini Index associated with the differences between mean
queue sizes: the value of the coefficient is close to zero. Extensive
numerical results are presented.

1School of Industrial Engineering and Management,
Afeka Tel-Aviv Academic College of Engineering, Tel-Aviv, Israel.

2Department of Statistics and Operations Research,
School of Mathematical Sciences, Tel-Aviv University, Tel-Aviv, Israel.

#### Keywords

polling systems; join the shortest queue; serve the longest queue; probability
generating functions; matrix geometric; Gini index

##  1 Introduction

Queueing systems with multiple queues have been studied extensively in the
literature. Noted are models where, upon arrival, customers join the shortest
queue (JSQ), along with other models where the server serves the longest queue
(SLQ). The main goal of each model-type is to balance between the queue sizes:
in the JSQ model, each customer’s action is directed into achieving the
balance, while in the SLQ model, the server’s actions are aimed at this goal.

Cohen, (1987) studied a single server system with two $M/G/1$-type queues,
under a non-preemptive Serve the Longest Queue (SLQ) regime. Each queue has
its own generally distributed service time and its own Poisson rate of
arrival, where upon each service completion, the server selects the longer
queue. If both queues are equal in size, the server chooses the next queue to
be served according to some probability. This model is analyzed by
constructing and solving a boundary value problem. The SLQ regime was further
studied by Flatto, (1989) who analyzed a Markovian system with two identical
queues and server’s preemptive switching policy, by calculating a
2-dimensional probability generating function of the queue lengths and
studying its asymptotic behavior. Houtum et al., (1997) analyzed a non-
preemptive Markovian system with $N$ symmetric queues, assuming that the
difference between the longest and the shortest queue is limited by some known
constant. As a result, it became possible to perform exact analysis of the
model and its variants via matrix geometric methods. Knessl and Yao, (2013)
derived asymptotic properties for a Markovian non-symmetric two-queue model in
heavy traffic, in which the server always attends the longer queue, and in
case when both queues are with equal sizes, the server splits its capacity
between the queues. An exact integral representation is derived for the joint
probability distribution of the queue sizes. Ravid et al., (2013) studied a
single repairman problem with two Markovian queues having non-identical
arrival rates, and pointed their model’s relation to a corresponding SLQ
model. The SLQ model was also studied in the context of wireless
communications and switching scheduling. For example, Maguluri et al., (2014)
analyzed a wireless network with SLQ mechanism, while Pedarsani and Walrand,
(2016) investigated the stability of an open multi-class network under SLQ
modeling.

The Join the Shortest Queue (JSQ) model has been also investigated widely in
the literature. Winston, (1977) showed that in a symmetric Markovian system
with identical servers, the JSQ policy maximizes the discounted number of jobs
completing service by some time $t$. Halfin, (1985) considered a service
center with two parallel single-server queues and identical Exponential
service rates, and derived, using linear programming techniques, bounds for
the probability distribution of the number of customers in the system. Hordijk
and Koole, (1990) showed, using Markov Decision Processes, that the JSQ policy
is optimal for queues with finite buffers and a general arrival process. Adan
et al., 1991a  studied a system with two non-symmetric parallel queues, and
showed, by using an elementary compensation procedure, that the equilibrium
distribution of the lengths of the two queues can be exactly represented by an
infinite sum of product form solutions. In a following paper Adan et al.,
1991b  considered the same system with jockeying between the queues, and
analyzed the model using a matrix geometric approach. Cohen, (1998) presented
an analytic solution for an asymmetric two-server queueing model with Poisson
arrivals, where arriving customers join the shorter queue and service times
are Exponentially distributed, by determining explicitly the bivariate
generating function of the stationary joint distribution of the queue lengths.
Foley and McDonald, (2001) derived asymptotic results for a Markovian system
with two servers, each of which receives a dedicated stream of customers in
addition to customers from a stream of smart customers who join the shorter
queue. Yao and Knessl, (2005, 2006) provided an asymptotic analysis for a
model with two $M/M/\infty$ queues where a new arrival is routed to the queue
with the smaller number of customers, by establishing and analyzing a bi-
variate generating function for the number of customers in each queue. Blanc,
(2009) studied a non-symmetric multi-server Markovian system where customers
follow the JSQ policy, and quantified the probability that a customers’s
service would have started earlier if one had joined another queue than the
queue that was actually chosen. By using the compensation procedure, Adan et
al., (2013) analyzed a system with two single-server queues, JSQ policy, where
customers inter-arrival times follow an Erlang distribution and service times
are Exponentially distributed. Adan et al., (2016) analyzed the JSQ polling
model by using the compensation method and by solving a boundary value
problem. Eschenfeldt and Gamarnik, (2018) and Braverman, (2020) analyzed the
steady-state properties of the JSQ model under the Halfin-Whitt regime. Dawson
et al., (2019) studied a queueing network with a large number of nodes, in
which each queue has a dedicated input stream, while, in addition, there is an
extra input stream, balancing the network load by directing its arrivals to
the shortest queue. A mean field interaction model is set up to study the
performance of this network in terms of limit results. Dimitriou, (2021)
studied a communication-related Markovian JSQ model with two infinite capacity
orbit queues. An arriving job that finds the server busy is forwarded to the
smaller-sized orbit queue. It is worth noting that in most of the JSQ models
described above, _each queue has its own server_. In contrast, in this work,
the combined JSQ-SLQ model analyzes a single-server polling-type system, where
the server switches from a shorter queue to a longer one, while arriving
customers join the shortest queue.

Recently, Perel et al., (2020) presented a combined unifying Join the shortest
queue - Serve the longest Queue (JSQ-SLQ) polling-type model under a
preemptive priority switching policy of the server. Such a combined policy
aims at equalizing the queue sizes. Specifically, the authors considered a
polling system comprised of two non-identical Markovian queues, attended by a
single server that alternates between them. An arriving customer always joins
the shortest queue, unless the queue lengths are equal, in which case the
customer joins queue $i$ with probability $p_{i}\geq 0$, $i=1,2$,
$p_{1}+p_{2}=1$. On the other hand, the server always serves the longest queue
while exercising a queue-size depending preemptive priority policy, i.e. the
server never resides in a shorter queue, giving priority to the longest queue.
This preemptive-type policy implies that at a moment when the number of
customers in an un-served queue exceeds the number of customers in the served
queue (either at a service completion or when arrival occurs), the server
immediately switches to the longest queue.

The current paper also combines the above two operating policies, i.e. JSQ and
SLQ, into a unified model, and extends the analysis to a 3-queue system.
Specifically, we study a polling system comprised of three asymmetrical
Markovian queues, denoted by $Q_{1}$, $Q_{2}$ and $Q_{3}$, and a single server
that alternates between them according to a preemptive SLQ policy. That is,
the server always attends the longest queue. As soon as the number of
customers in one of the non-attended queues rises above the number of
customers in the attended queue, the server stops serving the served customer
and instantaneously switches to the longest queue. The service of the
interrupted customer will resume anew when its turn comes again (within each
queue the service order is FCFS). In addition, if upon service completion the
two non-attendant queues ($Q_{i}$ and $Q_{j}$) have the (same) longest queue
length, the server will switch to $Q_{i}$ $(Q_{j})$ w.p.
$\frac{q_{i}}{q_{i}+q_{j}}$ $\left(\frac{q_{j}}{q_{i}+q_{j}}\right)$, where
$q_{i}>0$, $i=1,2,3$, is the weight the server assigns to $Q_{i}$. It is
assumed that service duration of an arbitrary customer in $Q_{i}$ is
Exponentially distributed with mean $1/\mu_{i}$, $i=1,2,3$. Customers arrive
at the system according to a Poisson process with rate $\lambda$. Each
arriving customer exercises the Join the Shortest Queue (JSQ) policy, whereas
if the lengths of all queues are equal, the customer joins $Q_{i}$ w.p.
$p_{i}$, where $p_{1}+p_{2}+p_{3}=1$. Furthermore, if only two queues, $Q_{i}$
and $Q_{j}$ ($i\neq j$), have the same shortest length, an arriving customer,
being consistent, joins $Q_{i}$ $(Q_{j})$ w.p. $\frac{p_{i}}{p_{i}+p_{j}}$
$\left(\frac{p_{j}}{p_{i}+p_{j}}\right)$ (note that according to this setting,
at most one of the $p_{i}^{\prime}s$ can be zero).

Even though we consider here a 3-queue system, where each queue is un-bounded
in its size, we are able to formulate this model as a 2-dimensional
continuous-time Markovian process with only one un-bounded size. In many
cases, solving a 2-dimensional un-bounded system may be done by applying a
boundary value problem analysis (see e.g. Flatto, (1989), Avrachenkov et al.,
(2014), Adan et al., (2016)), or by a truncation method as carried out in
Bright and Taylor, (1995). In contrast, using an un-conventional approach, as
was introduced in Perel et al., (2020), we are able to derive the equilibrium
joint probability mass function of the three queue lengths. For that purpose,
we apply both the two commonly used methods to solve QBD problems, namely, the
probability generating functions method and the matrix geometric analysis. The
main innovation of our analytical approach is to define only a finitely-many
number of probability generating functions representing the system’s states.

The paper is organized as follows. In Section 2 the model is described and
formulated. In Section 3 the steady-state equations are established and the
corresponding probability generating functions are derived. Various
performance measures are defined and calculated. In Section 4 we utilize the
matrix geometric method in order to derive the system’s stability condition,
as well as the equilibrium distribution of the system’s states. Extensive
numerical results are presented in Section 5, indicating that the balance
between the 3 queue sizes is very well achieved. This is further emphasized by
the corresponding Gini index that its value is close to zero in almost all
cases examined. Section 6 concludes the paper.

##  2 The Model

We consider a polling system with a single server and three non-identical
queues, denoted by $Q_{1}$, $Q_{2}$ and $Q_{3}$. Customers arrive at the
system according to a Poisson process with rate $\lambda$. Each arriving
customer exercises the ’Join the Shortest Queue’ (JSQ) policy, whereas, if the
lengths of all queues are equal, the customer joins $Q_{i}$ w.p. $p_{i}\geq
0$, where $p_{1}+p_{2}+p_{3}=1$. Furthermore, if only two queues, $Q_{i}$ and
$Q_{j}$ ($i\neq j$), have the same shortest length, an arriving customer,
being consistent, joins $Q_{i}$ $(Q_{j})$ w.p. $\frac{p_{i}}{p_{i}+p_{j}}$
$\left(\frac{p_{j}}{p_{i}+p_{j}}\right)$ (note that according to this setting,
at most one of the $p_{i}^{\prime}s$ can be zero). Service duration of an
arbitrary customer in $Q_{i}$ is Exponentially distributed with mean
$1/\mu_{i}$, $i=1,2,3$. The server alternates between the three queues
according to the preemptive ’Serve the Longest Queue’ (SLQ) switching policy.
That is, the server always attends (one of) the longest queue(s). As soon as
the number of customers in one of the non-attendant queues rises above the
number of customers in the attendant queue, the server stops serving the
served customer and instantaneously switches to the longest queue. The service
of the interrupted customer will resume anew when its turn comes again (within
each queue the service order is FCFS). In addition, if upon service completion
the two non-attendant queues ($Q_{i}$ and $Q_{j}$) have the (same) longest
queue length, the server will switch to $Q_{i}$ $(Q_{j})$ w.p.
$\frac{q_{i}}{q_{i}+q_{j}}$ $\left(\frac{q_{j}}{q_{i}+q_{j}}\right)$, where
$q_{i}>0$, $i=1,2,3$, is the weight assigned to $Q_{i}$. If, upon service
completion at $Q_{i}$ all queues are equal in size, the server stays in
$Q_{i}$.

At time $t>0$, let $L_{i}(t)$ denote the number of customers present in
$Q_{i}$, $i=1,2,3$, and, assuming stability (see in the sequel analysis of the
stability condition), let $L_{i}=\lim\limits_{t\rightarrow\infty}L_{i}(t)$.
The conventional approach would be to define a 3-dimensional un-bounded space
of the system’s states $\\{L_{1}(t),L_{2}(t),L_{3}(t)\\}$ which seems to be
untractable. Instead, we construct a non-conventional state space (following
Perel et al., (2020)) as follows: Let $I(t)$ denote the index of the queue
attended by the server at time $t$, and define the triple
$D(t)=\left(I(t);D_{2}(t),D_{3}(t)\right)$, where
$D_{j}(t)=L_{1}(t)-L_{j}(t)$, $j=2,3$. Let
$D=\lim\limits_{t\rightarrow\infty}D(t)$.

We now formulate the 3-dimensional intriguing polling system as a two-
dimensional continuous-time Markovian process, with state space $\\{(n,d)\\}$,
where $n\in\\{0,1,2,3,...\\}$, and

$\displaystyle d\in\mathfrak{D}=\\{$
$\displaystyle(1;1,1),(1;0,1),(1;1,0),(1;0,0),(2;0,1),(2;0,0),(2;-1,0),(2;-1,-1),$
$\displaystyle(3;1,0),(3;0,0),(3;0,-1),(3;-1,-1)\\}.$

In steady state, the joint probability distribution function is denoted by
$P_{n,d}=\mathbb{P}(L_{1}=n,D=d)$.

##  3 Steady-state analysis using probability generating functions

In this section we construct the steady-state joint probability distribution
function of the process $\\{(n,d)\\}$. We use the above mentioned non-
conventional construction (see Perel et al., (2020)) to obtain $12$ partial
probability generating functions (PGFs) defining the system’s states and
utilize their properties, as described below.

###  3.1 Balance equations and PGFs

Writing the balance equations for all $n$, we obtain:
For $d=(1;1,1)$,

$\displaystyle(\lambda+\mu_{1})P_{n,(1;1,1)}=$ $\displaystyle\lambda
p_{1}(P_{n-1,(1;0,0)}+P_{n-1,(2;0,0)}+P_{n-1,(3;0,0)})+\mu_{2}P_{n,(2;0,1)}+\mu_{3}P_{n,(3;1,0)},$
$\displaystyle n\geq 1.$ (3.1)

For $d=(1;0,1)$,

$\displaystyle(\lambda+\mu_{1})P_{n,(1;0,1)}=\lambda\frac{p_{2}}{p_{2}+p_{3}}P_{n,(1;1,1)}+\mu_{3}\frac{q_{1}}{q_{1}+q_{2}}P_{n,(3;0,0)},$
$\displaystyle n\geq 1.$ (3.2)

For $d=(1;1,0)$,

$\displaystyle(\lambda+\mu_{1})P_{n,(1;1,0)}=\lambda\frac{p_{3}}{p_{2}+p_{3}}P_{n,(1;1,1)}+\mu_{2}\frac{q_{1}}{q_{1}+q_{3}}P_{n,(2;0,0)},$
$\displaystyle n\geq 1.$ (3.3)

For $d=(1;0,0)$,

$\displaystyle\lambda P_{0,(1;0,0)}=\mu_{1}P_{1,(1;1,1)},\ $ $\displaystyle
n=0.$ (3.4)
$\displaystyle(\lambda+\mu_{1})P_{n,(1;0,0)}=\lambda\left(P_{n,(1;1,0)}+P_{n,(1;0,1)}\right)+\mu_{1}P_{n+1,(1;1,1)},$
$\displaystyle n\geq 1.$ (3.5)

For $d=(2;0,1)$,

$\displaystyle(\lambda+\mu_{2})P_{n,(2;0,1)}=\lambda\frac{p_{1}}{p_{1}+p_{3}}P_{n-1,(2;-1,0)}+\mu_{3}\frac{q_{2}}{q_{1}+q_{2}}P_{n,(3;0,0)},$
$\displaystyle n\geq 1.$ (3.6)

For $d=(2;0,0)$,

$\displaystyle\lambda P_{0,(2;0,0)}=\mu_{2}P_{0,(2;-1,0)},\ $ $\displaystyle
n=0.$ (3.7)
$\displaystyle(\lambda+\mu_{2})P_{n,(2;0,0)}=\lambda\left(P_{n-1,(2;-1,-1)}+P_{n,(2;0,1)}\right)+\mu_{2}P_{n,(2;-1,0)},$
$\displaystyle n\geq 1.$ (3.8)

For $d=(2;-1,0)$,

$\displaystyle(\lambda+\mu_{2})P_{n,(2;-1,0)}=$ $\displaystyle\lambda
p_{2}(P_{n,(1;0,0)}+P_{n,(2;0,0)}+P_{n,(3;0,0)})+\mu_{1}P_{n+1,(1;0,1)}+\mu_{3}P_{n,(3;-1,-1)},$
$\displaystyle n\geq 0.$ (3.9)

for $d=(2;-1,-1)$,

$\displaystyle(\lambda+\mu_{2})P_{n,(2;-1,-1)}=\lambda\frac{p_{3}}{p_{1}+p_{3}}P_{n,(2;-1,0)}+\mu_{1}\frac{q_{2}}{q_{2}+q_{3}}P_{n+1,(1;0,0)},$
$\displaystyle n\geq 0.$ (3.10)

For $d=(3;1,0)$,

$\displaystyle(\lambda+\mu_{3})P_{n,(3;1,0)}=\lambda\frac{p_{1}}{p_{1}+p_{2}}P_{n-1,(3;0,-1)}+\mu_{2}\frac{q_{3}}{q_{1}+q_{3}}P_{n,(2;0,0)},$
$\displaystyle n\geq 1.$ (3.11)

For $d=(3;0,0)$,

$\displaystyle\lambda P_{0,(3;0,0)}=\mu_{3}P_{0,(3;0,-1)},$ $\displaystyle
n=0.$ (3.12)
$\displaystyle(\lambda+\mu_{3})P_{n,(3;0,0)}=\lambda\left(P_{n-1,(3;-1,-1)}+P_{n,(3;1,0)}\right)+\mu_{3}P_{n,(3;0,-1)},$
$\displaystyle n\geq 1.$ (3.13)

For $d=(3;0,-1)$,

$\displaystyle(\lambda+\mu_{3})P_{n,(3;0,-1)}=\lambda
p_{3}(P_{n,(1;0,0)}+P_{n,(2;0,0)}+P_{n,(3;0,0)})+\mu_{1}P_{n+1,(1;1,0)}+\mu_{2}P_{n,(2;-1,-1)},$
$\displaystyle n\geq 0.$ (3.14)

Finally, for $d=(3;-1,-1)$,

$\displaystyle(\lambda+\mu_{3})P_{n,(3;-1,-1)}=\lambda\frac{p_{2}}{p_{1}+p_{2}}P_{n,(3;0,-1)}+\mu_{1}\frac{q_{3}}{q_{2}+q_{3}}P_{n+1,(1;0,0)},$
$\displaystyle n\geq 0.$ (3.15)

Note that in equations (3.1)–(3.15) above, equations (3.4), (3.7) and (3.12)
hold for $n=0$; equations (3.9), (3.10), (3.14) and (3.15) hold for $n\geq 0$,
while the rest of the equations hold for $n\geq 1$.

For each $d\in\mathfrak{D}$, define the conditional probability generating
function of the number of customers in $Q_{1}$ as:

$G_{d}(z)=\sum_{n}{P_{n,d}z^{n}},\ \ \ d\in\mathfrak{D}.$

Multiplying equation (3.1) by $z^{n}$ and summing over $n\geq 1$ results in

$\displaystyle(\lambda+\mu_{1})G_{(1;1,1)}(z)=$ $\displaystyle\lambda
p_{1}zG_{(1;0,0)}(z)+\lambda p_{1}zG_{(2;0,0)}(z)+\lambda
p_{1}zG_{(3;0,0)}(z)$
$\displaystyle+\mu_{2}G_{(2;0,1)}(z)+\mu_{3}G_{(3;1,0)}(z).$ (3.16)

Repeating this process for all $d\in\mathfrak{D}$, while using equations
(3.2)–(3.15), leads to

$\displaystyle(\lambda+\mu_{1})G_{(1;0,1)}(z)=\lambda\frac{p_{2}}{p_{2}+p_{3}}G_{(1;1,1)}(z)+\mu_{3}\frac{q_{1}}{q_{1}+q_{2}}G_{(3;0,0)}(z)-\mu_{3}\frac{q_{1}}{q_{1}+q_{2}}P_{0,(3;0,0)},$
(3.17)
$\displaystyle(\lambda+\mu_{1})G_{(1;1,0)}(z)=\lambda\frac{p_{3}}{p_{2}+p_{3}}G_{(1;1,1)}(z)+\mu_{2}\frac{q_{1}}{q_{1}+q_{3}}G_{(2;0,0)}(z)-\mu_{2}\frac{q_{1}}{q_{1}+q_{3}}P_{0,(2;0,0)},$
(3.18) $\displaystyle(\lambda+\mu_{1})zG_{(1;0,0)}(z)=\lambda
zG_{(1;1,0)}(z)+\lambda
zG_{(1;0,1)}(z)+\mu_{1}G_{(1;1,1)}(z)+\mu_{1}zP_{0,(1;0,0)},$ (3.19)
$\displaystyle(\lambda+\mu_{2})G_{(2;0,1)}(z)=\lambda
z\frac{p_{1}}{p_{1}+p_{3}}G_{(2;-1,0)}(z)+\mu_{3}\frac{q_{2}}{q_{1}+q_{2}}G_{(3;0,0)}(z)-\mu_{3}\frac{q_{2}}{q_{1}+q_{2}}P_{0,(3;0,0)},$
(3.20) $\displaystyle(\lambda+\mu_{2})G_{(2;0,0)}(z)=\lambda
zG_{(2;-1,-1)}(z)+\lambda
G_{(2;0,1)}(z)+\mu_{2}G_{(2;-1,0)}(z)+\mu_{2}P_{0,(2;0,0)},$ (3.21)
$\displaystyle(\lambda+\mu_{2})zG_{(2;-1,0)}(z)=\lambda
p_{2}zG_{(1;0,0)}(z)+\lambda p_{2}zG_{(2;0,0)}(z)+\lambda
p_{2}zG_{(3;0,0)}(z)$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ +\mu_{1}G_{(1;0,1)}(z)+\mu_{3}zG_{(3;-1,-1)}(z),$ (3.22)
$\displaystyle(\lambda+\mu_{2})zG_{(2;-1,-1)}(z)=\lambda
z\frac{p_{3}}{p_{1}+p_{3}}G_{(2;-1,0)}(z)+\mu_{1}\frac{q_{2}}{q_{2}+q_{3}}G_{(1;0,0)}(z)-\mu_{1}\frac{q_{2}}{q_{2}+q_{3}}P_{0,(1;0,0)},$
(3.23) $\displaystyle(\lambda+\mu_{3})G_{(3;1,0)}(z)=\lambda
z\frac{p_{1}}{p_{1}+p_{2}}G_{(3;0,-1)}(z)+\mu_{2}\frac{q_{3}}{q_{1}+q_{3}}G_{(2;0,0)}(z)-\mu_{2}\frac{q_{3}}{q_{1}+q_{3}}P_{0,(2;0,0)},$
(3.24) $\displaystyle(\lambda+\mu_{3})G_{(3;0,0)}(z)=\lambda
zG_{(3;-1,-1)}(z)+\lambda
G_{(3;1,0)}(z)+\mu_{3}G_{(3;0,-1)}(z)+\mu_{3}P_{0,(3;0,0)},$ (3.25)
$\displaystyle(\lambda+\mu_{3})zG_{(3;0,-1)}(z)=\lambda
p_{3}zG_{(1;0,0)}(z)+\lambda p_{3}zG_{(2;0,0)}(z)+\lambda
p_{3}zG_{(3;0,0)}(z)$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ +\mu_{1}G_{(1;1,0)}(z)+\mu_{2}zG_{(2;-1,-1)}(z),$ (3.26)
$\displaystyle(\lambda+\mu_{3})zG_{(3;-1,-1)}(z)=\lambda
z\frac{p_{2}}{p_{1}+p_{2}}G_{(3;0,-1)}(z)+\mu_{1}\frac{q_{3}}{q_{2}+q_{3}}G_{(1;0,0)}(z)-\mu_{1}\frac{q_{3}}{q_{2}+q_{3}}P_{0,(1;0,0)},$
(3.27)

The set of $12$ equations (3.1)-(3.27) can be written in a matrix form as

$A(z)\cdot\vec{G}(z)=\vec{P}(z),$ (3.28)

where

$A(z)=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}\lambda+\mu_{1}&0&0&-\lambda
p_{1}z&-\mu_{2}&-\lambda p_{1}z&0&0&-\mu_{3}&-\lambda
p_{1}z&0&0\\\\[6.40002pt] \frac{-\lambda
p_{2}}{p_{2}+p_{3}}&\lambda+\mu_{1}&0&0&0&0&0&0&0&\frac{-\mu_{3}q_{1}}{q_{1}+q_{2}}&0&0\\\\[6.40002pt]
\frac{-\lambda
p_{3}}{p_{2}+p_{3}}&0&\lambda+\mu_{1}&0&0&\frac{-\mu_{2}q_{1}}{q_{1}+q_{3}}&0&0&0&0&0&0\\\\[6.40002pt]
-\mu_{1}&-\lambda z&-\lambda
z&(\lambda+\mu_{1})z&0&0&0&0&0&0&0&0\\\\[6.40002pt]
0&0&0&0&\lambda+\mu_{2}&0&\frac{-\lambda
zp_{1}}{p_{1}+p_{3}}&0&0&\frac{-\mu_{3}q_{2}}{q_{1}+q_{2}}&0&0\\\\[6.40002pt]
0&0&0&0&-\lambda&\lambda+\mu_{2}&-\mu_{2}&-\lambda z&0&0&0&0\\\\[6.40002pt]
0&-\mu_{1}&0&-\lambda p_{2}z&0&-\lambda p_{2}z&(\lambda+\mu_{2})z&0&0&-\lambda
p_{2}z&0&-\mu_{3}z\\\\[6.40002pt]
0&0&0&\frac{-\mu_{1}q_{2}}{q_{2}+q_{3}}&0&0&\frac{-\lambda
zp_{3}}{p_{1}+p_{3}}&(\lambda+\mu_{2})z&0&0&0&0\\\\[6.40002pt]
0&0&0&0&0&\frac{-\mu_{2}q_{3}}{q_{1}+q_{3}}&0&0&\lambda+\mu_{3}&0&\frac{-\lambda
zp_{1}}{p_{1}+p_{2}}&0\\\\[6.40002pt]
0&0&0&0&0&0&0&0&-\lambda&\lambda+\mu_{3}&-\mu_{3}&-\lambda z\\\\[6.40002pt]
0&0&-\mu_{1}&-\lambda p_{3}z&0&-\lambda p_{3}z&0&-\mu_{2}z&0&-\lambda
p_{3}z&(\lambda+\mu_{3})z&0\\\\[6.40002pt]
0&0&0&\frac{-\mu_{1}q_{3}}{q_{2}+q_{3}}&0&0&0&0&0&0&\frac{-\lambda
zp_{2}}{p_{1}+p_{2}}&(\lambda+\mu_{3})z\end{pmatrix},$
$\displaystyle\vec{G}(z)=\Big{(}G$
${}_{(1;1,1)}(z),G_{(1;0,1)}(z),G_{(1;1,0)}(z),G_{(1;0,0)}(z),G_{(2;0,1)}(z),G_{(2;0,0)}(z),G_{(2;-1,0)}(z),$
$\displaystyle
G_{(2;-1,-1)}(z),G_{(3;1,0)}(z),G_{(3;0,0)}(z),G_{(3;0,-1)}(z),G_{(3;-1,-1)}(z)\Big{)}^{T}$

is a 12-dimensional column vector of the desired conditional PGFs, and

$\displaystyle\vec{P}(z)=\Big{(}0$
$\displaystyle,\frac{-\mu_{3}q_{1}}{q_{1}+q_{2}}P_{0,(3;0,0)},\frac{-\mu_{2}q_{1}}{q_{1}+q_{3}}P_{0,(2;0,0)},\mu_{1}zP_{0,(1;0,0)},\frac{-\mu_{3}q_{2}}{q_{1}+q_{2}}P_{0,(3;0,0)},\mu_{2}P_{0,(2;0,0)},$
$\displaystyle
0,\frac{-\mu_{1}q_{2}}{q_{2}+q_{3}}P_{0,(1;0,0)},\frac{-\mu_{2}q_{3}}{q_{1}+q_{3}}P_{0,(2;0,0)},\mu_{3}P_{0,(3;0,0)},0,\frac{-\mu_{1}q_{3}}{q_{2}+q_{3}}P_{0,(1;0,0)}\Big{)}^{T}$

is a 12-dimensional vector containing the three unknowns, so-called ’boundary
probabilities’, $P_{0,(1;0,0)}$, $P_{0,(2;0,0)}$ and $P_{0,(3;0,0)}$. Using
Cramer’s rule we have $G_{d}(z)=\frac{\mid A_{d}(z)\mid}{\mid A(z)\mid}$,
$d\in\mathfrak{D}$, where $|A|$ is the determinant of a matrix $A$, and
$A_{d}(z)$ is the matrix obtained from $A(z)$ by replacing the corresponding
column of the latter matrix with $\vec{P}(z)$. Each of the PGFs $G_{d}(z)$,
$d\in\mathfrak{D}$, is a function of the three unknown probabilities,
$P_{0,(1;0,0)}$, $P_{0,(2;0,0)}$ and $P_{0,(3;0,0)}$ appearing in
$\vec{P}(z)$. In order to calculate those probabilities, we use first the
normalization condition (3.29),

$\sum_{d\in\mathfrak{D}}G_{d}(1)=\sum_{d\in\mathfrak{D}}\lim\limits_{z\rightarrow
1}\frac{\mid A_{d}(z)\mid}{\mid A(z)\mid}=1,$ (3.29)

resulting in one equation in the above boundary probabilities, while two other
equations between the three boundary probabilities are derived as follows.
Since $G_{d}(z)$ is defined for all $|z|\leq 1$, each root of $|A(z)|$ is a
root of $|A_{d}(z)|$. As a result, we utilize the two roots of $|A(z)|$ lying
in $(-1,1)$, and solve $|A_{d}(z)|=0$, for arbitrary two different values of
$d$. Once the boundary probabilities are calculated, the 12 PGFs can be
obtained.

###  3.2 Performance measures

In this section we calculate the mean and variance of $L_{i}$, $i=1,2,3$, as
well as the covariances and correlation coefficients between each pair of
queue lengths. We also derive the effective arrival rate to each queue, the
probability that the server is idle, and the mean waiting time in each queue.

First

$\displaystyle\mathbb{E}[L_{1}]=\sum\limits_{d\in\mathfrak{D}}{G^{\prime}_{d}(1)},$
$\displaystyle\mathbb{E}[D_{2}]=\mathbb{E}[L_{1}-L_{2}]=G_{(1;1,1)}(1)+G_{(1;1,0)}(1)+G_{(3;1,0)}(1)-G_{(2;-1,0)}(1)-G_{(2;-1,-1)}(1)-G_{(3;-1,-1)}(1),$
$\displaystyle\mathbb{E}[D_{3}]=\mathbb{E}[L_{1}-L_{3}]=G_{(1;1,1)}(1)+G_{(1;0,1)}(1)+G_{(2;0,1)}(1)-G_{(2;-1,-1)}(1)-G_{(3;0,-1)}(1)-G_{(3;-1,-1)}(1),$
$\displaystyle\mathbb{E}[L_{i}]=\mathbb{E}[L_{1}]-\mathbb{E}[D_{i}],\ \
i=2,3.$

The above follows since the difference $L_{1}-L_{2}$ or $L_{1}-L_{3}$ can
assume only the values -1, 0, or 1. Next, for $j=2,3$,

$\displaystyle Cov(L_{1},L_{j})$
$\displaystyle=\mathbb{E}[L_{1}L_{j}]-\mathbb{E}[L_{1}]\mathbb{E}[L_{j}]=\mathbb{E}[L_{1}(L_{1}-D_{j})]-\mathbb{E}[L_{1}]\mathbb{E}[L_{j}]$
$\displaystyle=\mathbb{E}[L_{1}^{2}]-\mathbb{E}[L_{1}D_{j}]-\mathbb{E}[L_{1}]\mathbb{E}[L_{j}],$

where

$\displaystyle\mathbb{E}[L_{1}^{2}]=\sum\limits_{d\in\mathfrak{D}}{G^{\prime\prime}_{d}(1)}+\mathbb{E}[L_{1}],$
$\displaystyle\mathbb{E}[L_{1}D_{2}]=G^{\prime}_{(1;1,1)}(1)+G^{\prime}_{(1;1,0)}(1)+G^{\prime}_{(3;1,0)}(1)-G^{\prime}_{(2;-1,0)}(1)-G^{\prime}_{(2;-1,-1)}(1)-G^{\prime}_{(3;-1,-1)}(1),$
$\displaystyle\mathbb{E}[L_{1}D_{3}]=G^{\prime}_{(1;1,1)}(1)+G^{\prime}_{(1;0,1)}(1)+G^{\prime}_{(2;0,1)}(1)-G^{\prime}_{(2;-1,-1)}(1)-G^{\prime}_{(3;0,-1)}(1)-G^{\prime}_{(3;-1,-1)}(1).$

Also,

$\displaystyle Cov(L_{2},L_{3})$
$\displaystyle=\mathbb{E}[L_{2}L_{3}]-\mathbb{E}[L_{2}]\mathbb{E}[L_{3}]=\mathbb{E}[(L_{1}-D_{2})(L_{1}-D_{3})]-\mathbb{E}[L_{2}]\mathbb{E}[L_{3}]$
$\displaystyle=\mathbb{E}[L_{1}^{2}]-\mathbb{E}[L_{1}D_{2}]-\mathbb{E}[L_{1}D_{3}]+\mathbb{E}[D_{2}D_{3}]-\mathbb{E}[L_{2}]\mathbb{E}[L_{3}],$

where

$\mathbb{E}[D_{2}D_{3}]=G_{(1;1,1)}(1)+G_{(2;-1,-1)}(1)+G_{(3;-1,-1)}(1).$

Furthermore, the variance of $L_{i}$, for $i=1,2,3$, is given by

$\displaystyle
Var(L_{1})=\mathbb{E}[L_{1}^{2}]-\left(\mathbb{E}[L_{1}]\right)^{2},$
$\displaystyle
Var(L_{j})=\mathbb{E}[L_{j}^{2}]-\left(\mathbb{E}[L_{j}]\right)^{2}=\mathbb{E}[(L_{1}-D_{j})^{2}]-\left(\mathbb{E}[L_{j}]\right)^{2}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \
=\mathbb{E}[L_{1}^{2}]-2\mathbb{E}[L_{1}D_{j}]+\mathbb{E}[D_{j}^{2}]-\left(\mathbb{E}[L_{j}]\right)^{2},$

where

$\displaystyle\mathbb{E}[D_{2}^{2}]=G_{(1;1,1)}(1)+G_{(1;1,0)}(1)+G_{(2;-1,0)}(1)-G_{(2;-1,-1)}(1)-G_{(3;1,0)}(1)-G_{(3;-1,-1)}(1),$
$\displaystyle\mathbb{E}[D_{3}^{2}]=G_{(1;1,1)}(1)+G_{(1;0,1)}(1)+G_{(2;0,1)}(1)-G_{(2;-1,-1)}(1)-G_{(3;0,-1)}(1)-G_{(3;-1,-1)}(1).$

From all the above, the correlation coefficient between $L_{i}$ and $L_{j}$,
denoted by $Cor(L_{i},L_{j})$, can be explicitly calculated, using
$Cor(L_{i},L_{j})=\frac{Cov(L_{i},L_{j})}{\sqrt{Var(L_{i})Var(L_{j})}}$.

Let $\lambda^{i}_{eff}$ denote the effective arrival rate to $Q_{i}$,
$i=1,2,3$. Then,

$\displaystyle\lambda^{1}_{eff}$
$\displaystyle=\lambda\Big{(}p_{1}(G_{(1;0,0)}(1)+G_{(2;0,0)}(1)+G_{(3;0,0)}(1))+G_{(2;-1,-1)}(1)+G_{(3;-1,-1)}(1)+$
$\displaystyle\ \ \ \ \ \ \ \
\frac{p_{1}}{p_{1}+p_{2}}G_{(3;0,-1)}(1)+\frac{p_{1}}{p_{1}+p_{3}}G_{(2;-1,0)}(1)\Big{)},$
$\displaystyle\lambda^{2}_{eff}$
$\displaystyle=\lambda\Big{(}p_{2}(G_{(1;0,0)}(1)+G_{(2;0,0)}(1)+G_{(3;0,0)}(1))+G_{(1;1,0)}(1)+G_{(3;1,0)}(1)+$
$\displaystyle\ \ \ \ \ \ \ \
\frac{p_{2}}{p_{1}+p_{2}}G_{(3;0,-1)}(1)+\frac{p_{2}}{p_{2}+p_{3}}G_{(1;1,1)}(1)\Big{)},$
$\displaystyle\lambda^{3}_{eff}$
$\displaystyle=\lambda\Big{(}p_{3}(G_{(1;0,0)}(1)+G_{(2;0,0)}(1)+G_{(3;0,0)}(1))+G_{(1;0,1)}(1)+G_{(2;0,1)}(1)+$
$\displaystyle\ \ \ \ \ \ \ \
\frac{p_{3}}{p_{1}+p_{3}}G_{(2;-1,0)}(1)+\frac{p_{3}}{p_{2}+p_{3}}G_{(1;1,1)}(1)\Big{)}.$

Clearly, $\sum_{i=1}^{3}{\lambda^{i}_{eff}}=\lambda$. Defining
$\rho_{eff}^{i}=\frac{\lambda^{i}_{eff}}{\mu_{i}}$, we have that

$\mathbb{P}(\text{Server is
idle})=P_{0,(1;0,0)}+P_{0,(2;0,0)}+P_{0,(3;0,0)}=1-\sum_{i=1}^{3}{\rho_{eff}^{i}}.$

Furthermore, define $W_{i}$ as the sojourn time of a customer in $Q_{i}$, for
$i=1,2,3$. Then, by Little’s Law,

$\mathbb{E}[W_{i}]=\frac{\mathbb{E}[L_{i}]}{\lambda^{i}_{eff}}\ .$

Last, denote by $\gamma_{i}$ the probability that the sever resides in
$Q_{i}$, for $i=1,2,3$. Then,

$\displaystyle\gamma_{1}=G_{(1;1,1)}(1)+G_{(1;0,1)}(1)+G_{(1;1,0)}(1)+G_{(1;0,0)}(1),$
$\displaystyle\gamma_{2}=G_{(2;0,1)}(1)+G_{(2;0,0)}(1)+G_{(2;-1,0)}(1)+G_{(2;-1,-1)}(1),$
$\displaystyle\gamma_{3}=G_{(3;1,0)}(1)+G_{(3;0,0)}(1)+G_{(3;0,-1)}(1)+G_{(3;-1,-1)}(1).$

##  4 Matrix Geometric Analysis

###  4.1 Definitions and notations

In this section we present an alternative approach to analyze the 3-queue JSQ-
SLQ model by defining a Quasi Birth and Death (QBD) process, possessing 12
phases and infinite number of levels. With this analysis we derive the
stability condition of the system. Phases are represented by the states $d$,
$d\in\mathfrak{D}$, while levels are presented by the different values of
$L_{1}\geq 0$. For $n\geq 0$ define $\mathcal{S}_{n}$ to be the set of
possible states for level $n$, where

$\displaystyle\mathcal{S}_{0}=$
$\displaystyle\\{(0,(1;0,0)),(0,(2;0,0)),(0,(2;-1,0)),(0,(2;-1,-1)),(0,(3;0,0)),(0,(3;0,-1)),(0,(3;-1,-1))\\},$
$\displaystyle\mathcal{S}_{n}=$ $\displaystyle\\{(n,d)\mid
d\in\mathfrak{D}\\},\ \ n\geq 1,$

and arrange the system’s states in the order

$\displaystyle\mathcal{S}=\Big{\\{}\mathcal{S}_{0};\mathcal{S}_{1};\mathcal{S}_{2};\ldots;\mathcal{S}_{n}\ldots\Big{\\}}.$

The infinitesimal generator matrix of the QBD process, denoted by $Q$, is
given by

$Q=\begin{pmatrix}B_{1}&B_{0}&0&\cdots&\cdots&\cdots\\\\[8.00003pt]
B_{2}&A_{1}&A_{0}&0&\cdots&\cdots\\\\[8.00003pt]
0&A_{2}&A_{1}&A_{0}&0&\cdots\\\\[8.00003pt]
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots\end{pmatrix},$

where,

$B_{0}=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}\lambda
p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] \lambda
p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{3}}&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&\lambda
p_{1}&0&0&0&0&0&0\\\\[3.00003pt] \lambda
p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{2}}&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&\lambda&0&0\end{pmatrix}_{7\times 12},$

$B_{1}=\begin{pmatrix}-\lambda&0&\lambda p_{2}&0&0&\lambda
p_{3}&0\\\\[3.00003pt] 0&-\lambda&\lambda p_{2}&0&0&\lambda
p_{3}&0\\\\[3.00003pt] 0&\mu_{2}&-(\lambda+\mu_{2})&\frac{\lambda
p_{3}}{p_{1}+p_{3}}&0&0&0\\\\[3.00003pt]
0&0&0&-(\lambda+\mu_{2})&0&\mu_{2}&0\\\\[3.00003pt] 0&0&\lambda
p_{2}&0&-\lambda&\lambda p_{3}&0\\\\[3.00003pt]
0&0&0&0&\mu_{3}&-(\lambda+\mu_{3})&\frac{\lambda
p_{2}}{p_{1}+p_{2}}\\\\[3.00003pt]
0&0&\mu_{3}&0&0&0&-(\lambda+\mu_{3})\end{pmatrix}_{7\times 7},$

$B_{2}=\begin{pmatrix}\mu_{1}&0&0&0&0&0&0\\\\[3.00003pt]
0&0&\mu_{1}&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&\mu_{1}&0\\\\[3.00003pt]
0&0&0&\frac{\mu_{1}q_{2}}{q_{2}+q_{3}}&0&0&\frac{\mu_{1}q_{3}}{q_{2}+q_{3}}\\\\[3.00003pt]
0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0\end{pmatrix}_{12\times 7}.$

and

$A_{0}=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
\lambda p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] \lambda
p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{3}}&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&\lambda&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] \lambda
p_{1}&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{2}}&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&\lambda&0&0\end{pmatrix}_{12\times 12},$

$A_{1}=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}-(\lambda+\mu_{1})&\frac{\lambda
p_{2}}{p_{2}+p_{3}}&\frac{\lambda
p_{3}}{p_{2}+p_{3}}&0&0&0&0&0&0&0&0&0\\\\[1.50002pt]
0&-(\lambda+\mu_{1})&0&\lambda&0&0&0&0&0&0&0&0\\\\[1.50002pt]
0&0&-(\lambda+\mu_{1})&\lambda&0&0&0&0&0&0&0&0\\\\[1.50002pt]
0&0&0&-(\lambda+\mu_{1})&0&0&\lambda p_{2}&0&0&0&\lambda
p_{3}&0\\\\[1.50002pt]
\mu_{2}&0&0&0&-(\lambda+\mu_{2})&\lambda&0&0&0&0&0&0\\\\[1.50002pt]
0&0&\frac{\mu_{2}q_{1}}{q_{1}+q_{3}}&0&0&-(\lambda+\mu_{2})&\lambda
p_{2}&0&\frac{\mu_{2}q_{3}}{q_{1}+q_{3}}&0&\lambda p_{3}&0\\\\[1.50002pt]
0&0&0&0&0&\mu_{2}&-(\lambda+\mu_{2})&\frac{\lambda
p_{3}}{p_{1}+p_{3}}&0&0&0&0\\\\[1.50002pt]
0&0&0&0&0&0&0&-(\lambda+\mu_{2})&0&0&\mu_{2}&0\\\\[1.50002pt]
\mu_{3}&0&0&0&0&0&0&0&-(\lambda+\mu_{3})&\lambda&0&0\\\\[1.50002pt]
0&\frac{\mu_{3}q_{1}}{q_{1}+q_{2}}&0&0&\frac{\mu_{3}q_{2}}{q_{1}+q_{2}}&0&\lambda
p_{2}&0&0&-(\lambda+\mu_{3})&\lambda p_{3}&0\\\\[1.50002pt]
0&0&0&0&0&0&0&0&0&\mu_{3}&-(\lambda+\mu_{3})&\frac{\lambda
p_{2}}{p_{1}+p_{2}}\\\\[1.50002pt]
0&0&0&0&0&0&\mu_{3}&0&0&0&0&-(\lambda+\mu_{3})\end{pmatrix}_{12\times 12},$

$A_{2}=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}0&0&0&\mu_{1}&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&\mu_{1}&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&\mu_{1}&0\\\\[3.00003pt]
0&0&0&0&0&0&0&\frac{\mu_{1}q_{2}}{q_{2}+q_{3}}&0&0&0&\frac{\mu_{1}q_{3}}{q_{2}+q_{3}}\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt] 0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\\\\[3.00003pt]
0&0&0&0&0&0&0&0&0&0&0&0\end{pmatrix}_{12\times 12}.$

###  4.2 Stability condition

Define the matrix $A=A_{0}+A_{1}+A_{2}$. We get

$A=\setcounter{MaxMatrixCols}{12}\begin{pmatrix}-(\lambda+\mu_{1})&\frac{\lambda
p_{2}}{p_{2}+p_{3}}&\frac{\lambda
p_{3}}{p_{2}+p_{3}}&\mu_{1}&0&0&0&0&0&0&0&0\\\\[1.50002pt]
0&-(\lambda+\mu_{1})&0&\lambda&0&0&\mu_{1}&0&0&0&0&0\\\\[1.50002pt]
0&0&-(\lambda+\mu_{1})&\lambda&0&0&0&0&0&0&\mu_{1}&0\\\\[1.50002pt] \lambda
p_{1}&0&0&-(\lambda+\mu_{1})&0&0&\lambda
p_{2}&\frac{\mu_{1}q_{2}}{q_{2}+q_{3}}&0&0&\lambda
p_{3}&\frac{\mu_{1}q_{3}}{q_{2}+q_{3}}\\\\[1.50002pt]
\mu_{2}&0&0&0&-(\lambda+\mu_{2})&\lambda&0&0&0&0&0&0\\\\[1.50002pt] \lambda
p_{1}&0&\frac{\mu_{2}q_{1}}{q_{1}+q_{3}}&0&0&-(\lambda+\mu_{2})&\lambda
p_{2}&0&\frac{\mu_{2}q_{3}}{q_{1}+q_{3}}&0&\lambda p_{3}&0\\\\[1.50002pt]
0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{3}}&\mu_{2}&-(\lambda+\mu_{2})&\frac{\lambda
p_{3}}{p_{1}+p_{3}}&0&0&0&0\\\\[1.50002pt]
0&0&0&0&0&\lambda&0&-(\lambda+\mu_{2})&0&0&\mu_{2}&0\\\\[1.50002pt]
\mu_{3}&0&0&0&0&0&0&0&-(\lambda+\mu_{3})&\lambda&0&0\\\\[1.50002pt] \lambda
p_{1}&\frac{\mu_{3}q_{1}}{q_{1}+q_{2}}&0&0&\frac{\mu_{3}q_{2}}{q_{1}+q_{2}}&0&\lambda
p_{2}&0&0&-(\lambda+\mu_{3})&\lambda p_{3}&0\\\\[1.50002pt]
0&0&0&0&0&0&0&0&\frac{\lambda
p_{1}}{p_{1}+p_{2}}&\mu_{3}&-(\lambda+\mu_{3})&\frac{\lambda
p_{2}}{p_{1}+p_{2}}\\\\[1.50002pt]
0&0&0&0&0&0&\mu_{3}&0&0&\lambda&0&-(\lambda+\mu_{3})\end{pmatrix}_{12\times
12},$

The matrix $A$ is the infinitesimal generator matrix of the process describing
the evolution of $D$, given that $L_{1}\geq 1$. Let
$\vec{\pi}=\left(\pi_{1},\pi_{2},...,\pi_{12}\right)$ be the stationary vector
of the matrix $A$, i.e. $\vec{\pi}A=\vec{0}$ and $\vec{\pi}\cdot\vec{e}=1$
(where $\vec{e}$ is a 12-dimensional column vector with all its entries equal
to 1). From Neuts, (1981), we have that the stability condition is

$\vec{\pi}A_{0}\vec{e}<\vec{\pi}A_{2}\vec{e}.$

Hanukov and Yechiali, (2021) showed that if the 3 matrices $A_{0}$, $A_{1}$
and $A_{2}$ are all upper triangular, or all lower triangular, then there
exists a simple and direct stability condition. Unfortunately, this is not the
case here, so we have to use the above condition, which translates into

$\lambda\left(\pi_{8}+\pi_{12}+p_{1}(\pi_{4}+\pi_{6}+\pi_{10})+\frac{p_{1}}{p_{1}+p_{2}}\pi_{11}+\frac{p_{1}}{p_{1}+p_{3}}\pi_{7}\right)<\mu_{1}(\pi_{1}+\pi_{2}+\pi_{3}+\pi_{4}).$

Explicit numerical expressions for the elements of $\vec{\pi}$ can be
calculated by any mathematical software. We omit these expressions. However,
as an illustration, we assume values for the parameters $\mu_{i}$, $p_{i}$ and
$q_{i}$ for $i=1,2,3$, and define

$f(\lambda)=\lambda\left(\pi_{8}+\pi_{12}+p_{1}(\pi_{4}+\pi_{6}+\pi_{10})+\frac{p_{1}}{p_{1}+p_{2}}\pi_{11}+\frac{p_{1}}{p_{1}+p_{3}}\pi_{7}\right)-\mu_{1}(\pi_{1}+\pi_{2}+\pi_{3}+\pi_{4}).$

That is, the system is stable iff $f(\lambda)<0$. For example, consider the
symmetric case where $\mu_{i}=5$, $p_{i}=\frac{1}{3}$ and $q_{i}=1$ for
$i=1,2,3$. We then get

$\vec{\pi}=\left(\frac{1}{9},\frac{1}{18},\frac{1}{18},\frac{1}{9},\frac{1}{18},\frac{1}{9},\frac{1}{9},\frac{1}{18},\frac{1}{18},\frac{1}{9},\frac{1}{9},\frac{1}{18}\right),$

while the stability condition $f(\lambda)<0$ translates into
$\frac{\lambda-5}{3}<0$, so the system is stable, as expected, iff
$\lambda<5$. Further examples for the stability condition in terms of the
arrival rate $\lambda$ for various sets of parameters are given in Table 1.

Table 1: Stability condition in terms of $\lambda$. Parameters | Stability condition
---|---
$\mu_{1}=4$, $\mu_{2}=5$, $\mu_{3}=6$, $p_{i}=\frac{1}{3}$ and $q_{i}=1$ for $i=1,2,3$ | $\lambda<4.90215$
$\mu_{1}=4$, $\mu_{2}=5$, $\mu_{3}=6$, $p_{1}=0.2$, $p_{2}=0.3$, $p_{3}=0.5$ and $q_{i}=1$ for $i=1,2,3$ | $\lambda<5.0149$
$\mu_{1}=4$, $\mu_{2}=5$, $\mu_{3}=6$, $p_{1}=0.05$, $p_{2}=0.05$, $p_{3}=0.9$ and $q_{i}=1$ for $i=1,2,3$ | $\lambda<5.17335$
$\mu_{1}=4$, $\mu_{2}=5$, $\mu_{3}=6$, $p_{1}=0.9$, $p_{2}=0.05$, $p_{3}=0.05$ and $q_{i}=1$ for $i=1,2,3$ | $\lambda<4.63355$

###  4.3 Calculation of the system’s states steady state distribution

For $n\geq 0$ define the steady-state probability vector $\vec{P}_{n}$, as
follows:

$\vec{P}_{n}=\left\\{\begin{array}[]{ll}\left(P_{0,(1;0,0)},P_{0,(2;0,0)},P_{0,(2;-1,0)},P_{0,(2;-1,-1)},P_{0,(3;0,0)},P_{0,(3;0,-1)},P_{0,(3;-1,-1)}\right),\
n=0,\\\ \left(P_{n,(1;0,0)},P_{n,(1;0,1)},\ldots,P_{n,(3;-1,-1)}\right),\ \
n\geq 1.\\\ \end{array}\right.$

From Neuts, (1981),

$\vec{P}_{n}=\vec{P}_{1}R^{n-1},\ \ n\geq 1,$

where $R$ is the minimal non-negative solution of the matrix quadratic
equation

$A_{0}+RA_{1}+R^{2}A_{2}=0.$ (4.1)

The matrix $R$ can be calculated by using well-known algorithms, see e.g.
Neuts, (1981), Latouche and Ramaswami, (1999), Artalejo and Gómez-Corral,
(2008), and Harchol-Balter, (2013). The vectors $\vec{P}_{0}$ and
$\vec{P}_{1}$ are the solutions of the following linear system of equations:

$\displaystyle\vec{P}_{0}B_{1}+\vec{P}_{1}B_{2}=\vec{0},$
$\displaystyle\vec{P}_{0}B_{0}+\vec{P}_{1}(A_{1}+RA_{2})=\vec{0},$
$\displaystyle\vec{P}_{0}\vec{e}_{0}+\vec{P}_{1}[\textbf{I}-R]^{-1}\vec{e}=1,$
(4.2)

where $\vec{e}_{0}$ is a $7$-dimensional vector of 1’s. Equation (4.2) is the
normalization equation.

##  5 Numerical Results and the Gini Index

In this section we present numerical results for the performance measures
derived in Section 3.2. We calculate $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$,
$i=1,2,3$, for a wide range of parameters, as well as the correlation
coefficients $Cor(L_{i},L_{j})$ between each pair of queue lengths. We also
calculate $\lambda^{i}_{eff}$, i.e. the effective arrival rate to queue $i$,
as well as $\gamma_{i}$, which is the probability that the server resides in
$Q_{i}$, $i=1,2,3$. Tables 2-13 examine the influence of $\mu_{1}$ on the
various performance measures. Furthermore, we calculate the Gini coefficient,
(see e.g. (Dodge,, 2008, pp. 231–233)), which is a well known measure for
inequality. The Gini coefficient with respect to the differences between the
mean queue sizes in the JSQ-SLQ operating model is calculated as follows:

$GI=\frac{\sum\limits_{i=1}^{3}{\sum\limits_{i=1}^{3}{\big{|}\mathbb{E}[L_{i}]-\mathbb{E}[L_{j}]\big{|}}}}{2\cdot
3\cdot\sum\limits_{i=1}^{3}{\mathbb{E}[L_{i}]}}.$

It turns out, as it is shown in the tables below, that the Gini coefficient is
close to zero in almost all sets of parameters, even when the system is close
to saturation (resulting in large mean queue sizes). This outcome demonstrates
the effectiveness of the combined JSQ-SLQ policy in balancing the queue sizes.

In Tables 2-7 we set $\lambda=4$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$
and $q_{1}=q_{2}=q_{3}=1$, while $\mu_{2}$ varies between the tables (assuming
the values 3,4,5), and $\mu_{1}$ varies within each table (where its minimum
is close to the value ensuring stability). In Tables 8-13 we set $\lambda=4$,
$\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$,
while $\mu_{2}$ again varies between the tables assuming the values 3, 4 and
5. Tables 14-17 examine the impact of $q_{1}$, which is the weight assigned to
$Q_{1}$ from the server’s point of view.

Table 2: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=4.5$ | 48.64 | 48.69 | 48.62 | 35.46 | 39.90 | 34.53 | 0.9999 | 0.9999 | 0.9999 | 0.0003
$\mu_{1}=5$ | 8.31 | 8.37 | 8.31 | 5.96 | 6.91 | 5.96 | 0.9970 | 0.9970 | 0.9970 | 0.0016
$\mu_{1}=10$ | 1.33 | 1.46 | 1.40 | 0.88 | 1.25 | 1.05 | 0.9266 | 0.9260 | 0.9219 | 0.0207
$\mu_{1}=30$ | 0.67 | 0.85 | 0.79 | 0.42 | 0.76 | 0.62 | 0.8438 | 0.8390 | 0.8119 | 0.0519
Table 3: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=4.5$ | 1.37 | 1.22 | 1.41 | 0.307 | 0.409 | 0.284
$\mu_{1}=5$ | 1.395 | 1.21 | 1.395 | 0.293 | 0.414 | 0.293
$\mu_{1}=10$ | 1.51 | 1.16 | 1.33 | 0.232 | 0.438 | 0.33
$\mu_{1}=30$ | 1.60 | 1.12 | 1.28 | 0.204 | 0.446 | 0.350
Table 4: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=4$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=3.3$ | 88.46 | 88.44 | 88.41 | 70.28 | 66.44 | 62.69 | 0.9999 | 0.9999 | 0.9999 | 0.0001
$\mu_{1}=3.4$ | 22.85 | 22.83 | 22.80 | 18.04 | 17.20 | 16.22 | 0.9996 | 0.9996 | 0.9996 | 0.0005
$\mu_{1}=3.5$ | 13.26 | 13.25 | 13.22 | 10.41 | 10.01 | 9.43 | 0.9988 | 0.9988 | 0.9988 | 0.0007
$\mu_{1}=4$ | 4.49 | 4.49 | 4.47 | 3.44 | 3.44 | 3.23 | 0.9902 | 0.9902 | 0.9902 | 0.0010
$\mu_{1}=4.5$ | 2.82 | 2.84 | 2.81 | 2.11 | 2.19 | 2.05 | 0.9768 | 0.9769 | 0.9769 | 0.0024
$\mu_{1}=5$ | 2.11 | 2.14 | 2.11 | 1.56 | 1.67 | 1.56 | 0.9616 | 0.9618 | 0.9616 | 0.0031
$\mu_{1}=10$ | 0.80 | 0.89 | 0.86 | 0.54 | 0.72 | 0.66 | 0.8434 | 0.8434 | 0.8341 | 0.0235
$\mu_{1}=30$ | 0.45 | 0.59 | 0.56 | 0.29 | 0.50 | 0.45 | 0.7542 | 0.75154 | 0.7069 | 0.0583
Table 5: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=4$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=3.3$ | 1.26 | 1.33 | 1.41 | 0.383 | 0.334 | 0.283
$\mu_{1}=3.4$ | 1.27 | 1.33 | 1.40 | 0.377 | 0.337 | 0.286
$\mu_{1}=3.5$ | 1.27 | 1.32 | 1.41 | 0.372 | 0.339 | 0.289
$\mu_{1}=4$ | 1.31 | 1.31 | 1.38 | 0.349 | 0.349 | 0.302
$\mu_{1}=4.5$ | 1.34 | 1.29 | 1.37 | 0.332 | 0.357 | 0.311
$\mu_{1}=5$ | 1.36 | 1.28 | 1.36 | 0.318 | 0.364 | 0.318
$\mu_{1}=10$ | 1.47 | 1.23 | 1.30 | 0.262 | 0.389 | 0.349
$\mu_{1}=30$ | 1.56 | 1.19 | 1.25 | 0.237 | 0.399 | 0.364
Table 6: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=5$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=2.75$ | 55.11 | 55.03 | 55.03 | 46.49 | 39.10 | 39.10 | 0.9999 | 0.9999 | 0.9999 | 0.0003
$\mu_{1}=3$ | 8.37 | 8.31 | 8.31 | 6.91 | 5.95 | 5.95 | 0.9970 | 0.9970 | 0.9970 | 0.0016
$\mu_{1}=3.4$ | 3.73 | 3.68 | 3.68 | 2.993 | 2.67 | 2.67 | 0.9859 | 0.9859 | 0.9859 | 0.0030
$\mu_{1}=3.5$ | 3.29 | 3.25 | 3.25 | 2.63 | 2.37 | 2.37 | 0.9823 | 0.9823 | 0.9824 | 0.0027
$\mu_{1}=4$ | 2.14 | 2.11 | 2.11 | 1.67 | 1.56 | 1.56 | 0.9616 | 0.9616 | 0.9618 | 0.0031
$\mu_{1}=4.5$ | 1.63 | 1.61 | 1.61 | 1.24 | 1.2 | 1.2 | 0.9389 | 0.9389 | 0.9391 | 0.0027
$\mu_{1}=5$ | 1.33 | 1.33 | 1.33 | 1 | 1 | 1 | 0.9164 | 0.9164 | 0.9164 | 0
$\mu_{1}=10$ | 0.60 | 0.66 | 0.66 | 0.42 | 0.52 | 0.52 | 0.7744 | 0.7744 | 0.7624 | 0.0208
$\mu_{1}=30$ | 0.35 | 0.46 | 0.46 | 0.22 | 0.38 | 0.38 | 0.6870 | 0.6870 | 0.6316 | 0.0577
Table 7: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=5$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=2.75$ | 1.18 | 1.41 | 1.41 | 0.432 | 0.284 | 0.284
$\mu_{1}=3$ | 1.21 | 1.39 | 1.39 | 0.414 | 0.293 | 0.293
$\mu_{1}=3.4$ | 1.24 | 1.38 | 1.38 | 0.390 | 0.305 | 0.305
$\mu_{1}=3.5$ | 1.26 | 1.37 | 1.37 | 0.386 | 0.307 | 0.307
$\mu_{1}=4$ | 1.28 | 1.36 | 1.36 | 0.364 | 0.318 | 0.318
$\mu_{1}=4.5$ | 1.31 | 1.345 | 1.345 | 0.347 | 0.3265 | 0.3265
$\mu_{1}=5$ | $\frac{4}{3}$ | $\frac{4}{3}$ | $\frac{4}{3}$ | $\frac{1}{3}$ | $\frac{1}{3}$ | $\frac{1}{3}$
$\mu_{1}=10$ | 1.44 | 1.28 | 1.28 | 0.280 | 0.360 | 0.360
$\mu_{1}=30$ | 1.53 | 1.235 | 1.235 | 0.257 | 0.3715 | 0.3715

The following insights arise from the tables:

  1. 1.

In all tables, $\mathbb{E}[L_{i}]$ and $\mathbb{E}[W_{i}]$ ($i=1,2,3$) are all
decreasing functions of $\mu_{1}$.

  2. 2.

The correlation coefficient between $L_{i}$ and $L_{j}$ is always positive,
since customers always join the shortest queue. Furthermore, as $\mu_{1}$
increases, the correlation coefficient between $L_{i}$ and $L_{j}$ decreases.
This phenomenon occurs since for small values of $\mu_{1}$, each service in
$Q_{1}$ lasts for a long time and therefore, as a result of the JSQ policy,
the number of customers in all queues increases simultaneously. On the other
hand, for large values of $\mu_{1}$, the behavior of $L_{j}$ (for $j=2,3$) is
less affected by $L_{1}$, since the server resides in $Q_{1}$ a short amount
of time.

  3. 3.

The Gini coefficient increases monotonically when $\mu_{1}$ increases, but
remains very small (and close to zero is some cases). It is higher when
$\mu_{1}$ is large. This follows since, even if the mean queue sizes are
small, the relative differences between them are more significant.

  4. 4.

When $p_{1}$ is very close to 1, the expected sojourn time in $Q_{1}$ is
smaller than the expected sojourn times in $Q_{2}$ and $Q_{3}$. This occurs
since when $p_{1}$ approaches 1, $\lambda^{1}_{eff}$ increases, which, due to
the SLQ regime, increases the probability that the server resides in $Q_{1}$.
Clearly, the presence of the server in $Q_{1}$ reduces the waiting time there.

  5. 5.

When the value of $p_{1}$ approaches 1, the minimum value of $\mu_{1}$ that
ensures the system’s stability rises. For example, in Table 6, the minimal
value of $\mu_{1}$ for stability is 2.75, while in Table 12 it is 3.4.

  6. 6.

In all tables, as $\mu_{1}$ increases, the probability that the server resides
in $Q_{1}$ decreases.

  7. 7.

Tables 14 and 16 show that The correlation coefficient between $L_{i}$ and
$L_{j}$ remains unchanged as $q_{1}$ increases, while $\mathbb{E}[L_{i}]$ and
$\mathbb{E}[W_{i}]$ decrease. Furthermore, Tables 15 and 17 show that as
$q_{1}$ increases, the probability that the server resides in $Q_{1}$ slightly
increases, while the probabilities that the server resides in $Q_{2}$ and
$Q_{3}$ moderately decrease.

  8. 8.

The numerical results strongly demonstrate that the aim of the combined policy
JSQ+SLQ is very well achieved: in all tables, the mean queue sizes are almost
equal, as well as the mean sojourn times.

Table 8: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=4.2$ | 105.10 | 104.95 | 104.90 | 47.79 | 123.78 | 110.07 | 0.9999 | 0.9999 | 0.9999 | 0.0004
$\mu_{1}=4.5$ | 8.19 | 8.04 | 8.00 | 3.65 | 9.74 | 8.62 | 0.9967 | 0.9967 | 0.9974 | 0.0052
$\mu_{1}=5$ | 3.36 | 3.22 | 3.17 | 1.45 | 4.06 | 3.56 | 0.9814 | 0.9819 | 0.9858 | 0.0130
$\mu_{1}=10$ | 0.63 | 0.52 | 0.48 | 0.23 | 0.90 | 0.75 | 0.7225 | 0.7250 | 0.7910 | 0.0614
$\mu_{1}=30$ | 0.20 | 0.14 | 0.12 | 0.06 | 0.51 | 0.39 | 0.4698 | 0.4558 | 0.5629 | 0.1159
Table 9: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=4.2$ | 2.20 | 0.85 | 0.95 | 0.526 | 0.283 | 0.191
$\mu_{1}=4.5$ | 2.24 | 0.83 | 0.93 | 0.526 | 0.283 | 0.191
$\mu_{1}=5$ | 2.32 | 0.79 | 0.89 | 0.528 | 0.278 | 0.194
$\mu_{1}=10$ | 2.78 | 0.58 | 0.64 | 0.602 | 0.227 | 0.171
$\mu_{1}=30$ | 3.41 | 0.28 | 0.31 | 0.790 | 0.117 | 0.093
Table 10: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=4$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=3.7$ | 59.77 | 59.57 | 59.55 | 28.32 | 64.51 | 61.63 | 0.9999 | 0.9999 | 0.9999 | 0.0008
$\mu_{1}=4$ | 6.89 | 6.69 | 6.67 | 3.18 | 7.46 | 7.11 | 0.9954 | 0.9954 | 0.9963 | 0.0072
$\mu_{1}=4.5$ | 3.42 | 3.20 | 3.20 | 1.58 | 3.48 | 3.48 | 0.9824 | 0.9824 | 0.9860 | 0.0149
$\mu_{1}=5$ | 2.91 | 2.72 | 2.70 | 1.29 | 3.17 | 3.02 | 0.9757 | 0.9760 | 0.9811 | 0.0168
$\mu_{1}=10$ | 0.53 | 0.39 | 0.38 | 0.19 | 0.67 | 0.62 | 0.6530 | 0.6547 | 0.7349 | 0.0769
$\mu_{1}=30$ | 0.18 | 0.11 | 0.10 | 0.05 | 0.39 | 0.35 | 0.4146 | 0.4083 | 0.5083 | 0.1367
Table 11: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=4$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=3.7$ | 2.11 | 0.92 | 0.97 | 0.574 | 0.231 | 0.195
$\mu_{1}=4$ | 2.17 | 0.90 | 0.93 | 0.571 | 0.232 | 0.197
$\mu_{1}=4.5$ | 2.25 | 0.86 | 0.89 | 0.570 | 0.232 | 0.198
$\mu_{1}=5$ | 2.32 | 0.82 | 0.86 | 0.573 | 0.230 | 0.197
$\mu_{1}=10$ | 2.80 | 0.59 | 0.61 | 0.641 | 0.190 | 0.169
$\mu_{1}=30$ | 3.43 | 0.28 | 0.29 | 0.814 | 0.097 | 0.089
Table 12: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{2}=5$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$\mu_{1}=3.4$ | 46.31 | 46.08 | 46.08 | 22.60 | 47.24 | 47.24 | 0.9999 | 0.9999 | 0.9999 | 0.0011
$\mu_{1}=3.5$ | 14.46 | 14.23 | 14.23 | 6.99 | 14.74 | 14.74 | 0.9990 | 0.9990 | 0.9990 | 0.0036
$\mu_{1}=4$ | 3.42 | 3.20 | 3.20 | 1.58 | 3.48 | 3.48 | 0.9824 | 0.9824 | 0.9860 | 0.0149
$\mu_{1}=4.5$ | 2.02 | 1.81 | 1.81 | 0.90 | 2.06 | 2.06 | 0.9525 | 0.9525 | 0.9628 | 0.0248
$\mu_{1}=5$ | 1.47 | 1.26 | 1.26 | 0.63 | 1.49 | 1.49 | 0.9158 | 0.9158 | 0.9348 | 0.0351
$\mu_{1}=10$ | 0.48 | 0.33 | 0.33 | 0.17 | 0.55 | 0.55 | 0.6079 | 0.6079 | 0.6973 | 0.0877
$\mu_{1}=30$ | 0.17 | 0.09 | 0.09 | 0.05 | 0.32 | 0.32 | 0.3774 | 0.3774 | 0.4718 | 0.1524
Table 13: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{2}=5$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{1}=q_{2}=q_{3}=1$ $\mu_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$\mu_{1}=3.4$ | 2.05 | 0.975 | 0.975 | 0.607 | 0.1965 | 0.1965
$\mu_{1}=3.5$ | 2.07 | 0.965 | 0.965 | 0.605 | 0.1975 | 0.1975
$\mu_{1}=4$ | 2.17 | 0.915 | 0.915 | 0.60 | 0.20 | 0.20
$\mu_{1}=4.5$ | 2.25 | 0.875 | 0.875 | 0.60 | 0.20 | 0.20
$\mu_{1}=5$ | 2.32 | 0.84 | 0.84 | 0.601 | 0.1995 | 0.1995
$\mu_{1}=10$ | 2.81 | 0.595 | 0.595 | 0.665 | 0.1675 | 0.1675
$\mu_{1}=30$ | 3.45 | 0.275 | 0.275 | 0.83 | 0.085 | 0.085
Table 14: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{1}=4.5$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{2}=q_{3}=1$ $q_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$q_{1}=1$ | 48.64 | 48.69 | 48.62 | 35.46 | 39.90 | 34.53 | 0.9999 | 0.9999 | 0.9999 | 0.0003
$q_{1}=3$ | 38.08 | 38.16 | 38.10 | 26.88 | 31.83 | 27.53 | 0.9998 | 0.9998 | 0.9998 | 0.0005
$q_{1}=8$ | 34.01 | 34.10 | 34.04 | 23.59 | 28.73 | 24.82 | 0.9998 | 0.9998 | 0.9998 | 0.0006
$q_{1}=13$ | 33.00 | 33.10 | 33.04 | 22.79 | 27.96 | 24.16 | 0.9998 | 0.9998 | 0.9998 | 0.0007
$q_{1}=20$ | 32.43 | 32.53 | 32.46 | 22.33 | 27.52 | 23.77 | 0.9998 | 0.9998 | 0.9998 | 0.0007
$q_{1}=30$ | 32.07 | 32.17 | 32.10 | 22.04 | 27.25 | 23.54 | 0.9998 | 0.9998 | 0.9998 | 0.0007
$q_{1}=100$ | 31.56 | 31.66 | 31.60 | 21.63 | 26.86 | 23.20 | 0.9998 | 0.9998 | 0.9998 | 0.0007
Table 15: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{1}=4.5$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=p_{2}=p_{3}=\frac{1}{3}$, $q_{2}=q_{3}=1$ $q_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$q_{1}=1$ | 1.37 | 1.22 | 1.41 | 0.307 | 0.409 | 0.284
$q_{1}=3$ | 1.42 | 1.20 | 1.38 | 0.318 | 0.402 | 0.280
$q_{1}=8$ | 1.44 | 1.19 | 1.37 | 0.323 | 0.398 | 0.279
$q_{1}=13$ | 1.45 | 1.18 | 1.37 | 0.323 | 0.398 | 0.279
$q_{1}=20$ | 1.45 | 1.18 | 1.37 | 0.326 | 0.397 | 0.277
$q_{1}=30$ | 1.46 | 1.18 | 1.36 | 0.326 | 0.397 | 0.277
$q_{1}=100$ | 1.46 | 1.18 | 1.36 | 0.327 | 0.396 | 0.277
Table 16: Numerical results for $\mathbb{E}[L_{i}]$, $\mathbb{E}[W_{i}]$ and $Cor(L_{i},L_{j})$ when $\lambda=4$, $\mu_{1}=4.2$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{2}=q_{3}=1$ $q_{1}$ | $\mathbb{E}[L_{1}]$ | $\mathbb{E}[L_{2}]$ | $\mathbb{E}[L_{3}]$ | $\mathbb{E}[W_{1}]$ | $\mathbb{E}[W_{2}]$ | $\mathbb{E}[W_{3}]$ | $Cor(L_{1},L_{2})$ | $Cor(L_{1},L_{3})$ | $Cor(L_{2},L_{3})$ | $GI$
---|---|---|---|---|---|---|---|---|---|---
$q_{1}=1$ | 105.10 | 104.95 | 104.90 | 47.79 | 123.78 | 110.07 | 0.9999 | 0.9999 | 0.9999 | 0.0004
$q_{1}=3$ | 83.44 | 83.30 | 83.25 | 37.32 | 100.22 | 89.25 | 0.9999 | 0.9999 | 0.9999 | 0.0005
$q_{1}=8$ | 74.82 | 74.87 | 74.82 | 33.25 | 91.06 | 81.15 | 0.9999 | 0.9999 | 0.9999 | 0.00001
$q_{1}=13$ | 72.90 | 72.77 | 72.72 | 32.24 | 88.79 | 79.14 | 0.9999 | 0.9999 | 0.9999 | 0.0005
$q_{1}=20$ | 71.71 | 71.58 | 71.53 | 31.66 | 87.50 | 77.99 | 0.9999 | 0.9999 | 0.9999 | 0.0006
$q_{1}=30$ | 70.96 | 70.83 | 70.78 | 31.30 | 86.68 | 77.27 | 0.9999 | 0.9999 | 0.9999 | 0.0006
$q_{1}=100$ | 69.89 | 69.78 | 69.72 | 30.79 | 85.53 | 76.26 | 0.9999 | 0.9999 | 0.9999 | 0.0005
Table 17: Numerical results for $\lambda^{i}_{eff}$ and $\gamma_{i}$ when $\lambda=4$, $\mu_{1}=4.2$, $\mu_{2}=3$, $\mu_{3}=5$, $p_{1}=0.9999$, $p_{2}=p_{3}=0.00005$, $q_{2}=q_{3}=1$ $q_{1}$ | $\lambda^{1}_{eff}$ | $\lambda^{2}_{eff}$ | $\lambda^{3}_{eff}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$
---|---|---|---|---|---|---
$q_{1}=1$ | 2.20 | 0.85 | 0.95 | 0.526 | 0.283 | 0.191
$q_{1}=3$ | 2.24 | 0.83 | 0.93 | 0.535 | 0.278 | 0.187
$q_{1}=8$ | 2.26 | 0.82 | 0.92 | 0.540 | 0.275 | 0.185
$q_{1}=13$ | 2.26 | 0.82 | 0.92 | 0.541 | 0.274 | 0.185
$q_{1}=20$ | 2.26 | 0.82 | 0.92 | 0.542 | 0.274 | 0.184
$q_{1}=30$ | 2.27 | 0.82 | 0.91 | 0.543 | 0.273 | 0.184
$q_{1}=100$ | 2.27 | 0.82 | 0.91 | 0.543 | 0.273 | 0.184

##  6 Conclusions

This paper studies a queueing system combining two operating policies: Serve
the Longest Queue and Join the Shortest Queue, that were separately studied by
J.W. Cohen (Cohen, (1987) and Cohen, (1998)), respectively. The combined JSQ-
SLQ model (first introduced in Perel et al., (2020) for a preemptive 2-queue
system) is analyzed for a 3-queue, single-server, Markovian system. The system
is formulated in an un-conventional way, where, instead of defining an un-
bounded 3-dimensional state space $(L_{1},L_{2},L_{3})$, where $L_{i}$
represents the number of customers in $Q_{i}$, $i=1,2,3$, the system is
characterized by the couple $L_{1}$ and $D=(I,L_{1}-L_{2},L_{1}-L_{3})$, where
$I=1,2,3$ indicates the location of the server. This leads to a 2-dimensional
state space with infinite dimension for $L_{1}$ but finite dimension for $D$.
Consequently, one can apply the probability generating functions method, as
well as the matrix geometric analysis, for QBD processes to derive the
system’s joint 2-dimensional probability mass function, the system’s stability
condition, and consequently, calculate its performance measures. Numerical
results are presented in various tables with wide range of parameter values.
It is shown that the main purpose of the JSQ-SLQ policy, namely balancing
between the queue sizes, is greatly achieved. This is further demonstrated by
the Gini Index, which is close to zero in most cases, while its highest value,
$GI=0.1524$, is in the extreme case when $p_{1}=0.9999$ (Table 12). Note that
$GI=0$ in the symmetric case (Table 6).

A natural extension of the presented model is to consider a preemptive single-
server JSQ-SLQ system with $N\geq 4$ non-symmetrical queues. Such a system can
be modeled as a 2-dimensional continuous-time Markovian process, $(L_{1},D)$,
where $L_{1}$ denotes the (possibly infinite) size of $Q_{1}$, while $D$ is a
vector where its first coordinate indicates the position of the server, and
its $j$’th coordinate, $2\leq j\leq N$, represents the difference
$L_{1}-L_{j}$. Although $D$ assumes values belonging to a finite set
$\mathfrak{D}$, the size of $\mathfrak{D}$ increases exponentially with $N$.
Specifically, consider the case when $L_{1}=0$. Then, either all queues are
empty and the server may reside in any of the $N$ queues, or there are exactly
$i$ queues with a single customer, $1\leq i\leq N-1$, where the server may
reside in any of these $i$ queues. Therefore, the size of the set
$\mathfrak{D}$ when $L_{1}=0$ is

$N+\sum_{i=1}^{N-1}{i\binom{N-1}{i}}=N+(N-1)\cdot 2^{N-2}.$

For the case where $L_{1}\geq 1$, the number of states in $\mathfrak{D}$ is
$N\cdot 2^{N-1}$. This follows since the server may reside in any of the $N$
queues, say $Q_{i}$, while, the number of customers in all other $N-1$ queues
$Q_{j}$, $j\neq i$, can only be equal to $L_{i}$ or to $L_{i}-1$. Thus, the
resulting system becomes analytically intractable.

## References

  * Adan et al., (2016) Adan, I. J., Boxma, O. J., Kapodistria, S., and Kulkarni, V. G. (2016).  The shorter queue polling model.  Annals of Operations Research, 241(1-2):167–200.
  * Adan et al., (2013) Adan, I. J., Kapodistria, S., and van Leeuwaarden, J. S. (2013).  Erlang arrivals joining the shorter queue.  Queueing Systems, 74(2-3):273–302.
  * (3) Adan, I. J., Wessels, J., and Zijm, W. (1991a).  Analysis of the asymmetric shortest queue problem.  Queueing Systems, 8(1):1–58.
  * (4) Adan, I. J., Wessels, J., and Zijm, W. (1991b).  Analysis of the asymmetric shortest queue problem with threshold jockeying.  Communications in Statistics. Stochastic Models, 7(4):615–627.
  * Artalejo and Gómez-Corral, (2008) Artalejo, J. R. and Gómez-Corral, A. (2008).  Retrial Queueing Systems: A Computational Approach.  Springer.
  * Avrachenkov et al., (2014) Avrachenkov, K., Nain, P., and Yechiali, U. (2014).  A retrial system with two input streams and two orbit queues.  Queueing Systems, 77(1):1–31.
  * Blanc, (2009) Blanc, J. P. C. (2009).  Bad luck when joining the shortest queue.  European Journal of Operational Research, 195(1):167–173.
  * Braverman, (2020) Braverman, A. (2020).  Steady-state analysis of the join-the-shortest-queue model in the halfin–whitt regime.  Mathematics of Operations Research, 45(3):1069–1103.
  * Bright and Taylor, (1995) Bright, L. and Taylor, P. G. (1995).  Calculating the equilibrium distribution in level dependent quasi-birth-and-death processes.  Stochastic Models, 11(3):497–525.
  * Cohen, (1987) Cohen, J. W. (1987).  A two-queue, one-server model with priority for the longer queue.  Queueing Systems, 2(3):261–283.
  * Cohen, (1998) Cohen, J. W. (1998).  Analysis of the asymmetrical shortest two-server queueing model.  International Journal of Stochastic Analysis, 11(2):115–162.
  * Dawson et al., (2019) Dawson, D. A., Tang, J., Zhao, Y. Q., et al. (2019).  Performance analysis of joining the shortest queue model among a large number of queues.  Asia-Pacific Journal of Operational Research (APJOR), 36(04):1–23.
  * Dimitriou, (2021) Dimitriou, I. (2021).  Analysis of the symmetric join the shortest orbit queue.  Operations Research Letters, 49(1):23–29.
  * Dodge, (2008) Dodge, Y. (2008).  The concise encyclopedia of statistics.  Springer Science & Business Media.
  * Eschenfeldt and Gamarnik, (2018) Eschenfeldt, P. and Gamarnik, D. (2018).  Join the shortest queue with many servers. the heavy-traffic asymptotics.  Mathematics of Operations Research, 43(3):867–886.
  * Flatto, (1989) Flatto, L. (1989).  The longer queue model.  Probability in the Engineering and Informational Sciences, 3(4):537–559.
  * Foley and McDonald, (2001) Foley, R. D. and McDonald, D. R. (2001).  Join the shortest queue: stability and exact asymptotics.  Annals of Applied Probability, pages 569–607.
  * Halfin, (1985) Halfin, S. (1985).  The shortest queue problem.  Journal of Applied Probability, 22(4):865–878.
  * Hanukov and Yechiali, (2021) Hanukov, G. and Yechiali, U. (2021).  Explicit solutions for continuous-time QBD processes by using relations between matrix geometric analysis and the probability generating functions method.  Probability in the Engineering and Informational Sciences, 35:565–580.
  * Harchol-Balter, (2013) Harchol-Balter, M. (2013).  Performance modeling and design of computer systems: queueing theory in action.  Cambridge University Press.
  * Hordijk and Koole, (1990) Hordijk, A. and Koole, G. (1990).  On the optimality of the generalized shortest queue policy.  Probability in the Engineering and Informational Sciences, 4(4):477–487.
  * Houtum et al., (1997) Houtum, G.-J. V., Adan, I., and Der wal, J. (1997).  The symmetric longest queue system.  Stochastic Models, 13(1):105–120.
  * Knessl and Yao, (2013) Knessl, C. and Yao, H. (2013).  On the nonsymmetric longer queue model: joint distribution, asymptotic properties, and heavy traffic limits.  Advances in Operations Research, 2013.
  * Latouche and Ramaswami, (1999) Latouche, G. and Ramaswami, V. (1999).  Introduction to matrix analytic methods in stochastic modeling.  Siam.
  * Maguluri et al., (2014) Maguluri, S. T., Hajek, B., and Srikant, R. (2014).  The stability of longest-queue-first scheduling with variable packet sizes.  IEEE Transactions on Automatic Control, 59(8):2295–2300.
  * Neuts, (1981) Neuts, M. F. (1981).  Matrix-geometric solutions in stochastic models: an algorithmic approach.  The Johns Hopkins University Press.
  * Pedarsani and Walrand, (2016) Pedarsani, R. and Walrand, J. (2016).  Stability of multiclass queueing networks under longest-queue and longest-dominating-queue scheduling.  Journal of Applied Probability, 53(2):421–433.
  * Perel et al., (2020) Perel, E., Perel, N., and Yechiali, U. (2020).  A polling system with ’join the shortest - serve the longest’ policy.  Computers & Operations Research, 114:104809.
  * Ravid et al., (2013) Ravid, R., Boxma, O. J., and Perry, D. (2013).  Repair systems with exchangeable items and the longest queue mechanism.  Queueing Systems, 73(3):295–316.
  * Winston, (1977) Winston, W. (1977).  Optimality of the shortest line discipline.  Journal of Applied Probability, 14(1):181–189.
  * Yao and Knessl, (2005) Yao, H. and Knessl, C. (2005).  On the infinite server shortest queue problem: symmetric case.  Stochastic Models, 21(1):101–132.
  * Yao and Knessl, (2006) Yao, H. and Knessl, C. (2006).  On the infinite server shortest queue problem: Non-symmetric case.  Queueing Systems, 52(2):157–177.