<figure><img src="image_1.png" /></figure>

HAL
open science

# A Matrix-Analytic Solution for Randomized Load Balancing Models with PH Service Times

Quan-Lin Li, John Lui, Yang Wang

► To cite this version:

Quan-Lin Li, John Lui, Yang Wang. A Matrix-Analytic Solution for Randomized Load Balancing Models with PH Service Times. Performance Evaluation of Computer and Communication Systems (PERFORM), Oct 2010, Vienna, Austria. pp.240-253, 10.1007/978-3-642-25575-5_20. hal-01586888

HAL Id: hal-01586888

https://inria.hal.science/hal-01586888

Submitted on 13 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Matrix-Analytic Solution for Randomized Load Balancing Models with PH Service Times

Quan-Lin Li John C.S. Lui Yang Wang *

**Abstract.** In this paper, we provide a matrix-analytic solution for randomized load balancing models (also known as *supermarket models*) with phase-type (PH) service times. Generalizing the service times to the phase-type distribution makes analysis of the supermarket models more difficult and challenging than that of the exponential service time case which has been extensively discussed in the literature. We describe the supermarket model as a system of differential vector equations, provide a doubly exponential solution to the fixed point of the system of differential vector equations, and analyze the exponential convergence of the current location of the supermarket model to its fixed point.

## 1 Introduction

In the past few years, a number of companies (e.g., Amazon, Google and Microsoft) are offering the *cloud computing* service to enterprises, and many content publishers and application service providers are increasingly using *Data Centers* to host their services. This emerging computing paradigm allows service providers and enterprises to concentrate on developing and providing the own services/goods without worrying about computing system maintenance or upgrade, and thereby significantly reduce their operating cost. For companies that offer cloud computing service in the data centers, they can take advantage of the variation of computing workloads from these customers and achieve the computational multiplexing gain. One of the important technical challenges that they have to address is how to utilize these computing resources in the data centers efficiently since many of these servers can be used. There is a growing interest to examine simple and robust load balancing strategies to efficiently utilize the computing resource of the server farms.

Randomized load balancing is a simple and effective mechanism to fairly utilize computing resources. It can deliver surprisingly good performance measures such as reducing collisions, waiting times and backlogs. In a supermarket model, each arriving job randomly picks a small subset of servers and examines their instantaneous workload, and the job is routed to the least loaded server. When a job is committed to a server, jockeying is not allowed and each server uses the

* Quan-Lin Li is from the School of Economics and Management Sciences, Yanshan University. John C.S. Lui is from the Department of Computer Science & Engineering, The Chinese University of Hong Kong. Yang Wang is from the Department of Computer Science and Technology, Peking University.
first-come-first-service (FCFS) discipline to process all jobs, e.g., see Mitzenmacher [7]. For the supermarket models, most of recent research applied density dependent jump Markov processes to deal with the simple case with Poisson arrival processes and exponential service times, and illustrated that there exists a fixed point which decreases doubly exponentially. Readers may refer to, such as, a simple supermarket model by [7, 14]; simple variations by [8, 9, 11, 4]; load information by [2, 10]; fast Jackson network by Martin and Suhov [6, 5, 12]; and general service times by Bramson, Lu and Prabhakar [1]. When the arrival processes or the service times are more general, the available results of the supermarket models are few up to now. The purpose of this paper is to provide a novel approach for studying a supermarket model with PH service times, and show that the fixed point decreases doubly exponentially. Also, the PH approximation of order 2 gains new numerical insights on practical applications of the general service times to supermarket models.

The remainder of this paper is organized as follows. In the next section, we describe the supermarket model with the PH service times as a system of differential vector equations based on density dependent jump Markov processes. In Section 3, we set up a system of nonlinear equations satisfied by the fixed point, provide a doubly exponential solution to the system of nonlinear equations, and compute the expected sojourn time of a tagged arriving customer. In Section 4, we study exponential convergence of the current location of the supermarket model to its fixed point.

## 2 Supermarket Model

In this section, we describe a supermarket model with a Poisson arrival process and PH service times as a system of differential vector equations based on density dependent jump Markov processes.

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 arrival rate $n\lambda$ for $\lambda > 0$. The service time of each customer is of phase type with irreducible representation $(\alpha, T)$ of order $m$. Each arriving customer chooses $d \ge 1$ servers independently and uniformly at random from the $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 every server will be served in the FCFS manner. Please see Figure 1 for an illustration of the supermarket model.

For the supermarket models, the PH distribution allows us to model more realistic systems and understand their performance implication under the randomized load balancing strategy. As indicated in [3], the process lifetime of many parallel jobs, in particular, jobs to the data centers, tends to be non-exponential. For the PH service time distribution, we use the following irreducible representation: $(\alpha, T)$ of order $m$, the row vector $\alpha$ is a probability vector whose $j$th entry is the probability that a service begins in phase $j$ for $1 \le j \le m$; $T$ is an
<figure><img src="image_3.png" /><figcaption>Fig. 1. The supermarket model: each customer can probe the loading of d servers</figcaption></figure>

$m \times m$ matrix whose $(i, j)^{th}$ entry is denoted by $t_{i,j}$ with $t_{i,i} < 0$ for $1 \le i \le m$, and $t_{i,j} \ge 0$ for $1 \le i, j \le m$ and $i \ne j$. Let $T^0 = -Te \ge 0$, where $e$ is a column vector of ones with a suitable dimension in the context. The expected service time is given by $1/\mu = -\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 in the stable region $\rho = \lambda/\mu < 1$.

We define $n_k^{(i)}(t)$ as the number of queues with at least $k$ customers and the PH service time in phase $i$ at time $t \ge 0$. Clearly, $0 \le n_k^{(i)}(t) \le n$ for $k \ge 1$ and $1 \le i \le m$. Let

$$X_n^{(0)}(t) = \frac{n}{n} = 1,$$

and $k \ge 1$

$$X_n^{(k,i)}(t) = \frac{n_k^{(i)}(t)}{n},$$

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

$$X_n^{(k)}(t) = (X_n^{(k,1)}(t), X_n^{(k,2)}(t), \dots, X_n^{(k,m)}(t)), \quad k \ge 1,$$

$$X_n(t) = (X_n^{(0)}(t), X_n^{(1)}(t), X_n^{(2)}(t), \dots).$$

The state of the supermarket model may be described by the vector $X_n(t)$ for $t \ge 0$. Since the arrival process to the queueing system is Poisson and the service time of each server is of phase type, the stochastic process $\{X_n(t), t \ge 0\}$ is a Markov process whose state space is given by

$$\Omega_n = \{ (g_n^{(0)}, g_n^{(1)}, \dots) : g_n^{(0)} = 1, g_n^{(k)} \ge g_n^{(k+1)} \ge 0, \\ \text{and } ng_n^{(k)} \text{ is a vector of nonnegative integers for } k \ge 1 \}.$$

Let

$$s_0(n,t) = E[X_n^{(0)}(t)]$$

and $k \ge 1$

$$s_k^{(i)}(n,t) = E[X_n^{(k,i)}(t)].$$
Clearly, $s_0(n, t) = 1$. We write

$$S_k(n,t) = (s_k^{(1)}(n,t), s_k^{(2)}(n,t), \dots, s_k^{(m)}(n,t)), \quad k \ge 1.$$

As shown in Martin and Suhov [6] and Luczak and McDiarmid [4], the Markov process $\{X_n(t), t \ge 0\}$ is asymptotically deterministic as $n \to \infty$, this is due to the fact that from a block-structured point of view, the PH distribution can keep many excellent properties of the exponential distribution. Thus $\lim_{n \to \infty} E[X_n^{(0)}(t)]$ and $\lim_{n \to \infty} E[X_n^{(k,i)}]$ always exist by means of the law of large numbers. Based on this, we write

$$S_0(t) = \lim_{n \to \infty} s_0(n,t) = 1,$$

for $k \ge 1$

$$s_k^{(i)}(t) = \lim_{n \to \infty} s_k^{(i)}(n,t), \\ S_k(t) = (s_k^{(1)}(t), s_k^{(2)}(t), \dots, s_k^{(m)}(t))$$

and

$$S(t) = (S_0(t), S_1(t), S_2(t), \dots).$$

Let $X(t) = \lim_{n \to \infty} X_n(t)$. Then it is easy to see from the Poisson arrivals and the PH service times that $\{X(t), t \ge 0\}$ is also a Markov process whose state space is given by

$$\Omega = \left\{ (g^{(0)}, g^{(1)}, g^{(2)}, \dots) : g^{(0)} = 1, g^{(k)} \ge g^{(k+1)} \ge 0 \text{ for } k \ge 1 \right\}.$$

If the initial distribution of the Markov process $\{X_n(t), t \ge 0\}$ approaches the Dirac delta-measure concentrated at a point $g \in \Omega$, then its steady-state distribution is concentrated in the limit on the trajectory $S_g = \{S(t) : t \ge 0\}$. This indicates a law of large numbers for the time evolution of the fraction of queues of different lengths. Furthermore, the Markov process $\{X_n(t), t \ge 0\}$ converges weakly to the fraction vector $S(t) = (S_0(t), S_1(t), S_2(t), \dots)$, or for a sufficiently small $\varepsilon > 0$,

$$\lim_{n \to \infty} P \{ \|X_n(t) - S(t)\| \ge \varepsilon \} = 0,$$

where $\|a\|$ is the $L_\infty$-norm of vector $a$.

In what follows we provide a system of differential vector equations in order to determine fraction vector $S(t)$. To that end, we introduce the Hadamard Product of two matrices $A = (a_{i,j})$ and $B = (b_{i,j})$ as follows:

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

Specifically, for $k \ge 2$ we have

$$A^{\odot k} = \underbrace{A \odot A \odot \dots \odot A}_{k \text{ matrix } A}.$$
To determine the fraction vector $S(t)$, we need to set up a system of differential vector equations satisfied by $S(t)$. To that end, we provide a heuristic description in terms of a concrete example, and indicate how to derive the differential vector equations.

In the supermarket model with $n$ servers, we determine the expected change in the number of queues with at least $k$ customers over a small time period of length $dt$. The probability vector that during the time period $[0,dt)$, any arriving customer joins a queue with $k-1$ customers is given by

$$n [\lambda S_{k-1}^{\odot d}(n, t) - \lambda S_k^{\odot d}(n, t)] dt.$$

Similarly, the probability vector that during the time period $[0,dt)$, a customer leaves a server queued by $k$ customers is given by

$$n [S_k(n, t) T + S_{k+1}(n, t) T^0 \alpha] dt.$$

Therefore we can obtain

$$dE[n_k(n, t)] = n [\lambda S_{k-1}^{\odot d}(n, t) - \lambda S_k^{\odot d}(n, t)] dt + n [S_k(n, t) T + S_{k+1}(n, t) T^0 \alpha] dt,$$

which leads to

$$\frac{dS_k(n, t)}{dt} = \lambda S_{k-1}^{\odot d}(n, t) - \lambda S_k^{\odot d}(n, t) + S_k(n, t) T + S_{k+1}(n, t) T^0 \alpha. \quad (1)$$

Taking $n \to \infty$ in both sides of Equation (1), we have

$$\frac{dS_k(t)}{dt} = \lambda S_{k-1}^{\odot d}(t) - \lambda S_k^{\odot d}(t) + S_k(t)T + S_{k+1}(t)T^0\alpha. \quad (2)$$

Using a similar analysis to Equation (2), we obtain a system of differential vector equations for the fraction vector $S(t) = (S_0(t), S_1(t), S_2(t), \dots)$ as follows:

$$S_0(t) = 1,$$

$$\frac{d}{dt} S_0(t) = -\lambda S_0^d(t) + S_1(t) T^0, \quad (3)$$

$$\frac{d}{dt} S_1(t) = \lambda \alpha S_0^d(t) - \lambda S_1^{\odot d}(t) + S_1(t) T + S_2(t) T^0 \alpha, \quad (4)$$

and for $k \ge 2$,

$$\frac{d}{dt} S_k(t) = \lambda S_{k-1}^{\odot d}(t) - \lambda S_k^{\odot d}(t) + S_k(t) T + S_{k+1}(t) T^0 \alpha. \quad (5)$$

Note that the above derivations or the ordinary differential equations (ODE) descriptions are the approximation to the large scale stochastic system. One can be made the formal derivation using the “mean field” or the “stochastic density jump Markov process” techniques. While it is in form of ODE, we believe readers can easily understand the ODE of the above equations since the supplementary variable method are extensively used in the study of stochastic systems.
## 3 A Matrix-Analytic Solution

In this section, we provide a doubly exponential solution to the fixed point of the system of differential vector equations (3), (4) and (5).

A row vector $\pi = (\pi_0, \pi_1, \pi_2, ...)$ is called a fixed point of the fraction vector $S(t)$ if $\lim_{t \to +\infty} S(t) = \pi$. In this case, it is easy to see that

$$ \lim_{t \to +\infty} \left[ \frac{d}{dt} S(t) \right] = 0. $$

Therefore, as $t \to +\infty$ the system of differential vector equations (3), (4) and (5) can be simplified as

$$ -\lambda \pi_0^d + \pi_1 T^0 = 0, \qquad (6) $$

$$ \lambda \alpha \pi_0^d - \lambda \pi_1^{\odot d} + \pi_1 T + \pi_2 T^0 \alpha = 0, \qquad (7) $$

and for $k \ge 2$,

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

In general, it is more difficult and challenging to express the fixed point of the supermarket models with more general arrival processes or service times. Fortunately, we can derive a closed-form expression for the fixed point $\pi = (\pi_0, \pi_1, \pi_2, ...)$ for the supermarket model with PH service times by means of a novel matrix-analytic approach.

Noting that $S_0(t) = 1$ for all $t \ge 0$, it is easy to see that $\pi_0 = 1$. It follows from Equation (6) that

$$ \pi_1 T^0 = \lambda. \qquad (9) $$

To solve Equation (9), we denote by $\omega$ the stationary probability vector of the irreducible Markov chain $T + T^0\alpha$. Obviously, we have

$$ \begin{aligned} \omega T^0 &= \mu, \\ \frac{\lambda}{\mu} \omega T^0 &= \lambda. \end{aligned} \qquad (10) $$

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 (7) that

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

which leads to

$$ \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 (0, 1)$, and

$$ \pi_2 T^0 = \lambda \theta \rho^d. $$

Using a similar analysis to Equation (10), we have

$$ \pi_2 = \frac{\lambda \theta \rho^d}{\mu} \omega = \theta \rho^{d+1} \cdot \omega. \qquad (11) $$
Based on $\pi_1 = \rho \cdot \omega$ and $\pi_2 = \theta \rho^{d+1} \cdot \omega$, it follows from Equation (8) 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

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

Using a similar analysis on Equation (10), 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. \quad (12) $$

Based on Equations (11) and (12), we may infer that there is a structured expression $\pi_k = \theta^{d^{k-2}+d^{k-3}+\dots+d+1} \rho^{d^{k-1}+d^{k-2}+\dots+d+1} \cdot \omega$ for $k \ge 1$. To that end, the following theorem states this important result.

**Theorem 1** The fixed point $\pi = (\pi_0, \pi_1, \pi_2, \dots)$ is given by

$$ \pi_0 = 1, \quad \pi_1 = \rho \cdot \omega $$

and for $k \ge 2$,

$$ \pi_k = \theta^{d^{k-2}+d^{k-3}+\dots+d+1} \rho^{d^{k-1}+d^{k-2}+\dots+d+1} \cdot \omega, \quad (13) $$

or

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

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

It is clear that Equation (13) is correct for the cases with $l = 2, 3$ according to Equations (11) and (12). Now, we assume that Equation (13) is correct for the cases with $l = k$. Then it follows from Equation (8) that for $l = k + 1$, we have

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

which leads to

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

By a similar analysis to (10), we have

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

This completes the proof.

Now, we compute the expected sojourn time $T_d$ that a tagged arriving customer spends in the supermarket model. For the PH service times, a tagged
arriving customer is the kth customer in the corresponding queue with probability vector $\pi_{k-1}^{\odot d} - \pi_k^{\odot d}$. When $k \ge 1$, the head customer in the queue has been served, and so its service time is residual and is denoted as $X_R$. Let $X$ be of phase type with irreducible representation $(\alpha, T)$. Then $X_R$ is of phase type with irreducible representation $(\omega, T)$. Clearly, we have

$$E[X] = \alpha (-T)^{-1} e, \quad E[X_R] = \omega (-T)^{-1} e.$$

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

$$ 
\begin{aligned}
E[T_d] &= (\pi_0^d - \pi_1^{\odot d} e) E[X] + \sum_{k=1}^{\infty} (\pi_k^{\odot d} - \pi_{k+1}^{\odot d}) e \{E[X_R] + kE[X]\} \\
&= \rho^d \theta (\omega - \alpha) (-T)^{-1} e + \alpha (-T)^{-1} e \left(1 + \sum_{k=1}^{\infty} \theta^{\frac{d k - 1}{d-1}} \rho^{\frac{d k + 1 - d}{d-1}}\right). \quad (14)
\end{aligned}
$$

When the arrival process and the service time distribution are Poisson and exponential, respectively, it is clear that $\alpha = \omega = \theta = 1$ and $\alpha (-T)^{-1} e = 1/\mu$, thus we have

$$E[T_d] = \frac{1}{\mu} \sum_{k=0}^{\infty} \rho^{\frac{d k + 1 - d}{d-1}},$$

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

In what follows we consider an interesting problem: When using the PH approximation, how many moments of the service time distribution are needed to obtain a better accuracy for computing the fixed point or the expected sojourn time. It is well-known from the theory of probability distributions that the first three moments are basic for analyzing such an accuracy. Also, we can construct a PH distribution of order 2 by using the first three moments. Telek and Heindl [13] provided a fitting procedure for matching a PH distribution of order 2 in terms of the first three moments exactly. It is necessary to list the fitting procedure as follows:

For a nonnegative random variable $X$, let $m_n = E[X^n]$, $n \ge 1$. We take a PH distribution of order 2 with the canonical representation $(\alpha, T)$, where $\alpha = (\eta, 1 - \eta)$ and

$$T = \begin{pmatrix} -\xi_1 & \xi_1 \\ 0 & -\xi_2 \end{pmatrix},$$

$0 \le \eta \le 1$ and $0 < \xi_1 \le \xi_2$. Note that the three unknown parameters $\eta$, $\xi_1$ and $\xi_2$ can be obtained from the first three moments $m_1, m_2$ and $m_3$.

In Table 1, $c_X^2 = m_2/m_1^2 - 1$ which is the squared coefficient of variation. If the moments do not satisfy these conditions in Table 1, then we may analyze the following four cases:

(a.1) if $m_2 < 1.5m_1^2$, then we take $m_2 = 1.5m_1^2$;

(a.2) if $0.5 \le c_X^2 \le 1$ and $m_3 < 3m_1^3 (3c_X^2 - 1 + \sqrt{2}(1 - c_X^2)^{3/2})$, then we take $m_3 = 3m_1^3 (3c_X^2 - 1 + \sqrt{2}(1 - c_X^2)^{3/2})$;
<table><caption>**Table 1.** Specific Bounds of the First Three Moments</caption><thead><tr><th>Moment</th><th>Condition</th><th>Bounds</th></tr></thead><tbody><tr><td>m<sub>1</sub></td><td></td><td>0 &lt; m<sub>1</sub> &lt; &infin;</td></tr><tr><td>m<sub>2</sub></td><td></td><td>1.5m<sub>1</sub><sup>2</sup> &le; m<sub>2</sub></td></tr><tr><td>m<sub>3</sub></td><td>0.5 &le; c<sub>X</sub><sup>2</sup> &le; 1</td><td>3m<sub>1</sub><sup>3</sup>(3c<sub>X</sub><sup>2</sup> - 1 + &radic;2(1 - c<sub>X</sub><sup>2</sup>)<sup>2/3</sup>) &le; m<sub>3</sub> &le; 6m<sub>1</sub><sup>3</sup>c<sub>X</sub><sup>2</sup></td></tr><tr><td></td><td>1 &lt; c<sub>X</sub><sup>2</sup></td><td><sup>3</sup>&frasl;<sub>2</sub>m<sub>1</sub><sup>3</sup>(1 + c<sub>X</sub><sup>2</sup>)<sup>2</sup> &lt; m<sub>3</sub> &lt; &infin;</td></tr></tbody></table>

(a.3) if $0.5 \le c_X^2 \le 1$ and $m_3 > 6m_1^3c_X^2$, then we take $m_3 = 6m_1^3c_X^2$; and

(a.4) if $1 < c_X^2$ and $m_3 \le \frac{3}{2}m_1^3(1+c_X^2)^2$, then we take $m_3 = \frac{3}{2}m_1^3(1+c_X^2)^2$. Let $c = 3m_2^2 - 2m_1m_3$, $d = 2m_1^2 - m_2$, $b = 3m_1m_2 - m_3$ and $a = b^2 - 6cd$.

If the first three moments satisfy their specific bounds shown in Table 1 or the exceptive four cases, then the three unknown parameters $\eta$, $\xi_1$ and $\xi_2$ can be computed in the following three cases.

(1) If $c > 0$, then

$$ \eta = \frac{-b + 6m_1d + \sqrt{a}}{b + \sqrt{a}}, \quad \xi_1 = \frac{b - \sqrt{a}}{c}, \quad \xi_2 = \frac{b + \sqrt{a}}{c}. $$

(2) If $c < 0$, then

$$ \eta = \frac{b - 6m_1d + \sqrt{a}}{-b + \sqrt{a}}, \quad \xi_1 = \frac{b + \sqrt{a}}{c}, \quad \xi_2 = \frac{b - \sqrt{a}}{c}. $$

(3) If $c = 0$, then

$$ \eta = 0, \quad \xi_1 > 0, \quad \xi_2 = \frac{1}{m_1}. $$

From the above discussion, we can construct a PH distribution of order 2 to approximate an arbitrarily given general distribution under the same first three moments. In fact, such an approximation achieves a better accuracy in computations of the fixed point and the expected sojourn time.

For the PH distribution of order 2, we have

$$ T + T^0\alpha = \begin{pmatrix} -\xi_1 & \xi_1 \\ 0 & -\xi_2 \end{pmatrix} + \begin{pmatrix} 0 \\ \xi_2 \end{pmatrix} (\eta 1 - \eta) = \begin{pmatrix} -\xi_1 & \xi_1 \\ \xi_2\eta & -\xi_2\eta \end{pmatrix}, $$

which leads to

$$ \omega = \left( \frac{\xi_2\eta}{\xi_1 + \xi_2\eta}, \frac{\xi_1}{\xi_1 + \xi_2\eta} \right) $$

and

$$ \theta = \frac{\xi_1^d + \xi_2^d \eta^d}{(\xi_1 + \xi_2 \eta)^d}. $$

Note that the PH distributions are dense in the set of all nonnegative random variables, we can numerically provide necessary understanding for the role played
by the general service times in performance analysis of the supermarket model by means of the PH approximation of order 2. At the same time, all computations involved in the PH approximation of order 2 are very simple to implement.

In the remainder of this section, we provide two examples to illustrate that our approach is effective and efficient in analyzing the doubly exponential solution of supermarket models with non-exponential service requirements.

**Example one** (PH Distribution) Let $\lambda = 1, d = 5, m = 3, \alpha(1) = (1/3, 1/3, 1/3)$ and $\alpha(2) = (1/12, 7/12, 1/3)$,

$$T = \begin{pmatrix} -10 & 2 & 4 \\ 3 & -7 & 4 \\ 0 & 2 & -5 \end{pmatrix}.$$

Table 2 shows how the doubly exponential solution ($\pi_1$ to $\pi_4$) depends on the vectors $\alpha(1)$ and $\alpha(2)$, respectively.

<table><caption>**Table 2.** The doubly exponential solution depends on the vector $\alpha$</caption><thead><tr><th></th><th>&alpha; = (1/3, 1/3, 1/3)</th><th>&alpha; = (1/12, 7/12, 1/3)</th></tr></thead><tbody><tr><td>&pi;<sub>1</sub></td><td>(0.0741, 0.1358, 0.2346)</td><td>(0.0602, 0.1728, 0.2531)</td></tr><tr><td>&pi;<sub>2</sub></td><td>(5.619e-05, 1.030e-05, 1.779e-04)</td><td>(7.182e-05, 2.063e-04, 3.020e-04)</td></tr><tr><td>&pi;<sub>3</sub></td><td>(1.411e-20, 2.587e-20, 4.469e-20)</td><td>(1.739e-19, 4.993e-19, 7.311e-19)</td></tr><tr><td>&pi;<sub>4</sub></td><td>(1.410e-98, 2.586e-98, 4.466e-98)</td><td>(1.444e-92, 4.148e-92, 6.074e-92)</td></tr></tbody></table>

**Example two** (Expected Sojourn Time) We consider an *m*-order PH distribution with irreducible representation $(\alpha, T)$, where $m = 2, \alpha = (1/2, 1/2)$ and

$$T = \begin{pmatrix} -4 & 3 \\ 2 & -7 \end{pmatrix}.$$

Let $\lambda = 1, \mu = 2.7500$ and $d = 2$. Figure 2 provides a comparison for the expected sojourn times between the exponential service time and the PH service time when they have a same expected service time. This figure shows that the PH service time makes the lower expected sojourn time.

# 4 Exponential convergence to the fixed point

In this section, we study exponential convergence of the current location $S(t)$ of the supermarket model to its fixed point $\pi$.

We provide some notation for comparison of two vectors. Let $a = (a_1, a_2, a_3, ...)$ and $b = (b_1, b_2, b_3, ...)$.

We write $a \prec b$ if $a_k < b_k$ for some $k \ge 1$ and $a_l \le b_l$ for $l \ne k, l \ge 1$; and $a \preceq b$ if $a_k \le b_k$ for all $k \ge 1$.

**Proposition 1** If $S(0) \preceq \tilde{S}(0)$, then $S(t) \preceq \tilde{S}(t)$.
<figure><img src="image_4.png" /><figcaption>Fig. 2. The expected sojourn times correspond to the exponential and PH distributions</figcaption></figure>

It follows from the system of differential vector equations (3), (4) and (5) that for $d \ge 2$

$$ \frac{d}{dt} S(t) = S^{\odot d}(t) Q_1 + S(t) Q_2, \qquad (15) $$

where

$$ Q_1 = \begin{pmatrix} -\lambda I & \lambda I \\ -\lambda I & \lambda I \\ -\lambda I & \lambda I \\ \vdots & \vdots \\ \vdots & \vdots \end{pmatrix}, \quad Q_2 = \begin{pmatrix} 0 & & & \\ T^0 \alpha & T & & \\ & T^0 \alpha & T & \\ & & \ddots & \ddots \\ & & & T^0 \alpha & T \end{pmatrix}. $$

It follows from (15) that

$$ \frac{d}{dt} S^{\odot (1-d)}(t) = eQ_1 + S^{\odot (1-d)}(t) Q_2. $$

Let $W(t) = S^{\odot (1-d)}(t)$. Then

$$ \frac{d}{dt} W(t) = eQ_1 + W(t) Q_2. $$

Hence we obtain

$$ W(t) = \left[ W(0) + e \int_{0}^{t} \exp \{Q_2 x\} dx \right] \exp \{Q_1 t\}. $$

Since $W(t) = S^{\odot (1-d)}(t)$, we have

$$ W(0) = S^{\odot (1-d)}(0) $$

and

$$ S(t) = \left[ S^{\odot (1-d)}(t) \right]^{\odot \frac{1}{1-d}} = [W(t)]^{\odot \frac{1}{1-d}}. $$
Therefore, the solution to the system of differential vector equations (3), (4) and (5) is given by

$$S(t) = \left\{ \left[ S^{\odot(1-d)}(0) + e \int_{0}^{t} \exp\{Q_2x\} dx \right] \exp\{Q_1t\} \right\}^{\odot \frac{1}{1-d}} . \quad (16)$$

Note that $\int_{0}^{t} \exp\{Q_2x\}dx > 0$ and $\exp\{Q_1t\} > 0$, it is easy to see from (16) that if $S(0) \preceq \tilde{S}(0)$, then $S(t) \preceq \tilde{S}(t)$ for $t > 0$. This completes the proof. ■

Based on Proposition 1, the following theorem shows that the fixed point $\pi$ is an upper bound of the current location $S(t)$ for all $t \ge 0$.

**Theorem 2** For the supermarket model, if there exists some $k$ such that $S_k(0) = 0$, then the sequence $\{S_k(t)\}$ has an upper bound sequence which decreases doubly exponentially for all $t \ge 0$, that is, $S(t) \preceq \pi$ for all $t \ge 0$.

**Proof:** Let $\tilde{S}_k(0) = \pi_k$ for $k \ge 1$. Then for each $k \ge 1$, $\tilde{S}_k(t) = \tilde{S}_k(0) = \pi_k$ for all $t \ge 0$, since $\tilde{S}(0) = (\tilde{S}_1(0), \tilde{S}_2(0), \tilde{S}_3(0), \dots)$ is a fixed point in the supermarket model. If $S_k(0) = 0$ for some $k$, then $S_k(0) < \tilde{S}_k(0)$ and $S_j(0) \preceq \tilde{S}_j(0)$ for $j \ne k, j \ge 1$, thus $S(0) < \tilde{S}(0)$. It is easy to see from Proposition 1 that $S_k(t) \preceq \tilde{S}_k(t) = \pi_k$ for all $k \ge 1$ and $t \ge 0$. Thus we obtain that for all $k \ge 1$ and $t \ge 0$

$$S_k(t) \le \theta^{\frac{d^{k-1}-1}{d-1}} \rho^{\frac{d^{k-1}}{d-1}} \cdot \omega.$$

This completes the proof. ■

To show the exponential convergence, we define a Lyapunov function $\Phi(t)$ as

$$\Phi(t) = \sum_{k=1}^{\infty} w_k [\pi_k - S_k(t)] e$$

in terms of the fact that $S_k(t) \preceq \pi_k$ for $k \ge 1$ and $\pi_0 = S_0(t) = 1$, where $\{w_k\}$ is a positive scalar sequence with $w_{k+1} \ge w_k \ge w_1 = 1$ for $k \ge 2$.

The following theorem illustrates that the distance between the fixed point and the current location quickly comes close to zero with exponential convergence.

**Theorem 3** For $t \ge 0$, $\Phi(t) \le c_0 e^{-\delta t}$, where $c_0$ and $\delta$ are two positive constants. In this case, the potential function $\Phi(t)$ is exponentially convergent.

**Proof:** Note that

$$\Phi(t) = \sum_{k=1}^{\infty} w_k [\pi_k - S_k(t)] e,$$

we have

$$\frac{d}{dt}\Phi(t) = -\sum_{k=1}^{\infty} w_k \frac{d}{dt}S_k(t)e.$$
It follows from Equations (3) to (5) that

$$ \frac{d}{dt}\Phi(t) = -w_1[\lambda S_0^d(t)\alpha - \lambda S_1^{\odot d}(t) + S_1(t)T + S_2(t)T^0\alpha]e \\ - \sum_{k=1}^{\infty} w_k[\lambda S_{k-1}^{\odot d}(t) - \lambda S_k^{\odot d}(t) + S_k(t)T + S_{k+1}(t)T^0\alpha]e. $$

By means of $S_0(t) = 1$ and $Te = -T^0$, we can obtain

$$ \frac{d}{dt}\Phi(t) = -w_1[\lambda - \lambda S_1^{\odot d}(t) e - S_1(t)T^0 + S_2(t)T^0] \\ - \sum_{k=2}^{\infty} w_k[\lambda S_{k-1}^{\odot d}(t) e - \lambda S_k^{\odot d}(t) e - S_k(t)T^0 + S_{k+1}(t)T^0]. \quad (17) $$

We take some nonnegative constants $c_k(t)$ and $d_k(t)$ for $k \ge 1$ such that

$$ \lambda = f_1(t) S_1(t) T^0, $$

for $k \ge 1$

$$ \lambda S_k^{\odot d}(t) e = c_k(t) [\pi_k - S_k(t)] e $$

and

$$ S_k(t) T^0 = d_k(t) [\pi_k - S_k(t)] e. $$

Then it follows from (17) that

$$ \frac{d}{dt}\Phi(t) = -\{[(w_2 - w_1)] c_1(t) + w_1[f_1(t) - 1] d_1(t)\} \cdot [\pi_1 - S_1(t)] e \\ - \sum_{k=2}^{\infty} [(w_{k+1} - w_k) c_k(t) + (w_{k-1} - w_k) d_k(t)] \cdot [\pi_k - S_k(t)] e. $$

For a constant $\delta > 0$, we take

$$ w_1 = 1, $$

$$ [(w_2 - w_1)] c_1(t) + w_1[f_1(t) - 1] d_1(t) \ge \delta w_1 $$

and

$$ (w_{k+1} - w_k) c_k(t) + (w_{k-1} - w_k) d_k(t) \ge \delta w_k. $$

In this case, it is easy to see that

$$ w_2 \ge 1 + \frac{\delta + 1 - f_1(t)}{c_1(t)} $$

and for $k \ge 2$

$$ w_{k+1} \ge w_k + \frac{\delta w_k}{c_k(t)} + \frac{d_k(t)}{c_k(t)} (w_k - w_{k-1}). $$

Thus we have

$$ \frac{d}{dt}\Phi(t) \le -\delta \sum_{k=0}^{\infty} w_k [\pi_k - S_k(t)] e = -\delta \Phi(t), $$
which can leads to

$$ \Phi(t) \leq c_0 e^{-\delta t}. $$

This completes the proof.

<figure><img src="image_5.png" /></figure>

## Acknowledgements

The work of Quan-Lin Li was supported by the National Science Foundation of China under grant No. 10871114, John C.S. Lui was supported by the RGC grant, and Yang Wang was supported by the National Science Foundation of China under grant No. 61001075.

## References

1. Bramson, M., Lu, Y., Prabhakar, B.: Randomized load balancing with general service time distributions. In: Proceedings of the ACM SIGMETRICS international conference on Measurement and modeling of computer systems, pp. 275–286 (2010)

2. Dahlin, M.: Interpreting stale load information. IEEE Transactions on Parallel and Distributed Systems 11, 1033–1047 (1999)

3. Harchol-Balter, M., Downey, A.B.: Exploiting process lifetime distributions for dynamic load balancing. ACM Transactions on Computer Systems 15, 253–285 (1997)

4. Luczak, M., McDiarmid, C.: On the maximum queue length in the supermarket model. The Annals of Probability 34, 493–527 (2006)

5. Martin, J.B.: Point processes in fast Jackson networks. The Annals of Applied Probability 11, 650–663 (2001)

6. Martin, J.B., Suhov, Y.M.: Fast Jackson networks. The Annals of Applied Probability 9, 854–870 (1999)

7. Mitzenmacher, M.D.: The power of two choices in randomized load balancing. PhD thesis, University of California at Berkeley, Department of Computer Science, Berkeley, CA (1996).

8. Mitzenmacher, M.D.: Analyses of load stealing models using differential equations. In: Proceedings of the Tenth ACM Symposium on Parallel Algorithms and Architectures, pp. 212–221 (1998)

9. Mitzenmacher, M.D.: On the analysis of randomized load balancing schemes. Theory of Computing Systems 32, 361–386 (1999)

10. Mitzenmacher, M.D.: How useful is old information? IEEE Transactions on Parallel and Distributed Systems 11, 6–20 (2000)

11. Mitzenmacher, M.D., Richa, A., Sitaraman, R.: The power of two random choices: a survey of techniques and results. In: Handbook of Randomized Computing: Volume 1, Pardalos, P., Rajasekaran, S., Rolim, J. (eds.), pp. 255–312 (2001)

12. Suhov, Y.M., Vvedenskaya, N.D.: Fast Jackson Networks with Dynamic Routing. Problems of Information Transmission 38, 136–153 (2002)

13. Telek, M., Heindl, A.: Matching moments for acyclic discrete and continuous phase-type distributions of second order. International Journal of Simulation: Systems, Science & Technology 3, 47–57 (2002)

14. Vvedenskaya, N.D., Dobrushin, R.L., Karpelevich, F.I.: Queueing system with selection of the shortest of two queues: An asymptotic approach. Problems of Information Transmissions 32, 20–34 (1996)