# Asynchronous Load Balancing and Auto-scaling:
Mean-Field Limit and Optimal Design

J. Anselmi
Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG, 38000 Grenoble, France.
<EMAIL_ADDRESS>
###### Abstract

We introduce a Markovian framework for load balancing where classical
algorithms such as Power-of-$d$ are combined with asynchronous auto-scaling
features. These allow the net service capacity to scale up or down in response
to the current load within the same timescale of job dynamics. This is
inspired by serverless frameworks such as Knative, used among others by Google
Cloud Run, where servers are software functions that can be flexibly
instantiated in milliseconds according to user-defined scaling rules. In this
context, load balancing and auto-scaling are employed together to optimize
both user-perceived delay performance and energy consumption. In the
literature, these mechanisms are synchronous or rely on a central queue. The
architectural novelty of our work is to consider an _asynchronous_ and
_decentralized_ system, as in Knative, which takes scalability to the next
level.

Under a general assumption on the auto-scaling process, we prove a mean-field
limit theorem that provides an accurate and tractable approximation for the
system dynamics when the mean demand and nominal service capacity grow large
in proportion. We characterize the fixed points of the mean-field limit model
and provide a simple condition telling whether or not all the available
servers need to be turned on to handle the incoming demand. Then, we
investigate how to design optimal auto-scaling rules and find a general
condition able to drive the mean-field dynamics to delay and relative energy
optimality, a situation where the user-perceived delay and the relative energy
wastage induced by idle servers vanish. The proposed optimality condition
suggests to scale up capacity if and only if the mean demand exceeds the
overall rate at which servers become idle-on, i.e., idle and active. This
yields the definition of tractable optimization frameworks to trade off
between energy and performance, which we show as an application of our work.

##  1 Introduction

_Load balancing_ is the process of distributing work units (jobs) across a set
of distributed computational resources (servers). Exogenous jobs join the
system over time through one or more dispatchers, and these route each of them
to one out of $N$ parallel servers for processing. Then, each job leaves the
system upon service completion at its designated server. In large
architectures, each server has its own queue, as this enhances scalability,
and jobs are immediately dispatched upon their arrival. Given the stringent
latency requirements of modern applications, breaches of which can severely
impact revenue, load balancing techniques are usually designed to optimize
delay performance, and popular examples are Power-of-$d$ [24] and Join-the-
Idle-Queue (JIQ) [21]. Closely related to load balancing, _auto-scaling_ is a
term often used in cloud computing to refer to the process of adjusting the
available service capacity automatically in response to the current load [29].
In this context, auto-scaling techniques are meant to control capacity ($N$)
over time to avoid performance degradation, which yields unacceptably large
delays, and overprovisioning of resources, which yields high infrastructure
and energy costs. Google Cloud Run, Amazon Elastic Compute Cloud (EC2),
Microsoft Windows Azure and Oracle Cloud Platform are examples of platforms
that offer auto-scaling and load balancing features. Users of these platforms
deploy their applications defining how the system should scale up resources in
front of an increased load. Modern auto-scaling mechanisms are extremely
reactive in the sense that they adapt capacity relying on fresh observations
of the system state rather than historical data. This especially holds true in
_serverless computing platforms_ , or Function-as-a-Service, which nowadays
provide the convenient solution to deploy any type of application or backend
service [22].

Most of existing performance models for load balancing assume that the
available service capacity, $N$, remains constant over time even when designed
for cloud computing systems [37], i.e., auto-scaling is not taken into
account. Nonetheless, auto-scaling mechanisms are frequently employed in cloud
applications and affect delay performance. This does not mean that classic
load balancing models are inadequate for cloud systems but that they assume
that auto-scaling operates at a much slower timescale than load balancing.
Essentially, this means that jobs do not see any change in the overall
available capacity because they move much faster than servers. This makes
sense if servers are interpreted as physical or even virtual machines because
setup times are of the order of minutes if not longer [15] while in typical
applications hosted in cloud networks job service times are about ten
milliseconds. The large body of literature on load balancing, reviewed in
Section 5, is undoubtedly the proof that this timescale separation assumption
is well accepted by researchers and practitioners. In the context of
serverless computing however, a server is interpreted as a software function
that can be flexibly instantiated in milliseconds [40, 41], i.e., within a
time window that is comparable with the magnitude of job inter-arrival and
service times, and with negligible switching costs. Here, auto-scaling
mechanisms are extremely reactive and the decisions of turning servers on or
off are based on instantaneous observations of the current system state rather
than on the long-run equilibrium behavior or historical data. Therefore, the
timescale separation assumption above becomes arguable [22] because it would
mean to assume that job dynamics achieve stochastic equilibrium between
consecutive changes of $N$, i.e., in milliseconds.

While there is a large body of literature investigating load balancing and
auto-scaling _separately_ , very little has been done when both are applied
jointly within the same timescale. Most of existing works focused on
synchronous and centralized architectures where all servers share a common
queue [15, 22], where synchronous means that scale-up and dispatching
decisions are taken simultaneously. In typical cloud architectures however, no
central queue is maintained as it would affect scalability. A notable
exception is [26], where the authors consider a synchronous but decentralized
architecture where each server has its own queue; see also [16]. In
particular, they synchronize JIQ [21] with a specific auto-scaling strategy
and analyze the resulting performance in the mean-field limit; see Section 5
for further details. Also our framework considers a decentralized architecture
but the key novelty is that load balancing and auto-scaling operate
_asynchronously_ ; the term “asynchronous” is borrowed from the cloud
computing community [22]. In other words, the scale-up and dispatching
processes are decoupled. The asynchronous approach takes scalability to the
next level, for instance because it facilitates the deployment of systems with
a large number of dispatchers [21], and is one principle underlying real
systems such as Knative [1], a Kubernetes-based platform to deploy and manage
modern serverless workloads that is used among others by Google Cloud Run; see
Section 5.1 for further details. In addition, instead of focusing on a
specific load balancing and auto-scaling strategy, our work will provide a
framework to evaluate the impact of a general class of auto-scaling mechanisms
that are up to the platform user to choose or design. This will be one
advantage of our framework with respect to existing works because systems such
as Knative leave the platform user with such control.

###  1.1 Contributions

We propose a Markovian framework for load balancing that includes asynchronous
auto-scaling mechanisms. We refer to this framework as ‘Asynchronous Load
Balancing and Auto-scaling’ (ALBA). ALBA is inspired by and tailored to
serverless frameworks such as Knative; more details in Section 2. Two
mechanisms drive dynamics in ALBA:

  * i)

a _dispatching rule_ , or load balancing rule, which defines how jobs are
dispatched among the set of active servers as they join the system, and

  * ii)

a _scaling rule_ , which defines how the number of active servers scales up
and down over time, possibly as a function of the current system state.

In agreement with existing implementations of serverless frameworks [22], we
assume that dispatching and scaling rules operate _asynchronously_. This
implies that the capacity scaling decisions are not taken at the moments where
jobs arrive and represents a key difference with respect to the existing
approaches in the literature. The main motivation for us to consider an
asynchronous architecture is to provide the first performance model and
analysis of existing serverless platforms such as Knative. Their asynchronous
structure enhances scalability as it facilitates operations in large systems
with several dispatchers. The dispatching rules included in our framework are
Power-of-$d$ and Join-Below-Threshold-$d$ (JBT-$d$) because they are commonly
used in serverless frameworks [2]. A typical scale-down rule turns a server
off only if it remains idle during an expiration window [22, 1]. This rule,
also known as “delay-off” [15], is easy to implement and will be assumed
within the proposed framework. In contrast, we do not impose any particular
structure on scale-up rules as they are often defined by the user of the
serverless platform. Since we have fixed the scale-down rule, in the following
the term “scaling rule” will refer to a scale-up rule.

To the best of our knowledge, ALBA is the first framework that models the
interplay between load balancing and auto-scaling asynchronously. Our key
technical contributions hold under a general assumption on the structure of
scaling rules, see Assumption 2, and are as follows:

  * •

In Theorem 1, we prove a mean-field limit theorem that provides an accurate
deterministic approximation of the sample paths of the finite-dimensional
stochastic process underlying ALBA. This enables analytical tractability and
allows one to study dynamics easily. For instance, given an arbitrary scaling
rule, it may be used to estimate the resulting proportion of active servers,
and thus the energy costs, in both the transient and stationary regimes. We
prove Theorem 1 following the framework developed in [34, 8], which allows us
to handle discontinuities of the drift function of the underlying Markov
chain.

  * •

In Theorem 2, we characterize the _fixed points_ of the mean-field limit in
terms of a set of non-linear equations and provide a simple necessary and
sufficient condition on the scaling rule able to tell whether or not full
service capacity will be needed to handle the incoming demand. Within Power-
of-$d$, roughly speaking, there always exists a unique fixed point if the
scaling rule is “nice”. Within JBT-$d$, however, this is not the case: to
ensure uniqueness within a certain class of interesting scaling rules, it is
necessary that the scaling rule has access to the number of servers containing
exactly one job at any point in time (see Remark 1).

  * •

When the mean demand exceeds the current service capacity, one may expect that
a (strictly) positive scale-up rate at all times is enough to eventually drive
the system to the reversed situation where capacity exceeds, or matches,
demand. In Proposition 1, we provide a counterexample showing that this is not
necessarily the case. In addition, the mean-field dynamics may converge to
different fixed points, depending on the initial conditions. This poses the
challenge of designing scaling rules able to drive dynamics to a desired fixed
point regardless of the initial condition. This is addressed in our next
contribution.

  * •

In Theorem 3, we investigate how to design _fluid optimal_ scaling rules: we
provide a simple general condition ensuring that dynamics of the mean-field
limit converge to “delay and relative energy optimality”. Following [26], this
is a situation where the user-perceived delay and the relative energy wastage
induced by idle servers vanish. We show that delay and relative energy
optimality can only be achieved within JIQ, or equivalently JBT-$0$. The
proposed condition suggests to scale up capacity if and only if the mean
demand exceeds the overall rate at which servers become active _and_ idle.
From a practical standpoint, this requires the central controller to estimate
the mean job arrival rate, which can be done via machine learning techniques,
and to have access to the amount of both initializing, i.e., in set-up mode,
and busy servers containing exactly one job. In existing serverless platforms,
this information is easily accessible and implies only a constant
communication overhead per job as in standard implementations of JIQ (see
Remark 5). In addition, we find that being able to distinguish between busy
servers with one or more jobs is _necessary_ for fluid optimality (see Remark
3). This exposes a structural difference with respect to synchronous
architectures. In fact, delay and relative energy optimality can be achieved
even when JIQ is synchronized with the ad-hoc scale-up rule proposed in [26],
which turns a server on if all active servers are busy upon job arrivals. In
this respect, our contribution is to show that optimality can be achieved
under broader conditions and also within asynchronous architectures. In
existing serverless platforms, we stress that decentralized implementations of
load balancing and auto-scaling mechanisms are indeed asynchronous.

While we show that fluid optimality can only hold within JIQ, we observe that
this may not be the convenient choice within architectures with several
dispatchers. In this case, an exact implementation of JIQ implies an expensive
communication overhead per job and Power-of-$d$ may be the way to go. A
further advantage of Power-of-$d$ with respect to JIQ or JBT-$d$ is that is
_stateless_ , i.e., it does not require the dispatcher(s) to store information
about the server states. This stateless feature is generally required for load
balancing at the network layer [12]. Within these scenarios, the advantage of
the proposed framework is to allow one to be able to design and efficiently
evaluate the performance of different scaling rules. To the best of our
knowledge, our work is the first to provide such ability.

As application of our results, we develop a tractable optimization framework
allowing the system manager to trade off between performance and energy costs.
The main purpose of this optimization framework is to illustrate how the
results presented in this paper can be applied.

###  1.2 Organization

This paper is organized as follows. In Section 2, we build the ALBA framework
and define the structure of the scaling rules considered in our framework. In
Section 3, we present our main results: Theorems 1, 2 and 3. In Section 4, we
show an application example of our results to trade off between energy and
performance. In Section 5, we review the existing literature, and in Section
6, we draw the conclusions. Proofs of our main results are deferred to
Appendix 7.

##  2 Asynchronous Load Balancing and Auto-scaling (ALBA)

In this section, we first describe the main principles at the basis of
Asynchronous Load Balancing and Auto-scaling (ALBA). Since our aim is to
develop a model tailored to serverless computing, we will make several
references to Knative, a popular serverless framework for hosting Function-as-
a-Service processing that is used, among others, by Google Cloud Run. Then, we
propose two performance models for ALBA. The first captures the stochastic
nature of the underlying dynamics while the second is deterministic and will
serve to approximate the complex dynamics induced by the first. The advantage
of the deterministic model is its tractability. Finally, we formalize the
structure of the scaling rules investigated in this paper.

###  2.1 System Description

The proposed framework, ALBA, is composed of a system of $N$ parallel servers,
each with its own queue, that represent the nominal service capacity, i.e.,
the upper limit on the amount of resources that one user can have up and
running at the same time111In Knative, this upper limit is specified by the
max-scale-limit global key.. In the cloud computing community, servers are
also referred to as containers, cloud functions, instances or replicas. Public
serverless computing platforms usually require to specify such limit in order
to ensure service availability for other users. In the following, the terms
servers and queues will be used interchangeably. A server is said _warm_ if
turned on, _cold_ if turned off and _initializing_ if making the transition
from cold to warm. These are the possible server states [41, 23, 22]. An
initializing server performs basic startup operations such as connecting to
database, loading libraries, etc. This is the time to provision a new function
instance. Only warm servers are allowed to receive jobs. When no cold server
exists, we say that the system is saturated; in ALBA, a saturated system is
not representative of the average case. A server is also said _idle-on_ if
warm but not processing any job, and _busy_ if warm and processing some job.
Typically, billing policies charge per number of warm and initializing servers
used per time unit.

Jobs join the system from an exogenous source to receive service. Upon
arrival, each job is dispatched to a warm server according to some dispatching
rule. After dispatching, each job is processed by the selected server
according to the presumed scheduling discipline at that server. After
processing, each job leaves the system.

######  Assumption 1.

Jobs are dispatched to servers according to either Power-of-$d$ or Join-Below-
Threshold-$d$ (JBT-$d$).

We recall that Power-of-$d$ sends an incoming job to the shortest among $d\geq
1$ warm servers selected at random at the moment of its arrival and JBT-$d$
sends an incoming job to a warm server containing no more than $d\geq 0$ jobs
if one exists otherwise to a warm server selected at random. In all cases,
ties are broken randomly. If $d=0$, JBT-$d$ is also known as Join-the-Idle-
Queue (JIQ) [21]. We limit our framework to these types of schemes because
they involve a constant communication overhead per job (in architectures with
a single dispatcher) and because they are commonly used in practice. For
instance, Knative uses Power-of-$2$ if no limit is set on the queue length of
each server and JBT-$d$ if such limit is set to $d$ [2].

Alongside with the above job dynamics, the pools of warm/initializing/cold
servers change over time in the background and in an asynchronous manner.
Precisely, the platform monitors the system state at some epochs that we refer
to as _scaling times_. At such times, a cold server is selected, provided that
one exists, and becomes initializing according to the outcome of some scaling
rule. After some _initialization time_ , or coldstart latency, an initializing
server becomes idle-on. When a server becomes idle-on, it becomes cold after a
scale down delay, or _expiration time_ , if during such time the server
received no job; this scale-down rule is used in several serverless computing
platforms (including Knative) [41, 38] and also in other settings [15].
Therefore, the decision of activating a warm server is taken by a central
controller while the decision of making a server cold is taken by the server
itself. We observe that the number of warm servers fluctuates from 0 to $N$
over time. While in practice it may be possible to set a lower limit on the
number of warm servers, the scale down to zero (or one) servers configuration
is usually the default choice [3].

To a great extent, the scale up rule, the expiration rate and the scaling
times are under the control of the platform user, which may design them in a
way to optimize a trade-off between performance and energy. On the other hand,
several measurements indicate that initialization times are typically one
order of magnitude higher than jobs’ service times in serverless platforms
[22, 41].

###  2.2 Notation

We introduce some notation that will be used throughout the paper. Let
$B\in\mathbb{Z}_{+}\cup\\{+\infty\\}$ be a constant that will denote the
buffer size of each server. We use $\mathbb{I}_{\\{A\\}}$ to denote the
indicator function of $A$. If $a\in\mathbb{R}$ and $A$ denotes an interval,
$\mathbf{1}_{A}^{a}:=\mathbb{I}_{\\{a\in A\\}}$. We also let
$(\cdot)^{+}:=\max\\{\cdot,0\\}$ and $\|\cdot\|$ denotes the $L_{1}$ norm.
Unless specified otherwise, $(i,j)$ ranges over the set
$\\{0,\ldots,B\\}\times\\{0,1,2\\}$ if $B<\infty$ and over
$\mathbb{Z}_{+}\times\\{0,1,2\\}$ otherwise. The process of interest will take
values in

$\mathcal{S}:=\Big{\\{}(x_{i,j}\in\mathbb{R}_{+},\forall(i,j)):\sum_{i,j}x_{i,j}=1\Big{\\}}$
(1)

and our analysis holds under the distance function $d_{w}$ induced by the
weighted $\ell_{2}$ norm $\|\cdot\|_{w}$ on $\mathbb{R}^{\mathbb{Z}_{+}}$
defined by
$\|x-x^{\prime}\|_{w}^{2}:=\sum_{i,j}\frac{|x_{i,j}-x_{i,j}^{\prime}|^{2}}{2^{i+j}}.$
For $x\in\mathcal{S}$, let $y_{i}:=\sum_{k\geq i}x_{i,2}$. We also let
$\mathcal{S}_{1}:=\\{x\in\mathcal{S}:\sum_{i\geq 1}ix_{i,2}<\infty\\}$.

###  2.3 Markov Model

We model the dynamics induced by ALBA in terms of a continuous time Markov
chain. The exogenous arrival process of jobs is assumed to be Poisson with
rate $\lambda N$, with $0<\lambda<1$. Here, we may consider a time-varying
arrival rate of the form $\Lambda(t)N$, where $\Lambda(t)$ is a bounded
positive real-valued function, but we do not consider this case for
simplicity. The processing times, or service times, of jobs are independent
and exponentially distributed random variables with unit mean, and servers
process jobs according to any work-conserving discipline. Upon arrival, each
job is assigned to one of the available warm servers as specified in
Assumption 1. In the extreme case where no warm server exists, we assume that
the job is lost. We also assume that each server can contain at most $B>d$
jobs and a job that is sent to a server with $B$ jobs is rejected. If not
specified otherwise, $B$ is either finite or infinite. At each scaling time, a
cold server is selected uniformly at random, provided that one exists, and
becomes initializing with some probability $g$, referred to as _scaling
probability_ (or rule), that will depend on the system state; in the
Conclusions, we will discuss how our work adapts to the case where a random
number of cold servers is selected at each scaling time. Given that jobs
arrive with a rate proportional to $N$ and only one server can be added at
each scaling time, in our model we let the scaling frequency increase with $N$
as well. As it occurs in Knative, this implies that the number of servers
created in a time window of constant size is proportional to $N$ if within
such window the scaling probability is not zero. We let the inter-scaling,
initialization and expiration times be independent and exponentially
distributed with rate $\alpha N$, $\beta$ and $\gamma$, respectively.

The five sequences of inter-arrival, service, inter-scaling, initialization
and expiration times are also assumed mutually independent.

Let $Q^{N}(t):=(Q_{1}^{N}(t),\ldots,Q_{N}^{N}(t))$ be the vector of queue
lengths at time $t$, including the jobs in service, and let
$S^{N}(t):=(S_{1}^{N}(t),\ldots,S_{N}^{N}(t))$ be the vector of server states.
Specifically, $S_{k}^{N}(t)\in\\{0,1,2\\}$ indicates whether server $k$ is
cold ($S_{k}^{N}(t)=0$), initializing ($S_{k}^{N}(t)=1$) or warm
($S_{k}^{N}(t)=2$) at time $t$. Under the above assumptions, the stochastic
process $(Q^{N}(t),S^{N}(t))$ is a continuous-time Markov chain on state space
$\\{(n,s)\in\\{0,\ldots,B\\}^{N}\times\\{0,1,2\\}^{N}:n_{k}>0\Rightarrow
s_{k}=2,\,\forall k=1,\ldots,N\\}$.

It is convenient to describe dynamics in terms of the process
$X^{N}(t):=(X_{0,0}^{N}(t),X_{0,1}^{N}(t),X_{i,2}^{N}(t):i=0,\ldots,B)$ where

$\displaystyle
X_{i,j}^{N}(t):=\frac{1}{N}\sum_{k=1}^{N}\mathbb{I}_{\\{Q_{k}^{N}(t)=i,S_{k}^{N}(t)=j\\}}$
(2)

is the proportion of servers in state $j$ with $i$ jobs at time $t$. The
process $X^{N}(t)$ is still a Markov chain with values in some set
$\mathcal{S}^{(N)}$ that is a subset of $\mathcal{S}$. Let
$e_{i,j}:=\left(\delta_{i,i^{\prime}}\,\delta_{j,j^{\prime}}\in\\{0,1\\}:i^{\prime}\geq
0,j^{\prime}=0,1,2\right)$ where $\delta_{a,b}$ denotes the Kronecker delta
and let $x:=(x_{i,j})\in\mathcal{S}^{(N)}$ denote a generic state of
$X^{N}(t)$. For conciseness, the Markov chain $X^{N}(t)$ has the following
transitions:

$\begin{array}[]{lll}x\mapsto
x^{\prime}:=x+\frac{1}{N}(e_{i,2}-e_{i-1,2})&{\rm with\leavevmode\nobreak\
rate}\quad\lambda Nf_{i-1}(x)\\\ \vspace{-0.13cm}\\\ x\mapsto
x^{\prime}:=x+\frac{1}{N}(e_{i-1,2}-e_{i,2})&{\rm with\leavevmode\nobreak\
rate}\quad x_{i}N\\\ \vspace{-0.13cm}\\\ x\mapsto
x^{\prime}:=x+\frac{1}{N}(-e_{0,0},e_{0,1})&{\rm with\leavevmode\nobreak\
rate}\quad\alpha Ng\\\ \vspace{-0.13cm}\\\ x\mapsto
x^{\prime}:=x+\frac{1}{N}(-e_{0,1},e_{0,2})&{\rm with\leavevmode\nobreak\
rate}\quad\beta x_{0,1}N\\\ \vspace{-0.13cm}\\\ x\mapsto
x^{\prime}:=x+\frac{1}{N}(e_{0,0}-e_{0,2})&{\rm with\leavevmode\nobreak\
rate}\quad\gamma x_{0,2}N\end{array}$

for all $i=1,\ldots,B$, provided that $x,x^{\prime}\in\mathcal{S}^{(N)}$.
Here, $g:=g(x):\mathcal{S}\to[0,1]$ is the scaling probability, and
$f_{i}(x)$, which depends on the dispatching rule, represents the probability
of assigning an incoming job to a warm server containing exactly $i$ jobs. If
$y_{0}>0$, within Power-of-$d$ we have (assuming that server selections are
with replacement)

$f_{i}(x)=\frac{y_{i}^{d}-y_{i+1}^{d}}{y_{0}^{d}},\quad i=0,\ldots,B-1$ (3)

where $y_{i}:=y_{i}(x):=\sum_{j\geq i}x_{j,2}$, and within JBT-$d$ we have

$f_{i}(x)=\frac{x_{i,2}\,\mathbb{I}_{\\{\sum_{k=0}^{d}x_{k,2}=0\\}}}{y_{0}}+\frac{x_{i,2}\,\mathbb{I}_{\\{\sum_{k=0}^{d}x_{k,2}>0\\}}}{\sum_{k=0}^{d}x_{k,2}}\mathbb{I}_{\\{i\leq
d\\}},\quad i=0,\ldots,B-1$ (4)

where we have taken the convention that $0/0=0$. If $y_{0}=0$, then
$f_{i}(x)=0$ as no warm server exists.

###  2.4 Deterministic Model

We introduce the deterministic (or fluid, mean-field) model for the dynamics
of ALBA.

######  Definition 1.

A continuous function $x(t):\mathbb{R}_{+}\to\mathcal{S}$ is said to be a
_fluid model_ (or fluid solution) if for almost all $t\in[0,\infty)$

$\displaystyle\dot{{x}}_{0,0}$ $\displaystyle=\gamma{x}_{0,2}-\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (5a) $\displaystyle\dot{{x}}_{0,1}$ $\displaystyle=\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\beta{x}_{0,1}+\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (5b) $\displaystyle\dot{x}_{0,2}$ $\displaystyle=x_{1,2}-h_{0}(x)+\beta
x_{0,1}-\gamma x_{0,2}$ (5c) $\displaystyle\dot{x}_{i,2}$
$\displaystyle=x_{i+1,2}\mathbb{I}_{\\{i<B\\}}-x_{i,2}+h_{i-1}(x)-h_{i}(x)\mathbb{I}_{\\{i<B\\}},\quad
i=1,\ldots,B$ (5d)

where $g:=g(x):\mathcal{S}\to[0,1]$, and

$\displaystyle
h_{i}(x)=\left\\{\begin{array}[]{ll}\lambda\frac{y_{i}^{d}-y_{i+1}^{d}}{y_{0}^{d}}&{\rm
if}\,\,y_{0}>0\\\ \min\\{\beta x_{0,1},\lambda\\}&{\rm
otherwise}\end{array}\right.$ (8)

if Power-of-$d$ is applied, and

$\displaystyle
h_{i}(x)=\left\\{\begin{array}[]{ll}\lambda\frac{x_{i,2}}{\sum_{k=0}^{d}x_{k,2}}\mathbb{I}_{\\{i\leq
d\\}}&{\rm if}\,\,y_{0}>0,\,\,\sum_{k=0}^{d}x_{k,2}>0\\\
\left(\beta{x}_{0,1}+{x}_{d+1,2}\mathbb{I}_{\\{i=d\\}}\right)\mathbb{I}_{\\{{x}_{d+1,2}+(d+1)\beta{x}_{0,1}\leq\lambda\\}}&{\rm
if}\,\,y_{0}>0,\,\,\sum_{k=0}^{d}x_{k,2}=0,\,\,i\leq d,\\\
\frac{{x}_{i,2}}{y_{0}}(\lambda-{x}_{d+1,2}-(d+1)\beta{x}_{0,1})^{+}&{\rm
if}\,\,y_{0}>0,\,\,\sum_{k=0}^{d}x_{k,2}=0,\,\,i>d,\\\ \min\\{\beta
x_{0,1},\lambda\\}&{\rm otherwise}\end{array}\right.$ (13)

if JBT-$d$ is applied.

As for $X_{i,j}^{N}(t)$, $x_{i,j}(t)$ is interpreted as the proportion of
servers in state $j$ with $i$ jobs at time $t$.

Let us provide some intuition about the fluid model. First, when a strictly
positive fluid mass of warm server exists, i.e., $y_{0}>0$, the functions
$h_{i}$ are interpreted as the rate at which jobs are assigned to servers with
exactly $i$ jobs. When the amount of fluid of cold servers is strictly
positive, i.e., $x_{0,0}>0$, to some extent these equations may be interpreted
as the conditional expected change, or _drift_ , from state $x$ of the Markov
chain $X^{N}(t)$. In contrast, when $x_{0,0}=0$, there exists a term,
$-\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha g\\}}\gamma{x}_{0,2}$
(see (5a) and (5b)), that still drains the amount of cold servers down. This
is due to warm servers that become cold but immediately turn initializing and
it appears if the scaling rule is ‘greedy enough’, i.e., if the rate at which
new initializing servers can be created is greater than or equal to the rate
at which warm servers go cold. This term is due to fluctuations of order $1/N$
that appear when $X_{0,0}^{N}(t)=0$, which bring discontinuities in the drift
of $X^{N}(t)$, and will come out from the stochastic analysis developed in
Section 7.1.

Now, let us focus on (8) and (13), and let us assume that $y_{0}>0$. In the
case of Power-of-$d$, $h_{i}=\lambda f_{i}$ and $x(t)$ evolves following the
natural dynamics of Power-of-$d$ as in [24], though normalized on the variable
mass of warm servers $y_{0}(t)$. The case of JBT-$d$ is more delicate because
of the discontinuous structure of $f_{i}$ in (4). If a strictly positive
fraction of warm servers with no more than $d$ jobs exist, then $h_{i}=\lambda
f_{i}$ and $x(t)$ evolves following the natural dynamics of JBT-$d$, though
again normalized on a variable number of servers. On the other hand, when
$\sum_{k=0}^{d}x_{k,2}=0$, there is a flow of warm servers with at most $d$
jobs that are created but immediately used for dispatching jobs. Specifically,
there are two factors that come into play here: the first is due to
initializing servers that get warm with exactly $i$ jobs (with rate $\beta
x_{0,1}$), for all $i\leq d$, and the second is due service completions from
servers with exactly $d+1$ jobs (with rate $x_{d+1,2}$). The resulting rate
can not be greater than $\lambda$, the rate where jobs are assigned to
servers, and this justifies the
$\mathbb{I}_{\\{{x}_{d+1,2}+(d+1)\beta{x}_{0,1}\leq\lambda\\}}$ term. Then,
the excess of such rate, $(\lambda-{x}_{d+1,2}-(d+1)\beta{x}_{0,1})^{+}$, is
distributed uniformly over servers with $i>d$ jobs. In Theorem 3, we will show
that such rate is key for the design of fluid optimal scaling rules. Finally,
assume that no warm server exists, i.e., $y_{0}=0$. Here, initializing servers
get idle-on with rate $\beta x_{0,1}$ but all of them are immediately filled
by new arrivals if $\lambda\geq\beta x_{0,1}$, and in this case the mass of
idle-on servers remains zero. Otherwise, $x_{0,2}$ increases with surplus rate
$\beta x_{0,1}-\lambda$.

The existence of a fluid solution started in $x^{(0)}\in\mathcal{S}_{1}$ will
be direct from Theorem 1.

###  2.5 Scaling Rules

The scaling rule $g$ gives the probability to activate a new server at each
scaling time as a function of the system state. The following assumption,
which will hold throughout the paper, provides the structure of the scaling
rules investigated in this paper.

######  Assumption 2.

The scaling rule $g:\mathcal{S}\to[0,1]$ is Lipschitz continuous, and $g(x)>0$
if $x_{0,0}=1$.

The last technical condition is natural and will rule out the existence of
degenerate fixed points. We allow $g(x)$ to be greater than zero even when no
cold server exists, i.e., $x_{0,0}=0$. While this has no impact on the
dynamics of the stochastic model, it does affect the fluid model as there may
exist a flow of idle-on servers that go cold but instantly turn initializing
keeping the proportion of cold servers at zero. This situation can occur if
$\lambda$ is large enough and not only in the transient regime; see Theorem 2.

We propose two scaling rules that satisfy Assumption 2 and that will be used
for illustration purposes.

######  Definition 2.

At each scaling time, if the system state is $x$,

  * •

_Blind_ -$\theta$ activates a new server with probability $g(x)=\theta$,
$\theta\in(0,1]$;

  * •

_Rate-Idle_ activates a new server with probability
$g(x)=\frac{1}{\lambda}(\lambda-\beta x_{0,1}-x_{1,2})^{+}$.

Blind-$\theta$ is oblivious of the system state and thus highly scalable. On
the other hand, within Rate-Idle, the auto-scaler needs to know the amount of
initializing servers, the amount of busy servers with exactly one job at any
point in time and both the job arrival and server initialization rates, which
can be easily estimated. If $\lambda>\beta x_{0,1}+x_{1,2}$, then the
interpretation is that the scaling probability is proportional to the
difference between the mean demand, $\lambda$, and the rate at which servers
become idle-on, $\beta x_{0,1}+x_{1,2}$. Otherwise, the scale-up process is
turned off. Rate-Idle is interesting because it will induce fluid optimality
if combined with JIQ; see Theorem 3.

##  3 Main Results

In this section, we present our main results. First, in Theorem 1 we justify
the use of the deterministic model to approximate the behavior of the
stochastic. Then, we focus on properties of the deterministic model: we
characterize its fixed points in Theorem 2 and investigate the design of
optimal scaling rules in Theorem 3.

###  3.1 Connection between the Fluid and the Markov models

The following result shows that the fluid model can be seen as a first-order
approximation of the sample paths of the stochastic.

######  Theorem 1.

Let $T<\infty$, $x^{(0)}\in\mathcal{S}_{1}$ and assume that
$\|X^{N}(0)-x^{(0)}\|_{w}\to 0$ almost surely. Then, limit points of the
stochastic process $(X^{N}(t))_{t\in[0,T]}$ exist and almost surely satisfy
the conditions that define a fluid solution started at $x^{(0)}$.

###### Proof.

Given in Appendix 7. ∎

The stochastic and the deterministic models have some non-standard aspects
that prevent us to prove Theorem 1 by directly applying Kurtz’s theorem or
similar known results. Specifically, the technical issues are that the
trajectories of the deterministic model may cross or converge, as time
increases, to points of discontinuity of its drift function and that buffer
sizes are allowed to be infinite. To handle these issues, we follow the
general framework in [34, 8] developing ad-hoc arguments specific to the
structure of our problem.

In view of Theorem 1 and since typical and default maximum scale limit values
are 1000 or more, i.e., $N\geq 10^{3}$, we expect that the fluid model $x(t)$
provides an accurate approximation of the average behavior of $X^{N}(t)$. To
support this claim, we now present the results of numerical simulations; see
also Section 4. Figure 1 (left) plots some trajectories of $x(t)$ and
$X^{N}(t)$ when $N=10^{3}$ and $B=10^{2}$ along the coordinates of cold
($x_{0,0}$), initializing ($x_{0,1}$), idle-on ($x_{0,2}$) and busy ($y_{1}$)
servers. Also, Figure 1 (right) plots the average number of jobs per warm
server, which in state $x$ is given by $Q(x):=\frac{1}{y_{0}}\sum_{i\geq
1}ix_{i,2}$. The fluid (stochastic) trajectories are always represented by
dashed (continuous) lines and each curve is the average of ten simulations.
Each simulation is based on $10^{6}$ events. We have also set $\lambda=0.7$,
$\alpha=0.05$, $\beta=0.1$ and $\gamma=0.025$. If a time unit is 10
milliseconds, these are realistic choices [22, 10]. As scaling rule, we have
chosen Blind-$\theta$ where
$\theta=\frac{0.5}{\alpha}\frac{1-\lambda}{\frac{1}{\beta}+\frac{1}{\gamma}}$;
this choice will ensure that a strictly positive proportion of cold servers
exists in the long run (see Theorem 2). As dispatching algorithm, we have used
Power-of-$d$ with $d=1$. At time zero, we have assumed that the system is
minimally dimensioned for the average demand, i.e., $(1-\lambda)N$ servers are
cold and the remaining ones are idle-on.

Figure 1: Numerical convergence of the stochastic model $X^{N}(t)$ (continuous
lines), $N=10^{3}$, to the fluid model $x(t)$ (dashed lines) when combining
Power-of-$d$ with $d=1$ and Blind-$\theta$.

In both pictures, we observe that the fluid model captures the dynamics of
$X^{N}(t)$ accurately.

Let us provide some intuition about the dynamics in Figure 1. Initially, the
system is close to instability as capacity exactly matches demand. Here,
$Q(x(t))$ increases rapidly and as soon as a warm server is created, it is
filled with a job and as a result the proportion of idle-on servers decreases.
There is also another factor that makes the number of idle-on servers
decrease: Power-of-$d$, with $d=1$, is not so powerful in discovering idle-on
servers because a fraction of them go cold even in heavy load. This explains
why the number of cold servers (the blues lines) is increasing at the
beginning. Then, more warm servers are created to mitigate the effect of the
“close to instability” window on the accumulated overall number of jobs. Here,
the mass of busy servers ($y_{1}$) becomes greater than the average demand
$\lambda=0.7$ and $Q(x(t))$ decreases. Finally, dynamics stabilize and in
equilibrium there is a strictly positive fraction of servers that remain cold,
initializing and idle-on. This indicates that there is a flux of idle-on
servers that expires continuously even in equilibrium.

###  3.2 Characterization of Fixed Points

Theorem 1 justifies the use of the fluid model $x(t)$ to approximate the
behavior of the stochastic model $X^{N}(t)$. We now investigate properties of
the fluid model and in particular its fixed points when buffer sizes are
infinite. The fluid model has the form $\dot{x}=F(x)$; see Definition 1. We
say that $x^{*}\in\mathcal{S}_{1}$ is a _fixed point_ if $F(x^{*})=0$. Now,
for $x\in\mathcal{S}_{1}$, let us define the following conditions:

$\displaystyle x_{0,0}+x_{0,1}+x_{0,2}+\lambda=1$ (14a) $\displaystyle\beta
x_{0,1}=\gamma x_{0,2}$ (14b) $\displaystyle\gamma x_{0,2}\leq\alpha
g(x),\quad\mbox{ if }x_{0,0}=0$ (14c) $\displaystyle\gamma x_{0,2}=\alpha
g(x),\quad\mbox{ if }x_{0,0}>0$ (14d) $\displaystyle\mbox{ if Power-
of-}d\mbox{ is used: }$ $\displaystyle
x_{i,2}=(\lambda+x_{0,2})\left(\left(\tfrac{\lambda}{\lambda+x_{0,2}}\right)^{\frac{d^{i}-1}{d-1}}-\left(\tfrac{\lambda}{\lambda+x_{0,2}}\right)^{\frac{d^{i+1}-1}{d-1}}\right),$
$\displaystyle i\geq 1$ (14e) $\displaystyle\mbox{ if JBT-}d\mbox{ is used: }$
$\displaystyle\mbox{if }x_{0,2}=0:\quad x_{i,2}=0,\quad 0\leq i\leq d$ (14f)
$\displaystyle\,\qquad\qquad\quad\,\,\,x_{d+i,2}=x_{d+1,2}\left(1-\frac{{x}_{d+1,2}}{\lambda}\right)^{i-1},\quad
i\geq 2$ (14g) $\displaystyle\,\qquad\qquad\quad\,\,\,x_{d+1,2}\in(0,\lambda]$
(14h) $\displaystyle\,\qquad\qquad\quad\,\,\,g(x)=0$ (14i)
$\displaystyle\mbox{if }x_{0,2}>0:\quad
x_{i,2}=\left(\frac{\lambda}{z_{d}+x_{0,2}}\right)^{i}x_{0,2}\,\mathbb{I}_{\\{1\leq
i\leq d+1\\}},$ $\displaystyle i\geq 1$ (14j)

with $z_{d}\in[0,1]$ being the unique solution of

$\displaystyle
z_{d}+x_{0,2}=\frac{1-\left(\frac{\lambda}{z_{d}+x_{0,2}}\right)^{d+1}}{1-\frac{\lambda}{z_{d}+x_{0,2}}}x_{0,2}$
(15)

if $d\geq 1$ and $z_{d}=0$ if $d=0$. Here, $z_{d}$ is interpreted as the
proportion of busy servers with no more than $d$ jobs. Let us also introduce
the following assumption.

######  Assumption 3.

For any $x_{0,2}\in[0,1-\lambda]$, (14e)-(14j) uniquely determine $x_{i,2}$
for all $i\geq 1$.

Within Power-of-$d$, this assumption is clearly satisfied by (14e). Within
JBT-$d$, it is satisfied only if $x_{0,2}>0$, as if $x_{0,2}=0$, then
$x_{d+1,2}$ is only required to belong to $(0,\lambda]$. Under Assumption 3,
let $x^{\circ}=(x_{i,j}^{\circ})$ be the unique point in $\mathcal{S}_{1}$
such that $x_{0,0}^{\circ}=0$,
$x_{0,1}^{\circ}=\frac{\gamma}{\beta+\gamma}(1-\lambda)$,
$x_{0,2}^{\circ}=\frac{\beta}{\beta+\gamma}(1-\lambda)$.

The following result characterizes fixed points and provides a framework to
evaluate the equilibrium fluid performance induced by generic scaling rules;
see Section 7.2 for a proof.

######  Theorem 2.

If $x^{*}\in\mathcal{S}_{1}$ satisfies the conditions in (14)-(15), then it is
a fixed point of the fluid model with $B=+\infty$. In addition, under
Assumption 3

  1. 1.

$x_{0,0}^{*}>0$ if and only if

$\displaystyle\alpha
g(x^{\circ})<\frac{1-\lambda}{\frac{1}{\beta}+\frac{1}{\gamma}}.$ (16)

  2. 2.

If (16) does not hold, then $x^{\circ}$ is the unique fixed point.

###### Proof.

Given in Appendix 7. ∎

At the fluid scale and in a fixed point, Theorem 2 also provides the boundary
scaling probability that distinguishes between a saturated and a non-saturated
system. Specifically, if $g(x^{\circ})$ does not satisfy (16), then no cold
server exists in a fixed point but we observe that (14a) and (14b) imply that
a strictly positive fraction of servers remain initializing,
$\frac{\gamma}{\beta+\gamma}(1-\lambda)$. The interpretation is that there is
a mass of idle-on servers that go cold but instantly become initializing. This
corresponds to a waste of resources because initializing servers cannot
process jobs. In other words, a better performance may be obtained by keeping
the initializing servers warm at all times (no auto-scaling); recall also that
billing policies charge warm and initializing servers. If $g(x^{\circ})$
satisfies (16), then in a fixed point there still exists a fraction of idle-on
servers that go cold and instantly become initializing, provided that
$g(x^{\circ})>0$. However, the pool of cold servers remains non-empty and this
may give a benefit in terms of performance versus operational cost tradeoffs
because cold servers are not charged.

Within Blind-$\theta$, $g(x)=\theta$ and the conditions (14)-(15) easily
identify a unique fixed point, say $x^{*}$, with
$(x_{0,0}^{*},x_{0,1}^{*},x_{0,2}^{*})$ _not_ depending on the choice of the
load balancing algorithm. The following remark says that uniqueness is not
always guaranteed.

######  Remark 1 (Multiple Fixed Points).

Suppose that $g(x)=0$ whenever $y_{1}=\lambda=1-x_{0,0}$ and that JBT-$d$ is
used. Then, Theorem 2 implies that uncountably many fixed points exist. In
fact, while $x_{i,2}=0$ for all $i=0,\ldots,d$ and $x_{i,2}$ is uniquely
determined for all $i\geq d+2$ _once fixed_ $x_{d+1,2}$, the conditions
(14)-(15) do not tie $x_{d+1,2}\in(0,\lambda]$ to a specific value.

####  3.2.1 Blind-$\theta$ and Random Dispatching

For illustration purposes, let us consider Blind-$\theta$ with random
dispatching (Power-of-$1$). This combination does not involve any
communication overhead among the auto-scaler, dispatchers and servers, and for
this reason it is well suited for large systems with vast numbers of
dispatchers. Here, Theorem 2 identifies a unique fixed point, $x^{*}$. After
some algebra, we obtain
$x_{0,2}^{*}=\min\left\\{\frac{\alpha\theta}{\gamma},\,x_{0,2}^{\circ}\right\\}$
and for the mean queue length per warm server,
$Q(x)=\frac{1}{y_{0}}\sum_{i}ix_{i,2}$, we obtain (using also (14e))

$\displaystyle
Q(x^{*})=\frac{\lambda}{\min\left\\{\frac{\alpha\theta}{\gamma},\,x_{0,2}^{\circ}\right\\}}.$
(17)

As long as a strictly positive fraction of cold servers exists, or
equivalently $\frac{\alpha\theta}{\gamma}<x_{0,2}^{\circ}$, we remark that
$Q(x^{*})$ grows _linearly_ in $\lambda$.

###  3.3 Convergence to Multiple Fixed Points

As in the theory of dynamic systems, one is typically interested in whether or
not any fluid solution $x(t)$ converges to $x^{*}$, as $t\to\infty$,
regardless of the initial condition $x(0)$ and provided that $x^{*}$ is the
unique fixed point. In this case, $x^{*}$ is called “global attractor”.
Unfortunately, we have already shown in Remark 1 that uniqueness of a fixed
point is not always guaranteed. To guarantee uniqueness, the following
proposition shows that _it is not enough, as one may expect, to have a
strictly positive scaling probability whenever the current capacity of warm
servers is less than the average demand_ , i.e., $g(x)>0$ whenever
$y_{0}<\lambda$.

Let

$\displaystyle\mathcal{S}_{{\rm
subopt}}:=\\{x\in\mathcal{S}:x_{0,0}=1-\lambda,\,x_{0,1}=x_{0,2}=0,\,x_{1,2}<\lambda\mbox{
and \eqref{FP_PoD_conditions7aa} holds with }d=0\\}$ (18)

and let $\overline{Q}(x):=\sum_{i\geq 1}ix_{i,2}$ denote the average number of
jobs per server in state $x\in\mathcal{S}$; here, cold and initializing
servers are included in the counting.

######  Proposition 1.

Let $g(x)$ be any scaling rule such that

$\displaystyle g(x)=\frac{1}{\lambda}(x_{0,0}-1+\lambda)^{+},\qquad\forall
x\in\mathcal{S}:y_{0}<\lambda.$ (19)

Let $x(t)$ denote a fluid model induced by such $g(x)$ and JIQ such that

$\displaystyle x_{0,0}(0)>1-\lambda,\,\,x_{0,2}(0)=0,\,\,x_{1,2}(0)+\beta
x_{0,1}(0)<\lambda<\overline{Q}(x(0))<\infty.$ (20)

Suppose that $\beta<1$, $\alpha\neq\beta$ and $B=+\infty$. Then,

$\displaystyle g(x(t))$ $\displaystyle>0,\quad\forall t\geq 0$ (21)
$\displaystyle\lim_{t\to\infty}\overline{Q}(x(t))$
$\displaystyle=\overline{Q}(x(0))+\frac{\alpha+\beta}{\alpha\beta}(x_{0,0}(0)-1+\lambda)+\frac{x_{0,1}(0)}{\beta}>\lambda.$
(22)

In addition, $y_{0}(t)\uparrow\lambda$, and if $x_{1,2}(t)\to
x_{1,2}(\infty)$, then $x(t)\to x(\infty)$ with $x(\infty)\in\mathcal{S}_{{\rm
subopt}}$.

###### Proof.

Given in the Appendix. ∎

Thus, while the proportion of warm servers converges to $\lambda$, such
convergence may occur _from below_ even if there always exists a strictly
positive probability of creating new warm servers. In this case, the average
demand is greater than the current service capacity at any point in time and
this makes the mean queue length converge to a limit that depends on the
initial conditions.

Let us comment a little bit further and prepare the setting for our next
contribution. To create the underload situation above where
$y_{0}(t)\uparrow\lambda$, it is not necessary to assume that all warm servers
are initially busy ($x_{0,2}(0)=0$), though we have included this condition in
(20) to simplify our proof. In contrast, to avoid this situation, it may be
sufficient that $g(x)$ is bounded away from zero whenever $y_{0}<\lambda$. By
continuity, this implies that $g(x)>0$ as well whenever $y_{0}=\lambda$, but
in this case the resulting scaling rule will not possess the optimality
property stated in Theorem 3 below (as this will imply that $g(x^{\star})>0$).
On the other hand, one may consider a scaling rule that is discontinuous on
the set $\\{x:y_{1}=\lambda\\}$, a setting that does not satisfy Assumption 2.
Here, beyond revisiting Theorem 1 for justification of the fluid model, the
problem is that scale-up decisions would significantly depend on small
perturbations of the equilibrium system state, severely impacting robustness
from a practical standpoint.

###  3.4 Optimal (Stable) Design

Within Blind-$\theta$, Theorem 2 guarantees the existence of a unique fixed
point, say $x^{*}$, and all of our numerical simulations, which we omit,
indicate that it is a global attractor. Here, necessarily $x_{0,2}^{*}>0$, by
(14d), which means that a number of warm servers remain idle-on in
equilibrium. Clearly, this is not optimal for energy consumption because idle-
on servers consume energy. While it is mathematically interesting to find
conditions on the scaling rule that imply the existence of a global attractor,
our goal here is more ambitious: we are interested in how to drive the limit
behavior of any fluid solution $x(t)$ to the _optimal_ point
$x^{\star}\in\mathcal{S}$, uniquely defined by $x_{0,0}^{\star}=1-\lambda$ and
$x_{1,2}^{\star}=\lambda$, regardless of its initial condition
$x(0)\in\mathcal{S}_{1}$.

######  Remark 2 (Fluid Optimality).

Following [26], in $x^{\star}$ dynamics have achieved _“delay and relative
energy optimality”_ in the sense that both the waiting time of jobs and the
relative energy portion consumed by idle-on and initializing servers vanish in
the limit. Here, a possible intuition is that each job is always assigned to a
busy server with exactly one job but at the precise moment where it completes
the processing of its previous job. Therefore, service capacity perfectly
matches demand.

A direct consequence of Theorem 2 and (14d) is that it is necessary to impose
$g(x^{\star})=0$ to achieve fluid optimality. Within Power-of-$d$, this is
impossible as this condition would imply that $x_{0,2}=0$, and then (14e)
would imply $x_{i,2}=0$ for all $i$, contradicting that $\|x\|=1$. In fact,
Theorem 2 implies that the unique candidate is JIQ, though it leaves open the
possibility that $x(t)$ may converge to a fixed point in the sub-optimal set
$\mathcal{S}_{{\rm subopt}}$, see (18). Thus, it remains to understand what
additional structure the scaling rule $g(x)$ should satisfy to make
$x^{\star}$ a global attractor. Here, we observe that Remark 1 suggests that
even the knowledge of the amount of busy servers is not enough. More
precisely, it implies that one needs $g(x)>0$ for all $x\in\mathcal{S}_{{\rm
subopt}}$ as otherwise multiple fixed points exist. Therefore, given the
structure of $\mathcal{S}_{{\rm subopt}}$ and $x^{\star}$, we have the
following remark.

######  Remark 3.

A fluid optimal scaling rule needs the access to the amount of busy servers
with exactly one job, i.e., $x_{1,2}$.

The following result provides a general condition that yields fluid
optimality.

######  Theorem 3 (Optimal Design).

Let $\beta<1$ and let $x(t)$, with $x(0)\in\mathcal{S}_{1}$, denote a fluid
solution induced by JIQ and any scaling rule $g(x)$ that satisfies, beyond
Assumption 2,

$\displaystyle g(x)=0\leavevmode\nobreak\ {\rm if\leavevmode\nobreak\
and\leavevmode\nobreak\ only\leavevmode\nobreak\ if}\leavevmode\nobreak\
x_{1,2}+\beta x_{0,1}\geq\lambda.$ (23)

Then, $\lim_{t\to\infty}\|x(t)-x^{\star}\|_{w}=0$.

The key observation is that $x_{1,2}+\beta x_{0,1}$ is interpreted as the
overall rate at which servers become idle-on. Thus, our requirement for the
design of an optimal scaling rule is to ensure to scale up resources whenever
the excess of the mean demand over the rate at which servers become idle-on is
positive, as in this case JIQ is powerful enough to fill them up immediately
saturating the surplus service capacity. Otherwise, the scale-up process is
turned off ($g=0$), and in this case the natural dynamics induced by both JIQ
and the scale-down rule are enough to drive the system behavior to the
desirable configuration $x^{\star}$.

######  Remark 4.

_Rate-Idle_ , see Definition 2, satisfies (23).

As discussed in Section 2.1, the assumption $\beta<1$, i.e., the mean server
initialization rate is smaller than the mean job service rate, is largely
accepted in practice [22, 41]. From a mathematical standpoint, it is not
necessary for fluid optimality but simplifies our proof.

######  Remark 5 (Communication Overhead).

A scaling rule satisfying (23) requires the central controller to have access
to the amount of initializing and busy servers containing exactly one job,
i.e., $x_{0,1}$ and $x_{1,2}$. Since an initializing server informs the
platform as soon as it becomes warm, $x_{0,1}$ is easily obtained in practice.
For $x_{1,2}$, the auto-scaler can run a local memory with $N$ slots, where
the $n$-th slot indicates the state of server $n$, say ‘Cold’, ‘Init’, ‘Idle-
on’, ‘Busy1’ and ‘Busy≥2’, with obvious interpretations. Then, one way to
update the memory is by letting each server send a message to the auto-scaler
whenever the transitions ‘Busy≥2’ $\to$ ‘Busy1, ‘Busy1’ $\to$ ‘Idle-on’ and
‘Idle-on’ $\to$ ‘Busy1’ occur. As in standard implementations of JIQ, this
involves only a constant number of messages per job to be exchanged between
the auto-scaler and the servers.

Finally, the condition (23) is general enough to provide the platform user
with some flexibility when designing an optimal scaling rule. This opens the
way to a new level of optimizations as discussed in the next section.

##  4 Energy Optimization with Performance Guarantees

Given that the fluid model provides a tractable first-order approximation of
the stochastic dynamics within ALBA and that we have considered mild
assumptions on the structure of scaling rules (Assumption 2), it is natural to
use the fluid model in an optimization framework to trade off between
performance and energy costs. As proof of concept, the objective of this
section is to show a high-level application of our results that may serve as a
base for future more sophisticated techniques, which may integrate machine
learning mechanisms or consider switching costs.

We limit the focus on JIQ and on sets of scaling rules that satisfy the
assumptions in Theorem 3, say $\mathcal{G}$, as they automatically imply fluid
optimality (Remark 2) in the _stationary_ regime. Let us define the cost
function $\mathcal{J}_{g}$ as the long-run time average of a linear
combination between the power consumption
$P(x)=c_{0,1}x_{0,1}+c_{0,2}x_{0,2}+c_{1,2}y_{1}$, $c_{i,j}>0$, and the
average queue length per busy server $Q(x)=\frac{1}{y_{1}}\sum_{i}ix_{i,2}$
induced by the scaling rule $g$ (as in, e.g., [4]), i.e.,

$\displaystyle\mathcal{J}_{g}:=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\left(\kappa_{1}P(x(t))+\kappa_{2}Q(x(t))\right){\rm
d}t$ (24)

where $\kappa_{i}\geq 0$, $i=1,2$; one can think $\kappa_{1}$ in terms of
$/watt and $\kappa_{2}$ in terms of $/job. Then, Theorem 3 implies that

$\displaystyle\inf_{g}\mathcal{J}_{g}=\mathcal{J}_{g^{*}}=\kappa_{1}P(x^{\star})+\kappa_{2}Q(x^{\star})=\kappa_{1}c_{1,2}\lambda+\kappa_{2},\quad\forall
g^{*}\in\mathcal{G}.$ (25)

While all policies in $\mathcal{G}$ yield the same (optimal) cost, their
behavior is clearly different trajectory-wise. Depending on the application, a
platform user has several options to single out a policy in $\mathcal{G}$ that
satisfies a further level of optimization. For instance, a substantial portion
of the applications hosted in cloud networks have ultra-low delay
requirements, as this may have important consequences on e-commerce sales. On
the other hand, also energy bills are equally important from both financial
and environmental standpoints. Here, a system manager may want to look for a
scaling rule in $\mathcal{G}$ such that

$Q(x(t))\leq q,\quad\forall t\geq 0$ (26)

where $q$ is related to the desired user-perceived performance guarantee; by
Little’s law, (26) is equivalent to a constraint on the mean response time.
Towards this purpose, the parameterized subset of scaling rules taking the
form

$\displaystyle g(x)=\frac{1-\exp(-\frac{\eta}{\lambda}(\lambda-x_{1,2}-\beta
x_{0,1})^{+})}{1-\exp(-\eta)},\quad\eta>0,$ (27)

may be considered. Note that these satisfy both Assumption 2 and (23), that
Rate-Idle is recovered when $\eta\downarrow 0$, and that
$g(x)=\mathbb{I}_{\\{\lambda\geq x_{1,2}+\beta x_{0,1}\\}}$ when
$\eta\to\infty$. The control parameter $\eta>0$ indicates how aggressive the
scaling rule is. This provides intuition for considering scaling rule of the
form in (27). Then, the goal is to find $\eta$ such that (26) holds true,
provided that one exists. This problem can be easily addressed numerically
within the proposed deterministic model. Assume that the system is currently
in a light-load condition, say $\lambda=0.25$, and that, as a result, it is
dimensioned accordingly to save energy, say $x_{0,0}=1-\lambda-0.05$, with
$x_{0,2}=0.05$ and $x_{1,2}=\lambda$; the extra 0.05 is meant to keep a
reserve of idle-on servers ready to go. Then, at time zero, an unexpected
workload peak occurs, and $\lambda=0.5$. Here, the platform needs to
automatically adjust the service capacity while ensuring (26). Let us assume
$q=2$, $\alpha=0.35$, $\beta=0.1$ and $\gamma=0.025$. The dashed lines in
Figure 2 represent the dynamics of the fluid queue lengths $Q(x(t))$ and
scaling probabilities $g(x(t))$, for $\eta=1,10^{3}$. The corresponding
continuous lines represent the average of ten simulations of the stochastic
model $X^{N}(t)$ with $N=1000$.

Figure 2: Transient behavior of the queue lengths and scaling probabilities by
varying $\eta$, see (27), for the both the fluid ($x(t)$) and stochastic
($X^{N}(t)$) models with $N=1000$.

First, let us remark that the fluid model approximation accurately captures
the dynamics of $X^{N}(t)$, though it slightly underestimates queue lengths
and scaling probabilities. Now, let us consider $\eta=1$. Initially, queue
lengths increase as expected due to the surge of demand and the scaling
probability is large enough to drive the proportion of cold servers to zero.
This explains the non-differentiability point of the trajectory of the scaling
rule because the amount of initializing servers stops to grow. Then, the
system has enough capacity to drain the load and at some point the rate at
which servers become idle-on overflows the mean demand, i.e., $x_{1,2}+\beta
x_{0,1}>\lambda$, so that eventually $g(x(t))=0$. Finally, queue lengths
assess to their asymptotic value $Q(x^{\star})=1$ We conclude that $\eta=1$ is
enough to make (26) holds true. We also notice that the choice $\eta=10^{3}$,
which essentially means to scale up resources at the maximum available rate
$\alpha$ whenever $x_{1,2}+\beta x_{0,1}<\lambda$, has little impact on
performance. Nonetheless, it should be clear that the larger the value of
$\eta$, the larger the resulting time-average power consumption. A deeper
analysis is out of the scope of the present work.

##  5 Literature Review

The existing literature related to load balancing and auto-scaling is huge and
our goal here is to provide the necessary background highlighting the
difference with respect to our work.

The most popular examples of load balancing techniques that work well when
servers are homogeneous, i.e., all servers have the same processing speed, are
Random, Round-Robin (RR) [20], Power-of-$d$ [24], Join-the-Idle-Queue (JIQ)
[21], Least-Left-Workload (LLW) and Size Interval Task Allocation (SITA) [19,
18]. Random sends each job to random server, RR sends jobs to servers in a
cyclic manner, Power-of-$d$ sends an incoming job to the least loaded server
among $d$ selected uniformly at random. JIQ sends an incoming job to a random
idle server if an idle server exists and to a random one otherwise, LLW sends
an incoming jobs to the queue having the shortest workload, and SITA sends a
job to a given server if its size belong to a given interval. In all cases,
ties are broken randomly. In general, it is not possible to identify which of
these algorithms is the best because the general answer depends on the
underlying architecture, load conditions, service time distribution and on the
amount of information available to the dispatcher [37]. Recently, a number of
works attempted to understand under which conditions the mean waiting time can
be driven down to zero in the mean field limit, i.e., in the limiting regime
where the arrival rate grows linearly with the number of servers while keeping
the average load below one. It has been shown that this is possible within
different load balancing algorithms and architectures. Examples include JIQ
[32], Power-of-$d$ with $d\to\infty$ as the network size grows to infinity
[25], Power-of-$d$ with memory [6], SITA combined with RR [5] and the pull-
based policies developed in [13, 36]. To some extent, the fundamental limits
of load balancing are described in [13], where the authors investigate trade-
offs between performance, communication overhead and memory within a certain
class of symmetric architectures. There are several variants and
generalization of the above algorithms and these are meant to handle i)
heterogeneous servers [33, 9], ii) multiple dispatchers [35, 33], iii)
increased information to the dispatcher [17, 6], or iv) the case where model
parameters such as arrival and service rates are unknown [7].

The dynamics of the above dispatching algorithms have been analyzed under the
assumption that jobs are dispatched among a _constant_ number of parallel
servers. To the best of our knowledge, the only exceptions are [26, 16], which
provide the first performance model where load balancing and auto-scaling
operate jointly within the same timescale. Here, JIQ is combined with a
specific scaling strategy that turns a server off if that server remains idle
for a certain amount of time, and turns a server on at the moment of a job
arrival if all active servers are busy upon arrival of that job. In a regime
where the traffic demand and nominal service capacity grow large in
proportion, this mechanism yields zero user-perceived delay as in the
classical version of JIQ with no auto-scaling [32] but also deactivates any
surplus idle servers, ensuring that the relative energy wastage vanishes in
the limit as well. In this sense, this strategy induces asymptotic _delay and
relative energy optimality_. These properties have been strengthened in [27],
where the authors relax some finite buffer assumptions.

Let us comment on the differences between our work and [26, 16]. First, while
in these references the load balancing and auto-scaling processes are
_synchronous_ , in our framework they operate _asynchronously_ ; see below for
more details about the synchronous and asynchronous approach. Second, we do
not focus on a specific load balancing and auto-scaling strategy: our
framework allows for the design of general scaling rules (see Assumption 2)
and considers a wide class of load balancing schemes that is not limited to
JIQ. In fact, Power-of-$d$ may be the best option when the underlying
architecture is composed of multiple dispatchers because it still involves a
constant communication overhead per job, which is not the case within JIQ.
Third, we show that delay and relative energy optimality can be achieved under
broad assumptions (see Theorem 3). In turn, these yield the definition of
optimization frameworks allowing the system manager to further trade off
between energy and performance according to his needs.

Another reference that is close to our work is [22], where the authors propose
a performance model for auto-scaling that is tailored to the serverless
computing platform AWS Lambda. Here, the auto-scaling strategy is similar to
[26] in the sense that it is synchronized with job arrivals but different
because an incoming job that finds all servers busy is sent anyway, after a
coldstart latency, to an idle server that is turned on at the moment of its
arrival. If this is not possible because no further server is available, then
the job is rejected. Within this approach, there is no queueing because jobs
are always assigned to idle servers.

We also mention [11], where the authors propose a synchronous and
decentralized approach that is an hybrid between [22] and [26]. Upon arrival,
a job visit the active servers in a sequential order until a idle one is met.
This server will process the job if such server exists, otherwise the job
queues up in the last active server and this situation triggers a scale-up
decision if the number of jobs in such queue overflows a given threshold. This
architecture is meant to operate entirely within the network layer (Layer 4)
using IPv6 Segment Routing. As in [22], their performance analysis assumes a
timescale separation between the load balancing and auto-scaling processes.

There is a line of research works developing performance models for auto-
scaling in a different vein. One approach that has been widely considered in
the literature is based on the assumption that jobs join a global queue before
being dispatched to one out of a number of parallel servers. This type of
approach yields models that are variations of M/M/K or M/M/K/N queues [15, 43,
28, 4]. In contrast, the focus of this paper is on decentralized architectures
where each server has its own queue, which enhance scalability. A further
class of performance models for auto-scaling in cloud systems assume that job
dynamics reach stochastic equilibrium between consecutive changes of the
overall processing capacity [39, 31]; see also Section 8.4 in [29]. This
approach is arguable in the context of serverless computing platforms because
the auto-scaling process operates at a fast timescale [40, 41, 22]. Finally,
we mention that auto-scaling can be performed at the level of a single
centralized processor. This type of auto-scaling is referred to as “speed
scaling” [42, 4].

###  5.1 Synchronous vs Asynchronous in Serverless Computing

The load balancing and auto-scaling processes of existing implementations of
public serverless computing platforms can be “synchronous” or “asynchronous”;
this terminology is borrowed from the cloud computing community [22].

  * •

In platforms such as AWS Lambda, Azure Functions, IBM Cloud Functions and
Apache OpenWhisk, a new server (instance) is turned on at the arrival time of
a job if the job itself finds all servers busy [30, 41]. Then, once the server
is turned on, the job is executed on that server. The extra time that an
incoming job may incur if it finds all servers busy (usually referred to as
coldstart latency [22]) can pay off because of the order of $10^{2}$
milliseconds [41]. Sporadic excesses of such duration are acceptable for a
substantial number of applications hosted in these systems. Within these
technologies, the auto-scaling process is thus _synchronized_ with job arrival
and dispatching decision times, and all the above implementations require the
incoming jobs to wait in a _central_ queue before dispatching. This type of
auto-scaling is also referred to as _scale-per-request_ [22].

  * •

In platforms such as Google Cloud Run and Knative, the load balancing and
auto-scaling processes are decoupled in the sense that servers are turned on
_asynchronously_ at certain user-defined scaling times. Specifically, upon
arrival, a job is immediately dispatched to some running server according to
some load balancing algorithm such as Power-of-$d$ [2]. Independently of this,
the auto-scaling process decides whether the current processing capacity
should increase as a function of user-defined metrics, or scaling rules, that
may depend on instantaneous observations of the current system state. Because
of this decoupling, jobs may be assigned to busy servers. Asynchronous
implementations of serverless platforms such as Knative are _decentralized_ ,
i.e., each server has its own queue, and thus highly scalable. This type of
auto-scaling is also referred to as _concurrency-value-scaling_ [22].

######  Remark 6.

Existing implementations of decentralized serverless architectures are
asynchronous. In this setting, no dynamic performance model is available in
the literature, to the best of our knowledge. The main motivation underlying
this work has been to fill this gap.

##  6 Conclusions

In cloud systems, load balancing and auto-scaling are key mechanisms to
optimize both user-perceived delay performance and energy consumption. The
focus of the existing literature has been on architectures where these
mechanisms are synchronous or rely on a central queue. The novelty of our work
is to consider an asynchronous and decentralized architecture. In the limit
where the network demand grows in proportion to the nominal service capacity,
our work provides a tractable framework to evaluate the impact of existing or
new auto-scaling algorithms that are up to the platform user to choose or
design. In addition, it provides a general condition under which, with a
minimal communication overhead, it is possible to drive dynamics to the ideal
situation where the user-perceived delay and the relative energy wastage
induced by idle servers vanish. While this property also holds within the TABS
scheme recently proposed in [26], our contribution consists in showing that it
holds as well in asynchronous architectures provided that resources are scaled
up if and only if the mean demand is greater than the rate at which servers
become idle-on (Theorem 3). We have discussed in Remark 5 how such rate,
$x_{0,1}+\beta x_{1,2}$, can be estimated over time but here one may employ
other approaches, for instance based on machine learning techniques.

We discuss some perspectives and generalizations of our work:

  * •

The load-balancing algorithms considered within our framework are Power-of-$d$
and JBT-$d$. However, one may consider some of the algorithms mentioned in
Section 5 as well to see, for example, whether fluid optimality can be
obtained even without using JIQ (while keeping a comparable communication
overhead).

  * •

We have assumed that only one server at a time can be activated at each
scaling time. Our framework and results generalize trivially to the case where
a random number $C$ of cold servers is selected instead of one, provided that
the distribution of $C$ does not depend on the overall number of servers $N$.
Mutatis mutandis, it is enough to replace $\alpha$ by $\alpha\,\mathbb{E}[C]$.

  * •

Blind-$\theta$, see Definition 2, is the simplest scaling rule and it may be
the best option in systems with a large number of dispatchers because in this
setting the communication overhead induced by algorithms such as JIQ becomes a
limiting factor. While Theorem 2 implies the existence of a unique fixed
point, we leave as future work whether it is a global attractor. This is not
trivial as one may check that natural adaptations of classical Lyapunov
functions from queueing theory do not work.

  * •

Since we have proposed a class, instead of a specific instance, of fluid
optimal scaling rules, our work provides the system manager with some degrees
of freedom when designing a suitable scale-up strategy. This opens the way to
a new level of optimization problems such as the one presented in Section 4 as
an application of our results.

## References

  * [1] Knative docs v1.3.  https://knative.dev/docs/, 2022.  Online; accessed: 2022-04-04.
  * [2] Knative Load balancing.  https://knative.dev/docs/serving/load-balancing/, 2022.  Online; accessed: 2022-04-04.
  * [3] Knative scale bounds.  https://knative.dev/docs/serving/autoscaling/scale-bounds/, 2022.  Online; accessed: 2022-04-04.
  * [4] L. L. Andrew, M. Lin, and A. Wierman.  Optimality, fairness, and robustness in speed scaling designs.  In Proceedings of the ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems, SIGMETRICS ’10, page 37–48, New York, NY, USA, 2010. Association for Computing Machinery.
  * [5] J. Anselmi.  Combining size-based load balancing with round-robin for scalable low latency.  IEEE Transactions on Parallel and Distributed Systems, 31(4):886–896, 2020.
  * [6] J. Anselmi and F. Dufour.  Power-of-_d_ -choices with memory: Fluid limit and optimality.  Math. Oper. Res., 45(3):862–888, 2020.
  * [7] R. Atar, I. Keslassy, G. Mendelson, A. Orda, and S. Vargaftik.  Persistent-idle load-distribution.  Stochastic Systems, 10(2):152–169, 2020.
  * [8] M. Bramson.  State space collapse with application to heavy traffic limits for multiclass queueing networks.  Queueing Syst. Theory Appl., 30(1/2):89–148, June 1998.
  * [9] H. Chen and H.-Q. Ye.  Asymptotic optimality of balanced routing.  Oper. Res., 60(1):163–179, Jan. 2012.
  * [10] J. Dean and L. A. Barroso.  The tail at scale.  Commun. ACM, 56(2):74–80, Feb. 2013.
  * [11] Y. Desmouceaux, M. Enguehard, and T. H. Clausen.  Joint monitorless load-balancing and autoscaling for zero-wait-time in data centers.  IEEE Transactions on Network and Service Management, 18(1):672–686, 2021.
  * [12] Y. Desmouceaux, P. Pfister, J. Tollet, M. Townsley, and T. Clausen.  Srlb: The power of choices in load balancing with segment routing.  In 2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS), pages 2011–2016, 2017.
  * [13] D. Gamarnik, J. N. Tsitsiklis, and M. Zubeldia.  Delay, memory, and messaging tradeoffs in distributed service systems.  In Proceedings of the 2016 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Science, SIGMETRICS ’16, pages 1–12, New York, NY, USA, 2016. ACM.
  * [14] D. Gamarnik, J. N. Tsitsiklis, and M. Zubeldia.  Delay, memory, and messaging tradeoffs in distributed service systems.  Stochastic Systems, 8(1):45–74, 2018.
  * [15] A. Gandhi, S. Doroudi, M. Harchol-Balter, and A. Scheller-Wolf.  Exact analysis of the m/m/k/setup class of markov chains via recursive renewal reward.  In Proceedings of the ACM SIGMETRICS/International Conference on Measurement and Modeling of Computer Systems, SIGMETRICS ’13, page 153–166, New York, NY, USA, 2013. Association for Computing Machinery.
  * [16] D. Goldsztajn, A. Ferragut, F. Paganini, and M. Jonckheere.  Controlling the number of active instances in a cloud environment.  SIGMETRICS Perform. Eval. Rev., 45(3):15–20, Mar. 2018.
  * [17] V. Gupta and N. Walton.  Load balancing in the nondegenerate slowdown regime.  Operations Research, 67(1):281–294, 2019.
  * [18] M. Harchol-Balter, M. E. Crovella, and C. D. Murta.  On choosing a task assignment policy for a distributed server system.  Journal of Parallel and Distributed Computing, 59(2):204 – 228, 1999.
  * [19] M. Harchol-Balter, A. Scheller-Wolf, and A. R. Young.  Surprising results on task assignment in server farms with high-variability workloads.  SIGMETRICS ’09, pages 287–298, New York, NY, USA, 2009. ACM.
  * [20] Z. Liu and R. Righter.  Optimal load balancing on distributed homogeneous unreliable processors.  Operations Research, 46(4):563–573, 1998.
  * [21] Y. Lu, Q. Xie, G. Kliot, A. Geller, J. R. Larus, and A. Greenberg.  Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services.  Perform. Eval., 68(11):1056–1071, Nov. 2011.
  * [22] N. Mahmoudi and H. Khazaei.  Performance modeling of serverless computing platforms.  IEEE Transactions on Cloud Computing, pages 1–1, 2020.
  * [23] N. Mahmoudi, C. Lin, H. Khazaei, and M. Litoiu.  Optimizing serverless computing: Introducing an adaptive function placement algorithm.  In Proceedings of the 29th Annual International Conference on Computer Science and Software Engineering, CASCON ’19, page 203–213, USA, 2019\. IBM Corp.
  * [24] M. Mitzenmacher.  The power of two choices in randomized load balancing.  IEEE Trans. Parallel Distrib. Syst., 12(10):1094–1104, Oct. 2001\.
  * [25] D. Mukherjee, S. C. Borst, J. S. H. van Leeuwaarden, and P. A. Whiting.  Asymptotic Optimality of Power-of-$d$ Load Balancing in Large-Scale Systems.  ArXiv e-prints, Dec. 2016.
  * [26] D. Mukherjee, S. Dhara, S. C. Borst, and J. S. van Leeuwaarden.  Optimal service elasticity in large-scale distributed systems.  Proc. ACM Meas. Anal. Comput. Syst., 1(1), June 2017.
  * [27] D. Mukherjee and A. Stolyar.  Join idle queue with service elasticity: Large-scale asymptotics of a nonmonotone system.  Stochastic Systems, 9(4):338–358, 2019.
  * [28] H. Qian, D. Medhi, and K. Trivedi.  A hierarchical model to evaluate quality of experience of online services hosted by cloud computing.  In 12th IFIP/IEEE International Symposium on Integrated Network Management (IM 2011) and Workshops, pages 105–112, 2011.
  * [29] C. Qu, R. N. Calheiros, and R. Buyya.  Auto-scaling web applications in clouds: A taxonomy and survey.  ACM Comput. Surv., 51(4), July 2018.
  * [30] M. Shahrad, R. Fonseca, I. Goiri, G. Chaudhry, P. Batum, J. Cooke, E. Laureano, C. Tresness, M. Russinovich, and R. Bianchini.  Serverless in the wild: Characterizing and optimizing the serverless workload at a large cloud provider.  In 2020 USENIX Annual Technical Conference (USENIX ATC 20), pages 205–218. USENIX Association, July 2020.
  * [31] U. Sharma, P. Shenoy, and D. F. Towsley.  Provisioning multi-tier cloud applications using statistical bounds on sojourn time.  In Proceedings of the 9th International Conference on Autonomic Computing, ICAC ’12, page 43–52, New York, NY, USA, 2012. Association for Computing Machinery.
  * [32] A. L. Stolyar.  Pull-based load distribution in large-scale heterogeneous service systems.  Queueing Syst. Theory Appl., 80(4):341–361, Aug. 2015.
  * [33] A. L. Stolyar.  Pull-based load distribution among heterogeneous parallel servers: The case of multiple routers.  Queueing Syst. Theory Appl., 85(1-2):31–65, Feb. 2017.
  * [34] J. N. Tsitsiklis and K. Xu.  On the power of (even a little) resource pooling.  Stoch. Syst., 2(1):1–66, 2012.
  * [35] M. van der Boor, S. Borst, and J. van Leeuwaarden.  Load balancing in large-scale systems with multiple dispatchers.  In IEEE INFOCOM 2017 - IEEE Conference on Computer Communications, pages 1–9, May 2017.
  * [36] M. van der Boor, S. C. Borst, and J. van Leeuwaarden.  Hyper-scalable JSQ with sparse feedback.  Proc. ACM Meas. Anal. Comput. Syst., 3(1):4:1–4:37, 2019.
  * [37] M. van der Boor, S. C. Borst, J. S. van Leeuwaarden, and D. Mukherjee.  Scalable load balancing in networked systems: A survey of recent advances.  arXiv preprint arXiv:1806.05444, 2018.
  * [38] E. van Eyk, A. Iosup, C. L. Abad, J. Grohmann, and S. Eismann.  A spec rg cloud group’s vision on the performance challenges of faas cloud architectures.  In Companion of the 2018 ACM/SPEC International Conference on Performance Engineering, ICPE ’18, page 21–24, New York, NY, USA, 2018. Association for Computing Machinery.
  * [39] D. Villela, P. Pradhan, and D. Rubenstein.  Provisioning servers in the application tier for e-commerce systems.  In Twelfth IEEE International Workshop on Quality of Service, 2004\. IWQOS 2004., pages 57–66, 2004.
  * [40] L. Wang, M. Li, Y. Zhang, T. Ristenpart, and M. Swift.  Peeking behind the curtains of serverless platforms.  In Proceedings of the 2018 USENIX Conference on Usenix Annual Technical Conference, USENIX ATC ’18, page 133–145, USA, 2018. USENIX Association.
  * [41] L. Wang, M. Li, Y. Zhang, T. Ristenpart, and M. Swift.  Peeking behind the curtains of serverless platforms.  In Proceedings of the 2018 USENIX Conference on Usenix Annual Technical Conference, USENIX ATC ’18, page 133–145, USA, 2018. USENIX Association.
  * [42] M. Weiser, B. Welch, A. Demers, and S. Shenker.  Scheduling for reduced cpu energy.  In Proceedings of the 1st USENIX Conference on Operating Systems Design and Implementation, OSDI ’94, page 2–es, USA, 1994. USENIX Association.
  * [43] B. Yang, F. Tan, Y.-S. Dai, and S. Guo.  Performance evaluation of cloud service considering fault recovery.  In Proceedings of the 1st International Conference on Cloud Computing, CloudCom ’09, page 571–576, Berlin, Heidelberg, 2009. Springer-Verlag.

##  7 Proofs of Theorems 1, 2 and 3

We provide proofs for Theorems 1, 2 and 3.

###  7.1 Theorem 1: Connection between the Fluid and the Markov models

To prove Theorem 1, we follow two main steps. First, we use a common coupling
technique to define the processes $(X^{N}(t))_{t\in[0,T]}$, for all
$N\in\mathbb{Z}_{+}$, on a common probability space and show that limit
trajectories exist and are Lipschitz continuous with probability one. The
arguments used in this step are quite standard [34, 8, 13]. Then, we prove
that limit trajectories are fluid solutions, which is the main technical
difficulty, and here we develop arguments specific to our model.

####  7.1.1 Coupled Construction of Sample Paths

Let $\mathcal{N}_{c}(t)$ denote a Poisson process of rate $c$. We construct a
probability space where the stochastic processes
$\\{(X^{N}(t))_{t\in[0,T]}\\}_{N\geq 1}$ are coupled. All the processes of
interest can be constructed in terms of the following mutually independent
primitive processes:

  * •

$\mathcal{N}_{\phi}(t)$, a Poisson process of rate
$\phi:=\lambda+1+\alpha+\beta+\gamma$. This process is defined on
$(\Omega_{E},\mathcal{A}_{E},\mathbb{P}_{E})$ and each jump of
$\mathcal{N}_{\phi}(t)$ denotes the occurrence of an _event_.

  * •

$(W_{n})_{n}$, where the random variables $W_{n}$ are $\\{0,1,2,3,4\\}$-valued
i.i.d. and such that $\mathbb{P}(W_{n}=0)=\lambda/\phi$,
$\mathbb{P}(W_{n}=1)=1/\phi$, $\mathbb{P}(W_{n}=2)=\alpha/\phi$
$\mathbb{P}(W_{n}=3)=\beta/\phi$ and $\mathbb{P}(W_{n}=4)=\gamma/\phi$. This
process is defined on $(\Omega_{W},\mathcal{A}_{W},\mathbb{P}_{W})$ and will
identify the _type_ of the $n$-th event. Specifically, $W_{n}=0$ indicates a
job arrival, $W_{n}=1$ a potential job departure, $W_{n}=2$ a scaling time,
$W_{n}=3$ a potential server initialization, i.e., a server completed the
initialization phase, and $W_{n}=4$ a potential server expiration.

  * •

$(A_{n}^{p})_{n}$, $p=1,\ldots,d$, $(D_{n})_{n}$, $(I_{n})_{n}$, $(E_{n})_{n}$
and $(R_{n})_{n}$, where the random variables $A_{n}^{p}$, $D_{n}$, $I_{n}$,
$E_{n}$ and $R_{n}$, for all $n$, are all i.i.d. and uniform over the interval
$[0,1]$. The rvs $A_{n}^{p}$, $D_{n}$, $I_{n}$, $E_{n}$ will be respectively
used to select a server that i) will process an arriving job, ii) fires a
departure, iii) fires an initialization and iv) fires an expiration. The rv
$R_{n}$ is related to the scaling rule and will decide whether a new server
will be activated. These processes are defined on
$(\Omega_{S},\mathcal{A}_{S},\mathbb{P}_{S})$;

  * •

$(X^{N}(0))_{N}$, the process of the initial conditions, where each random
variable $X^{N}(0)$ takes values in $\mathcal{S}_{N}$. This process is defined
on $(\Omega_{0},\mathcal{A}_{0},\mathbb{P}_{0})$.

Using that $\mathcal{N}_{\phi}(Nt)$ and $\mathcal{N}_{\phi N}(t)$ are equal in
distribution and the well-known fact that thinnings of a Poisson process
produce independent Poisson processes, each process
$\\{(X^{N}(t))_{t\in[0,T]}\\}$, $N\geq 1$, can be constructed on the product
space, say $(\Omega,\mathcal{A},\mathbb{P})$.

Now, let $t_{n}$ be the time of the $n$-th jump of $\mathcal{N}_{\phi}(Nt)$.
Let also $X^{N}(t^{-}):=\lim_{s\uparrow t}X^{N}(s)$,
$Y_{i}^{N}(t):=\sum_{j=0}^{i}X_{j,2}^{N}(t)$ for all $i\geq 0$,
$Y_{-1}^{N}(t)=0$ and $\mathbf{1}_{A}^{x}=1$ if $x\in A$ and 0 otherwise. The
coordinates of $X^{N}(t)$ are then given by

$\displaystyle X_{0,0}^{N}(t)=$ $\displaystyle
X_{0,0}^{N}(0)+\frac{1}{N}\sum_{n=1}^{\mathcal{N}_{\phi}(Nt)}\left(\mathbb{I}_{\\{W_{n}=4\\}}\mathbf{1}_{(0,X_{0,2}^{N}(t_{n}^{-})]}^{E_{n}}-\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N}(t_{n}^{-}))]}^{R_{n}}\right)$
(28a) $\displaystyle X_{0,1}^{N}(t)=$ $\displaystyle
X_{0,1}^{N}(0)+\frac{1}{N}\sum_{n=1}^{\mathcal{N}_{\phi}(Nt)}\left(\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N}(t_{n}^{-}))]}^{R_{n}}-\mathbb{I}_{\\{W_{n}=3\\}}\mathbf{1}_{(0,X_{0,1}^{N}(t_{n}^{-})]}^{I_{n}}\right)$
(28b) $\displaystyle X_{0,2}^{N}(t)=$ $\displaystyle
X_{0,2}^{N}(0)+\frac{1}{N}\sum_{n=1}^{\mathcal{N}_{\phi}(Nt)}\left(\mathbb{I}_{\\{{W_{n}=1}\\}}\mathbf{1}_{[Y_{0}^{N}(t_{n}^{-}),Y_{1}^{N}(t_{n}^{-})]}^{D_{n}}-\mathbb{I}_{\\{{W_{n}=0}\\}}H_{0}(t_{n}^{-})\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\,\,\left.+\mathbb{I}_{\\{W_{n}=3\\}}\mathbf{1}_{(0,X_{0,1}^{N}(t_{n}^{-})]}^{I_{n}}-\mathbb{I}_{\\{W_{n}=4\\}}\mathbf{1}_{(0,X_{0,2}^{N}(t_{n}^{-})]}^{E_{n}}\right)$
(28c) $\displaystyle X_{i,2}^{N}(t)=$ $\displaystyle
X_{i,2}^{N}(0)+\frac{1}{N}\sum_{n=1}^{\mathcal{N}_{\phi}(Nt)}\bigg{(}\mathbb{I}_{\\{W_{n}=0\\}}\left(H_{i-1}(t_{n}^{-})-H_{i}(t_{n}^{-})\mathbb{I}_{\\{i<B\\}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\,\,\,\left.+\mathbb{I}_{\\{W_{n}=1\\}}\left(\mathbf{1}_{(Y_{i}^{N}(t_{n}^{-}),Y_{i+1}^{N}(t_{n}^{-})]}^{D_{n}}-\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{D_{n}}\right)\right)$
(28d)

for all $i\geq 1$, where within Power-of-$d$ (servers are selected with
replacement)

$\displaystyle
H_{i}(t_{n}^{-}):=\prod_{p=1}^{d}\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),1]}^{A_{n}^{p}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}-\prod_{p=1}^{d}\mathbf{1}_{(Y_{i}^{N}(t_{n}^{-}),1]}^{A_{n}^{p}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\in\\{0,1\\}$
(29)

and within JBT-$d$

$\displaystyle
H_{i}(t_{n}^{-}):=\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\mathbb{I}_{\\{Y_{d}^{N}(t_{n}^{-})=0\\}}+\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}Y_{d}^{N}(t_{n}^{-})}\mathbb{I}_{\\{i\leq
d\\}}\mathbb{I}_{\\{Y_{d}^{N}(t_{n}^{-})>0\\}}\in\\{0,1\\}.$ (30)

These expressions follow by uniformization of $X^{N}(t)$. For instance,
$X_{0,0}^{N}(t)$ has an upward jump of size $1/N$ at time $t_{n}$ if the event
occurring at that time is of type 4 (potential server expiration) and an idle-
on server is actually selected at time $t_{n}^{-}$ by the uniformized process.
Analogously, $X_{0,0}^{N}(t)$ decreases by $1/N$ at time $t_{n}$ if the event
occurring at that time is of type 2, provided that at time $t_{n}^{-}$ the
cold servers pool is not empty and the scaling rule applies. Similar
interpretations hold along the other coordinates of $X^{N}(t)$.

####  7.1.2 Tightness of Sample Paths and Lipschitz Property

We now prove tightness of sample paths. The lemmas in this section are
equivalent to the lemmas in [14, Section 5.2].

Let us introduce the following formulas for quick reference.

######  Lemma 1.

Let $T>0$. There exists $\mathcal{C}\subseteq\Omega$ such that
$\mathbb{P}(\mathcal{C})=1$ and for all $\omega\in\mathcal{C}$:

$\displaystyle\lim_{N\to\infty}\sup_{t\in[0,T]}|\frac{1}{N}\mathcal{N}_{\phi}(Nt,\omega)-\phi
t|=0,$ (31)
$\displaystyle\lim_{N\to\infty}\sup_{t\in[0,T]}\Big{|}\frac{1}{N}\sum_{n=1}^{\mathcal{N}_{\phi}(Nt,\omega)}\mathbb{I}_{\\{W_{n}(\omega)=k\\}}-\mathbb{P}(W_{1}=k)\,\phi\,t\,\Big{|}=0,\qquad
k\in\\{0,\ldots,4\\}$ (32)
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{n=1}^{N}\prod_{p=1}^{d}\mathbf{1}_{(a_{p},b_{p}]}^{c_{p}A_{n}^{p}}=\prod_{p=1}^{d}\frac{b_{p}-a_{p}}{c_{p}},\quad\forall
a_{p},b_{p},c_{p}\in[0,1],c_{p}>0,p=1,\ldots,d.$ (33)

###### Proof.

This lemma directly follows by applying the functional strong law of large
numbers for the Poisson process (for (31)), the fact that thinnings of a
Poisson process produce independent Poisson processes (for (32)) and the
strong law of the large numbers (for (33)). ∎

We will work on a fixed $\omega$ that belongs to $\mathcal{C}$.

Let $x^{0}\in[0,1]$, sequences $A_{N}\downarrow 0$ and $B_{N}\downarrow 0$ be
given. Let also $D[0,T]$ denote the Skorokhod space endowed with the uniform
metric $d(x,y):=\sup_{t\in[0,T]}|x(t)-y(t)|$, for all $x,y\in D[0,T]$. For
$N\geq 1$, let also

$\displaystyle\mathcal{E}_{N}(B_{N},A_{N},x^{0})$
$\displaystyle:=\big{\\{}x\in D[0,T]:|x(0)-x^{0}|\leq
B_{N},\,|x(a)-x(b)|\leq\phi|a-b|+A_{N},\,\forall a,b\in[0,T]\big{\\}}$
$\displaystyle\mathcal{E}_{c}(x^{0})$ $\displaystyle:=\big{\\{}x\in
D[0,T]:x(0)=x^{0},\,|x(a)-x(b)|\leq\phi|a-b|,\,\forall a,b\in[0,T]\big{\\}}.$

The next lemma says that the sample paths along any coordinate is
approximately Lipschitz continuous. The proof is omitted because follows
exactly the same standard arguments used in Lemma 5.2 of [14], which basically
use the fact that the jumps of the Markov chain of interest are of the order
of $1/N$ and that the evolution of such Markov chain on a given coordinate
only depends on the evolution of such Markov chain on a finite number of other
coordinates.

######  Lemma 2.

Fix $T>0$, $\omega\in\mathcal{C}$, and some $x^{0}\in\mathcal{S}_{1}$. Suppose
that $\|X^{N}(\omega,0)-x^{0}\|_{w}\leq\tilde{B}_{N}$, for some sequence
$\tilde{B}_{N}\downarrow 0$. Then, there exists sequences
$\left\\{B_{N}^{(i,j)}\downarrow 0\right\\}_{i,j}$ and $A_{N}\downarrow 0$
such that

$\displaystyle
X_{i,j}^{N}(\omega,\cdot)\in\mathcal{E}_{N}(B_{N}^{(i,j)},A_{N},x^{0}),\quad\forall(i,j),\,\forall
N.$ (34)

The next proposition shows that any sequence of sample paths $X^{N}(\omega,t)$
contains a further subsequence that converges in $D^{\infty}[0,T]$, endowed
with the metric $d^{\mathbb{Z}_{+}}(x,y):=\sup_{t\in[0,T]}\|x(t)-y(t)\|_{w}$,
to a coordinate-wise Lipschitz continuous trajectory $x(t)$, as long as
$\omega\in\mathcal{C}$. The proof is routine and omitted because it is a
repetition of the argument used in the proof of Proposition 11 in [34]
(equivalently, see also Proposition 5.3 in [14]).

######  Proposition 2.

Fix $T>0$, $\omega\in\mathcal{C}$, and some $x^{0}\in\mathcal{S}_{1}$. Suppose
that $\|X^{N}(\omega,0)-x^{0}\|_{w}\leq\tilde{B}_{N}$, for some sequence
$\tilde{B}_{N}\downarrow 0$. Then, every subsequence of
$\\{X^{N}(\omega,\cdot)\\}_{N=1}^{\infty}$ contains a further subsequence
$\\{X^{N_{k}}(\omega,\cdot)\\}_{k=1}^{\infty}$ such that

$\displaystyle\lim_{k\to\infty}d^{\mathbb{Z}_{+}}(X^{N_{k}},x)=0$ (35)

where $x(0)=x^{0}$ and $x_{i,j}\in\mathcal{E}_{c}(x^{0})$, for all $i$ and
$j$.

Since Lipschitz continuity implies absolute continuity, we have obtained that
limit points of $X^{N}(t)$ exist and are absolutely continuous. Since all
sample paths of $X^{N}(t)$ take values in $\mathcal{S}$, these limit points
must belong as well to $\mathcal{S}$ because $\mathcal{S}$ is a closed set.
Therefore, to conclude the proof of Theorem 1 it remains to show that the
derivative of $x_{i,j}(t)$ is as in Definition 1 for all $i$ and $j$, provided
that $t$ is a regular time. This is done in the next subsection and will also
prove that a fluid solution started in $x^{(0)}\in\mathcal{S}_{1}$ exists.

####  7.1.3 Limit Trajectories are Fluid Solutions

Fix $\omega\in\mathcal{C}$ and let $\\{X^{N_{k}}(\omega,t)\\}_{k=1}^{\infty}$
be a subsequence that converges to $\overline{x}$ (by Proposition 2), i.e.

$\lim_{k\to\infty}\sup_{t\in[0,T]}\|X^{N_{k}}(\omega,t)-\overline{x}(t)\|_{w}=0.$
(36)

In the remainder, we fix such $\omega\in\mathcal{C}$ such that (36) holds and
for simplicity we drop the dependency on $\omega$. Since $\overline{x}$ must
be Lipschitz continuous (by Proposition 2), it is also absolutely continuous
and to conclude the proof of Theorem 1, it remains to show that
$\overline{x}(t)$ satisfies the conditions on the derivatives given in
Definition 1 whenever $\overline{x}_{i,j}(t)$ is differentiable, for all
$i,j$.

We say that $t$ is a point of differentiability (of $\overline{x}$) if
$x_{i,j}(t)$ is differentiable for all $i,j$.

We will (implicitly) use several times the following elementary lemma, which
holds true because $\overline{x}$ is a non-negative absolutely continuous
function.

######  Lemma 3.

If $\overline{x}_{i,j}(t)=0$ and $t$ is a point of differentiability of
$\overline{x}_{i,j}$, then $\dot{\overline{x}}_{i,j}(t)=0$.

Let $\epsilon>0$. By Lemma 2, there exists a sequence $A_{N_{k}}\downarrow 0$
such that $X_{i,j}^{N_{k}}(\omega,u)\in[\overline{x}_{i,j}(t)-\epsilon\phi-
A_{N_{k}},\overline{x}_{i,j}(t)+\epsilon\phi+A_{N_{k}}]$, for all
$u\in[t,t+\epsilon]$. Thus, for all $k$ sufficiently large,
$X_{i,j}^{N_{k}}(\omega,u)\in[\overline{x}_{i,j}(t)-2\epsilon\phi,\overline{x}_{i,j}(t)+2\epsilon\phi]$,
for all $u\in[t,t+\epsilon]$. Thus, we have

$\displaystyle|X_{i,j}^{N_{k}}(u)-\overline{x}_{i,j}(t)|\leq
2\phi\epsilon,\quad\forall u\in[t,t+\epsilon]$ (37)

for all $k$ sufficiently large. In addition, using (37) and that $g$ is
Lipschitz, we obtain

$\displaystyle|g(X^{N_{k}}(u))-g(\overline{x}(u))|$ $\displaystyle\leq
L\|X^{N_{k}}(u))-\overline{x}(u)\|_{w}$ (38a) $\displaystyle\leq 2\phi\epsilon
L\sqrt{\sum_{i,j}\frac{1}{2^{i+j}}}=2\phi\epsilon L\sqrt{2},\quad\forall
u\in[t,t+\epsilon]$ (38b)

where $L$ is the Lipschitz constant of the scaling rule $g$.

We will refer to the following lemma, which is a straightforward consequence
of (37) and of the strong law of the large numbers. In points where the fluid
drift function is continuous, it will provide an expression for terms related
to job departures, server initializations/departures and, in some cases,
dispatching decisions.

######  Lemma 4.

Fix $\omega\in\mathcal{C}$ and let (36) hold. Then,

$\displaystyle\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=1\\}}\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{D_{n}}$
$\displaystyle\,=\overline{x}_{i,2}(t)$ $\displaystyle\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=3\\}}\mathbf{1}_{(0,X_{0,1}^{N}(t_{n}^{-})]}^{I_{n}}$
$\displaystyle\,=\beta\overline{x}_{0,1}(t)$
$\displaystyle\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=4\\}}\mathbf{1}_{(0,X_{0,2}^{N}(t_{n}^{-})]}^{E_{n}}$
$\displaystyle\,=\gamma\overline{x}_{0,2}(t).$

In addition,

$\displaystyle\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}Y_{d}^{N}(t_{n}^{-})}\mathbb{I}_{\\{i\leq
d\\}}$ $\displaystyle\,=\frac{\lambda\mathbb{I}_{\\{i\leq
d\\}}\overline{x}_{i}(t)}{\sum_{j=0}^{d}\overline{x}_{j,2}(t)}$

provided that $\sum_{j=0}^{d}\overline{x}_{j,2}>0$, and

$\displaystyle\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}H_{i}(X^{N}(t_{n}^{-}))$
$\displaystyle\,=\lambda h_{i}(\overline{x})$

provided that Power-of-$d$ is used.

###### Proof.

Given in Section 8. ∎

The next proposition proves the desired condition on the amount of fluid of
cold and initializing servers.

######  Proposition 3.

Fix $\omega\in\mathcal{C}$, let (36) hold and assume that $t$ is a point of
differentiability. Then,

$\displaystyle\dot{\overline{x}}_{0,0}$
$\displaystyle=\gamma\overline{x}_{0,2}(t)-\alpha\mathbb{I}_{\\{\overline{x}_{0,0}(t)>0\\}}g(\overline{x}(t))-\gamma\overline{x}_{0,2}(t)\,\mathbb{I}_{\\{\overline{x}_{0,0}(t)=0,\,\gamma\overline{x}_{0,2}(t)\leq\alpha
g(\overline{x}(t))\\}}$ (39) $\displaystyle\dot{\overline{x}}_{0,1}$
$\displaystyle=\alpha
g(\overline{x}(t))\mathbb{I}_{\\{\overline{x}_{0,0}(t)>0\\}}-\beta\overline{x}_{0,1}(t)+\gamma\overline{x}_{0,2}(t)\,\mathbb{I}_{\\{\overline{x}_{0,0}(t)=0,\,\gamma\overline{x}_{0,2}(t)\leq\alpha
g(\overline{x}(t))\\}}.$ (40)

###### Proof.

Assume that $\overline{x}_{0,0}(t)>0$ and let
$\epsilon\in(0,\frac{\overline{x}_{0,0}(t)}{2\phi})$. Given that

$\displaystyle t_{n}\in(t,t+\epsilon]\mbox{ whenever
}n\in\\{\mathcal{N}_{\phi}({N_{k}}t)+1,\ldots,\mathcal{N}_{\phi}({N_{k}}(t+\epsilon))\\},$
(41)

(37) implies that for all $k$ sufficiently large,
$|X_{0,0}^{N_{k}}(t_{n}^{-})-\overline{x}_{0,0}(t)|\leq
2\phi\epsilon<\overline{x}_{0,0}(t)$ and thus $X_{0,0}^{N_{k}}(t_{n}^{-})>0$.
We have shown that

$\mathbb{I}_{\\{X_{0,0}^{N_{k}}(t_{n}^{-})>0\\}}=1,\quad\forall
n\in\\{\mathcal{N}_{\phi}({N_{k}}t)+1,\ldots,\mathcal{N}_{\phi}({N_{k}}(t+\epsilon))\\}$
(42)

for all $k$ sufficiently large. Using (28), Lemma 4 and (42), we have

$\displaystyle\dot{\overline{x}}_{0,0}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\left(X_{0,0}^{N_{k}}(t+\epsilon)-X_{0,0}^{N_{k}}(t)\right)$
$\displaystyle=\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\left(\mathbb{I}_{\\{W_{n}=4\\}}\mathbf{1}_{(0,X_{0,2}^{N_{k}}(t_{n}^{-})]}^{E_{n}}-\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N_{k}}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}\right)$
(43a) $\displaystyle=\gamma\overline{x}_{0,2}(t)-\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=2\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}.$
(43b)

Since $t$ is a point of differentiability, the double limit in the RHS of
(43b) exists. Then, (38) implies that given $\epsilon>0$ small enough,
$g(X^{N_{k}}(t_{n}^{-}))\in[g(\overline{x}(t))-2\phi\epsilon
L\sqrt{2},g(\overline{x}(t))+2\phi\epsilon L\sqrt{2}]$ for all $k$
sufficiently large. Combining these bounds with Lemma 1 and letting
$\epsilon\downarrow 0$ (as in the proof of Lemma 4), we obtain

$\displaystyle\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=2\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}=\alpha
g(\overline{x}(t)).$ (44)

Similarly, on coordinates (0,1), we obtain

$\displaystyle\dot{\overline{x}}_{0,1}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\left(X_{0,1}^{N_{k}}(t+\epsilon)-X_{0,1}^{N_{k}}(t)\right)$
$\displaystyle=\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\left(\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N_{k}}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}-\mathbb{I}_{\\{W_{n}=3\\}}\mathbf{1}_{(0,X_{0,1}^{N_{k}}(t_{n}^{-})]}^{I_{n}}\right)$
$\displaystyle=\alpha g(\overline{x})-\beta\overline{x}_{0,1}(t).$

Now, let us assume that $\overline{x}_{0,0}(t)=0$. First, we notice that

$\displaystyle\dot{\overline{x}}_{0,0}(t)$
$\displaystyle=\gamma\overline{x}_{0,2}(t)-\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}$
(45) $\displaystyle\geq\gamma\overline{x}_{0,2}(t)-\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=2\\}}\,\mathbf{1}_{(0,g(X^{N_{k}}(t_{n}^{-}))]}^{R_{n}}=\gamma\overline{x}_{0,2}(t)-\alpha
g(\overline{x})$

where the first equality follows by (43a) and Lemma 4, and the last equality
follows by (44). Thus, if $\overline{x}_{0,0}(t)=0$ and
$\gamma\overline{x}_{0,2}(t)>\alpha g(\overline{x})$, then by the previous
inequality $\dot{\overline{x}}_{0,0}(t)>0$, which is not possible because if
$t$ is a point of differentiability and $\overline{x}_{0,0}(t)=0$ then
necessarily $\dot{\overline{x}}_{0,0}(t)=0$ as $\overline{x}_{0,0}$ is a non-
negative absolutely continuous function. Thus, in a point of differentiability
$t$ where $\overline{x}_{0,0}(t)=0$, we must have
$\gamma\overline{x}_{0,2}(t)\leq\alpha g(\overline{x})$. and, necessarily,
$\dot{\overline{x}}_{0,0}(t)=0$. In this case, (45) gives

$\displaystyle\gamma\overline{x}_{0,2}(t)$
$\displaystyle=\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=2\\}}\mathbb{I}_{\\{X_{0,0}^{N}(t_{n}^{-})>0\\}}\,\mathbf{1}_{(0,g(X^{N}(t_{n}^{-}))]}^{R_{n}}.$
(46)

This term is interpreted as the amount of idle-on servers that become cold but
instantly turn initializing. Substituting (46) in the previous equalities
within the conditions $\gamma\overline{x}_{0,2}(t)\leq\alpha g(\overline{x})$
and $\overline{x}_{0,0}(t)=0$, we obtain (39) and (40). ∎

On the coordinates associated to warm servers, it remains to prove that

$\displaystyle\dot{\overline{x}}_{0,2}(t)=\overline{x}_{1,2}(t)-\lambda
h_{0}(\overline{x}(t))+\beta\overline{x}_{0,1}(t)-\gamma\overline{x}_{0,2}(t)$
(47)
$\displaystyle\dot{\overline{x}}_{i,2}(t)=\overline{x}_{i+1,2}(t)\mathbb{I}_{\\{i<B\\}}-\overline{x}_{i,2}(t)+\lambda(h_{i-1}(\overline{x}(t))-h_{i}(\overline{x}(t))\mathbb{I}_{\\{i<B\\}}),\quad
i\geq 1,$ (48)

whenever $t$ is a point of differentiability of $\overline{x}$. Let

$\displaystyle\mathcal{H}_{i}(t):=\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}H_{i}(t_{n}^{-})\geq
0,$ (49)

which is interpreted as the rate at which jobs are assigned to warm servers
with exactly $i$ jobs. Using Lemma 4 and (28), we have

$\displaystyle\dot{\overline{x}}_{0,2}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\left(X_{0,2}^{N_{k}}(t+\epsilon)-X_{0,2}^{N_{k}}(t)\right)=\overline{x}_{1,2}(t)-\mathcal{H}_{0}(t)+\beta\overline{x}_{0,1}(t)-\gamma\overline{x}_{0,2}(t)$
(50a) $\displaystyle\dot{\overline{x}}_{i,2}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\left(X_{i,2}^{N_{k}}(t+\epsilon)-X_{i,2}^{N_{k}}(t)\right)=\overline{x}_{i+1,2}(t)\mathbb{I}_{\\{i<B\\}}-\overline{x}_{i,2}(t)+\mathcal{H}_{i-1}(t)-\mathcal{H}_{i}(t)\mathbb{I}_{\\{i<B\\}}.$
(50b)

In the following, we need to show that
$\mathcal{H}_{i}(t)=h_{i}(\overline{x}(t))$ where the $h_{i}$’s are as in
Definition 1. We treat the cases of Power-of-$d$ and JBT-$d$ separately.

######  Lemma 5.

Assume that Power-of-$d$ is applied. Then, (47) and (48) hold true.

###### Proof.

If $\overline{x}_{0,0}+\overline{x}_{0,1}<1$, then the structure of the
$H_{i}$’s in (29) and Lemma 4 immediately give (47) and (48). Now, let us
assume that $\overline{x}_{0,0}+\overline{x}_{0,1}=1$. On coordinate (0,2), in
a point of differentiability we necessarily have $\dot{\overline{x}}_{0,2}=0$.
Using Lemma 4 and (28), we obtain

$\displaystyle\dot{\overline{x}}_{0,2}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\left(X_{0,2}^{N_{k}}(t+\epsilon)-X_{0,2}^{N_{k}}(t)\right)=\beta\overline{x}_{0,1}(t)-\mathcal{H}_{0}(t)=0.$
(51)

Similarly, on coordinate $(1,2)$, Lemma 4 and (51) imply that in a point of
differentiability we have
$\dot{\overline{x}}_{1,2}(t)=\mathcal{H}_{0}(t)-\mathcal{H}_{1}(t)=0$ and thus
$\mathcal{H}_{1}(t)=\mathcal{H}_{0}(t)=\beta\overline{x}_{0,1}(t)$. Then, on
coordinate $(i,2)$ by induction we obtain
$\mathcal{H}_{i}(t)=\mathcal{H}_{i-1}(t)=\beta\overline{x}_{0,1}(t)$. On the
other hand, we also have

$\displaystyle\dot{\overline{x}}_{0,2}(t)$
$\displaystyle=\beta\overline{x}_{0,1}(t)-\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{{W_{n}=0}\\}}H_{0}(t_{n}^{-})\geq\beta\overline{x}_{0,1}(t)-\lambda$
(52)

where in the last inequality we have just used that $H_{0}(t_{n}^{-})\leq 1$.
Thus, if $\beta\overline{x}_{0,1}(t)>\lambda$, we get a contradiction and $t$
can not be a point of differentiability. Substituting
$\mathcal{H}_{i}(t)=\beta\overline{x}_{0,1}(t)$ in (50) when
$\beta\overline{x}_{0,1}(t)\leq\lambda$, we obtain (47) and (48). ∎

The case of JBT-$d$ is more delicate than Power-of-$d$ because of the
discontinuous structure of the $H_{i}$’s when
$\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})=0$, see (30). In addition to a more
involved argument than the one presented in the proof of Lemma 5, which we
will develop in Lemma 7 below, we need the following lemma, which we will use
to determine an expression for $\mathcal{H}_{i}$ when
$\sum_{j=0}^{d}\overline{x}_{j,2}(t)=0$.

######  Lemma 6.

Assume that $\overline{x}(t)$ satisfies
$\overline{x}_{0,0}(t)+\overline{x}_{0,1}(t)<1$. Then, for all $i$

$\lim_{\epsilon\downarrow 0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbf{1}_{(Y_{i-1}^{N_{k}}(t_{n}^{-}),Y_{i}^{N_{k}}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N_{k}}(t_{n}^{-})-X_{0,1}^{N_{k}}(t_{n}^{-}))}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N_{k}}(t_{n}^{-})>0\\}}\\\
=\frac{\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N_{k}}(t_{n}^{-})>0\\}}.$
(53)

###### Proof.

Given in Section 8. ∎

The following lemma proves the desired property in the case of JBT-$d$.

######  Lemma 7.

Assume that JBT-$d$ is applied. Then, (47) and (48) hold true.

###### Proof.

We analyze $\mathcal{H}_{i}$ and the resulting expression will be substituted
in (50). This will give (47) and (48).

First, if $\overline{x}_{0,0}+\overline{x}_{0,1}=1$, the argument in the proof
of Lemma 5 gives i) $\mathcal{H}_{i}(t)=\beta\overline{x}_{0,1}(t)$ when
$\beta\overline{x}_{0,1}(t)\leq\lambda$ and ii) $t$ not a point of
differentiability when $\beta\overline{x}_{0,1}(t)>\lambda$. This gives (47)
and (48) (when $\overline{x}_{0,0}+\overline{x}_{0,1}=1$) and in the remainder
we assume that $\overline{x}_{0,0}+\overline{x}_{0,1}<1$.

Let us now assume that $\sum_{j=0}^{d}\overline{x}_{j,2}(t)>0$ and let
$\epsilon\in(0,\frac{\sum_{j=0}^{d}\overline{x}_{j,2}(t)}{2\phi(d+1)})$. Since
$t_{n}\in(t,t+\epsilon]$ whenever
$n\in\\{\mathcal{N}_{\phi}({N_{k}}t)+1,\ldots,\mathcal{N}_{\phi}({N_{k}}(t+\epsilon))\\}$,
(37) and the triangular inequality imply that for all $k$ sufficiently large
$|\sum_{j=0}^{d}X_{j,2}^{N_{k}}(t_{n})-\overline{x}_{j,2}(t)|\leq
2(d+1)\phi\epsilon<\sum_{j=0}^{d}\overline{x}_{j,2}(t)$ and thus
$\sum_{j=0}^{d}X_{j,2}^{N_{k}}(t_{n})>0$. We have shown that

$\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N_{k}}(t_{n}^{-})>0\\}}=1,\quad\forall
n\in\\{\mathcal{N}_{\phi}({N_{k}}t)+1,\ldots,\mathcal{N}_{\phi}({N_{k}}(t+\epsilon))\\}$
(54)

for all $k$ sufficiently large, given $\epsilon>0$ sufficiently small.
Substituting (54) in (30) and applying Lemma 4, we obtain (47) and (48) (under
the conditions $\overline{x}_{0,0}+\overline{x}_{0,1}<1$ and
$\sum_{j=0}^{d}\overline{x}_{j,2}(t)>0$).

It remains to understand the terms $\mathcal{H}_{i}$ in the case where
$\sum_{j=0}^{d}\overline{x}_{j,2}(t)=0$, which we assume in the remainder of
the proof.

Suppose that $t$ is a point of differentiability. Then, by applying Lemma 4 to
$X_{0,2}^{N}$ (see (28)), we obtain

$\displaystyle\dot{\overline{x}}_{0,2}(t)=\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\lim_{k\to\infty}X_{0,2}^{N_{k}}(t+\epsilon)-X_{0,2}^{N_{k}}(t)=\overline{x}_{1,2}(t)\mathbb{I}_{\\{d=0\\}}+\beta\overline{x}_{0,1}(t)-\mathcal{H}_{0}(t),$
(55)

and given that necessarily $\dot{\overline{x}}_{0,2}(t)=0$, we obtain

$\displaystyle\mathcal{H}_{0}(t)=\overline{x}_{1,2}(t)\mathbb{I}_{\\{d=0\\}}+\beta\overline{x}_{0,1}(t).$
(56)

Similarly, on coordinate $(i,2)$, with $0<i\leq d$, we obtain

$\displaystyle\dot{\overline{x}}_{i,2}(t)=\overline{x}_{i+1,2}(t)\mathbb{I}_{\\{i=d\\}}+\mathcal{H}_{i-1}(t)-\mathcal{H}_{i}(t)=0.$
(57)

By induction, this gives
$\mathcal{H}_{i}(t)=\mathcal{H}_{0}(t)=\beta\overline{x}_{0,1}(t)$ for all
$i<d$ and
$\mathcal{H}_{d}(t)=\beta\overline{x}_{0,1}(t)+\overline{x}_{d+1,2}(t)$, that
is,

$\displaystyle\mathcal{H}_{i}(t)=\beta\overline{x}_{0,1}(t)+\overline{x}_{d+1,2}(t)\mathbb{I}_{\\{i=d\\}},\qquad
i\leq d.$ (58)

We have proven (58) under the hypothesis that $t$ was a point of
differentiability but now we show that $\overline{x}(t)$ is not differentiable
if $\lambda<\overline{x}_{d+1,2}+(d+1)\beta\overline{x}_{0,1}$. Towards this
purpose, first we notice that

$\displaystyle\sum_{i=0}^{d}\mathcal{H}_{i}(t)=\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\sum_{i=0}^{d}H_{i}(t_{n}^{-})$
$\displaystyle\qquad=\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}\leq\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}=\lambda.$

Here, the first equality follows because the limits $\mathcal{H}_{i}(t)$ exist
and the second inequality follows by the fact that (recall the definition of
$\mathcal{H}_{i}$ in (49))

$\displaystyle\sum_{i=0}^{d}H_{i}(t_{n}^{-})=\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}+\mathbf{1}_{(0,Y_{d}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})=0\\}}$
(59)

and by Lemma 4 because $\sum_{j=0}^{d}\overline{x}_{j,2}(t)=0$. Then, using
(58), we necessarily have

$\displaystyle\sum_{i=0}^{d}\mathcal{H}_{i}(t)=\overline{x}_{d+1,2}(t)+(d+1)\beta\overline{x}_{0,1}(t)\leq\lambda$
(60)

and, given that necessarily $\mathcal{H}_{i}\geq 0$, we conclude that $t$ can
not be a point of differentiability whenever (60) does not hold true.

Now, we investigate $\mathcal{H}_{i}$ when $i>d$ and assuming that (60) holds
as otherwise $\overline{x}(t)$ would not be differentiable. We have

$\displaystyle\mathcal{H}_{i}(t)$
$\displaystyle=\frac{\lambda\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}-\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}$
$\displaystyle=\frac{\lambda\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}-\frac{\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=0\\}}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}$
$\displaystyle=\frac{\lambda\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}-\frac{\overline{x}_{i,2}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}\sum_{i=0}^{d}\mathcal{H}_{i}(t)$
$\displaystyle=\overline{x}_{i,2}(t)\frac{\lambda-\overline{x}_{d+1,2}(t)-(d+1)\beta\overline{x}_{0,1}(t)}{1-\overline{x}_{0,0}(t)-\overline{x}_{0,1}(t)}.$

In the first equality, we have used (30) and applied Lemma 4 to the definition
of $\mathcal{H}_{i}$ in (49); In the second, we have applied Lemma 6. In the
third, we have used (59) and that

$\displaystyle 0\leq\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\mathbf{1}_{(0,Y_{d}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})=0\\}}\leq\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\mathbf{1}_{(0,Y_{d}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}=0$

with the last inequality following by Lemma 4 as
$\sum_{j=0}^{d}\overline{x}_{j,2}(t)=0$; in the fourth, we have substituted
(58). This concludes the proof. ∎

Thus, we have shown that $\overline{x}$ is a fluid solution.

###  7.2 Proof of Theorem 2: Fixed points

We now prove Theorem 2. By definition, $x\in\mathcal{S}_{1}$ is a fixed point
if and only if

$\displaystyle 0$ $\displaystyle=\gamma{x}_{0,2}-\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (61a) $\displaystyle 0$ $\displaystyle=\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\beta{x}_{0,1}+\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (61b) $\displaystyle 0$ $\displaystyle=x_{1,2}-h_{0}(x)+\beta
x_{0,1}-\gamma x_{0,2}$ (61c) $\displaystyle 0$
$\displaystyle=x_{i+1,2}-x_{i,2}+h_{i-1}(x)-h_{i}(x),\quad i\geq 1.$ (61d)

Together with $\|x\|=1$, we now show that these conditions coincide with (14).

If i) $x_{0,0}=0$ and $\gamma{x}_{0,2}>\alpha g$, or if ii)
$x_{0,0}+x_{0,1}=1$, then we easily observe that $x$ cannot be a fixed point.
Therefore, in the following we exclude these conditions. Now, summing (61a)
and (61b), we obtain

$\beta x_{0,1}=\gamma x_{0,2}$ (62)

which gives (14b). Then, (14c) and (14d) directly follow from (61a) and (61b).

Substituting (62) in (61c), the conditions (61c)-(61d) become

$\displaystyle 0$ $\displaystyle=x_{1,2}-h_{0}(x)$ (63a) $\displaystyle 0$
$\displaystyle=x_{i+1,2}-x_{i,2}+h_{i-1}(x)-h_{i}(x),\quad i\geq 1,$ (63b)

and taking summations

$\displaystyle x_{i,2}=h_{i-1}(x),\quad i\geq 1.$ (64)

The equations in (63) are interpreted as the mean-field fixed-point equations
associated to Power-of-$d$ and JBT-$d$ when the number of servers is $Ny_{0}$
instead of $N$; we recall that $y_{i}=\sum_{i\geq 0}x_{i,2}$ is the proportion
of warm servers with at least $i$ jobs. Within Power-of-$d$, one can directly
check that for any given $x_{0,2}$, (64) holds if and only if $x_{i,2}$ is
given by (14e) and that, after a substitution, this gives $\sum_{i\geq
1}x_{i,2}=\lambda$ so that (14a) must hold true. The following lemma, given in
Section 8, handles the more delicate case of JBT-$d$.

######  Lemma 8.

Within JBT-$d$, for any given $x_{0,2}$, (63) holds if and only if $x_{i,2}$
satisfies (14f)-(14j). In addition, (14a) holds true.

Therefore, the conditions in (61) are equivalent to (14). This proves the
first statement of Theorem 2.

Now, under Assumption 3, $x_{i,2}$ is a function of $x_{0,2}$, for all $i\geq
1$, and we write $x_{i,2}$ as a shorthand notation for $x_{i,2}(x_{0,2})$.
Using (62), we can then focus only on the following conditions:

$\displaystyle x_{0,0}+\left(\frac{\gamma}{\beta}+1\right)x_{0,2}$
$\displaystyle=1-\lambda$ (65a) $\displaystyle\gamma x_{0,2}$
$\displaystyle\leq\alpha g,\quad\mbox{ if }x_{0,0}=0$ (65b)
$\displaystyle\gamma x_{0,2}$ $\displaystyle=\alpha g,\quad\mbox{ if
}x_{0,0}>0.$ (65c)

Here, we notice that
$(x_{0,0}^{\circ},x_{0,2}^{\circ})=\big{(}0,\frac{\beta}{\beta+\gamma}(1-\lambda)\big{)}$
uniquely solves (65) if

$\displaystyle\left(\frac{1}{\beta}+\frac{1}{\gamma}\right)\alpha
g(x^{\circ})\geq 1-\lambda$ (66)

where $x^{\circ}$ is uniquely determined by
$(x_{0,0}^{\circ},x_{0,2}^{\circ})$. So, let us assume that (66) does _not_
hold true. Then, if a point $(x_{0,0},x_{0,2})\in[0,1)^{2}$ that solves (65)
exists, then necessarily $x_{0,0}>0$ as otherwise $x_{0,2}=x_{0,2}^{\circ}$
(by (65a)) and (66) would hold, contradicting the hypothesis. This proves the
second part of Theorem 2.

###  7.3 Proof of Theorem 3: Fluid Optimality

We give a proof for Theorem 3. The non-linear structure taken by the $h_{i}$’s
when $x_{0,2}=0$, see (13), complicates the analysis and the identification of
a Lyapunov function. For this reason, our strategy is based on a divide-and-
conquer approach. This will actually provide insights about the dynamics
followed by fluid solutions. For simplicity, we provide a proof assuming that
$B<\infty$, which is essentially equivalent to assume that $x_{i,2}(0)=0$ for
all $i$ large enough; this is not critical as $x_{i,2}^{\star}=0$ for all
$i\geq 2$.

Let $\overline{Q}(x):=\sum_{i=1}^{B}ix_{i,2}$, i.e., the overall number of
jobs in the system in state $x$. The following lemma gives a property on the
time derivative of $\overline{Q}(x(t))$.

######  Lemma 9.

Let $x(t)$ be a fluid solution induced by JIQ such that
$x(0)\in\mathcal{S}_{1}$ and $B<\infty$. If $t$ is a point of
differentiability, then

$\dot{\overline{Q}}(x(t))=\lambda-y_{1}(t).$ (67)

###### Proof.

First, we notice that

$\displaystyle\dot{\overline{Q}}(x(t))=\sum_{i\geq
1}i\dot{x}_{i,2}(t)=-y_{1}+\sum_{i=0}^{B-1}h_{i}(x(t))$ (68)

where the second equality follows by applying Definition 1. Now, we treat the
cases $x_{0,2}(t)>0$ and $x_{0,2}(t)=0$ separately. Suppose that
$x_{0,2}(t)>0$. Then, $h_{i}(x(t))=\lambda\mathbb{I}_{\\{i=0\\}}$ (by (13))
and substituting in (68) we immediately get $\dot{\overline{Q}}(x(t))=\lambda-
y_{1}(t)$ as desired. Thus, suppose in the remainder that $x_{0,2}(t)=0$. Now,
assume that $y_{0}(t)>0$. Then, using again (13),

$\displaystyle\dot{\overline{Q}}(x(t))$
$\displaystyle=-y_{1}+\left(\beta{x}_{0,1}+{x}_{1,2}\right)\mathbb{I}_{\\{{x}_{1,2}+\beta{x}_{0,1}\leq\lambda\\}}+\frac{y_{1}}{y_{0}}(\lambda-{x}_{1,2}-\beta{x}_{0,1})^{+}$
$\displaystyle=-y_{1}+\left(\beta{x}_{0,1}+{x}_{1,2}\right)\mathbb{I}_{\\{{x}_{1,2}+\beta{x}_{0,1}\leq\lambda\\}}+(\lambda-{x}_{1,2}-\beta{x}_{0,1})^{+}$

and the statement follows immediately if
${x}_{1,2}(t)+\beta{x}_{0,1}(t)\leq\lambda$. On the other hand, if
${x}_{1,2}(t)+\beta{x}_{0,1}(t)>\lambda$, then, since $x_{0,2}(t)=0$ and $t$
is supposed to be a point of differentiability, we get (by (5c)) the
contradiction that $0=\dot{x}_{0,2}(t)=x_{1,2}(t)-h_{0}(x(t))+\beta
x_{0,1}(t)-\gamma x_{0,2}(t)=x_{1,2}(t)+\beta x_{0,1}(t)>\lambda$; the first
equality holds because $x_{0,2}(t)$ is a non-negative absolutely continuous
function. This shows that $t$ cannot be a point of differentiability. Finally,
if $y_{0}(t)=0$, then the differentiability at $t$ and the normalizing
condition $\|x\|=1$ give
$0=\dot{y}_{0}(t)=-\dot{x}_{0,0}(t)-\dot{x}_{0.1}(t)=\beta x_{0,1}$ and thus
$x_{0,0}(t)=1$. Assumption 2 requires that $g(x)>0$ when $x_{0,0}=1$, so (5a)
implies that $\dot{{x}}_{0,0}<0$. This contradicts that $t$ is a point of
differentiability because $x_{0,0}(t)$ is uniformly bounded by one and
absolutely continuous. ∎

We now prove Theorem 3 by showing that $\|x(t)-x^{\star}\|\to 0$ in each of
the following complete and mutually exclusive cases. For each case, we show
that $x(t)$ follows a unique trajectory that stays in $\mathcal{S}_{1}$.

_Case i)_. Suppose that $x_{0,2}(t)=0$ for all $t\geq 0$. This rules out the
possibility that $x_{0,0}(t)$ stays on zero for all $t$ large enough because
(5a) and (5b), together with the normalizing condition $\|x\|=1$, would imply
that $y_{1}(t)\to 1$ as $t\to\infty$, and in this case Lemma 9 yields the
contradiction that $\overline{Q}(x(t))$ is eventually negative. Thus, without
loss of generality, let us assume that $x_{0,0}(0)>0$. Then, using (5), $x(t)$
satisfies

$\displaystyle\dot{{x}}_{0,0}$ $\displaystyle=-\alpha g(x)$ (69a)
$\displaystyle\dot{{x}}_{0,1}$ $\displaystyle=\alpha g(x)-\beta{x}_{0,1}$
(69b) $\displaystyle\dot{x}_{0,2}$ $\displaystyle=0,\quad x_{1,2}+\beta
x_{0,1}\leq\lambda.$ (69c)

Note that $\lim_{t\to\infty}x_{0,0}(t)$ exists, say $x_{0,0}(\infty)$, because
$\dot{x}_{0,0}(t)\leq 0$ and $x_{0,0}(t)$ is uniformly bounded. Thus, as
$t\to\infty$, $\dot{x}_{0,0}(t)=-\alpha g(x(t))\to 0$. Given the assumptions
on $g$, $(\lambda-x_{0,1}(t)-\beta x_{1,2}(t))^{+}\to 0$ and since
$x_{1,2}(t)+\beta x_{0,1}(t)\leq\lambda$ for all $t$, by (69c), we obtain that
$x_{1,2}(t)+\beta x_{0,1}(t)\to\lambda$. Then, (69b) and $g(x(t))\to 0$ imply
that $x_{0,1}(t)\to 0$ and thus $x_{1,2}(t)\to\lambda$. In turn, (5d) gives
$x_{i,2}(t)\to 0$ for all $i\geq 2$, and the normalizing condition $\|x\|=1$
implies that necessarily $x_{0,0}(t)\to 1-\lambda$. Thus,
$\|x(t)-x^{\star}\|\to 0$.

_Case ii)_. Suppose that $x_{0,2}(t)>0$ for all $t$. Then, $x(t)$ satisfies
the following conditions (using Definition 1)

$\displaystyle\dot{{x}}_{0,0}$ $\displaystyle=\gamma{x}_{0,2}-\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (70a) $\displaystyle\dot{{x}}_{0,1}$ $\displaystyle=\alpha
g\mathbb{I}_{\\{{x}_{0,0}>0\\}}-\beta{x}_{0,1}+\gamma{x}_{0,2}\,\mathbb{I}_{\\{{x}_{0,0}=0,\,\gamma{x}_{0,2}\leq\alpha
g\\}}$ (70b) $\displaystyle\dot{x}_{0,2}$ $\displaystyle=x_{1,2}-\lambda+\beta
x_{0,1}-\gamma x_{0,2}$ (70c) $\displaystyle\dot{x}_{1,2}$
$\displaystyle=x_{2,2}-x_{1,2}+\lambda$ (70d) $\displaystyle\dot{x}_{i,2}$
$\displaystyle=x_{i+1,2}\mathbb{I}_{\\{i<B\\}}-x_{i,2},\quad i\geq 2.$ (70e)

The ODE system (70d)-(70e) is an autonomous linear ODE system with constant
coefficients and, developing the matrix-exponential general solution of such
ODE system, for all $i\geq 1$ we obtain

$\displaystyle
x_{i,2}(t)=\lambda\mathbb{I}_{\\{i=1\\}}+e^{-t}\sum_{k=i}^{B}\frac{t^{k-i}}{(k-i)!}(x_{k,2}(0)-\lambda\mathbb{I}_{\\{k=1\\}})$
(71)

and thus $x_{i,2}(t)\to\lambda\mathbb{I}_{\\{i=1\\}}$ as $t\to\infty$. In
turn, $\lim_{t\to\infty}(\lambda-x_{1,2}(t)-\beta
x_{0,1}(t))^{+}=\lim_{t\to\infty}(-\beta x_{0,1}(t))^{+}=0$ and therefore
$g(x(t))\to 0$. Since $g(x(t))\to 0$, $x_{0,1}(t)\to 0$ necessarily by (70b),
and using this in (70c) we obtain $x_{0,2}(t)\to 0$. Since $\|x\|=1$,
necessarily $x_{0,0}(t)\to 1-\lambda$ and we have shown that
$\|x(t)-x^{\star}\|\to 0$.

_Case iii)_. If the conditions in cases _i)_ and _ii)_ are not met, then there
exists $t_{0}$, $t_{1}$, with $t_{0}\leq t_{1}<\infty$, and $\delta>0$ such
that

  1. 1.

$x_{0,2}(t)=0$ for all $t\in[t_{0},t_{1}]$

  2. 2.

$x_{0,2}(t)>0$ and $\dot{x}_{0,2}(t)<0$ for all $t\in[t_{0}-\delta,t_{0})$,
and

  3. 3.

$x_{0,2}(t)>0$ and $\dot{x}_{0,2}(t)>0$ for all $t\in(t_{1},t_{1}+\delta]$.

On $[t_{0}-\delta,t_{0})$, $h_{0}(x(t))=\lambda$ (by (13)) and using (5c), we
obtain $\dot{x}_{0,2}(t)=x_{1,2}(t)-\lambda+\beta x_{0,1}(t)-\gamma
x_{0,2}(t)<0$ and thus by continuity

$\displaystyle 0\geq\lim_{t\uparrow t_{0}}x_{1,2}(t)-\lambda+\beta
x_{0,1}(t)-\gamma x_{0,2}(t)=x_{1,2}(t_{0})-\lambda+\beta x_{0,1}(t_{0}).$
(72)

Since $x_{0,2}(t_{0})=0$ on $[t_{0},t_{1}]$, (5c) implies that (72) holds as
well on $[t_{0},t_{1}]$. On $(t_{1},t_{1}+\delta]$, $h_{0}(x(t))=\lambda$ (by
(13)) and using again (5c), we obtain

$\displaystyle 0<\dot{x}_{0,2}(t)=x_{1,2}(t)-h_{0}(x(t))+\beta
x_{0,1}(t)-\gamma x_{0,2}(t)<x_{1,2}(t)-\lambda+\beta x_{0,1}(t)$ (73)

and therefore $g(x(t))=0$. By continuity of fluid solutions,
$x_{1,2}(t_{1})+\beta x_{0,1}(t_{1})=\lambda$. In addition, on
$(t_{1},t_{1}+\delta]$, $x(t)$ is uniquely defined by

$\displaystyle\dot{{x}}_{0,0}$ $\displaystyle=\gamma{x}_{0,2}$ (74a)
$\displaystyle\dot{{x}}_{0,1}$ $\displaystyle=-\beta{x}_{0,1}$ (74b)
$\displaystyle\dot{x}_{0,2}$ $\displaystyle=x_{1,2}-\lambda+\beta
x_{0,1}-\gamma x_{0,2}$ (74c) $\displaystyle\dot{x}_{1,2}$
$\displaystyle=x_{2,2}-x_{1,2}+\lambda$ (74d) $\displaystyle\dot{x}_{i,2}$
$\displaystyle=x_{i+1,2}\mathbb{I}_{\\{i<B\\}}-x_{i,2},\quad i\geq 2,$ (74e)

and we also know that $\dot{x}_{0,2}(t)>0$. As long as _a)_ $g(x(t))=0$ and
_b)_ $x_{0,2}(t)>0$, on $[t_{1},\infty)$ the fluid solution under
investigation $x(t)$ is indeed uniquely given by the trajectory induced by
(74) on $[t_{1},\infty)$. In the remainder, we show that both _a)_ and _b)_
hold true for all $t$. This will conclude the proof because $x^{\star}$ is the
unique fixed point of (74) and because (74) is a linear ODE system with
constant coefficients. For simplicity of notation, let us shift time and
assume that $t_{1}=0$. Now, since $x_{0,1}(t)=x_{0,1}(0)e^{-\beta t}$ (by
(74b)) and since $x_{1,2}(t)$ takes the form given in (71), substituting in
(74c) we obtain

$\displaystyle\dot{x}_{0,2}(t)$ $\displaystyle=\beta x_{0,1}(0)e^{-\beta
t}-\gamma
x_{0,2}(t)+e^{-t}(x_{1,2}(0)-\lambda)+e^{-t}\sum_{k=2}^{B}\frac{t^{k-1}}{(k-1)!}x_{k,2}(0)$
$\displaystyle=\beta x_{0,1}(0)e^{-\beta t}-\gamma x_{0,2}(t)-\beta
x_{0,1}(0)e^{-t}+e^{-t}\sum_{k=2}^{B}\frac{t^{k-1}}{(k-1)!}x_{k,2}(0)$
$\displaystyle\geq\beta x_{0,1}(0)(e^{-\beta t}-e^{-t})-\gamma x_{0,2}(t).$

Thus, $x_{0,2}(t)\geq z(t)$ where $z(t)$ is uniquely defined by
$\dot{z}(t)=\beta x_{0,1}(0)(e^{-\beta t}-e^{-t})-\gamma z(t)$ with
$z(0)=x_{0,2}(0)$. The solution of this differential equation is

$\displaystyle z(t)=\beta x_{0,1}(0)e^{-\gamma
t}\left(\frac{1-e^{-t\left(\beta-\gamma\right)}}{\beta-\gamma}-\frac{1-e^{-t\left(1-\gamma\right)}}{1-\gamma}\right)$

and now we notice that $z(t)>0$ if $\beta>1$, for all $t$. This proves
property _b)_. To prove property _a)_ , we use again (71) and
$x_{1,2}(t_{1})+\beta x_{0,1}(t_{1})=\lambda$ to obtain

$\displaystyle x_{1,2}(t)+\beta x_{0,1}(t)-\lambda=\beta
x_{0,1}(0)\left(e^{-\beta
t}-e^{-t}\right)+e^{-t}\sum_{k=2}^{B}\frac{t^{k-1}}{(k-1)!}x_{k,2}(0)>0,$ (75)

where the last inequality follows because $\beta<1$. Given (23), (75) implies
$g(x(t))=0$.

##  8 Proofs of Technical Lemmas

###  8.1 Proof of Lemma 4

We give a proof for the first limit because the argument used for the others
is identical.

Since $t_{n}\in(t,t+\epsilon]$ whenever
$n\in\\{\mathcal{N}_{\phi}({N_{k}}t)+1,\ldots,\mathcal{N}_{\phi}({N_{k}}(t+\epsilon))\\}$,
(37) implies that for all $k$ sufficiently large
$|Y_{i}^{N_{k}}(t_{n})-\sum_{j=0}^{i}\overline{x}_{j,2}(t)|\leq C\epsilon$,
for some constant $C$, i.e.,

$\displaystyle\mathbf{1}_{(\sum_{j=0}^{i}\overline{x}_{j,2}(t)+C\epsilon,\sum_{j=0}^{i}\overline{x}_{j,2}(t)-C\epsilon]}^{D_{n}}\leq\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{D_{n}}\leq\mathbf{1}_{(\sum_{j=0}^{i}\overline{x}_{j,2}(t)-C\epsilon,\sum_{j=0}^{i}\overline{x}_{j,2}(t)+C\epsilon]}^{D_{n}}$
(76)

Let $\Gamma$ denote the LHS of the first equation in Lemma 4. Applying Lemma
1, we obtain

$\displaystyle\Gamma$ $\displaystyle\leq\lim_{\epsilon\downarrow
0}\lim_{k\to\infty}\frac{1}{\epsilon
N_{k}}\sum_{n=\mathcal{N}_{\phi}(N_{k}t)+1}^{\mathcal{N}_{\phi}(N_{k}(t+\epsilon))}\mathbb{I}_{\\{W_{n}=1\\}}\mathbf{1}_{(\sum_{j=0}^{i}\overline{x}_{j,2}(t)-C\epsilon,\sum_{j=0}^{i}\overline{x}_{j,2}(t)+C\epsilon]}^{D_{n}}=\overline{x}_{i,2}(t)$

and using (76) in the other direction we obtain $\Gamma=\overline{x}_{i,2}(t)$
as desired.

###  8.2 Proof of Lemma 6

We recall that we have analyzed $\overline{x}$ along a fixed
$\omega\in\mathcal{C}$, where $\mathbb{P}(\mathcal{C})=1$. We now explicit the
dependence on $\omega$ and treat quantities $\overline{x}(t)$ and $X^{N}(t)$
as random variables. Let

$\displaystyle
Z_{n}^{N}:=\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N_{k}}(t_{n}^{-}))}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}.$
(77)

For all $n$, the random variable $Z_{n}^{N}$ is $\mathcal{F}_{n}$-measurable
where $\mathcal{F}_{n}:=\\{X^{N}(t_{n}^{N,\lambda-}),A_{n}^{1},W_{n}\\}$, and

$\displaystyle\mathbb{E}[Z_{n}^{N}|\mathcal{F}_{n}\setminus
A_{n}^{1}]=\frac{X_{i,2}^{N}(t_{n}^{-})}{1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-})}\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}$
(78)

where the set $\mathcal{F}_{n}\setminus W_{n}$ denotes the set
$\mathcal{F}_{n}$ with $A_{n}^{1}$ removed. Now, let
$\Delta_{n}^{N}:=Z_{n}^{N}-\mathbb{E}[Z_{n}^{N}|\mathcal{F}_{n}\setminus
W_{n}]$. Then, $\mathbb{E}[\Delta_{n}^{N}|\mathcal{F}_{n}\setminus W_{n}]=0$
and $|\Delta_{n}^{N}|\leq 2$, and applying the Azuma–Hoeffding inequality, we
get

$\mathbb{P}\left(\frac{1}{N}\left|\sum_{n=1}^{N}\Delta_{n}^{N}\right|>\delta\right)\leq
2\exp\left(-\frac{(N\delta)^{2}}{8N}\right)$ (79)

for any $\delta>0$. Since $\sum_{N}\exp\left(-{N\delta^{2}/8}\right)<\infty$,
an application of the Borel–Cantelli lemma shows that
$\frac{1}{N}\sum_{n=1}^{N}\Delta_{n}^{N}\to 0$ almost surely. In particular,

$\displaystyle\lim_{N\to\infty}\frac{1}{\epsilon N}\sum_{{\scriptstyle
n=\mathcal{N}_{\phi}(Nt)+1:}\atop{\scriptstyle
W_{n}=0}}^{\mathcal{N}_{\phi}(N(t+\epsilon))}\mathbf{1}_{(Y_{i-1}^{N}(t_{n}^{-}),Y_{i}^{N}(t_{n}^{-})]}^{A_{n}^{1}(1-X_{0,0}^{N}(t_{n}^{-})-X_{0,1}^{N}(t_{n}^{-}))}\left(\mathbb{I}_{\\{\sum_{j=0}^{d}X_{j,2}^{N}(t_{n}^{-})>0\\}}-\mathbb{E}[Z_{n}^{N}|\mathcal{F}_{n}\setminus
W_{n}]\right)=0$ (80)

almost surely. We now come back to work on a given trajectory $\omega$. In
view of the previous equality, we may redefine $\mathcal{C}$ in Lemma 1 to be
a subset of $\mathcal{C}^{\prime}$ where $\mathbb{P}(\mathcal{C}^{\prime}=1)$
and (80) holds for all $\omega\in\mathcal{C}^{\prime}$. Therefore, we fix
$\omega\in\mathcal{C}$ and use (37) and (78) to obtain that

$\lim_{k\to\infty}\frac{1}{\epsilon N_{k}}\sum_{{\scriptstyle
n=\mathcal{N}_{\phi}(N_{k}t)+1:}\atop{\scriptstyle