Exploiting Data Locality to Improve Performance of Heterogeneous Server Clusters

Motivated by this, in the current paper, we consider a bipartite graph model for large-scale load balancing systems, which has recently gained popularity in the research community. In this model, a bipartite graph between the servers and task types describes the compatibility between the two, where an edge represents the server’s ability to process the corresponding task type. This encompasses the classical full-flexibility models as those having a complete bipartite compatibility graph. An immediate difficulty of the new model is that when the graph is non-trivial (i.e., not a collection of isolated pairs or a complete bipartite graph), the mean-field techniques break down. This is because, the queues no longer remain exchangeable, making the aggregate processes such as the vector of number of servers with queue length $i$ with $i=0,1,2,\ldots$ non-Markovian. In addition, we also consider that each dispatcher handles the arrival flow of one of $K$ possible task types and that there are $M$ server types. The rate of service at a server depends on its type. Throughout the paper, the key quantity of interest will be the global occupancy process $\mathbf{q}^{N}(t)=(q^{N}_{m,l}(t),m=1,\ldots,M,l\geq 1)$, where $q^{N}_{m,l}(t)$ represents the fraction of servers of type $m$ with queue-length at least $l$ at time $t$ in the $N$-th system with $N$ servers, and we will look at the large-system asymptotic regime: $N\to\infty$.

Due to the compatibility constraints, the servers become non-exchangeable, even if they belong to the same type. This causes most of the existing frameworks [8, 17, 24] to break down. To characterize the process-level limit of the queue length process, we take resort to the theory of weakly interacting particle systems and asymptotically couple the evolution of the $N$-dimensional vector of queue lengths with an appropriately defined infinite system of independent McKean-Vlasov processes [26, 15]. We also show the asymptotic independence of any finite number of queue length processes, also know as the propagation of chaos property. This convergence of the queue length processes (in $L_{2}$ sense) is then used to establish the transient convergence of the occupancy process. One downside of the above convergence is that it depends on the assumption that the initial queue lengths within each set of servers of the same type are independent and identically distributed (i.i.d.) and are independent across the the set of servers of different types. Due to this assumption, this convergence result cannot be used to establish the interchange of $t\to\infty$ and $N\to\infty$ limits, which is crucial in studying the limit of steady states.

To overcome this issue, we use the framework of [22], recently introduced in the context of homogeneous systems. Here, a notion called proportional sparsity for graph sequences was introduced, which ensures that the empirical queue length distribution within the set of compatible servers of any dispatcher is close to the empirical queue length distribution of the entire system. This was used in [22] to construct conditions on graphs that match the performance of a fully flexible system. In the current setup, however, this notion is inadequate, since our goal is not to match the performance of the fully flexible system (which is usually poor under heterogeneity). That is why, we extend this notion to what we call the clustered proportional sparsity for a sequence of graphs with increasing size, to accommodate the heterogeneous systems. The clustered proportional sparsity property allows us to construct a stochastic coupling between the system and another intermediate system whose task allocation is done by a carefully constructed algorithm called GWSQ($d$) (Algorithm 1). This coupling with the intermediate system, along with clustered proportional sparsity, helps us establish that if the initial occupancy of two systems are close, then the distance (in the $\ell_{1}$-norm) between their global occupancy remains small uniformly over any finite time interval. In turn, it implies that their limits of the global occupancy systems are the same. As a consequence, we can remove the i.i.d. assumption of the initial queue lengths, since the above guarantees that under clustered proportional sparsity, the convergence of the occupancy process depends only on the initial occupancy and not on how the individual queues are distributed.

The above process-level limit result shows that the transient limit of the occupancy process can be described as a system of ODEs that depend on various graph parameters. Next, we also show that the interchange of limits holds and that the sequence of occupancy states in stationarity, converges weakly to the unique fixed point of the ODE. One celebrated feature of the classical JSQ($d$) policy for homogeneous systems under full flexibility is that the steady-state queue length decays doubly exponentially as $\lambda^{(d^{i}-1)/(d-1)}$, where $\lambda\in(0,1)$ is the load per server [17, 32]. We establish this double-exponential decay property for the heterogeneous system.

It is worthwhile to note that the strength of the above results lie in that they hold for arbitrary deterministic sequence of graphs satisfying certain properties. However, we show that all these properties are satisfied almost surely by a sequence of inhomogeneous random graphs with parameters prescribed by the theorems. This makes it easy to design graphs with the desired favorable properties.

The research on task allocation systems with limited flexibility can be traced back to the works of Turner [30] and Foss and Chernova [9]. Of particular importance to the current work, Foss and Chernova [9] considered stability properties of the system using the fluid model. Later, Bramson [4] generalized some parts of results in [9] to a broad class of Join-Shortest-Queue (JSQ)-type systems, including the JSQ($d$) policy, via the Lyapunov function approach. Stolyar [23] considered optimal routing in output-queued flexible server system, which is essentially the bipartite graph model for the load balancing system. Here the author considered a system with a fixed number of servers and dispatchers in the conventional heavy traffic regime and proposed a routing policy that is optimal in terms of server workload. Recently, Cruise et al. [7] considered load balancing problems on hypergraphs and proved its stability conditions. The above works, however, did not aim to precisely characterize the system performance in the large-scale scenario.

The analysis in the large-scale scenario became prominent in the last decade, with the emergence of its applications to load balancing in data centers and cloud networks. In the full-flexibility setup, the analysis of heterogeneous-server systems gained some attention. In this case, Stolyar [24, 25] studied the zero-queueing property of the Join-Idle-Queue (JIQ) policy, Mukhopadhyay et al. [21] and Mukhopadhyay and Mazumdar [20] analyzed the JSQ($d$) policy in heterogeneous systems with processor-sharing service discipline, Hurtado-Lange and Maguluri [13] studied the throughput and delay optimality properties of JSQ($d$), and Bhambay and Mukhopadhyay [3] studied a speed-aware JSQ policy. The above works on the JSQ($d$) policy observe that the stability region shrinks if the dispatcher applies the JSQ($d$) policy blindly. One way to mitigate this performance degradation is to take the server speeds into consideration while sampling servers or while assigning tasks to the sampled servers. Such a ‘hybrid JSQ($d$)’ scheme is able to recover the stability region. The current work can be contrasted with this approach. First, in the presence of data locality, both the server speeds and the underlying compatibility constraints need to be taken into account during the sampling procedure, and the approach becomes significantly more complicated. Second, we show how exploiting the data locality, the blind JSQ($d$) policy can recover the stability region and even achieve the double-exponential decay of tail probabilities of the stead-state queue length distribution. One advantage of the latter approach is that the dispatchers can be oblivious to the server speeds, which reduces the implementation complexity and also, makes it robust against changes to the servers (e.g., when servers are add/removed).

Recently, Allmeier and Gast [2] studied the application of (refined) mean-field approximations for heterogeneous systems. Their method is using an ODE to approximate the evolution of each server, and the error vanishes as the system scales. However, this method cannot be directly used in our case. Due to the bipartite compatibility graph structure, it is hard to capture the interactions between two servers, which means that we cannot write the transition rates of the underlying Markov chain as [2] does. Also, one important assumption in their work is the finite buffer, but we consider the infinite buffer case here.

The aspect of task-server compatibility constraints in large-scale load balancing and scheduling gained popularity only recently, as the data locality became prominent in data centers and cloud networks. This led to many works in this area [11, 18, 5, 22, 33, 29, 28]. All these works consider homogeneous processing speeds at the servers. The initial works [30, 11] focused on certain fixed-degree graphs and showed that the flexibility to forward tasks to even a few neighbors with possibly shorter queues may significantly improve the waiting time performance as compared to dedicated arrival streams or a collection of independent M/M/1 queues. Tsitsiklis and Xu [28, 29] considered asymptotic optimality properties of the bipartite graph topology in an input-queued, dynamic scheduling framework. Later, in the (output-queued) load balancing setup, Mukherjee et al. [18] considered the JSQ policy and Budhiraja et al. [5] considered the transient analysis of the JSQ($d$) policy on non-bipartite graphs. The goal in these papers was to provide sufficient condition on the graph sequence to asymptotically match the performance of a complete graph. Here we should mention that the non-bipartite graph model cannot be used to capture the data locality constraints. In the presence of data locality constraints, the analysis of the JSQ($d$) policy for homogeneous systems, including both transient and interchange of limits, was performed by Rutten and Mukherjee [22]. Weng et al. [33] is the first to consider the large-scale heterogeneous-server model under data locality. They showed that the Join-the-Fastest-Shortest-Queue (JFSQ) and Join-the-Fastest-Idle-Queue (JFIQ) policies achieve asymptotic optimality for minimizing mean steady-state waiting time when the bipartite graph is sufficiently well connected. However, these results fall in the category of JSQ-type policies where the asymptotic behavior is degenerate in the sense that the queue lengths at servers can be either 0 or 1. Naturally, the results and their analysis are very different from the JSQ($d$)-type policies where queues of any length is possible.

Let $\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$. For a set $S$, its cardinality is denoted as $|S|$. For a polish space $\mathcal{S}$, the space of right continuous functions with left limits from $[0,\infty)$ to $\mathcal{S}$ is denoted as $\mathbb{D}([0,\infty),\mathcal{S})$, endowed with the Skorokhod topology. The distribution of $\mathcal{S}$-valued random variable $X$ will be denoted as $\mathcal{L}(X)$. For a function $f:\ [0,\infty)\rightarrow\mathbb{R}$, let $\left\lVert f\right\rVert_{*,t}\coloneqq\sup_{0\leq s\leq t}|f(s)|$. The distribution of $\mathcal{S}$-valued random variable $X$ will be denoted as $\mathcal{L}(X)$. For $x\in\mathcal{S}$, the Dirac measure at the point $x$ is denoted as $\delta_{x}$. $\left\lVert\cdot\right\rVert_{p}$ represents the $\ell_{p}$-norm. Define ${X\choose Y}=\frac{X(X-1)\cdots(X-Y+1)}{Y!}$ if $X\geq Y$ and is 0, otherwise. RHS is the acronym of Right Hand side.

As discussed earlier, when the server speeds are heterogeneous, the fully flexible systems (with the complete bipartite compatibility graph) may not be stable under the JSQ($d$) policy, even if we assume that the sufficient service capacity in (2.1) is satisfied. The next lemma provides a necessary and sufficient condition for ergodicity of the queue length process. Recall $\delta^{N}_{i}=|\mathcal{N}^{N}_{w}(i)|$. For any fixed $N$, define

$$\rho^{N}\coloneqq\max_{\begin{subarray}{c}U\subseteq V^{N}\\ U\neq\emptyset\end{subarray}}\Big{\{}\Big{(}\sum_{j\in U}\sum_{m\in\mathcal{M}}\mathds{1}_{(j\in V^{N}_{m})}u_{m}\Big{)}^{-1}\sum_{i\in W^{N}}\Big{(}\mathds{1}_{(\delta^{N}_{i}\geq d)}\sum_{\begin{subarray}{c}S\subseteq(U\cap\mathcal{N}^{N}_{w}(i)):\\ |S|=d\end{subarray}}\frac{\lambda}{{\delta^{N}_{i}\choose d}}+\mathds{1}_{(\delta^{N}_{i}<d)}\frac{|U\cap\mathcal{N}^{N}_{w}(i)|}{\delta^{N}_{i}}\Big{)}\Big{\}}.$$

(3.1)

The above lemma is an immediate consequence of [9, Theorem 2.5]; see also [4]. We omit its proof. Intuitively, $\rho^{N}<1$ means that in the $N$-th system, for any subset $U$ of servers with possibly long queues (compared to the rest servers), the total rate at which tasks are assigned to some server in this set must be less than the rate of departure from this set.

Since we are interested in large-$N$ behavior, we will assume certain asymptotic version of the above stability criterion. This is fairly standard in the large-system analysis, as one would want to avoid the ‘heavy-traffic’ regime when $\rho^{N}\uparrow 1$ as $N\to\infty$. The behavior in the latter scenario is typically qualitatively different from the so-called ‘subcritical’ regime as defined below.

Throughout this paper we will assume that the sequence of systems under consideration is in subcritical regime. From Lemma 3.1, it is immediate that if a sequence of systems is in subcritical regime, then its queue length process is ergodic for all large enough $N$. The potential non-ergodicity of fully flexible, heterogeneous server clusters brings us to the question that when the sufficient service capacity in (2.1) is satisfied, whether we can design the underlying compatibility structure carefully so that the queue length process is ergodic. In other words, can we regain the stability region? Proposition 3.3 below shows that this is indeed the case. In some sense, this highlights the first-order improvements (i.e., in terms of stability properties) of a careful compatibility structure design in contrast to a fully flexible system.

In the following sections, we will demonstrate, in addition to the first-order improvements, how asymptotic queue length distribution can be improved as well, for example, in terms of having a double-exponential decay of tail probabilities.

The proof of Proposition 3.3 is provided in Appendix A. It relies on first building a simple criteria involving the system parameters, which, for sequence of systems satisfying Condition 2.1, ensures stability for all large-enough $N$ (Lemma 3.4). Then we show that given other parameters, a value of $(p_{k,m})_{k\in\mathcal{K},m\in\mathcal{M}}$ satisfying this criteria can be found by checking the feasibility region defined by $M$ inequalities.

We end this subsection by presenting the above-mentioned simple asymptotic criteria for subcriticality. The proof is given in Appendix A. Denote $\delta_{k}\coloneqq\sum_{m\in\mathcal{M}}p_{k,m}v_{m}$ for each $k\in\mathcal{K}$.

Our first main result characterizes the process-level limit of the queue-length process $\big{(}X^{N}_{j}$, $j\in V\big{)}$, as $N\to\infty$, when the starting states $\big{\{}X^{N}_{j}(0):j\in V^{N}_{m}\big{\}}$ are i.i.d. for all $m\in\mathcal{M}$ and independent across different $m$-values. When the sequence of graphs $\{G^{N}\}_{N}$ satisfies a stronger condition, called clustered proportional sparsity (Definition 3.8), the i.i.d. condition can be removed. This is the content of Section 3.3.

Now, note that for a fixed $N\geq 1$, $\{X^{N}_{j}:j\in V_{N}\}$ is a system of $N$ stochastic processes with mean-field type interactions. Exploiting tools from the theory of weakly interacting particles, we show in Theorem 3.5 below that as the system size becomes large, queue-length processes converge weakly to those of an infinite system of independent McKean-Vlasov processes $\{X_{j}:j\in\mathbb{N}\}$ (see e.g. [26, 15]). In fact, using a suitable coupling to be described in more details in Section 4.1, the convergence holds in $L_{2}$. For ease of describing such processes and coupling, although we only assumed that certain fractions of servers are of certain task types in the model description, it will be convenient to fix the type of each server $j\in\mathbb{N}$ in this subsection, by defining a membership map $\mathbf{M}:\mathbb{N}\rightarrow\mathcal{M}$, so that $V^{N}_{m}=\{j\in V^{N}:\mathbf{M}(j)=m\}$ with $\lim_{N\rightarrow\infty}\frac{|V^{N}_{m}|}{N}=v_{m}$ and $V_{m}=\lim_{N\rightarrow\infty}V^{N}_{m}$ for each $m\in\mathcal{M}$. With such fixed server types and $X_{j}^{N}(0)\equiv X_{j}(0)$, let

	$\displaystyle X_{j}(t)$	$\displaystyle=X_{j}(0)-\int_{0}^{t}\mathds{1}_{\big{(}X_{j}(s-)>0\big{)}}D_{j}(ds)+\int_{[0,t]\times\mathbb{R}_{+}}\mathds{1}_{\big{(}0\leq y\leq C_{j}(s-)\big{)}}A_{j}(dsdy),$		(3.3)
	$\displaystyle C_{j}(t)$	$\displaystyle=d\xi\sum_{k\in\mathcal{K}}\frac{p_{k,m}w_{k}}{\delta_{k}}\sum_{(M_{2},...,M_{d})\in\mathcal{M}^{d-1}}h_{t}(j,M_{2},...,M_{d}),$		(3.4)

where $\mathbf{M}(j)=m$ and

	$\displaystyle h_{t}(j,M_{2},...,M_{d})$	$\displaystyle=\prod_{h=2}^{d}\frac{v_{M_{h}}p_{k,M_{h}}}{\delta_{k}}\int_{\mathbb{N}^{d-1}}b\big{(}X_{j}(t),x_{j_{2}},...,x_{j_{d}}\big{)}\mu^{M_{2}}_{t}(dx_{j_{2}})\cdots\mu^{M_{d}}_{t}(dx_{j_{d}}),$
	$\displaystyle b(\mathbf{x})$	$\displaystyle=b(x_{1},...,x_{d})\coloneqq\sum_{r=1}^{d}\frac{1}{r}\mathds{1}_{\big{(}x_{1}=\min_{j\in[d]}\mathbf{x},|\operatorname*{arg\,min}\mathbf{x}|=r\big{)}},\quad\mathbf{x}=(x_{1},...,x_{d})\in\mathbb{N}_{0}^{d},$		(3.5)
	$\displaystyle\mu^{m}_{t}$	$\displaystyle=\mathcal{L}\big{(}X_{i}(t)\big{)},\quad\forall\,i\in V_{m},m\in\mathcal{M},t\geq 0.$

Here, $\{D_{j}:j\in V_{m}\}$ are i.i.d. Poisson processes with rate $u_{m}$ for each $m\in\mathcal{M}$, $\{A_{j}:j\in\mathbb{N}\}$ are i.i.d. Poisson random measures on $[0,\infty)\times\mathbb{R}_{+}$ with intensity $\lambda dsdy$, and all $D_{j}$’s and $A_{j}$’s are independent. Loosely speaking, $A_{j}$ corresponds to the arrival processes and $D_{j}$ corresponds to the departure processes at servers. We note that the existence and uniqueness of solutions to (3.3) and (3.4) can be proved by standard arguments (see e.g., [26, 15]) using the boundedness and Lipschitz property of the functions $b$ and $x\mapsto\mathds{1}_{(x>0)}$ on $\mathbb{N}_{0}$.

Theorem 3.5 gives us the limit law of all individual queues. Next, in Theorem 3.6, we will show how such a server-level convergence can be used to obtain a convergence result for the global occupancy process $\mathbf{q}^{N}(\cdot)$ to a deterministic dynamical system, which was our primary goal. The proofs of Theorem 3.5 and Theorem 3.6 are provided in Section 4.

Theorem 3.6 requires the strong assumption that for each $m\in\mathcal{M}$, $X^{N}_{j}(0)$, $j\in V^{N}_{m}$, are i.i.d.. In order to argue the interchange of limits, we need to relax this assumption on initial states. This is because the arguments for the interchange of limits involves initiating the prelimit system at the steady state and then showing that as $N\to\infty$, the system must converge to the unique fixed point of the limiting ODE. The above requires us to characterize the (process-level) limiting trajectory of the system starting from arbitrary occupancy state. We achieve this in this section.

Intuitively, the assumption of i.i.d. in Theorems 3.5 and 3.6 ensures that the local occupancy observed by any dispatcher $i\in W^{N}_{k}$, $k\in\mathcal{K}$ is ‘close’, in suitable sense, to the average occupancy at the entire system. This phenomenon can be ensured asymptotically, even without the i.i.d. assumption if the graph sequence satisfies a property we call the clustered proportional sparsity. This notion was first introduced for the homogeneous systems in [22]. The definition below is a modified notion that is suitable for the current heterogeneous setting.

The proof of Theorem 3.10 is given in Section 4.4

In the last section, we showed the process-level convergence of global occupancy process $\mathbf{q}^{N}(\cdot)$ to a mean-field limit $\mathbf{q}(\cdot)$. In this section, we will establish the convergence of the sequence of stationary distributions to the unique fixed point of the mean-field limit by establishing the interchange of large-$N$ and large-$t$ limits: $\lim_{t\rightarrow\infty}\lim_{N\rightarrow\infty}\mathbf{q}^{N}(t)=\lim_{N\rightarrow\infty}\lim_{t\rightarrow\infty}\mathbf{q}^{N}(t)$. Throughout this section, we will assume that the sequence of systems is in subcritical regime (recall Definition 3.2). The first result below states that the limiting system of ODEs have a unique fixed point $\mathbf{q}^{*}$ and it satisfies the global stability property, i.e., for any initial point $\mathbf{q}(0)\in\mathcal{S}$, $\lim_{t\rightarrow\infty}\mathbf{q}(t)=\mathbf{q}^{*}$.

The proof of Theorem 3.11 is given in Section 5. It relies on a monotonicity property of the system, which ensures that for two processes $\mathbf{q}^{1}(\cdot)$ and $\mathbf{q}^{2}(\cdot)$, if $\mathbf{q}^{1}(0)\leq\mathbf{q}^{2}(0)$, then $\mathbf{q}^{1}(t)\leq\mathbf{q}^{2}(t)$ for all $t\geq 0$ (see ref. [24, 14]).

The last ingredient that we need in order to prove the interchange of limits is to establish tightness of the sequence of random variables $\{\mathbf{q}^{N}(\infty)\}_{N\geq 1}$ under suitable metric, where $\mathbf{q}^{N}(\infty):=\lim_{t\to\infty}\mathbf{q}^{N}(t)$. Here, as before, we should note that the process $(\mathbf{q}^{N}(t))_{t\geq 0}$ is not Markovian. That is why, the random variable $\mathbf{q}^{N}(\infty)$ should be interpreted as the functional applied to the steady-state system. The tightness result is stated in the next theorem.

Theorem 3.12 is proved in Section 5. The key idea is to use Lyapunov function approach to bound the expected sum of tails $q^{N}_{m,l}(\infty)$. Combining Theorems 3.10, 3.11, and 3.12 we can prove the following interchange of limits result.

One major discovery about the JSQ($d$) policy for the classical, homogeneous, fully flexible system is that the limit of the stationary distribution (which, in our case, is given by $\mathbf{q}^{*}$) has a double-exponential decay of tail [17, 32] for any $d\geq 2$. This is in sharp contrast with the (single) exponential decay of the corresponding tail for random routing or $d=1$. In fact, in this case, for any $d\geq 2$, $\mathbf{q}^{*}$ can be characterized explicitly as: $q^{*}_{l}=\lambda^{\frac{d^{l}-1}{d-1}}$, where $q_{l}^{*}$ is the (limiting) steady-state fraction of servers with queue length at least $l=1,2,\ldots$. In the current case of heterogeneous systems, it is intractable to characterize the fixed point $\mathbf{q}^{*}$ explicitly. However, as stated in the next theorem, we can still prove that the doubly exponential decay of the tails $q^{*}_{m,l}$ for each $m\in\mathcal{M}$ holds.

Sections 3.1–3.4 characterize the performance of the occupancy process for arbitrary deterministic sequence of systems where the underlying graph sequence satisfies certain properties. In particular, Condition 2.1 and Definition 3.8 provide a sufficient criteria under which both the process-level convergence (Theorem 3.10) and interchange of limits (Theorem 3.13) hold. In this section, we show that graphs satisfying the above required criteria can be obtained easily if the compatibility graph is designed suitably randomly. Given the asymptotic edge-density parameters in Condition 2.1, we define a certain sequence of inhomogeneous random graphs or irg as follows.

For any $\mathbf{p}$ for which the asymptotic stability criterion holds, we have the following result for the sequence of irg($\mathbf{p}$).

The proof of Theorem 3.16 is provided in Appendix I. It relies on verifying that the sequence of irg($\mathbf{p}$) graphs satisfies Condition 2.1 and the property of clustered proportional sparsity almost surely. The verification involves using concentration of measure arguments to establish structural properties of the compatibility graphs.

First, we will need a characterization of the evolution of the queue length process at each server. To describe this evolution, let us introduce the following notations:

	$\displaystyle\texttt{set}^{N}(j)$	$\displaystyle\coloneqq\Big{\{}(j_{2},...,j_{d})\in[N]^{d-1}:(j,j_{2},\ldots,j_{d})\text{ are distinct}\Big{\}},$		(4.1)
	$\displaystyle\texttt{sett}^{N}(j)$	$\displaystyle\coloneqq\Big{\{}(j_{2},\ldots,j_{d},j^{\prime}_{2},\ldots,j^{\prime}_{d})\in[N]^{2d-2}:(j_{2},...,j_{d})\in\texttt{set}^{N}{(j)},(j^{\prime}_{2},...,j^{\prime}_{d})\in\texttt{set}^{N}{(j)},$
		$\displaystyle\hskip 199.16928pt(j_{2},...,j_{d})\cap(j^{\prime}_{2},...,j^{\prime}_{d})\neq\emptyset\Big{\}}.$		(4.2)

To represent the graph, define the edge occupancy $\xi^{N}_{i,j}$ to be the binary variable:

$$\xi^{N}_{i,j}=\begin{cases}1,&\text{if}\ (i,j)\in E^{N},\\ 0,&\text{otherwise},\end{cases}\quad\mbox{for all }i\in W^{N},j\in V^{N}.$$

Recall the function $b$, Poisson processes $\{D_{j}\}$ and Poisson random measures $\{A_{j}\}$ in and below (3.5). By Condition 2.1, for all large enough $N$, all dispatchers in the $N$-th system have at least $d$ neighbors. Hence, WLOG, in the rest of this section, we will only consider the case $\delta^{N}_{i}\geq d$, $\forall i\in W^{N}$. In that case, due to the Poisson thinning property, note that we can write $X^{N}_{j}(t)$ as follows:

$$X^{N}_{j}(t)=X^{N}_{j}(0)-\int_{0}^{t}\mathds{1}_{\big{(}X^{N}_{j}(s-)>0\big{)}}D_{j}(ds)+\int_{[0,\infty)\times\mathbb{R}_{+}}\mathds{1}_{\big{(}0\leq y\leq C^{N}_{j}(s-)\big{)}}A_{j}(dsdy),$$

(4.3)

where

	$\displaystyle C^{N}_{j}(s)$	$\displaystyle=\sum_{i\in W^{N}}\xi^{N}_{i,j}\sum_{(j_{2},...,j_{d})\in\texttt{set}^{N}(j)}\frac{\xi^{N}_{i,j_{2}}\times\cdots\times\xi^{N}_{i,j_{d}}}{{\delta^{N}_{i}\choose d}(d-1)!}b\big{(}X^{N}_{j}(s),X^{N}_{j_{2}}(s),...,X^{N}_{j_{d}}(s)\big{)}$		(4.4)
		$\displaystyle=\sum_{k\in\mathcal{K}}\sum_{i\in W^{N}_{k}}\xi^{N}_{i,j}\sum_{(j_{2},...,j_{d})\in\texttt{set}^{N}(j)}\frac{\xi^{N}_{i,j_{2}}\times\cdots\times\xi^{N}_{i,j_{d}}}{{\delta^{N}_{i}\choose d}(d-1)!}b\big{(}X^{N}_{j}(s),X^{N}_{j_{2}}(s),...,X^{N}_{j_{d}}(s)\big{)}.$

The RHS of the first summation in (4.4) represents the probability that a job arriving at the dispatcher $i\in W^{N}$ will be assigned to the server $j\in V^{N}$ given the state $\big{(}X^{N}_{j},j\in V^{N}\big{)}$. Moreover, by Condition 2.1, the term $C^{N}_{j}$ for all $j\in V^{N}$ can be upper bounded, uniformly for all $t$, by a constant for all large enough $N$, which is stated in Lemma 4.2 below.

When we do some estimation, like bounding the term $C^{N}_{j}$, we need to uniformly bound the number of the neighbors of servers or dispatchers. Such uniformity is stated in Lemma 4.1 and is a direct result of Condition 2.1. Recall $\delta^{N}_{i}=|\mathcal{N}^{N}_{w}(i)|$ and $\delta_{k}=\sum_{m\in\mathcal{M}}p_{k,m}v_{m}$.