Large-System Insensitivity of Zero-Waiting Load Balancing Algorithms

1. Introduction

Zero-wait load balancing refers to a load balancing algorithm in which a job is routed to an idle server and is processed immediately on its arrival. The zero-waiting property has become increasingly important as the amount of modern machine learning (ML) and artificial intelligence (AI) applications running on large-scale data centers explodes. Although increasing the number of servers and the processing speed of each server is a critical step to meet the increasing demand, the design of load balancing algorithms that can efficiently use available resources to minimize or even eliminate the waiting time of incoming jobs is equally important, especially when a minor increase of latency (e.g., 100 milliseconds) can lead to a significant drop in a cloud-computing provider’s revenue (7% drop in sales according to a recent Akamai report; Akamai 2017).

Significant progress has been made over the past few years on achieving asymptotic zero-waiting (as the system size approaches infinity) in a large-scale data center with distributed queues, including the classic supermarket model (Stolyar 2015; Mukherjee et al. 2016; Banerjee and Mukherjee 2018; Eschenfeldt and Gamarnik 2018; Liu and Ying 2018, 2020, 2021; Gupta and Walton 2019; Braverman 2020; Liu et al. 2021), models with data locality (Weng et al. 2020, Rutten and Mukherjee 2022), and models where each job consists of parallel tasks (Weng and Wang 2020, Hong and Wang 2021, Wang et al. 2021).

However, almost all these results assume exponential service time distributions. Although each of these results provided important insights into achieving zero-waiting in a practical system, theoretically, it is not clear whether these principles hold for general service times. This is a very important question to answer because it is well known that service time distributions in real-world systems are not exponential. Understanding a queueing system’s performance with general service times remains one of the most important and intensively studied problems in stochastic networks (Schelling 1978, Whittle 1985, Zachary 2007). A concept that excited many theorists in the area is “insensitivity” (Burman 1981). A queueing system is called insensitive if the steady-state distribution of queue lengths is invariant to the service time distribution. Therefore, any conclusion drawn from exponential service time distributions can be applied to general service time distributions. A result that is insensitive is robust and is expected to be widely applicable in practical systems. Unfortunately, insensitivity results are rare and often hold only under some special queueing disciplines such as processor sharing (PS) (Burman 1981, Bonald et al. 2004, Jonckheere and Prabhu 2018, Kielanski and Van Houdt 2021). One of the reasons is that insensitivity, although appealing, is a very strong notion of “robustness.” It requires the steady-state distribution under a general service time distribution to be exactly the same as that under the exponential distribution. Some recent studies started to relax it to weaker notions such as insensitivity in the heavy-traffic regime (Wang et al. 2018) or the large-system regime (Bramson et al. 2012, Stolyar 2015), that is, insensitivity in the limiting regimes. In the light of these recent developments, this paper addresses the following important question:

Are the zero-waiting algorithms insensitive and if so, for which notion of insensitivity?

1.1. Main Contributions

This paper provides some positive answers to the previous question. First, it is well known that most zero-waiting algorithms, such as join-the-shortest-queue (JSQ) (Winston 1977) and join-the-idle-queue (JIQ) (Lu et al. 2011), are not insensitive (according to its original definition). However, we prove that in the sub-Halfin-Whitt regime, LB-zero identified in Liu and Ying (2020) in fact achieves asymptotic zero-waiting for jobs whose service time follows a Coxian distribution with a finite number of phases. (Informally speaking, LB-zero refers to a class of load balancing algorithms that assigns an incoming job to an idle server with a high probability and achieves asymptotic zero-waiting. The formal definition is in Section 3.) This result establishes the large-system insensitivity of LB-zero for Coxian service time distributions with a finite number of phases in the sub-Halfin-Whitt regime. Because the Coxian family is dense in the class of positive-valued distributions, our result justifies load balancing algorithms in the LB-zero family will be able to minimize unnecessary waiting in large-scale data centers for a large class of job size or service time distributions.

To prove this result, this paper develops a new technique, called “iterative state-space peeling” (ISSP). ISSP first identifies an iterative relation between upper and lower bounds on the queue states. Then by iteratively “peeling off” the low-probability states, it proves that the system “lives” near the fixed point of the iterative bounds with a high probability. Based on ISSP, the steady-state distribution of the system can be further analyzed by using Stein’s method in a small neighborhood of the fixed point, where in the large-system regime, the queue lengths under LB-zero behave similarly as an $M/M/\infty$ system because almost all incoming jobs are immediately assigned to an idle server, so we can expect some insensitivity. ISSP, like the state-space collapse in the heavy-traffic analysis, is a general technique that may be used to study other complex stochastic systems, for example, large-system insensitivity of load balancing algorithms for other models such as Mukherjee et al. (2018), Weng and Wang (2020), Weng et al. (2020), and Wang et al. (2021).

We remark that this paper does not establish the large-system insensitivity for an arbitrary service time distribution, for which we need to show that our results continue to hold for a large but finite N when the number of phases of the Coxian distribution goes to infinity. This requires an interchange limit argument or a continuity argument and is an interesting open problem. It is also worth mentioning that a Coxian representation of a probability distribution is nonunique. The choice of the Coxian representation is out of the scope of this paper.

1.2. Related Work

Steady-state analysis of distributed queueing systems has been an active research topic since the seminal work on power-of-two-choices (Mitzenmacher 1996, Vvedenskaya et al. 1996). The most popular approach to study a large-scale distributed queueing system is the mean-field approach where the system is approximated using a deterministic dynamical system (a set of ordinary differential equations), called a mean-field model. In the large-system limit (as the number of servers approaches infinity), the steady state of the stochastic system can often be shown to converge to the equilibrium point of the mean-field model using the interchange of limits (Vvedenskaya et al. 1996, Xie et al. 2015, Ying et al. 2015) or Stein’s method (Ying 2016, Gast 2017).

Although most studies on this topic assume exponential service time distributions for tractability, the approach has been used to study nonexponential service time distributions theoretically or numerically (Mitzenmacher 1996, Gupta et al. 2007, Bramson et al. 2012, Stolyar 2015, Aghajani et al. 2017, Hellemans and Van Houdt 2018, Van Houdt 2019, Vasantam et al. 2019, Kielanski and Van Houdt 2021). Gupta et al. (2007) studied JSQ load balancing with processing sharing service discipline and established its insensitivity property. When the service time distribution has a decreasing hazard rate (DHR), the system often exhibits a monotonicity property such that the system starting from the empty state is dominated by the system starting from any other state. Leveraging this monotonicity, Bramson et al. (2012) proved the convergence of power-of-d-choices (pod), and Stolyar (2015) proved the convergence of JIQ to the corresponding mean-field limit, respectively. Recently, Van Houdt (2019) studied load balancing policies under the hyper-exponential service time distribution in the light traffic regime (or in a critical traffic regime). By transforming the hyper-exponential distribution to a Coxian distribution with DHR, the monotonicity property holds in a partial order and the global stability of the mean-field model was established. Kielanski and Van Houdt (2021) studied Pod with PS servers for a hyper-exponential distribution of order 2 in the light traffic regime. They also established the global stability result in a spirit similar to Van Houdt (2019). Kielanski and Van Houdt (2021) considered Pod where d is a constant independent of the number of servers so the system has a nonvanishing delay in the large-system limit. For the Pod algorithm to be in the LB-zero family, d is a function of N and the algorithm achieves asymptotic zero waiting in the heavy-traffic regime without the DHR assumption.

Without DHR, the results are very limited. Bramson et al. (2012) proved the convergence of join the least loaded of d queues (LL(d)) for general service time distributions and that of Pod when the load of the system is small (less than 1/4). Foss and Stolyar (2017) proved the asymptotic optimality of JIQ under general service time distributions when the normalized load is less than 0.5 (asymptotic optimality means as $N\to \infty ,$ the steady-state probability of an arriving customer being routed to a nonidle server vanishes). Because JIQ is an LB-zero policy, our result confirms the conjecture made in Foss and Stolyar (2017) that JIQ is asymptotically optimal for any load less than one and not just less than 0.5. Another significant result is Bonald et al. (2004), which identifies a set of policies that are insensitive in many-server load-balancing systems and are optimal in the class of insensitive load balancing algorithms. The asymptotic blocking probability of this class of insensitive algorithms in a finite buffer system was later studied in Jonckheere and Prabhu (2018). Many LB-zero such as JSQ and Pod algorithms are sensitive, so the results in Bonald et al. (2004) and Jonckheere and Prabhu (2018) do not apply. The waiting probability in our paper includes both blocking (we consider server systems with finite buffer size and jobs are possible to be blocked) and being queued in the buffer, so our result implies asymptotic zero blocking of LB-zero in the sub-Halfin-Whitt regime.

Liu et al. (2021) is the work most related to this paper, which considers the Coxian-2 distribution and shows that LB-zero achieves asymptotic zero-waiting in the sub-Halfin-Whitt regime. Inspired by Liu et al. (2021), this paper develops the ISSP technique for general Coxian distributions with a finite number of phases and establishes its large-system insensitivity. We remark that Liu et al. (2021) used a key property of Coxian-2 service time distribution that a job in the first phase (phase 1) either departs or enters the last phase (phase 2) immediately, which does not hold under a general Coxian distribution which may have many phases because a job starting in phase 1 requires going through the intermediate phases before entering the last phase. The intermediate phases complicate the dependence among states and make the analysis much more challenging. Technically speaking, it requires a dedicated partition and aggregation of the states for a general Coxian distribution such that the upper/lower bounds on the steady-state distribution can be established in an iterative way. Such an iterative analysis is not necessary in Coxian-2.

In terms of the proof, each step in ISSP uses the tail bound in Bertsimas et al. (2001) to “peel off” a low probability subspace. The tail bound is based on the Lyapunov drift analysis and is a critical step to prove state-space collapse in the traditional heavy traffic regime with a fixed number of servers (Eryilmaz and Srikant 2012, Maguluri and Srikant 2016, Wang et al. 2018). The key difference is that Eryilmaz and Srikant (2012) use the tail bound only once, whereas ISSP repeatedly uses the tail bound guided by an iterative relation between the upper and lower bounds.

2. Model and State Representation

We consider a many-server system with N homogeneous servers, where job arrival follows a Poisson process with rate λN with $\lambda =1-N^{-\alpha },0<\alpha <0.5$, that is, the system is in the sub-Halfin-Whitt regime. We assume the service times follow the Coxian distribution with M phases as shown in Figure 1, where ${\mu }_m>0$ is the rate a job finishes phase m when in service and $0\le p_i<1,1\le i<M$ is the probability that a job enters phase i + 1 after completing phase i and $p_M=0.$ We assume $\lambda =1-N^{-\alpha }$ for the ease of exposition, but our results can be easily extended to the case that $\lambda =1-\beta N^{-\alpha }$ with any positive constant $\beta >0$ independent of N. As convention, we define ${\sum }_{i=a}^b{x}_i=0$ if a > b and ${\prod }_{i=a}^b{x}_i=1$ if a > b for the series $\{x_i,i\ge 1\}.$

Without loss of generality, we normalize the mean service time to be one, that is,

where v_m is viewed as the average time spent in phase m for a job. Given the unit service rate, λ is the normalized load of the system and $\lambda v_m$ is the normalized load of jobs in phase m.

Taking a Coxian-3 distribution as an example (Figure 2), a job is colored in black if it waits in the buffer, and colored in light red, blue, and green when it is in phases 1, 2, and 3, respectively. Jobs are served with the first in first out (FIFO) discipline, and we assume each server has a buffer of size $b-1,$ so can hold at most b jobs (b − 1 in the buffer and one in service). The assumption of finite buffer is imposed due to a technical reason and will be explained later in the paper (specifically, it is required when proving state-space collapse argument in Lemma F.2). Relaxing the finite-buffer assumption is not trivial, but we conjecture that our results hold without this assumption.

To represent the system, define $Q_{j,m}(t)$ ($m=1,2,\dots ,M$) to be the fraction of servers that have j jobs at time t and are serving a job in phase m. Because an idle server does not have a phase, we define $Q_{0,1}(t)$ to be the fraction of servers that are idle at time t and $Q_{0,m}(t)=0,\forall 2\le m\le M$ for convenience. Let $Q(t)\in R^{b\times M}$ be a matrix such that the (j, m)th entry of the matrix is $Q_{j,m}(t).$ We further define $S_{i,m}(t)={\sum }_{j\ge i}{Q}_{j,m}(t)$ and $S_i(t)={\sum }_{m=1}^M{S}_{i,m}(t).$ Therefore, $S_{i,m}(t)$ is the fraction of servers that have at least i jobs and are serving a job in phase m at time t and $S_i(t)$ is the fraction of servers with at least i jobs at time t. Stack $S_{j,m}(t)$ to be a matrix S(t) such that the (j, m)th entry of the matrix is $S_{j,m}(t).$ Because there is a one-to-one mapping between Q(t) and S(t), we focus on S(t) throughout the paper. We consider load balancing policies which dispatch jobs to servers based on S(t) and under which the finite-state continuous-time Markov chain (CTMC) $\{S(t),t\ge 0\}$ is irreducible, and so it has a unique stationary distribution. This includes well-known load balancing policies such as JSQ (Winston 1977, Eschenfeldt and Gamarnik 2018, Braverman 2020), JIQ (Lu et al. 2011, Stolyar 2015), I1F (Gupta and Walton 2019), and Pod (Mitzenmacher 1996, Vvedenskaya et al. 1996).

Let $Q_{j,m}$ be a random variable that has the distribution of $Q_{j,m}(t)$ at steady state. Correspondingly, define $S_{i,m}={\sum }_{j\ge i}{Q}_{j,m}$ and $S_i={\sum }_m{S}_{i,m}.$ In other words, $S_{i,m}$ is the fraction of servers that have at least i jobs and are serving a job in phase m, and S_i is the fraction of servers with at least i jobs, both at steady state. Consider a system with 10 servers and Coxian-3 service time distribution. A realization of state representation $S_{i,m}$ is shown in Figure 3. Define $S\in R^{b\times M}$ to be a matrix such that the (i, m)th entry is $S_{i,m}$ and $s\in {\mathbb{R}}^{b\times M}$ to be a realization of S. Define ${\mathcal{S}}^{(N)}$ to be a set of s as follows:

that is, ${\mathcal{S}}^{(N)}$ is the set of all possible s in a system with N servers. Note ${\mathcal{S}}^{(N)}$ is a set defined for the scaled states.

3. Main Results

We introduce LB-zero load balancing policies and main results in this section. Let $A_1(s)$ be the probability that an incoming job is routed to a busy server given that the system is in state $s\in {\mathcal{S}}^{(N)};$ that is,

We now consider the set of zero-waiting load balancing policies, named as LB-zero,

This class of policies is similar to the one considered in Liu and Ying (2020). Several well-known policies satisfy this condition, as summarized in Table 1.

Table 1. Examples of LB-Zero Policies: Join-the-Shortest-Queue, Join-the-Idle-Queue, Idle-One-First, and Power-of-d-Choices with a Carefully Chosen d

To prove the large-system insensitivity of LB-zero, we first show that $S_{1,m}$ is “close” to $s_{1,m}^{*}=\lambda v_m,$ which is the normalized load from phase m of the jobs and is also the equilibrium point of the mean-field system assuming zero-waiting (details can be found in Sections 4 and 5). We call $s^{*}$ the zero-waiting equilibrium. Theorem 1 shows that at the steady state, $S_{1,m}$ concentrates around the zero-waiting equilibrium $s_{1,m}^{*}$ for large N under any LB-zero policy. The proof of this theorem can be found in Section 6.

Theorem 1.

There exist positive constants θ_m with $1\le m\le M$ independent of N. For any LB-zero policy in $\mathrm{\Pi },$ the following result holds for a large N:

$$\mathbb{P}\left(s_{1,m}^{*}-\frac{{\theta }_m \log N}{\sqrt{N}}\le S_{1,m}\le s_{1,m}^{*}+\frac{1}{N^{\alpha }}-\frac{{\sum }_{r\ne m}{\theta }_r \log N}{\sqrt{N}}\right)\ge 1-\frac{M}{N^3}.$$

Remark 1.

Theorem 1 shows that $S_{1,m}$ differs from $s_{1,m}^{*}$ by at most $\max \left\{\frac{{\theta }_m \log N}{\sqrt{N}},\frac{1}{N^{\alpha }}\right\}$ with a probability at least $1-M/N^3,$ which implies the convergence to the steady-state $S_{1,m},\forall m$ to the zero-waiting equilibrium as $N\to \infty$ in probability and mean-square senses, that is, ($\mathbb{E}[{\Vert S_{1,m}-s_{1,m}^{*}\Vert }^2]\to 0$ as $N\to \infty$). To be best of our knowledge, this is the first result to establish such a steady-state convergence for a load balancing system under Coxian-M service time distributions in the sub-Halfin-Whitt regime. Because the high probability bound holds for a large but finite N, it also provides a rate of convergence. The specific requirement on N and values of constants θ_m can be found in Section 6 and Appendix A, respectively.

From Theorem 1, we have the following corollary on the distance between the steady-state $S_{1,m}$ and its mean-field equilibrium $s_{1,m}^{*}.$

Corollary 1.

When N is sufficiently large, LB-zero policies in Π achieve

$$\mathbb{P}\left(|S_{1,m}-s_{1,m}^{*}|\le \frac{1}{N^{\alpha }}\right)\ge 1-\frac{M}{N^3} \mathit{\text{and}} \mathbb{E}\left[|S_{1,m}-s_{1,m}^{*}|\right]=O\left(\frac{1}{N^{\alpha }}\right).$$

From Theorem 1 and Corollary 1, it is not clear whether the probability of waiting approaches zero under an LB-zero policy because the average number of waiting jobs (per server) ${\sum }_{i=2}^b{S}_i$ is not clear yet. This will be studied in the next theorem and its key step is to establish the upper bound on the average queue length $\mathbb{E}[{\sum }_{m=1}^b{S}_i]$ in Lemma 6. Let $\mathcal{W}$ denote the event that an incoming job is routed to a busy server in the system, and $\mathbb{P}(\mathcal{W})$ denote the probability of this event at steady state. We have the following result on the waiting probability. The proof of Theorem 2 can be found in Section 7.

Theorem 2.

For any LB-zero policy in $\mathrm{\Pi },$ the following result holds for a large $N:$

$$\mathbb{P}(\mathcal{W})\le \frac{1}{\sqrt{N}}+\frac{10{\mu }_{\text{max}}+4}{N^{0.5-\alpha } \log N}.$$

Remark 2.

Theorem 2 shows that the waiting probability is $O(1/N^{0.5-\alpha })$ for a large but finite N, which implies the asymptotic zero waiting, that is, $\mathbb{P}(\mathcal{W})\to 0$ as $N\to \infty$ in the sub-Halfin-Whitt regime. This asymptotic result implies LB-zero is large-system insensitive to Coxian-M distributions with a finite number of phases. The specific requirement on N can be found in Section 7.

Next, we will establish these two main results. We first introduce the system dynamic of LB-zero in Section 4 and present ISSP in Section 5, which is used to prove that the system lives near a limiting regime in Theorem 1 in Section 6 and to prove the large system insensitivity in Theorem 2 in Section 7. A roadmap can be found in Figure 4, which demonstrates the relationship of the key lemmas and theorems. The condition in LB-zero is used when we analyze the dynamics of $S_{1,1}$ or $S_1,$ including Lemma 1 in ISSP (the lower bound on $S_{1,1}$) and “zero-waiting” in Theorem 2 and its associated intermediate Lemma F.2.

4. System Dynamics

Define $e_{i,m}\in {\mathbb{R}}^{b\times M}$ to be a b × M-dimensional matrix with the (i, m)th entry being $1/N$ and all other entries being zero. Furthermore, define $A_{i,m}(s)$ to be the probability that an incoming job is routed to a server with at least i jobs and the job in service is in phase m given the system state s, that is,

Given state s (or the corresponding q) of the CTMC, each of the following three events triggers a state transition, which is illustrated individually in Figure 5.

• Event 1: A job arrives and is routed to a server that has i − 1 jobs and the job in service is in phase m as in the left figure in Figure 5. When this occurs, $q_{i,m}$ increases by $1/N,$ and $q_{i-1,m}$ decreases by $1/N$ (note m = 1 if i = 1 since we define the fraction of idle servers to be $q_{0,1}$). Therefore, the CTMC has the following transition:

  $$\begin{array}{l}q\to q+e_{i,m}-e_{i-1,m},\hfill \\ s\to s+e_{i,m},\hfill \end{array}$$
  where the transition of s can be verified according to the definition $s_{i,m}={\sum }_{j\ge i}{q}_{j,m}$ so only $s_{i,m}$ increasing by $1/N.$ This event occurs with rate
  $$\lambda N(A_{i-1,m}(s)-A_{i,m}(s)),$$
  where $A_{i-1,m}(s)-A_{i,m}(s)$ is the probability that an incoming job is routed to a server that has i − 1 jobs and the job in service is in phase m.

• Event 2: A server with i jobs finishes serving a job in phase m, and the job departs from the system without entering into the next phase as in the middle figure in Figure 5. When this event occurs, $q_{i,m}$ decreases by $1/N$ and $q_{i-1,1}$ increases by $1/N,$ so the CTMC has the following transition:

  $$\begin{array}{l}q\to q-e_{i,m}+e_{i-1,1},\hfill \\ s\to s-{\displaystyle \sum _{j=1}^i{e}_{j,m}}+{\displaystyle \sum _{j=1}^{i-1}{e}_{j,1}},\hfill \end{array}$$
  where the transition of s can be verified based on the definition $s_{i,m}={\sum }_{j\ge i}{q}_{j,m}$ so $s_{j,m}$ decreases by $1/N$ for any $j\le i$ and $s_{j,1}$ increases by $1/N$ for any $j<i.$ This event occurs with rate
  $${\mu }_{m}N{q}_{i,m}(1-p_m),$$
  where ${\mu }_{m}N{q}_{i,m}$ is the rate at which a job in phase m finishes the service and $(1-p_m)$ is the probability that a job finishes phase m and departs from the system immediately.

• Event 3: A server with i jobs finishes serving a job in phase m, and the job enters the next phase m + 1 as shown in the right figure in Figure 5. When this event occurs, a server in state (i, m) transits to state $(i,m+1),$ so $q_{i,m}$ decreases by $1/N$ and $q_{i,m+1}$ increases by $1/N.$ Therefore, the CTMC has the following transition:

  $$\begin{array}{l}q\to q-e_{i,m}+e_{i,m+1},\hfill \\ s\to s-{\displaystyle \sum _{j=1}^i{e}_{j,m}}+{\displaystyle \sum _{j=1}^i{e}_{j,m+1}},\hfill \end{array}$$
  where the transition of s holds because $s_{i,m}$ decreases by $1/N$ for any $j\le i$ and $s_{j,m+1}$ increases by $1/N$ for any $j\le i.$ This event occurs with rate
  $${\mu }_{m}N{q}_{i,m}{p}_m,$$
  where ${\mu }_{m}N{q}_{i,m}$ is the rate at which a job in phase m finishes the service and p_m is the probability that a job enters phase m + 1 after finishing phase m.

Define G to be the generator of CTMC $(S(t):t\ge 0).$ Given function $f:{\mathcal{S}}^{(N)}\to \mathbb{R},$ we have

To understand the dynamics better, we write down the mean-field model (MFM) according to the generator:

This mean-field model is nonlinear in s because $A_{i,m}(s)$ is a nonlinear function in s ($A_{i,m}(s)=0,\forall i\le 0.$) and its equilibrium point is difficult to calculate in general. However, suppose zero-waiting occurs, that is, $s_{i,m}=0,\forall i\ge 2,$ and the fraction of jobs dropped is negligible, then we can obtain the following equilibrium:

We call $s^{*}$ zero-waiting equilibrium because it is a conjectured equilibrium by assuming zero-waiting. Note zero-waiting equilibrium corresponds to the mean-field model in (5)–(6) with only $A_{0,1}=1$ and $A_{i,m}=0$ for any $(i,m)\ne (0,1)$. In this following analysis, we will not solve Mean-Field Model (5)–(6) to check whether its equilibrium is close to $s^{*}$. Instead, we will directly prove $S_{1,m}$ concentrates around $s_{1,m}^{*}$ and zero waiting occurs at the steady state with a high probability.

5. ISSP

In this section, we illustrate the key idea of ISSP, which will be applied to prove Theorem 1. Intuitively, the original stochastic system $S_{i,m}(t)$ and the steady-state $S_{i,m}$ would be close to the MFM $s_{i,m}(t)$ and the zero-waiting equilibrium $s_{i,m}^{*},$ respectively. However, because of the “nonmonotonicity” of the system (the “monotonicity” means a partial order of two systems starting from two initial conditions will be maintained over time), it is extremely challenging to justify this argument. To tackle the challenge, we develop a new technique, called ISSP, which first identifies an iterative relation between upper and lower bounds on the original stochastic queue states and then proves that the system lives in a regime concentrated around the “fixed point” of the iterative bounds with a high probability. We outline the general steps to apply ISSP. The key step is to identify the lower and upper bound functions according to the drift (or dynamics) of the mean-field model such that if the system is in region 1, then it leaves region 1 and moves toward region 2 with a rate (or drift) of x. The drift will provide a high probability upper or lower bound on the steady-state distribution. In particular, if the mean-field system leaves a region with a certain rate, then the probability that the steady-state system is in this region is low.

In particular, we focus on $S_{1,m}$ the fraction of busy servers with the job in service in phase m because we hypothesize $S_{i,m}\to 0,\forall i\ge 2$ (in other words, the fraction of servers with any waiting jobs is negligible in a large system).

Let n be the nth step of ISSP. Let $L_{1,m}(n)$ denote a high probability lower bound on $S_{1,m}$ and $U_m(n)$ be a high probability upper bound on ${\sum }_{r=2}^m{S}_{1,r},$ established at the nth step of ISSP, that is,

are close to one. Our goal is to show that as n increases, $L_{1,m}(n)$ and $U_m(n)$ approach the zero-waiting equilibrium $s_{1,m}^{*}$ and ${\sum }_{r=2}^m{s}_{1,r}^{*}$. Taking a Coxian-3 distribution as an example, we need to show that as n increases $L_{1,m}(n)\to s_{1,m}^{*},\forall m,$ $U_2\to s_{1,2}^{*}$ and $U_3\to s_{1,2}^{*}+s_{1,3}^{*}.$

5.1. ISSP: Hand-Waving Explanation

In this section, we present a hand-waving explanation of our ISSP, which is based on the following iterative relation between the upper and lower bounds on the queue states, denoted by $L_{1,m}$ and $U_m,\forall m\ge 2:$

where the initial condition $L_{1,m}(0)=0,\forall m,$ and $U_m(0)=1,\forall m\ge 2,$ and “$\approx$” is used because we ignore diminishing terms $($e.g., $\frac{\log N}{\sqrt{N}})$ in Equations (7)–(9) for explaining the intuition.

We next provide a hand-waiving argument to show these equations iteratively provide upper and lower bounds. The exact definitions of $L_{1,m}$ and $U_m,\forall m\ge 2$ and their relations can be found in (15)–(17). The formal proofs that establish these are indeed high-probability upper and lower bounds can be found in Section 6. The right-hand sides of (7)–(9) can be regarded as the iterative update to refine the upper and lower bounds. Considering (7), given ${\sum }_{r=2}^m{S}_{1,r}\le U_M(n)$ holds with a high probability, we can establish $S_{1,1}\ge L_{1,1}(n+1)$ holds with a high-probability and the formula of $L_{1,1}(n+1)$ is given in the right-hand side in (7). Considering (8), given $S_{1,m-1}\ge L_{1,m-1}(n+1)$ holds with a high-probability, we can establish $S_{1,m}\ge L_{1,m}(n+1)$ holds with a high probability, where the formula of $L_{1,1}(n+1)$ is given in the right-hand side in (8). Considering (9), given $S_{1,1}\ge L_{1,1}(n+1)$ and ${\sum }_{r=2}^{m-1}{S}_{1,r}\le U_{1,m}(n+1)$ holds with a high probability, we can establish ${\sum }_{r=2}^m{S}_{1,r}\le U_{1,m}(n+1)$ with a high probability, where the formula of $L_{1,1}(n+1)$ is given in the right-hand side in (9).

From (7)–(9), we can obtain a recursive equation for $L_{1,1}:$

Therefore, as $n\to \infty ,$ $L_{1,m}(n)\to s_{1,m}^{*}$ and $U_m(n)\to {\sum }_{r=2}^m{s}_{1,r}^{*}.$

To provide the intuition behind (7)–(9), we consider the mean-field model under JSQ as an example and focus on $s_{1,m}$ in (5)–(6) by ignoring $s_{i,m},\forall i\ge 2,$ that is,

where an equilibrium can be verified to be $s_{1,m}^{*}=\lambda v_m,\forall m.$

We next carefully analyze (11)–(12) to establish (7)–(9). To derive the lower and upper bounds on the equilibrium point of a dynamical system x(t) (or $s_{1,m}(t)$), we use the following straightforward ideas. We need to emphasize that these are intuitive ideas and informal arguments to guide us to formalize ISSP in Section 6.

• If

  $$\dot{x}(t)>L-x(t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{with} x(0)<L,$$(13)
  then x(t) eventually becomes greater than or equal to L, that is, $x(\infty )\ge L,$ because otherwise x(t) continues to increase.

• If

  $$\dot{x}(t)<U-x(t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{with} x(0)>U,$$(14)
  then x(t) eventually becomes smaller than or equal to U, that is, $x(\infty )\le U,$ because otherwise, x(t) continues to decrease.

In the following, we will explain (7)–(9) based on the previous ideas. The explanation is not a rigorous proof. The detailed proof will be presented later in Section 6. We will ignore the iteration index n occasionally when confusion does not arise.

The intuition to obtain (7): We start with the dynamic of $s_{1,1}$ in (11), which is

given ${\sum }_{m=2}^M{s}_{1,m}<U_M,$ we havewhen $s_{1,1}<1-U_M,$ which implies ${\dot{s}}_{1,1}>0$ when $s_{1,1}<\min \{1-U_M,s_{1,1}^{*}\}.$ Therefore, at the equilibrium point, $s_{1,1}\ge L_{1,1}\triangleq \min \{1-U_M,s_{1,1}^{*}\}$ because otherwise, $s_{1,1}$ will continue to increase because ${\dot{s}}_{1,1}>0.$

The intuition to obtain (8): Consider the dynamic of $s_{1,m},\forall m\ge 2$ in (12):

Given $s_{1,m-1}\ge L_{1,m-1},$ we have

which implies at the equilibrium point,because otherwise, $s_{1,m}$ will continue to increase because ${\dot{s}}_{1,m}>0.$

Note $\frac{v_m}{v_{m-1}}=\frac{s_{1,m}^{*}}{s_{1,m-1}^{*}},$ so we have $\frac{L_{1,m}}{L_{1,m-1}}=\frac{s_{1,m}^{*}}{s_{1,m-1}^{*}},$ which means that the ratio of the lower bounds is the same as that of the corresponding equilibrium points.

The intuition to obtain (9): We focus on the dynamic of ${\sum }_{r=2}^m{s}_{1,r}$ that

Given $s_{1,m}\ge L_{1,m},\forall m$ and ${\sum }_{r=2}^{m-1}{s}_{1,r}\le U_{m-1},$ we have

which implies at the equilibrium point,because otherwise, ${\sum }_{r=2}^m{s}_{1,r}$ will continue to decrease because ${\sum }_{r=2}^m{\dot{s}}_{1,r}<0.$

By invoking $L_{1,m}=\frac{v_m}{v_{m-1}}{L}_{1,m-1},$ we eventually have

5.2. ISSP: Illustrative Example

To demonstrate ISSP, we consider JSQ with Erlang-3 distribution and no buffer, that is, Erlang-3 with $b=1.$

The mean-field model under JSQ with Erlang-3 and $b=1.$

With the Erlang-3 service time distribution, we have

Therefore, the MFM in this case is

ISSP for JSQ with Erlang-3 and $b=1.$ The values of the key parameters in iterative equations in this case are $a_1=a_2=a_3=\frac{1}{2},$ $b_1=1,b_2=\frac{1}{2},b_3=0,$ and $v_1=v_2=v_3=\frac{1}{3}.$ The corresponding iterative equations are

The iterative relation in terms of $L_{1,1}(n)$ is

which implies that $L_{1,1}(n)\to \frac{\lambda }{3}$ as $n\to \infty ,$ and

The iterative procedure has been visualized in Figure 6. In each iteration n, we first establish $L_{1,1}(n+1)$ (shown in light red) based on $U_3(n)$ from the last iteration, and then obtain $L_{1,2}(n+1)\approx L_{1,1}(n+1)$ (shown in light blue) based on $L_{1,1}(n+1);$ and $L_{1,3}(n+1)\approx L_{1,2}(n+1)$ (shown in light green) based on $L_{1,2}(n+1).$ Given the new lower bounds on $L_{1,m}(n+1),$ we first refine $U_2(n+1)$ (shown in light purple) and finally establish $U_3(n+1)$ given all the new lower bounds $L_{1,m}(n+1)$ and upper bound $U_2(n+1),$ which in turn will improve $L_{1,1}$ in the next iteration.

In the following sections, we formalize the ISSP, which is then combined with Stein’s method to prove Theorems 1 and 2.

6. Proof of Theorem 1 Based on ISSP

In this section, we formalize ISSP to prove Theorem 1. We first present the formal definitions of the iterative Equations (7)–(9) in the following:

These iterative lower and upper bounds are proved to hold with high probability in the following lemmas, that is,

The probability $\{{ϵ}_m(n),{\sigma }_m(n)\}$ are diminishing quantities for a large N. The specific values of $\{{ϵ}_m(n),{\sigma }_m(n)\}$ and positive constants $\{C,\xi ,a_m,b_m,c_m\}$ can be found in Appendix A. We only state the lemmas and their proofs can be found in Appendix C.

A lower bound on $L_{1,1}(n+1)$ given ${\sum }_{m=2}^M{S}_{1,m}\le U_M(n):$ The following lemma is the rigorous statement of (7).

Lemma 1.

Given

$$\mathbb{P}\left({\displaystyle \sum _{m=2}^M{S}_{1,m}}>U_M(n)\right)\le {\sigma }_M(n),$$

we have

$$\mathbb{P}\left(S_{1,1}<L_{1,1}(n+1)\right)\le {ϵ}_1(n+1),$$

where

$$L_{1,1}(n+1)=\min \left\{s_{1,1}^{*}-\frac{6\mathrm{\Delta }}{C},1-U_M(n)-\frac{1-\xi }{2{\mu }_1{N}^{\alpha }}-\frac{6\mathrm{\Delta }}{C}\right\}.$$

A lower bound on $S_{1,m}$ given $S_{1,m-1}\ge L_{1,m-1}$: The following lemma is the rigorous statement of (8).

Lemma 2.

Consider $m\ge 2.$ Given

$$\mathbb{P}\left(S_{1,m-1}<L_{1,m-1}(n+1)\right)\le {ϵ}_{m-1}(n+1),$$

we have

$$\mathbb{P}\left(S_{1,m}<L_{1,m}(n+1)\right)\le {ϵ}_m(n+1),$$

where

$$L_{1,m}(n+1)=\frac{v_m}{v_{m-1}}{L}_{1,m-1}(n+1)-\frac{5v_m}{C}\mathrm{\Delta }.$$

An upper bound on ${\sum }_{r=2}^m{S}_{1,r}$ given ${\sum }_{r=2}^{m-1}{S}_{1,r}\le U_{m-1}$ and $S_{1,m-1}\ge L_{1,m-1}:$ The following lemma is the rigorous statement of (9).

Lemma 3.

Consider $m\ge 2.$ Given

$$\mathbb{P}\left({\displaystyle \sum _{r=2}^{m-1}{S}_{1,r}}\ge U_{m-1}(n+1)\right)\le {\sigma }_{m-1}(n+1),$$

$$\mathbb{P}\left(S_{1,r}<L_{1,r}(n+1)\right)\le {ϵ}_r(n+1)\hspace{1em}1\le r\le M,$$

we have

$$\mathbb{P}\left({\displaystyle \sum _{r=2}^m{S}_{1,r}}\ge U_m(n+1)\right)\le {\sigma }_m(n+1),$$

where

$$U_m(n+1)=1-a_m-b_m{L}_{1,1}(n+1)+a_m{U}_{m-1}(n+1)+\frac{c_m}{C}\mathrm{\Delta }.$$

6.1. Convergence of $L_{1,1}(n).$

Based on Lemmas 1–3, we will show ${\{L_{1,1}(n)\}}_n$ is an increasing sequence and approaches $s_{1,1}^{*}.$

Lemma 4.

Given $\mathbb{P}(S_{1,1}<L_{1,1}(n))\le {ϵ}_1(n)$ and $L_{1,1}(n)\le s_{1,1}^{*}-\frac{6\mathrm{\Delta }}{C},$ we have

$$\mathbb{P}(S_{1,1}<L_{1,1}(n+1))\le {ϵ}_1(n+1),$$

where

$$L_{1,1}(n+1)=\frac{1}{{\mu }_1}-\frac{C_M\mathrm{\Delta }}{C(1-\xi )}-\frac{1}{2{\mu }_1{N}^{\alpha }}+\xi \left(L_{1,1}(n)-\frac{1}{{\mu }_1}+\frac{C_M\mathrm{\Delta }}{C(1-\xi )}+\frac{1}{2{\mu }_1{N}^{\alpha }}\right)$$

and ${ϵ}_1(n+1)=\max \{{ϵ}_1(n)(M^2+2){\left(\frac{2C}{\overline{v}\mathrm{\Delta }}+1\right)}^{2M},e^{-\frac{{\overline{v}}^2 {\log }^2 N}{C^2}}\}.$ Furthermore, $L_{1,1}(n+1)>L_{1,1}(n)$ holds when $L_{1,1}(n)\le s_{1,1}^{*}-\frac{6\mathrm{\Delta }}{C}.$

6.2. Proving Theorem 1

Because $\{L_{1,1}(n)\}$ is an increasing sequence before reaching its equilibrium, we can apply Lemma 4 a sufficient number of times so that $L_{1,1}(n)$ is close to $s_{1,1}^{*}-\frac{6\mathrm{\Delta }}{C},$ which is formalized in the following lemma. We can keep applying Lemma 4 because both the lower bound $L_{1,1}(n+1)$ and the probability ${ϵ}_1(n+1)$ only depend on the previous iteration $L_{1,1}(n)$ and ${ϵ}_1(n)$.

Lemma 5.

$$\mathbb{P}\left(S_{1,1}<s_{1,1}^{*}-\frac{6\mathrm{\Delta }}{C}\right)\le {\left(\frac{1}{\sqrt{N}}\right)}^{M+8}.$$

Proof.

To prove this lemma, we apply Lemma 4 n times iteratively with $n=\lceil \frac{\log N}{2 \log (1/\xi )}\rceil$ such that ${\xi }^n\le \mathrm{\Delta }.$ We obtain

$$\begin{array}{l}\mathbb{P}\left(S_{1,1}<\frac{1}{{\mu }_1}-\frac{C_M\mathrm{\Delta }}{(1-\xi )C}-\frac{1}{2{\mu }_1{N}^{\alpha }}-\frac{\mathrm{\Delta }}{{\mu }_1}\right)\le {ϵ}_1(1){(M^2+2)}^n{\left(\frac{2C}{\overline{v}\mathrm{\Delta }}+1\right)}^{2Mn}\hfill \\ \le {ϵ}_1(1){(M^2+2)}^{\frac{\log N}{\log (1/\xi )}+1}{N}^{\frac{2M \log N}{\log (1/\xi )}+2M}\hfill \\ \le e^{-\frac{{\overline{v}}^2 {\log }^2 N}{C^2}}{(M^2+2)}^{\frac{\log N}{\log (1/\xi )}+1}{N}^{\frac{2M \log N}{\log (1/\xi )}+2M}\hfill \\ =N^{-\frac{{\overline{v}}^2 \log N}{C^2}+\frac{2M \log N}{\log (1/\xi )}+2M}{(M^2+2)}^{\frac{\log N}{\log (1/\xi )}+1}\hfill \\ \le N^{-\frac{{\overline{v}}^2 \log N}{C^2}+\frac{3M \log N}{\log (1/\xi )}+2M+1},\hfill \end{array}$$

where the third inequality holds because ${ϵ}_1(n+1)=\max \{{ϵ}_1(n)(M^2+2){\left(\frac{2C}{\overline{v}\mathrm{\Delta }}+1\right)}^{2M},e^{-\frac{{\overline{v}}^2 {\log }^2 N}{C^2}}\}$ in Lemma 4 implies ${ϵ}_1(1)\le e^{-\frac{{\overline{v}}^2 {\log }^{2}N}{C^2}}$ and the last inequality holds because $N\ge M^2+2.$

Recalling $C=\sqrt{\frac{2{\overline{v}}^2 \log (1/\xi )}{3M+(3M+4) \log (1/\xi )}}$ and noting $N\ge 1/\xi ,$ we have

$$\begin{array}{l}\mathbb{P}\left(S_{1,1}\le \frac{1}{{\mu }_1}-\frac{C_M\mathrm{\Delta }}{(1-\xi )C}-\frac{\mathrm{\Delta }}{{\mu }_1}-\frac{1}{2{\mu }_1{N}^{\alpha }}\right)\hfill \\ =\mathbb{P}\left(S_{1,1}\le \frac{\lambda }{{\mu }_1}-\frac{C_M\mathrm{\Delta }}{(1-\xi )C}-\frac{\mathrm{\Delta }}{{\mu }_1}+\frac{1}{2{\mu }_1{N}^{\alpha }}\right)\le {\left(\frac{1}{\sqrt{N}}\right)}^{M+8},\hfill \end{array}$$

which implies

$$\mathbb{P}\left(S_{1,1}<s_{1,m}^{*}-\frac{6\mathrm{\Delta }}{C}\right)\le {\left(\frac{1}{\sqrt{N}}\right)}^{M+8},$$

because $N^{0.5-\alpha }\ge \frac{2{\mu }_1{C}_M}{C(1-\xi )}\log N.$ □

The result in Lemma 5 established that

with probability $1-{\left(\frac{1}{\sqrt{N}}\right)}^{M+8}.$

Combining with Lemma 2, we next prove that $S_{1,m}\ge s_{1,m}^{*}-\mathrm{\Omega }(\mathrm{\Delta }),\forall m\ge 2$ holds with a high probability. Applying Lemma 2 iteratively for $S_{1,m},\forall m\ge 2,$ we have

and $S_{1,m}\ge L_{1,m}(n)$ holds with the probability ${ϵ}_m(n).$ Note ${ϵ}_m(n)\le {ϵ}_M(n),\forall m$ and ${ϵ}_M(n)$ is bounded as follows:

Therefore, we have proved the lower bound in the theorem.

Define the event $\mathcal{K}=\{S_{1,m}\ge s_{1,m}^{*}+\frac{{\theta }_m \log N}{\sqrt{N}},\forall m\}.$ We have $\mathbb{P}({\mathcal{K}}^c)\le \frac{M}{N^3}$ according to the union bound. We now establish the upper bound in Theorem 1 as follows:

The proof is completed because $1-{\sum }_{r\ne m}^M{s}_{1,r}^{*}=s_{1,m}^{*}+\frac{1}{N^{\alpha }}.$

7. Proof of Theorem 2

Theorem 1 shows that $S_{1,m}$ is “close” to $s_{1,m}^{*}$ with a high probability. However, it is not clear whether the average total queue length or the waiting probability is small under an LB-zero policy because Theorem 2 alone does not quantify the expected total queue length $\mathbb{E}[{\sum }_{m=1}^b{S}_i]$ and the waiting probability. To establish Theorem 2, we first prove an important lemma on the upper bound of the average total queue, which is used to establish the waiting probability in Theorem 2 by the Markov inequality. The proof of this lemma can be found in Sections 7.1–7.3.

Lemma 6.

Define $w_m=(1-p_m){\mu }_m,$ $w_u={\max }_m{w}_m,$ $w_l={\min }_m{w}_m,$ ${\mu }_{\text{max}}={\max }_m{\mu }_m,$ $\zeta =\frac{4w_{u}b}{w_l}((\frac{1}{w_l}-\frac{1}{w_u}){\sum }_m{\theta }_m{w}_m+\frac{1}{w_l}+6)$, and $k=\frac{{\sum }_m{\theta }_m{w}_m}{w_u}+(1+\frac{w_l}{4w_{u}b})\zeta -{\sum }_m{\theta }_m.$ Under a load balancing policy in LB-zero, the following bound holds:

$$\mathbb{E}\left[\max \left\{{\displaystyle \sum _{i=1}^b{S}_i}-\lambda -\frac{k \log N}{\sqrt{N}},0\right\}\right]\le \frac{5{\mu }_{\text{max}}+2}{\sqrt{N} \log N},$$

when N satisfies

$$\min \left\{2k{\mu }_1,{\displaystyle \sum _{m=1}^M{\theta }_m},\frac{C(1-\xi )}{2{\mu }_1{C}_M}\right\}{N}^{0.5-\alpha }\ge \log N\ge \max \left\{\log \left(\frac{1}{\xi }\right),\frac{2{\mu }_1}{1-\xi },\frac{4b}{w_l\zeta },\frac{C}{{\mu }_1}\right\}.$$

Lemma 6 establishes an upper bound on the average total queue length. Recall Theorem 1 indicates the service rate under an LB-zero policy in Π is close to arrival rate λN at steady state because ${\sum }_{m=1}^M(1-p_m)S_{1,m}\approx {\sum }_{m=1}^M(1-p_m)s_{1,m}^{*}=\lambda .$ Therefore, it is reasonable to couple the distributed load balancing system with a simple centralized server system with a similar arrival rate and service rate. This coupling will be done via Stein’s method (Braverman et al. 2016, Ying 2016). Stein’s method allows us to understand the key performance metrics (e.g., average queue length) of a complicated load balancing system from the performance of a simple fluid system (to be introduced in the next section). Formally, we study the generator difference between the distributed load balancing system and a simple centralized system within a small state space identified by ISSP. This idea of coupling a simple (and almost trivial) fluid model, without ISSP, has been also used in Liu and Ying (2020, 2021) and Liu et al. (2021). We will introduce it next so the paper is self-contained. Generator coupling in Section 7.1 and state space peeling in Section 7.2 are used to prove Lemma 6.

7.1. Generator Coupling with a Single Server System

Denote $\mathrm{\Delta }=\frac{\log N}{\sqrt{N}}.$ We consider a single server queue with arrival rate λ and service rate $\lambda +\mathrm{\Delta }.$ The fluid model with respect to the queue length x is

Let function g(x) be the solution of the following Stein’s equation or Poisson equation (Ying 2016) for the fluid system above and distance function h(x) such that

where $g\prime (x)=\frac{dg(x)}{dx}$. Because (19) has a very simple form, both $g\prime$ and $g″$ can be easily solved, which is different from other applications of Stein’s method, for example, Braverman et al. (2016), where establishing the gradient bounds is a key difficulty.

To analyze the total queue length at steady state under an LB-zero policy in $\mathrm{\Pi },$ we choose a truncated distance function:

where k is a positive constant. The distance $h({\sum }_{i=1}^b{S}_i)$ can be viewed as a proxy to measure the total queue length $(N{\sum }_{i=1}^b{S}_i)$ at steady state.

To couple the one-dimensional fluid system in (18) with the b × M-dimensional stochastic system, we define

Note f(s) is bounded for $s\in {\mathcal{S}}^{(N)},$ and we imply the generator of stochastic system G defined in (2)–(4) on function f and have the basic adjoint relationship for the stationary distribution S such that

Combining (19) and (21), we connect the performance metric $h(·)$ with the generator difference between the simple single-server system and G as follows:

In the following lemma, we provide an upper bound on (22) that includes two terms: the first term is from the gradient bounds and the second term is from ISSP. The proof of this lemma can be found in Appendix D.

Lemma 7.

Define the regions $\mathcal{T}=\left\{x|x>\eta +\frac{1}{N}\right\}$ and denote the normalized service rate $D_1={\sum }_{m=1}^M(1-p_m){\mu }_m{S}_{1,m},$ then we have

$$\mathbb{E}\left[h\left({\displaystyle \sum _{i=1}^b{S}_i}\right)\right]\le J_1+\frac{5{\mu }_{\text{max}}+\lambda }{\sqrt{N} \log N},$$(23)

with

$$J_1=\mathbb{E}\left[g\prime \left({\displaystyle \sum _{i=1}^b{S}_i}\right)\left(\lambda A_b(S)-\lambda -\mathrm{\Delta }+D_1\right){\mathbb{I}}_{{\sum }_{i=1}^b{S}_i\in \mathcal{T}}\right].$$(24)

To establish Lemma 6 and Theorem 2 based on the previous lemma, we need to provide the upper bounds on (24), which are related to the difference between the normalized arrival rate and the normalized service rate, which will be bounded based on the ISSP result that shows the normalized service rate at the steady state is close to the zero-waiting equilibrium value.

7.2. State Space Peeling on ${\sum }_{m=1}^M(1-p_m)S_{1,m}$

We analyze J₁ in (24) by splitting it into two cases and show it is $O\left(\frac{1}{N^{1.5} \log N}\right)$ in Section 7.3:

where the equality is due to the definition of $g\prime$ in Stein’s Equation (19), and the inequality holds because $\frac{1}{\mathrm{\Delta }}h\left({\sum }_{i=1}^b{S}_i\right){\mathbb{I}}_{{\sum }_{i=1}^b{S}_i>\eta +\frac{1}{N}}\ge 0.$ We focus onwhere we recall $\eta =\lambda +k\mathrm{\Delta }$, and $d_1={\sum }_{m=1}^M(1-p_m){\mu }_m{s}_{1,m}$ is the total service rate when the system is in the state s. Although we have established $S_{1,m}\approx s_{1,m}^{*}$ in Theorem 1, we might have a trivial upper bound of (25), that is, $\mathbb{E}[h\left({\sum }_{i=1}^b{S}_i\right){\mathbb{I}}_{{\sum }_{i=1}^b{S}_i>\eta +\frac{1}{N}}],$ because $\lambda +\mathrm{\Delta }-D_1$ could be greater than Δ. In fact, we need one more state space peeling to show (25) is nonpositive.

We define two regions ${\mathcal{S}}_{{\mathit{\text{ssp}}}_1}$ and ${\mathcal{S}}_{{\mathit{\text{ssp}}}_2}$:

where ${\mathcal{S}}_{{\mathit{\text{ssp}}}_1}$ is the region with sufficient many busy servers and ${\mathcal{S}}_{{\mathit{\text{ssp}}}_2}$ is the region with bounded total queue length. We further define a regionand consider two cases: $s\in {\mathcal{S}}_{\mathit{\text{ssp}}}$ and $s\notin {\mathcal{S}}_{\mathit{\text{ssp}}},$ where

• Case 1: In Lemma E.1, we show any $s\in {\mathcal{S}}_{{\mathit{\text{ssp}}}_1}$ satisfies

  $${\displaystyle \sum _{m=1}^M(1-p_m)}{\mu }_m{s}_{1,m}\ge \lambda +\mathrm{\Delta }.$$

  It implies $(\lambda +\mathrm{\Delta }-{\sum }_{m=1}^M(1-p_m){\mu }_m{s}_{1,m}){\mathbb{I}}_{{\sum }_{i=1}^b{s}_i>\eta +\frac{1}{N}}\le 0$ for any $s\in {\mathcal{S}}_{{\mathit{\text{ssp}}}_1}.$ For any $s\in {\mathcal{S}}_{{\mathit{\text{ssp}}}_2},$ we have

  $${\mathbb{I}}_{{\sum }_{i=1}^b{s}_i>\eta +\frac{1}{N}}=0.$$

  It implies $(\lambda +\mathrm{\Delta }-{\sum }_{m=1}^M(1-p_m){\mu }_m{s}_{1,m}){\mathbb{I}}_{{\sum }_{i=1}^b{s}_i>\eta +\frac{1}{N}}=0$ for any $s\in {\mathcal{S}}_{{\mathit{\text{ssp}}}_2}.$

• Case 2: In Lemma F.1, we show that

  $$\mathbb{P}\left(S\notin {\mathcal{S}}_{\mathit{\text{ssp}}}\right)\le \frac{2}{N^2}$$
  using an ISSP approach on S₁ and ${\sum }_{i=2}^b{S}_i.$

Therefore, (25) are nonpositive under $s\in {\mathcal{S}}_{\mathit{\text{ssp}}}$ and $O(\frac{1}{N^{1.5} \log N})$ under $s\notin {\mathcal{S}}_{\mathit{\text{ssp}}}$, respectively, as shown in Section 7.3.

7.3. Proving Lemma 6

Recall (25) $=\mathbb{E}[\frac{1}{\mathrm{\Delta }}h\left({\sum }_{i=1}^b{S}_i\right)\left(\lambda +\mathrm{\Delta }-D_1\right){\mathbb{I}}_{{\sum }_{i=1}^b{S}_i>\eta +\frac{1}{N}}].$ Based on the two regions $s\in S_{\mathit{\text{ssp}}}$ and $s\notin S_{\mathit{\text{ssp}}}$, we split (25) into these two regions as follows:

The first item in (27) is negative because Lemma E.1 shows $\lambda +\mathrm{\Delta }-D_1\le 0$ for any $S\in {\mathcal{S}}_{\mathit{\text{ssp}}1}$ and $\left(\lambda +\mathrm{\Delta }-D_1\right){\mathbb{I}}_{{\sum }_{i=1}^b{S}_i>\eta +\frac{1}{N}}=0$ for any $S\in {\mathcal{S}}_{{\mathit{\text{ssp}}}_2}.$ Therefore, we have

where the last inequality in (28) holds because the average total number of jobs per server is at most b and $\left(\lambda +\mathrm{\Delta }-D_1\right){\mathbb{I}}_{S\notin {\mathcal{S}}_{\mathit{\text{ssp}}}}{\mathbb{I}}_{{\sum }_{i=1}^b{S}_i>\eta +\frac{1}{N}}<1.$

By combining (23) and (28), we can now establish Lemma 6 in the following:

7.4. Proving Theorem 2

Once we have Lemma 6, we can prove Theorem 2 with the property of LB-zero and the Markov inequality. For an LB-zero policy in $\mathrm{\Pi },$ the waiting probability satisfies

where the first term is bounded by $\frac{1}{\sqrt{N}}$ because of the definition of LB-zero. For the second term, we havewhere the second inequality holds because of the Markov inequality; the third inequality holds because $\log N\ge 2$ and $\frac{N^{0.5-\alpha }}{\log N}\ge 2k;$ and the last inequality holds because of Lemma 6.

Finally, we remark that we choose $k=\mathrm{\Omega }(b)$ to prove Lemma F.2, which is the technical reason we assumed b is finite. We, however, believe our results hold even for $b=\infty .$

8. Simulations

In this section, we confirm our theoretical results of ISSP in Theorem 1 and the large system insensitivity in Theorem 2 with simulations. We considered two policies, JSQ and JIQ, and $\lambda =1-N^{\alpha }.$ We conjecture that our results hold even for $\alpha =0.5$ because it holds for any $\alpha <0.5.$ To confirm this, we used $\alpha =0.5$ in our experiments (the additional experiments with $\alpha =0.3$ and $\alpha =0.4$ can be found in Appendix G). We set the buffer size b = 10, and in our simulations, no blocking events were observed.

8.1. Large System Insensitivity Under JSQ and JIQ

We first studied the average total queue length (per server) $\mathbb{E}[{\sum }_{i=1}^b{S}_i]$ and the waiting probability $\mathbb{P}(\mathcal{W})$ under a Coxian-4 service time distribution with the parameters $p=[0.5,0.5,0.5,1.0]$ and $\mu =[1.875,1.875,1.875,1.875].$ We plotted $\mathbb{E}[{\sum }_{i=1}^b{S}_i]$ and $\mathbb{P}(\mathcal{W})$ versus the number of servers N. The results are obtained with 10 trials, and each trial has 10⁷ steps. From Figure 7, the waiting probability tends to zero when N increases as we expected and JIQ almost has the identical performance as JSQ. We then investigated $\mathbb{E}[{\sum }_{i=1}^b{S}_i]$ and $\mathbb{P}(\mathcal{W})$ versus Coxian-M with various number of phases M and fixed $N={10}^4.$ In particular, we consider Coxian-M with $p=[p_1,p_2,\dots ,1]$ and $\mu =[{\mu }_1,{\mu }_2,\dots ,{\mu }_M],$ where $p_m=0.5,\forall 1\le m\le M-1$ and $\overline{\mu }={\mu }_m,\forall m$ (identical service times). We plotted $\mathbb{E}[{\sum }_{i=1}^b{S}_i]$ and $\mathbb{P}(\mathcal{W})$ versus M. From Figure 8, we can observe that the average queue length and the waiting probability remain roughly the same under different Ms, which confirms the insensitivity.

8.2. ISSP Under JSQ and JIQ

In this section, we investigated ISSP by studying the trajectory of $S_{1,m}(t),\forall m$ under JSQ and JIQ. We considered n = 10,000 and a Coxian-4 service time distribution with $p=[0.5,0.5,0.5,1.0]$ and $\mu =[1.875,1.875,1.875,1.875].$ We plotted $S_{1,m}(t)$ of JSQ and JIQ in Figure 9. The results are obtained with 10 trials, and each trial has 10⁶ steps. We observed that $S_{1,m(t)}$ under both policies concentrates around dash lines $s_{1,m}^{*},\forall m,$ which confirms the high probability bounds in one. The system was initialized with the zero-waiting equilibrium instead of the empty state in our simulations.

We also provide additional experiments in Appendix G. These additional results justify the finite buffer assumption and demonstrate the performance of LB-zero under different traffic intensities and the service time distributions with different variations.

9. Conclusions and Discussions

In this paper, we studied a distributed queueing system under Coxian service time distributions in the sub-Halfin-Whitt regime. We established that a set of load balancing policies, named LB-zero, achieves asymptotic zero-waiting, that is, insensitive in the large-system regime. To tackle the nonmonotonicity under general service time distributions, we developed a technique, called iterative state-space peeling (ISSP), which iteratively removes the low-probability states, and results in a small state-space that can be analyzed using a simple mean-field model. This ISSP approach may be used for other problems as well. One possible application is to study load-balancing in many server systems with heterogeneous servers or jobs belonging to multiple priority classes. For heterogeneous servers, we can use ISSP to identify the “typical” load of each type of the servers; and for jobs with different priorities, we can use ISSP to identify the “typical” distribution of job types in the system. In the reduced state space based on the typical load of the typical distribution at the steady state, the steady-state performance of the many server system may become tractable like in this paper.