# Doubly Exponential Solution for Randomized Load Balancing Models with
General Service Times

Quan-Lin Li
School of Economics and Management Sciences
Yanshan University, Qinhuangdao 066004, P.R. China

(September 5, 2010)

###### Abstract

The randomized load balancing model (also called _supermarket model_) is now
being applied to the study of load balancing in data centers and multi-core
servers systems. It is very interesting to analyze the general service times
in the supermarket model, and specifically understand influence of the heavy-
tailed service times on the doubly exponential solution. Since the supermarket
model is a complex queueing system, the general service times make its
analysis more challenging than the exponential or PH service case. Up to now,
it still is an open problem whether or how the heavy-tailed service times can
disrupt the doubly exponential structure of the fixed point in the supermarket
model.

In this paper, we provide a novel and simple approach to study the supermarket
model with general service times. This approach is based on the supplementary
variable method used in analyzing stochastic models extensively. We organize
an infinite-size system of integral-differential equations by means of the
density dependent jump Markov process, and obtain a close-form solution:
doubly exponential structure, for the fixed point satisfying the system of
nonlinear equations, which is always a key in the study of supermarket models.
The fixed point is decomposited into two groups of information under a product
form: the arrival information and the service information. Based on this, we
indicate two important observations: the fixed point for the supermarket model
is different from the tail of stationary queue length distribution for the
ordinary M/G/1 queue, and the doubly exponential solution to the fixed point
can extensively exist even if the service time distribution is heavy-tailed.
Furthermore, we analyze the exponential convergence of the current location of
the supermarket model to its fixed point, and study the Lipschitz condition in
the Kurtz Theorem under general service times. Based on these analysis, one
can gain a new understanding how workload probing can help in load balancing
jobs with general service times such as heavy-tailed service.

Keywords: Randomized load balancing, supermarket model, density dependent jump
Markov process, fixed point, doubly exponential solution, heavy-tailed
distribution, exponential convergence, Lipschitz condition.

##  1 Introduction

Randomized load balancing, where a job is assigned to a server from a small
subset of randomly chosen servers, is very simple to implement, and can
surprisingly deliver better performance (for example reducing collisions,
waiting times, backlogs) in a number of applications, such as, data center,
hash tables, distributed memory machines, path selection in networks, and task
assignment at web servers. One useful model that has been extensively used to
study the randomized load balancing schemes is the supermarket model. In the
supermarket model, a key result by Vvedenskaya, Dobrushin and Karpelevich [31]
indicated that when each Poisson arriving job is assigned to the shortest one
of $d\geq 2$ randomly chosen queues with exponential service times, the
equilibrium queue length can decay doubly exponentially in the limit as the
population number $n\rightarrow\infty$, and the stationary fraction of queues
with at least $k$ customers is $\rho^{\frac{d^{k}-1}{d-1}}$, which is a
substantially exponential improvement over the case for $d=1$, where the tail
of stationary queue length distribution in the corresponding M/M/1 queue is
$\rho^{k}$.

The distributed load balancing strategies, in which individual job decisions
are based on information on a limited number of other processors, have been
studied analytically by Eager, Lazokwska and Zahorjan [5, 6, 7] and through
trace-driven simulations by Zhou [33]. Based on this, the supermarket model is
developed by queueing theory and Markov processes. Most of recent research
applied the density dependent jump Markov processes to deal with a simple
supermarket model with Poisson arrival processes and exponential service
times, a key result of which illustrates that there exists a unique fixed
point which is decreasing doubly exponentially. That approach used in the
literature relies on determining the behavior of the supermarket model as its
size grows to infinity, and its behavior is naturally described as a system of
differential equations, which leads to a closed form solution: doubly
exponential structure, of the fixed point. Readers may refer to, such as,
analyzing a basic and simple supermarket model by Azar, Broder, Karlin and
Upfal [2], Vvedenskaya, Dobrushin and Karpelevich [31], Mitzenmacher [19, 20].

Certain generalization of the supermarket model have been explored, for
example, simple variations by Mitzenmacher and Vöcking [27], Mitzenmacher [21,
22, 25], Vöcking [30], Mitzenmacher, Richa, and Sitaraman [26] and Vvedenskaya
and Suhov [32]; and analyzing load information by Mirchandaney, Towsley, and
Stankovic [28], Dahlin [3], Mitzenmacher [24, 26]. Furthermore, Martin and
Suhov [18], Martin [17], Suhov and Vvedenskaya [29] studied the supermarket
mall model by means of the fast Jackson network, where each node in a Jackson
network is replaced by $N$ parallel servers, and a job joins the shortest of
$d$ randomly chosen queues at the node to which it is directed. Luczak and
McDiamid [15, 16] studied the maximum queue length of the original supermarket
model with exponential service times when the service speed scales linearly
with the number of jobs in the queue. Li, Lui and Wang [11, 12] discussed the
supermarket model with PH service times and the supermarket model with
Markovian arrival processes, respectively. Readers may refer to an excellent
overview by Mitzenmacher, Richa, and Sitaraman [26].

This paper is interested in analyzing the supermarket model with general
service times, which is an open problem for determining whether or how the
heavy-tailed service times can disrupt the doubly exponential structure of the
fixed point. On the other hand, note that the supermarket model is a complex
queueing system and has much different characteristics from the ordinary
queueing systems, thus the general service times make its analysis more
challenging than the exponential or PH service case. Up to now, there has not
been an effective method to be able to deal with the supermarket model with
general service times.

The main contributions of the paper are threefold. The first one is to provide
a novel and simple approach to study the supermarket model with general
service times. This approach is based on the supplementary variable method but
is described as a new integral-differential structure for expressing and
computing the fraction of queues efficiently. Using the new approach, we setup
an infinite-size system of integral-differential equations, which makes
applications of the density dependent jump Markov processes to be able to deal
with the general distributions, such as general service times, involved in the
supermarket model. The second one is to obtain a close-form solution: doubly
exponential structure, for the fixed point satisfying the system of nonlinear
equations, which is always a key in the study of supermarket models.
Furthermore, this paper analyzes the exponential convergence of the current
location of the supermarket model to its fixed point, and studies the
Lipschitz condition in the Kurtz Theorem under general service times. Also,
this paper provides numerical examples to illustrate the effectiveness of our
approach in analyzing the randomized load balancing schemes with the non-
exponential service requirements. The third one is to obtain that the fixed
point is decomposited into two groups of information under a product form: the
arrival information and the service information. Based on this, we indicate
three important observations:

(a)

The fixed point for the supermarket model is different from the tail of
stationary queue length distribution for the ordinary M/G/1 queue, because the
fixed point is light-tailed but the stationary queue length is heavy-tail if
the service times are heavy-tailed. Note that such a difference is illustrated
in this paper for the first time, while it can not be observed in the
literature for the supermarket model with Poisson arrivals and exponential
service times, e.g., see Mitzenmacher, Richa, and Sitaraman [26].

(b)

The doubly exponential solution to the fixed point can extensively exist even
if the service time distribution is heavy-tailed. This is an answer of the
above open problem to illustrate the role played by the heavy-tailed service
time distribution in the doubly exponential solution to the fixed point.

(c)

The doubly exponential solution to the fixed point is not unique for a more
general supermarket model. Note that we give three different doubly
exponential solutions in the supermarket model with Poisson arrivals and PH
service times, thus it is very interesting to provide all the doubly
exponential solutions for a more general supermarket model.

Based on this, one can gain the new and important understanding how the
workload probing can help in load balancing jobs with general service times
such as heavy-tailed service.

The remainder of this paper is organized as follows. In Section 2, we first
describe a supermarket model with general service times, which is always
useful in the study of randomized load balancing schemes. Then the supermarket
model is expressed as a systems of integral-differential equations in terms of
the density dependent jump Markov processes. In Section 3, we first introduce
a fixed point of the system of integral-differential equations, and set up a
system of nonlinear equations satisfied by the fixed point. Then we provide a
close-form solution: doubly exponential structure, to the system of nonlinear
equations. In Section 4, we provide a necessary discussion on the key
parameter $\theta$ used in the doubly exponential structure, and indicate that
the doubly exponential solution to the fixed point extensively exists even if
the service time distribution is heavy-tailed. In Section 5, we give three
methods to analyze the supermarket model with Poisson arrivals and PH service
times, and provide three different ways to determine the key parameter
$\theta$ and compute the doubly exponential solution to the fixed point. We
show that the doubly exponential solution to the fixed point is not unique for
a more general supermarket. In Section 6, we study the exponential convergence
of the current location of the supermarket model to its fixed point. Not only
does the exponential convergence indicates the existence of the fixed point,
but it also explains such a convergent process is very fast. In Section 7, we
apply the Kurtz Theorem to study the supermarket model with the general
service times, and analyze the Lipschitz condition with respect to general
service times. Some concluding remarks are given in Section 8.

##  2 Supermarket Model

In this section, we first describe a supermarket model with general service
times. Then we provide a novel and simple approach to setup an infinite-size
system of integral-differential equations based on the density dependent jump
Markov processes. Note that this approach is based on the supplementary
variable method but is described as a new integral-differential structure so
that the corresponding boundary conditions are written in a different version.

Let us formally describe the supermarket model, which is abstracted as a
multi-server multi-queue stochastic system. Customers arrive at a queueing
system of $n>1$ servers as a Poisson process with an average arrival rate
$n\lambda$ for $\lambda>0$. The service time $\chi_{k}$ of the $k$th customer
is general with the distribution function

$G\left(x\right)=P\left\\{\chi_{k}\leq
x\right\\}=1-\exp\left\\{-\int_{0}^{x}\mu\left(y\right)dy\right\\},$

where all the random variables $\chi_{k}$ for $k\geq 1$ are i.i.d. with the
mean $E\left[\chi_{k}\right]=1/\mu$. Each arriving customer chooses $d\geq 1$
servers independently and uniformly at random from these $n$ servers, and
waits for service at the server which currently contains the fewest number of
customers. If there is a tie, servers with the fewest number of customers will
be chosen randomly. All customers in any service center will be served in the
First-Come-First-Served (FCFS) manner. Figure 1 simply shows such a
supermarket model.

Figure 1: working structure of the supermarket model

In the study of supermarket models, it is necessary for us to study general
service time distributions, for example, heavy-tailed distributions. Not only
because the general distribution makes analysis of the supermarket models more
difficult and challenging than those in the literature for the exponential or
PH service case, but it also allows us to model more realistic systems and
understand their performance implication under the randomized load balancing
strategy. As indicated in [8], the process times of many parallel jobs, in
particular, jobs to data centers, tend to be non-exponential. Unless we state
otherwise, we assume that all the random variables defined above are
independent, and that the system is operating under the condition:
$\rho=\lambda/\mu<1$.

###### Lemma 1

The supermarket model with general service times is stable if
$\rho=\lambda/\mu<1.$

Proof: Let $Q_{O}\left(t\right)$ and $Q_{S}\left(t\right)$ be the queue
lengths of the ordinary M/G/1 queue and of an arbitrary server in the
supermarket model at time $t$, respectively. Note that in the supermarket
model, each customer chooses $d$ servers independently and uniformily at
random, and the queue length of the entering server is currently shorten, it
is easy to see that for each $t\geq 0$,

$0\leq
Q_{S}\left(t\right)\underset{\text{st}}{\leq}Q_{O}\left(t\right)\text{.}$ (1)

Since the ordinary M/G/1 queue is stable if $\rho=\lambda/\mu<1$, it follows
from (1) that the supermarket model with general service times is stable if
$\rho=\lambda/\mu<1$. This completes the proof.

For $k\geq 1$, we define $n_{k}\left(t,x\right)$d$x$ as the number of queues
with at least $k$ customers and the residual service time of each server be in
the interval $[x,x+$d$x)$ at time $t\geq 0$. Clearly, $0\leq
n_{k}\left(t,x\right)\leq n$ for $x\geq 0$ and $1\leq k\leq n$. Let

$s_{k,n}\left(t,x\right)=\frac{n_{k}\left(t,x\right)}{n},$

which is the density function of the fraction of queues with at least $k$
customers and the residual service time of each server be $x$. We write

$S_{k}\left(t,x\right)=\lim_{n\rightarrow\infty}s_{k,n}\left(t,x\right),\text{
\ for }k\geq 1.$

We define $n_{0,n}\left(t\right)$ as the number of queues with at least $0$
customers at time $t\geq 0$. Clearly, $n_{0,n}\left(t\right)=n$. Let

$s_{0,n}\left(t\right)=\frac{n_{0,n}\left(t\right)}{n}.$

Then $s_{0,n}\left(t\right)=1$ for all $t\geq 0$ and

$S_{0}\left(t\right)=\lim_{n\rightarrow\infty}s_{0,n}\left(t\right)=1.$

Let $V\left(t\right)$ be the fraction of queues with zero customer at time
$t$. Then

$S_{0}\left(t\right)=V\left(t\right)+\int_{0}^{+\infty}S_{1}\left(t,x\right)\text{d}x.$

Thus we have

$S\left(t,x\right)=\left(S_{0}\left(t\right),S_{1}\left(t,x\right),S_{2}\left(t,x\right),\ldots\right).$

The following proposition shows that the sequence
$\left\\{S_{k}\left(t,x\right)\right\\}$ is monotone increasing for $k\geq 1$,
while its proof is easily by means of the definition of
$S_{k}\left(t,x\right)$ for $k\geq 1$.

###### Proposition 1

For $1\leq k<l$

$S_{l}\left(t,x\right)<S_{k}\left(t,x\right)<S_{0}\left(t\right)=1$

and

$\int_{0}^{+\infty}S_{l}\left(t,x\right)\text{d}x<\int_{0}^{+\infty}S_{k}\left(t,x\right)\text{d}x<S_{0}\left(t\right)=1.$

Now, we setup a system of integral-differential equations by means of the
density dependent jump Markov process. To that end, we provide an example with
$k\geq 2$ to indicate how to derive the system of integral-differential
equations.

Consider the supermarket model with $n$ queues, and determine the expected
change in the number of servers with at least $k$ customers and the residual
service time of each server be $x$ over a small time period of length d$t$.
The probability that a customer arriving during this time period is
$n\lambda$d$t$, and the probability that an arriving customer joins a queue of
size $k-1$ is given by
$\int_{0}^{+\infty}s_{k-1,n}^{d}\left(t,x\right)$d$x-\int_{0}^{+\infty}s_{k,n}^{d}\left(t,x\right)$d$x$.
Thus, the probability that during this time period, any arriving customer
joins a queue of size $k-1$ is given by

$n\lambda\text{d}t\cdot\left[\int_{0}^{+\infty}s_{k-1,n}^{d}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}s_{k,n}^{d}\left(t,x\right)\text{d}x\right].$

Similarly, the probability that a customer leaves a server of size $k$ is
given by

$ndt\cdot\left[\int_{0}^{+\infty}\mu\left(x\right)s_{k,n}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)s_{k+1,n}\left(t,x\right)\text{d}x\right].$

Therefore we can obtain

$\displaystyle\frac{\text{d}\int_{0}^{+\infty}n_{k}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle
n\lambda\left[\int_{0}^{+\infty}s_{k-1,n}^{d}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}s_{k,n}^{d}\left(t,x\right)\text{d}x\right]$
$\displaystyle+n\left[\int_{0}^{+\infty}\mu\left(x\right)s_{k,n}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)s_{k+1,n}\left(t,x\right)\text{d}x\right],$

which leads to

$\displaystyle\frac{\text{d}\int_{0}^{+\infty}s_{k,n}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda\left[\int_{0}^{+\infty}s_{k-1,n}^{d}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}s_{k,n}^{d}\left(t,x\right)\text{d}x\right]$
$\displaystyle+\left[\int_{0}^{+\infty}\mu\left(x\right)s_{k,n}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)s_{k+1,n}\left(t,x\right)\text{d}x\right].$
(2)

Taking $n\rightarrow\infty$ in the both sides of (2), we have

$\displaystyle\frac{\text{d}\int_{0}^{+\infty}S_{k}\left(t\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda\left[\int_{0}^{+\infty}S_{k-1}^{d}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}S_{k}^{d}\left(t,x\right)\text{d}x\right]$
$\displaystyle+\left[\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(t,x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)S_{k+1}\left(t,x\right)\text{d}x\right].$
(3)

Using a similar analysis to that for deriving Equation (3), we can easily
obtain a system of integral-differential equations for the fraction density
vector $S\left(t,x\right)$ as follows:

$S_{0}\left(t\right)=1\text{ for all }t\geq 0,$ (4)
$\frac{\mathtt{d}}{\text{d}t}S_{0}\left(t\right)=-\lambda
S_{0}^{d}\left(t\right)+\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(t,x\right)\text{d}x,$
(5)
$\displaystyle\frac{\mathtt{d}\int_{0}^{+\infty}S_{1}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda
S_{0}^{d}\left(t\right)-\lambda\int_{0}^{+\infty}S_{1}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{2}\left(t,x\right)\text{d}x,$
(6)

and for $k\geq 2$,

$\displaystyle\frac{\mathtt{d}\int_{0}^{+\infty}S_{k}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda\int_{0}^{+\infty}S_{k-1}^{d}\left(t,x\right)\text{d}x-\lambda\int_{0}^{+\infty}S_{k}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{k+1}\left(t,x\right)\text{d}x.$
(7)

###### Remark 1

When there are $n$ servers in the supermarket model, it is necessary to give a
finite-size system of integral-differential equations for the fraction density
vector
$S^{\left(n\right)}\left(t,x\right)=\left(s_{0,n}\left(t\right),s_{1,n}\left(t,x\right),\ldots,s_{n,n}\left(t,x\right)\right)$
as follows:

$s_{0,n}\left(t\right)=1\text{ for all }t\geq 0,$
$\frac{\mathtt{d}}{\text{d}t}s_{0,n}\left(t\right)=-\lambda
s_{0,n}^{d}\left(t\right)+\int_{0}^{+\infty}\mu\left(x\right)s_{1,n}\left(t,x\right)\text{d}x,$
$\displaystyle\frac{\mathtt{d}\int_{0}^{+\infty}s_{1,n}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda
s_{0,n}^{d}\left(t\right)-\lambda\int_{0}^{+\infty}s_{1,n}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)s_{1,n}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)s_{2,n}\left(t,x\right)\text{d}x,$

and for $n\geq k\geq 2$,

$\displaystyle\frac{\mathtt{d}\int_{0}^{+\infty}s_{k,n}\left(t,x\right)\text{d}x}{\text{d}t}=$
$\displaystyle\lambda\int_{0}^{+\infty}s_{k-1,n}^{d}\left(t,x\right)\text{d}x-\lambda\int_{0}^{+\infty}s_{k,n}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)s_{k,n}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)s_{k+1,n}\left(t,x\right)\text{d}x.$

##  3 Doubly Exponential Solution

In this section, we discuss the fixed point of the system of integral-
differential equations in Equations (4) to (7), and set up a system of
nonlinear equations satisfied by the fixed point. Also, we provide a closed-
form solution: doubly exponential structure, to the system of nonlinear
equations.

A row vector
$\pi\left(x\right)=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is called a fixed point of the fraction density vector
$S\left(t,x\right)=\left(S_{0}\left(t\right),S_{1}\left(t,x\right),S_{2}\left(t,x\right),\ldots\right)$
if there exists a $t_{0}\geq 0$ such that $S_{0}\left(t\right)=\pi_{0}$ and
$S_{k}\left(t,x\right)=\pi_{k}\left(x\right)$ for all $t\geq t_{0}$ and $k\geq
1$. It is easy to see that if $\pi\left(x\right)$ is a fixed point of the
fraction density vector $S\left(t,x\right)$ for all $t\geq t_{0}$, then

$\frac{\mathtt{d}}{\text{d}t}S_{0}\left(t\right)_{|t\geq t_{0}}=0$

and for $k\geq 1$

$\frac{\mathtt{d}}{\text{d}t}S_{k}\left(t,x\right)_{|t\geq t_{0}}=0$

which leads to

$\int_{0}^{+\infty}\frac{\mathtt{d}}{\text{d}t}S_{k}\left(t,x\right)_{|t\geq
t_{0}}\text{d}x=0.$ (8)

Since for $k\geq 1$

$0\leq S_{k}\left(t,x\right)\leq 1,$

using the Dominated Convergence Theorem we obtain

$\frac{\mathtt{d}}{\text{d}t}\int_{0}^{+\infty}S_{k}\left(t,x\right)_{|t\geq
t_{0}}\text{d}x=0.$

Therefore, if
$\pi\left(x\right)=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is a fixed point of the fraction density vector
$S\left(t,x\right)=\left(S_{0}\left(t\right),S_{1}\left(t,x\right),S_{2}\left(t,x\right),\ldots\right)$
for all $t\geq t_{0}$, then the system of integral-differential equations (4)
to (7) can be simplified as

$\pi_{0}=1$ (9)
$-\lambda\pi_{0}^{d}+\int_{0}^{+\infty}\mu\left(x\right)\pi_{1}\left(x\right)\text{d}x=0,$
(10)
$\lambda\pi_{0}^{d}\left(t\right)-\lambda\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)\pi_{1}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)\pi_{2}\left(x\right)\text{d}x=0,$
(11)

and for $k\geq 2$,

$\lambda\int_{0}^{+\infty}\pi_{k-1}^{d}\left(x\right)\text{d}x-\lambda\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x-\int_{0}^{+\infty}\mu\left(x\right)\pi_{k}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)\pi_{k+1}\left(x\right)\text{d}x=0.$
(12)

In what follows we derive a closed-form expression for
$\pi\left(x\right)=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$.
It follows from Equations (9) and (10) that

$\int_{0}^{+\infty}\mu\left(x\right)\pi_{1}\left(x\right)\text{d}x=\lambda.$
(13)

To solve Equation (13), using the fact that
$\int_{0}^{+\infty}\mu\left(x\right)\overline{G}\left(x\right)$d$x=1$ we have

$\pi_{1}\left(x\right)=\lambda\overline{G}\left(x\right)=\rho\cdot\mu\overline{G}\left(x\right).$
(14)

Based on the fact that $\pi_{0}=1$ and
$\pi_{1}\left(x\right)=\rho\cdot\mu\overline{G}\left(x\right)$, it follows
from Equations (11) and (13) that

$-\lambda\rho^{d}\cdot\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)\pi_{2}\left(x\right)\text{d}x=0.$

Let
$\theta=\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}$d$x$,
and we assume that $0<\theta<+\infty$. Then

$\int_{0}^{+\infty}\mu\left(x\right)\pi_{2}\left(x\right)\text{d}x=\lambda\theta\rho^{d}.$
(15)

Using a similar analysis on Equation (15), we have

$\pi_{2}\left(x\right)=\lambda\theta\rho^{d}\overline{G}\left(x\right)=\theta\rho^{d+1}\cdot\mu\overline{G}\left(x\right).$
(16)

Based on $\pi_{1}\left(x\right)=\rho\cdot\mu\overline{G}\left(x\right)$ and
$\pi_{2}\left(x\right)=\theta\rho^{d+1}\cdot\mu\overline{G}\left(x\right)$, we
can compute

$\lambda\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}x=\lambda\theta\rho^{d},$
$\int_{0}^{+\infty}\mu\left(x\right)\pi_{2}\left(x\right)\text{d}x=\theta\rho^{d+1}\mu\int_{0}^{+\infty}\mu\left(x\right)\overline{G}\left(x\right)\text{d}x=\lambda\theta\rho^{d}$

and

$\lambda\int_{0}^{+\infty}\pi_{2}^{d}\left(x\right)\text{d}x=\lambda\theta^{d}\rho^{d^{2}+d}\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}\text{d}x=\lambda\theta^{d+1}\rho^{d^{2}+d},$

thus it follows from Equation (12) that for $k=2$,

$\int_{0}^{+\infty}\mu\left(x\right)\pi_{3}\left(x\right)\text{d}x=\lambda\theta^{d+1}\rho^{d^{2}+d},$

which leads to

$\pi_{3}\left(x\right)=\theta^{d+1}\rho^{d^{2}+d}\overline{G}\left(x\right)=\theta^{d+1}\rho^{d^{2}+d+1}\cdot\mu\overline{G}\left(x\right).$
(17)

Based on the above analysis for the simple expressions $\pi_{k}\left(x\right)$
for $k=1,2$ and $3$, we can summarize the following theorem.

###### Theorem 1

The fixed point
$\pi=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is given by

$\pi_{0}=1,$ $\pi_{1}\left(x\right)=\rho\cdot\mu\overline{G}\left(x\right)$

and for $k\geq 2,$

$\pi_{k}\left(x\right)=\theta^{d^{k-2}+d^{k-3}+\cdots+1}\rho^{d^{k-1}+d^{k-2}+\cdots+1}\cdot\mu\overline{G}\left(x\right),$
(18)

or

$\pi_{k}\left(x\right)=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\mu\overline{G}\left(x\right).$
(19)

Proof By induction, one can easily derive the above result.

It is clear from (16) and (17) that Equation (18) or (19) is correct for the
cases with $l=2,3$. Now, we assume that Equation (19) is correct for the cases
with $l=k$. Then

$\lambda\int_{0}^{+\infty}\pi_{k-1}^{d}\left(x\right)\text{d}x=\lambda\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-d}{d-1}},$
$\int_{0}^{+\infty}\mu\left(x\right)\pi_{k}\left(x\right)\text{d}x=\lambda\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-d}{d-1}}$

and

$\lambda\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x=\lambda\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k+1}-d}{d-1}},$

it follows from Equation (12) that

$\int_{0}^{+\infty}\mu\left(x\right)\pi_{k+1}\left(x\right)\text{d}x=\lambda\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k+1}-d}{d-1}}.$

Thus, for $l=k+1$ we have

$\pi_{k+1}\left(x\right)=\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k+1}-1}{d-1}}\cdot\mu\overline{G}\left(x\right).$

This completes the proof.

Let
$\widetilde{\theta}=\int_{0}^{+\infty}\left[\overline{G}\left(x\right)\right]^{d}$d$x$.
Then $\theta=\mu^{d}\widetilde{\theta}$. The following corollary provides
another expression for the fixed point.

###### Corollary 2

$\pi_{0}=1$

and for $k\geq 1$

$\pi_{k}\left(x\right)=\lambda^{\frac{d^{k}-1}{d-1}}\cdot\left\\{\widetilde{\theta}^{\frac{d^{k-1}-1}{d-1}}\overline{G}\left(x\right)\right\\}.$

It is easy to see from Corollary 2 that the fixed point is decomposited into
two groups of information under a product form: the arrival information and
the service information. At the same time, the service information indicates
that the doubly exponential solution to the fixed point must exist for
$0<\mu<+\infty$, even if the service times are heavy-tailed.

The following corollary provides an upper bound for the fixed point.

###### Corollary 3

For $k\geq 1$ and $x\geq 0,$

$\pi_{k}\left(x\right)<\int_{0}^{+\infty}\pi_{k}\left(x\right)\text{d}x<\rho^{\frac{d^{k-1}-1}{d-1}}\frac{\lambda^{d^{k}}}{\mu}.$

Proof: Note that $0\leq\overline{G}\left(x\right)\leq 1$, we have

$\widetilde{\theta}=\int_{0}^{+\infty}\left[\overline{G}\left(x\right)\right]^{d}\text{d}x<\int_{0}^{+\infty}\overline{G}\left(x\right)\text{d}x=\frac{1}{\mu}.$

It follows from Corollary 2 that

$\pi_{k}\left(x\right)<\int_{0}^{+\infty}\pi_{k}\left(x\right)\text{d}x<\int_{0}^{+\infty}\lambda^{\frac{d^{k}-1}{d-1}}\widetilde{\theta}^{\frac{d^{k-1}-1}{d-1}}\overline{G}\left(x\right)\text{d}x=\rho^{\frac{d^{k-1}-1}{d-1}}\frac{\lambda^{d^{k}}}{\mu}.$

This completes the proof.

Now, we compute the expected sojourn time $T_{d}$ which a tagged arriving
customer spends in the supermarket model. For the general service times, a
tagged arriving customer is the $k$th customer in the corresponding queue with
the following probability

$\int_{0}^{+\infty}\pi_{k-1}^{d}\left(x\right)\text{d}x-\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x=\theta^{\frac{d^{k-2}-1}{d-1}}\rho^{\frac{d^{k-1}-1}{d-1}}-\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}.$

When $k\geq 1$, the head customer in the queue has been served, and so its
service time is residual and is denoted as $X_{R}$. Under the stationary
setting, we have

$P\left\\{X_{R}\leq
x\right\\}=\int_{0}^{x}\left[\mu\overline{G}\left(y\right)\right]\text{d}y$

with

$E\left[X_{R}\right]=\int_{0}^{+\infty}\int_{x}^{+\infty}\left[\mu\overline{G}\left(y\right)\right]\text{d}y\text{d}x.$

Thus it is easy to see that the expected sojourn time of the tagged arriving
customer is given by

$\displaystyle E\left[T_{d}\right]=$ $\displaystyle\left[\pi_{0}^{\odot
d}-\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}x\right]E\left[X\right]$
$\displaystyle+\sum_{k=1}^{\infty}\left[\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x-\int_{0}^{+\infty}\pi_{k+1}^{d}\left(x\right)\text{d}x\right]\left[E\left[X_{R}\right]+kE\left[X\right]\right]$
$\displaystyle=$
$\displaystyle\left[1-\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}x\right]E\left[X\right]+\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}xE\left[X_{R}\right]$
$\displaystyle+E\left[X\right]\sum_{k=1}^{\infty}k\left[\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x-\int_{0}^{+\infty}\pi_{k+1}^{d}\left(x\right)\text{d}x\right]$
$\displaystyle=$
$\displaystyle\left\\{E\left[X_{R}\right]-E\left[X\right]\right\\}\int_{0}^{+\infty}\pi_{1}^{d}\left(x\right)\text{d}x+E\left[X\right]\left\\{1+\sum_{k=1}^{\infty}\int_{0}^{+\infty}\pi_{k}^{d}\left(x\right)\text{d}x\right\\}$
$\displaystyle=$
$\displaystyle\theta\rho^{d}\left\\{E\left[X_{R}\right]-E\left[X\right]\right\\}+E\left[X\right]\left[\sum_{k=1}^{\infty}\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k}-d}{d-1}}\right].$

If the service times are exponential, then
$E\left[X_{R}\right]=E\left[X\right]$, thus we obtain

$E\left[T_{d}\right]=\frac{1}{\mu}\left[\sum_{k=1}^{\infty}\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k}-d}{d-1}}\right],$

which is the same as Corollary 3.8 in Mitzenmacher [20].

We consider a computational example for the expected sojourn time in the
supermarket model with an Erlang service time distribution
$E\left(m,\mu\right)$, where $m=2,\mu=1,d=2$. Figure 3 shows how the the
expected sojourn time depends on the arrival rate.

Figure 2: the expected sojourn time $E\left[T_{d}\right]$ depends on the
arrival rate $\lambda$

With the results from Equation (18) or (19), let us now provide some useful
discussions on the _asymptotic behavior_ of the fixed point
$\pi\left(x\right)=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$.
Note that we express $a_{k}\backsim O\left(b_{k}\right)$ if
$\lim_{k\rightarrow\infty}a_{k}/b_{k}=c\in\left(-\infty,0\right)\cup\left(0,+\infty\right)$.

###### Remark 2

If the general distribution $G\left(x\right)$ and its mean $1/\mu$ are given,
then
$\theta=\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}$d$x$
is a deterministic factor. We have

$\frac{\pi_{k}\left(x\right)}{\theta^{\frac{d^{k-1}-1}{d-1}}}\sim
O\left(\rho^{\frac{d^{k}-1}{d-1}}\right)\mu\overline{G}\left(x\right),\text{ \
as }k\rightarrow\infty.$

In this case, the heavy traffic should have a bigger influence on the
asymptotic behavior of the sequence
$\left\\{\frac{\pi_{k}\left(x\right)}{\theta^{\frac{d^{k-1}-1}{d-1}}}\right\\}$.

###### Remark 3

If $\rho$ is given, then

$\frac{\pi_{k}\left(x\right)}{\rho^{\frac{d^{k}-1}{d-1}}}\sim
O\left(\theta^{\frac{d^{k-1}-1}{d-1}}\right)\mu\overline{G}\left(x\right),\text{
\ as }k\rightarrow\infty.$

In this case, the maximal value $\theta_{\max}$ of the positive number
$\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}$d$x$ should
have a bigger influence on the asymptotic behavior of the sequence
$\left\\{\frac{\pi_{k}\left(x\right)}{\rho^{\frac{d^{k}-1}{d-1}}}\right\\}$.

##  4 A discussion for the key parameter $\theta$

In this section, we provide a necessary discussion for the key parameter
$\theta$ in the doubly exponential solution of Theorem 1. Based on this, for
the fixed point we give a new and important observation: the doubly
exponential solution to the fixed point can extensively exist for
$0<\mu<+\infty$, even if the service time distribution is heavy-tailed.

Note that

$\displaystyle\theta$
$\displaystyle=\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}\text{d}x$
$\displaystyle=\frac{\int_{0}^{+\infty}\left[\overline{G}\left(x\right)\right]^{d}\text{d}x}{\left[\int_{0}^{+\infty}\overline{G}\left(x\right)\text{d}x\right]^{d}},$
(20)

it is easy to see that $\theta=1$ if $d=1$. Thus, we need to analyze the case
for $k\geq 2$ as follows.

Since $0\leq\overline{G}\left(x\right)\leq 1$, we get that
$0\leq\left[\overline{G}\left(x\right)\right]^{d}\leq\overline{G}\left(x\right)\leq
1$, which leads to

$\int_{0}^{+\infty}\left[\overline{G}\left(x\right)\right]^{d}dx\leq\int_{0}^{+\infty}\overline{G}\left(x\right)dx=1/\mu.$

It is easy to see that $0<\theta<\mu^{d-1}$, and thus if $0<\mu<+\infty$, then
$0<\theta<+\infty$.

In what follows we analyze five simple and useful examples. In first two
examples, the service time distribution is light-tailed; while in the last
three examples, the service time distribution is heavy-tailed. Specifically,
the examples with heavy-tailed service times illustrate two important
observations: the first one indicates that the fixed point for the supermarket
model is different from the tail of stationary queue length distribution for
the ordinary M/G/1 queue, and the second one is to show that the doubly
exponential solution to the fixed point can exist extensively if the service
time mean is non-zero and finite.

Example one: Exponential distribution. Let $\overline{G}\left(x\right)=e^{-\mu
x}$. Then $\theta=\mu^{d-1}/d$. It is easy to see that when
$\mu>\sqrt[d-1]{d}$, $\theta>1$; when $\mu=\sqrt[d-1]{d}$, $\theta=1$; and
when $\mu<\sqrt[d-1]{d}$, $0<\theta<1$. If $d=2$, then $\theta$ is a linear
function of $\mu$, and if $d=3$, then $\theta$ is a nonlinear function of
$\mu$. Figures 3 and 4 show the functions $\theta=\mu/2$ and
$\theta=\mu^{2}/3$, respectively.

Figure 3: $\theta$ is a linear function of $\mu$ for $d=2$ Figure 4: $\theta$
is a nonlinear function of $\mu$ for $d=3$

Example two: Erlang distribution $E\left(m,\mu\right)$. Let
$\overline{G}\left(x\right)=e^{-\mu x}\sum_{k=0}^{m}\frac{\left(\mu
x\right)^{k}}{k!}$. Then $\theta$ is given by

$\theta=\left(\frac{\mu}{m}\right)^{d}\int_{0}^{+\infty}e^{-\mu
dx}\left[\sum_{k=0}^{m}\frac{\left(\mu x\right)^{k}}{k!}\right]^{d}\text{d}x.$

Let $\mu=1$. Table 1 lists how $\theta$ depends on the parameter pair
$\left(m,d\right)$. As seen from Table 1, $\theta$ is decreasing for each of
the two parameters $m$ and $d$.

Table 1: $\theta$ depends on the Erlang parameter pair $(m,d)$ $(m,d)$ | (2, 2) | (2, 5) | (2, 10) | (5, 2) | (10, 2) | (5, 5) | (10, 10)
---|---|---|---|---|---|---|---
$\theta$ | 0.52 | 0.19 | $9.15\times 10^{-2}$ | $4.13\times 10^{-2}$ | $9.48\times 10^{-4}$ | $1.11\times 10^{-3}$ | $6.51\times 10^{-10}$

Example three: Weibull distribution $W\left(\tau,\mu\right)$. Let
$\overline{G}\left(x\right)=\exp\left\\{-\left(\mu t\right)^{\tau}\right\\}$.
It is easy to check that the mean of the Weibull distribution is given by

$\frac{1}{\mu}\Gamma\left(1+\frac{1}{\tau}\right),$

which follows that $\theta$ is given by

$\theta=\frac{\mu^{d-1}}{d^{\frac{1}{\tau}}\left[\Gamma\left(1+\frac{1}{\tau}\right)\right]^{d-1}},$

where $\Gamma(\alpha)=\int_{0}^{+\infty}x^{\alpha-1}e^{-x}$d$x$. Obviously,
the Weibull distribution $W\left(\tau,\mu\right)$ is heavy-tailed if
$0<\tau<1$; and the Weibull distribution $W\left(\tau,\mu\right)$ is light-
tailed if $\tau>1$. To indicate the role played by the heavy-tailed parameter
$\tau$ for $0<\tau<1$, taking $\mu=5$ and $d=2$ we have

$\theta=\frac{5}{2^{\frac{1}{\tau}}\Gamma\left(1+\frac{1}{\tau}\right)}.$

Table 2 indicates how $\theta$ depends on the heavy-tailed parameter $\tau$,
such as, $0<\theta<1$ if $\tau=0.2$; $\theta>1$ if $\tau=0.9$. This example,
together with Theorem 1, illustrates an important observation that the fixed
point
$\pi=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is doubly exponential (clearly, it is light-tailed) even if the service time
distribution is heavy-tailed. Based on this, the the fixed point is different
from the tail of stationary queue length distribution of the ordinary M/G/1
queue, since for the ordinary M/G/1 queue, the stationary queue length
distribution is heavy-tailed if the service time distribution is heavy-tailed,
e.g., see Adler, Feldman and Taqqu [1].

Table 2: $\theta$ depends on the heavy-tailed parameter $\tau$ for $0<\tau<1$ $\tau$ | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---
$\theta$ | $1.3\times 10^{-3}$ | $5.3\times 10^{-2}$ | 0.27 | 0.63 | 1.05 | 1.47 | 1.86 | 2.19

Example four: Power law distribution. Let
$\overline{G}\left(x\right)=\left(\mu+x\right)^{-\alpha}$. If $0<\alpha\leq
1$, then the power law distribution does not exist the finite mean. In this
case, we can not setup the system of integral-differential equations for the
fraction density vector
$S\left(t,x\right)=\left(S_{0}\left(t\right),S_{1}\left(t,x\right),S_{2}\left(t,x\right),\ldots\right)$
which leads to the analysis for the fixed point. Thus we only deal with the
case with $\alpha>1$. Note that for each $\alpha>1$

$\int_{0}^{+\infty}\overline{G}\left(x\right)\text{d}x=\frac{1}{\mu}$

and

$\int_{0}^{+\infty}\left[\overline{G}\left(x\right)\right]^{d}\text{d}x=\frac{1}{\mu},$

thus we obtain $\theta=\mu^{d-1}$. It is easy to see that when $\mu>1$,
$\theta>1$; when $\mu=1$, $\theta=1$; and when $0<\mu<1$, $0<\theta<1$. It
follows from Theorem 1 that for $k\geq 1$

$\pi_{k}\left(x\right)=\mu^{d^{k-1}-1}\rho^{\frac{d^{k}-1}{d-1}}\cdot\mu\overline{G}\left(x\right).$

This indicates that the fixed point
$\pi=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is doubly exponential (of course, it is light-tailed) if the service time
distribution is power law.

Example five: Almost exponential distribution. Let
$\overline{G}\left(x\right)=\exp\left\\{-x\left(\ln
x\right)^{-\alpha}\right\\}$. Then it is easy to see that the almost
exponential distribution is heavy-tailed if $\alpha>0$

$\theta=\frac{\int_{0}^{+\infty}\exp\left\\{-dx\left(\ln
x\right)^{-\alpha}\right\\}\text{d}x}{\left[\int_{0}^{+\infty}\exp\left\\{-x\left(\ln
x\right)^{-\alpha}\right\\}\text{d}x\right]^{d}}$

Table 3 lists how $\theta$ depends on the parameter pair $(d,\alpha)$. As seen
from Table 3, $\theta$ is decreasing for each of the two parameters $d$ and
$\alpha$.

Table 3: $\theta$ depends on the parameter pair $(d,\alpha)$ $(d,\alpha)$ | (2, 2) | (4, 2) | (2, 4) | (4, 4)
---|---|---|---|---
$\theta$ | $2.24\times 10^{-2}$ | $2.01\times 10^{-4}$ | $3.44\times 10^{-5}$ | $1.18\times 10^{-13}$

##  5 The key parameter $\theta$ for PH Service Times

In this section, as an important example we provide three methods to analyze a
supermarket model with Poisson arrivals and PH service times. Our purpose is
to provide three different ways to determine the key parameter $\theta$ and
compute the doubly exponential solution to the fixed point. Also, we indicate
that the doubly exponential solution to the fixed point is not unique for a
more general supermarket model.

The supermarket model with Poisson arrivals and PH service times is described
as follows. Customers arrive at a queueing system of $n>1$ servers as a
Poisson process with arrival rate $n\lambda$ for $\lambda>0$. The service
times of these customers are of phase type with irreducible representation
$\left(\alpha,T\right)$ of order $m$. Each arriving customer chooses $d\geq 1$
servers independently and uniformly at random from these $n$ servers, and
waits for service at the server which currently contains the fewest number of
customers. If there is a tie, servers with the fewest number of customers will
be chosen randomly. All customers in any service center will be served in the
FCFS manner. For the PH service time distribution, we use the irreducible
representation $\left(\alpha,T\right)$ of order $m$, where the row vector
$\alpha$ is a probability vector whose $j$th entry is the probability that a
service begins in phase $j$ for $1\leq j\leq m$; and $T$ is a matrix of order
$m$ whose $\left(i,j\right)^{\text{th}}$ entry is denoted by $t_{i,j}$ with
$t_{i,i}<0$ for $1\leq i\leq m$, and $t_{i,j}\geq 0$ for $1\leq i,j\leq m$ and
$i\neq j$. Let $T^{0}=-Te\gvertneqq 0$, where $e$ is a column vector of ones
with a suitable dimension in the context. When a PH service time is in phase
$i$, the transition rate from phase $i$ to phase $j$ is $t_{i,j}$, the service
completion rate is $t_{i}^{0}$. At the same time, the mean service rate is
given by

$\mu=-\frac{1}{\alpha T^{-1}e}.$

Unless we state otherwise, we assume that all the random variables defined
above are independent, and that the system is operating at the stable region:
$\rho=\lambda/\mu<1$.

We introduce some useful notation. Let $n_{k}^{\left(i\right)}\left(t\right)$
be the number of queues with at least $k$ customers and the service time in
phase $i$ at time $t\geq 0$. Clearly, $0\leq
n_{k}^{\left(i\right)}\left(t\right)\leq n$ for $1\leq i\leq m$ and $0\leq
k\leq n$. We define

$s_{k}^{\left(i\right)}\left(t\right)=\frac{n_{k}^{\left(i\right)}\left(t\right)}{n},$

which is the fraction of queues with at least $k$ customers and the service
time in phase $i$. We write

$S_{0}\left(t\right)=\left(s_{0}\left(t\right)\right)$

and for $k\geq 1$,

$S_{k}\left(t\right)=\left(s_{k}^{\left(1\right)}\left(t\right),s_{k}^{\left(2\right)}\left(t\right),\ldots,s_{k}^{\left(m\right)}\left(t\right)\right),$
$S\left(t\right)=\left(S_{0}\left(t\right),S_{1}\left(t\right),S_{2}\left(t\right),\ldots\right).$

We now introduce Hadamard Product of two matrices $A=\left(a_{i,j}\right)$ and
$B=\left(b_{i,j}\right)$ as follows:

$A\odot B=\left(a_{i,j}b_{i,j}\right).$

Specifically, for $k\geq 2$ we have

$A^{\odot k}=\underset{k\text{ matrix }A}{\underbrace{A\odot A\odot\cdots\odot
A}}.$

Let $a=\left(a_{1},a_{2},a_{3},\ldots\right)$. We write

$a^{\odot\frac{1}{d}}=\left(a_{1}^{\frac{1}{d}},a_{2}^{\frac{1}{d}},a_{3}^{\frac{1}{d}},\ldots\right).$

Using a similar analysis to that in Equations (4) to (7), we can obtain the
following systems of differential vector equations for the fraction density
vector
$S\left(t\right)=\left(S_{0}\left(t\right),S_{1}\left(t\right),S_{2}\left(t\right),\ldots\right)$.

$S_{0}\left(t\right)=1,\text{ \ for }t\geq 0,$ (21)
$\frac{\mathtt{d}}{\text{d}t}S_{0}\left(t\right)=-\lambda S_{0}^{\odot
d}\left(t\right)+S_{1}\left(t\right)T^{0},$ (22)
$\frac{\mathtt{d}}{\text{d}t}S_{1}\left(t\right)=\lambda\alpha S_{0}^{\odot
d}\left(t\right)-\lambda S_{1}^{\odot
d}\left(t\right)+S_{1}\left(t\right)T+S_{2}\left(t\right)T^{0}\alpha,$ (23)

and for $k\geq 2$,

$\frac{\mathtt{d}}{\text{d}t}S_{k}\left(t\right)=\lambda S_{k-1}^{\odot
d}\left(t\right)-\lambda S_{k}^{\odot
d}\left(t\right)+S_{k}\left(t\right)T+S_{k+1}\left(t\right)T^{0}\alpha.$ (24)

If $\pi=\left(\pi_{0},\pi_{1},\pi_{2},\ldots\right)$ is a fixed point of the
fraction density vector $S\left(t\right)$, then the system of differential
vector equations (21) to (24) can be simplified as

$\pi_{0}=1$ (25) $-\lambda\pi_{0}^{\odot d}+\pi_{1}T^{0}=0,$ (26)
$\lambda\alpha\pi_{0}^{\odot d}-\lambda\pi_{1}^{\odot
d}+\pi_{1}T+\pi_{2}T^{0}\alpha=0,$ (27)

and for $k\geq 2$,

$\lambda\pi_{k-1}^{\odot d}-\lambda\pi_{k}^{\odot
d}+\pi_{k}T+\pi_{k+1}T^{0}\alpha=0.$ (28)

In what follows we provide three methods to solve the system of nonlinear
equations (25) to (28), and give three different doubly exponential solutions
to the fixed point.

###  5.1 The first method

The first method is based on Theorem 1 given in this paper. For the PH service
time distribution

$\overline{G}\left(x\right)=\alpha\exp\left\\{Tx\right\\}e$

Let
$\theta=\int_{0}^{+\infty}\left[\mu\overline{G}\left(x\right)\right]^{d}$d$x$,
and we assume that $0<\theta<+\infty$. Then the fixed point
$\pi=\left(\pi_{0},\pi_{1}\left(x\right),\pi_{2}\left(x\right),\ldots\right)$
is given by

$\pi_{0}=1,$

and for $k\geq 1$

$\pi_{k}\left(x\right)=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\mu\overline{G}\left(x\right).$
(29)

###  5.2 The second method

The second method is proposed in Li, Wang and Liu [11], and the key parameter
$\theta$ is based on the stationary probability vector $\omega$ of the
irreducible Markov chain $T+T^{0}\alpha$, that is, $\theta=\omega^{\odot d}e$.

It follows from Equation (26) that

$\pi_{1}T^{0}=\lambda.$

Note that

$\omega T^{0}=\mu,$ $\frac{\lambda}{\mu}\omega T^{0}=\lambda.$ (30)

Thus, we obtain

$\pi_{1}=\frac{\lambda}{\mu}\omega=\rho\cdot\,\omega.$

Based on the fact that $\pi_{0}=1$ and $\pi_{1}=\rho\cdot\omega$, it follows
from Equation (27) that

$\lambda\alpha-\lambda\rho^{d}\cdot\omega^{\odot d}+\rho\cdot\omega
T+\pi_{2}T^{0}\alpha=0,$

which leads to

$\lambda-\lambda\rho^{d}\cdot\omega^{\odot d}e+\rho\cdot\omega
Te+\pi_{2}T^{0}=0.$

Note that $\omega Te=-\mu$ and $\rho=\lambda/\mu$, we obtain

$\pi_{2}T^{0}=\lambda\rho^{d}\omega^{\odot d}e.$

Let $\theta=\omega^{\odot d}e$. Then it is easy to see that
$\theta\in\left(0,1\right)$, and

$\pi_{2}T^{0}=\lambda\theta\rho^{d}.$

Using a similar analysis on Equation (30), we have

$\pi_{2}=\frac{\lambda\theta\rho^{d}}{\mu}\omega=\theta\rho^{d+1}\cdot\omega.$

Based on $\pi_{1}=\rho\,\omega$ and $\pi_{2}=\theta\rho^{d+1}\cdot\omega$, it
follows from Equation (28) that for $k=2$,

$\lambda\rho^{d}\cdot\omega^{\odot
d}-\lambda\theta^{d}\rho^{d^{2}+d}\cdot\omega^{\odot
d}+\theta\rho^{d+1}\cdot\omega T+\pi_{3}T^{0}\alpha=0,$

which leads to

$\lambda\theta\rho^{d}-\lambda\theta^{d+1}\rho^{d^{2}+d}+\theta\rho^{d+1}\cdot\omega
Te+\pi_{3}T^{0}=0,$

thus we obtain

$\pi_{3}T^{0}=\lambda\theta^{d+1}\rho^{d^{2}+d}.$

Using a similar analysis on Equation (30), we have

$\pi_{3}=\frac{\lambda\theta^{d+1}\rho^{d^{2}+d}}{\mu}\omega=\theta^{d+1}\rho^{d^{2}+d+1}\cdot\omega.$

Now, we assume that
$\pi_{k}=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\omega$
is correct for the cases with $l=k$. Then it follows from Equation (28) that
for $l=k+1$, we have

$\displaystyle\lambda$
$\displaystyle\theta^{d^{k-2}+d^{k-3}+\cdots+d}\rho^{d^{k-1}+d^{k-2}+\cdots+d}\cdot\omega^{\odot
d}-\lambda\theta^{d^{k-1}+d^{k-2}+\cdots+d}\rho^{d^{k}+d^{k-1}+\cdots+d}\cdot\omega^{\odot
d}$
$\displaystyle+\theta^{d^{k-2}+d^{k-3}+\cdots+1}\rho^{d^{k-1}+d^{k-2}+\cdots+1}\cdot\omega
T+\pi_{k+1}T^{0}\alpha=0,$

which leads to

$\displaystyle\lambda$
$\displaystyle\theta^{d^{k-2}+d^{k-3}+\cdots+d+1}\rho^{d^{k-1}+d^{k-2}+\cdots+d}-\lambda\theta^{d^{k-1}+d^{k-2}+\cdots+d+1}\rho^{d^{k}+d^{k-1}+\cdots+d}$
$\displaystyle+\theta^{d^{k-2}+d^{k-3}+\cdots+1}\rho^{d^{k-1}+d^{k-2}+\cdots+1}\cdot\omega
Te+\pi_{k+1}T^{0}=0,$

thus we obtain

$\pi_{k+1}T^{0}=\lambda\theta^{d^{k-1}+d^{k-2}+\cdots+d+1}\rho^{d^{k}+d^{k-1}+\cdots+d}.$

By a similar analysis to (30), we have

$\displaystyle\pi_{k+1}$
$\displaystyle=\frac{\lambda\theta^{d^{k-1}+d^{k-2}+\cdots+d+1}\rho^{d^{k}+d^{k-1}+\cdots+d}}{\mu}\omega$
$\displaystyle=\theta^{d^{k-1}+d^{k-2}+\cdots+d+1}\rho^{d^{k}+d^{k-1}+\cdots+d+1}\cdot\omega.$

Therefore, by induction the fixed point
$\pi=\left(\pi_{0},\pi_{1},\pi_{2},\ldots\right)$ is given by

$\pi_{0}=1,$

and for $k\geq 1$

$\pi_{k}=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\omega.$
(31)

###  5.3 The third method

The third method is based on the matrix computation for the system of
nonlinear equations (25) to (28), and shows that the key parameter $\theta$ is
based on the initial probability vector $\alpha$ in the PH service time
distribution, that is, $\theta=1/\alpha^{\odot\frac{1}{d}}e.$

It follows from (25) to (28) that

$\displaystyle\left(\pi_{1}^{\odot d},\pi_{2}^{\odot d},\pi_{3}^{\odot
d},\ldots\right)\left(\begin{array}[c]{ccccc}-\lambda&\lambda&&&\\\
&-\lambda&\lambda&&\\\ &&-\lambda&\lambda&\\\
&&&\ddots&\ddots\end{array}\right)+\left(\pi_{1},\pi_{2},\pi_{3},\ldots\right)\left(\begin{array}[c]{cccc}T&&&\\\
T^{0}\alpha&T&&\\\ &T^{0}\alpha&T&\\\ &&\ddots&\ddots\end{array}\right)$
$\displaystyle=-\left(\lambda\alpha,0,0,\ldots\right),$

which leads to

$\displaystyle\left(\pi_{1},\pi_{2},\pi_{3},\ldots\right)=$
$\displaystyle\left(\pi_{1}^{\odot d},\pi_{2}^{\odot d},\pi_{3}^{\odot
d},\ldots\right)\left(\begin{array}[c]{ccccc}R&V&&&\\\ &R&V&&\\\ &&R&V&\\\
&&&\ddots&\ddots\end{array}\right)$
$\displaystyle+\left(\lambda\alpha\left(-T\right)^{-1},0,0,\ldots\right),$

where

$V=\lambda\left(-T\right)^{-1}$

and

$R=\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}.$

Thus we obtain

$\pi_{1}=\lambda\alpha\left(-T\right)^{-1}+\pi_{1}^{\odot
d}\left[\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}\right]$ (32)

and for $k\geq 2$

$\pi_{k}=\pi_{k-1}^{\odot
d}\left[\lambda\left(-T\right)^{-1}\right]+\pi_{k}^{\odot
d}\left[\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}\right].$ (33)

To omit the term $\pi_{k}^{\odot
d}\left[\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}\right]$ for $k\geq
1$, we assume that $\left\\{\pi_{k},k\geq 1\right\\}$ has the following
expression

$\pi_{k}=r\left(k\right)\alpha^{\odot\frac{1}{d}}.$

In this case, we have

$\pi_{k}^{\odot
d}\left[\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}\right]=r^{d}\left(k\right)\alpha\left[\lambda\left(-I+e\alpha\right)\left(-T\right)^{-1}\right]=0,$

thus it follows from (32) and (33) that

$\pi_{1}=\lambda\alpha\left(-T\right)^{-1}$ (34)

and for $k\geq 2$

$\pi_{k}=\pi_{k-1}^{\odot d}\left[\lambda\left(-T\right)^{-1}\right].$ (35)

It follows from (34) that

$r\left(1\right)\alpha^{\odot\frac{1}{d}}=\lambda\alpha\left(-T\right)^{-1},$

which follows that

$r\left(1\right)=\theta\rho,$

where

$\theta=\frac{1}{\alpha^{\odot\frac{1}{d}}e}.$

It follows from (35) that

$r\left(k\right)\alpha^{\odot\frac{1}{d}}=r^{d}\left(k-1\right)\alpha\left[\lambda\left(-T\right)^{-1}\right],$

which follows that

$r\left(k\right)=r^{d}\left(k-1\right)\theta\rho=\left(\theta\rho\right)^{\frac{d^{k}-1}{d-1}}.$

Therefore, we can obtain

$\pi_{0}=1$

and for $k\geq 1$

$\pi_{k}=\left(\theta\rho\right)^{\frac{d^{k}-1}{d-1}}\cdot\alpha^{\odot\frac{1}{d}}.$
(36)

###  5.4 Non-uniqueness

Based on the above three methods, we can summarize the key parameter and the
doubly exponential solution to the fixed point in the following table.

Table 4: Comparison for the three methods | Key Parameter | Fixed point
---|---|---
Method 1 | $\theta=\frac{\int_{0}^{+\infty}\left[\alpha\exp\left\\{Tx\right\\}e\right]^{d}\text{d}x}{\left[-\alpha T^{-1}e\right]^{d}}$ | $\pi_{k}=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\mu\overline{G}\left(x\right)$
Method 2 | $\theta=\omega^{\odot d}e$ | $\pi_{k}=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\omega$
Method 3 | $\theta=1/\alpha^{\odot\frac{1}{d}}e$ | $\pi_{k}=\left(\theta\rho\right)^{\frac{d^{k}-1}{d-1}}\cdot\alpha^{\odot\frac{1}{d}}$

When the PH service time is an $m$-order Erlang distribution with the
irreducible representation $(\alpha,T)$, where

$\alpha=\left(1,0,\ldots,0\right)$

and

$T=\left(\begin{array}[c]{ccccc}-\eta&\eta&&&\\\ &-\eta&\eta&&\\\
&&\ddots&\ddots&\\\ &&&-\eta&\eta\\\ &&&&-\eta\end{array}\right),\text{ \ \
}T^{0}=\left(\begin{array}[c]{c}0\\\ 0\\\ \vdots\\\ 0\\\
\eta\end{array}\right).$

We have

$\alpha^{\odot\frac{1}{d}}=\left(1,0,\ldots,0\right)$

and

$\theta=\frac{1}{\alpha^{\odot\frac{1}{d}}e}=1.$

Thus the doubly exponential solution by the third method is given by

$\pi_{k}=\rho^{\frac{d^{k}-1}{d-1}}\cdot\left(1,0,\ldots,0\right),\text{ \
}k\geq 1.$ (37)

It is clear that

$T+T^{0}\alpha=\left(\begin{array}[c]{ccccc}-\eta&\eta&&&\\\ &-\eta&\eta&&\\\
&&\ddots&\ddots&\\\ &&&-\eta&\eta\\\ \eta&&&&-\eta\end{array}\right),$

which leads to the stationary probability vector of the Markov chain
$T+T^{0}\alpha$ as follows:

$\omega=\left(\frac{1}{m},\frac{1}{m},\ldots,\frac{1}{m}\right),$ $\mu=\omega
T^{0}=\frac{\eta}{m},$ $\rho=\frac{\lambda}{\mu}=\frac{m\lambda}{\eta}$

and

$\theta=\omega^{\odot d}e=m\left(\frac{1}{m}\right)^{d}=m^{1-d}.$

Thus the doubly exponential solution by the second method is given by

$\displaystyle\pi_{k}$
$\displaystyle=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\left(\frac{1}{m},\frac{1}{m},\ldots,\frac{1}{m}\right)$
$\displaystyle=\rho^{\frac{d^{k}-1}{d-1}}\cdot\left(m^{-d^{k-1}},m^{-d^{k-1}},\ldots,m^{-d^{k-1}}\right),\text{
\ }k\geq 1.$ (38)

It is clear that the three doubly exponential solutions (29), (37) and (38)
are different for $m\geq 2$.

###### Remark 4

For the supermarket model with Poisson arrivals and PH service times, we have
obtained three different doubly exponential solutions to the fixed point. It
is interesting but difficult how to be able to find another new doubly
exponential solution. We believe that it is an open problem how to give all
the doubly exponential solutions to the fixed point for a more general
supermarket model including the case with MAP arrivals, PH service times or
general service times.

##  6 Exponential convergence to the fixed point

In this section, we study the exponential convergence of the current location
$S\left(t,x\right)$ of the supermarket model to its fixed point
$\pi\left(x\right)$ for $t\geq 0$ and $x\geq 0$. Not only does the exponential
convergence indicates the existence of the fixed point, but it also explains
such a convergent process is very fast.

For the supermarket model, the initial point $S\left(0,x\right)$ can affect
the current location $S\left(t,x\right)$ for each $t>0$, since the service
process in the supermarket model is under a unified structure. Here, we
provide notation for comparison of two vectors. Let
$a=\left(a_{1},a_{2},a_{3},\ldots\right)$ and
$b=\left(b_{1},b_{2},b_{3},\ldots\right)$. We write $a\prec b$ if
$a_{k}<b_{k}$ for all $k\geq 1$; $a\preceq b$ if $a_{k}\leq b_{k}$ for all
$k\geq 1$.

Now, we can obtain the following useful proposition whose proof is clear from
a sample path analysis and is omitted here.

###### Proposition 2

If $S\left(0,x\right)\preceq\widetilde{S}\left(0,x\right)$, then
$S\left(t,x\right)\preceq\widetilde{S}\left(t,x\right)$.

Based on Proposition 2, the following theorem shows that the fixed point
$\pi\left(x\right)$ is an upper bound of the current location
$S\left(t,x\right)$ for all $t\geq 0$ and $x\geq 0$.

###### Theorem 4

For the supermarket model, if there exists some $k$ such that
$S_{k}\left(0,x\right)=0$, then the sequence
$\left\\{S_{k}\left(t,x\right)\right\\}$ has a upper bound sequence which
decreases doubly exponentially for all $t\geq 0$ and $x\geq 0$, that is,
$S\left(t,x\right)\preceq\pi\left(x\right)$ for all $t\geq 0$ and $x\geq 0$.

Proof Let

$\widetilde{S}_{0}\left(0\right)=\pi_{0}$
$\widetilde{S}_{k}\left(0,x\right)=\pi_{k}\left(x\right),\text{ \ }k\geq 1.$

Then
$\widetilde{S}_{0}\left(t\right)=\widetilde{S}_{0}\left(0\right)=\pi_{0}$, and
for each $k\geq 1$,
$\widetilde{S}_{k}\left(t,x\right)=\widetilde{S}_{k}\left(0,x\right)=\pi_{k}\left(x\right)$
for all $t\geq 0$ and $x\geq 0$, since $\widetilde{S}\left(0\right)$ is a
fixed point in the supermarket model. If $S_{k}\left(0,x\right)=0$ for some
$k$, then $S_{k}\left(0,x\right)<\widetilde{S}_{k}\left(0,x\right)$ and
$S_{j}\left(0,x\right)\leq\widetilde{S}_{j}\left(0,x\right)$ for all $j\geq 1$
and $j\neq k$, thus $S\left(0,x\right)\preceq\widetilde{S}\left(0,x\right)$.
It is easy to see from Proposition 2 that
$S_{k}\left(t,x\right)\leq\widetilde{S}_{k}\left(t,x\right)=\pi_{k}\left(x\right)$
for all $k\geq 1$, $t\geq 0$ and $x\geq 0$. Thus we obtain that for all $k\geq
1$, $t\geq 0$ and $x\geq 0$

$S_{k}\left(t,x\right)\leq\pi_{k}\left(x\right)=\theta^{\frac{d^{k-1}-1}{d-1}}\rho^{\frac{d^{k}-1}{d-1}}\cdot\left[\mu\overline{G}\left(x\right)\right].$

This completes the proof.

To show the exponential convergence, we use Theorem 4 to define a potential
function (or Lyapunov function) $\Phi\left(t\right)$ as follows:

$\Phi\left(t\right)=\sum_{k=1}^{\infty}w_{k}\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x,$

where $\left\\{w_{k}\right\\}$ is a positive scalar sequence with
$w_{k}>w_{k-1}\geq w_{1}=1$ for $k\geq 2$. Note that
$\pi_{0}=S_{0}\left(t\right)=1$. It is easy to see from Proposition 2 that
$\Phi\left(t\right)\geq 0$ for all $t\geq 0$.

When
$\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]$d$x>0$
for $k\geq 1$, we write

$\frac{\int_{0}^{+\infty}S_{k}^{d}\left(t,x\right)\text{d}x}{\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x}=c_{k}\left(t\right)$

and

$\frac{\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(t,x\right)\text{d}x}{\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x}=d_{k}\left(t\right).$

The following lemma provide a method to determine the positive scalar sequence
$\left\\{w_{k}\right\\}$ with $w_{k}>w_{k-1}\geq w_{1}=1$ for $k\geq 2$. This
proof is easy by means of some simple computation.

###### Lemma 2

If $\delta$ is a positive constant,

$w_{1}=1,$ $\lambda\left(w_{1}-w_{2}\right)c_{1}\left(t\right)=-\delta w_{1}$

and for $k\geq 2$

$\lambda\left(w_{k}-w_{k+1}\right)c_{k}\left(t\right)+\left(w_{k}-w_{k-1}\right)d_{k}\left(t\right)=-\delta
w_{k},$

then

$w_{2}=1+\frac{\delta}{\lambda c_{1}\left(t\right)},$

and for $k\geq 3$

$w_{k}=w_{k-1}+\frac{\delta
w_{k-1}+\left(w_{k-1}-w_{k-2}\right)d_{k-1}\left(t\right)}{\lambda
c_{k-1}\left(t\right)}.$

The following theorem measures the distance $\Phi\left(t\right)$ of the
current location $S\left(t,x\right)$ for $t\geq 0$ and the fixed point
$\pi\left(x\right)$ for $x\geq 0$, and illustrates that the distance
$\Phi\left(t\right)$ to the fixed point from the current location is very go
to zero with exponential convergence. Hence, it shows that from any suitable
starting point, the supermarket model can be quickly close to the fixed point,
that is, there always exists a fixed point in the supermarket model.

###### Theorem 5

For $t\geq 0$,

$\Phi\left(t\right)\leq c_{0}e^{-\delta t},$

where $c_{0}$ and $\delta$ are two positive constants, and they possibly
depend on time $t\geq 0$. In this case, the potential function
$\Phi\left(t\right)$ is exponentially convergent.

Proof Note that

$\Phi\left(t\right)=\sum_{k=1}^{\infty}w_{k}\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x,$

we have

$\displaystyle\frac{\text{d}}{\text{d}t}\Phi\left(t\right)$
$\displaystyle=\frac{\text{d}}{\text{d}t}\sum_{k=1}^{\infty}w_{k}\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x$
$\displaystyle=-\sum_{k=1}^{\infty}w_{k}\frac{\text{d}}{\text{d}t}\int_{0}^{+\infty}S_{k}\left(t,x\right)\text{d}x$

by means of the Dominated Convergence Theorem. It follows from (4) to (7) that

$\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(t,x\right)\text{d}x=\lambda,$
(39)

and using (39) we obtain

$\displaystyle\frac{d}{dt}\Phi\left(t\right)=$
$\displaystyle-\sum_{k=1}^{\infty}w_{k}\frac{\mathtt{d}}{\text{d}t}\int_{0}^{+\infty}S_{k}\left(t,x\right)\text{d}x$
$\displaystyle=$ $\displaystyle-
w_{1}[\lambda-\lambda\int_{0}^{+\infty}S_{1}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{2}\left(t,x\right)\text{d}x]$
$\displaystyle-\sum_{k=2}^{\infty}w_{k}[\lambda\int_{0}^{+\infty}S_{k-1}^{d}\left(t,x\right)\text{d}x-\lambda\int_{0}^{+\infty}S_{k}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{k+1}\left(t,x\right)\text{d}x]$
$\displaystyle=$ $\displaystyle-
w_{1}[-\lambda\int_{0}^{+\infty}S_{1}^{d}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{2}\left(t,x\right)\text{d}x]$
$\displaystyle-\sum_{k=2}^{\infty}w_{k}[\lambda\int_{0}^{+\infty}S_{k-1}^{d}\left(t,x\right)\text{d}x-\lambda\int_{0}^{+\infty}S_{k}^{d}\left(t,x\right)\text{d}x$
$\displaystyle-\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(t,x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{k+1}\left(t,x\right)\text{d}x],$

which follows

$\displaystyle\frac{d}{dt}\Phi\left(t\right)=$
$\displaystyle\lambda\left(w_{1}-w_{2}\right)c_{1}\left(t\right)\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x$
$\displaystyle+\sum_{k=2}^{\infty}\left[\lambda\left(w_{k}-w_{k+1}\right)c_{k}\left(t\right)+\left(w_{k}-w_{k-1}\right)d_{k}\left(t\right)\right]$
$\displaystyle\cdot\int_{0}^{+\infty}\left[\pi_{k}\left(x\right)-S_{k}\left(t,x\right)\right]\text{d}x.$

Using Lemma 2 we can easily choose a parameter $\delta>0$ and a suitable
positive scalar sequence $\left\\{w_{k}\right\\}$ with $w_{k}>w_{k-1}\geq
w_{1}=1$ for $k\geq 2$ such that

$\lambda\left(w_{1}-w_{2}\right)c_{1}\left(t\right)\leq-\delta w_{1}$

and for $k\geq 2$

$\left(w_{k}-w_{k-1}\right)d_{k}\left(t\right)-\lambda\left(w_{k+1}-w_{k}\right)c_{k}\left(t\right)\leq-\delta
w_{k},$

thus we can obtain

$\frac{d}{dt}\Phi\left(t\right)\leq-\delta\Phi\left(t\right),$

which leads to

$\Phi\left(t\right)\leq c_{0}e^{-\delta t}.$

This completes the proof.

###### Remark 5

We have provided an algorithm for computing the positive scalar sequence
$\left\\{w_{k}\right\\}$ with $1=w_{1}\leq w_{k-1}<w_{k}$ for $k\geq 2$ as
follows:

Step one:

$w_{1}=1.$

Step two:

$w_{2}=1+\frac{\delta}{\lambda c_{1}\left(t\right)}.$

Step three: for $k\geq 2$

$w_{k}=w_{k-1}+\frac{\delta
w_{k-1}+\left(w_{k-1}-w_{k-2}\right)d_{k-1}\left(t\right)}{\lambda
c_{k-1}\left(t\right)}.$

This illustrates that $w_{k}$ is a function of time $t$. Note that
$\lambda,\delta,c_{k}\left(t\right),d_{l}\left(t\right)>0$, it is clear that
for $k\geq 2$

$1=w_{1}\leq w_{k-1}<w_{k}.$

##  7 Lipschitz Condition

In this section, we apply the Kurtz Theorem to study the supermarket model
with general service times, and analyze the Lipschitz condition with respect
to general service times.

The supermarket model can be analyzed by a density dependent jump Markov
process, where the density dependent jump Markov process is a Markov process
with a single parameter $n$ which corresponds to the population size. Kurtz’s
work provides a basis for the density dependent jump Markov processes in order
to relate the infinite-size system of differential equations to the
corresponding finite-size system of differential equations. Readers may refer
to Kurtz [9] for more details.

In the supermarket model, the states of density dependent jump Markov process
can be normalized and interpreted as measuring population densities, so that
the transition rates depend only on these densities. Hence, the infinite-size
system of differential equations can be regarded as the limiting model of the
corresponding finite-size system of differential equations as the population
size grows arbitrarily large. When the population size is $n$, we write

$E_{n}=\left\\{k:k=0,1,\ldots,n\right\\}.$

For $k\geq 1$ and $x\geq 0$, we write

$s_{k}^{\left(n\right)}\left(x\right)=\left(\frac{k}{n},x\right),$

where $x$ is the residual service time of each server, and

$s_{k}^{\left(n\right)}=\int_{0}^{+\infty}s_{k}^{\left(n\right)}\left(x\right)\text{d}x.$

Let

$S_{0}=\lim_{n\rightarrow\infty}s_{0}^{\left(n\right)}$

and for $k\geq 1$

$S_{k}=\lim_{n\rightarrow\infty}\int_{0}^{+\infty}s_{k}^{\left(n\right)}\left(x\right)\text{d}x.$

Let $\left\\{\widehat{X}_{n}\left(t\right):t\geq 0\right\\}$ be a density
dependent jump Markov process on the state space $E_{n}$ whose transition
rates are given by

$q_{k,k+l}^{\left(n\right)}=n\beta_{l}\left(\frac{k}{n}\right)=n\beta_{l}\left(s_{k}^{\left(n\right)}\right).$

In this supermarket model, $\widehat{X}_{n}\left(t\right)$ is the unscaled
process which records the number of servers with at least $k$ customers for
$0\leq k\leq n$.

Let $a$ and $b$ denote an arrival and a service completion, respectively.
Hence taking $l=a$ or $b$ for $a>b>0$, we write

$\beta_{a}\left(s_{0}^{\left(n\right)}\right)=-\lambda,$
$\beta_{b}\left(s_{0}^{\left(n\right)}\right)=\int_{0}^{+\infty}\mu\left(x\right)s_{1}^{\left(n\right)}\left(x\right)\text{d}x;$
$\beta_{a}\left(s_{1}^{\left(n\right)}\right)=\lambda-\lambda\int_{0}^{+\infty}\left[s_{1}^{\left(n\right)}\left(x\right)\right]^{d}\text{d}x,$
$\beta_{b}\left(s_{1}^{\left(n\right)}\right)=-\int_{0}^{+\infty}\mu\left(x\right)s_{1}^{\left(n\right)}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)s_{2}^{\left(n\right)}\left(x\right)\text{d}x;$

and for $n\geq k\geq 2,$

$\beta_{a}\left(s_{k}^{\left(n\right)}\right)=\lambda\int_{0}^{+\infty}\left[s_{k-1}^{\left(n\right)}\left(x\right)\right]^{d}\text{d}x-\lambda\int_{0}^{+\infty}\left[s_{k}^{\left(n\right)}\left(x\right)\right]^{d}\text{d}x,$
$\beta_{b}\left(s_{k}^{\left(n\right)}\right)=-\int_{0}^{+\infty}\mu\left(x\right)s_{k}^{\left(n\right)}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)s_{k+1}^{\left(n\right)}\left(x\right)\text{d}x.$

Using Chapter 7 in Kurtz [9] or Subsection 3.4.1 in Mitzenmacher [20], the
Markov process $\left\\{\widehat{X}_{n}\left(t\right):t\geq 0\right\\}$ with
transition rates $q_{k,k+l}^{\left(n\right)}$ is given by

$\widehat{X}_{n}\left(t\right)=\widehat{X}_{n}\left(0\right)+\sum_{l=a,b}lY_{l}\left(n\int_{0}^{t}\beta_{l}\left(\frac{\widehat{X}_{n}\left(u\right)}{n}\right)\text{d}u\right),$
(40)

where $Y_{l}\left(x\right)$ for $l=a$ and $b$ are two independent standard
Poisson processes. Clearly, the jump Markov process by Equation (40) at time
$t$ is determined by the starting point and the transition rates which are
integrated over its history.

Let

$F\left(y\right)=a\beta_{a}\left(y\right)+b\beta_{b}\left(y\right).$ (41)

Taking $X_{n}\left(t\right)=n^{-1}\widehat{X}_{n}\left(t\right)$ which is an
appropriate scaled process, we have

$X_{n}\left(t\right)=X_{n}\left(0\right)+\sum_{l=a,b}ln^{-1}\widehat{Y}_{l}\left(n\int_{0}^{t}\beta_{l}\left(X_{n}\left(u\right)\right)\text{d}u\right)+\int_{0}^{t}F\left(X_{n}\left(u\right)\right)\text{d}u,$
(42)

where $\widehat{Y}_{l}\left(y\right)=Y_{l}\left(y\right)-y$ is a Poisson
process centered at its expectation. Note that in (42), the function
$F\left(y\right)$ given in (41) is for $y=s_{k}^{\left(n\right)},0\leq k\leq
n$.

Taking $X\left(t\right)=\lim_{n\rightarrow\infty}X_{n}\left(t\right)$ and
$x_{0}=\lim_{n\rightarrow\infty}X_{n}\left(0\right)$, we obtain

$X\left(t\right)=x_{0}+\int_{0}^{t}F\left(X\left(u\right)\right)\text{d}u,\text{
\ }t\geq 0,$ (43)

due to the fact that

$\lim_{n\rightarrow\infty}\frac{1}{n}\widehat{Y}_{l}\left(n\int_{0}^{t}\beta_{l}\left(X_{n}\left(u\right)\right)\text{d}u\right)=0$

by means of the law of large numbers. Note that in (43), the function
$F\left(y\right)$ given in (41) is for $y=S_{k},k\geq 1$. In the supermarket
model, the deterministic and continuous process $\left\\{X\left(t\right),t\geq
0\right\\}$ is described by the infinite-size system of integral-differential
equations (4) to (7), or simply in the below

$\frac{d}{dt}X\left(t\right)=F\left(X\left(t\right)\right)$ (44)

with the initial condition

$X\left(0\right)=x_{0}.$ (45)

Now, we consider the uniqueness of the limiting deterministic process
$\left\\{X\left(t\right),t\geq 0\right\\}$ with (44) to (45), or the
uniqueness of solution to the infinite-size system of integral-differential
equations (4) to (7). To that end, a sufficient condition is Lipschitz, that
is, for some constant $M>0,$

$|F\left(y\right)-F\left(z\right)|\leq M|y-z|.$

In general, the Lipschitz condition is standard and sufficient for the
uniqueness of solution to the finite-size system of differential equations;
while for the countable infinite-size case, readers may refer to Theorem 3.2
in Deimling [4] and Subsection 3.4.1 in Mitzenmacher [20] for some
generalization.

To check the Lipschitz condition, as $n\rightarrow\infty$ we have

$\beta_{a}\left(S_{0}\right)=-\lambda,$
$\beta_{b}\left(S_{0}\right)=\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(x\right)\text{d}x;$
$\beta_{a}\left(S_{1}\right)=\lambda-\lambda\int_{0}^{+\infty}\left[S_{1}\left(x\right)\right]^{d}\text{d}x,$
$\beta_{b}\left(S_{1}\right)=-\int_{0}^{+\infty}\mu\left(x\right)S_{1}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{2}\left(x\right)\text{d}x;$

and for $k\geq 2,$

$\beta_{a}\left(S_{k}\right)=\lambda\int_{0}^{+\infty}\left[S_{k-1}\left(x\right)\right]^{d}\text{d}x-\lambda\int_{0}^{+\infty}\left[S_{k}\left(x\right)\right]^{d}\text{d}x,$
$\beta_{b}\left(S_{k}\right)=-\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(x\right)\text{d}x+\int_{0}^{+\infty}\mu\left(x\right)S_{k+1}\left(x\right)\text{d}x.$

Let

$\zeta_{k}=\frac{\int_{0}^{+\infty}\left[S_{k}\left(x\right)\right]^{d}\text{d}x}{\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x}$

and

$\eta_{k}=\frac{\int_{0}^{+\infty}\mu\left(x\right)S_{k}\left(x\right)\text{d}x}{\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x}.$

Then $\zeta_{k},\eta_{k}>0$ for $k\geq 1$.

The following theorem shows that the supermarket model with general service
times satisfies the Lipschitz condition for the infinite-size system of
integral-differential equations (4) to (7).

###### Theorem 6

The supermarket model with general service times satisfies the Lipschitz
condition.

Proof Let

$\Omega=\left\\{S_{k}:k\geq 0\right\\}.$

For two arbitrary entries $y,z\in\Omega$, we have

$|F\left(y\right)-F\left(z\right)|\leq
a|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|+b|\beta_{b}\left(y\right)-\beta_{b}\left(z\right)|.$

Now, we analyze the following four cases for the function
$\beta_{a}\left(y\right)$, while the function $\beta_{b}\left(y\right)$ can be
analyzed similarly.

Case one: $y=S_{0},z=S_{1}$. In this case, we have

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle=|-\lambda-\lambda+\lambda\int_{0}^{+\infty}\left[S_{1}\left(x\right)\right]^{d}\text{d}x|$
$\displaystyle=\lambda|2-\zeta_{1}\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x|$
$\displaystyle=\lambda\left[2-\zeta_{1}\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x\right],$

since $0<\zeta_{1},\int_{0}^{+\infty}S_{1}\left(x\right)$d$x<1$. Taking

$M_{a}\left(0,1\right)\geq\frac{2-\zeta_{1}\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x}{2-\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x},$

it is clear that

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle\leq
M_{a}\left(0,1\right)\lambda\left[2-\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x\right]$
$\displaystyle=M_{a}\left(0,1\right)|y-z|.$

Case two: $y=S_{0},z=S_{k}$ for $k\geq 2$. In this case, we have

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle=|-\lambda-\lambda\int_{0}^{+\infty}\left[S_{k-1}\left(x\right)\right]^{d}\text{d}x+\lambda\int_{0}^{+\infty}\left[S_{k}\left(x\right)\right]^{d}\text{d}x|$
$\displaystyle=\lambda|-1-\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|$
$\displaystyle=\lambda\left[1+\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x-\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x\right]$

due to that $0<\zeta_{k},\int_{0}^{+\infty}S_{k}\left(x\right)$d$x<1$. Let

$M_{a}\left(0,k\right)\geq\frac{1+\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x-\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x}{1+\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x-\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x.}$

Then

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle\leq
M_{a}\left(0,k\right)\lambda\left[1+\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x-\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x\right]$
$\displaystyle=M_{a}\left(0,k\right)|y-z|.$

Case three: $y=S_{1},z=S_{k}$ for $k\geq 2$. In this case, we have

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle=|\lambda-\lambda\int_{0}^{+\infty}\left[S_{1}\left(x\right)\right]^{d}\text{d}x-\lambda\int_{0}^{+\infty}\left[S_{k-1}\left(x\right)\right]^{d}\text{d}x+\lambda\int_{0}^{+\infty}\left[S_{k}\left(x\right)\right]^{d}\text{d}x|$
$\displaystyle=\lambda|1-\zeta_{1}\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x-\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|$

Let

$M_{a}\left(1,k\right)\geq\frac{|1-\zeta_{1}\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x-\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|}{|1-\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x-\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|}.$

Then

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|$
$\displaystyle\leq
M_{a}\left(1,k\right)\lambda|1-\int_{0}^{+\infty}S_{1}\left(x\right)\text{d}x-\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|$
$\displaystyle=M_{a}\left(1,k\right)|y-z|.$

Case four: $y=S_{l},z=S_{k}$ for $k>l\geq 2$. In this case, we have

$\displaystyle|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|=$
$\displaystyle|\lambda\int_{0}^{+\infty}\left[S_{l-1}\left(x\right)\right]^{d}\text{d}x-\lambda\int_{0}^{+\infty}\left[S_{l}\left(x\right)\right]^{d}\text{d}x$
$\displaystyle-\lambda\int_{0}^{+\infty}\left[S_{k-1}\left(x\right)\right]^{d}\text{d}x+\lambda\int_{0}^{+\infty}\left[S_{k}\left(x\right)\right]^{d}\text{d}x|$
$\displaystyle=$
$\displaystyle\lambda|\zeta_{l-1}\int_{0}^{+\infty}S_{l-1}\left(x\right)\text{d}x+\zeta_{l}\int_{0}^{+\infty}S_{l}\left(x\right)\text{d}x$
$\displaystyle-\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|.$

Let

$M_{a}\left(l,k\right)\geq\frac{|\zeta_{l-1}\int_{0}^{+\infty}S_{l-1}\left(x\right)\text{d}x+\zeta_{l}\int_{0}^{+\infty}S_{l}\left(x\right)\text{d}x-\zeta_{k-1}\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\zeta_{k}\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|}{|\int_{0}^{+\infty}S_{l-1}\left(x\right)\text{d}x+\int_{0}^{+\infty}S_{l}\left(x\right)\text{d}x-\int_{0}^{+\infty}S_{k-1}\left(x\right)\text{d}x+\int_{0}^{+\infty}S_{k}\left(x\right)\text{d}x|}.$

Then

$|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|\leq
M_{a}\left(l,k\right)|y-z|.$

Based on the above four cases, taking

$M_{a}=\max\left\\{M_{a}\left(l,k\right):k>l\geq 0\right\\}$

we obtain that for two arbitrary entries $y,z\in\Omega,$

$|\beta_{a}\left(y\right)-\beta_{a}\left(z\right)|\leq M_{a}|y-z|.$

Similarly, we can choose a positive number $M_{b}$ such that for two arbitrary
entries $y,z\in\Omega,$

$|\beta_{b}\left(y\right)-\beta_{b}\left(z\right)|\leq M_{b}|y-z|.$

Let $M=\max\left\\{aM_{a},bM_{b}\right\\}$. Then for two arbitrary entries
$y,z\in\Omega,$

$|F\left(y\right)-F\left(z\right)|\leq M|y-z|.$

This completes the proof.

Based on Theorem 6, the following theorem easily follows from Theorem 3.13 in
Mitzenmacher [20].

###### Theorem 7

In the supermarket model with general service times,
$\left\\{X_{n}\left(t\right)\right\\}$ and $\left\\{X\left(t\right)\right\\}$
are respectively given by (42) and (43), we have

$\lim_{n\rightarrow\infty}\sup_{u\leq
t}|X_{n}\left(u\right)-X\left(u\right)|=0,\text{ \ }a.s.$

Proof It is seen from that in the supermarket model with general service
times, the function $F\left(y\right)$ for $y\in\Omega$ satisfies the Lipschitz
condition. At the same time, it is easy to take a subset
$\Omega^{\ast}\subset\Omega$ such that

$\left\\{X\left(u\right):u\leq t\right\\}\subset\Omega^{\ast}$

and

$a\sup_{y\in\Omega^{\ast}}\beta_{a}\left(y\right)+a\sup_{y\in\Omega^{\ast}}\beta_{a}\left(y\right)<+\infty.$

Thus, this proof can easily be completed by means of Theorem 3.13 in
Mitzenmacher [20]. This completes the proof.

Using Theorem 3.11 in Mitzenmacher [20] and Theorem 7, we can obtain the
following theorem for the expected sojourn time that a customer spends in an
initially empty supermarket model with general service times over the time
interval $\left[0,T\right]$.

###### Theorem 8

In the supermarket model with general service times, the expected sojourn time
that a customer spends in an initially empty system over the time interval
$\left[0,T\right]$ is bounded above by

$\theta\rho^{d}\left\\{E\left[X_{R}\right]-E\left[X\right]\right\\}+E\left[X\right]\left[\sum_{k=1}^{\infty}\theta^{\frac{d^{k}-1}{d-1}}\rho^{\frac{d^{k}-d}{d-1}}\right]+o\left(1\right),$

where $o\left(1\right)$ is understood as $n\rightarrow\infty$.

##  8 Concluding remarks

In this paper, we provide a novel and simple approach to study the randomized
load balancing model with general service times, which is described as an
infinite-size system of integral-differential equations. This approach is
based on the supplementary variable method, which is always applied in dealing
with stochastic models of M/G/1 type, e.g., see Li and Zhao [13, 14] and Li
[10]. We organize an infinite-size system of integral-differential equations
by means of the density dependent jump Markov process, and obtain a close-form
solution: doubly exponential structure, for the fixed point satisfying the
system of nonlinear equations, which is always a key in the study of
supermarket models. Since the fixed point is decomposited into two groups of
information under a product form, we indicate three important observations:

  1. 1.

the fixed point for the supermarket model is different from the tail of
stationary queue length distribution for the ordinary M/G/1 queue;

  2. 2.

the doubly exponential solution to the fixed point can exist extensively for
$0<\mu<+\infty$ even if the service time distribution is heavy-tailed; and

  3. 3.

the doubly exponential solution to the fixed point is not unique for a more
general supermarket model.

Furthermore, we analyze the exponential convergence of the current location of
the supermarket model to its fixed point, and study the Lipschitz condition in
the Kurtz Theorem under general service times. Finally, we present numerical
examples to illustrate the effectiveness of our approach in analyzing the
randomized load balancing schemes with the non-exponential service
requirements. Based on this analysis, one can gain a new and important
understanding how workload probing can help in load balancing jobs with
general service times such as heavy-tailed service.

The approach of this paper is useful in analyzing the randomized load
balancing schemes in resource allocation in computer networks. We expect that
this approach will be applicable to the study other randomized load balancing
schemes with general service times, for example, generalizing the arrival
process to non-Poisson: the renewal arrival process or the Markovian arrival
process.

## Acknowledgements

The work of Q.L. Li was supported by the National Science Foundation of China
under grant No. 10871114 and the National Grand Fundamental Research 973
Program of China under grant No. 2006CB805901.

## References

  * [1] R. Adler, R. Feldman and M.S. Taqqu (1998). A Practical Guide to Heavy Tails: Statistical Techniques for Analyzing Heavy Tailed Distributions. Birkhäuser: Boston.
  * [2] Y. Azar, A.Z. Broder, A.R. Karlin and E. Upfal (1999). Balanced allocations. SIAM Journal on Computing 29, 180–200. A preliminary version of this paper appeared in Proceedings of the Twenty-Sixth Annual ACM Symposium on the Theory of Computing, 1994.
  * [3] M. Dahlin (1999). Interpreting stale load information. In Proceedings of the Nineteenth Annual IEEE International Conference on Distributed Computing Systems.
  * [4] K. Deimling (1977). Ordinary Differential Equations in Banach Spaces. Springer-Verlag, Lecture Notes in Math., Vol. 96.
  * [5] D.L. Eager, E.D. Lazokwska and J. Zahorjan (1986). Adaptive load sharing in homogeneous distributed systems. IEEE Transactions on Software Engineering 12, 662–675.
  * [6] D.L. Eager, E. D. Lazokwska and J. Zahorjan (1986). A comparison of receiver-initiated and sender-initiated adaptive load sharing. Performance Evaluation Review 6, 53–68.
  * [7] D.L. Eager, E.D. Lazokwska and J. Zahorjan (1988). The limited performance benefits of migrating active processes for load sharing. Performance Evaluation Review 16, 63–72.
  * [8] M. harchol-Balter and A.B. Downey. Exploiting process lifetime distributions for dynamic load balancing. ACM Transactions on Computer Systems 15, 253–285, 1997.
  * [9] T.G. Kurtz (1981). Approximation of Population Processes. SIAM.
  * [10] Q.L. Li (2010). Constructive Computation in Stochastic Models with Applications: The RG-Factorizations. Tsinghua Press and Springer.
  * [11] Q.L. Li, John C.S. Lui and Y. Wang (2010). A matrix-analytic solution for randomized load balancing models with PH service times. Submited for publication.
  * [12] Q.L. Li, John C.S. Lui and Y. Wang (2010). Doubly exponential asymptotics for randomized load balancing with Markovian arrival processes. Submited for publication.
  * [13] Q.L. Li and Y.Q. Zhao (2004). A MAP/G/1 queue with negative customers. Queueing Systems 47, 5-43.
  * [14] Q.L. Li, Y. Ying and Y.Q. Zhao (2006). A BMAP/G/1 retrial queue with a server subject to breakdowns and repairs. Annals of Operations Research 141, 233-270.
  * [15] M. Luczak and C. McDiarmid (2005). On the power of two choices: Balls and bins in continuous time. The Annals of Applied Probability 15, 1733–1764.
  * [16] M. Luczak and C. McDiarmid (2006). On the maximum queue length in the supermarket model. The Annals of Probability 34, 493–527.
  * [17] J.B. Martin (2001). Point processes in fast Jackson networks. Annals of Applied Probability 11, 650-663.
  * [18] J.B. Martin and Y.M Suhov (1999). Fast Jackson networks. Annals of Applied Probability 9, 854–870.
  * [19] M.D. Mitzenmacher (1996). Load balancing and density dependent jump Markov processes. In Proceedings of the Thirty-Seventh Annual Symposium on Foundations of Computer Science, pages 213–222.
  * [20] M.D. Mitzenmacher (1996). The power of two choices in randomized load balancing. PhD thesis, University of California at Berkeley, Department of Computer Science, Berkeley, CA, 1996.
  * [21] M. Mitzenmacher (1998). Analyses of load stealing models using differential equations. In Proceedings of the Tenth ACM Symposium on Parallel Algorithms and Architectures, pages 212–221.
  * [22] M. Mitzenmacher (1999). On the analysis of randomized load balancing schemes. Theory of Computing Systems 32, 361–386.
  * [23] M. Mitzenmacher (1999). Studying balanced allocations with differential equations. Combinatorics, Probability, and Computing 8, 473–482.
  * [24] M. Mitzenmacher (2000). How useful is old information? IEEE Transactions on Parallel and Distributed Systems 11, 6–20.
  * [25] M. Mitzenmacher (2001). The power of two choices in randomized load balancing. IEEE Transactions on Parallel and Distributed Computing 12, 1094-1104.
  * [26] M. Mitzenmacher, A. Richa, and R. Sitaraman (2001). The power of two random choices: a survey of techniques and results. In Handbook of Randomized Computing: volume 1, edited by P. Pardalos, S. Rajasekaran and J. Rolim, pp. 255-312.
  * [27] M. Mitzenmacher and B. Vöcking (1998). The asymptotics of selecting the shortest of two, improved. In Proceedings of the 37th Annual Allerton Conference on Communication, Control, and Computing, pages 326–327. A full version is available as Harvard Computer Science TR-08-99.
  * [28] R. Mirchandaney, D. Towsley, and J.A. Stankovic (1989). Analysis of the e ects of delays on load sharing. IEEE Transactions on Computers 38, 1513–1525.
  * [29] Y.M. Suhov and N.D. Vvedenskaya (2002). Fast Jackson Networks with Dynamic Routing. Problems of Information Transmission 38, 136{153.
  * [30] B. Vöcking (1999). How asymmetry helps load balancing. In Proceedings of the Fortieth Annual Symposium on Foundations of Computer Science, pages 131–140.
  * [31] N.D. Vvedenskaya, R.L. Dobrushin and F.I. Karpelevich (1996). Queueing system with selection of the shortest of two queues: An asymptotic approach. Problems of Information Transmissions 32, 20–34.
  * [32] N.D. Vvedenskaya and Y.M. Suhov (1997). Dobrushin’s mean-field approximation for a queue with dynamic routing. Markov Processes and Related Fields 3, 493–526.
  * [33] S. Zhou (1988). A trace-driven simulation study of dynamic load balancing. IEEE Transactions on Software Engineering, 1327–1341.