Utility Maximizing Load Balancing Policies

1. Introduction

We consider a service system where incoming tasks are instantaneously assigned to one out of many heterogeneous server pools. All the tasks sharing a server pool are executed in parallel and the execution times do not depend on the class of the server pool or the number of tasks currently contending for service. Nevertheless, associated with each server pool is a not necessarily increasing concave utility function that does depend on the class of the server pool and the number of tasks currently sharing it. These features are characteristic of streaming and online gaming services, where the duration of tasks is mainly determined by the application, but still congestion can have a strong impact on the experienced performance (e.g., video resolution and fluency).

The goal is to assign tasks to maximize the overall utility of the system, defined as the aggregate utility of all the server pools normalized by the number of server pools. We derive an upper bound for its stationary mean through an optimization problem where the optimization variable is a sequence that describes the distribution of a fractional number of tasks across the server pools; the objective of the problem is the overall utility function, and the main constraint is that the total number of tasks must be equal to the offered load of the system. We construct an optimal (fractional) task assignment that solves this problem and has a particularly insightful structure, and we formulate the upper bound for the mean stationary overall utility in terms of this solution.

Armed with the previous insight, we propose and analyze two assignment policies that maintain the occupancy state of the system aligned with an optimal task assignment. Specifically, we examine a policy that assigns every new task to a server pool with the largest marginal utility; this policy is dubbed join the largest marginal utility (JLMU). We also introduce a multithreshold policy that follows the same greedy principle but only approximately and uses significantly less state information. The optimal threshold values depend on the typically unknown offered load of the system and are adjusted over time through an inbuilt learning scheme; thus, we name this policy self-learning threshold assignment (SLTA). Assuming exponential service times, we characterize the asymptotic transient and stationary behavior of both policies on the scale of the number of server pools, and we prove that both policies achieve the upper bound for the mean stationary overall utility as the number of server pools grows large.

A fundamental difference between JLMU and SLTA is that the former is naturally agnostic to the offered load, whereas for the latter, the optimal thresholds depend on the offered load. However, we show that the online learning scheme of SLTA is capable of finding the optimal threshold values without any prior knowledge of the offered load, which makes it possible to deploy SLTA if the offered load is not known in advance.

1.1. Main Contributions

The main contribution of this paper is an upper bound for the mean stationary overall utility that is asymptotically tight for exponentially distributed service times and thereby serves as a crucial performance benchmark. The asymptotic tightness of the upper bound is proved by studying the stationary behavior of JLMU and SLTA in the regime where the number of server pools grows large and by establishing that both assignment policies achieve the upper bound in the latter regime.

The analysis of JLMU is based on a fluid limit given by an infinite system of differential equations with a discontinuous right-hand side. We prove that the associated initial value problem always has a unique solution by making a connection with a system of integral equations, expressed in terms of Skorokhod one-dimensional reflection mappings and using a uniqueness result for certain Kolmogorov backward equations. Moreover, we show that the fluid limit holds with respect to an ${\ell }_1$ norm and that the system of differential equations is globally asymptotically stable with respect to this norm. These results are used to prove that the stationary distribution of the process that describes the occupancy state of the system converges in ${\ell }_1$ to an optimal task assignment for the offered load of the system. The asymptotic optimality of JLMU is then established by proving that the stationary overall utilities form a convergent and uniformly integrable sequence of random variables; the proof of the latter properties exploits a representation of the overall utility as a linear functional on ${\ell }_1$ and our convergence results with respect to the ${\ell }_1$ norm.

Although SLTA is simple to implement, its analysis is inherently challenging due to the complex interdependence between two components of the policy. Namely, the dispatching rule, which depends on the multiple thresholds, and the online learning scheme, which adjusts the thresholds over time. Furthermore, an additional technical difficulty is that the learning scheme is triggered by excursions of the occupancy state of the system that asymptotically vanish on the scale of the number of server pools.

To analyze the large-scale transient behavior of SLTA, we use a methodology of Goldsztajn et al. (2021a) that allows overcoming of the aforementioned challenges by means of a nontraditional fluid limit analysis. In this paper, we extend the latter methodology to also prove weak convergence of the stationary distribution of the occupancy process and thresholds. Here our contributions are proofs of ergodicity and tightness of stationary distributions through a careful drift analysis, as well as a suitably adapted interchange of limits argument designed to leverage the large-scale transient result obtained with the methodology of Goldsztajn et al. (2021a). Equipped with the convergence results for the stationary distributions, we prove the asymptotic optimality of SLTA in a similar way as for JLMU by showing that all our limit theorems hold with respect to the ${\ell }_1$ norm and exploiting the linear representation of the overall utility function.

1.2. Related Work

Load balancing and task assignment in parallel-server systems has received immense attention in the past decades; some relevant papers are Winston (1977), Vvedenskaya et al. (1996), Mitzenmacher (2001), Lu et al. (2011), Stolyar (2015), and Eschenfeldt and Gamarnik (2018). Traditionally the focus used to be on performance, but more recently the implementation overhead has emerged as an equally important issue. In large-scale deployments, this overhead has two main sources: the communication burden of messaging between the dispatcher and the servers and the operational cost of storing and managing state information at the dispatcher (Gamarnik et al. 2018, 2020). We refer to Van der Boor et al. (2022) for an extensive survey on scalable load balancing.

Whereas the load balancing literature has been predominantly concerned with systems of parallel single-server queues, the present paper considers an infinite-server setting where the service times of tasks do not depend on the number of competing tasks. As noted earlier, this feature is characteristic of streaming applications, where the level of congestion does not significantly affect the duration of tasks. The level of congestion has, however, a strong impact on the amount of resources received by individual streaming sessions, and thus on the experienced quality-of-service, which can be modeled through utility functions. Infinite-server dynamics have been commonly adopted as a natural paradigm for modeling streaming sessions on flow-level (Benameur et al. 2002, Key et al. 2004), and the problem of managing large data centers serving streaming sessions has been recently addressed in Mukherjee et al. (2020). Systems with infinite-server dynamics have also been analyzed in Mukhopadhyay et al. (2015a, b), Karthik et al. (2017), and Xie et al. (2015), which concern loss models that are different in nature from the setting considered in the present paper.

When the server pools are homogeneous, the overall utility is a Schur-concave function of the vector describing the number of tasks at each server pool. In this case, maximizing the aggregate utility of the system boils down to equalizing the number of tasks across the various server pools. Join the shortest queue (JSQ) maximizes the mean stationary overall utility of the system for exponential service times and, in fact, has stronger stochastic optimality properties (Menich and Serfozo 1991, Sparaggis et al. 1993). In the homogeneous setting, JLMU reduces to JSQ and is thus optimal for exponential service times. In addition, SLTA reduces to the policy considered in Goldsztajn et al. (2021a, b), which asymptotically matches the performance of JSQ on the fluid and diffusion scales for exponentially distributed service times. Although the policy considered in Zhou et al. (2017, 2018) and Horváth et al. (2019) is similar to SLTA in name, this policy does not equalize the queue lengths.

The problem of maximizing the overall utility of the system is more challenging if the server pools are heterogeneous as in this paper. Heterogeneity is the norm in data centers, where servers from different generations coexist because old machines are only gradually replaced by more powerful versions; as shown in Figure 1, this feature has been recently addressed in the load balancing literature for single-server models (Gardner and Stephens 2019, Jaleel et al. 2020, Gardner et al. 2021) but not in the infinite-server context. When the server pools are heterogeneous, it is no longer optimal to maintain an evenly balanced distribution of the load; in fact, it is not even obvious at all how tasks should be distributed to maximize the overall utility function, and the optimal distribution of tasks across the server pools depends on this function. Another striking difference with the homogeneous setting is that JLMU is generally not optimal in the prelimit for exponentially distributed service times; we establish that, in general, the optimality is only achieved asymptotically in the heterogeneous case.

From a theoretical perspective, one of the most interesting features of SLTA is its capacity to learn the offered load of the system. The problem of adaptation to unknown demands was previously addressed in Mukherjee et al. (2017) and Goldsztajn et al. (2018, 2022) in the context of single-server models by assuming that the number of servers can be right-sized on the fly to match the load of the system. However, in the latter papers the dispatching rule remains the same at all times because the right-sizing mechanism alone is sufficient to maintain small queues, by adjusting the number of servers. Different from these right-sizing mechanisms, the learning scheme of SLTA modifies the parameters of the dispatching rule over time to maximize the overall utility of the system.

1.3. Outline of the Paper

In Section 2, we introduce some of the notation used throughout the paper and we prove the upper bound for the mean stationary overall utility. In Section 3, we specify the JLMU and SLTA policies, and we state their asymptotic optimality with respect to the mean stationary overall utility. In Section 4, we present several results that pertain to the asymptotic transient behavior of these two policies and that are used to establish their asymptotic optimality. To characterize the asymptotic behavior of JLMU and SLTA, we construct systems of different sizes on a common probability space in Section 5, where we also prove relative compactness results. Limit theorems for the transient regimes of JLMU and SLTA are proved in Sections 6 and 7, respectively, and the asymptotic optimality of these policies is established in Section 8. Some proofs are deferred to Online Appendices A, B and C.

2. Performance Upper Bound

In this section, we define some of the notation used throughout the paper and we prove the upper bound for the mean stationary overall utility. In Section 2.1, we introduce two descriptors for specifying the state of the system, and we define the overall utility function. In Section 2.2, we present the optimization problem used to derive the upper bound for the mean stationary overall utility. In Section 2.3, we construct a solution of this problem explicitly, and in Section 2.4, we use the constructed solution to formulate the upper bound for the mean stationary overall utility, which we prove in Section 2.5.

2.1. Basic Notation

Consider a system with m classes of server pools. All the tasks sharing a server pool are executed in parallel and the execution times do not depend on the class of the server pool or the number of tasks currently contending for service. Nevertheless, associated with each server pool is a concave utility function which does depend on the class of the server pool and the number of tasks sharing it. For example, these functions can be used to model the overall quality-of-service provided to streaming tasks sharing an underlying resource with a fixed capacity. The objective is to assign the incoming tasks to the various server pools to maximize the aggregate utility of all the server pools in stationarity.

The number of server pools is denoted by n, and the number and fraction of server pools of class i are denoted by $A_n(i)$ and ${\alpha }_n(i)=A_n(i)/n$, respectively. We assume that tasks arrive as a Poisson process of intensity $n\lambda$ with independent and identically distributed service times of mean $1/\mu$, and we define ${\mathit{X}}_n(i,k)$ as the number of tasks in server pool k of class i; bold symbols are used in the paper to indicate time dependence. Server pools of the same class that have the same number of tasks are exchangeable; thus, we usually consider a different state descriptor. Specifically, we let

denote the fraction of server pools which are of class i and have at least j tasks. The values of ${\mathit{X}}_n$ and ${\mathit{q}}_n$ at a given time are referred as the occupancy state or task assignment.

The concave utility function associated with server pools of class i is denoted by u_i, and the overall utility of the system is defined as the aggregate utility of all the server pools normalized by the number of server pools. More precisely, we let

Note that ${\mathit{q}}_n(i,j)-{\mathit{q}}_n(i,j+1)$ is the fraction of server pools of class i with j tasks. Thus, the overall utility may equivalently be expressed as

Although the overall utility function is generally not linear as a function of ${\mathit{X}}_n$, it is always linear as a function of ${\mathit{q}}_n$, as shown by the previous expression.

The total number of tasks in the system, normalized by the number of server pools, can be expressed in terms of the occupancy state ${\mathit{q}}_n$ as follows:

The quantity $j\left[{\mathit{q}}_n(i,j)-{\mathit{q}}_n(i,j+1)\right]$ represents the number of tasks in server pools of class i with exactly j tasks, normalized by the total number of server pools. Hence, ${\mathit{s}}_n$ indeed corresponds to the normalized total number of tasks.

Throughout the paper, we write $P$ and $E$ to denote the probability and expectation with respect to a given probability measure. If $\rho ≔\lambda /\mu$ denotes the normalized offered load, then the stationary distribution of the total number of tasks is Poisson with mean $n\rho$ due to the infinite-server dynamics of the system. Thus, $E\left[{\mathit{s}}_n\right]=\rho$ in stationarity, for any task assignment policy.

2.2. Optimization Problem

Based on the previous statements, we now formulate an optimization problem that yields an upper bound for the mean stationary overall utility:

To see that the optimum of (2) yields an upper bound for the mean stationary overall utility, consider any policy such that ${\mathit{q}}_n$ has a stationary distribution. We assume that the policy is such that the evolution of the system over time can be described by a Markov process, with a possibly uncountable state space, that has a stationary distribution. Let q_n be a random variable with the stationary distribution of ${\mathit{q}}_n$ and define $E\left[q_n\right]$ as the sequence whose (i, j) element is $E\left[q_n(i,j)\right]$. Observe that

Indeed, the utility functions are concave, so for each i there exists $j_i\in \mathbb{N}$ such that $u_i(j)$ and $u_i(k)$ have the same sign if $j,k>j_i$. Therefore, (3) follows from Tonelli’s theorem. In addition, $E\left[q_n\right]$ satisfies the constraints of (2) because the total number of tasks in stationarity has mean $n\rho$. Thus, $E\left[u(q_n)\right]$ is upper bounded by the optimum of (2).

2.3. Structure of an Optimal Solution

For brevity, we refer to an optimizer of (2) as an optimal task assignment; the term optimal fractional task assignment would be more appropriate since the offered load $n\rho$ may not be integral. In this section we define a ranking of the server pools that can be used to construct an optimal task assignment. For this purpose, consider the sets

A server pool has coordinates $(i,j)\in {\mathcal{I}}_{+}$ if its class is i and it has precisely $j-1$ tasks; for example, in Figure 2, server pool A of class 1 has coordinates (1, 4), and both server pools of class 3 have coordinates (3, 1). Because server pools with the same coordinates are statistically identical, we may focus on ranking coordinates rather than server pools.

Formally, we define a total order on ${\mathcal{I}}_{+}$ that gives precedence to coordinates associated with larger marginal utilities. The marginal utility of a server pool of class i with j tasks is denoted by $\mathrm{\Delta }(i,j)≔u_i(j+1)-u_i(j)$ and represents the change in the utility function of such a server pool if it receives an additional task. The marginal utility of the coordinates (i, j) is just $\mathrm{\Delta }(i,j-1)$, the marginal utility of server pools of class i with $j-1$ tasks.

Consider the dictionary order $\prec$ on $\mathcal{I}$, defined by

We obtain a total order $⊲$ on ${\mathcal{I}}_{+}$ by writing $(i_1,j_1)⊲(i_2,j_2)$ if and only if one of the following conditions holds:

In particular, the marginal utility of server pools with coordinates (i₁, j₁) is smaller than or equal to that of server pools with coordinates $(i_2,j_2)⊳(i_1,j_1)$. The dictionary order is used to break the tie when both coordinates are associated with the same marginal utility, but a different tie breaking rule could be used instead.

Consider the task assignment $q_n^{*}$ defined by

In Section 2.5, we prove that $q_n^{*}$ constitutes an optimal task assignment if ${\sigma }_n^{*}$ is defined as the unique element of ${\mathcal{I}}_{+}$ that satisfies

For the uniqueness of ${\sigma }_n^{*}$, the number of terms in the summations on both sides of (4) increases as the ranking of ${\sigma }_n^{*}$ becomes worse and that the summation on the right has exactly one more term than the summation on the left.

2.3.1. Numerical Examples.

Figure 3 depicts the optimal task assignments obtained through (4) for two sets of utility functions and different values of ρ. If n is such that $n{\alpha }_n(i)$ and $n\rho$ are integers for all i, then the plots can be interpreted as sets of n adjacent columns, where each column represents a server pool and the colored portion of a column indicates the number of tasks sharing the server pool, as in the diagram of Figure 2. The thick vertical lines separate the server pool classes and the quantities $q_n^{*}(i,j)$ can be read off by rotating the plots.

The left plot corresponds to utility functions of the form $u_i(x)=xg(r(i)/x)$, with g a concave and increasing function. These utility functions can be used to model the aggregate quality-of-service provided to streaming tasks sharing a single server pool. The quantity r(i) represents the total amount of resources in a server pool of class i and $g(r(i)/x)$ models the quality-of-service provided to a single task when the server pool is shared by x tasks and each task gets a fraction $r(i)/x$ of the total resources. For a given ρ, the left plot of Figure 3 depicts server pools with roughly $r(i)/c(r,\rho )$ tasks, with $c(r,\rho )$ a normalizing constant that does not depend on g; for some values of ρ, all server pools have exactly $r(i)/c(r,\rho )$ tasks, but in other cases some of these numbers are rounded. This behavior is explained by noting that the derivative of $u_i(x)$ can be expressed as a function of $r(i)/x$, thus the occupancy levels $r(i)/c(r,\rho )$ equalize the marginal utilities.

In the left plot of Figure 3, the occupancy levels of the various server pools maintain approximately fixed ratios as ρ increases. The right plot shows a completely different behavior: as ρ increases from 0 to 1.25, only the occupancy of server pools of class 2 grows, but from 1.25 to 6.25, only the occupancy of server pools of class 3 grows, and eventually exceeds the occupancy of server pools of class 2. Furthermore, as ρ increases beyond 7.5, only the occupancy of server pools of class 1 increases.

2.4. Performance Upper Bound

We now state the upper bound for the mean stationary overall utility. The proof uses the remarks of Sections 2.2 and 2.3 and is given in Section 2.5.

Theorem 1.

Consider any task assignment policy such that the occupancy process ${\mathit{q}}_n$ has a stationary distribution and let q_n be a random variable distributed as this stationary distribution. Then

$$E\left[u(q_n)\right]\le u\left(q_n^{*}\right).$$

In Section 3.4, we establish that the upper bound is asymptotically achievable when service times are exponentially distributed. In particular, we will see that JLMU achieves the upper bound of Theorem 1 as the number of server pools grows large; recall that JLMU is generally not optimal in the prelimit, not even for exponential service times. Moreover, we will establish that SLTA also achieves the upper bound asymptotically while relying on considerably less state information.

2.5. Proof of the Upper Bound

Recall from Section 2.2 that the optimum of (2) is an upper bound for the mean stationary overall utility. Therefore, we only need to prove that the task assignment $q_n^{*}$ defined in (4) is an optimizer of (2). For this purpose, we introduce the following definition. We say that a sequence $q\in {\mathbb{R}}^{\mathcal{I}}$ is eventually zero if there exists k > 0 such that $q(i,j)=0$ for all i and j > k. The following lemma implies that $q_n^{*}$ is an optimizer of (2) if we impose the additional constraint that the solution must be eventually zero.

Lemma 1.

If q satisfies the constraints of (2) and is eventually zero, then $u(q)\le u(q_n^{*})$.

Proof.

Because q is eventually zero, it is possible to write

$$\begin{array}{ll}u(q)\hfill & ={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{j=0}^{\infty }{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^m{u}_i}(0)q(i,0)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^{\infty }\mathrm{\Delta }}}(i,j-1)q(i,j)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^m{u}_i}(0){\alpha }_n(i)+{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathrm{\Delta }}(i,j-1)q(i,j).\hfill \end{array}$$(6)

In the last expression, the terms of the summation are ordered with respect to $⊲$ and in particular in nonincreasing order of the marginal utilities $\mathrm{\Delta }(i,j-1)$. The task assignment $q_n^{*}$ is obtained by choosing the coefficients q(i, j) so that the first coefficients are maximal while all the coefficients add up to ρ. Thus, $q_n^{*}$ maximizes the right-hand side of (6). □

We now provide a solution of (2), without imposing any additional constraints.

Proposition 1.

The task assignment $q_n^{*}$ defined in (4) is an optimizer of (2).

Proof.

By Lemma 1, it suffices to prove that $u(q)\le u(q_n^{*})$ for each q that is not eventually zero and satisfies the constraints of (2). Next, we fix one such sequence q and we construct an eventually zero sequence z such that $u(q)\le u(z)$ and z satisfies the constraints of (2). Then $u(q)\le u(z)\le u(q_n^{*})$ by Lemma 1, as desired.

Choose $k\in \mathbb{N}$ such that ${\alpha }_n(i)k>\rho$ for all i. For each i, we define z(i, j) iteratively, by

$$z(i,j)≔\{\begin{array}{ll}\min \left\{q(i,j)+\sum _{l=k+1}^{\infty }q(i,j)-\sum _{l=0}^{j-1}\left[z(i,l)-q(i,l)\right],{\alpha }_n(i)\right\}\hfill & \text{if}j\le k,\hfill \\ 0\hfill & \text{if}j>k.\hfill \end{array}$$

Informally, each coefficient q(i, j) can be regarded as a container with capacity ${\alpha }_n(i)$, as shown in Figure 4. For each i, the sequence $q(i,·)$ is transformed into $z(i,·)$ in two steps: first we remove all the mass from the coefficients q(i, j) with j > k, and then we place this mass on the coefficients q(i, j) with $j\le k$. In the latter step, we start with the first coefficient, placing as much mass as possible without exceeding the capacity ${\alpha }_n(i)$. The remainder of mass is placed in the following coefficients in the same fashion, and in increasing order of j. Observe that all the mass will have been placed right after the coefficient q(i, k) is done because q satisfies the constraints of (2) and ${\alpha }_n(i)k>\rho$. Furthermore, the following property holds:

$${\displaystyle \sum _{j=0}^k\left[z(i,j)-q(i,j)\right]}={\displaystyle \sum _{j=k+1}^{\infty }q}(i,j)\text{for}\text{all}i.$$(7)

The overall utility of the task assignment q satisfies

$$\begin{array}{l}u(q)\hfill \\ ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^{\infty }{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=k+1}^{\infty }{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ \le {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ +{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=k+1}^{\infty }\left[u_i(k)+(j-k)\mathrm{\Delta }(i,k)\right]}}\left[q(i,j)-q(i,j+1)\right]\hfill \\ ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]+{\displaystyle \sum _{i=1}^m{u}_i}(k)q(i,k+1)+{\displaystyle \sum _{i=1}^m\mathrm{\Delta }}(i,k){\displaystyle \sum _{j=k+1}^{\infty }q}(i,j)\hfill \\ ={\displaystyle \sum _{i=1}^m{u}_i}(0){\alpha }_n(i)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k\mathrm{\Delta }}}(i,j-1)q(i,j)+{\displaystyle \sum _{i=1}^m\mathrm{\Delta }}(i,k){\displaystyle \sum _{j=k+1}^{\infty }q}(i,j)\hfill \\ =u(z)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k\mathrm{\Delta }}}(i,j-1)\left[q(i,j)-z(i,j)\right]+{\displaystyle \sum _{i=1}^m\mathrm{\Delta }}(i,k){\displaystyle \sum _{j=k+1}^{\infty }q}(i,j).\hfill \end{array}$$

For the last step, recall that z is eventually zero, so u(z) can be computed as in (6). Note that $q(i,j)\le z(i,j)$ and $\mathrm{\Delta }(i,j-1)\ge \mathrm{\Delta }(i,k)$ if $j\le k$. Therefore,

$$u(q)\le u(z)+{\displaystyle \sum _{i=1}^m\mathrm{\Delta }}(i,k){\displaystyle \sum _{j=1}^k\left[q(i,j)-z(i,j)\right]}+{\displaystyle \sum _{i=1}^m\mathrm{\Delta }}(i,k){\displaystyle \sum _{j=k+1}^{\infty }q}(i,j)=u(z)\le u(q_n^{*}).$$

The middle equality and last inequality follow from (7) and Lemma 1, respectively. □

The proof of Theorem 1 follows easily from Proposition 1.

Proof of Theorem 1.

As indicated at the end of Section 2.2, the optimum of (2) upper bounds the mean stationary overall utility $E\left[u(q_n)\right]$. Thus, it follows from Proposition 1 that $u\left(q_n^{*}\right)$ upper bounds the mean stationary overall utility. □

3. Load Balancing Policies

In this section we describe the load balancing policies considered in the paper, and we state their asymptotic optimality with respect to the mean stationary overall utility when service times are exponential. In Sections 3.1 and 3.2, we specify JLMU and SLTA, respectively. In Section 3.3, we define stochastic models, based on continuous-time Markov chains, for the analysis of both policies. In Section 3.4, we state the asymptotic optimality result.

Before proceeding, it is illustrative to draw an analogy between the setting considered in this paper and the load balancing literature for systems of parallel single-server queues, where the primary objective is to minimize queueing delay. The natural policy for the setting considered in this paper is JLMU, while the natural policy for minimizing queueing delay in systems of parallel single-server queues is JSQ. The deployment of these policies involves a considerable communication overhead, or storing and managing a significant amount of state information. In the setting considered in this paper, SLTA provides a asymptotically optimal performance for exponential service times and uses substantially less state information than JLMU. From this perspective, SLTA is the counterpart of JIQ in the load balancing literature for systems of parallel single-server queues (Lu et al. 2011, Stolyar 2015).

3.1. Join the Largest Marginal Utility (JLMU)

JLMU assigns every new task to a server pool that currently has the best ranked coordinates, thus also the largest marginal utility. Formally, define

for each occupancy state q. The maximum is taken with respect to $⊲$, and the condition $q(i,j-1)>q(i,j)$ implies that some server pool of class i has precisely $j-1$ tasks. If q_n is the occupancy state right before a task arrives, then JLMU assigns the task to a server pool of class ${\sigma }_i\left(q_n\right)$ with exactly ${\sigma }_j\left(q_n\right)-1$ tasks.

The coordinates obtained through (8) correspond to server pools with the largest marginal utility by definition of $⊲$. In addition, observe that the dictionary order is used to break ties between coordinates associated with the same marginal utility. If two server pools have the same coordinates, then it does not matter which of them is assigned the new task since they are statistically identical. For definiteness, we postulate that the tie is broken uniformly at random.

If all the server pools have the same utility function, then JLMU reduces to JSQ, and the overall utility is a Schur-concave function of ${\mathit{X}}_n$. If in addition the service times are exponential, then the stochastic optimality properties proved in Sparaggis et al. (1993) and Menich and Serfozo (1991) for JSQ imply that JLMU maximizes the mean stationary overall utility in this homogeneous setting. It might be natural to expect that the optimality with respect to the mean stationary overall utility extends to the heterogeneous setting. We refute this, however, in Section 8.3, where we construct a heterogeneous system for which JLMU is strictly suboptimal. The constructed example also hints at the underlying reasons for the suboptimality in the heterogeneous case. Essentially, instead of always assigning incoming tasks greedily, such that the increase in the overall utility is maximal, it is sometimes advantageous to dispatch the new tasks conservatively, to hedge against pronounced drops of the overall utility that may be caused by a quick succession of departures. The right balance between greedy and conservative actions depends intricately on the utility functions, but we prove that JLMU is always asymptotically optimal for exponential service times, regardless of the specific set of utility functions; this result is stated formally in Section 3.4.

3.2. Self-Learning Threshold Assignment (SLTA)

JLMU relies on complete information about the number of tasks per server pool, which could be impractical in large-scale deployments. In contrast, SLTA only requires to store at most two bits per server pool, which is considerably less state information. To specify this policy, we need to describe its two components. Namely, the dispatching rule, for assigning the incoming tasks to the server pools, and the learning scheme, for dynamically adjusting a set of thresholds that the dispatching rule uses.

Consider the unique enumeration $\{(i_k,j_k):k\ge 1\}$ of ${\mathcal{I}}_{+}$ such that $(i_k,j_k)⊳(i_{k+1},j_{k+1})$ for all k. Given $r\ge 1$, we define a set of thresholds by

Recall the optimal task assignment defined at the end of Section 2.3. The learning scheme keeps an estimate $(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$ of the coordinates ${\sigma }_n^{*}$, which depend on the typically unknown offered load of the system. The index ${\mathit{r}}_n$ determines this estimate and is used to compute thresholds from the previous expression, which are in turn used to assign tokens to the server pools. Specifically, a server pool of class i with exactly $j-1$ tasks has

• A green token if $j-1<{\ell }_i({\mathit{r}}_n)$,

• A yellow token if $i=i_{{\mathit{r}}_n}$ and $j-1\le {\ell }_i({\mathit{r}}_n)$.

The first condition is equivalent to ${\mathit{r}}_n>1$ and $(i,j)⊳(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$, and the second condition is equivalent to $i=i_{{\mathit{r}}_n}$ and $(i,j)\mathrm{⊵}(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$. Also, a server pool can have both a green and a yellow token at the same time. As indicated in Figure 5, a larger increase in the overall utility is obtained by dispatching tasks to server pools with green tokens first, then to server pools with yellow tokens and only afterward to server pools without tokens.

The tokens are used by the dispatching rule. Specifically, when a task arrives, it is assigned to a server pool according to the following criteria.

• In the presence of green tokens of class $i\ne i_{{\mathit{r}}_n-1}$, the dispatcher picks one of these green tokens uniformly at random, and if only green tokens of class $i=i_{{\mathit{r}}_n-1}$ remain, then one of these is picked. Then the task is sent to the corresponding server pool.

• In the presence of only yellow tokens, the dispatcher picks a yellow token uniformly at random and sends the task to the corresponding server pool.

• Otherwise, the task is sent to a server pool chosen uniformly at random.

If $(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$ are the coordinates ${\sigma }_n^{*}$ defined in (5), then this dispatching rule drives the occupancy state of the system toward the optimal task assignment $q_n^{*}$ specified in (4).

The learning scheme aims at finding the coordinates ${\sigma }_n^{*}$, which depend on the typically unknown offered load. The learning scheme is parameterized by ${\beta }_n>0$ and adjusts the value of ${\mathit{r}}_n$ at certain arrival epochs, in steps of one unit. Specifically, when a task arrives, the learning scheme acts only under the following circumstances.

• If the system has at least $n{\beta }_n$ green tokens and at least one belongs to a server pool of class $i_{{\mathit{r}}_n-1}$, then ${\mathit{r}}_n$ is decremented by one after the task is dispatched.

• If the number of yellow tokens is smaller than or equal to one and there are no other tokens, then ${\mathit{r}}_n$ is incremented by one after the task is dispatched.

Observe that exactly one of the thresholds changes when the value of ${\mathit{r}}_n$ is modified and that this threshold changes by one unit. Also, note that

are the number of green and yellow tokens, respectively.

3.2.1. Comparison with the Homogeneous Case.

When all the server pools are of the same class, SLTA reduces to the load balancing policy studied in Goldsztajn et al. (2021a). In this case, there is a single threshold whose optimal value is simply $\lfloor \rho \rfloor$. When the threshold has this value, the number of green tokens and the total number of tokens are typically small and positive, respectively. On the other hand, when the threshold is below optimal, the total number of tokens tends to be zero, and when the threshold is larger than optimal, the number of green tokens tends to be relatively large. These properties are used to adjust the threshold in an online manner when the offered load is unknown. In few words, the threshold is increased in the absence of tokens, and it is decreased if the number of green tokens is large enough.

In the general case, there are as many thresholds as the number of server pool classes, and the optimal threshold values depend intricately on the utility functions and the offered load. However, the ranking $⊲$ introduced in Section 2.3 makes it possible to express all the thresholds as a function of the coordinates ${\sigma }_n^{*}$. Hence, the optimal thresholds can still be found through a one-directional search, but now in the totally ordered space ${\mathcal{I}}_{+}$. Moreover, because this space is countable, the search can be carried out by adjusting the integral parameter ${\mathit{r}}_n$ until $(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$ reaches the optimal value ${\sigma }_n^{*}$.

The learning scheme of SLTA operates so that all the thresholds are at their optimal values if and only if ${\mathit{r}}_n$ is at its optimal value, and when this happens, the number of green tokens and the total number of tokens are typically small and positive, respectively. When ${\mathit{r}}_n$ is below optimal, all the thresholds are smaller than or equal to their optimal values, and at least one of the thresholds is strictly smaller than optimal; in this case, the total number of tokens tends to be zero. Similarly, when ${\mathit{r}}_n$ is above optimal, all the thresholds are larger than or equal to their optimal values and at least one of the thresholds is above optimal; as a result, the number of green tokens tends to be relatively large. As in the homogeneous case, these observations are used to adjust ${\mathit{r}}_n$ over time. Loosely speaking, in the absence of tokens, ${\mathit{r}}_n$ is increased by one unit, which implies that one of the thresholds is increased by one unit, and when the number of green tokens is large enough, ${\mathit{r}}_n$ is decreased by one unit, and thus one of the thresholds is decreased by one unit.

The dispatching rule and the online learning scheme of SLTA make distinctions between green tokens of class $i_{{\mathit{r}}_n-1}$ and green tokens of any other class. The rationale is that the marginal utility of server pools of class i with ${\ell }_i({\mathit{r}}_n)-1$ tasks is the lowest when $i=i_{{\mathit{r}}_n-1}$; thus, it makes sense to give green tokens of this class the lowest priority for receiving new tasks. Although this may slightly improve performance, it is not crucial. Nevertheless, the differential treatment of class $i_{{\mathit{r}}_n-1}$ simplifies the mathematical analysis of SLTA. In particular, the distinction made in the description of the learning scheme ensures that if the system had yellow tokens, then it will continue to have yellow tokens after ${\mathit{r}}_n$ is decreased, which is used in Remark 1. In addition, the differential treatment by the dispatching rule simplifies the proof of Proposition 6.

3.3. Stochastic Models

If service times are exponentially distributed, then ${\mathit{q}}_n$ and $({\mathit{q}}_n,{\mathit{r}}_n)$ are continuous-time Markov chains when the load balancing policies are JLMU and SLTA, respectively. In either case, the process ${\mathit{s}}_n$ that describes the normalized total number of tasks is defined by (1). Due to the infinite-server dynamics of the system, $n{\mathit{s}}_n$ has the law of an $M/M/\infty$ queue with arrival rate $n\lambda$ and service rate μ.

Let ${\ell }_1$ be the space of absolutely summable sequences in ${\mathbb{R}}^{\mathcal{I}}$, equipped with the norm

Throughout we assume that ${\mathit{s}}_n(0)$ is finite, so ${\mathit{q}}_n(0)$ takes values in ${\ell }_1$. As a result, if we let $F_n≔\{k/n:0\le k\le n\}$, then ${\mathit{q}}_n$ takes values in the set

If the load balancing policy is JLMU, then the state space of ${\mathit{q}}_n$ is defined as the subset S_n of Q_n that is reachable from an empty occupancy state. If the load balancing policy is SLTA, then the state space of $({\mathit{q}}_n,{\mathit{r}}_n)$ is the subset S_n of $Q_n\times \{r\in \mathbb{N}:r\ge 1\}$ that is reachable from an empty occupancy state with ${\mathit{r}}_n=1$.

The notation used for the processes ${\mathit{s}}_n$ and ${\mathit{q}}_n$, for the state space S_n, and for some other objects that will be defined later, is exactly the same for JLMU and SLTA, but we always indicate which policy is being considered.

3.4. Asymptotic Optimality

Throughout the rest of the paper, we assume that there exist constants $\alpha (i)\in (0,1)$ and a random variable q₀ such that the following limits hold:

In analogy with (5) and (4), we consider the unique ${\sigma }_{*}=\left(i_{r_{*}},j_{r_{*}}\right)\in {\mathcal{I}}_{+}$ such that

and we define an occupancy state $q_{*}$ in terms of ${\sigma }_{*}$ by

In Section 8.1, we establish that ${\mathit{q}}_n$ and $({\mathit{q}}_n,{\mathit{r}}_n)$ have a unique stationary distribution for all n when the assignment policies are JLMU and SLTA, respectively. The following theorem is proved in Section 8.2 and implies that both policies are asymptotically optimal if the marginal utilities are bounded, and the service times are exponentially distributed; we also require that (9) holds, and we impose some mild technical assumptions, to be stated in Section 4.2.1. The condition on the marginal utilities always holds when the utility functions are nondecreasing due to the concavity of these functions.

Theorem 2.

Suppose that service times are exponentially distributed. Also, if the load balancing policy is SLTA, assume that the assumptions of Section 4.2.1 hold and that

$$\underset{n\to \infty }{\lim }{\beta }_n=0\mathit{and}\underset{n\to \infty }{\lim \inf }{n}^{{\gamma }_0}{\beta }_n>0\mathit{for}\mathit{some}{\gamma }_0\in (0,1/2).$$

The following statements hold in the asymptotic regime described by (9).

• (a) Suppose that JLMU is used and let q_n have the stationary distribution of ${\mathit{q}}_n$. Then the random variables q_n converge weakly in ${\ell }_1$ to $q_{*}$.

• (b) Suppose that SLTA is used and let (q_n, r_n) have the stationary distribution of $({\mathit{q}}_n,{\mathit{r}}_n)$. The random variables (q_n, r_n) converge weakly in ${\ell }_1\times \mathbb{N}$ to $(q_{*},r_{*})$.

• (c) Furthermore, if the load balancing policy is either JLMU or SLTA and the marginal utilities are bounded, then the random variables $u(q_n)$ are uniformly integrable and

  $$\underset{n\to \infty }{\lim }E\left[u(q_n)\right]=u(q_{*})=\underset{n\to \infty }{\lim }u\left(q_n^{*}\right).$$

The claims concerning the stationary overall utilities are proved using (a) and (b), as well as the fact that u(q) is a bounded linear functional of $q\in {\ell }_1$. To establish (a) and (b), we first use drift analysis to prove that the random variables in (a) and (b) are tight in ${\ell }_1$ and ${\ell }_1\times \mathbb{N}$, respectively. Then (a) is established through an interchange of limits argument based on a fluid limit and a global asymptotic stability result for the fluid dynamics; these two results are stated in Theorems 5 and 4, respectively. A different type of argument is used to prove (b). Namely, the fluid limit step is circumvented, and Theorem 6 serves as the counterpart of Theorems 4 and 5, as illustrated in Figure 6.

Deriving a fluid limit for a SLTA system would be inherently difficult due to the intricate interdependence between the dispatching rule and the learning scheme and because the actions of the learning scheme are triggered by excursions of the occupancy process that have vanishing size. We deal with these challenges using a methodology of Goldsztajn et al. (2021a) to derive the fluid approximation of Theorem 6, which consists of asymptotic bounds, over arbitrarily long intervals of time, for the occupancy state and the thresholds. As noted earlier, this fluid approximation serves as a counterpart of both the fluid limit and the global asymptotic stability results for JLMU.

3.4.1. Simulation Experiments.

The asymptotic optimality of JLMU and SLTA is illustrated by Table 1, which shows estimates of the mean stationary overall utility $E[u(q_n)]$ for simulation experiments with different values of n. All the estimates correspond to systems with two server pool classes of equal size and utility functions of the form $u_i(x)=x\log (r(i)/x)$ for $r=(20,30)$. Also, two different values of ρ are considered, so that the optimal task assignment defined in (4) takes two distinct forms. For $\rho =9.75$, the optimal task assignment is such that all server pools of class 1 have 8 tasks, half of the server pools of class 2 have 11 tasks, and the other half have 12 tasks. For ρ = 10, the optimal task assignment is such that all server pools of class 1 have 8 tasks and all server pools of class 2 have 12 tasks.

Table 1. Results of Simulation Experiments with Different Values of n

All the server pools of the same class have the same number of tasks when ρ = 10, and thus we say that the optimal task assignment $q_n^{*}$ is integral; that is, the optimal task assignment is integral if $q_n^{*}({\sigma }_n^{*})=0$. In contrast, the optimal task assignment is fractional when $\rho =9.75$ because server pools of class 2 may have either 11 or 12 tasks. Server pools of class 1 behave similarly in the fractional and integral settings: almost all of the time all server pools of class 1 have precisely eight tasks when n is moderately large. However, the behavior of server pools of class 2 depends on the setting. In the fractional case, server pools typically have 11 or 12 tasks and the fractions of server pools with 11 and 12 tasks oscillate around a half. In the integral case, server pools typically have 12 tasks but a small number of server pools sometimes have 11 or 13 tasks instead. The aggregate utility of the system decreases by $\mathrm{\Delta }(2,11)$ whenever a server pool of class 2 goes from 12 to 11 tasks and increases by the same quantity when the server pool goes from 11 to 12 tasks. Therefore, the contributions to the average overall utility of the oscillations observed in the fractional case roughly balance each other. In the integral case, the aggregate utility increases by $\mathrm{\Delta }(2,12)$, instead of $\mathrm{\Delta }(2,11)$, if a server pool of class 2 goes from 12 to 13 tasks. Therefore, in the integral case, the contributions to the average overall utility of class 2 server pools that drop to 11 tasks or reach 13 tasks are amplified by different marginal utilities. As a result, the mean stationary overall utility $E[u(q_n)]$ is closer to the upper bound $u(q_n^{*})$ in the fractional case, as reflected by the estimates in Table 1.

Although there is a difference between the fractional and integral settings, in both settings, the empirical mean of the overall utility $u({\mathit{q}}_n)$ is extremely close to the upper bound $u(q_n^{*})$ across all the values of n listed in Table 1. Furthermore, the deviation of the empirical mean from the upper bound approaches zero as n increases in both cases. We also observe that the empirical mean of $u({\mathit{q}}_n)$ is almost the same for JLMU and SLTA in all the experiments and particularly for the largest values of n.

A final remark on the simulation experiments is that the empirical mean of $u({\mathit{q}}_n)$ is slightly larger than $u(q_n^{*})$ in a few of the experiments within the fractional setting: both for JLMU and SLTA when n = 350 and just for JLMU when n = 450. It may be checked that the statement of Theorem 1 still holds if the stationary expectation sign is replaced by a time average and the upper bound is computed through (4) and (5) but with ρ replaced by the time average of ${\mathit{s}}_n$; the proof does not change. The value of the upper bound when ρ is replaced by the time average of ${\mathit{s}}_n$ is displayed in Table 1, and in all the experiments, the empirical mean of $u({\mathit{q}}_n)$ is indeed smaller than this empirical upper bound. Thus, the experiments where the empirical mean of $u({\mathit{q}}_n)$ slightly exceeds the upper bound $u(q_n^{*})$ are an indication of how close the performance of JLMU and SLTA is to optimal.

4. Approximation Theorems

In this section, we assume exponential service times and we state several results used to prove Theorem 2. In Section 4.1, we specify a fluid model of a JLMU system, based on differential equations, and we state some properties of this model. In Section 4.2, we state limit theorems that characterize the asymptotic transient behavior of JLMU and SLTA.

4.1. Fluid Model of JLMU

Consider a large-scale system where the load balancing policy is JLMU and assume that $\alpha (i)$ is the fraction of server pools of class i. Then the occupancy state of the system remains within the set

The evolution of the occupancy state of this large-scale system can be modeled through the system of differential equations introduced in the following definition.

Definition 1.

We say that $\mathit{q}:[0,\infty )\to Q$ is a fluid trajectory if the coordinate functions $\mathit{q}(i,j)$ are absolutely continuous for all $(i,j)\in \mathcal{I}$ and the following conditions hold almost everywhere with respect to the Lebesgue measure:

$$\dot{\mathit{q}}(i,j)=\mathrm{\Lambda }(\mathit{q},i,j)-\mu j\left[\mathit{q}(i,j)-\mathit{q}(i,j+1)\right],$$(12a)

$$\mathrm{\Lambda }(\mathit{q},i,j)\ge 0\text{for}\text{all}(i,j)\in {\mathcal{I}}_{+},$$(12b)

where $\mathrm{\Lambda }:Q\times {\mathcal{I}}_{+}\to \mathbb{R}$ is defined by

$$\mathrm{\Lambda }(q,i,j)≔\{\begin{array}{ll}0\hfill & \text{if}(i,j)⊲\sigma (q),\hfill \\ \mu j\left[\alpha (i)-q(i,j+1)\right]\hfill & \text{if}(i,j)⊳\sigma (q),\hfill \\ \lambda -\sum _{(k,l)⊳\sigma (q)}\mu l\left[\alpha (k)-q(k,l+1)\right]\hfill & \text{if}(i,j)=\sigma (q).\hfill \end{array}$$

In the latter definition, λ is the arrival rate of tasks normalized by the number of server pools, μ is the service rate of tasks, and $\mathit{q}(i,j)$ represents the fraction of server pools that are of class i and have at least j tasks. Thus, the system of Differential Equations (12) has a simple interpretation. The right-most term of (12a) corresponds to the departure rate of tasks from server pools of class i with exactly j tasks, and $\mathrm{\Lambda }(\mathit{q},i,j)$ represents the arrival rate of tasks to server pools that belong to class i and have precisely $j-1$ tasks. The definition of Λ is motivated by the following remarks.

• Server pools of class i with exactly $j-1$ tasks are not assigned additional tasks if $(i,j)⊲\sigma (\mathit{q})$. Hence, we should have $\mathrm{\Lambda }(\mathit{q},i,j)=0$ in this case.

• All server pools of class i have at least j tasks if $(i,j)⊳\sigma (\mathit{q})$. Therefore, $\mathrm{\Lambda }(\mathit{q},i,j)$ should be equal to the last term of (12a) in this case, because $\mathit{q}(i,j)$ is at its maximum value and thus its derivative should be zero.

• The total arrival rate of tasks normalized by the number of server pools is equal to λ, and this determines the value of $\mathrm{\Lambda }(\mathit{q},i,j)$ when $(i,j)=\sigma (\mathit{q})$.

4.1.1. Properties of Fluid Trajectories.

The two results stated here are proved in Section 6.1. The first one is a uniqueness theorem for the solutions of (12). Existence is ensured by Theorem 5 of Section 4.2.

Theorem 3.

Fix an initial condition $q\in Q$. If there exists a fluid trajectory $\mathit{q}$ such that $\mathit{q}(0)=q$, then this fluid trajectory is unique.

To prove this theorem, we first show that all fluid trajectories satisfy an infinite system of integral equations, stated using Skorokhod one-dimensional reflection mappings. The theorem is then proved using a Lipschitz property of these mappings and a uniqueness result for certain Kolmogorov backward equations.

Besides uniqueness of solutions of (12), we also establish that there exists a unique equilibrium point and that this equilibrium point is globally asymptotically stable; that is, all fluid trajectories converge to the unique equilibrium over time.

Theorem 4.

Let $q_{*}$ be as in (11). Then $q_{*}$ is the unique equilibrium of (12). Furthermore, all fluid trajectories converge to $q_{*}$ in ${\ell }_1$ over time.

Recall that (11) is the counterpart of (4), which is used to formulate the upper bound for the mean stationary overall utility provided in Theorem 1. It is not difficult to check that $u\left(q_n^{*}\right)\to u\left(q_{*}\right)$ as n grows large, which hints at the asymptotic optimality of JLMU.

4.2. Limit Theorems

In Section 5.1, we construct the processes defined in Section 3.3 on a common probability space $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ for all n, in such a way that the sample paths of the occupancy processes lie in the space $D_{{\ell }_1}[0,\infty )$ of càdlàg functions with values in ${\ell }_1$, which we endow with the topology of uniform convergence over compact sets. This construction is used to prove limit theorems that characterize the asymptotic transient behavior of JLMU and SLTA. Before stating these theorems, we introduce some mild technical assumptions.

4.2.1. Technical Assumptions.

As indicated earlier, we assume that (9) holds with q₀ a random variable that takes values in Q and represents the limiting initial occupancy state. The initial number of tasks in the limit, normalized by the number of server pools, is defined as

If the load balancing policy is SLTA, then we assume that the first inequality in (10) is strict and that there exists a constant ${\gamma }_0\in (0,1/2)$ such that

These assumptions are used to prove that the learning scheme reaches an equilibrium in all large enough systems with probability one. Finally, we adopt the following assumptions about the initial state of the system: There exists a random variable $R\ge 1$ such that

for all n with probability one. We impose (15b) just to simplify the analysis; this property always holds after a certain time, which depends on the initial state of the system.

Remark 1.

Property (15b) is preserved by arrivals and departures; thus, it holds at all times provided that it holds at time zero. Furthermore, every new task is sent to a server pool with coordinates $(i,j)\mathrm{⊵}(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$ if the number of tokens is positive right before the arrival. Hence, (15b) implies that tasks are sent to server pools with coordinates $(i,j)\mathrm{⊵}(i_{{\mathit{r}}_n},j_{{\mathit{r}}_n})$ at all times and for all n with a probability of one.

4.2.2. Statements of the Theorems.

First we state a fluid limit for JLMU, proven in Section 6.2. In view of Theorem 4, this fluid limit implies that, as n grows large, the occupancy processes of JLMU approach functions that converge over time to the unique equilibrium of (12).

Theorem 5.

Suppose that the load balancing policy is JLMU. Then there exists a set of probability one Γ with the following property. If $\omega \in \mathrm{\Gamma }$, then ${\mathit{q}}_n(\omega )$ converges in $D_{{\ell }_1}[0,\infty )$ to the unique fluid trajectory with initial condition $q_0(\omega )$.

Because q₀ is arbitrary, the previous theorem implies that solutions to (12) exist for all initial conditions. Therefore, Theorems 3 and 5 imply that for each initial condition $q\in Q$ there exists a unique fluid trajectory with initial condition q.

The proof of Theorem 5 uses a methodology of Bramson (1998) to prove that, with a probability of one, every subsequence of $\{{\mathit{q}}_n:n\ge 1\}$ has a further subsequence that converges uniformly over compact sets with respect to a metric for the product topology of ${\mathbb{R}}^{\mathcal{I}}$. Then we show that this convergence in fact holds with respect to ${\Vert ·\Vert }_1$ and that the limits of convergent subsequences are fluid trajectories, also with a probability of one.

The counterpart of Theorems 4 and 5 for SLTA is the following result. The proof is provided in Section 7.2 and is based on a methodology of Goldsztajn et al. (2021a).

Theorem 6.

Suppose that the load balancing policy is SLTA and let ${\sigma }_{*}$ and $r_{*}$ be as in (10). There exist ${\tau }_{\text{eq}}:[0,\infty )\to \mathbb{R}$ and a set of probability one Γ with the following property. If $\omega \in \mathrm{\Gamma }$ and $T\ge \tau >{\tau }_{\text{eq}}(s_0(\omega ))$, then the next limits hold:

$$\underset{n\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }|{\mathit{r}}_n(\omega )-r_{*}|=0,$$(16a)

$$\underset{n\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }{n}^{\gamma }\left|\alpha (i)-{\mathit{q}}_n(\omega ,t,i,j)\right|=0\mathit{if}(i,j)⊳{\sigma }_{*}\mathit{and}\gamma \in [0,1/2),$$(16b)

$$\underset{n\to \infty }{\lim \sup }\underset{t\in [\tau ,T]}{\sup }{\displaystyle \sum _{(i,j)⊲{\sigma }^{*}}{\mathit{q}}_n}(\omega ,t,i,j){\text{e}}^{\mu (t-\tau )}\le c(\omega ,\tau ),$$(16c)

where $c(\omega ,\tau )$ can be expressed in terms of ρ, τ, and $s_0(\omega )$.

5. Strong Approximations

In this section, we construct the processes defined in Section 3.3 on a common probability space for all n. In addition, we prove that $\{{\mathit{q}}_n:n\ge 1\}$ is almost surely relatively compact in $D_{{\ell }_1}[0,\infty )$ both for JLMU and SLTA. The construction of the processes is carried out in Section 5.1, and the relative compactness results are provided in Section 5.2.

5.1. Coupled Construction of Sample Paths

Consider the following stochastic processes and random variables.

• Driving Poisson processes: A collection $\left\{\mathcal{N}\right\}\cup \{{\mathcal{N}}_{\nu }:\nu \in {\mathcal{I}}_{+}\}$ of independent Poisson processes with unit rate, for counting arrivals and departures. These processes are defined on a common probability space $({\mathrm{\Omega }}_D,{\mathcal{F}}_D,{\mathbb{P}}_D)$.

• Selection variables: A family $\{U_k:k\ge 1\}$ of independent random variables, uniformly distributed on $[0,1)$ and defined on a common probability space $({\mathrm{\Omega }}_S,{\mathcal{F}}_S,{\mathbb{P}}_S)$.

• Initial conditions: Sequences $\{{\mathit{q}}_n(0):n\ge 1\}$ and $\{{\mathit{r}}_n(0):n\ge 1\}$ for SLTA of random variables for describing the initial states of the systems, defined on a common probability space $({\mathrm{\Omega }}_I,{\mathcal{F}}_I,{\mathbb{P}}_I)$ and satisfying the assumptions of Section 4.2.1.

Denote the completion of the product probability space of $({\mathrm{\Omega }}_D,{\mathcal{F}}_D,{\mathbb{P}}_D),\hspace{0.17em}({\mathrm{\Omega }}_S,{\mathcal{F}}_S,{\mathbb{P}}_S)$ and $({\mathrm{\Omega }}_I,{\mathcal{F}}_I,{\mathbb{P}}_I)$ by $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$. The processes introduced in Section 3.3 are constructed on the latter space as deterministic functions of the stochastic primitives.

5.1.1. Construction for JLMU.

Let ${\mathcal{N}}_n^{\lambda }(t)≔\mathcal{N}(n\lambda t)$ for each $t\ge 0$ and each n. This quantity will be used to count the number of tasks arriving to the system with n server pools during the interval $[0,t]$. Also, denote the jump times of ${\mathcal{N}}_n^{\lambda }$ by $\{{\tau }_{n,k}:k\ge 1\}$ and define ${\tau }_{n,0}≔0$. For each function $\mathit{q}:[0,\infty )\to Q_n$ and each n, we define two counting processes, for arrivals and departures, denoted ${\mathcal{A}}_n(\mathit{q})$ and ${\mathcal{D}}_n(\mathit{q})$, respectively. The coordinates $(i,j)\in \mathcal{I}$ of these processes are identically zero if j = 0, whereas the other coordinates are defined as follows:

For each n, the functional equation

has a unique solution with a probability of one. More precisely, there exists a set of probability of one ${\mathrm{\Gamma }}_0$ with the following property: For each $\omega \in {\mathrm{\Gamma }}_0$ and each n, there exists a unique càdlàg function ${\mathit{q}}_n(\omega ):[0,\infty )\to Q_n$ that solves (18). This solution can be constructed by forward induction on the jump times of the driving Poisson processes. The assumption ${\mathit{q}}_n(\omega ,0)\in Q_n$ implies that ${\mathit{q}}_n(\omega ,0)$ has finitely many nonzero coordinates and ensures that the constructed solution is defined on $[0,\infty )$ with probability one; that is, the constructed solution does not explode in finite time.

The occupancy processes are defined by extending the previous solutions to Ω, setting ${\mathit{q}}_n(\omega ,t)=0$ for all $t\ge 0$ and all $\omega \notin {\mathrm{\Gamma }}_0$. In addition, we let

and we note that the sample paths of ${\mathcal{A}}_n,\hspace{0.17em}{\mathcal{D}}_n$ and ${\mathit{q}}_n$ lie in $D_{{\ell }_1}[0,\infty )$. Functional Equation (18) can now be rewritten as follows:

This construction endows the processes ${\mathit{q}}_n$ with the intended statistical behavior. The processes ${\mathcal{A}}_n(i,j)$ count the arrivals to server pools of class i with precisely $j-1$ tasks and the processes ${\mathcal{D}}_n(i,j)$ count the departures from server pools of class i with exactly j tasks. Indeed, ${\mathcal{A}}_n(i,j)$ has a jump at the arrival epoch ${\tau }_{n,k}$ if and only if the incoming task should be assigned to a server pool of class i with $j-1$ tasks under the JLMU policy. In addition, the intensity of ${\mathcal{D}}_n(i,j)$ equals the total number of tasks in server pools of class i with exactly j tasks times the rate at which tasks are executed, and this totals the departure rate from server pools of class i with precisely j tasks.

5.1.2. Construction for SLTA.

The processes $({\mathit{q}}_n,{\mathit{r}}_n)$ are constructed to a large extent as in Section 5.1.1 when the load balancing policy is SLTA. The only differences arise in (17a) and (18). Namely, (17a) must be modified to capture the dispatching rule of SLTA, and (18) must be accompanied by another equation for describing the evolution of ${\mathit{r}}_n$.

The counterpart of (17a), with an extra argument $\mathit{r}:[0,\infty )\to \{r\in \mathbb{N}:r\ge 1\}$, is

The functions η_k are defined in Online Appendix B using the selection variables U_k, so that they have the following property. If (q, r) is the value of $({\mathit{q}}_n,{\mathit{r}}_n)$ when the $k^{\text{th}}$ task arrives, then SLTA sends this task to a server pool of class i with precisely $j-1$ tasks if and only if ${\eta }_k(q,r,i,j)=1$. Moreover, ${\eta }_k(q,r,h,l)=0$ for all $(h,l)\ne (i,j)$.

The analog of the functional Equation (18) is

where the sets $I_{n,k}$ and $D_{n,k}$ are defined formally in Online Appendix A. The former set indicates that $(\mathit{q},\mathit{r})$ corresponds to a system with no green tokens and at most one yellow token right before the $k^{\text{th}}$ arrival. The latter set indicates that the number of green tokens is larger than or equal to $n{\beta }_n$ and that at least one of these tokens belongs to a server pool of class $i_{\mathit{r}-1}$ right before the $k^{\text{th}}$ arrival. Also, ${\mathcal{D}}_n(\mathit{q})$ is defined as in (17b).

As in Section 5.1.1, there exists a set of probability of one ${\mathrm{\Gamma }}_0$ with the next property. For each $\omega \in {\mathrm{\Gamma }}_0$ and each n, there exists a unique pair of càdlàg functions ${\mathit{q}}_n(\omega ):[0,\infty )\to Q_n$ and ${\mathit{r}}_n(\omega ):[0,\infty )\to \{r\in \mathbb{N}:r\ge 1\}$ that solve (20); these functions can be constructed by forward induction on the jumps of the driving Poisson processes. The processes $({\mathit{q}}_n,{\mathit{r}}_n)$ are defined by extending the above solutions to Ω, setting ${\mathit{q}}_n(\omega ,t)=0$ and ${\mathit{r}}_n(\omega ,t)=0$ for all $t\ge 0$ and all $\omega \notin {\mathrm{\Gamma }}_0$. In addition, we define

and we note that the sample paths of ${\mathcal{A}}_n,\hspace{0.17em}{\mathcal{D}}_n$ and ${\mathit{q}}_n$ lie in $D_{{\ell }_1}[0,\infty )$. Functional Equation (20) can now be rewritten as follows:for all $\omega \in {\mathrm{\Gamma }}_0$ and all $t\ge 0$.

5.2. Relative Compactness Results

Let $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ denote the space of càdlàg functions on $[0,\infty )$ with values in ${\mathbb{R}}^{\mathcal{I}}$. We endow the space ${\mathbb{R}}^{\mathcal{I}}$ with the metric defined in Online Appendix B, which is compatible with the product topology, and we equip $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ with the topology of uniform convergence over compact sets. The following proposition is proved in Online Appendix B.

Proposition 2.

Suppose that the load balancing policy is JLMU or SLTA. There exists a set of probability of one ${\mathrm{\Gamma }}_{\infty }$, where

$$\underset{n\to \infty }{\lim }{\Vert {\mathit{q}}_n(0)-q_0\Vert }_1=0,$$(22a)

$$\underset{n\to \infty }{\lim }\underset{t\in [0,T]}{\sup }{n}^{\gamma }\left|\frac{1}{n}{\mathcal{N}}_n^{\lambda }(t)-\lambda t\right|=0,$$(22b)

$$\underset{n\to \infty }{\lim }\underset{t\in [0,\mu jT]}{\sup }{n}^{\gamma }\left|\frac{1}{n}{\mathcal{N}}_{(i,j)}(nt)-t\right|=0\mathit{for}\mathit{all}(i,j)\in {\mathcal{I}}_{+},$$(22c)

for all $T\ge 0$ and $\gamma \in [0,1/2)$. Also, $\{{\mathcal{A}}_n(\omega ):n\ge 1\},\hspace{0.17em}\{{\mathcal{D}}_n(\omega ):n\ge 1\}$ and $\{{\mathit{q}}_n(\omega ):n\ge 1\}$ are relatively compact subsets of $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ for all $\omega \in {\mathrm{\Gamma }}_{\infty }$ and satisfy that the limit of every convergent subsequence is a function with locally Lipschitz coordinates. If the load balancing policy is SLTA, then there exists a random variable $R\ge 1$ such that, apart from the previous properties, we also have on ${\mathrm{\Gamma }}_{\infty }$ that

$${\mathit{r}}_n(0)\le R\mathit{and}{\mathit{q}}_n(0,i,j)<{\alpha }_n(i)\mathit{for}\mathit{all}(i,j)\mathrm{⊴}\left(i_{{\mathit{r}}_n(0)},j_{{\mathit{r}}_n(0)}\right)\mathit{and}n.$$(23)

The product topology of ${\mathbb{R}}^{\mathcal{I}}$ is coarser than the topology of ${\ell }_1$; thus, convergence in $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ does not imply convergence in $D_{{\ell }_1}[0,\infty )$. The following technical lemma is used to demonstrate that $\{{\mathcal{A}}_n:n\ge 1\},\hspace{0.17em}\{{\mathcal{D}}_n:n\ge 1\}$ and $\{{\mathit{q}}_n:n\ge 1\}$ are relatively compact in $D_{{\ell }_1}[0,\infty )$ with a probability of one; the proof is given in Online Appendix A.

Lemma 2.

Suppose that the load balancing policy is JLMU or SLTA. There exists a set of probability of one $\mathrm{\Gamma }\subset {\mathrm{\Gamma }}_{\infty }$ with the following property. For each $\omega \in \mathrm{\Gamma }$ and $T\ge 0$, there exist $j_T(\omega )$ and $n_T(\omega )$ such that

$${\mathcal{A}}_n(\omega ,t,i,j)=0\mathit{for}\mathit{all}t\in [0,T],1\le i\le m,j>j_T(\omega )\mathit{and}n\ge n_T(\omega ).$$

Also, if the load balancing policy is SLTA, then there exists $R_T(\omega )$ such that

$${\mathit{r}}_n(\omega ,t)\le R_T(\omega )\mathit{for}\mathit{all}t\in [0,T]\mathit{and}n\ge n_T(\omega ).$$

Proposition 3.

The sequences $\{{\mathcal{A}}_n(\omega ):n\ge 1\},\hspace{0.17em}\{{\mathcal{D}}_n(\omega ):n\ge 1\}$ and $\{{\mathit{q}}_n(\omega ):n\ge 1\}$ are relatively compact in $D_{{\ell }_1}[0,\infty )$ for all $\omega \in \mathrm{\Gamma }$.

Proof.

We fix some $\omega \in \mathrm{\Gamma }$, which we omit from the notation. For $\{{\mathcal{A}}_n:n\ge 1\}$, the claim is a straightforward consequence of Proposition 2 and Lemma 2. Here we prove the claim for $\{{\mathit{q}}_n:n\ge 1\}$. If the load balancing policy is JLMU, then (19) and (22a) imply that the claim also holds for $\{{\mathcal{D}}_n:n\ge 1\}$. If the load balancing policy is SLTA, then we must invoke (21a) instead of (19).

Consider any increasing sequence of natural numbers. By Proposition 2, there exists a subsequence $\mathcal{K}$ such that $\{{\mathit{q}}_k:k\in \mathcal{K}\}$ converges in $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ to a function $\mathit{q}\in D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ that satisfies $\mathit{q}(0)=q_0$ and has locally Lipschitz coordinates. Therefore, it suffices to prove that the latter limit in fact holds in $D_{{\ell }_1}[0,\infty )$. More specifically, we have to demonstrate that $\mathit{q}\in D_{{\ell }_1}[0,\infty )$ and that

$$\underset{k\to \infty }{\lim }\underset{t\in [0,T]}{\sup }{\Vert {\mathit{q}}_k(t)-\mathit{q}(t)\Vert }_1=0\text{for}\text{all}T\ge 0.$$

For this purpose, fix arbitrary $T\ge 0$ and $\epsilon >0$. In addition, let j_T and n_T be as in the statement of Lemma 2, which implies that ${\mathit{q}}_k(i,j)$ and $\mathit{q}(i,j)$ are nonincreasing on $[0,T]$ provided that $j>j_T$ and $k\ge n_T$. The coordinates of $\mathit{q}$ are continuous; thus, we may conclude from the monotone convergence theorem that ${\Vert \mathit{q}(s)-\mathit{q}(t)\Vert }_1\to 0$ as $s\to t\in [0,T]$ monotonically, from above or below. Because T is arbitrary, $\mathit{q}$ is continuous with respect to ${\Vert ·\Vert }_1$ and, in particular, $\mathit{q}\in D_{{\ell }_1}[0,\infty )$.

For all $t\in [0,T],\hspace{0.17em}l\ge j_T$ and $k\ge n_T$, we have

$$\begin{array}{ll}{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}\left|{\mathit{q}}_k(t,i,j)-\mathit{q}(t,i,j)\right|}}\hfill & \le {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}{\mathit{q}}_k}}(t,i,j)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}\mathit{q}}}(t,i,j)\hfill \\ \hfill & \le {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}{\mathit{q}}_k}}(0,i,j)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}{q}_0}}(i,j)\hfill \\ \hfill & \le {\Vert {\mathit{q}}_k(0)-q_0\Vert }_1+2{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>l}{q}_0}}(i,j).\hfill \end{array}$$

Because $q_0\in {\ell }_1$ and ${\Vert {\mathit{q}}_k(0)-q_0\Vert }_1\to 0$ with k, there exist $j_{\epsilon }\ge j_T$ and $k_0\ge n_T$ such that the following inequality holds for all $t\in [0,T]$ and $k\ge k_0$:

$${\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>j_{\epsilon }}\left|{\mathit{q}}_k(t,i,j)-\mathit{q}(t,i,j)\right|}}\le {\Vert {\mathit{q}}_k(0)-q_0\Vert }_1+2{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j>j_{\epsilon }}{q}_0}}(i,j)\le \frac{\epsilon }{2}.$$

Convergence in $D_{{\mathbb{R}}^{\mathcal{I}}}[0,\infty )$ implies uniform convergence over compact sets of the coordinate functions. In particular, there exists $k_{\epsilon }\ge k_0$ such that

$$\underset{t\in [0,T]}{\sup }{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^{j_{\epsilon }}\left|{\mathit{q}}_k(t,i,j)-\mathit{q}(t,i,j)\right|}}\le \frac{\epsilon }{2}\text{for}\text{all}k\ge k_{\epsilon }.$$

Therefore, we conclude that

$$\underset{t\in [0,T]}{\sup }{\Vert {\mathit{q}}_k(t)-\mathit{q}(t)\Vert }_1\le \epsilon \text{for}\text{all}k\ge k_{\epsilon },$$

which completes the proof because T and ε are arbitrary. □

6. Limiting Behavior of JLMU

In this section we assume that JLMU is used and we prove Theorems 3, 4, and 5. The first two theorems are proved in Section 6.1, which is devoted to the study of fluid trajectories. The proof of Theorem 5 is provided in Section 6.2.

6.1. Properties of Fluid Trajectories

To prove the uniqueness of fluid trajectories, we show that every fluid trajectory satisfies a system of equations involving one-dimensional Skorokhod reflection mappings. Then we use a Lipschitz property of these mappings to prove that the system of equations cannot have multiple solutions for a given initial condition. The proof strategy is inspired by a fluid limit derived in Bhamidi et al. (2022) using Skorokhod reflection mappings; this fluid limit corresponds to a system of parallel single-server queues with a JSQ policy.

Consider the space $D[0,\infty )$ of all real càdlàg functions defined on $[0,\infty )$ and let

The next lemma introduces the one-dimensional Skorokhod mappings with upper reflecting barrier; a proof is provided in Online Appendix A.

Lemma 3.

Fix $\alpha \in \mathbb{R}$ and suppose that $\mathit{x}\in D[0,\infty )$ is such that $\mathit{x}(0)\le \alpha$. Then there exist unique $\mathit{y},\mathit{z}\in D[0,\infty )$ such that the following statements hold.

• (a) We have $\mathit{z}(t)=\mathit{x}(t)-\mathit{y}(t)\le \alpha$ for all $t\ge 0$.

• (b) The function $\mathit{y}$ is nondecreasing, thus absolutely continuous, and $\mathit{y}(0)=0$.

• (c) The function $\mathit{y}$ is flat off $\{t\ge 0:\mathit{z}(t)=\alpha \}$, i.e., $\stackrel{.}{\mathit{y}}(t){\mathrm{𝟙}}_{\{\mathit{z}(t)<\alpha \}}=0$ almost everywhere.

The map $({\mathrm{\Psi }}_{\alpha },{\mathrm{\Phi }}_{\alpha })$ such that ${\mathrm{\Psi }}_{\alpha }(\mathit{x})=\mathit{y}$ and ${\mathrm{\Phi }}_{\alpha }(\mathit{x})=\mathit{z}$ is called the one-dimensional Skorokhod mapping with upper reflecting barrier at α and satisfies

$${\mathrm{\Psi }}_{\alpha }(\mathit{x})(t)=\underset{s\in [0,t]}{\sup }{\left[\mathit{x}(s)-\alpha \right]}^{+}\mathit{and}{\mathrm{\Phi }}_{\alpha }(\mathit{x})(t)=\mathit{x}(t)-{\mathrm{\Psi }}_{\alpha }(\mathit{x})(t).$$(24)

In addition, if $\mathit{x},\mathit{y}\in D[0,\infty )$ are any two functions such that $\mathit{x}(0),\mathit{y}(0)\le \alpha$, then for each $T\ge 0$ we have the following Lipschitz properties:

$${\Vert {\mathrm{\Psi }}_{\alpha }(\mathit{x})-{\mathrm{\Psi }}_{\alpha }(\mathit{y})\Vert }_T\le {\Vert \mathit{x}-\mathit{y}\Vert }_T\mathit{and}{\Vert {\mathrm{\Phi }}_{\alpha }(\mathit{x})-{\mathrm{\Phi }}_{\alpha }(\mathit{y})\Vert }_T\le 2{\Vert \mathit{x}-\mathit{y}\Vert }_T.$$

Consider càdlàg functions $\mathit{x}:[0,\infty )\to {\mathbb{R}}^{\mathcal{I}}$ such that $\mathit{x}(0,i,j)\le \alpha (i)$ for all $(i,j)\in {\mathcal{I}}_{+}$ and families of càdlàg functions ${\mathit{v}}_k:[0,\infty )\to \mathbb{R}$ such that ${\mathit{v}}_k(0)=0$ for all $k\in \mathbb{N}$. Define for each $k\ge 1$ a mapping Θ_k as follows:

here recall the enumeration of ${\mathcal{I}}_{+}$ introduced in Section 3.2. The next lemma establishes that $(\mathit{q},\mathit{w})$ satisfies the following set of conditions if $\mathit{q}$ is a fluid trajectory and $\mathit{w}$ is defined suitably in terms of $\mathit{q}$:

Lemma 4.

Let $\mathit{q}$ be a fluid trajectory and define $\mathit{r}$ such that

$$\sigma \left(\mathit{q}(t)\right)=\left(i_{\mathit{r}(t)},j_{\mathit{r}(t)}\right)\mathit{for}\mathit{all}t\ge 0.$$

Also, consider the absolutely continuous functions ${\mathit{w}}_k$ such that ${\mathit{w}}_k(0)=0$ and

$${\stackrel{.}{\mathit{w}}}_k(t)=\left[\lambda -{\displaystyle \sum _{s=1}^k\mu }{j}_s\left[\alpha (i_s)-\mathit{q}(t,i_s,j_s+1)\right]\right]{\mathrm{𝟙}}_{\{k<\mathit{r}(t)\}}\mathit{for}\mathit{all}k\ge 0.$$

Then $(\mathit{q},\mathit{w})$ satisfies (25).

Proof.

It is clear that $(\mathit{q},\mathit{w})$ satisfies (25c) and (25d), so we only need to verify that (25a) and (25b) hold as well. By Lemma 3, it is enough to check the following properties.

• (a) We have $\mathit{q}(t,i_k,j_k)={\mathrm{\Theta }}_k(\mathit{q},\mathit{w})(t)-{\mathit{w}}_k(t)\le \alpha (i_k)$ for all $t\ge 0$ and $k\ge 1$.

• (b) We have ${\mathit{w}}_k(0)=0$ and ${\stackrel{.}{\mathit{w}}}_k(t)\ge 0$ almost everywhere for all $k\ge 1$.

• (c) We have ${\stackrel{.}{\mathit{w}}}_k(t){\mathrm{𝟙}}_{\{\mathit{q}(t,i_k,j_k)<\alpha (i_k)\}}=0$ almost everywhere for all $k\ge 1$.

To establish (a), it suffices to show that

$$\stackrel{.}{\mathit{q}}(i_k,j_k)={\dot{\mathrm{\Theta }}}_k(\mathit{q},\mathit{w})-{\stackrel{.}{\mathit{w}}}_k={\stackrel{.}{\mathit{w}}}_{k-1}-{\stackrel{.}{\mathit{w}}}_k-\mu j_k\left[\mathit{q}(i_k,j_k)-\mathit{q}(i_k,j_k+1)\right]$$(26)

holds almost everywhere. Indeed, (a) holds at t = 0 and $\mathit{q}(i_k,j_k)\le \alpha (i_k)$ at all times. Because ${\stackrel{.}{\mathit{w}}}_{k-1}-{\stackrel{.}{\mathit{w}}}_k=\mathrm{\Lambda }(\mathit{q},i_k,j_k)$ and $\mathit{q}$ is a fluid trajectory, we obtain (26) from (12a).

Property (b) is a consequence of the definition of ${\mathit{w}}_k$ and (12b), and (c) follows from the following observation. If $\mathit{q}(t,i_k,j_k)<\alpha (i_k)$, then $(i_k,j_k)\mathrm{⊴}\sigma \left(\mathit{q}(t)\right)$ and thus $k\ge \mathit{r}(t)$, which implies that ${\stackrel{.}{\mathit{w}}}_k(t)=0$. □

Next we prove Theorem 3. As noted earlier, the proof relies on the Lipschitz property of the Skorokhod reflection mappings ${\mathrm{\Psi }}_{\alpha }$ and ${\mathrm{\Phi }}_{\alpha }$. In addition, a uniqueness of solutions result for certain Kolmogorov backward equations is used.

Proof of Theorem 3.

Suppose that there exist two fluid trajectories $\mathit{x}$ and $\mathit{y}$ such that $\mathit{x}(0)=\mathit{y}(0)=q$. Define $\mathit{v}$ in terms of $\mathit{x}$ and $\mathit{w}$ in terms of $\mathit{y}$ as in Lemma 4. It follows from the same lemma that $(\mathit{x},\mathit{v})$ and $(\mathit{y},\mathit{w})$ satisfy (25). Next we fix T > 0 and we prove that $\mathit{x}(t)=\mathit{y}(t)$ and $\mathit{v}(t)=\mathit{w}(t)$ for all $t\in [0,T]$.

As a first step, we demonstrate that there exists M > 0 such that

$${\mathit{v}}_k(t)={\mathit{w}}_k(t)=0\text{for}\text{all}t\in [0,T]\text{and}k\ge M.$$(27)

Because ${\stackrel{.}{\mathit{v}}}_{k-1}-{\stackrel{.}{\mathit{v}}}_k=\mathrm{\Lambda }(\mathit{x},i_k,j_k)\ge 0$, we conclude that ${\stackrel{.}{\mathit{v}}}_k\le {\stackrel{.}{\mathit{v}}}_{k-1}\le {\stackrel{.}{\mathit{v}}}_0=\lambda$ almost everywhere and for all $k\ge 1$. It follows from (12a), or equivalently from (25a) and (25b), that the following inequalities hold almost everywhere:

$$\stackrel{.}{\mathit{x}}(i_k,j_k)\le {\stackrel{.}{\mathit{v}}}_{k-1}-{\stackrel{.}{\mathit{v}}}_k\le {\stackrel{.}{\mathit{v}}}_{k-1}\le \lambda \text{for}\text{all}k\ge 1.$$(28)

Define ${\alpha }_{\min }≔\min \{\alpha (i):1\le i\le m\}$ and note that there exist $\epsilon >0$ and k₀ such that $q(i_k,j_k)\le \epsilon <{\alpha }_{\min }$ for all $k\ge k_0$ because $q\in Q\subset {\ell }_1$. Property (c) of Lemma 3 implies that ${\mathit{v}}_k$ is zero until $\mathit{x}(i_k,j_k)$ reaches $\alpha (i_k)\ge {\alpha }_{\min }$ for the first time. Using this remark and (28), it is possible to prove by induction on $k\ge k_0$ that

$${\mathit{v}}_k(t)=0\text{for}\text{all}t\in \left[0,(k+1-k_0)\frac{{\alpha }_{\min }-\epsilon }{\lambda }\right).$$

The same property holds if $(\mathit{x},\mathit{v})$ is replaced by $(\mathit{y},\mathit{w})$; thus, (27) holds.

Next we show that $\mathit{x}(t,i,j)=\mathit{y}(t,i,j)$ for all $t\in [0,T]$ and $(i,j)⊲(i_M,j_M)$. For this purpose, fix an arbitrary i and let $J_i≔\min \{j\ge 1:(i,j)⊲(i_M,j_M)\}$. By (27), both $\mathit{x}(i)$ and $\mathit{y}(i)$ satisfy the following initial value problem:

$$\stackrel{.}{\mathit{z}}(j)=-\mu j\left[\mathit{z}(j)-\mathit{z}(j+1)\right]\text{and}\mathit{z}(0,j)=q(i,j)\text{for}\text{all}j\ge J_i.$$(29)

The previous system of differential equations are the backward Kolmogorov equations of the pure birth process with state space $E_i≔\{j\in \mathbb{N}:j\ge J_i\}$ that has birth rate ${\lambda }_j≔\mu j$ at state j. This process is nonexplosive because ${\sum }_{j\ge J_i}1/{\lambda }_j=\infty$; hence, it follows from Fontes (1989) that the initial value problem (29) has a unique solution $\mathit{z}$ such that $\mathit{z}$ is bounded on $[0,t]\times E_i$ for all $t\ge 0$. Both $\mathit{x}(i)$ and $\mathit{y}(i)$ satisfy the latter condition because fluid trajectories take values in Q; thus, $\mathit{x}(t,i,j)=\mathit{y}(t,i,j)$ for all $t\ge 0$ and $j\ge J_i$.

We conclude by proving that $\mathit{x}(i_k,j_k)=\mathit{y}(i_k,j_k)$ and ${\mathit{v}}_k={\mathit{w}}_k$ along the interval $[0,T]$ for all $k\le M$. Let $\mathit{f}(i,j)≔\mathit{x}(i,j)-\mathit{y}(i,j)$ and ${\mathit{g}}_k≔{\mathit{v}}_k-{\mathit{w}}_k$. The subsequent arguments are analogous to those in (Bhamidi et al. 2022, section 4.1).

For all $t\in [0,T]$ and $k\le M$, we have

$$\begin{array}{ll}\hfill & {\Vert \mathit{f}(i_k,j_k)\Vert }_t\le 2{\Vert {\mathit{g}}_{k-1}\Vert }_t+2\mu j_k{\displaystyle {\int }_0^t{\Vert \mathit{f}(i_k,j_k)\Vert }_s}ds+2\mu j_k{\displaystyle {\int }_0^t{\Vert \mathit{f}(i_k,j_k+1)\Vert }_s}ds,\hfill \\ \hfill & {\Vert {\mathit{g}}_k\Vert }_t\le {\Vert {\mathit{g}}_{k-1}\Vert }_t+\mu j_k{\displaystyle {\int }_0^t{\Vert \mathit{f}(i_k,j_k)\Vert }_s}ds+\mu j_k{\displaystyle {\int }_0^t{\Vert \mathit{f}(i_k,j_k+1)\Vert }_s}ds;\hfill \end{array}$$

these inequalities follow from the Lipschitz properties of ${\mathrm{\Psi }}_{\alpha (i_k)}$ and ${\mathrm{\Phi }}_{\alpha (i_k)}$. Let us define $J≔\max \{j_k:k\le M\}$ and $\mathit{h}(t)≔\max \{{\Vert \mathit{f}(i_k,j_k)\Vert }_t:k\le M\}$, then

$${\Vert {\mathit{g}}_k\Vert }_t\le {\Vert {\mathit{g}}_{k-1}\Vert }_t+2\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds\le 2k\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds$$

for all $t\in [0,T]$ and $k\le M$. For the first inequality, $\mathit{f}(i_k,j_k+1)$ is identically zero along the interval $[0,T]$ if $(i_k,j_k+1)=(i_l,j_l)$ for some l > M, and for the last inequality, observe that ${\mathit{g}}_0$ is identically zero by (25c). In addition, we have

$${\Vert \mathit{f}(i_k,j_k)\Vert }_t\le 4k\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds+4\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds\le 4(M+1)\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds$$

for all $t\in [0,T]$ and $k\le M$. We conclude that

$$\mathit{h}(t)\le 4(M+1)\mu J{\displaystyle {\int }_0^t\mathit{h}}(s)ds\text{for}\text{all}t\in [0,T].$$

Therefore, Grönwall’s inequality yields $\mathit{h}(t)=0$ for all $t\in [0,T]$, and this in turn implies that $(\mathit{x},\mathit{v})=(\mathit{y},\mathit{w})$ along the interval $[0,T]$. □

We conclude this section by establishing that (12) has a unique equilibrium point and that all fluid trajectories converge to this equilibrium point over time.

Proof of Theorem 4.

First, we verify that $q_{*}$ is an equilibrium of (12). To this end, $\sigma \left(q_{*}\right)={\sigma }_{*}$, and the right-hand side of (12a) equals zero if $(i,j)\ne {\sigma }_{*}$ and $\mathit{q}=q_{*}$. It only remains to be shown that this also holds for $(i_{*},j_{*})≔{\sigma }_{*}$.

If $(i,j)=\sigma (q)$, then

$$\mathrm{\Lambda }(q,i,j)-\mu j\left[q(i,j)-q(i,j+1)\right]=\lambda -{\displaystyle \sum _{(r,s)\mathrm{⊵}(i,j)}\mu }s\left[q(r,s)-q(r,s+1)\right].$$

Define $J^i≔\max \{j\ge 1:(i,j)\mathrm{⊵}{\sigma }_{*}\}$ for each i. If $(i,j)=(i_{*},j_{*})$ and $\mathit{q}$ is replaced by $q_{*}$, then the right-hand side of (12a) equals

$$\begin{array}{ll}\lambda -{\displaystyle \sum _{(i,j)\mathrm{⊵}{\sigma }_{*}}\mu }j\left[q_{*}(i,j)-q_{*}(i,j+1)\right]\hfill & =\lambda -{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^{J^i}\mu }}j\left[q_{*}(i,j)-q_{*}(i,j+1)\right]\hfill \\ \hfill & =\lambda -{\displaystyle \sum _{i\ne i_{*}}\mu }{J}^i\alpha (i)-\mu (j_{*}-1)\alpha (i_{*})-\mu q_{*}(i_{*},j_{*})\hfill \\ \hfill & ={\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\mu }\alpha (i)-{\displaystyle \sum _{i\ne i_{*}}\mu }{J}^i\alpha (i)-\mu (j_{*}-1)\alpha (i_{*}).\hfill \end{array}$$

The expression in the last line equals zero by definition of Jⁱ, and therefore we conclude that $q_{*}$ is indeed an equilibrium point of (12).

Next, we prove that all fluid trajectories converge coordinatewise to $q_{*}$ over time; this implies, in particular, that $q_{*}$ is the unique equilibrium of (12). Afterward we prove that all fluid trajectories in fact converge to $q_{*}$ in ${\ell }_1$.

If $\mathit{q}$ is a fluid trajectory and $T\ge 0$, then there exists j_T such that $\mathrm{\Lambda }(\mathit{q},i,j)=0$ in $[0,T]$ for all i and $j\ge j_T$. This can be established directly from (12), but also using Lemma 2, Theorem 5, and the uniqueness of solutions. For each i and $k\ge j_T$,

$${\displaystyle \sum _{j=j_T}^k\stackrel{.}{\mathit{q}}}(i,j)=-{\displaystyle \sum _{j=j_T}^k\mu }j\left[\mathit{q}(i,j)-\mathit{q}(i,j+1)\right]=-\mu (j_T-1)\mathit{q}(i,j_T)-{\displaystyle \sum _{j=j_T}^k\mathit{q}}(i,j)+\mu k\mathit{q}(i,k+1).$$

The last term vanishes as $k\to \infty$ because $\mathit{q}$ takes values in ${\ell }_1$. Also, $\mathit{q}(i,j)$ does not increase in $[0,T]$ for all $j\ge j_T$ because $\mathrm{\Lambda }(\mathit{q},i,j)=0$ along $[0,T]$. Therefore,

$$\begin{array}{ll}\left|{\displaystyle \sum _{j=j_T}^k\stackrel{.}{\mathit{q}}}(t,i,j)-{\displaystyle \sum _{j=j_T}^{\infty }\stackrel{.}{\mathit{q}}}(t,i,j)\right|\hfill & ={\displaystyle \sum _{j=k+1}^{\infty }\mathit{q}}(t,i,j)+\mu k\mathit{q}(t,i,k+1)\hfill \\ \hfill & \le {\displaystyle \sum _{j=k+1}^{\infty }\mathit{q}}(0,i,j)+\mu k\mathit{q}(0,i,k+1)\text{for}\text{all}t\in [0,T].\hfill \end{array}$$

The right-hand side vanishes as $k\to \infty$ because $\mathit{q}(0)\in {\ell }_1$; thus, the left-hand side converges uniformly to zero over $[0,T]$ and, by Rudin (1976, theorem 7.17), the derivative of

$${\displaystyle \sum _{j=j_T}^{\infty }\mathit{q}}(t,i,j)\text{is}{\displaystyle \sum _{j=j_T}^{\infty }\stackrel{.}{\mathit{q}}}(t,i,j)\text{for}\text{all}i\text{and}t\in (0,T).$$

This allows for the interchanges of summation and differentation that appear later.

Consider the function

$$\mathit{f}≔{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathit{q}}(i,j).$$

It follows from (12a) that $\stackrel{.}{\mathit{f}}=\lambda -\mu \mathit{f}$. Thus, $\mathit{f}(t)=\rho +\left[\mathit{f}(0)-\rho \right]{\text{e}}^{-\mu t}$ for all $t\ge 0$. We conclude from this identity and (10) that

$$t_0≔\min \left\{t\ge 0:\mathit{f}(t)\le {\displaystyle \sum _{(i,j)\mathrm{⊵}{\sigma }_{*}}\alpha }(i)\right\}$$

exists and is finite. If $t>t_0$, then the inequality inside the minimum sign is strict, and this implies that $\sigma \left(\mathit{q}(t)\right)\mathrm{⊵}{\sigma }_{*}$ because $\mathit{q}(t,i,j)=\alpha (i)$ for all $(i,j)⊳\sigma \left(\mathit{q}(t)\right)$.

Consider the function

$$\mathit{g}≔{\displaystyle \sum _{(i,j)\mathrm{⊴}{\sigma }_{*}}\mathit{q}}(i,j).$$

As noted previously, if $t>t_0$, then $\sigma \left(\mathit{q}(t)\right)\mathrm{⊵}{\sigma }_{*}$, and thus we have

$$\begin{array}{ll}\dot{\mathit{g}}(t)\hfill & =\left[\lambda -{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\mu }j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\right]{\mathrm{𝟙}}_{\{\sigma (\mathit{q}(t))={\sigma }_{*}\}}\hfill \\ \hfill & -{\displaystyle \sum _{(i,j)\mathrm{⊴}{\sigma }_{*}}\mu }j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right],\hfill \end{array}$$

because $\mathrm{\Lambda }(\mathit{q}(t),i,j)=0$ for all $(i,j)\mathrm{⊴}{\sigma }_{*}$ except perhaps for $(i,j)={\sigma }_{*}$. Hence,

$$\begin{array}{ll}\dot{\mathit{g}}(t)\hfill & =\left[\lambda -{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mu }j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\right]{\mathrm{𝟙}}_{\{\sigma (\mathit{q}(t))={\sigma }_{*}\}}\hfill \\ \hfill & -\left[{\displaystyle \sum _{(i,j)\mathrm{⊴}{\sigma }_{*}}\mu }j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\right]{\mathrm{𝟙}}_{\{\sigma (\mathit{q}(t))\ne {\sigma }_{*}\}}\hfill \\ \hfill & =\left[\lambda -{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mu }\mathit{q}(t,i,j)\right]{\mathrm{𝟙}}_{\{\sigma (\mathit{q}(t))={\sigma }_{*}\}}\hfill \\ \hfill & -\left[{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=J_i}^{\infty }\mu }}j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\right]{\mathrm{𝟙}}_{\{\sigma (\mathit{q}(t))\ne {\sigma }_{*}\}},\hfill \end{array}$$

where $J_i≔\min \{j\ge 1:(i,j)\mathrm{⊴}{\sigma }_{*}\}$. Note that $\mathit{q}(t,i,j)=\alpha (i)$ if $(i,j)⊳\sigma \left(\mathit{q}(t)\right)$. Thus,

$$\begin{array}{ll}\dot{\mathit{g}}(t)\hfill & =\left[\lambda -\mu {\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i)-\mu \mathit{g}(t)\right]{\mathrm{𝟙}}_{\{\sigma \left(\mathit{q}(t)\right)={\sigma }_{*}\}}\hfill \\ \hfill & -\left[{\displaystyle \sum _{i=1}^m\left(\mu (J_i-1)\mathit{q}(t,i,J_i)+{\displaystyle \sum _{j=J_i}^{\infty }\mu }\mathit{q}(t,i,j)\right)}\right]{\mathrm{𝟙}}_{\{\sigma \left(\mathit{q}(t)\right)\ne {\sigma }_{*}\}}\hfill \\ \hfill & =\left[\lambda -\mu {\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i)-\mu \mathit{g}(t)\right]{\mathrm{𝟙}}_{\{\sigma \left(\mathit{q}(t)\right)={\sigma }_{*}\}}\hfill \\ \hfill & -\mu \mathit{g}(t){\mathrm{𝟙}}_{\{\sigma \left(\mathit{q}(t)\right)\ne {\sigma }_{*}\}}-\mu {\displaystyle \sum _{i=1}^m(J_i-1)}\mathit{q}(t,i,J_i){\mathrm{𝟙}}_{\{\sigma \left(\mathit{q}(t)\right)\ne {\sigma }_{*}\}}\hfill \\ \hfill & \le {\left[\lambda -\mu {\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i)\right]}^{+}-\mu \mathit{g}(t).\hfill \end{array}$$

By definition of ${\sigma }_{*}$, we have $\theta ≔{\sum }_{(i,j)⊳{\sigma }_{*}}\alpha (i)\le \rho$, and from the above bound for $\dot{\mathit{g}}$, we get

$$\mathit{g}(t)\le \rho -\theta +\left[\mathit{g}(t_0)-(\rho -\theta )\right]{\text{e}}^{-\mu (t-t_0)}\text{for}\text{all}t\ge t_0.$$

We conclude that

$$\underset{t\to \infty }{\lim \inf }{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\mathit{q}}(t,i,j)=\underset{t\to \infty }{\lim \inf }\left[\mathit{f}(t)-\mathit{g}(t)\right]\ge \rho -(\rho -\theta )={\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i).$$

This proves that $\mathit{q}(t,i,j)\to \alpha (i)=q_{*}(i,j)$ over time for all $(i,j)⊳{\sigma }_{*}$.

Consider now the function

$$\mathit{h}≔{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\mathit{q}}(i,j).$$

Recall that $\sigma \left(\mathit{q}(t)\right)\mathrm{⊵}{\sigma }_{*}$ for all $t>t_0$. Therefore,

$$\dot{\mathit{h}}(t)=-{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\mu }j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\le -\mu \mathit{h}(t)\text{for}\text{all}t>t_0.$$

It follows that $\mathit{h}(t)\le \mathit{h}(t_0){\text{e}}^{-\mu (t-t_0)}$ for all $t\ge t_0$, and thus $\mathit{q}(t,i,j)\to 0=q_{*}(i,j)$ over time for all $(i,j)⊲{\sigma }_{*}$. Furthermore, we have

$$\underset{t\to \infty }{\lim }\mathit{q}(t,{\sigma }_{*})=\underset{t\to \infty }{\lim }\left[\mathit{f}(t)-\mathit{h}(t)-{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\mathit{q}}(t,i,j)\right]=\rho -{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i)=q_{*}({\sigma }_{*}),$$

and thus $\mathit{q}(t,i,j)\to q_{*}(i,j)$ for all $(i,j)\in \mathcal{I}$.

Finally, observe that

$$\underset{t\to \infty }{\lim }{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\left|\mathit{q}(t,i,j)-q_{*}(i,j)\right|}=\underset{t\to \infty }{\lim }\mathit{h}(t)=0.$$

Consequently, $\mathit{q}(t)\to q_{*}$ over time not only coordinatewise but also in ${\ell }_1$. □

6.2. Proof of the Fluid Limit

To prove Theorem 5, it suffices to demonstrate, for each $\omega \in \mathrm{\Gamma }$, that every subsequence of $\{{\mathit{q}}_n(\omega ):n\ge 1\}$ has a further subsequence with a limit in $D_{{\ell }_1}[0,\infty )$ and that this limit is the unique fluid trajectory starting at $q_0(\omega )$. The first part is covered by Proposition 3, every subsequence of $\{{\mathit{q}}_n(\omega ):n\ge 1\}$ has a further subsequence with a limit in $D_{{\ell }_1}[0,\infty )$. Next, we characterize the limits of the convergent subsequences.

Let us fix an arbitrary $\omega \in \mathrm{\Gamma }$, which we omit from the notation for brevity, and an increasing sequence $\mathcal{K}\subset \mathbb{N}$ such that $\{{\mathcal{A}}_k:k\in \mathcal{K}\},\hspace{0.17em}\{{\mathcal{D}}_k:k\in \mathcal{K}\}$ and $\{q_k:k\in \mathcal{K}\}$ converge in $D_{{\ell }_1}[0,\infty )$ to certain functions $\mathit{a},\hspace{0.17em}\mathit{d}$ and $\mathit{q}$, respectively, which have locally Lipschitz coordinates by Proposition 2. To characterize these three limits, it suffices to just characterize $\mathit{a}$ and $\mathit{d}$ because (19) and (22a) imply that

Because $\mathit{a}$ and $\mathit{d}$ have locally Lipschitz coordinates, there exists $\mathcal{R}\subset (0,\infty )$ such that ${\mathcal{R}}^c$ has zero Lebesgue measure and the derivatives of $\mathit{a}(i,j)$ and $\mathit{d}(i,j)$ exist for all $(i,j)\in \mathcal{I}$ at all points in $\mathcal{R}$. These derivatives are zero if j = 0 by the definitions of ${\mathcal{A}}_k$ and ${\mathcal{D}}_k$. The following lemma computes the derivatives for $(i,j)\in {\mathcal{I}}_{+}$.

Lemma 5.

Fix an arbitrary $t_0\in \mathcal{R}$, we have

$$\dot{\mathit{d}}(t_0,i,j)=\mu j\left[\mathit{q}(t_0,i,j)-\mathit{q}(t_0,i,j+1)\right]\mathit{for}\mathit{all}(i,j)\in {\mathcal{I}}_{+}.$$(31)

Furthermore, $\mathit{q}(t_0)$ and the derivatives $\dot{\mathit{a}}(t_0,i,j)$ satisfy

$$\mathrm{\Lambda }(\mathit{q}(t_0),i,j)=\dot{\mathit{a}}(t_0,i,j)\ge 0\mathit{for}\mathit{all}(i,j)\in {\mathcal{I}}_{+}.$$(32)

Proof.

The sequences $\{{\mathcal{D}}_k(i,j):k\in \mathcal{K}\}$ and $\{{\mathit{q}}_k(i,j):k\in \mathcal{K}\}$ converge uniformly over compact sets to $\mathit{d}(i,j)$ and $\mathit{q}(i,j)$, respectively, for all $(i,j)\in \mathcal{I}$. This remark, the definition of ${\mathcal{D}}_k$ and (22c) imply that

$$\mathit{d}(t,i,j)={\displaystyle {\int }_0^t\mu }j\left[\mathit{q}(s,i,j)-\mathit{q}(s,i,j+1)\right]ds\text{for}\text{all}t\ge 0.$$

It is clear that this identity establishes (31).

We now prove that

$${\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\dot{\mathit{a}}}(t_0,i,j)=\lambda \text{and}\dot{\mathit{a}}(t_0,i,j)\ge 0\text{for}\text{all}(i,j)\in {\mathcal{I}}_{+}.$$(33)

The derivatives $\dot{\mathit{a}}(t_0,i,j)$ are nonnegative because the processes ${\mathcal{A}}_k(i,j)$ are nondecreasing, so we only need to show that the derivatives add up to λ. For this purpose, note that

$${\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}{\mathcal{A}}_k}(t,i,j)={\mathcal{N}}_k^{\lambda }(t)\text{for}\text{all}t\ge 0\text{and}k\in \mathcal{K}.$$

Fix $T>t_0$, and let j_T and n_T be as in Lemma 2. The left-hand side has at most $m{j}_T$ nonzero terms for all $k\ge n_T$ and $t\in [0,T]$. It follows from (22b) that

$${\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathit{a}}(t,i,j)=\lambda t\text{for}\text{all}t\in [0,T].$$

This yields (33) because the left-hand side has at most mj_T nonzero terms.

It follows from (30) that the derivative of $\mathit{q}(i,j)$ exists at t₀ and

$$\stackrel{.}{\mathit{q}}(t_0,i,j)=\dot{\mathit{a}}(t_0,i,j)-\dot{\mathit{d}}(t_0,i,j)\text{for}\text{all}(i,j)\in \mathcal{I}.$$

Note that $\mathit{q}(i,j)$ is upper bounded by $\alpha (i)$, so $\mathit{q}(t_0,i,j)=\alpha (i)$ implies $\stackrel{.}{\mathit{q}}(t_0,i,j)=0$. Thus,

$$\dot{\mathit{a}}(t_0,i,j)=\dot{\mathit{d}}(t_0,i,j)=\mu j[\alpha (i)-\mathit{q}(t_0,i,j+1)]\text{if}\mathit{q}(t_0,i,j)=\alpha (i).$$(34)

To prove the equality in (32), define ${\sigma }_0=(i_0,j_0)≔\sigma \left(\mathit{q}(t_0)\right)$. If $(i,j)⊳{\sigma }_0$, then $\mathit{q}(t_0,i,j)=\mathit{q}(t_0,i,j-1)$ by (8). Moreover, for a fixed i, the marginal utility $\mathrm{\Delta }(i,j)$ does not increase with j because u_i is a concave function. Therefore, $(i,j)⊳{\sigma }_0$ and j > 1 imply that $(i,j-1)⊳{\sigma }_0$. We conclude that

$$\mathit{q}(t_0,i,j)=\mathit{q}(t_0,i,0)=\alpha (i)\text{for}\text{all}(i,j)⊳{\sigma }_0.$$

The last property and (34) imply that (32) holds if $(i,j)⊳{\sigma }_0$.

Note that $\mathit{q}(t_0,{\sigma }_0)<\alpha (i_0)$ by (8). Because $\mathit{q}({\sigma }_0)$ is continuous and ${\mathit{q}}_k({\sigma }_0)$ converges uniformly over compact sets to $\mathit{q}({\sigma }_0)$, there exist $\epsilon >0$ and $k_{\epsilon }\in \mathcal{K}$ such that

$${\mathit{q}}_k(t,{\sigma }_0)<{\alpha }_k(i_0)\text{for}\text{all}t\in (t_0-\epsilon ,t_0+\epsilon )\text{and}k\ge k_{\epsilon }.$$

It follows from (8) and the last statement that

$${\sigma }_0\mathrm{⊴}\sigma \left({\mathit{q}}_k(t)\right)\text{for}\text{all}t\in (t_0-\epsilon ,t_0+\epsilon )\text{and}k\ge k_{\epsilon }.$$

Thus, ${\mathcal{A}}_k(i,j)$ is constant over $(t_0-\epsilon ,t_0+\epsilon )$ if $k\ge k_{\epsilon }$ and $(i,j)⊲{\sigma }_0$. Indeed, server pools of class i with exactly $j-1$ tasks are not assigned incoming tasks in the system with k server pools if $(i,j)⊲\sigma \left({\mathit{q}}_k\right)$. This proves (32) for $(i,j)⊲{\sigma }_0$, and we conclude from (33) that (32) must also hold in the case $(i,j)={\sigma }_0$. □

Here we complete the proof of Theorem 5.

Proof of Theorem 5.

As previously, we fix some $\omega \in \mathrm{\Gamma }$ that we omit from the notation. Every subsequence of $\{{\mathit{q}}_n:n\ge 1\}$ has a further subsequence that converges in $D_{{\ell }_1}[0,\infty )$ by Proposition 3. It follows from (22a) and Lemma 5 that the limit $\mathit{q}$ of this convergent subsequence is a fluid trajectory with $\mathit{q}(0)=q_0$, which determines $\mathit{q}$ by Theorem 3. □

7. Limiting Behavior of SLTA

In this section, we assume that the load balancing policy is SLTA, and we leverage a methodology developed in Goldsztajn et al. (2021a) to prove Theorem 6. The first steps of the proof are carried out in Section 7.1, where we establish that certain dynamical properties of the system hold asymptotically with probability one. The proof is completed in Section 7.2, where we analyze the evolution of the learning scheme over time. Although the arguments used here are more involved due to the heterogeneity of the system, most of the proofs are conceptually similar to those in Goldsztajn et al. (2021a) and hence are deferred to Online Appendix C.

7.1. Asymptotic Dynamical Properties

In this section, we establish asymptotic dynamical properties pertaining to the total and tail mass processes, which are defined as

respectively. Recall that the total mass process was introduced in (1) and represents the total number of tasks in the system, normalized by the number of server pools. The following proposition is proved in Online Appendix C.

Proposition 4.

For each $\omega \in \mathrm{\Gamma }$, the sequence $\{{\mathit{s}}_n(\omega ):n\ge 1\}$ converges uniformly over compact sets to the unique function $\mathit{s}(\omega )$ such that

$$\dot{\mathit{s}}(\omega )=\lambda -\mu \mathit{s}(\omega )\mathit{and}\mathit{s}(\omega ,0)=s_0(\omega ),$$(35)

where s₀ is as defined in (13). Explicitly, $\mathit{s}(\omega ,t)=\rho +\left[s_0(\omega )-\rho \right]{\text{e}}^{-\mu t}$.

Although the previous law of large numbers is known to hold weakly, it is not straightforward that it holds with probability one under the coupled construction of sample paths adopted in Section 5.1; this fact is established in Proposition 4.

The next result is also proved in Online Appendix C, and it provides an asymptotic upper bound for certain tail mass processes, under specific conditions concerning ${\mathit{q}}_n$ and ${\mathit{r}}_n$. The upper bound implies at least an exponentially fast decay over time.

Proposition 5.

Suppose that the next conditions hold for a given $\omega \in \mathrm{\Gamma }$ and a given increasing sequence $\mathcal{K}$ of natural numbers.

• (a) The sequence $\{{\mathit{q}}_k(\omega ):k\in \mathcal{K}\}$ converges in $D_{{\ell }_1}[0,\infty )$ to some function $\mathit{q}$.

• (b) There exist r > 1 and $0\le t_0<t_1$ such that

  $${\mathit{r}}_k(\omega ,t)\le r\mathit{and}{\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\mathit{q}}_k}(\omega ,t,i,j)<{\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\alpha }_k}(i)$$
  for all $t\in [t_0,t_1]$ and $k\in \mathcal{K}$.

Then $\mathit{q}(i,j)$ is differentiable on (t₀, t₁) for all $(i,j)\mathrm{⊴}(i_r,j_r)$ and satisfies

$$\stackrel{.}{\mathit{q}}(t,i,j)=-\mu j\left[\mathit{q}(t,i,j)-\mathit{q}(t,i,j+1)\right]\mathit{for}\mathit{all}t\in (t_0,t_1).$$

Furthermore, the sequence of tail mass processes $\{{\mathit{v}}_k(\omega ,r):k\in \mathcal{K}\}$ converges uniformly over compact sets to a function $\mathit{v}(\omega ,r)$ that satisfies

$$\mathit{v}(\omega ,t,r)<\mathit{s}(\omega ,t_0){\text{e}}^{-\mu \left(t-t_0\right)}\mathit{for}\mathit{all}t\in [t_0,t_1].$$

7.2. Evolution of the Learning Scheme

In this section, we complete the proof of Theorem 6. In Section 7.2.1, we establish that there exists a neighborhood of zero outside of which ${\mathit{r}}_n$ is asymptotically upper bounded by $r_{*}$ with a probability of one. This property partially proves (16a) and is used to obtain (16c). The proof of (16a) is finished in Section 7.2.2, where we also establish (16b).

7.2.1. Preliminary Results.

The following proposition states that ${\mathit{r}}_n$ is asymptotically upper bounded by $r_{*}$ outside of a neighborhood of zero with probability one; the proof is deferred to Online Appendix C.

Proposition 6.

There exists a function ${\tau }_{\text{bd}}:[0,\infty )\to \mathbb{R}$ with the following property. If $\omega \in \mathrm{\Gamma }$ and $T\ge \tau >{\tau }_{\text{bd}}(s_0(\omega ))$, then there exists $n_{\text{bd}}^{\tau ,T}(\omega )$ such that

$${\mathit{r}}_n(\omega ,t)\le r_{*}\mathit{for}\mathit{all}t\in [\tau ,T]\mathit{and}n\ge n_{\text{bd}}^{\tau ,T}(\omega ).$$

The following corollary establishes (16c).

Corollary 1.

For each $\omega \in \mathrm{\Gamma }$ and $T\ge \tau >{\tau }_{\text{bd}}(s_0(\omega ))$, we have

$$\underset{n\to \infty }{\lim \sup }\underset{t\in [\tau ,T]}{\sup }{\mathit{v}}_n(\omega ,t,r_{*}+1){\text{e}}^{\mu (t-\tau )}\le \mathit{s}(\omega ,\tau ).$$

In particular, (16c) holds.

Proof.

We fix $\omega \in \mathrm{\Gamma }$ and $T\ge \tau >{\tau }_{\text{bd}}(s_0(\omega ))$, and we omit ω from the notation for brevity. Suppose that the statement of the corollary does not hold; then, there exist $\epsilon >0$ and an increasing sequence $\mathcal{K}$ of natural numbers such that

$$\underset{t\in [\tau ,T]}{\sup }{\mathit{v}}_k(t,r_{*}+1){\text{e}}^{\mu (t-\tau )}>\mathit{s}(\tau )+\epsilon \text{for}\text{all}k\in \mathcal{K}.$$

By Propositions 3 and 6, we may assume that $\{{\mathit{q}}_k:k\in \mathcal{K}\}$ has a limit in $D_{{\ell }_1}[0,\infty )$ and that ${\mathit{r}}_k(t)\le r_{*}$ for all $t\in [\tau ,T]$ and all $k\in \mathcal{K}$. The latter property implies that

$${\displaystyle \sum _{(i,j)\mathrm{⊵}(i_{r_{*}},j_{r_{*}})}{\mathit{q}}_k}(t,i,j)<{\displaystyle \sum _{(i,j)\mathrm{⊵}(i_{r_{*}},j_{r_{*}})}{\alpha }_k}(i)\text{for}\text{all}t\in [\tau ,T],$$

because otherwise the number of tokens would be zero, which cannot occur by Remark 1. Therefore, Proposition 5 holds with $r=r_{*}+1$ along the interval $[\tau ,T]$, and in particular, there exists a function $\mathit{v}(r_{*}+1)$ such that

$$\begin{array}{ll}\hfill & \mathit{v}(t,r_{*}+1)<\mathit{s}(\tau ){\text{e}}^{-\mu (t-\tau )}\text{for}\text{all}t\in [\tau ,T],\hfill \\ \hfill & \underset{k\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }|{\mathit{v}}_k(t,r_{*}+1)-\mathit{v}(t,r_{*}+1)|=0.\hfill \end{array}$$

This leads to a contradiction, so the statement of the corollary must hold. □

7.2.2. Proof of Theorem 6.

Here we complete the proof of Theorem 6. For this purpose, let

for all $t\ge 0$ and $r\ge 1$. It follows from (22b) and (22c) that

The following two technical lemmas are proved in Online Appendix C.

Lemma 6.

Fix $\omega \in \mathrm{\Gamma },\hspace{0.17em}T\ge 0$ and r > 1. Suppose that there exist an increasing sequence $\mathcal{K}$ of natural numbers and random times $0\le {\tau }_{k,1}\le {\tau }_{k,2}\le T$ such that

$${\mathit{r}}_k(\omega ,t)\le r\mathit{and}{\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\mathit{q}}_k}(\omega ,t,i,j)<{\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\alpha }_k}(i)$$

for all $t\in \left[{\tau }_{k,1}(\omega ),{\tau }_{k,2}(\omega )\right)$ and $k\in \mathcal{K}$. Then

$$\begin{array}{ll}\hfill & {\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\mathit{q}}_k}(\omega ,t,i,j)-{\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\mathit{q}}_k}\left(\omega ,{\tau }_{k,1}(\omega ),i,j\right)\ge \left[t-{\tau }_{k,1}(\omega )\right]\left[\lambda -\mu {\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\alpha }_k}(i)\right]-2\underset{s\in [0,T]}{\sup }|{\delta }_k(\omega ,s,r)|\hfill \end{array}$$

for all $t\in \left[{\tau }_{k,1}(\omega ),{\tau }_{k,2}(\omega )\right]$ and $k\in \mathcal{K}$.

Lemma 7.

Fix $\omega \in \mathrm{\Gamma },\hspace{0.17em}T\ge 0$ and $1\le r\le r_{*}$. Assume that there exist an increasing sequence $\mathcal{K}$ of natural numbers and random times $0\le {\zeta }_{k,1}\le {\zeta }_{k,2}\le T$ such that

$${\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\mathit{q}}_k}\left(\omega ,{\zeta }_{k,1}(\omega ),i,j\right)={\displaystyle \sum _{(i,j)⊳(i_r,j_r)}{\alpha }_k}(i)\mathit{and}{\mathit{r}}_k(\omega ,t)\le r$$

for all $t\in \left[{\zeta }_{k,1}(\omega ),{\zeta }_{k,2}(\omega )\right]$ and $k\in \mathcal{K}$. For each $\gamma \in [0,1/2)$, we have

$${\mathit{r}}_k(\omega ,t)=r\mathit{and}{\mathit{q}}_k(\omega ,t,i,j)\ge {\alpha }_k(i)-k^{-\gamma }$$

for all $(i,j)⊳(i_r,j_r),\hspace{0.17em}t\in \left[{\zeta }_{k,1}(\omega ),{\zeta }_{k,2}(\omega )\right]$ and all large enough $k\in \mathcal{K}$.

These lemmas are used to complete the proof of Theorem 6.

Proof of Theorem 6.

We define ${\tau }_{\text{eq}}$ as follows:

$${\tau }_{\text{eq}}(s)≔{\tau }_{\text{bd}}(s)+\frac{1}{\mu }\log \left(\frac{\rho }{\rho -\sum _{(i,j)⊳{\sigma }_{*}}\alpha (i)}\right)\text{for}\text{all}s\ge 0.$$

Fix $\omega \in \mathrm{\Gamma }$ and $T\ge \tau >{\tau }_{\text{eq}}(s_0(\omega ))$ as in the statement of the theorem; we omit ω from the notation for brevity. Given $0<\epsilon <\rho -{\sum }_{(i,j)⊳{\sigma }_{*}}\alpha (i)$, we define

$$\tau (\epsilon )≔\min \left\{t\ge 0:\rho \left(1-{\text{e}}^{-\mu t}\right)-\epsilon \ge {\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}\alpha }(i)\right\}=\frac{1}{\mu }\log \left(\frac{\rho }{\rho -\sum _{(i,j)⊳{\sigma }_{*}}\alpha (i)-\epsilon }\right).$$

Fix ${\tau }_0>{\tau }_{\text{bd}}(s_0)$ and ε such that $\tau ={\tau }_0+\tau (\epsilon )+\epsilon$. This is possible because ${\tau }_{\text{bd}}(s_0)+\tau (\epsilon )+\epsilon$ decreases to ${\tau }_{\text{eq}}(s_0)$ as $\epsilon \to 0$. In addition, consider the random times

$${\xi }_n≔\inf \left\{t\ge {\tau }_0:{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}{\mathit{q}}_n}(t,i,j)={\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}{\alpha }_n}(i)\right\}.$$

The proofs of (16a) and (16b) will be completed if we demonstrate that ${\xi }_n\le \tau$ for all large enough n. Indeed, if this is established, then (16a) and (16b) follow from Lemma 7 with $r≔r_{*},\hspace{0.17em}{\zeta }_{n,1}≔{\xi }_n$ and ${\zeta }_{n,2}≔T$. The hypotheses of the lemma hold because ${\mathit{r}}_n(t)\le r_{*}$ for all $t\in [{\tau }_0,T]$ and all large enough n by the choice of τ₀ and Proposition 6.

To prove that ${\xi }_n\le \tau$ for all large enough n, we show that

$$\underset{n\to \infty }{\lim \sup }{\xi }_n\le {\tau }_0+\tau (\epsilon )<{\tau }_0+\tau (\epsilon )+\epsilon =\tau .$$(38)

If $r_{*}=1$, then ${\xi }_n={\tau }_0$ for all n and the above inequality holds, so suppose that $r_{*}>1$.

Assume that (38) does not hold. Then there exists an increasing sequence $\mathcal{K}$ of natural numbers such that ${\xi }_k>{\tau }_0+\tau (\epsilon )$ for all $k\in \mathcal{K}$. Moreover, by Propositions 3 and 6, this sequence may be chosen so that the next two properties hold.

• (i) The sequence $\{{\mathit{q}}_k:k\in \mathcal{K}\}$ converges in $D_{{\ell }_1}[0,\infty )$.

• (ii) We have ${\mathit{r}}_k(t)\le r_{*}$ for all $t\in [{\tau }_0,T]$ and $k\in \mathcal{K}$.

The definition of ξ_k implies that

$${\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}{\mathit{q}}_k}(t,i,j)<{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}{\alpha }_k}(i)\text{for}\text{all}t\in \left[{\tau }_0,{\tau }_0+\tau (\epsilon )\right]\subset \left[{\tau }_0,{\xi }_k\right)\text{and}k\in \mathcal{K}.$$

The hypotheses of Proposition 5 hold with $r≔r_{*},\hspace{0.17em}{t}_0≔{\tau }_0$ and $t_1≔{\tau }_0+\tau (\epsilon )$, by the previous remark and properties (i) and (ii). Let $\mathit{v}(r_{*})$ be the function defined in this proposition, as the uniform limit of the tail processes ${\mathit{v}}_k(r_{*})$ over $[0,T]$. It follows from Propositions 4 and 5 that

$$\underset{t\in [0,T]}{\sup }|{\mathit{s}}_k(t)-{\mathit{v}}_k(t,r_{*})-\left[\mathit{s}(t)-\mathit{v}(t,r_{*})\right]|\le \frac{\epsilon }{2}$$

for all sufficiently large $k\in \mathcal{K}$. For each of these k, we have

$$\begin{array}{ll}{\displaystyle \sum _{(i,j)⊳{\sigma }_{*}}{\mathit{q}}_k}\left({\tau }_0+\tau (\epsilon ),i,j\right)\hfill & ={\mathit{s}}_k\left({\tau }_0+\tau (\epsilon )\right)-{\mathit{v}}_k({\tau }_0+\tau (\epsilon ),r_{*})\hfill \\ \hfill & \ge \mathit{s}\left({\tau }_0+\tau (\epsilon )\right)-\mathit{v}({\tau }_0+\tau (\epsilon ),r_{*})-\frac{\epsilon }{2}\hfill \\ \hfill & >\rho +\left[\mathit{s}({\tau }_0)-\rho \right]{\text{e}}^{-\mu \tau (\epsilon )}-\mathit{s}({\tau }_0){\text{e}}^{-\mu \tau (\epsilon )}-\frac{\epsilon }{2}\hfill \\ \hfill & =\rho \left(1-{\text{e}}^{-\mu \tau (\epsilon )}\right)-\frac{\epsilon }{2}={\displaystyle \sum _{(i,j)⊳(i_r,j_r)}\alpha }(i)+\frac{\epsilon }{2}.\hfill \end{array}$$

The third inequality follows from Proposition 5 and the last equality from the definition of $\tau (\epsilon )$. It follows from (9) that the right-hand side is strictly larger than ${\sum }_{(i,j)⊳{\sigma }_{*}}{\alpha }_k(i)$ for all large enough $k\in \mathcal{K}$, which is a contradiction.

We conclude that (38) holds, which proves (16a) and (16b). We had already proved (16c) in Corollary 1; thus, the proof of the theorem is complete. □

8. Asymptotic Optimality

In this section, we prove Theorem 2. Specifically, in Section 8.1, we use drift analysis to demonstrate that the continuous-time Markov chains introduced in Section 3.3 are irreducible and positive-recurrent and to derive upper bounds for certain expectations and tail probabilities. In Section 8.2, we use these upper bounds to establish that the stationary distributions of the latter Markov chains are tight, and then we complete the proof of Theorem 2 using the results of Sections 6 and 7. Finally, in Section 8.3, we demonstrate that JLMU is not optimal in general, although it is asymptotically optimal.

8.1. Drift Analysis

Denote the state space and the generator matrix of the continuous-time Markov chains defined in Section 3.3 by S_n and A_n, respectively. We use exactly the same notation for JLMU and SLTA, but we always indicate which policy is being considered. The drift of a function $f:S_n\to [0,\infty )$ is the function A_nf defined by

The proof of the following proposition uses a Foster-Lyapunov argument, which is based on the drift of certain suitably chosen functions.

Proposition 7.

For each given n, the two continuous-time Markov chains introduced in Section 3.3 are irreducible and positive-recurrent. In particular, each of these Markov chains has a unique stationary distribution π_n.

Proof.

Suppose first that the load balancing policy is JLMU. Any occupancy state can reach the empty occupancy state after a finite number of consecutive departures. By the definition of S_n provided in Section 3.3, the latter remark implies that ${\mathit{q}}_n$ is irreducible. Moreover, ${\mathit{q}}_n$ is the empty occupancy state if and only if ${\mathit{s}}_n=0$, which implies that the empty occupancy state is positive-recurrent, because the $M/M/\infty$ queue ${\mathit{s}}_n$ is irreducible and positive-recurrent. Thus, ${\mathit{q}}_n$ is positive-recurrent.

Suppose now that the load balancing policy is SLTA. Any state $(q,r)\in S_n$ can reach the empty occupancy sate with ${\mathit{r}}_n=r$ after a finite number of consecutive departures. Moreover, the latter state can reach the empty occupancy state with ${\mathit{r}}_n=1$ after a finite number of alternate arrivals and departures. We conclude from the definition of S_n provided in Section 3.3 that $({\mathit{q}}_n,{\mathit{r}}_n)$ is irreducible.

Next we use a Foster-Lyapunov argument to prove the positive recurrence. Consider the functions $f,g:S_n\to [0,\infty )$ defined by

$$f(q,r)≔{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}q}(i,j)\text{and}g(q,r)≔r\text{for}\text{all}(q,r)\in S_n.$$(39)

All server pools together form an infinite-server system; thus, $A_{n}f(q,r)=\lambda -\mu f(q,r)$ for all $(q,r)\in S_n$. In addition, we have

$$A_{n}g(q,r)=\lambda \left[{\mathrm{𝟙}}_I(q,r)-{\mathrm{𝟙}}_D(q,r)\right]\text{for}\text{all}(q,r)\in S_n.$$

Here I corresponds to those states (q, r) such that ${\mathit{r}}_n$ increases if $({\mathit{q}}_n,{\mathit{r}}_n)=(q,r)$ and the next event is an arrival. Specifically,

$$I≔\left\{(q,r)\in S_n:q(i,j)={\alpha }_n(i)\hspace{0.17em}\text{for}\text{all}\hspace{0.17em}(i,j)⊳(i_r,j_r),\hspace{0.17em}nq(i_r,j_r)=n{\alpha }_n(i_r)-1\right\}.$$

Also, D corresponds to those states (q, r) such that ${\mathit{r}}_n$ decreases if $({\mathit{q}}_n,{\mathit{r}}_n)=(q,r)$ and the next event is an arrival. Specifically,

$$D≔\left\{(q,r)\in S_n:r>1,\hspace{0.17em}n-{\displaystyle \sum _{i=1}^{m}n}q\left(i,{\ell }_i(r)\right)\ge n{\beta }_n,\hspace{0.17em}q(i_{r-1},j_{r-1})<{\alpha }_n(i_{r-1})\right\}.$$

Consider the function $h≔f+2g$ and let F be the set of those $(q,r)\in S_n$ that satisfy the following two conditions.

• (i) We have $f(q,r)\le 4\rho$.

• (ii) We have r = 1 or r > 1 and $\left[{\alpha }_n(i_{r-1})-{\beta }_n\right]j_{r-1}\le 4\rho$.

The first condition holds for finitely many $q\in Q_n$ and the second condition holds for finitely many $r\ge 1$, thus F is finite. Next we establish that $A_{n}h\le -\lambda +4\lambda {\mathrm{𝟙}}_F$. Note that $({\mathit{q}}_n,{\mathit{r}}_n)$ is nonexplosive because the infinite-server queue ${\mathit{s}}_n$ has this property. Therefore, it follows from Hajek (2006, proposition 2.20.1) that $({\mathit{q}}_n,{\mathit{r}}_n)$ is positive-recurrent.

The latter inequality holds for all $(q,r)\in F$ since $A_{n}h\le 3\lambda$. Hence, let us assume that $(q,r)\notin F$. Suppose that $(q,r)\notin F$ violates (i). Then

$$A_{n}h(q,r)\le 3\lambda -\mu f(q,r)<3\lambda -4\lambda =-\lambda .$$

Assume now that $(q,r)\notin F$ satisfies (i). Then (q, r) satisfies (i) and violates (ii), which implies that r > 1 and $\left[{\alpha }_n(i_{r-1})-{\beta }_n\right]j_{r-1}>4\rho \ge f(q,r)$. From this we conclude that $q(i_{r-1},j_{r-1})<{\alpha }_n(i_{r-1})-{\beta }_n$, because otherwise $f(q,r)\ge j_{r-1}q(i_{r-1},j_{r-1})>4\rho$. Thus,

$$n-{\displaystyle \sum _{i=1}^{m}n}q\left(i,{\ell }_i(r)\right)={\displaystyle \sum _{i=1}^{m}n}\left[{\alpha }_n(i)-q\left(i,{\ell }_i(r)\right)\right]\ge n\left[{\alpha }_n(i_{r-1})-q(i_{r-1},j_{r-1})\right]>n{\beta }_n.$$

It follows that $(q,r)\in D$; thus, $A_{n}h(q,r)=\lambda -\mu f(q,r)-2\lambda \le -\lambda$. □

Next we provide upper bounds for certain expectations and tail probabilities, which are used in the following section to demonstrate that the sequence of stationary distributions $\{{\pi }_n:n\ge 1\}$ is tight, both for JLMU and SLTA. First we state a technical lemma; the proof follows from Fubini’s theorem and is provided in Online Appendix A.

Lemma 8.

Let x_n have the stationary distribution π_n, where x_n = q_n if JLMU is used and $x_n=(q_n,r_n)$ if SLTA is used. If $f:S_n\to [0,\infty )$ satisfies

$$E\left[{\displaystyle \sum _{y\in S_n}|A_n(x_n,y)|}f(y)\right]<\infty ,\mathit{then}E\left[A_{n}f(x_n)\right]=0.$$

Consider the quantities

$${\theta }_n^k≔{\displaystyle \sum _{(i,j)\mathrm{⊵}(i_k,j_k)}{\alpha }_n}(i)\mathit{for}\mathrm{all}k\ge 1.$$(40)

This lemma is used to prove the following two propositions.

Proposition 8.

Suppose that the load balancing policy is JLMU, fix n and consider the functions $f_k:S_n\to [0,\infty )$ defined by

$$f_k(q)≔{\displaystyle \sum _{(i,j)⊲(i_k,j_k)}q}(i,j)\mathit{for}\mathrm{all}q\in S_n\mathit{and}k\ge 1.$$

If q_n has the stationary distribution π_n, then

$$E\left[f_k(q_n)\right]\le \rho {\text{e}}^{-n\left({\theta }_n^k-2\rho \right)}\mathit{for}\mathit{all}k\ge 1.$$

Proof.

Fix some $k\ge 1$ and define $J_i≔\min \{j\ge 1:(i,j)⊲(i_k,j_k)\}$. The drift of f_k with respect to ${\mathit{q}}_n$ satisfies

$$\begin{array}{ll}{A}_n{f}_k(q)\hfill & =\lambda {\mathrm{𝟙}}_{\{\sigma (q)⊲(i_k,j_k)\}}-{\displaystyle \sum _{(i,j)⊲(i_k,j_k)}\mu }j\left[q(i,j)-q(i,j+1)\right]\hfill \\ \hfill & =\lambda {\mathrm{𝟙}}_{\{\sigma (q)⊲(i_k,j_k)\}}-{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=J_i}^{\infty }\mu }}j\left[q(i,j)-q(i,j+1)\right]\hfill \\ \hfill & =\lambda {\mathrm{𝟙}}_{\{\sigma (q)⊲(i_k,j_k)\}}-{\displaystyle \sum _{i=1}^m\left[\mu \left(J_i-1\right)q\left(i,J_i\right)+\mu {\displaystyle \sum _{j=J_i}^{\infty }q}(i,j)\right]}\hfill \\ \hfill & =\lambda {\mathrm{𝟙}}_{\{\sigma (q)⊲(i_k,j_k)\}}-\mu f_k(q)-\mu {\displaystyle \sum _{i=1}^m\left(J_i-1\right)}q\left(i,J_i\right)\hfill \\ \hfill & \le \lambda {\mathrm{𝟙}}_{\{f(q)\ge {\theta }_n^k\}}-\mu f_k(q)\text{for}\text{all}q\in S_n,\hfill \end{array}$$(41)

where f(q) is as in (39). For the last step, $(i,j)⊳\sigma (q)$ implies that $q(i,j)={\alpha }_n(i)$.

Observe that $|A_n(x,x)|=n\lambda +n\mu f(x)$ because $nf(x)$ is the total number of tasks at the occupancy state x and $n\lambda$ is the arrival rate of tasks. Hence,

$$\begin{array}{ll}{\displaystyle \sum _{y\in S_n}|A_n(x,y)|}{f}_k(y)\hfill & =|A_n(x,x)|f_k(x)+{\displaystyle \sum _{y\ne x}{A}_n}(x,y)f_k(y)\hfill \\ \hfill & =2|A_n(x,x)|f_k(x)+{\displaystyle \sum _{y\in S_n}{A}_n}(x,y)f_k(y)\hfill \\ \hfill & =2\left[n\lambda +n\mu f(x)\right]f_k(x)+A_n{f}_k(x)\le 2n\left[\lambda +\mu f(x)\right]f(x)+\lambda .\hfill \end{array}$$

The right-hand side has a finite mean with respect to π_n because $nf(q_n)$ is the total number of tasks in the system in stationarity, which is Poisson distributed with mean $n\rho$. Thus, we conclude that $E\left[A_n{f}_k(q_n)\right]=0$ by Lemma 8.

Taking expectations with respect to π_n on both sides of (41) and recalling that $nf(q_n)$ is Poisson distributed with mean $n\rho$, we obtain

$$E\left[f_k\left(q_n\right)\right]\le \rho E\left[{\mathrm{𝟙}}_{\{f(q_n)\ge {\theta }_n^k\}}\right]=\rho P\left(nf(q_n)\ge n{\theta }_n^k\right)\le \rho {\text{e}}^{-n\left({\theta }_n^k-2\rho \right)},$$

where the last inequality follows from a Chernoff bound. □

Proposition 9.

Suppose that the load balancing policy is SLTA, fix n and consider the functions $f_k:S_n\to [0,\infty )$ defined by

$$f_k(q,r)≔{\displaystyle \sum _{(i,j)⊲(i_k,j_k)}q}(i,j)\mathit{for}\mathit{all}(q,r)\in S_n\mathit{and}k\ge 1.$$

Let (q_n, r_n) have the stationary distribution π_n. For each $k\ge 1$, we have

$$E\left[f_k(q_n,r_n)\right]\le \rho P\left(r_n>k\right),$$(42a)

$$P\left(r_n>k\right)\le {\text{e}}^{-n\left({\theta }_n^k-2\rho \right)}+{\text{e}}^{-n\left[{\theta }_n^k-L(k){\beta }_n-2\rho \right]}+{\text{e}}^{-n\left[\ell (k){\alpha }_n^{\min }-2\rho \right]},$$(42b)

where ${\theta }_n^k$ is defined as in (40), ${\alpha }_n^{\min }≔\min \{{\alpha }_n(i):1\le i\le m\}$,

$$L(k)≔\max \{{\ell }_i(k+1):1\le i\le m\}\mathit{and}\ell (k)≔\min \{{\ell }_i(k+1):1\le i\le m\}.$$

Proof.

Fix $k\ge 1$. As in the proof of Proposition 8, we see that

$$A_n{f}_k(q,r)\le \lambda {\mathrm{𝟙}}_{\{r>k\}}-\mu f_k(q,r)\text{for}\text{all}(q,r)\in S_n,$$

and that f_k satisfies the hypothesis of Lemma 8. Then we obtain (42a) by taking the expectation with respect to π_n on both sides of the latter inequality.

Consider the sets I and D defined in the proof of Proposition 7 and let

$$I_k≔I\cap \{(q,r)\in S_n:r\ge k\}\text{and}{D}_k≔D\cap \{(q,r)\in S_n:r>k\}.$$

The first step of the proof of (42b) is to establish that

$$P\left((q_n,r_n)\in I_k\right)=P\left((q_n,r_n)\in D_k\right).$$(43)

Fix l > k and consider the function $g_k^l:S_n\to [0,\infty )$ defined by

$$g_k^l(q,r)={(r-k)}^{+}{\mathrm{𝟙}}_{\{r<l\}}+(l-k){\mathrm{𝟙}}_{\{r\ge l\}}\text{for}\text{all}(q,r)\in S_n.$$

As in the proof of Proposition 7, we obtain

$$A_n{g}_k^l(q,r)=\lambda \left[{\mathrm{𝟙}}_{I_k^l}(q,r)-{\mathrm{𝟙}}_{D_k^l}(q,r)\right],$$(44)

where the sets $I_k^l$ and $D_k^l$ are defined by

$$I_k^l≔I\cap \{(q,r)\in S_n:k\le r<l\}\text{and}{D}_k^l≔D\cap \{(q,r)\in S_n:k<r\le l\}.$$

Define f(q, r) as in (39) and note that $|A_n\left((x,r),(x,r)\right)|=n\lambda +n\mu f(x,r)$. Thus,

$$\begin{array}{ll}{\displaystyle \sum _{(y,s)\in S_n}|A_n\left((x,r),(y,s)\right)|}{g}_k^l(y,s)\hfill & =2|A_n\left((x,r),(y,s)\right)|g_k^l(x,r)+A_n{g}_k^l(x,r)\hfill \\ \hfill & \le 2n\left[\lambda +\mu f(x,r)\right](l-k)+\lambda .\hfill \end{array}$$

The right-hand side has a finite mean with respect to π_n because $nf(q_n,r_n)$ is the total number of tasks in stationarity, which is Poisson distributed with mean $n\rho$. Therefore, it follows from Lemma 8 and (44) that $P\left((q_n,r_n)\in I_k^l\right)=P\left((q_n,r_n)\in D_k^l\right)$. The sets $I_k^l$ and $D_k^l$ increase to I_k and D_k, respectively, as $l\to \infty$. This implies (43) because

$$P\left((q_n,r_n)\in I_k\right)=\underset{l\to \infty }{\lim }P\left((q_n,r_n)\in I_k^l\right)=\underset{l\to \infty }{\lim }P\left((q_n,r_n)\in D_k^l\right)=P\left((q_n,r_n)\in D_k\right).$$

Now we may write

$$\begin{array}{ll}P\left(r_n>k\right)\hfill & =P\left(r_n>k,(q_n,r_n)\in D\right)+P\left(r_n>k,(q_n,r_n)\notin D\right)\hfill \\ \hfill & =P\left((q_n,r_n)\in D_k\right)+P\left(r_n>k,(q_n,r_n)\notin D\right)\hfill \\ \hfill & =P\left((q_n,r_n)\in I_k\right)+P\left(r_n>k,(q_n,r_n)\notin D\right).\hfill \end{array}$$(45)

Using the definition of I_k, we can bound the first term on the last line by

$$P\left(f(q_n,r_n)\ge {\displaystyle \sum _{(i,j)⊳(i_{r_n},j_{r_n})}{\alpha }_n}(i),r_n>k\right)\le P\left(f(q_n,r_n)\ge {\theta }_n^k\right).$$(46)

The second term on the last line of (45) can be bounded by

$$\begin{array}{ll}\hfill & P\left({\displaystyle \sum _{i=1}^m{q}_n}\left(i,{\ell }_i(r_n)\right)>1-{\beta }_n,r_n>k\right)+P\left(q_n(i_{r_n-1},j_{r_n-1})={\alpha }_n(i_{r_n-1}),r_n>k\right)\le \hfill \\ \hfill & P\left(f(q_n,r_n)\ge {\theta }_n^k-L(k){\beta }_n\right)+P\left(f(q_n,r_n)\ge \ell (k){\alpha }_n^{\min }\right).\hfill \end{array}$$(47)

For the last inequality, observe that the condition inside the first probability sign of the left-hand side of (47) implies that

$${\displaystyle \sum _{i=1}^m{q}_n}\left(i,{\ell }_i(k+1)\right)\ge {\displaystyle \sum _{i=1}^m{q}_n}\left(i,{\ell }_i(r_n)\right)>1-{\beta }_n,$$

and this in turn implies that

$$f(q_n,r_n)\ge {\displaystyle \sum _{(i,j)\mathrm{⊵}(i_k,j_k)}{\alpha }_n}(i)-L(k){\beta }_n={\theta }_n^k-L(k){\beta }_n,$$

because $q_n(i,j)$ is nonincreasing in j for all i. In addition, the condition inside the second probability sign of the left-hand side of (47) implies that

$$f(q_n,r_n)\ge j_{r_n-1}{q}_n(i_{r_n-1},j_{r_n-1})=j_{r_n-1}{\alpha }_n(i_{r_n-1})\ge \ell (k){\alpha }_n^{\min }.$$

We obtain (42b) from (46) and (47), recalling that $nf(q_n,r_n)$ is Poisson distributed with mean $n\rho$ and applying Chernoff bounds. □

8.2. Proof of the Asymptotic Optimality

In this section, we prove Theorem 2. As a first step, we establish that the sequence of stationary distributions $\{{\pi }_n:n\ge 1\}$ is tight both for JLMU and SLTA.

Proposition 10.

If the load balancing policy is JLMU, then $\{{\pi }_n:n\ge 1\}$ is tight in ${\ell }_1$. If the load balancing policy is SLTA, then $\{{\pi }_n:n\ge 1\}$ is tight in ${\ell }_1\times \mathbb{N}$.

Proof.

Suppose first that the load balancing policy is JLMU and let q_n have the stationary distribution π_n for each n. The sequence $\{q_n:n\ge 1\}$ is tight with respect to the product topology because the random variables q_n take values in ${[0,1]}^{\mathcal{I}}$, which is compact with respect to the product topology. Therefore, as in Mukherjee et al. (2018, lemma 2), the tightness in ${\ell }_1$ of $\{q_n:n\ge 1\}$ will follow if we establish that

$$\underset{k\to \infty }{\lim }\underset{n\to \infty }{\lim \sup }P\left({\displaystyle \sum _{(i,j)⊲(i_k,j_k)}{q}_n}(i,j)>\epsilon \right)=0\text{for}\text{all}\epsilon >0.$$(48)

By Proposition 8 and Markov’s inequality, we have

$$P\left({\displaystyle \sum _{(i,j)⊲(i_k,j_k)}{q}_n}(i,j)>\epsilon \right)\le \frac{\rho }{\epsilon }{\text{e}}^{-n\left({\theta }_n^k-2\rho \right)}\text{for}\text{all}k,n\ge 1\text{and}\epsilon >0.$$

For all sufficiently large k, the exponent on the right-hand side converges to minus infinity as n grows large and k is held fixed. Thus, $\{q_n:n\ge 1\}$ is tight in ${\ell }_1$.

Suppose now that the load balancing policy is SLTA and let (q_n, r_n) have the stationary distribution π_n for each n. To prove that $\{(q_n,r_n):n\ge 1\}$ is tight in ${\ell }_1\times \mathbb{N}$, it suffices to show that $\{q_n:n\ge 1\}$ and $\{r_n:n\ge 1\}$ are tight in ${\ell }_1$ and $\mathbb{N}$, respectively. Indeed, if the latter properties hold, then for each $\epsilon >0$ there exist compact sets $K_q\subset {\ell }_1$ and $K_r\subset \mathbb{N}$ such that

$$\max \{P\left(q_n\notin K_q\right),P\left(r_n\notin K_r\right)\}\le \frac{\epsilon }{2}\text{for}\text{all}n.$$

Therefore, the compact set $K_q\times K_r\subset {\ell }_1\times \mathbb{N}$ satisfies

$$P\left((q_n,r_n)\notin K_q\times K_r\right)\le P\left(q_n\notin K_q\right)+P\left(r_n\notin K_r\right)\le \epsilon \text{for}\text{all}n.$$

By Proposition 9 and Markov’s inequality, we have

$$\begin{array}{ll}\hfill & P\left({\displaystyle \sum _{(i,j)⊲(i_k,j_k)}{q}_n}(i,j)>\epsilon \right)\le \frac{\rho }{\epsilon }P\left(r_n>k\right)\text{for}\text{all}k,n\ge 1\text{and}\epsilon >0,\hfill \\ \hfill & P\left(r_n>k\right)\le {\text{e}}^{-n\left({\theta }_n^k-2\rho \right)}+{\text{e}}^{-n\left[{\theta }_n^k-L(k){\beta }_n-2\rho \right]}+{\text{e}}^{-n\left[\ell (k){\alpha }_n^{\min }-2\rho \right]}\text{for}\text{all}k,n\ge 1.\hfill \end{array}$$

For all large enough k, the right-hand side of the second inequality is summable over n and in particular vanishes with n. This implies that the sequences $\{q_n:n\ge 1\}$ and $\{r_n:n\ge 1\}$ are tight in ${\ell }_1$ and $\mathbb{N}$, respectively. □

We also need the following technical lemma.

Lemma 9.

Fix $q\in {\ell }_1$ and suppose that the marginal utilities are bounded. Then

$$u(q)={\displaystyle \sum _{i=1}^m{u}_i}(0)q(i,0)+{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathrm{\Delta }}(i,j-1)q(i,j).$$

Proof.

Note that

$$\begin{array}{ll}u(q)\hfill & ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^{\infty }{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ \hfill & =\underset{k\to \infty }{\lim }{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=0}^k{u}_i}}(j)\left[q(i,j)-q(i,j+1)\right]\hfill \\ \hfill & =\underset{k\to \infty }{\lim }{\displaystyle \sum _{i=1}^m\left[u_i(0)q(i,0)+{\displaystyle \sum _{j=1}^k\mathrm{\Delta }}(i,j-1)q(i,j)-u_i(k)q(i,k+1)\right]},\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^m{u}_i}(0)q(i,0)+{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathrm{\Delta }}(i,j-1)q(i,j)-\underset{k\to \infty }{\lim }{\displaystyle \sum _{i=1}^m{u}_i}(k)q(i,k+1).\hfill \end{array}$$

The second term in the last line is absolutely convergent because the marginal utilities are bounded and $q\in {\ell }_1$. Moreover, there exists $M\ge 0$ such that

$$\begin{array}{l}\underset{k\to \infty }{\lim }|u_i(k)|q(i,k+1)\hfill \\ \le \underset{k\to \infty }{\lim }\left|u_i(0)+{\displaystyle \sum _{j=0}^{k-1}\mathrm{\Delta }}(i,j)\right|q(i,k+1)\hfill \\ \le \underset{k\to \infty }{\lim }\left(|u_i(0)|+Mk\right)q(i,k+1)\text{for}\text{all}i,\end{array}$$

and the latter limit is zero for all i because $q\in {\ell }_1$. □

Now we are ready to prove Theorem 2.

Proof of Theorem 2.

Suppose first that the load balancing policy is JLMU. It follows from Prokhorov’s theorem and Proposition 10 that the stationary distributions $\{{\pi }_n:n\ge 1\}$ are relatively compact in ${\ell }_1$; thus, every subsequence has a further subsequence that converges in distribution. To establish (a), it suffices to prove the following statement: if $\mathcal{K}$ is an increasing sequence of natural numbers such that $\{{\pi }_k:k\in \mathcal{K}\}$ converges weakly to π, then π is the Dirac measure concentrated at $q_{*}$. Similarly, the sequence $\{{\pi }_n:n\ge 1\}$ is relatively compact in ${\ell }_1\times \mathbb{N}$ if the load balancing policy is SLTA, and to prove (b), it suffices to establish the following statement: if $\mathcal{K}$ is an increasing sequence of natural number such that $\{{\pi }_k:k\in \mathcal{K}\}$ converges weakly to π, then π is the Dirac measure at $(q_{*},r_{*})$. Here we prove (a) and (b) in parallel, proceeding as indicated previously.

Fix an arbitrary increasing sequence of natural numbers $\mathcal{K}$ such that $\{{\pi }_k:k\in \mathcal{K}\}$ converges weakly to a certain probability measure π. The following constructions use Skorokhod’s representation theorem. If the load balancing policy is JLMU, then there exist random variables q_k and q, distributed as π_k and π, respectively, that are defined on a common probability space $({\mathrm{\Omega }}_I,{\mathcal{F}}_I,{\mathbb{P}}_I)$ and satisfy

$$\underset{k\to \infty }{\lim }{\Vert q_k(\omega )-q(\omega )\Vert }_1=0\text{for}\text{all}\omega \in {\mathrm{\Omega }}_I.$$

If the load balancing policy is SLTA, then there exist random variables (q_k, r_k) and (q, r), with distributions given by the probability measures π_k and π, respectively, that are defined on some common probability space $({\mathrm{\Omega }}_I,{\mathcal{F}}_I,{\mathbb{P}}_I)$ and satisfy

$$\underset{k\to \infty }{\lim }{\Vert q_k(\omega )-q(\omega )\Vert }_1=0\hspace{0.17em}\text{and}\hspace{0.17em}\underset{k\to \infty }{\lim }{r}_k(\omega )=r(\omega )\text{for}\text{all}\omega \in {\mathrm{\Omega }}_I.$$

The second limit implies that there exists a random variable R such that $r_k(\omega )\le R(\omega )$ for all $k\in \mathcal{K}$ and $\omega \in {\mathrm{\Omega }}_I$. Moreover, $q_k(i,j)<{\alpha }_k(i)$ for all $(i,j)\mathrm{⊴}(i_{r_k},j_{r_k})$ on Ω_I by the definition of the state space S_k for SLTA. Indeed, recall from Remark 1 that the latter property is preserved by arrivals and departures and observe that it holds for the empty occupancy state with ${\mathit{r}}_k=1$.

If the load balancing policy is JLMU, then we may construct occupancy processes ${\mathit{q}}_k$ on a common probability space $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ as in Section 5.1.1, such that ${\mathit{q}}_k(0)=q_k$ and (9) holds with $q_0=q$. If the load balancing policy is SLTA, then we may construct processes $({\mathit{q}}_k,{\mathit{r}}_k)$ on a common probability space $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ as in Section 5.1.2, such that $\left({\mathit{q}}_k(0),{\mathit{r}}_k(0)\right)=(q_k,r_k)$ and the assumptions of Section 4.2.1 hold with $q_0=q$. For both load balancing policies, we may assume by Proposition 3 that ${\mathit{q}}_k$ converges in $D_{{\ell }_1}[0,\infty )$ to a process $\mathit{q}$ with a probability of one; this may require to replace $\mathcal{K}$ by a subsequence.

If the load balancing policy is JLMU, then Theorem 5 implies that $\mathit{q}(\omega )$ is the unique fluid trajectory such that $\mathit{q}(\omega ,0)=q(\omega )$ for each $\omega \in \mathrm{\Omega }$. Moreover, by Theorem 4,

$$\underset{t\to \infty }{\lim }{\left|\left|\mathit{q}(\omega ,t)-q_{*}\right|\right|}_1=0\text{for}\text{all}\omega \in \mathrm{\Omega }.$$(49)

If the load balancing policy is SLTA, then (16b) and (16c) hold on a set of probability of one by Theorem 6. Fix any $\omega \in \mathrm{\Omega }$ such that ${\mathit{q}}_k(\omega )$ converges to $\mathit{q}(\omega )$ in $D_{{\ell }_1}[0,\infty )$ and such that (16b) and (16c) hold. Also, choose $\tau >{\tau }_{\text{eq}}(s_0(\omega ))$. It follows from (16b) that

$$\begin{array}{ll}|\alpha (i)-\mathit{q}(\omega ,T,i,j)|\hfill & \le \underset{t\in [\tau ,T]}{\sup }|\alpha (i)-\mathit{q}(\omega ,t,i,j)|\hfill \\ \hfill & \le \underset{k\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }|\alpha (i)-{\mathit{q}}_k(\omega ,t,i,j)|+\underset{k\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }{\Vert {\mathit{q}}_k(\omega ,t)-\mathit{q}(\omega ,t)\Vert }_1=0,\hfill \end{array}$$

for all $(i,j)⊳{\sigma }_{*}$ and $T\ge \tau$. Thus, $\mathit{q}(\omega ,t,i,j)=q_{*}(i,j)$ for all $(i,j)⊳{\sigma }_{*}$ and $t\ge \tau$. Similarly, we conclude from (16c) that

$$\begin{array}{ll}{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\mathit{q}}(\omega ,T,i,j){\text{e}}^{\mu (T-\tau )}\hfill & \le \underset{t\in [\tau ,T]}{\sup }{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\mathit{q}}(\omega ,t,i,j){\text{e}}^{\mu (t-\tau )}\hfill \\ \hfill & \le \underset{k\to \infty }{\lim \sup }\underset{t\in [\tau ,T]}{\sup }{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}{\mathit{q}}_k}(\omega ,t,i,j){\text{e}}^{\mu (t-\tau )}\hfill \\ \hfill & +\underset{k\to \infty }{\lim }\underset{t\in [\tau ,T]}{\sup }{\Vert {\mathit{q}}_k(\omega ,t)-\mathit{q}(\omega ,t)\Vert }_1{\text{e}}^{\mu (t-\tau )}\le c(\omega ,\tau )\hfill \end{array}$$

for all $T\ge \tau$; for the last inequality, in the last line, ${\text{e}}^{\mu (t-\tau )}\le {\text{e}}^{\mu (T-\tau )}$ for all $t\in [\tau ,T]$, so the last limit equals zero. Because T can be arbitrarily large, we have

$$\underset{t\to \infty }{\lim }{\displaystyle \sum _{(i,j)⊲{\sigma }_{*}}\mathit{q}}(\omega ,t,i,j)\le \underset{t\to \infty }{\lim }c(\omega ,\tau ){\text{e}}^{-\mu (t-\tau )}=0,$$

and in particular $\mathit{q}(\omega ,t,i,j)\to q_{*}(i,j)$ as $t\to \infty$ for all $(i,j)⊲{\sigma }_{*}$. This also holds for $(i,j)={\sigma }_{*}$ by Proposition 4, so we conclude that

$$\mathbb{P}\left(\underset{t\to \infty }{\lim }{\Vert \mathit{q}(t)-q_{*}\Vert }_1=0\right)=1.$$(50)

It follows from the stationarity of π_k that ${\mathit{q}}_k(t)$ is distributed as q_k for all $t\ge 0$ in the case of JLMU and that $\left({\mathit{q}}_k(t),{\mathit{r}}_k(t)\right)$ has the same distribution as (q_k, r_k) for all $t\ge 0$ if the load balancing policy is SLTA. Furthermore, recall that in either case, we have

$$\mathbb{P}\left(\underset{k\to \infty }{\lim }{\Vert {\mathit{q}}_k(t)-\mathit{q}(t)\Vert }_1=0\right)=1\text{for}\text{all}t\ge 0\text{and}{q}_k\Rightarrow q\text{as}k\to \infty ,$$

which implies that $\mathit{q}(t)$ is distributed as q for all $t\ge 0$. Moreover, $\mathit{q}(t)$ converges weakly to $q_{*}$ as $t\to \infty$ by (49) and (50). Therefore, q corresponds to the Dirac probability measure concentrated at $q_{*}$ both for JLMU and SLTA.

This completes the proof of (a). To finish the proof of (b), observe that $s_0=\rho$ with a probability of one because $q_0=q$ is equal to $q_{*}$ with a probability of one. It follows from (16a) that

$$\mathbb{P}\left(\underset{k\to \infty }{\lim }{\mathit{r}}_k(t)=r_{*}\right)=1\text{for}\text{all}t>{\tau }_{\text{eq}}(\rho ).$$

Hence, ${\mathit{r}}_k(t)$ converges weakly to $r_{*}$ for all $t>{\tau }_{\text{eq}}(\rho )$. Because ${\mathit{r}}_k(t)$ has the same distribution as r_k for all $t\ge 0$, we conclude that r_k converges weakly to $r_{*}$. Note that $q_{*}$ and $r_{*}$ are deterministic; thus, q_k and r_k converge to $q_{*}$ and $r_{*}$, respectively, in probability, which implies that (q_k, r_k) converges to $(q_{*},r_{*})$ in probability. This completes the proof of (b).

Next we prove the statements about the stationary overall utilities, and here exactly the same arguments apply both for JLMU and SLTA. Suppose that q_n has the stationary distribution. By (a) and (b), q_n converges in probability to $q_{*}$ in ${\ell }_1$, and by Lemma 9,

$$u(x)={\displaystyle \sum _{i=1}^m{u}_i}(0)x(i,0)+{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}\mathrm{\Delta }}(i,j-1)x(i,j)\text{for}\text{all}x\in {\ell }_1.$$

Because the marginal utilities are bounded, there exists a constant a such that

$$|u(x)-u(y)|\le a{\Vert x-y\Vert }_1\text{for}\text{all}x,y\in {\ell }_1.$$

This implies that $u(q_n)$ converges in probability to $u(q_{*})$.

Finally, observe that $|u(q_n)|\le b+c{s}_n$ with

$$b≔{\displaystyle \sum _{i=1}^m|}{u}_i(0)|,c≔\underset{(i,j)\in \mathcal{I}}{\max }|\mathrm{\Delta }(i,j)|\text{and}{s}_n≔{\displaystyle \sum _{(i,j)\in {\mathcal{I}}_{+}}{q}_n}(i,j).$$

Because $n{s}_n$ is Poisson distributed with mean $n\rho$, we have

$$\begin{array}{ll}E\left[u{(q_n)}^2\right]\hfill & \le b^2+2bcE\left[s_n\right]+c^{2}E\left[s_n^2\right]\le b^2+2bc\rho +c^2{\rho }^2+\frac{c^2\rho }{n}\le b^2+2bc\rho +c^2{\rho }^2+c^2\rho .\hfill \end{array}$$

The last expression does not depend on n; therefore, we conclude that $\{u(q_n):n\ge 1\}$ is uniformly integrable. As a result, we have

$$\underset{n\to \infty }{\lim }E\left[u(q_n)\right]=u(q_{*}).$$

This completes the proof. □

8.3. Suboptimality Result

In this section, we prove that JLMU is generally not optimal in the prelimit, although it becomes optimal as the number of server pools grows large. For this purpose, we construct an example in which JLMU is strictly outperformed by another policy.

Specifically, consider a system with two heterogeneous server pools and assume that the associated utility functions are of the following form:

where a > 0 and $\epsilon \in (0,1)$ are constants. Note that

As a result, JLMU sends tasks to server pool 2 if and only if this server pool is empty; thus, the number of tasks in server pool 2 is zero or one in stationarity. This assignment rule guarantees the largest increase in the aggregate utility of the system at each arrival epoch. However, this increase can be very small when tasks are sent to server pool 1, whereas the decrease in the aggregate utility can be comparatively large when a departure leaves server pool 2 with zero tasks. Particularly, this is the case when ε is small.

Suppose that tasks arrive as a Poisson process of intensity λ with exponential service times of mean $1/\mu$. In addition, let $\mathit{X}$ and $\mathit{Y}$ denote the number of tasks in server pools 1 and 2, respectively. Next we provide an upper bound for the mean stationary overall utility of JLMU when the utility functions are as defined previously.

As already noted, a new task is sent to server pool 2 if and only if $\mathit{Y}=0$. Therefore, the following statements hold.

• We have $\mathit{Y}(t)\in \{0,1\}$ for all sufficiently large t.

• Conditional on $\mathit{Y}(0)\in \{0,1\}$, the process $\mathit{Y}$ alone is a birth-death process with state space {0, 1}, birth rate λ, and death rate μ.

By Proposition 7, the process $(\mathit{X},\mathit{Y})$ has a unique stationary distribution (X, Y), and by the previous observations, we have

If we let $U≔u_1(X)+u_2(Y)$ denote the aggregate utility in stationarity, then

Consider now the policy that sends all tasks to server pool 2. Whereas JLMU dispatches the new tasks in a greedy fashion, this other policy is conservative because it tries to avoid drops of $u_2(\mathit{Y})$ from a to zero by keeping a positive number of tasks in server pool 2.

The mean stationary aggregate utility can be computed explicitly for the policy that we described previously because $\mathit{Y}$ is now an $M/M/\infty$ queue and $\mathit{X}(t)=0$ for all sufficiently large t. Let (X, Y) be the stationary distribution of $(\mathit{X},\mathit{Y})$ and let $V≔u_1(X)+u_2(Y)$ denote the aggregate utility of the system in stationarity, we have

It is not difficult to verify that ε can be chosen so that

for all ρ in some open interval contained in $(0,\infty )$. If ε and ρ are chosen so that the latter inequality holds, then $E\left[U\right]<E\left[V\right]$, and thus JLMU is strictly suboptimal.

In general, the right balance between greedy and conservative actions is difficult to determine and depends intricately on the set of utility functions. However, the benefits of conservative actions, which prevent the number of tasks in server pools of specific classes from dropping below certain occupancy levels, diminish as the scale of the system grows. Indeed, the average fraction of these server pools that have less tasks than the optimal quantity decreases as the number of server pools grows, even if the assignment policy is purely greedy as JLMU; essentially, this is a consequence of the increase in the number of server pools per class and the decrease in the coefficient of variation of the total number of tasks in the system. Therefore, the associated loss in mean stationary overall utility asymptotically vanishes as the scale of the system grows.