Learning and balancing unknown loads in large-scale systems

Consider a system of identical server pools where tasks with exponentially distributed service times arrive as a time-inhomogenenous Poisson process. An admission threshold is used in an inner control loop to assign incoming tasks to server pools while, in an outer control loop, a learning scheme adjusts this threshold over time to keep it aligned with the unknown offered load of the system. In a many-server regime, we prove that the learning scheme reaches an equilibrium along intervals of time where the normalized offered load per server pool is suitably bounded, and that this results in a balanced distribution of the load. Furthermore, we establish a similar result when tasks with Coxian distributed service times arrive at a constant rate and the threshold is adjusted using only the total number of tasks in the system. The novel proof technique developed in this paper, which differs from a traditional fluid limit analysis, allows to handle rapid variations of the first learning scheme, triggered by excursions of the occupancy process that have vanishing size. Moreover, our approach allows to characterize the asymptotic behavior of the system with Coxian distributed service times without relying on a fluid limit of a detailed state descriptor.


Introduction
Consider a service system where tasks are instantaneously dispatched to one of n identical server pools. All the tasks sharing a server pool are executed in parallel and the execution times do not depend on the number of tasks contending for service. Nevertheless, the portion of shared resources available to individual tasks does depend on the number of tasks sharing a server pool and determines the performance perceived by tasks. In order to optimize the overall experienced performance, the system needs to maintain a balanced distribution of the number of tasks across the server pools by dispatching the incoming tasks suitably.
The latter load balancing problem is particularly challenging in large-scale systems, where it is computationally prohibitive to store and manage complete state information. Motivated by this remark, we consider a threshold-based dispatching rule which maintains a balanced distribution of the load while it involves only limited state information and knowledge of the offered load of the system. This knowledge has to be acquired in an online manner to set the threshold optimally since the offered load is typically not exactly known or even time-varying. This paper examines the load balancing policies that result from blending the latter dispatching rule with one of two learning schemes for dynamically tracking the offered load. Namely, a basic scheme, which keeps track of the total number of tasks in the system and uses this quantity to estimate the offered load, and a refined scheme, which uses more detailed state information. An orthogonal dimension of the problem, which is also addressed in this paper, concerns the underlying traffic scenario; i.e., the distribution of the service times and the characteristics of the arrival process. We focus on two settings: first we assume that tasks with exponentially distributed service times arrive as a timeinhomogeneous Poisson process of intensity nλ, and that the refined learning scheme is used to adjust the threshold over time; then we assume that tasks with Coxian distributed service times arrive at a constant rate, and that the basic learning scheme is used to estimate the offered load. The analysis of the former setting emphasizes the learning aspect of the problem, whereas the latter setting shifts the spotlight towards the service time distribution. We also indicate how the results established in this paper can be extended to other possible combinations of the underlying traffic scenario and the learning scheme.
Altogether, the load balancing policies considered in this paper integrate two control loops. In the inner control loop, the dispatching rule aims at maintaining the number of tasks at each of the server pools between two given values; these values are always consecutive multiples of a positive integer ∆, which is a parameter of the policy, the lowest of these values is referred to as the threshold and is denoted n . In the outer control loop, the learning scheme adjusts the threshold over time to maintain the normalized offered load per server pool between n and n + ∆. Hence, the parameter ∆ regulates a tradeoff between the degree of load balance and the stability of the learning scheme. Specifically, as ∆ decreases, the dispatching rule is more stringent in terms of the load balancing objective and the learning scheme is more sensitive to demand variations.
The control actions of the dispatching rule do not influence the basic learning scheme since the total number of tasks in the system does not depend on these actions. In contrast, there is a complex interdependence between the inner and outer control actions when the refined learning scheme is used. Namely, the actions of the refined learning scheme are triggered by changes in the occupancy state of the system, while at the same time, the evolution of the occupancy state is governed by the dispatching rule, which depends on the output of the refined learning scheme.

Main contributions
Our main result for the refined scheme, in the time-inhomogeneous exponential scenario, is that the threshold reaches an equilibrium along intervals of time where the offered load is suitably bounded. Specifically, assume that there exist an integer m and an interval of time [a, b] where the offered load per server pool is bounded between m∆ and (m + 1)∆. We prove that there exists σ such that n = m∆ along [a + σ, b] for all large enough n with probability one. We also demonstrate that the fraction of server pools with less than m∆ tasks vanishes uniformly along the interval [a + σ, b] and that the fraction of server pools with more than (m + 1)∆ tasks decays at least exponentially fast with time; if ∆ = 1, then this means that the total number of tasks becomes evenly balanced across the server pools as n grows large. Analogous results are established for the basic learning scheme in the time-homogeneous Coxian setting; namely, we prove that the basic learning scheme reaches an equilibrium in all large enough systems and that the system asymptotically achieves a balanced distribution of the load, as described above.
In order to prove our main results, we develop a novel methodology which is fundamentally different from a traditional fluid limit analysis. In particular, our approach does not decouple the analysis across the scale parameter and the time variable as in a traditional fluid limit approach, where first a sequence of scaled processes is proven to converge to a solution of a certain system of differential equations and then stability results are established for this limiting dynamical system. Instead, our methodology leverages the tractability of specific system dynamics to identify key dynamical properties of the threshold and occupancy processes, and uses strong approximations to prove that these dynamical properties hold asymptotically as n grows large; these asymptotic dynamical properties are then employed to establish our main results.
An important motivation for our novel approach is the nature of the refined learning scheme, which allows for rapid changes of the threshold, triggered by small excursions of the occupancy process. Specifically, these excursions trip changes of the threshold if their size exceeds a parameter of the learning scheme which vanishes on the fluid scale as n grows large, making a traditional fluid limit analysis especially challenging, if not impossible. Furthermore, our approach allows to derive our main result for the Coxian setting in a direct way which avoids proving a fluid limit for a detailed state descriptor of the system; if tasks can undergo at most r phases, then such a state descriptor could be given by a multidimensional sequence in which the element (j 1 , . . . , j r ) indicates the fraction of server pools with exactly j m tasks in phase m.
The analysis of the refined learning scheme in the time-inhomogeneous exponential scenario is more challenging from a learning perspective, whereas the analysis of the basic learning scheme in the time-homogeneous Coxian scenario involves a higher complexity in terms of the service time distribution, since the residual service time distribution of the tasks in the system is heterogeneous. Nevertheless, we demonstrate that our methodology can be applied in both of these substantially distinct settings. Moreover, we believe that the principle underlying this methodology, of leveraging process-level dynamical properties through strong approximations, is of broader applicability, and can be particularly useful in situations where fluid limits are unavailable or difficult to derive.

Related work
The problem of balancing the load in parallel-server systems has received immense attention in the past few decades; a recent survey is provided in [4]. While traditionally the focus in this literature used to be on performance, more recently the implementation overhead has emerged as an equally important issue. This overhead has two sources: the communication burden of exchanging messages between the dispatcher and the servers, and the cost in large-scale deployments of storing and managing state information at the dispatcher [10,11].
As noted earlier, the present paper concerns an infinite-server setting where the execution times of tasks do not depend on the number of competing tasks. In contrast, the load balancing literature has been predominantly concerned with single-server models, where performance is measured in terms of queueing delays. In the latter scenario, the Join the Shortest Queue (JSQ) policy [8,31] minimizes the mean delay for exponentially distributed service times, but a naive implementation of JSQ involves an excessive communication burden for large-scale systems. So-called power-of-d strategies [20,21,30] involve substantially less communication overhead and yet provide significant improvements in delay performance over purely random routing. Another alternative are pull-based policies [3,27], which reduce the communication burden by maintaining state information at the dispatcher. Particularly, the Join the Idle Queue (JIQ) policy [18,27] asymptotically matches the optimality of JSQ and involves only one message per task, by keeping a list of idle queues.
With the exception of [27], the papers listed above only consider exponentially distributed service times, and their results have only been extended partially or under stringent assumptions to more general service times. In particular, [27,28] treat a JIQ system where service times are assumed to have a decreasing hazard rate function, whereas [9] considers a JIQ system with general service times but only under the assumption that the load is strictly below 1/2. Also, mean-field limits for power-of-d policies with generally distributed service times are examined in [6,7]. Assuming that a certain propagation of chaos property holds, the latter papers establish power-of-choice benefits similar to those originally reported in [20,30] for exponential service times. Power-of-d strategies were also studied in [1,2] through measure-valued processes, proving that a suitably scaled sequence of processes converges to a hydrodynamic limit. More closely related to the Coxian setting considered in the present paper is [29], which examines mean-field limits for power-of-d policies in loss systems with mixed-Erlang service times.
As mentioned above, the infinite-server setting considered in this paper has received only limited attention in the load balancing literature, even for exponentially distributed service times. While queueing delays are hardly meaningful in this setting, load balancing still plays a crucial role in optimizing different performance measures. Besides loss probabilities in finite-capacity scenarios, further relevant performance measures are the utility metrics introduced in [22], associated with quality-of-service as perceived by streaming applications. Specifically, if x i denotes the number of tasks sharing server pool i, then a utility metric of the form ϕ(x i ) := ψ(1/x i ) can be used as a proxy for measuring the quality-of-service provided to tasks assigned to server pool i, as a function of their resource share. Provided that ψ is concave and increasing, the overall utility n i=1 x i ϕ(x i ) is a Schur-concave function of the vector (x 1 , . . . , x n ) and is thus maximized by distributing the total number of tasks evenly across the various server pools. As in the single-server setting, JSQ is optimal for balancing the load [19,26], and it was proven in [22] that the performance of JSQ can be asymptotically matched on the diffusion scale by certain power-of-d strategies which considerably reduce the associated communication overhead.
Just like the dispatching algorithms studied in [22], the policies considered in this paper aim at optimizing the overall quality-of-service experienced by tasks. The refined version of our policy was previously considered in [12], assuming exponentially distributed service times and a constant arrival rate of tasks. Resorting to a traditional fluid limit analysis, [12] proves fluid-scale optimality provided that a parameter of the learning scheme satisfies a certain condition, stated in terms of the normalized offered load per server pool. The present paper establishes how this parameter of the refined learning scheme should scale with the number of server pools, and strengthens the latter result by proving, through a novel methodology, that fluid-scale optimality is always attained with the proposed parameter scaling, even if the arrival rate of tasks is time-varying.
As alluded to above, one of the most interesting features of our policies is their capability of adapting the threshold to unknown time-varying loads. The problem of adaptation to uncertain demand patterns was addressed in the context of single-server queueing models in [13,14,23,24], which assume that the number of servers can be right-sized on the fly to match the load. However, in all these papers the dispatching algorithm remains the same at all times since the right-sizing mechanism alone is sufficient to maintain small queueing delays when the demand for service changes. Different from these right-sizing mechanisms, the learning schemes considered in this paper modify the dispatching rule of the system over time to maintain a balanced distribution of the load while the number of server pools remains constant.

Outline of the paper
In Section 2 we specify the dispatching rule, the learning schemes and the traffic scenarios considered in this paper. The notation used throughout the paper and some technical assumptions are introduced in Section 3, and our main results are stated in Section 4, where we also outline our methodology and we discuss some extensions. The proof of our main result for the refined learning scheme, in the time-inhomogenenous exponential scenario, is carried out in Section 5, and the proof of our main result for Coxian distributed service times, in the time-homogeneous scenario with the basic learning scheme, is provided in Section 6. The proofs of some intermediate results are provided in Appendices A and B.

Model description
Consider a system of n identical server pools. Each server pool has infinitely many servers, tasks arrive as a (possibly time-varying) Poisson process of intensity nλ and service times are independent and identically distributed with mean 1/µ. In order to optimize the overall experienced performance of tasks, it is necessary to maintain an even distribution of the load. Specifically, if the total number of tasks at a given time is nρ, then the Schurconcave utility function defined in Section 1.2 is maximal when all server pools have either ρ or ρ tasks at that time; here · and · denote the floor and ceiling functions.

Dispatching rule
Suppose briefly that the arrival rate of tasks is constant over time and let ρ := λ/µ denote the offered load of the system; strictly speaking, ρ is the offered load per server pool. The total number of tasks in stationarity is Poisson distributed with mean nρ since all server pools together constitute an infinite-server system. To optimize the overall experienced performance of tasks, it is desirable that the stationary fraction of server pools with fewer than ρ or more than ρ tasks vanishes with n. For this purpose, we consider a dispatching algorithm that uses an integral parameter ∆ ≥ 1 and a threshold n ∈ N to distribute the incoming tasks among the server pools. Specifically, the dispatching policy aims at maintaining the number of tasks at each server pool between n and h n := n + ∆, by dispatching every incoming task as follows.
If some server pool has strictly less than n tasks right before the arrival, then the new task is sent to a server pool having strictly less than n tasks, chosen uniformly at random.
If all server pools have at least n tasks right before the arrival, and some server pool has strictly less than h n tasks, then the new task is sent to a server pool chosen uniformly at random among those having at least n and strictly less than h n tasks.
If all server pools have at least h n tasks right before the arrival, then the new task is sent to a server pool having at least h n tasks, chosen uniformly at random.
It was proved in [12] that the latter dispatching policy results in an evenly balanced distribution of the load as n grows large, assuming exponentially distributed service times and provided that ∆ = 1 and n = ρ for all n. In addition, this policy admits a tokenbased implementation which involves at most two messages per task and storage of two bits of information per server pool, as indicated in [12]. We observe that the insensitive load balancing policy studied in [15], in the context of parallel processor sharing servers, can also be implemented using tokens to determine routing probabilities based on a threshold.

Learning schemes
The dispatching rule described above relies on knowledge of the offered load for selecting the optimal threshold value, the one that yields a balanced distribution of the load as n grows large. In practice, the offered load is typically unknown and may experience changes over time due to demand fluctuations, which introduces the challenge of keeping the threshold aligned with an unknown and possibly time-dependent optimal value.

Basic learning scheme
A basic online procedure for learning the optimal threshold value relies on keeping track of the total number of tasks in the system N n and estimating the offered load as this quantity divided by the number of server pools. The optimal threshold value is then estimated as N n (t)/n for all t, as in Section 2.1. If the arrival rate of tasks is constant over time, then the total number of tasks in stationarity is Poisson distributed with mean nρ, so the latter procedure is asymptotically effective by a law of large numbers. Moreover, if the arrival rate of tasks is time-varying, then it is reasonable to expect a similar effectiveness along intervals of time where the arrival rate of tasks is roughly constant.
One caveat is that the threshold determined through this scheme could present quick fluctuations if the arrival rate of tasks has large enough high-frequency oscillations. Due to implementation considerations, it might be desirable to sacrifice load balance, to a certain extent, in order to gain stability of the threshold. This compromise is attained by increasing ∆ and setting and k ∈ Z for all t.

Refined learning scheme
The dispatching rule described in Section 2.1 aims at maintaining the number of tasks at each server pool between n and h n . An alternative learning scheme for finding the optimal threshold value consists of adjusting the threshold in the appropriate direction whenever the objective of the dispatching rule is violated. This more refined scheme is parameterized by α n ∈ (0, 1) and uses the same state information as the dispatching rule. In this scheme, the time-dependent threshold takes values on the non-negative multiples of ∆ and is adjusted in steps of ∆ units at certain arrival epochs. Specifically, when an arrival occurs, the refined learning scheme acts under the following circumstances.
If the fraction of server pools having at least n tasks is smaller than or equal to α n right before the arrival, then the threshold is decreased ∆ units once the new task has been dispatched.
If all the server pools have at least n tasks and the number of server pools with at least h n tasks is larger than or equal to n − 1 right before the arrival, then the threshold is increased ∆ units right after dispatching the task that has just arrived.
Observe that, eventually, the threshold is such that the number of server pools with at least h n tasks is smaller than or equal to n − 1 at all times; depending on the initial state of the system, this property may hold from time zero or after a certain number of arrivals. Once this property is in force, incoming tasks are always sent to server pools with at most h n − 1 tasks.
As for the basic learning scheme, the parameter ∆ regulates a tradeoff between load balance and stability of the threshold. Naturally, changes in the threshold will occur if the arrival rate of tasks exhibits large variations, but we will prove that the threshold remains constant and nearly optimal in a many-server regime during intervals of time where these variations are bounded by ∆. In fact, the threshold is optimal during these intervals provided that ∆ = 1, and the degree of optimality decreases with ∆ while the stability against demand fluctuations increases.

Traffic scenarios
The adaptive load balancing policies that result from the combination of the dispatching rule described in Section 2.1 and any of the two learning schemes introduced in Section 2.2 are analyzed in the context of two different traffic scenarios. Specifically, we first assume that tasks with exponentially distributed service times arrive as a time-inhomogeneous Poisson process, and we then asssume that tasks with Coxian distributed service times arrive as a time-homogenenous Poisson process. In the former scenario, our analysis focuses on the refined learning scheme, although our results easily extend to the basic learning scheme. In the latter scenario, our analysis focuses on the basic learning scheme, but most of our results extend to the refined learning scheme, as discussed in Section 4.3.
Note that the class of Coxian distributions is dense in the set of non-negative distributions with respect to the topology of weak convergence. In particular, any service time distribution can be approximated with an arbitrary precision by a Coxian distribution with respect to any metric that is compatible with this topology (e.g., the Lévy-Prokhorov metric). The denseness of the class of Coxian distributions was proved in [25] and can also be established as in [16,

Notation and assumptions
We consider a finite interval of time [0, T ], we assume that tasks arrive at a bounded rate λ : [0, T ] −→ [0, ∞) and we denote the instantaneous offered load per server pool by ρ(t) := λ(t)/µ. All tasks undergo at most r phases before they leave the system, they always enter the system in the first phase and they complete phase m after an exponentially distributed time of mean 1/µ m . After completing phase m, a task moves to the next phase with probability p m or leaves with probability 1 − p m , where p r = 0. Note that the case r = 1 corresponds to exponentially distributed service times and that the mean service time of a task is For technical reasons, we assume that p m ∈ (0, 1) for all m = r. The subclass of Coxian distributions that satisfy this property is dense in the class of all Coxian distributions with respect to the topology of weak convergence, hence dense in the set of all non-negative distributions. The pair formed by the threshold and the vector-valued process which describes the number of tasks in each phase at each server pool constitutes a Markov process. Instead of considering each server pool individually, we adopt an aggregate state description, denoting by s n (j 1 , . . . , j r ) the fraction of server pools with exactly j m tasks in phase m. In view of symmetry, the pair formed by the threshold n and the occupancy process i nqn(i) 1 · · · n · · · · · · · · · · · · Figure 1: Schematic representation of the occupancy state of the system for exponentially distributed service times. White circles represent servers and black circles represent tasks. Each row corresponds to a server pool, and these are arranged so that the number of tasks increases from top to bottom. The number of tasks in column i is nq n (i). s n := {s n (j 1 , . . . , j r ) : j 1 , . . . , j r ∈ N} also constitutes a Markov process. Sequences of these Markov processes, indexed by n, will be constructed on a common probability space (Ω, F, P) in Sections 5.1 and 6.1, for exponentially distributed service times and general Coxian distributed service times, respectively. We write s n (ω, t, j 1 , . . . , j r ) for the value of s n (j 1 , . . . , j r ) at (ω, t) ∈ Ω × [0, T ], although ω and t are often omitted.
The fraction of server pools with at least i tasks, in any phase, is denoted by and we let q n := {q n (i) : i ≥ 0}. If service times are exponentially distributed, then all the tasks in the system are identical, in the sense that they have the same residual service time distribution. In this case, s n can be recovered from q n , and we use the latter process to describe the occupancy state of the system, which can be represented as in the diagram of Figure 1. The latter diagram is also useful when service times are not exponentially distributed, although it only provides partial information about the occupancy state of the system since it ignores the phases of tasks. The occupancy process s n takes values in while the process q n takes values in Consider the product topology in R N r and assume that there exists s(0) ∈ S such that lim n→∞ s n (0) = s(0) almost surely.
We also assume that there exists B ∈ N such that q n (0, B + 1) = 0 and n (0) ≤ B for all n with probability one. This technical assumption precludes degenerate sequences of initial conditions. For example, occupancy states with a (possibly vanishing) fraction of server pools having a number of tasks that approaches infinity with n, or initial thresholds which approach infinity with n. Let G m denote the complementary cumulative distribution function of the residual service time of a task in phase m, and define u : The quantity u m (0) may be interpreted as the normalized total number of tasks in phase m for the limiting initial occupancy state s(0). The integral expression corresponds to the functional law of large numbers of an infinite-server system where the arrival rate of tasks scales proportionally to the time-varying function λ and the service times have the Coxian distribution specified above. In the case of exponential service times, we may also write where q(0) is defined in terms of s(0) as in (1). If λ is a regular enough function, then the latter integral equation can be stated in differential form asu(t) = λ(t) − µu(t).

Main results
In this section we state our main results, we discuss several extensions and we outline our methodology. In Section 4.1 we state our main result for the refined learning scheme, in the time-inhomogeneous exponential scenario, and we provide an overview of our proof technique. In Section 4.2 we state our main result for Coxian distributed service times, in the time-homogeneous scenario with the basic learning scheme, and we indicate the main steps of the proof. Extensions to our main results are discussed in Section 4.3.

Refined learning scheme
Assume that λ is time-varying, that service requirements are exponentially distributed and that the refined learning scheme is used. We consider the asymptotic behavior of the threshold and the occupancy state of the system over intervals of time where the instantaneous offered load is sufficiently well-behaved, as specified in the next definition.

Definition 1. Consider an integer m ≥ 0 and an interval
Moreover, assume that the length of [a, b] is strictly larger than We say that the offered load is (m, ∆)-bounded on [a, b] if the latter conditions hold, and we say that the offered load is ∆-bounded on [a, b] if it is (m, ∆)-bounded for some m.
Our main result regarding the refined learning scheme is presented below. Loosely speaking, it states that the threshold reaches an equilibrium during ∆-bounded intervals in all sufficiently large systems with probability one, provided that the control parameter α n approaches one with n at a suitable rate. Furthermore, the occupancy states are asymptotically nearly balanced under the latter conditions, with the degree of balance increasing as ∆ decreases.
If the offered load is (m, ∆)-bounded on an interval [a, b] and σ(a, b, m, ∆) < σ < b − a, then there exist c > 0 and a set of probability one where the following statements hold: The constant c can be taken equal to u(a + σ).
Since the threshold takes discrete values, (6a) implies that n (t) = m∆ for all t ∈ [a + σ, b] and all large enough n. Also, it follows from (6b) and (6c) that the number of server pools with less than m∆ tasks is uniformly o (n γ ) and the fraction of server pools with more than (m+1)∆ tasks is upper bounded by The constant σ(a, b, m, ∆) provides a closed-form upper bound for the asymptotic time until the threshold reaches an equilibrium. When the system is initially empty and ρ is constant over time, this upper bound is tight and corresponds to the asymptotic time until the total scaled number of tasks reaches m∆ for the first time.    14 16

Figure 2:
Numerical experiments with two different choices of ∆. The evolution of the thresholds is plotted in (a), the occupancy state of the system with ∆ = 3 along the 3-bounded interval [3,12] is depicted in (b) and the occupancy state of the system with ∆ = 1 along the 1-bounded interval [14,23] is plotted in (c). Both experiments correspond to initially empty systems with µ = 1, n = 300, α n = 1 − 1/n 0.48 and the arrival rate function plotted in (a); the superscripts in the legends indicate the value of ∆.
As ∆ decreases, the ∆-boundedness condition becomes more stringent, but the degree of load balance attained when the threshold reaches an equilibrium increases. In particular, if ∆ = 1, then m∆ = ρ(t) for all t ∈ [a, b] and (6a) implies that the threshold is optimal on the interval [a + σ, b] for all sufficiently large systems. Furthermore, (6b) and (6c) imply that the number of server pools with less than ρ(t) tasks is o (n γ ) and the fraction of server pools with more than ρ(t) tasks is upper bounded by The numerical experiments of Figure 2 illustrate the tradeoff between optimality and stability regulated by ∆. The experiments correspond to two systems operating under the same conditions, except for the choice of ∆. The time-varying offered load, the number of server pools and α n are the same in both systems, but while ∆ = 1 in one of the systems, ∆ = 3 in the other.
The offered load and threshold processes of the systems are plotted in Figure 2a. On a slower time scale, we observe that the offered load switches between two modes: the first one corresponds to a period of larger demand, during the interval [3,12], while the second one corresponds to a period of smaller demand, during [14,23]. On a faster time scale, the offered load exhibits rapid oscillations, which have larger amplitude along [3,12] and become smaller during [14,23]. The threshold of the system with ∆ = 3 is constant along both intervals, as one would expect from (6a) since both are 3-bounded. In contrast, the threshold of the system with ∆ = 1 is constant during the interval [14,23] but presents oscillations along [3,12], which is not 1-bounded.
The plots on the bottom of Figure 2 depict the occupancy state of the systems along intervals of time where the threshold remains constant. Figure 2b shows the occupancy state of the system with ∆ = 3 along the interval [3,12]. In consonance with (6b) and (6c), the fraction of server pools with fewer than (3) n tasks or more than (3) n + 3 tasks is nearly zero. However, the load is not perfectly balanced since the fraction of server pools with exactly (3) n + i tasks remains positive for all 0 ≤ i ≤ 3. In contrast, the load is perfectly balanced in the system with ∆ = 1 along [14,23], as Figure 2c shows; i.e., at all times virtually all server pools have either (1) n or (1) n + 1 tasks. As explained above, the system with ∆ = 1 manages to distribute the load evenly among the server pools during the interval [14,23], where the offered load only exhibits small oscillations. However, the larger demand fluctuations during the interval [3,12] result in an unstable threshold, which exhibits frequent oscillations along this interval. In contrast, the system with ∆ = 3 holds the threshold fixed along the intervals [3,12] and [14,23], although this comes at the expense of the load being only coarsely balanced along each of these intervals.

Outline of the proof
We now present an outline of the proof of Theorem 2, which is illustrated in Figure 3. This theorem concerns a sequence of occupancy and threshold processes constructed on a common probability space in Section 5.1, where it is also established that the occupancy processes have relatively compact sample paths with probability one. The latter property is proven in Proposition 4 through a technique developed in [5], which relies on strong approximations; namely, the strong law of large numbers for the Poisson process is used.
A slight refinement of the latter law of large numbers is derived in Section 5.2, where we define a set of sample paths for which this refined strong law of large numbers, and a few additional technical properties, hold. Propositions 8 and 9 prove that certain dynamical properties hold asymptotically within this set of sample paths, which has probability one. These properties concern the total and tail mass processes, defined as respectively. The former represents the total number of tasks in the system, normalized by the number of server pools. The latter corresponds to the total number of tasks in all but the first j − 1 columns of the diagram depicted in Figure 1, also normalized by the number of server pools. Proposition 8 establishes that u n converges uniformly over [0, T ] to the Proposition 4 Relative compactness of occupancy processes

Theorem 2 Proposition 10
Part of the proof of (6a)

Lemmas 12 and 13
Leveraging refined strong approximations

Proposition 9
Asymptotic dynamical properties of tail mass

Proposition 8
Asymptotic dynamical properties of total mass

Proposition 7
Nice sample paths (refined strong approximations) Section 5.1 Section 5.3 Section 5.2 function u defined by (4) with probability one, which is the strong law of large numbers of an infinite-server system. Proposition 9 is conceived to be used in conjunction with Proposition 4 since it assumes that a convergent subsequence of occupancy process sample paths is given. This proposition provides an asymptotic upper bound for a certain tail mass process under suitable conditions on the threshold and occupancy processes. Essentially, the conditions ensure that server pools with at least a certain number j of tasks are not assigned any further tasks, and the asymptotic upper bound implies that v n (j + 1) decays at least exponentially fast in a many-server regime. The proof of Theorem 2 is completed in Section 5.3 through several steps. First Proposition 8 is used to prove part of (6a). Specifically, it is established in Proposition 10 that n is almost surely upper bounded by m∆ along [a, b] minus a right neighborhood of a, for all large enough n. Next Propositions 4, 9 and 10 are used to establish (6c). As mentioned above, Propositions 4 and 9 are used in conjunction. Loosely speaking, under the assumption that (6c) does not hold, it is possible to construct a convergent subsequence of occupancy processes violating (6c), but Proposition 9 leads to a contradiction. Finally, Lemmas 12 and 13 leverage the refined strong approximations derived in Section 5.2 to analyze the occupancy processes on the scale of n γ 0 . This analysis and the asymptotic dynamical properties derived in Section 5.2 are then used to finish the proof of (6a) and to also prove (6b), completing the proof of Theorem 2.

Coxian service times
Assume that λ is constant over time, that service times are Coxian distributed and that the basic learning scheme is used. Our main result in this setting is that the threshold eventually reaches an equilibrium value in a many-server regime with probability one, and that the load becomes nearly balanced after the threshold settles.
Then there exist constants σ, η, c > 0 and a set of probability one where the following three statements hold: The constants σ, η and c do not depend on T .
This result has similar consequences as Theorem 2. Namely, (8a) implies that n (t) = ρ ∆ for all t ∈ [σ, T ] and all large enough n. Also, (8b) and (8c) imply that the total number of tasks in the system is asymptotically nearly balanced across the server pools during the interval [σ, T ]. Furthermore, the load is evenly balanced if ∆ = 1; i.e., the fraction of server pools with less than ρ tasks vanishes with n and the fraction of server pools with more than ρ tasks is upper bounded by [c + o(1)] e −η(t−σ) . The technical condition ρ / ∈ {i∆ : i ∈ N} ensures that the threshold does not oscillate between ρ − ∆ and ρ when the offered load is integral and a multiple of ∆. In addition, the constant η is a lower bound for the hazard rate function of the service time distribution, and is equal to µ for exponentially distributed service requirements. The proof of Theorem 3 is based on the methodology outlined in Section 4.1.1. Specifically, we first construct a sequence of occupancy and threshold processes on a common probability space from a set of stochastic primitives, and we establish that the sequence of occupancy processes is relatively compact with probability one; this is done in Section 6.1. Then we identify certain dynamical properties of the system and prove that they hold asymptotically in Section 6.2; these properties concern aggregate quantities analogous to the total and tail mass processes introduced in (7). Finally, the latter properties are used in Section 6.3 to prove Theorem 3.

Extensions
Our proof of Theorem 2 can be easily adapted to prove an analogous result for the basic learning scheme. Specifically, assume that the hypotheses of Theorem 2 hold, except that the basic learning scheme is used instead of the refined learning scheme, and let It follows from (4) that m∆ < u(t) < (m + 1)∆ for all a +σ eq (a, b, m, ∆) < t < b.
Using a strong law of large numbers for the total number of tasks in the system, it can be proved that (6a) holds for all σ >σ eq (a, b, m, ∆); the law of large numbers is provided in then it is possible to prove that, with probability one, for all sufficiently large n there exists a +σ eq (a, b, m, ∆) < t < a + σ such that q n (t, m∆) = 1. For any given γ ∈ [0, 1/2), it can then be established, as in Lemma 13, that q n (t, m∆) ≥ 1 − n −γ for all t ∈ [a + σ, b] and all large enough n with probability one. The latter statement implies (6b), and (6c) can be proven as in Corollary 11. It is also possible to prove that (8b) and (8c) hold if the hypotheses of Theorem 3 hold except that the basic learning scheme is replaced by the refined learning scheme. Indeed, as in the proof of Proposition 10, it can be shown that there exists σ > 0 such that n (t) ≤ ρ ∆ for all t ≥ σ and all large enough n with probability one. Then (8b) and (8c) can be established as in the proof of Theorem 3. Nevertheless, proving that the threshold eventually reaches an equilibrium in all sufficiently large systems would require non-trivial additional arguments. In particular, the arguments used in Section 5.3 to establish this property in the case of exponentially distributed service times cannot be easily extended to a scenario with general Coxian service times.
As a final remark, we observe that the statements in Theorem 2 are only meaningful because the constant σ(a, b, m, ∆) can be evaluated and is reasonably small in many situations. The derivation of this constant heavily relies on the simplicity of (4) for exponential service times, whereas the more complex equation (3) holds for general Coxian service times. Deriving a suitable analog of σ(a, b, m, ∆) would be challenging in a timeinhomogeneous scenario with Coxian distributed service times, even for the basic learning scheme. Furthermore, a closed-form expression as the one provided in Definition 1 is not likely to be available in the latter setting.

Refined learning scheme
In this section we prove Theorem 2. For this purpose, we assume that the refined learning scheme is used, that λ is time-varying and that service times are exponentially distributed. We proceed as indicated in Section 4.1.1.

Relative compactness of occupancy processes
In Section 5.1.1 we construct the sample paths of the occupancy states q n and the thresholds n on a common probability space for all n.
In Section 5.1.2 we use a methodology developed in [5] to establish that {q n : n ≥ 1} is a relatively compact subset of D R N[0, T ] with probability one.

Coupled construction of sample paths
The sample paths of the occupancy states and thresholds are defined as deterministic functions of the following stochastic primitives. Consider the product probability space of (Ω D , F D , P D ), (Ω S , F S , P S ) and (Ω I , F I , P I ); denote its completion by (Ω, F, P). The occupancy states and thresholds are defined on the latter space from a set of equations involving the fundamental processes and random variables introduced above. In order to write these equations, it is convenient to introduce some notation.
For each occupancy state q ∈ Q, we define the intervals These intervals form a partition of [0, 1) such that the length of I i (q) is the fraction of server pools with precisely i − 1 tasks. If q(j) < 1 for some j ∈ N, then we also define These intervals yield another partition of [0, 1). In this case, the length of J i (q, j) is the fraction of server pools with exactly i − 1 tasks, but only among those server pools with at most j − 1 tasks; we adopt the convention that J i (q, j) := ∅ for all 1 ≤ i ≤ j if q(j) = 1.
If ∈ N represents a threshold and q( ) < 1, then the length of J i (q, ) is equal to the probability of picking a server pool with precisely i − 1 tasks uniformly at random among those with strictly less than tasks. Similarly, if h := + ∆, q( ) = 1 and q(h) < 1, then the length of J i (q, h) is equal to the probability of picking a server pool with exactly i − 1 tasks uniformly at random among those server pools with at least tasks and strictly less than h tasks. We define for all i, j ≥ 1.
Note that r ij (q, ) ∈ {0, 1} and that for a fixed j there exists a unique i such that r ij (q, ) = 1. This family of random variables will be used to describe the dispatching decisions in the systems. Specifically, if q and are the occupancy state and the threshold, respectively, when the j th task arrives to the system, then this task is sent to a server pool with exactly i − 1 tasks if and only if r ij (q, ) = 1; this coincides with the dispatching rule described in Section 2.1. We postulate that and D n (q, ), respectively. The i = 0 coordinates of these two counting processes are identically zero, whereas the other coordinates are defined as follows: for all i ≥ 1.
Consider the following set of implicit equations: Recall that q n (0, B + 1) = 0 for all n almost surely, as indicated in Section 3. In particular, the initial number of tasks in the system is finite with probability one. Using this property, and proceeding by forward induction on the jumps of the driving Poisson processes, it is possible to prove that there exists a set of probability one Γ 0 with the next property. The sample paths of the occupancy states and the thresholds are defined by extending these solutions to Ω, setting q n (ω) ≡ 0 and n (ω) ≡ 0 for all ω / ∈ Γ 0 . In addition, we define A n := A n (q n , n ), D n := D n (q n , n ) and h n := n + ∆, and we rewrite the implicit equations as follows: The above construction endows the processes q n and n with the intended statistical behavior. Indeed, the time-inhomogeneous Poisson process N λ n has instantaneous rate nλ(t) for each n, and the random functions r ij apply the dispatching rule described in Section 2.1 when the occupancy state and threshold are passed as arguments. In particular, A n (i) models the arrivals to server pools with precisely i−1 tasks in a system with n server pools. In addition, the instantaneous intensity of D n (i) is given by nµi[q n (t, i)−q n (t, i+1)], which is equal to the total number of tasks in server pools with exactly i tasks times the rate at which tasks are executed. Therefore, D n (i) models the departures from server pools with precisely i tasks, and thus (9a) corresponds to the evolution of the occupancy state in a system with n server pools. Furthermore, the refined learning scheme, for adjusting the threshold, is captured by (9b).

Relative compactness of sample paths
Below we state the relative compactness result mentioned at the start of Section 5.1, deferring the proof to Appendix A.

Proposition 4.
There exists a set of probability one Γ T ⊂ Γ 0 with the following property. The sequences {A n (ω) : n ≥ 1}, {D n (ω) : n ≥ 1} and {q n (ω) : n ≥ 1} are relatively compact for each ω ∈ Γ T and have the property that the limit of every convergent subsequence is a function with Lipschitz coordinates, also for each ω ∈ Γ T .
Essentially, the proof of the above proposition relies on the decomposition (9a) and the fact that the coordinates of A n and D n have intensities which are O(n). In particular, observe that the proof provided in Appendix A does not depend on the learning scheme used for adjusting the threshold over time, or even on the rule used to take the dispatching decisions.

Asymptotic dynamical properties
In this section we prove asymptotic dynamical properties concerning the total and tail mass processes. These were defined in (7) as q n (i) and v n (j) := ∞ i=j q n (i) for all j ≥ 1, respectively. The total mass process u n corresponds to the total number of tasks in the system, normalized by the number of server pools, as the diagram of Figure 1 suggests. The tail mass process v n (j) represents the total number of tasks located in column j of the latter diagram, or further to the right, also normalized by the number of server pools. Section 5.2.1 contains some preliminary results, and the asymptotic dynamical properties are established in Section 5.2.2.

Set of nice sample paths
We begin by introducing a set of probability one consisting of sample paths which are well-behaved, in a sense that will be made clear below. For this purpose, we need two technical lemmas. The first one establishes asymptotic upper bounds, which are uniform over time, for the number of tasks sharing a server pool and the threshold. Proof. First let us prove that there exists k ∈ N such that q n (t, k∆+1) = 0 for all t ∈ [0, T ] and all sufficiently large n with probability one. Note that the set m≥1 n≥m {ω ∈ Ω : q n (ω, t, k∆ + 1) > 0 for some t ∈ [0, T ]} is measurable for each k ∈ N since the occupancy processes have right-continuous sample paths. We must establish that this set has probability zero for some k. Consider the sets E k n := {ω ∈ Ω : q n (ω, t, k∆ + 1) > 0 for some t ∈ [0, T ]} .
By the Borel-Cantelli lemma, it is enough to prove that there exists k ∈ N such that Recall from Section 3 that q n (0, B + 1) = 0 and n (0) ≤ B for all n with probability one, and suppose that k > B/∆. Observe that q n (k∆ + 1) remains equal to zero while the threshold is smaller than or equal to (k − 1)∆. Furthermore, the threshold remains strictly smaller than k∆ while nq n (k∆) < n − 1. Therefore, The total number of tasks in the system at time t is upper bounded by N λ n (t) + nu n (0), the number of arrivals during [0, t] plus the initial number of tasks in the system. Hence, The second inequality follows from q n (0, B + 1) = 0, which implies that u n (0) ≤ B. Applying a Chernoff bound, we conclude that Set k such that A (e − 1) + B − k∆ < 0, then (10) holds, so q n (t, k∆ + 1) = 0 for all t ∈ [0, T ] and all large enough n with probability one. Note that q n (t, (k + 1)∆) = 0 for all t ∈ [0, T ] implies that the threshold can never reach (k + 1)∆ during [0, T ], so we may define B T := (k + 1)∆.
The following lemma is a refinement of the functional strong law of large numbers for the Poisson process; a proof is provided in Appendix B. The last two lemmas are used to prove the next proposition, which defines a set of probability one consisting of well-behaved sample paths.  Proof. By Proposition 19 of Appendix A, we have (12a) on Γ T . Also, it follows from Lemma 5 that there exists a positive constant B T ∈ N such that property (a) holds on a subset of Γ T which has probability one. By Lemma 6, this subset can be chosen so that (12b) and (12c) hold.

Properties of total and tail mass processes
Below we prove the aforementioned asymptotic dynamical properties of the total and tail mass processes, starting with the following proposition; the proof is provided in Appendix B. As mentioned in Section 3, the last proposition is simply the functional law of large numbers for an infinite-server system. However, it is not straightforward that this law holds with probability one under the coupled construction of sample paths adopted in Section 5.1.1; Proposition 8 establishes this fact. The next proposition provides an asymptotic upper bound for certain tail mass processes, under specific conditions concerning the occupancy and threshold processes. Proposition 9. Suppose that the following two conditions hold for a given ω ∈ Ω T and a given increasing sequence K of natural numbers.
(a) {q k (ω) : k ∈ K} converges to a function q ∈ D R N[0, T ] with respect to .
Then the coordinate functions q(i) are differentiable on (t 0 , t 1 ) for all i > j, and they satisfẏ Proof. It follows from (b) that A k (t, i) = A k (t 0 , i) for all t ∈ [t 0 , t 1 ] and all i > j, because tasks are only assigned to server pools with strictly fewer than j tasks during [t 0 , t 1 ] in all the systems corresponding to the sequence K. Thus, (9a) implies that for all t ∈ [t 0 , t 1 ] and all i > j. (13) As indicated in Appendix A, convergence of a sequence of functions with respect to implies uniform convergence over [0, T ] of the coordinate functions. By (a) and (12c), By Proposition 7, there exists B T > 0 such that q k (i) is identically zero for all i > B T and all large enough k ∈ K; this implies, in particular, that q(i) is identically zero as well for all i > B T . Therefore, we conclude from (a), (13) and (14) that Proceeding by backward induction, starting from i = B T , it is possible to conclude that q(i) is differentiable on (t 0 , t 1 ) for all j < i ≤ B T , and such thatq(i) = −µi [q(i) − q(i + 1)]. As observed, q(i) is identically zero for all i > B T , so this proves the first claim of the proposition. The convergence of {v k (j + 1) : k ∈ K} follows from (a) and the observation that q k (i) is identically zero for all i > B T and all large enough k ∈ K. Moreover, v(j +1) = B T i=j+1 q(i), and in particular, v(j + 1) is differentiable on (t 0 , t 1 ). It follows thaṫ Note that q(j + 1) is a non-negative function, thus The last inequality is strict since j ≥ 1.

Evolution of the threshold
In this section we prove Theorem 2. For this purpose, we fix an integer m ≥ 0 and an interval [a, b] ⊂ [0, T ]. As in the statement of Theorem 2, we assume that the control parameters α n satisfy (5), that the offered load is (m, ∆)-bounded on [a, b] and that the length of this interval is strictly larger than σ(a, b, m, ∆).
The first step towards establishing Theorem 2 is to prove that all large enough constants σ have the next property: the threshold is smaller than or equal to m∆ along the interval [a+σ, b] in all large enough systems with probability one. This step is carried out in Section 5.3.1, which also provides an upper bound for the infimum of the constants σ having the latter property, and in addition establishes (6c). The next step towards proving Theorem 2 is to demonstrate that σ > σ(a, b, m, ∆) implies that the threshold is in fact equal to m∆ along the interval [a + σ, b] in all sufficiently large systems with probability one, which is done in Section 5.3.2; the latter statement corresponds to (6a), and (6b) is also proven in Section 5.3.2.

Preliminary results
As mentioned above, we begin by establishing that there exists a constant σ such that the threshold is upper bounded by m∆ along the interval [a + σ, b] in all sufficiently large systems with probability one. This is done in the following proposition, which also provides an upper bound for the infimum of the constants σ having the latter property.

Proposition 10. Suppose the hypotheses of Theorem 2 hold, and fix a constant
with ρ max as in Definition 1. For each ω ∈ Ω T , there exists n bd (ω) such that n (ω, t) ≤ m∆ for all t ∈ [a + σ, b] and all n ≥ n bd (ω).
Proof. It follows from (4) that the next inequality holds for all t ∈ [a, b]: We fix an arbitrary ω ∈ Ω T , which is omitted from the notation for brevity. By Definition 1, we may consider positive constants δ and ε such that Note that ε can be expressed in terms of δ. In other words, if δ is given, then ε is determined. For each positive δ, satisfying the previous inequality, and each θ ∈ (0, 1), we define Observe that σ bd (a, b, m, ∆) is the infimum of δ + σ(δ, θ). Therefore, we can choose δ and θ such that σ > δ + σ(δ, θ). Below we assume that δ and θ enforce the latter condition. By Proposition 7 and (5), there exists n 0 such that the next statements hold if n ≥ n 0 .
(ii) N λ n has at least B T − m∆ jumps on each subinterval of [0, T ] of length at least δ. Suppose that j ≥ (m + 1)∆ and n ≥ n bd . If t ∈ [a, b], then The first inequality follows from the fact that q n (t) ∈ Q is a non-increasing sequence, and the last inequality follows from (15). If t ∈ [a + σ(δ, θ), b], then t − a ≥ σ(δ, θ) and thus (δ, θ), b] and n ≥ n bd .
For each n ≥ n bd , the last inequality implies that threshold increases cannot occur beyond m∆ along the interval [a + σ(δ, θ), b]. Moreover, the second to last inequality implies the next property: if n (t) ≥ (m + 1)∆ for some t ∈ [a + σ(δ, θ), b] and some n ≥ n bd , then the threshold decreases with each arrival until it reaches m∆, which occurs in at most δ units of time by (i) and (ii). These two observations imply that n (t) ≤ m∆ for all t ∈ [a + σ, b] and all n ≥ n bd ; recall that σ > δ + σ(δ, θ) and that the threshold takes values on the multiples of ∆.
The next result is a consequence of Propositions 9 and 10, and it implies (6c). Proof. We fix an arbitrary ω ∈ Ω T . For brevity, ω is omitted from the notation and we set r := (m + 1)∆. Suppose that the statement of the corollary does not hold. Then there exist ε > 0 and an increasing sequence of natural numbers K such that sup t∈ [a+σ,b] v k (t, r + 1)e µ[t−(a+σ)] > u(a + σ) + ε for all k ∈ K.
By Proposition 4, we may assume without loss of generality that {q k : k ∈ K} has a limit with respect to ; this may require to replace K by a subsequence. Furthermore, by Proposition 10 we may also assume that k (t) ≤ m∆ for all t ∈ [a + σ, b] and all k ∈ K.
Note that k (t) ≤ m∆ for all t ∈ [a + σ, b] implies that q k (t, r) ≤ 1 − 1/k for all t ∈ [a + σ, b]. Indeed, no arrival can result in q k (r) going from 1 − 1/k to one because such an arrival would also result in the threshold increasing to r = (m + 1)∆. Consequently, Proposition 9 holds with j = r along the interval [a + σ, b], and in particular, there exists a function v(r + 1) such that
This leads to a contradiction, so the statement of the corollary must hold.

Proof of Theorem 2
We have already established that σ > σ bd (a, b, m, ∆) implies that the threshold is upper bounded by m∆ along the interval [a + σ, b] in all large enough systems with probability one, and we have shown that (6c) holds. It remains to prove that σ > σ(a, b, m, ∆) implies that the threshold is in fact equal to m∆ in all sufficiently large systems with probability one, and to establish (6b). To this end, we introduce the processes δ n (j), defined as for all t ∈ [0, T ] and all j ≥ 1. Observe that (12b) and (12c) imply that The following two technical lemmas make use of the latter definition and property.
Then the next inequality holds for all t ∈ [τ k,1 (ω), τ k,2 (ω)] and all k ∈ K: Proof. We omit ω from the notation for brevity. It follows from (9a) that Note that all tasks are assigned to server pools with at most j − 1 tasks during (τ k,1 , t], so the first term on the right-hand side can be expressed as follows: Using the process δ k (j) defined in (16), we may write The last inequality follows from the fact that This completes the proof.
For all sufficiently large k ∈ K and all γ ∈ [0, 1/2), we have Proof. We omit ω from the notation for brevity, and we assume that j > 0 since otherwise the statement holds trivially. Fix γ ∈ [0, 1/2) and consider the random times defined as Observe that k (t) = j for all t ∈ [ζ k,1 , ζ k,2 ∧ ξ k,2 ] and all k ∈ K since k (ζ k,1 ) = j, the threshold is upper bounded by j along [ζ k,1 , ζ k,2 ] and the threshold can only decrease if q k (j) ≤ α k right before an arrival epoch. Thus, it suffices to demonstrate that ξ k,2 ≥ ζ k,2 for all large enough k ∈ K.
In order to prove this, we show that and all large enough k ∈ K.
To prove (18), let τ k,1 := ζ k,2 ∧ ξ k,1 and τ k, For the first inequality, note that the processes q k (i) have jumps of size 1/k. The second inequality follows from Lemma 12, since k (t) = j and q k (t, j) < 1 for all t ∈ [τ k,1 , τ k,2 ) and all k ∈ K. For the third inequality, recall that the offered load is (m, We conclude from (5) and (17) that If 1 − α k is replaced by k −γ on the left-hand side, then the limit is also equal to zero by (17). Consequently, for all sufficiently large k ∈ K, we have This implies (18) for all large enough k ∈ K; (18) holds trivially when ξ k,1 > ζ k,2 .
Proof of Theorem 2. As in the proof of Proposition 10, it follows from (4) that the following inequality holds for all times a ≤ a + σ 0 ≤ t ≤ b: In this equation, ρ min is as in Definition 1. Fix σ(a, b, m, ∆) < σ < b−a as in the statement of the theorem, and for each 0 < ε < ρ min − m∆ let Recall the definitions of σ(a, b, m, ∆) and σ bd (a, b, m, ∆) provided in Definition 1 and Proposition 10, respectively. We fix σ 0 > σ bd (a, b, m, ∆) and ε > 0 such that σ = σ 0 + σ(ε) + ε. This is possible since σ bd (a, b, m, ∆) + σ(ε) + ε decreases to σ(a, b, m, ∆) as ε goes to zero. We fix some ω ∈ Ω T , which is omitted from the notation, and we define The proofs of (6a) and (6b) will be completed if we show that ξ n ≤ a + σ for all large enough n. Indeed, if this claim is established, then (6a) and (6b) become a consequence of Lemma 13 with j := m∆, ζ n,1 := ξ n and ζ n,2 := b. In order to see why, observe that n (t) ≤ m∆ for all t ∈ [a + σ 0 , b] and all sufficiently large n by Proposition 10, so the hypothesis of Lemma 13 holds.
To prove that ξ n ≤ a + σ for all large enough n, we will establish that lim sup If m = 0, then ξ n = a + σ 0 for all n and the above inequality holds. Thus, we assume that m ≥ 1. In order to prove (20), we assume that (20) does not hold and we arrive to a contradiction. Suppose then that (20) is false. This implies that there exists an increasing sequence K of natural numbers such that ξ k > a+σ 0 +σ(ε) for all k ∈ K. Furthermore, by Propositions 4 and 10, this sequence may be chosen so that the next two properties hold.
(i) {q k : k ∈ K} converges to some function q ∈ D R N[0, T ] with respect to .
Let v(m∆+1) be the function defined in the statement of Proposition 9, as the uniform limit of the processes v k (m∆ + 1) over the interval [0, T ]. It follows from Propositions 8 and 9 that sup for all sufficiently large k ∈ K. For each of these k ∈ K, we have The third inequality follows from (19) and Proposition 9, and the last equality follows from the definition of σ(ε). This is a contradiction, so we conclude that (20) holds, and this establishes (6a) and (6b). Recall that (6c) holds by Corollary 11, thus the proof of Theorem 2 is complete.

Coxian service times
In this section we prove Theorem 3. For this purpose, we assume that the basic learning scheme is used, that λ is constant over time and that service times are Coxian distributed. Below we proceed as indicated in Section 4.2.

Relative compactness of occupancy processes
In Section 6.1.1 we construct the sample paths of the occupancy states s n and the thresholds n on a common probability space for all n. The sample paths of the occupancy states lie in the space D R N r [0, T ] of all càdlàg functions on [0, T ] with values in R N r , equipped with the topology of uniform convergence. In Section 6.1.2 we establish that the sequence {s n : n ≥ 1} is almost surely relatively compact.
Before we continue, we introduce some notation and we define a metric on D R N r [0, T ]. A server pool is said to have occupancy profile ν = (ν(1), . . . , ν(r)) ∈ N r if it has exactly ν(m) tasks in phase m for each m. Fix an enumeration {ν i : i ≥ 0} of N r and consider the metric d defined by for all x, y ∈ R N r .
The latter metric generates the product topology of R N r . The topology of uniform convergence equipped to D R N r [0, T ] is generated by the following metric: Note that the previous definitions generalize the ones introduced in Section 5.1, where r = 1. Also, the specific choice of the enumeration {ν i : i ≥ 0} is not particularly relevant since convergence with respect to is equivalent to uniform convergence of the coordinate functions.

Coupled construction of sample paths
The sample paths of the occupancy states and thresholds are defined as deterministic functions of the following stochastic primitives. satisfying the assumptions introduced in Section 3.
Consider the product probability space of (Ω D , F D , P D ), (Ω S , F S , P S ) and (Ω I , F I , P I ); denote its completion by (Ω, F, P). The occupancy states and thresholds are defined on the latter space from a set of equations involving the stochastic primitives. In order to write these equations, we introduce some notation which is analogous to that of Section 5.1.1.
For each occupancy state s ∈ S, we define the intervals These intervals form a partition of [0, 1) such that the length of I i (s) is the fraction of server pools with occupancy profile ν i . Let Σ i := r m=1 ν i (m) denote the total number of tasks in a server pool with profile ν i , and define q ∈ Q as in (1), but in terms of s. If q(j) < 1, then we let These intervals yield another partition of [0, 1). In this case, the length of J i (s, j) is the fraction of server pools with occupancy profile ν i but only among those server pools with at most j − 1 tasks in total; we adopt the convention J i (s, j) := ∅ for all i ≥ 0 when q(j) = 1.
If ∈ N represents the threshold and q( ) < 1, then the length of J i (s, ) is equal to the probability of picking a server pool with occupancy profile ν i uniformly at random among those with strictly less than tasks in total. Similarly, if h := + ∆, q( ) = 1 and q(h) < 1, then the length of J i (s, h) is equal to the probability of picking a server pool with occupancy profile ν i uniformly at random among those with at least and strictly fewer than h tasks in total. We define for all i, j ≥ 1.
Note that r ij (s, ) ∈ {0, 1} and that for a fixed j there exists a unique i such that r ij (s, ) = 1. This family of random variables will be used to describe the dispatching decisions, as in Section 5.1.1. Specifically, if s and are the occupancy state and the threshold, respectively, when the j th incoming task arrives, then this task is sent to a server pool with occupancy profile ν i if and only if r ij (s, ) = 1; this coincides with the dispatching rule described in Section 2.1. We postulate that N λ n (t) := N (0,0,0) (nλt) is the number of tasks that arrive to the system with n server pools during the interval [0, t], we denote the jump times of N λ n by {τ n,k : k ≥ 1} and we define τ n,0 := 0. Note that N λ n is a Poisson process with intensity nλ, as required by our model.
As in Section 5.1.1, it is possible to construct càdlàg processes s n and n such that the following equations hold on a set of probability one Γ 0 for all n and all t ∈ [0, T ]. s n (t) = s n (0) +Ā n (t) +P n (t) +D n (t), (21a) In the above equations, u n is defined in terms of q n as in (7), whereas q n is defined in terms of s n using (1). In addition, the processesĀ n ,P n andD n satisfy: Here {e 1 , . . . , e r } denotes the canonical basis of R r . Also, A n (t, ν) and P n (t, ν, k) are defined as zero if the vector ν has a negative component. A coupled construction of all the above processes can be performed by forward induction on the jumps of the driving Poisson processes since the initial number of tasks in the system is almost surely finite for all n; this precludes an infinite number of events in finite time with probability one. The above construction endows the processes s n and n with the intended statistical behavior. In particular, the process A n (ν i ) counts the arrivals to server pools with occupancy profile ν i , the process P n (ν i , m) counts transitions from phase m to phase m + 1 of tasks in server pools with profile ν i , and the process D n (ν i , m) counts departures of tasks in phase m from server pools with occupancy profile ν i . Observe that (21a) is equivalent to (9a) if r = 1, whereas (21b) now corresponds to the basic learning scheme.

Relative compactness of sample paths
Below we state the relative compactness result mentioned at the start of Section 6.1. The proof essentially relies on the decomposition (21a) and the fact that the coordinate processes A n (ν), P n (ν, m) and D n (ν, m) have O(n) intensities for all ν and m. In particular, the proof is analogous to that of Proposition 4 and is thus omitted.

Proposition 14.
There exists a set of probability one Γ T ⊂ Γ 0 with the next property. The sequences {A n (ω) : n ≥ 1}, {P n (ω, m) : n ≥ 1}, {D n (ω, m) : n ≥ 1} and {s n (ω) : n ≥ 1} are relatively compact with respect to for each m and each ω ∈ Γ T . In addition, they have the property that the limit of every convergent subsequence is a function with Lipschitz coordinates.

Asymptotic dynamical properties
In this section we establish a few asymptotic dynamical properties concerning the following four processes: q n (i) and z n := ∞ i= ρ ∆ +∆+1 q n (i). (22) At any given time, the server pools can be arranged as in the diagram of Figure 1, ignoring the phases of tasks. The mass processes x n , y n and z n correspond to the number of tasks in the columns of the latter diagram covered by the respective summations in (22). The total mass process u n , considered in Section 5.2 as well, represents the total number of tasks in the system, normalized by n, or equivalently, the total mass in the diagram of Figure 1.

Set of nice sample paths
As in Section 5.2.1, we begin by introducing a set of probability one consisting of well-behaved sample paths. For this purpose, we need the following result, counterpart of Lemma 5; the proof is provided in Appendix B. Proof. It follows from Lemma 15 that there exists a positive B T ∈ N such that property (a) holds on a set Ω T ⊂ Γ T of probability one. By (2) and the strong law of large numbers for the Poisson process, the latter set can be chosen so that property (b) holds as well.

Properties of mass processes
Next we prove the asymptotic dynamical properties of the mass processes. The proof of the next proposition provided in Appendix B. While the above law of large numbers is known to hold weakly, and even for more general service times, it is not straightforward that it holds with probability one under the coupled construction of sample paths of Section 6.1.1; Proposition 17 proves this fact.
The following proposition pertains to the dynamics of the mass processes x n , y n and z n along intervals of time where n = := ρ ∆ . Before stating it formally, we provide some intuition using the diagram of Figure 1. Namely, at any given time, the server pools, represented by the rows in said diagram, can be arranged monotonically with respect to the total number of tasks that they currently have. The tasks on each row may be arranged arbitrarily, ignoring the phase in which they are, but regardless of how tasks are arranged within each row, every task leaves the system at a rate that is lower bounded by If n = , then at least one of the first + ∆ columns of the diagram is not full, which implies that new tasks are always dispatched to one of these columns. Similarly, tasks are sent to one of the first columns if x n < = n . Recall that z n represents the number of tasks outside of the first + ∆ columns, hence z n decreases at a rate larger than ηz n while n = . Likewise, y n decreases at a rate that is lower bounded by η [y n − ∆q n ( + ∆ + 1)] ≥ η (y n − ∆z n ) while x n < = n .
Proposition 18. Let us fix ω ∈ Ω T and an increasing sequence K of natural numbers, in the sequel we omit ω for brevity. Suppose that the following conditions hold.
(a) The sequence {s k : k ∈ K} converges to a function s ∈ D R N r [0, T ] with respect to .
Then the sequences {x k : k ∈ K}, {y k : k ∈ K} and {z k : k ∈ K} converge to absolutely continuous functions x, y and z, respectively, uniformly over [0, T ]. Also, the next properties hold for t ∈ (t 0 , t 1 ).
Proof. The uniform convergence of the sequences {x k : k ∈ K}, {y k : k ∈ K} and {z k : k ∈ K} follows from (a) and (23); the latter equation implies that s k (ν i ) is identically zero if Σ i > B T for all large enough k ∈ K. Furthermore, it follows from Proposition 14 that s(ν i ) is Lipschitz continuous for all i, therefore x, y and z are absolutely continuous, in fact Lipschitz.
Assume that x(t) < ρ ∆ for some t ∈ (t 0 , t 1 ) and write := ρ ∆ for brevity. Since x is continuous and x k converges uniformly to x, there exist a < t < b such that Below we establish (i) by proving that Suppose τ and k are such that (23) and (25) hold. It follows from (1) and (21a) that For the last identity, first observe that A k (τ, ν i ) − A k (a, ν i ) = 0 if Σ i ≥ since q k ( ) < 1 along the interval [a, b] by (25). In addition, note that As a result, it is possible to write Here the last identity in (26a) uses (23), which implies that D k (σ) = 0 if σ > B T . Using (24c), we conclude that Therefore, taking limits on both sides of (26a), we obtain For the last inequality, note that Clearly, the inequality between the extremes also holds if s k , y k and z k are replaced by their limits s, y and z. This completes the proof (i). The proof of (ii) follows from similar arguments. Specifically, k = implies that u k < + ∆ and therefore q k ( + ∆) < 1. Consequently, A k (b, ν i ) − A k (a, ν i ) = 0 for all t 0 ≤ a < b ≤ t 1 and k ∈ K if Σ i ≥ + ∆. Then one may proceed as in the proof of (i).

Proof of Theorem 3
In this section we use the asymptotic dynamical properties of the mass processes to complete the proof of Theorem 3.
Proof of Theorem 3. Let := ρ ∆ for brevity, and choose ε > 0 such that + 3ε < ρ < + ∆. It follows from (3) that u is differentiable. Furthermore, u(t) → ρ andu(t) → 0 as t → ∞. Therefore, there exists a time t 0 ≥ 0 such that The previous statement is trivial if t 0 > T , but if T is large enough, then the statement holds for some t 0 < T . Let us fix some ω ∈ Ω T , which we omit from the notation for brevity. Proposition 17 states that u n converges uniformly over [0, T ] to u. Therefore, < u n (t) < + ∆, and thus n (t) = , for all t ∈ [t 0 , T ] and all large enough n. In particular, this proves (8a).
Fix an increasing sequence of natural numbers. By Proposition 14, there exists a subsequence K such that {s k : k ∈ K} converges to a function s ∈ D R N r [0, T ] with respect to . Also, it follows from Proposition 18 that {x k : k ∈ K}, {y k : k ∈ K} and {z k : k ∈ K} converge uniformly over [0, T ] to absolutely continuous functions x, y and z, respectively. Moreover, we have by property (ii) of Proposition 18. In particular, Suppose that x(t) < for some t ∈ (t 1 , T ), then the next statements hold almost everywhere on a neighborhood of t. Since y = u − x − z > u − − z, Proposition 18 and (27) imply thatẏ As a result, we haveẋ =u −ẏ −ż ≥u −ẏ > ηε, where the last inequality uses (27). Thus, Here σ is an upper bound for the time until x reaches if x(t 1 ) < , computed froṁ x > ηε. Also, x cannot decrease once it reaches sinceẋ > ηε almost everywhere on {t ∈ (t 1 , T ) : We conclude from (28) and (29) that (8b) and (8c) hold with c = u(t 0 ) along the sequence K. In fact, (8b) and (8c) hold without restricting the corresponding limits to the sequence K since t 0 , t 1 , σ and c do not depend on K, they only depend on ρ, ∆, η, ε and u. Specifically, if (8b) and (8c) do not hold, then one may construct an increasing sequence of natural numbers which violates the latter equations. Nevertheless, reasoning as above, one obtains a subsequence K along which (8b) and (8c) hold, which leads to a contradiction.
We now establish the relative compactness of the occupancy state sample paths.
Proof of Proposition 4. As above, we fix an arbitrary ω ∈ Γ T and we omit it from the notation for brevity. It follows from (9a) and (30a) that it is enough to establish the result for the sequences such that the limits of their convergent subsequences are Lipschitz; the sufficiency was explained above. We fix i ≥ 0 and we do this just for {A n (i) : n ≥ 1} since the same arguments can be used if A n is replaced by D n . Let M i and {ε i,n : n ≥ 1} be as in the statement of Lemma 22. It follows from Lemma 21 that for each n there exists x n (i) ∈ L M i such that ||A n (i) − x n (i)|| T ≤ 4ε i,n . Recall that L M i is compact, thus every increasing sequence of natural numbers has a subsequence K such that {x k (i) : k ∈ K} converges to some function x ∈ L M i . Also, note that lim sup the limits are taken along K. This proves that every subsequence of {A n (i) : n ≥ 1} has a further subsequence which converges to a Lipschitz function.

Appendix B Auxiliary results
Below are the proofs of some auxiliary results.
For each k ≥ 1, the process [N (t) − t] 2k is a submartingale. It follows from Doob's inequality that where µ i (x) denotes the i th central moment of the Poisson distribution with mean x. Below we prove that there exists k such that the sum over n of the expression on the right-hand side is finite.
The central moments µ i (x) are polynomials of x such that for the last identity see [17,Equation (5.22)]. This recursion implies that µ 2i and µ 2i+1 have degree i. Therefore, there exist coefficients a k,j such that µ 2k (nT ) = k j=0 a k,j n j .
Fix k such that k − 2kγ > 1 and observe that To conclude, note that 2k − j − 2kγ ≥ k − 2kγ > 1 for all 0 ≤ j ≤ k. Hence, This completes the proof.
Proof of Proposition 8. We fix an arbitrary ω ∈ Ω T , which is omitted from the notation for brevity. It follows from Proposition 7 that q n (i), A n (i) and D n (i) are identically zero for all large enough n and all i > B T := B T + ∆ 1 . Consequently, we have D n (t, i) for t ∈ [0, T ] and large enough n.
To prove the proposition, it suffices to establish that every increasing sequence of natural numbers has a subsequence K such that {u k : k ∈ K} converges uniformly over [0, T ] to u. For this purpose, we first observe that (12a) and (12b) imply that The unique solution to this integral equation is u. Therefore, the sequence {u n : n ≥ 1} converges uniformly over [0, T ] to u.
Proof of Lemma 15. It suffices to establish that there exists k ∈ N such that u n (t) < (k + 1)∆ for all t ∈ [0, T ] and all sufficiently large n. Specifically, if u n (t) < (k + 1)∆, then it follows from (21b) that n (t) ≤ k∆, and clearly q n (t, (k + 1)∆) < 1. In particular, incoming tasks are never sent to server pools with (k + 1)∆ tasks or more during [0, T ]. Also, recall from Section 3 that q n (0, B + 1) = 0 for all n with probability one, so we can take B T = max{B, (k + 1)∆}. Therefore, it is enough to show that The total number of tasks in the system at time t is upper bounded by N λ n (t) + nu n (0), the number of arrivals during [0, t] plus the initial number of tasks in the system. Hence, The second inequality follows from q n (0, B + 1) = 0, which implies that u n (0) ≤ B.
Applying a Chernoff bound, we conclude that of (32b) follows from the following arithmetic computations: For the second and last identities, note that P n (ν, j) is identically zero if ν(j) = 0. For all sufficiently large n, the i th term in each of the summations appearing in (32a) and (32b) is zero if Σ i > B T , a condition that is met by all but finitely many i. Therefore, as in the proof of Proposition 8, it follows from (24a), (24b) and (24c) that every increasing sequence of natural numbers has a subsequence K such that {(u 1 k , . . . , u r k ) : k ∈ K} converges uniformly over [0, T ] to a function (u 1 , . . . , u r ) that solveṡ This system of differential equations has a unique solution with initial conditions u m (0) defined as in the statement of the proposition. Since u m n (0) converges to u m (0) for all m by (24a), it follows that {(u 1 n , . . . , u r n ) : n ≥ 1} converges uniformly over [0, T ] to this unique solution.
The above arguments imply that {u n : n ≥ 1} converges uniformly over [0, T ] on Ω T . By the weak law of large numbers of an infinite-server system, {u n : n ≥ 1} converges weakly to the deterministic process defined by (3). Since Ω T has probability one, it follows that {u n (ω) : n ≥ 1} converges uniformly over [0, T ] to the function defined by (3) for all ω ∈ Ω T .