(

# Near-Optimal Stochastic Bin-Packing in Large Service Systems with Time-
Varying Item Sizes

Yige Hong1 Qiaomin Xie2 Weina Wang1
1Carnegie Mellon University 2University of Wisconsin-Madison
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Motivated by the virtual machine scheduling problem in today’s computing
systems, we propose a new setting of stochastic bin-packing in service systems
that allows the item sizes (job resource requirements) to vary over time. In
this setting, items (jobs) arrive to the system, vary their sizes, and depart
from the system following certain Markovian assumptions. We focus on
minimizing the expected number of non-empty bins (active servers) in steady
state, where the expectation in steady state is equal to the long-run time-
average with probability $1$ under the Markovian assumptions. Our main result
is a policy that achieves an optimality gap of $O(\sqrt{r})$ in the objective,
where the optimal objective value is $\Theta(r)$ and $r$ is a scaling factor
such that the item arrival intensity scales linearly with it. When specialized
to the setting where the item sizes do not vary over time, our result improves
upon the state-of-the-art $o(r)$ optimality gap. Our technical approach
highlights a novel policy conversion framework that reduces the policy design
problem to that in a single-bin (single-server) system.

##  1 Introduction

The bin-packing problem is a classical problem that has been extensively
studied in theoretical computer science and operations research (see, e.g.,
the survey paper Coffman et al. 1984). In a bin-packing problem, we are given
$T$ items with sizes $s_{1},s_{2},\dots,s_{T}$, and the goal is to pack them
into as few as possible bins where each bin can hold items whose total size
does not exceed $1$. A stochastic version of the problem, referred to as
_stochastic bin-packing (SBP)_ , assumes that items arrive one by one and
their sizes are random variables that are independent and identically
distributed (i.i.d.) (Courcoubetis and Weber 1986, Courcobetis and Weber 1990,
Csirik et al. 2006, Gupta and Radovanović 2020, Freund and Banerjee 2019,
Ayyadevara et al. 2022). For stochastic bin-packing, an online algorithm is
usually compared against the optimal offline algorithm OPT that sees the
realizations of all the item sizes. The optimality gap of the online algorithm
is then measured by the difference in the expected numbers of bins used under
the algorithm and under OPT. A fundamental question for stochastic bin-packing
is: what optimality gap is achievable as $T$ becomes large? The answer to this
question can depend on the regime that the setting falls into. It has been
shown that there are three possible regimes based on the expected wasted space
under OPT (Courcoubetis and Weber 1986, Courcobetis and Weber 1990): Linear
Waste (LW) for $\Theta(T)$ waste, Perfectly Packable (PP) for
$\Theta(\sqrt{T})$ waste, and Bounded Waste (BW) for $O(1)$ waste. However,
the recent work of Freund and Banerjee (2019) develops an algorithm that
achieves an $O(1)$ optimality gap across all three regimes.

Motivated by the problem of scheduling virtual machines in large-scaling
computing systems, stochastic bin-packing in large _service systems_ has been
proposed (Stolyar 2013, Stolyar and Zhong 2013, 2015, Stolyar 2017, Stolyar
and Zhong 2021, Ghaderi et al. 2014). Due to the application background, this
formulation is usually presented in terms of _jobs and servers, where jobs
correspond to items, whose sizes are the amounts of resources they need, and
servers correspond to bins_. Since our paper focuses on service systems, we
will use the terminology of jobs and servers from this point onward. In
service systems, jobs not only arrive over time but also _depart_ over time.
Moreover, a further deviation from the conventional SBP setting is the
Markovian assumptions on job dynamics in service systems. In particular, jobs
are assumed to arrive according to Poisson processes, i.e., the interarrival
times between two consecutive jobs are all i.i.d. following an exponential
distribution. For the departure dynamics, each job is assumed to stay in the
system for an exponentially distributed service time. The service time of a
job is unknown until the job departs. These assumptions make the Markovian
representation of the system state succinct. They have a long history in
stochastic modeling of service systems such as call centers and computing
clusters (e.g., see the literature for stochastic systems and queueing
theory); and they strike a good balance between providing relevant insights
for practical systems and being mathematically tractable. We note that job
departures are also handled in the extensively studied dynamic bin-packing
problem (e.g. Coffman et al. 1983, Li et al. 2014, Buchbinder et al. 2021, and
references therein); however, there, no stochastic assumptions are made and
the emphasis is on worst-case analysis.

For SBP in service systems, when a job arrives, the scheduler uses a policy to
determine where to place the job, either to an active server (non-empty bin)
or an idle server (empty bin), based on all the system state information up to
the current time. The goal is to design a policy $\sigma$ to minimize the
expected number of active servers in steady state, denoted as $N(\sigma)$.
Here a steady state exists because of the Markovian assumptions, and the
expectation of a random quantity in steady state is equal to its long-run
time-average with probability $1$. With this formulation, it is more natural
to compare a policy $\sigma$ with the _optimal policy $\sigma^{*}$_, rather
than with algorithms that are offline or have access to information in the
future. The optimality gap of a policy $\sigma$ is $N(\sigma)-N(\sigma^{*})$.

Since SBP in service systems typically considers an infinite time horizon and
aims to model today’s large-scale computing systems, the optimality gap of a
policy is usually studied in the regime where the arrival rate/intensity of
jobs becomes large. Let the total arrival rate of jobs scale linearly with a
scaling factor $r$. Then the optimal value $N(\sigma^{*})$ can be shown to be
$\Theta(r)$. Therefore, we say a policy is asymptotically optimal if its
optimality gap is $o(r)$. Stolyar (2013), Stolyar and Zhong (2013) proposed
greedy policies that are asymptotically optimal, but the scheduler that
executes these policies needs to know detailed state information, which is in
a high-dimensional space. Stolyar and Zhong (2015, 2021) later developed
policies that use much less state information and achieve $\Theta(r)$ (with
arbitrarily small constant) and $o(r)$ optimality gaps, respectively.

As one can see from the literature above, we have developed a good
understanding of the problem of SBP in large service systems thus far.
However, when examining the application of scheduling virtual machines (or
containers in many modern cloud computing systems) (Delimitrou and Kozyrakis
2014, Cloud 2023, Foundation 2023), we find that there is an essential
characteristic of practical systems that cannot be captured by the current
formulation. That is, the actual resource requirement of a job (item size in
bin-packing) _varies over time_. For example, in computing clusters that host
user-facing services, a virtual machine job serves requests from users. The
actual CPU and memory usage of such a job depends on the instantaneous demand
of requests, which varies widely during a day (Delimitrou and Kozyrakis 2014).
To achieve better server utilization under time-varying resource requirements,
the system does not schedule jobs based on their peak requirements, but rather
allows temporary resource overcommitment/contention, i.e., allows the total
resource requirement on a server to temporarily exceed its capacity. Given the
above flexibility in the capacity constraint, a key question for the
scheduling problem is how to best use it to maximize server utilization.
Currently, data centers typically suffer from a low utilization level, which
results in increased operational costs and energy footprint, and has become
one of the key limiting factors for the continued scaling of today’s computing
systems (Reiss et al. 2012, Delimitrou and Kozyrakis 2014, Lo et al. 2015,
Bashir et al. 2021).

Motivated by this background, in this paper, we propose a _new setting of
stochastic bin-packing in service systems_ that allows the resource
requirement of a job (i.e., size of an item) to vary over time. We then
capture the flexibility in the capacity constraint by setting a budget for the
cost associated with resource contention. Most of the remainder of this
section will be devoted to an overview of our work. We will begin by
describing our formulation, providing a simplified version for ease of
exposition. We will then state our main result: a proposed policy that
achieves an optimality gap of $O(\sqrt{r})$. This result can be specialized to
the setting where the resource requirement of a job does not vary over time,
improving upon the state-of-the-art $o(r)$ optimality gap result. We will also
provide an overview of our technical approach, which highlights a novel policy
conversion framework.

###  1.1 Overview of Our Work

#### Problem Formulation: Simplified Version

(a) A simplified version of our job model. Each job in service is in either an
$L$ phase or an $H$ phase, associated with low and high resource requirements,
respectively. When the job is completed, it is said to be in the state
$\perp$. The job transitions between the two phases while in service until it
is completed, following a continuous-time Markov chain with rates
$\mu_{ii^{\prime}}$, $i,i^{\prime}\in\\{L,H,\perp\\}$.

(b) A system model with an infinite number of identical servers. As soon as a
job arrives to the system, the job needs to be dispatched to a server to start
service immediately. The configuration of each server is the number of jobs in
each phase on the server.

Figure 1: Job model and system model.

We first describe our job model that features time-varying resource
requirements. For ease of exposition, here we present a setting where each job
in service can be in one of the two _phases_ , $L$ and $H$, associated with
_low_ and _high_ resource requirements, respectively. Our full model,
presented in Section 2, allows more than two phases and more than one type of
resources. To model the temporal variation in the resource requirement, we
assume that each job transitions between the two phases while in service until
it is completed, following a continuous-time Markov chain illustrated in
Figure 1(a). We use an absorbing state $\perp$ to denote that the job is
completed. A job can initialize in either phase $L$ or phase $H$, and they are
referred to as _type_ $L$ and _type_ $H$ jobs, respectively. Note that the
setting where a job’s resource requirement does not vary over time is a
special case of our job model where the transition rates between phases are
$0$.

We consider a system with an infinite number of identical servers, illustrated
in Figure 1(b). For job arrivals, we assume the same Poisson assumption as in
the existing work on SBP in service systems. In particular, we assume that the
two types of jobs arrive to the system following two independent Poisson
processes, with rates $\Lambda_{L}$ and $\Lambda_{H}$, respectively; i.e., the
interarrival times of type $L$ and type $H$ jobs are i.i.d., following
exponential distributions with means $1/\Lambda_{L}$ and $1/\Lambda_{H}$,
respectively. Upon arrival, a job needs to be dispatched to a server according
to a _dispatch policy_ , and the job enters service immediately. The goal is
to design a policy to minimize the expected number of _active servers_
(servers currently serving a positive number of jobs) in steady state, denoted
as $N(\sigma)$.

To quantify the cost of resource contention on a server, we first represent
the state of a server by its _configuration_ , a vector
${\bm{k}}=(k_{L},k_{H})$ where $k_{L}$ and $k_{H}$ are the numbers of jobs in
phase $L$ and phase $H$, respectively. Then a _cost rate function_ $h(\cdot)$
maps a server’s configuration to a rate of cost. Typically, the cost rate is
positive when the total resource requirement of all jobs on a server exceeds
the server’s resource capacity, although our assumptions allow more general
forms of $h(\cdot)$. Here we assume that the resource contention does not
affect the transition rates in the job model nor prompt jobs to be killed,
suitable for the application scenarios where contention level is low and
manageable. Let $C(\sigma)$ denote the average expected cost rate per server.

Now our bin-packing problem can be formulated as follows:

$\displaystyle\underset{\sigma}{\text{minimize}}$ $\displaystyle N(\sigma)$
(1) subject to $\displaystyle C(\sigma)\leq\epsilon,$

where $\epsilon$ is a budget for the cost of resource contention. We are
interested in solving this problem in the asymptotic regime where the arrival
rates $(\Lambda_{L},\Lambda_{H})$ scale to infinity (Stolyar and Zhong 2013,
2015, Stolyar 2017, Stolyar and Zhong 2021), motivated by the ever-increasing
computing demand that drives today’s computing systems to be large-scale.
Specifically, we assume
$(\Lambda_{L},\Lambda_{H})=(\lambda_{L}r,\lambda_{H}r)$ for some fixed
coefficients $\lambda_{L}$ and $\lambda_{H}$ and a scaling factor
$r\to\infty$.

#### Main Result

We propose a meta-policy named Join-Requesting-Server (JRS) for dispatching
jobs in the considered system. JRS can take different subroutines as input,
and we develop a specific subroutine named Single-OPT. We give a brief
description of JRS and Single-OPT in the Technical Approach below, and leave
the detailed description to Section 4 and Appendix C. We prove that JRS with
the Single-OPT subroutine is _asymptotically optimal_ , in the sense that the
expected number of active servers is at most
$\left(1+O\left(r^{-0.5}\right)\right)$ times the optimal value of the
optimization problem in (1), while the cost rate incurred is at most
$\left(1+O\left(r^{-0.5}\right)\right)\cdot\epsilon$. It can be shown that the
optimal value of (1) is $\Theta(r)$, so the resultant additive optimality gap
is $O(\sqrt{r})$. Note that the policy does violate the cost constraint by a
diminishing fraction. This main result will be formally presented in Theorem
1.

When we specialize our result to the setting where the resource requirement of
a job does not vary over time, we remove the constraint
$C(\sigma)\leq\epsilon$ and consider a configuration to be feasible if the
total resource requirement is within the capacity. Then we will get a policy
that has an $O(\sqrt{r})$ optimality gap in the objective value.

#### Technical Approach

Figure 2: A single-server system with an infinite supply of jobs. A single-
server policy decides when to request jobs and how many jobs of each type to
request.

We develop a novel _policy-conversion_ framework that reduces the optimal-
policy design problem in (1) to an optimal-policy design problem in a single-
server system. Specifically, this single-server system, as illustrated in
Figure 2, has an infinite supply of jobs of all types. As a result, the server
can start the service of any number of new jobs of any type at any time. When
the server starts the service of a new job, we say that the server _requests_
the job (from the infinite supply). We assume that the server is also subject
to a resource contention cost, following the same cost model as that in the
original infinite-server system. Then the optimal-policy design problem in
this single-server system is as follows. We aim to find a policy that decides
when to request jobs of each type, such that the system throughput is equal to
$(\Lambda_{L}/\overline{N},\Lambda_{H}/\overline{N})$ for some number
$\overline{N}$ as small as possible, while the cost rate is kept within a
budget $\epsilon$. The objective $\overline{N}$ can be interpreted as the
copies of the single-server system needed so that when they run in parallel,
their total throughput matches the arrival rates $(\Lambda_{L},\Lambda_{H})$
in the infinite-server system. For ease of notation, we simply write
$\overline{N}$ to mean $\lceil\overline{N}\rceil$, the smallest integer number
larger than or equal to $\overline{N}$. Let $\overline{N}^{*}$ be the optimal
value for $\overline{N}$. Our approach has the following three steps.

##### Lower Bound.

We first show that the expected number of active servers in the original
infinite-server system is lower bounded by $\overline{N}^{*}$. This result is
formally presented in Theorem 2. The proof is carried out using a coupling
argument. Intuitively, the single-server system can always be made more
efficient than an “average” server in the infinite-server system, since the
former can mimic the dynamics of any server in the latter system.

##### Converting a Single-Server Policy to an Infinite-Server Policy.

Next, we come to the _key_ part of our framework, where we use the proposed
meta-policy JRS to convert a policy $\overline{\sigma}$ in the single-server
system to a dispatch policy in the original infinite-server system. We say
that JRS takes $\overline{\sigma}$ as a subroutine. Under this conversion, we
prove that the $\overline{N}$ under $\overline{\sigma}$ corresponds to the
expected number of active servers under the converted policy in the infinite-
server system, and that the converted policy conforms to the budget
$\epsilon$, all in an asymptotic sense. This result is presented in Theorem 3,
which we refer to as the _conversion theorem_. With the conversion theorem,
designing an asymptotically optimal policy in the original infinite-server
system _reduces_ to designing a policy in the single-server system that
achieves the optimal value $\overline{N}^{*}$.

We outline this conversion procedure of JRS below. Fix a policy
$\overline{\sigma}$ for the single-server system. We let each server $\ell$
with $\ell\in\\{1,2,\dots,\overline{N}\\}$ in the infinite-server system adopt
$\overline{\sigma}$, based on the current configuration of server $\ell$.
However, when $\overline{\sigma}$ requests a job in the single-server policy,
server $\ell$ cannot obtain a job immediately in the infinite-server system.
In this case, server $\ell$ generates a _token_ and sends it to the dispatcher
to store the request. The type of the token is set to be the type of the
requested job. Then when a job of type $i$ (e.g., $i\in\\{L,H\\}$) arrives to
the system, the dispatcher dispatches the job based on the tokens of type $i$.
In particular, the dispatcher chooses a type $i$ token uniformly at random,
sends the job to the server that generates the token, and removes the token;
when there are no type $i$ tokens, the dispatcher sends the job to an idle
server with an index greater than $\overline{N}$.

Under this conversion, although each server
$\ell\in\\{1,2,\dots,\overline{N}\\}$ in the infinite-server system cannot
obtain a requested job immediately as in the single-server system, it should
be able to receive the job with a diminishing delay since the arrival rates of
jobs scale up to infinity. Therefore, the dynamics of server $\ell$ can be
approximated by that of a single-server system. We formally prove this
approximation through Stein’s method (Braverman et al. 2017, Braverman and Dai
2017, Braverman 2022).

##### Optimal Single-Server Policy.

Now the last piece of the proposed conversion framework is to develop a
single-server policy $\overline{\sigma}$ that achieves the optimal value
$\overline{N}^{*}$. This task is typically much simpler than finding an
optimal policy in the infinite-server system. We design a single-server
policy, Single-OPT, through a linear program based on state-action
frequencies. The optimality is presented as Theorem 4.

Combining the three steps above, we establish the asymptotic optimality of the
meta-policy JRS with the Single-OPT subroutine. We comment that as stated in
the conversion theorem, our policy-conversion framework allows us to convert
not only Single-OPT but rather a general class of single-server policies to
infinite-server policies with performance guarantees. This framework also
avoids the typically tedious work involved in fluid-limit based approaches
such as justifying the validity of fluid limit for each policy.

#### Proof Organization

We present the general problem formulation in Section 2 (in contrast to the
simplified version in Introduction for ease of presentation). In Section 3, we
give a more detailed overview of our main result and approach, with a short
proof of Theorem 1 (main result) based on Theorems 2–4 at the end of the
section. In Section 4, we describe our meta-policy, Join-Requesting-Server
(JRS), in full detail, and prove its performance guarantee (Theorem 3)
assuming an irreducibility assumption. The rest of the proofs, including the
proofs of Theorem 2 (lower bound), Theorem 3 (conversion theorem) without the
irreducibility assumption, and Theorem 4 (single-server analysis) are deferred
to the appendices. We conclude the paper and discuss some future directions in
Section 5.

##  2 Problem Formulation

##### Job Model.

As described in Section 1, we consider a job model where each job in service
can be in one of multiple _phases_ , each phase associated with a different
resource requirement. Here the resource requirement can be a multi-dimensional
vector, with each coordinate specifying the requirement of one type of
resources. To model the temporal variation in the resource requirement, we
assume that each job transitions between phases while in service until it is
completed. The phase transition process is described by a continuous-time
Markov chain on the state space $\mathcal{I}\cup\\{\perp\\}$, where
$\mathcal{I}$ is the set of phases and $\perp$ is the absorbing state that
denotes the completion of the job. We call a transition between two phases in
$\mathcal{I}$ an internal transition, and let $\mu_{ii^{\prime}}$ denote the
transition rate from phase $i$ to phase $i^{\prime}$; the departure of a job
then corresponds to a transition from a phase $i\in\mathcal{I}$ to $\perp$,
whose transition rate is denoted as $\mu_{i\perp}$. The phase transitions of
different jobs are assumed to be independent of each other.

We classify a job as a type $i$ job if it starts from phase $i\in\mathcal{I}$
when entering service. Jobs of each type $i$ arrive to the system according to
an independent Poisson process with rate $\Lambda_{i}$.

##### Server Model.

We consider an infinite-server system with identical servers. As soon as a job
arrives to the system, the job needs to be dispatched to a server to start
service immediately. Note that this is always feasible because there are an
infinite number of servers in the system. We assume that jobs cannot be
preempted or migrated. To describe the state of a server, we define the
configuration of a server as an
$\left\lvert\mathcal{I}\right\rvert$-dimensional vector
${\bm{k}}=(k_{i})_{i\in\mathcal{I}}$, whose $i$-th entry $k_{i}$ is the number
of jobs in phase $i$ on the server. Each server has a limit on the total
number of jobs in service at the same time. This limit is denoted as
$K_{\max}$ and referred to as the service limit. Then the set of feasible
server configurations is
$\mathcal{K}\triangleq\\{{\bm{k}}\colon\sum_{i\in\mathcal{I}}k_{i}\leq
K_{\max}\\}$.

##### System Dynamics.

The system state can be represented by the concatenation of the configurations
of all servers. Specifically, we index the servers by positive integer numbers
and denote the configuration of server $\ell$ at time $t$ as
${\bm{K}}^{\ell}(t)$. Then the state of the entire system can be represented
by the infinite vector $({\bm{K}}^{\ell}(t))_{\ell\in\mathbb{Z}_{+}}$.

Suppose that the system is in state
$({\bm{k}}^{\ell})_{\ell\in\mathbb{Z}_{+}}$. Let $\bm{e}_{i}$ be an
$\left\lvert\mathcal{I}\right\rvert$-dimensional vector whose $i$-th entry is
$1$ and all other entries are $0$. Then the following state transitions can
happen:

  * •

${\bm{k}}^{\ell}\rightarrow{\bm{k}}^{\ell}+\bm{e}_{i}$,
${\bm{k}}^{\ell^{\prime}}\rightarrow{\bm{k}}^{\ell^{\prime}}$
$\forall\ell^{\prime}\neq\ell$: a type $i$ job arrives and is dispatched to
server $\ell$;

  * •

${\bm{k}}^{\ell}\rightarrow{\bm{k}}^{\ell}+\bm{e}_{i^{\prime}}-\bm{e}_{i}$,
${\bm{k}}^{\ell^{\prime}}\rightarrow{\bm{k}}^{\ell^{\prime}}$
$\forall\ell^{\prime}\neq\ell$: a job on server $\ell$ transits from phase $i$
to phase $i^{\prime}$;

  * •

${\bm{k}}^{\ell}\rightarrow{\bm{k}}^{\ell}-\bm{e}_{i}$,
${\bm{k}}^{\ell^{\prime}}\rightarrow{\bm{k}}^{\ell^{\prime}}$
$\forall\ell^{\prime}\neq\ell$: a job on server $\ell$ departs the system from
phase $i$.

The specifics of the system dynamics depend on the employed _dispatch policy_
that decides which server to dispatch to when a job arrives.

##### Active Servers.

We are interested in the number of _active servers_ , i.e., servers currently
serving a positive number of jobs. Note that given the arrival rates of jobs,
the smaller the number of active servers, the better the system is utilized.
Let $X_{{\bm{k}}}(t)$ be the number of servers in configuration ${\bm{k}}$ at
time $t$, i.e.,
$X_{{\bm{k}}}(t)=\sum_{\ell=1}^{\infty}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(t)={\bm{k}}\\}}.$
Then the number of active servers can be written as
$\sum_{{\bm{k}}\neq\bm{0}}X_{{\bm{k}}}(t),$ where
$\bm{0}\in\mathbb{R}^{|\mathcal{I}|}$ is the zero vector.

##### Cost of Resource Contention.

Recall that the cost rate function $h(\cdot)$ maps a server’s configuration to
a rate of cost. We assume that $h(\cdot)$ is any function that is
$\Gamma$-Lipschitz continuous with respect to the $L^{1}$ distance for some
constant $\Gamma>0$ and satisfies $h(\bm{0})=0$.

##### Performance Goal.

Our high-level goal is to design dispatch policies that minimize the number of
active users while keeping the cost rate of resource contention within a
certain budget. Specifically, we consider policies that are allowed to be
randomized and non-Markovian (i.e., the policies can make history-dependent
decisions). We further focus on policies that induce a unique stationary
distribution on the configuration process
$\\{({\bm{K}}^{\ell}(t))_{\ell\in\mathbb{Z}_{+}}\\}$, assuming that the
configuration process is embedded in a Markov chain that has a stationary
distribution. We are interested in such policies because the resulting time
averages of quantities related to the configurations are equal to the
corresponding expectations under the unique stationary distribution regardless
of the initial state. Let $\sigma$ be a policy of interest,
$({\bm{K}}^{\ell})_{\ell\in\mathbb{Z}_{+}}$ be a random element that follows
the stationary distribution of the system state induced by $\sigma$, and
$X_{{\bm{k}}}$ be the corresponding number of servers in configuration
${\bm{k}}$ in steady state under $\sigma$. Then the expected number of active
servers is given by

$N(\sigma)\triangleq\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}],$

and the average expected cost rate per active server is given by

$C(\sigma)\triangleq\frac{\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}]}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}]}.$

Now our goal can be formulated as the following optimization problem, referred
to as problem $\mathcal{P}((\Lambda_{i})_{i\in\mathcal{I}},\epsilon)$:

$\displaystyle\underset{\sigma}{\text{minimize}}$ $\displaystyle N(\sigma)$
(2) subject to $\displaystyle C(\sigma)\leq\epsilon,$

where $\epsilon$ is a budget for the cost of resource contention.

##### Asymptotic Optimality.

We focus on the asymptotic regime where for all $i\in\mathcal{I}$, the arrival
rate $\Lambda_{i}=\lambda_{i}r$ for some constant coefficient $\lambda_{i}$
and a positive scaling factor $r\to+\infty$. To define asymptotic optimality,
we first define the following notion of approximation to the optimization
problem $\mathcal{P}((\Lambda_{i})_{i\in\mathcal{I}},\epsilon)$ in (2): a
policy $\sigma$ is said to be $(\alpha,\beta)$-optimal if
$N(\sigma)\leq\alpha\cdot N^{*}((\Lambda_{i})_{i\in\mathcal{I}},\epsilon)$ and
$C(\sigma)\leq\beta\cdot\epsilon$, where
$N^{*}((\Lambda_{i})_{i\in\mathcal{I}},\epsilon)$ is the optimal objective
value in (2). Now consider a family of policies $\sigma^{(r)}$ indexed by the
scaling factor $r$. We say that the policy $\sigma^{(r)}$ is _asymptotically
optimal_ if it is $\left(\alpha^{(r)},\beta^{(r)}\right)$-optimal to the
optimization problem $\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$
with $\alpha^{(r)},\beta^{(r)}\to 1$ as $r\to\infty$. We will suppress the
superscript $(r)$ for simplicity when it is clear from the context.

We note that under any policy $\sigma$, $N(\sigma)=\Theta(r)$. This can be
proven using the renowned Little’s Law (Kleinrock 1975) in the following way.
The total job arrival rate is $\Theta(r)$ and the expected time that a job
spends in the system is $O(1)$. So by Little’s Law, the expected number of
jobs in the system in steady state is $\Theta(r)$. Since each server can
accommodate a constant number of jobs, the expected number of active servers
is $\Theta(r)$. Given this, $\left(\alpha^{(r)},\beta^{(r)}\right)$-optimality
implies an optimality gap of $\alpha^{(r)}\cdot r$ in the objective value.

##  3 Main Result and Our Approach

###  3.1 Main Result

Our main result, Theorem 1, is the asymptotic optimality of our proposed
policy Join-Requesting-Server (JRS), with a subroutine we call Single-OPT, as
briefly discussed in Section 1. This asymptotic optimality result implies an
$O(\sqrt{r})$ optimality gap in the expected number of active servers. We
defer the detailed descriptions of JRS and Single-OPT to Section 4 and
Appendix C. Theorem 1 follows immediately from Theorems 2–4 to be introduced
in Section 3.2; a short proof is included at the end of Section 3.2 for
clarity.

######  Theorem 1 (Asymptotic Optimality).

Consider a stochastic bin-packing problem in service systems with time-varying
job resource requirements. Let the arrival rates be
$(\lambda_{i}r)_{i\in\mathcal{I}}$ and the cost rate budget be $\epsilon$.
Then the policy Join-Requesting-Server (JRS) with the subroutine Single-OPT is
$\left(1+O\left(r^{-0.5}\right),1+O\left(r^{-0.5}\right)\right)$-optimal. That
is, the expected number of active servers under JRS with Single-OPT is at most
$\left(1+O\left(r^{-0.5}\right)\right)$ times the optimal value of the problem
$\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, while the cost rate
incurred is at most $\left(1+O\left(r^{-0.5}\right)\right)\cdot\epsilon$.

##### Specialization to Non-Time-Varying Resource Requirements.

As mentioned in Section 1, we can specialize this result to the setting where
the resource requirement of a job does not vary over time. To do that, we
remove the cost constraint in
$\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, and redefine the set
of feasible server configurations, $\mathcal{K}$, to also incorporate hard
capacity constraints for each type of resources. The rest of the analysis is
almost identical to that of the analysis for time-varying resource
requirements; we will omit the details due to the space limit. This
specialization will result in a policy that is
$\left(1+O\left(r^{-0.5}\right)\right)$-optimal in the expected number of
active servers.

###  3.2 Our Approach

In a nutshell, our approach is to reduce the original optimization problem in
an infinite-server system to an optimization problem in a single-server
system, which is defined below.

##### A Single-Server System.

Consider a single-server system serving jobs with time-varying resource
requirements. The system has an infinite supply of jobs of all types. As a
result, the server can request any number of new jobs of any type at any time.
Once a job is requested, it immediately enters service.

We represent the server configuration at time $t$ using a vector
$\overline{{\bm{K}}}(t)=(\overline{K}_{i}(t))_{i\in\mathcal{I}}$, whose $i$-th
entry denotes the number of jobs in phase $i$. We assume that the single-
server system has the same service limit $K_{\max}$ and cost rate function
$h(\cdot)$ as a server in the original infinite-server system. Therefore, the
server configuration $\overline{{\bm{K}}}(t)$ is also in the set
$\mathcal{K}=\\{{\bm{k}}\colon\sum_{i\in\mathcal{I}}k_{i}\leq K_{\max}\\}$,
and the cost rate at time $t$ is $h(\overline{{\bm{K}}}(t))$.

A single-server policy $\overline{\sigma}$ determines when and how many jobs
of each type to request. We allow the policy to be randomized and make
decisions based on the current configuration and the history. Let
$(\pi({\bm{k}}))_{{\bm{k}}\in\mathcal{K}}$ be a stationary distribution of the
server configuration under the policy $\overline{\sigma}$, and let
$\overline{{\bm{K}}}(\infty)$ be a random variable following the stationary
distribution. When we consider a policy $\overline{\sigma}$ and its stationary
distribution $\pi$, we assume that the system is initialized from $\pi$. The
policy $\overline{\sigma}$ together with $\pi$ defines the request rate of
type $i$ jobs $\overline{\lambda}_{i}$, which is the expected number of type
$i$ jobs requested per unit time in steady state. Note that
$\overline{\lambda}_{i}$ is the throughput of type $i$ jobs since the system
has a finite state space.

We consider the following single-server optimization problem, denoted as
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$:

$\displaystyle\underset{\overline{N},\mspace{10.0mu}\overline{\sigma},\mspace{10.0mu}\pi}{\text{minimize}}$
$\displaystyle\overline{N}$ (3) subject to
$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]\leq\epsilon,$
$\displaystyle\overline{N}\cdot\overline{\lambda}_{i}=\lambda_{i}r,\quad\forall
i\in\mathcal{I}.$

The single-server optimization problem can be interpreted as follows. We can
think of $\overline{N}$ as the number of copies of the single-server system
under $\overline{\sigma}$ needed to support the arrival rates
$(\lambda_{i}r)_{i\in\mathcal{I}}$ in the infinite-server system. To minimize
$\overline{N}$, it is equivalent to maximizing the throughput
$(\overline{\lambda}_{i})_{i\in\mathcal{I}}$ in each single-server system,
while maintaining their proportions.

We remark that for the problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, we only
need to consider policies that do not depend on the scaling factor $r$. To see
this, we can replace the decision variable $\overline{N}$ with
$\overline{n}\triangleq\overline{N}/r$ and the optimization problem can be
equivalently formulated as follows, which does not involve $r$:

$\displaystyle\underset{\overline{n},\mspace{10.0mu}\overline{\sigma},\mspace{10.0mu}\pi}{\text{minimize}}$
$\displaystyle\overline{n}$ (4) subject to
$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]\leq\epsilon,$
$\displaystyle\overline{n}\cdot\overline{\lambda}_{i}=\lambda_{i},\quad\forall
i\in\mathcal{I}.$

##### Lower Bound.

The single-server optimization problem gives a lower bound to the original
problem in (2) as stated in the following theorem. The proof is given in
Appendix A.

######  Theorem 2 (Lower Bound).

Consider a stochastic bin-packing problem in service systems with time-varying
job resource requirements. Let the arrival rates be
$(\lambda_{i}r)_{i\in\mathcal{I}}$ and the cost rate budget be $\epsilon$. Let
$N^{*}$ be the optimal value of the original infinite-server problem in (2),
and let $\overline{N}^{*}$ be the optimal value of the single-server
optimization problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, then
$N^{*}\geq\overline{N}^{*}$.

##### Converting From the Single-Server System to the Infinite-Server System.

Having established a lower bound on the infinite-server problem
$\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in terms of the
optimal value of the single-server problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, next we
focus on finding an asymptotically optimal policy. In particular, we will
characterize the performance guarantee of a class policies and then show that
the best policy within the class is asymptotically optimal. Specifically, we
consider a meta-policy called Join-Requesting-Server (JRS), which converts a
Markovian single-server policy $\overline{\sigma}$ into an infinite-server
policy. We call the policy resulting from the conversion a JRS policy with a
subroutine $\overline{\sigma}$. Through analyzing the meta-policy JRS, we show
that the performance of each JRS policy can be characterized by the
performance of its subroutine, as stated in Theorem 3 below. The proof of
Theorem 3 under an irreducibility assumption is given in Section 4, and the
proof for the full version is given in Appendix B.4.

######  Theorem 3 (Conversion Theorem).

Consider a stochastic bin-packing problem in service systems with time-varying
job resource requirements. Let the arrival rates be
$(\lambda_{i}r)_{i\in\mathcal{I}}$ and the cost rate budget be $\epsilon$. Let
$(\overline{N},\overline{\sigma},\pi({\bm{k}}))$ be a solution feasible to the
single-server optimization problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$. In
addition, we assume that the policy $\overline{\sigma}$ is Markovian. Let the
infinite-server policy $\sigma$ be JRS with a subroutine $\overline{\sigma}$.
Then under $\sigma$, we have

$\displaystyle\left\lvert\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}\left[X^{(r)}_{{\bm{k}}}\right]-\left\lceil\overline{N}\right\rceil\cdot\mathbb{P}\left(\overline{{\bm{K}}}\neq\bm{0}\right)\right\rvert$
$\displaystyle=O\left(r^{0.5}\right),$ (5)
$\displaystyle\left\lvert\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}\left[X^{(r)}_{{\bm{k}}}\right]-\left\lceil\overline{N}\right\rceil\cdot\mathbb{E}\left[h(\overline{{\bm{K}}})\right]\right\rvert$
$\displaystyle=O\left(r^{0.5}\right).$ (6)

As a result,

$\displaystyle N(\sigma)$
$\displaystyle\leq\left(1+O\left(r^{-0.5}\right)\right)\cdot\overline{N},$ (7)
$\displaystyle C(\sigma)$
$\displaystyle\leq\left(1+O\left(r^{-0.5}\right)\right)\cdot\epsilon.$ (8)

We briefly describe JRS here. Let each server with index
$\ell\leq\lceil\overline{N}\rceil$ request new jobs based on its
configuration, according to the input single-server policy
$\overline{\sigma}$. A request for a type $i$ job is sent to the centralized
dispatcher in the form of a type $i$ token. When a type $i$ job arrives, the
dispatcher randomly chooses a type $i$ token, removes the token, and sends the
job to the server that generates the token; when there are no type $i$ tokens,
the dispatcher sends the job to a server with index
$\ell>\lceil\overline{N}\rceil$. The full description of JRS is presented in
Section 4, which includes additional details that we omit here.

Under JRS with a subroutine $\overline{\sigma}$, the configurations of the
first $\lceil\overline{N}\rceil$ servers
$({\bm{K}}^{1}(t),\dots,{\bm{K}}^{\lceil\overline{N}\rceil}(t))$ have
approximately the same steady-state distribution as $\lceil\overline{N}\rceil$
independent and identical copies of the single-server system under
$\overline{\sigma}$. Moreover, the number of active servers with index greater
than $\lceil\overline{N}\rceil$ is diminishing with $r$. These two facts allow
us to connect the performance metrics of the infinite-server problem (2) and
the single-server problem (3). The intuition behind these facts is that when
the arrival rate goes to infinity, each server should be able to receive the
requested jobs with diminishing delays, making them behave like independent
copies of the single-server system.

##### Optimal Single-Server Policy.

Theorem 3 together with the lower bound in Theorem 2 reduces the infinite-
server optimization problem
$\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (2) to the single-
server optimization problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (3). To
solve the single-server optimization problem, we consider the linear program
in (9) below, denoted as
$\overline{\mathcal{LP}}((\lambda_{i})_{i\in\mathcal{I}},\epsilon)$ with
decision variables $\Phi\in\mathbb{R}$, $\bm{\pi}\in\mathbb{R}^{\mathcal{K}}$,
and $\bm{u}_{i}\in\mathbb{R}^{\mathcal{K}}$ for $i\in\mathcal{I}$. This linear
program is derived as a relaxation to the optimization problem (3), and
$\bm{\pi}\in\mathbb{R}^{\mathcal{K}}$ and
$\bm{u}_{i}\in\mathbb{R}^{\mathcal{K}}$ in the solution will be used to
construct a single-server policy:

$\displaystyle\underset{\begin{subarray}{c}\Phi,\mspace{10.0mu}\bm{\pi},\mspace{8.0mu}(\bm{u}_{i})_{i\in\mathcal{I}}\end{subarray}}{\text{maximize}}$
$\displaystyle\Phi$ (9) subject to
$\displaystyle\bm{h}^{T}\bm{\pi}\leq\epsilon(1-\pi(\bm{0}))$
$\displaystyle\bm{1}_{o}^{T}\bm{u}_{i}=\Phi\cdot\lambda_{i}\quad\forall i$
$\displaystyle
A\bm{\pi}+\sum\nolimits_{i\in\mathcal{I}}B_{i}\bm{u}_{i}=\bm{0}$
$\displaystyle\bm{1}^{T}\bm{\pi}=1$ $\displaystyle\bm{\pi}\geq
0,\bm{u}_{i}\geq 0\quad\forall i\in\mathcal{I}$

where $\bm{h}$ is the vector form of the cost rate function $h$;
$\bm{1}_{o}^{T}$ is a $\left\lvert\mathcal{K}\right\rvert$-dimensional vector
with one in all entries except those with
$\sum_{i\in\mathcal{I}}k_{i}=K_{\max}$. The linear equation
$A\bm{\pi}+\sum_{i\in\mathcal{I}}B_{i}\bm{u}_{i}=0$ is explicitly given by

$\displaystyle\mspace{22.0mu}\sum_{i}u_{i}({\bm{k}}-\bm{e}_{i})\mathbbm{1}_{\\{k_{i}\geq
1\\}}+\sum_{i}(k_{i}+1)\mu_{i\perp}\pi({\bm{k}}+\bm{e}_{i})$
$\displaystyle\mspace{22.0mu}+\sum_{i,i^{\prime}\colon i\neq
i^{\prime}}(k_{i}+1)\mu_{ii^{\prime}}\pi({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}})\mathbbm{1}_{\\{k_{i^{\prime}}\geq
1\\}}$
$\displaystyle=\sum_{i}u_{i}({\bm{k}})+\Big{(}\sum_{i}k_{i}\mu_{i\perp}+\sum_{i,i^{\prime}\colon
i\neq
i^{\prime}}k_{i}\mu_{ii^{\prime}}\Big{)}\pi({\bm{k}}),\quad\forall{\bm{k}}\in\mathcal{K}.$

The decision variables and constraints of the problem have the following
interpretations: $\Phi$ is a factor that scales inversely with the request
rates of the system; $\bm{\pi}$ is a column vector representing the stationary
distribution of configuration; $\bm{u}_{i}$ will be referred to as transition
frequency, whose ${\bm{k}}$-th entry $u_{i}({\bm{k}})$ represents the number
of transitions from configuration ${\bm{k}}$ to configuration
${\bm{k}}+\bm{e}_{i}$ per unit time in steady state, indicating how often type
$i$ jobs should be requested when the system is in configuration ${\bm{k}}$.
The first two constraints of the linear program correspond to the first two
constraints in
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, and the
last three constraints characterize the fact that $\bm{\pi}$ is a stationary
distribution.

Given a solution $(\Phi,\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ feasible to
the linear program in (9), we construct a single-server policy using timers
whose durations follow exponential distributions. In particular, if the system
enters a configuration ${\bm{k}}$ with $\pi({\bm{k}})\neq 0$, the policy
starts a timer with an exponentially distributed duration with rate
$\frac{u_{i}({\bm{k}})}{\pi({\bm{k}})}$ for each $i\in\mathcal{I}$, and
requests a type $i$ job when the corresponding timer ticks. Note that the rate
of the timer gives us the desired transition frequency in steady-state because
$\frac{u_{i}({\bm{k}})}{\pi({\bm{k}})}\cdot\pi({\bm{k}})=u_{i}({\bm{k}})$. The
case where $\pi({\bm{k}})=0$ is dealt with separately. We refer to this
single-server policy as
$\overline{\sigma}(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$. We call such
single-server policies LP-based policies, and we call an LP-based policy
obtained from an optimal LP solution Single-OPT. Theorem 4 below states that
Single-OPT is an optimal policy for
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$. More
details on the construction of the policy and all the analysis of the single-
server system are given in Appendix C.

######  Theorem 4 (Optimality of Single-OPT).

Given an optimal solution
$(\Phi^{*},\bm{\pi}^{*},(\bm{u}_{i})_{i\in\mathcal{I}})$ to the linear program
$\overline{\mathcal{LP}}((\lambda_{i})_{i\in\mathcal{I}},\epsilon)$, we can
solve the single-server optimization problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (3) to
achieve an optimal value $r/\Phi^{*}$, using the optimal policy
$\overline{\sigma}(\bm{\pi}^{*},(\bm{u}_{i})_{i\in\mathcal{I}}^{*})$ and the
optimal stationary distribution $\bm{\pi}^{*}$. Moreover, the optimal policy
$\overline{\sigma}(\bm{\pi}^{*},(\bm{u}_{i})_{i\in\mathcal{I}}^{*})$ is a
Markovian policy.

###### Proof of the Theorem 1 Based on Theorems 2–4..

Because Single-OPT (along with an optimal stationary distribution) optimally
solves $\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, it
achieves the optimal value $\overline{N}^{*}$. Let $\sigma$ be JRS with a
subroutine Single-OPT, then according to Theorem 3, we have
$N(\sigma)\leq\left(1+O\left(r^{-0.5}\right)\right)\cdot\overline{N}^{*}$ and
$C(\sigma)\leq\left(1+O\left(r^{-0.5}\right)\right)\cdot\epsilon$. By Theorem
2, we also have $N^{*}\geq\overline{N}^{*}.$ So we conclude that JRS with a
subroutine Single-OPT is $\left(1+O(r^{-0.5}),1+O(r^{-0.5})\right)$-optimal. ∎

##  4 Join-Requesting-Server and Performance Analysis

In this section, we first describe our meta-policy, Join-Requesting-Server
(JRS), in full detail. Then we prove Theorem 3 to establish the performance
guarantee for JRS. For ease of presentation, we focus on the case where the
subroutine policy $\overline{\sigma}$ for JRS is ${\bm{k}}^{0}$-irreducible,
i.e., under $\overline{\sigma}$, there exists a configuration ${\bm{k}}^{0}$
such that the single-server system can return to ${\bm{k}}^{0}$ from any other
configurations (which is equivalent to assuming that the configuration of the
single-server system under policy $\overline{\sigma}$ forms a _unichain_). The
algorithm and analysis for the general case are in Appendix B.4. The proofs of
all the lemmas presented in this section are given in Appendix B.

###  4.1 How the Single-Server Policy Request Jobs

Before going into the definition of JRS, we first describe how the Markovian
single-server policy $\overline{\sigma}$ requests jobs. Note that the
description of a single-server policy in this section, referred to as a _non-
impulse view_ , is different from the description specified by the solution to
the LP in (9) in Section 3.2, which is referred to as an _impulse view_ (see
Remark 2 in Appendix C.2 for a more detailed discussion). But one can verify
that for any policy specified from the impulse view, we can describe it using
the non-impulsive view in this section. Let $a_{i}$ denote the number of type
$i$ jobs requested, and let $\bm{a}\triangleq(a_{i})_{i\in\mathcal{I}}$. We
say $\bm{a}$ is feasible if the total number of jobs on the server does not
exceed $K_{\max}$ after adding the jobs. The policy $\overline{\sigma}$
performs one of the following two types of requests based on the current
configuration.

  * •

Reactive requests. A _reactive request_ is triggered by _either an internal
transition or a departure_. The changes of the configuration when a reactive
request is made can be represented by the diagram

${\bm{k}}\to{\bm{k}}^{\prime}\to{\bm{k}}^{\prime}+\bm{a},$

where ${\bm{k}}\to{\bm{k}}^{\prime}$ is due to the internal transition or
departure, and ${\bm{k}}^{\prime}\to{\bm{k}}^{\prime}+\bm{a}$ happens since
the policy immediately requests $\bm{a}$ jobs. The policy $\overline{\sigma}$
specifies a probability distribution over all feasible $\bm{a}$ when it
decides to perform reactive requests for the configuration
${\bm{k}}^{\prime}$.

  * •

Proactive requests. A _proactive request_ happens on its own, and it happens
at a _finite rate_ depending on the current configuration of the server. The
change of the configuration when a proactive request happens can be
represented by the diagram

${\bm{k}}\to{\bm{k}}+\bm{a}.$

More specifically, suppose the policy $\overline{\sigma}$ decides to perform
proactive requests for a configuration ${\bm{k}}$. Then for each feasible
$\bm{a}$, the policy $\overline{\sigma}$ specifies a rate and runs a timer
with an exponentially distributed duration with the specified rate. When the
timer ticks, $\bm{a}$ jobs are requested. When the configuration changes, all
the timers are canceled and restarted with new rates based on the new
configuration.

###  4.2 Definition of Join-Requesting-Server

##### Policy description.

The input of JRS are: (1) the single-server policy $\overline{\sigma}$, (2)
the objective value of $\overline{\sigma}$ in the single-server problem (3),
denoted as $\overline{N}$, and (3) the transition rates of the job model.

We first divide the infinite server pool into two sets based on the server
index $\ell$. Let $L=\lceil\overline{N}\rceil$. We call servers with index
$\ell\leq L$ normal servers; we call servers with index $\ell>L$ backup
servers. The normal servers are responsible for serving most of the jobs,
while the backup servers are active only when there are too many jobs on the
normal servers.

The Join-Requesting-Server is specified in two steps.

  * •

Step 1 (Job Requesting on a Normal Server): We let each normal server request
jobs using its subroutine, the single-server policy $\overline{\sigma}$. The
input to the policy $\overline{\sigma}$ is what we refer to as the _observed
configuration_ of the server, which will be further explained later. When
$\overline{\sigma}$ requests $\bm{a}=(a_{i})_{i\in\mathcal{I}}$ jobs, $a_{i}$
type $i$ tokens are generated for each $i\in\mathcal{I}$ to store the job
requests. The server _pauses_ this job requesting process if it already has
tokens on it, and resumes when all the tokens are removed.

  * •

Step 2 (Arrival Assignment): When a type $i$ job arrives, the dispatcher
chooses a type $i$ token uniformly at random, removes the token, and assigns
the job to the corresponding server. There are two exceptions:

    * –

When there are no type $i$ tokens, the dispatcher sends the job to an idle
server with an index greater than $\overline{N}$.

    * –

When the total number of type $i$ tokens throughout the system exceeds the
limit $\eta_{\max}=\left\lceil L^{0.5}\right\rceil$ (we call it the token
limit), a type $i$ virtual arrival is triggered, which causes the dispatcher
to choose a type $i$ token uniformly at random, remove the token, and assign a
virtual job to the corresponding server. A virtual job has the same transition
dynamics as a real job, but does not use resources.

The _observed configuration_ of a normal server, which is the input to the
policy $\overline{\sigma}$, is a vector whose $i$-th entry represents the
total number of real and virtual jobs in phase $i$ on this server. The
observed configuration changes when there is a new job or a new virtual job
arrival assigned to the server, or when a job or a virtual job on the server
has a phase transition or departs. We update the input to the policy
$\overline{\sigma}$ when the observed configuration changes. Whenever the
observed configuration changes, the policy $\overline{\sigma}$ cancels the
exponential timers in progress; but a reactive request from the policy
$\overline{\sigma}$ can only be triggered when a job or a virtual job on the
server has a phase transition or departs.

##### Intuition for why JRS works.

Our intuition for why each normal server’s configuration in the infinite-
server system is close to the configuration in a single-server system is based
on the following observation. In a single-server system, a job request gets
fulfilled immediately; meanwhile, in an infinite-server system, we can show
that after a job request is generated, the request gets fulfilled quickly,
either by a real job or a virtual job.

Because each job request is stored as a token, showing that each job request
gets fulfilled quickly is equivalent to showing that each token gets removed
quickly. To see this, observe that tokens are generated at a rate proportional
to the number of normal servers, $L$, and there are at most $\eta_{\max}$
tokens. By Little’s Law, for each $i\in\mathcal{I}$, a type $i$ token exists
for $O(\eta_{\max}/L)$ expected time. When the arrival rate is of order
$\Theta\left(r\right)$, we have that $L=\Theta\left(r\right)$ and
$\eta_{\max}=\left\lceil L^{0.5}\right\rceil=\Theta\left(r^{0.5}\right)$.
Therefore, a type $i$ token only exists for $O(r^{-0.5})$ expected time.

###  4.3 Proof of Theorem 3 Assuming Irreducibility

##### Preliminaries.

Consider an infinite-server system under the JRS policy. For each normal
server $\ell$, we describe its status at time $t$ using the following
variables: configuration of real jobs ${\bm{K}}^{\ell}(t)$ (referred to as
“configuration” in previous sections), tokens $\bm{\eta}^{\ell}(t)$,
configuration of virtual jobs $\bm{\zeta}^{\ell}(t)$, and observed
configuration
$\widehat{{\bm{K}}}^{\ell}(t)\triangleq{\bm{K}}^{\ell}(t)+\bm{\zeta}^{\ell}(t)$.
We use the superscript “${1:L}$” to refer to a certain descriptor of all
normal servers, for example,
$\widehat{{\bm{K}}}^{1:L}(t)\triangleq\big{(}\widehat{{\bm{K}}}^{\ell}(t)\big{)}_{\ell=1,2,\dots
L}$. The system under JRS is a Markov chain with a unique Markovian
representation
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$.
The following lemma shows that the system has a unique stationary distribution
(the proof is provided in Appendix B.1).

######  Lemma 1 (Unique Stationary Distribution).

Consider an infinite-server system under the JRS policy with
$\overline{\sigma}$ as its subroutine, where $\overline{\sigma}$ is a single-
server policy that is Markovian and ${\bm{k}}^{0}$-irreducible. Then the state
of the system
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$
has a unique stationary distribution.

Let
$\overline{{\bm{K}}}^{1:L}(t)\triangleq\big{(}\overline{{\bm{K}}}^{\ell}(t)\big{)}_{\ell=1,2,\dots,L}$
be the configuration vector of $L$ i.i.d. copies of the single-server system
under $\overline{\sigma}$. As discussed in Section 4.2, we will show that
${\bm{K}}^{1:L}(\infty)$ can be approximated by
$\overline{{\bm{K}}}^{1:L}(\infty)$. In the remainder of this section, we omit
the steady-state symbol $(\infty)$ for simplicity.

To rigorously discuss the approximation of the steady-state random variables,
we define some metrics. Recall that
$\mathcal{K}\triangleq\\{{\bm{k}}\colon\sum_{i\in\mathcal{I}}k_{i}\leq
K_{\max}\\}$ is the set of feasible single-server configurations. Let
$\mathcal{K}^{L}\triangleq\\{{\bm{k}}^{1:L}\colon{\bm{k}}^{\ell}\in\mathcal{K},\forall\ell\\}$
be the set of feasible configurations for all normal servers. We use
$\lVert\cdot\rVert$ to denote the $L^{1}$ norm in both space $\mathcal{K}$ and
space $\mathcal{K}^{L}$:

$\displaystyle\lVert{\bm{k}}-{\bm{k}}^{\prime}\rVert$
$\displaystyle=\textstyle\sum_{i\in\mathcal{I}}\left\lvert
k_{i}-k_{i}^{\prime}\right\rvert,\quad\text{for
}{\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K},$
$\displaystyle\lVert{\bm{k}}^{1:L}-{\bm{k}}^{\prime{1:L}}\rVert$
$\displaystyle=\textstyle\sum_{\ell=1}^{L}\lVert{\bm{k}}^{\ell}-{\bm{k}}^{\prime\ell}\rVert,\quad\text{for
}{\bm{k}}^{1:L},{\bm{k}}^{\prime{1:L}}\in\mathcal{K}^{L}.$

For any two random variables $\bm{U}^{a},\bm{U}^{b}\in\mathcal{K}^{L}$, their
closeness will be measured in terms of Wasserstein distance as follows:

$d(\bm{U}^{a},\bm{U}^{b})\triangleq\sup_{f\in\text{Lip}(1)}\left\\{\mathbb{E}\left[f\left(\bm{U}^{a}\right)\right]-\mathbb{E}[f(\bm{U}^{b})]\right\\},$

where the supremum is taken over the all Lipschitz-$1$ functions from
$\mathcal{K}^{L}$ to $\mathbb{R}$.

##### Proof Outline.

Observe that the process $(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t))$
itself forms a Markov chain and governs the dynamics of other state variables,
including the number of virtual jobs $\bm{\zeta}^{1:L}(t)$ and the number of
jobs on the backup servers. We can thus decompose our analysis into two steps.
In Step 1, we focus on the dynamics of
$(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t))$ and prove that
$\widehat{{\bm{K}}}^{1:L}$ is close to $\overline{{\bm{K}}}^{1:L}$ in terms of
the Wasserstein distance. In Step 2, we use the result of the first step to
show that the total number of virtual jobs
$\sum_{i\in\mathcal{I}}\sum_{\ell=1}^{L}\zeta_{i}^{\ell}$ and the total number
of jobs on the backup servers are both small. Recall that
${\bm{K}}^{1:L}=\widehat{{\bm{K}}}^{1:L}-\bm{\zeta}^{1:L}$. Step 1 and Step 2
together imply that ${\bm{K}}^{1:L}$ and $\overline{{\bm{K}}}^{1:L}$ are close
in the Wasserstein distance, and the backup servers are almost all idle. Since
$\overline{{\bm{K}}}^{1:L}$ are i.i.d. copies of the single-server system, we
thus establish a connection between the real-job configurations of the
infinite-server system $({\bm{K}}^{\ell})_{\ell=1,2,\dots}$ and the stationary
distribution of the single-server system. Then the theorem follows.

###### Step 1.

Lemma 2 below bounds the Wasserstein distance between
$\widehat{{\bm{K}}}^{1:L}$ and $\overline{{\bm{K}}}^{1:L}$.

######  Lemma 2.

Under the conditions of Theorem 3 and $\overline{\sigma}$ being
${\bm{k}}^{0}$-irreducible, we have

$d\left(\overline{{\bm{K}}}^{1:L},\widehat{{\bm{K}}}^{1:L}\right)=O\left(r^{0.5}\right).$
(10)

The proof of Lemma 2 follows the framework of Stein’s method (see, for example
Braverman and Dai 2017, Braverman et al. 2017, Braverman 2022). Stein’s method
usually consists of three steps: generator comparison, Stein factor bounds,
and moment bounds. In our case, due to the finiteness of the state space
$\mathcal{K}$, we only need to do the generator comparison and the Stein
factor bounds. The generator comparison step justifies the intuition that most
normal servers behave like single-server systems under the policy
$\overline{\sigma}$, and the Stein’s factor bounds corresponds to the
intuition that the small fraction of servers whose dynamics deviate from
single-server systems do not contribute much to the overall distance. We defer
the detailed proof to Appendix B.2.

###### Step 2.

We establish Lemma 3 and Lemma 4 below, which bound the steady-state expected
number of virtual jobs and the number of jobs on backup servers, respectively.
The proofs of them are provided in Appendix B.3.

######  Lemma 3.

Under the conditions of Theorem 3 and $\overline{\sigma}$ being
${\bm{k}}^{0}$-irreducible, for each $i\in\mathcal{I}$, the steady-state
expected number of virtual jobs of type $i$ is of the order
$O\left(r^{0.5}\right)$, i.e.,

$\mathbb{E}\left[\textstyle\sum_{\ell=1}^{L}\zeta_{i}^{\ell}\right]=O\left(r^{0.5}\right).$

######  Lemma 4.

Under the conditions of Theorem 3 and $\overline{\sigma}$ being
${\bm{k}}^{0}$-irreducible, for each $i\in\mathcal{I}$, the steady-state
expected number of type $i$ jobs on backup servers is of the order
$O\left(r^{0.5}\right)$, i.e.,

$\mathbb{E}\left[\textstyle\sum_{\ell=L+1}^{\infty}K_{i}^{\ell}\right]=O\left(r^{0.5}\right).$

###### Proof of Theorem 3 with ${\bm{k}}^{0}$-irreducible $\overline{\sigma}$
Based on Lemmas 1–4..

First we show that Lemmas 2 and 3 imply the closeness between ${\bm{K}}^{1:L}$
and $\overline{{\bm{K}}}^{1:L}$. From the Wasserstein distance bound in Lemma
2, we have that for any $f:\mathcal{K}^{L}\to\mathbb{R}$ that is 1-Lipschitz
continuous,

$\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}\big{)}\right]=O\left(r^{0.5}\right).$

Recall that $\widehat{{\bm{K}}}^{\ell}={\bm{K}}^{\ell}+\bm{\zeta}^{\ell}$, and
$\mathbb{E}\left[\sum_{\ell=1}^{L}\zeta_{i}^{\ell}\right]=O\left(r^{0.5}\right)$
by Lemma 3. Then by the Lipschitz continuity of $f$, we have

$\mathbb{E}\left[f(\overline{{\bm{K}}}^{1:L})\right]-\mathbb{E}\left[f({\bm{K}}^{1:L})\right]=O\left(r^{0.5}\right).$
(11)

Now we prove $\eqref{eq:sat:active-servers}$, the bound on expected number of
the active servers, by leveraging the Wasserstein distance bound above. We
write the difference in the expected number of active servers between the two
systems as follows:

$\displaystyle\mspace{23.0mu}\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}\Big{[}X_{{\bm{k}}}\Big{]}-L\cdot\mathbb{P}(\overline{{\bm{K}}}\neq
0)$
$\displaystyle=\mathbb{E}\Big{[}\sum_{\ell=1}^{L}\mathbbm{1}_{\\{{\bm{K}}^{\ell}\neq\bm{0}\\}}\Big{]}+\mathbb{E}\Big{[}\sum_{\ell=L+1}^{\infty}\mathbbm{1}_{\\{{\bm{K}}^{\ell}\neq\bm{0}\\}}\Big{]}-\mathbb{E}\Big{[}\sum_{\ell=1}^{L}\mathbbm{1}_{\\{\overline{{\bm{K}}}^{\ell}\neq\bm{0}\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}f_{1}({\bm{K}}^{1:L})\Big{]}-\mathbb{E}\Big{[}f_{1}(\overline{{\bm{K}}}^{1:L})\Big{]}+\mathbb{E}\Big{[}\sum_{\ell=L+1}^{\infty}\mathbbm{1}_{\\{K_{i}^{\ell}\neq\bm{0}\\}}\Big{]},$
(12)

where $f_{1}$ is given by
$f_{1}({\bm{k}}^{1:L})=\sum_{\ell=1}^{L}\mathbbm{1}_{\\{{\bm{k}}^{\ell}\neq\bm{0}\\}}.$
Note that for any ${\bm{k}}^{1:L},{\bm{k}}^{\prime{1:L}}\in\mathcal{K}^{L}$,
we have

$\displaystyle\left\lvert
f_{1}({\bm{k}}^{1:L})-f_{1}({\bm{k}}^{\prime{1:L}})\right\rvert=\left\lvert\textstyle\sum_{\ell=1}^{L}\left(\mathbbm{1}_{\\{{\bm{k}}^{\ell}\neq\bm{0}\\}}-\mathbbm{1}_{\\{{\bm{k}}^{\prime\ell}\neq\bm{0}\\}}\right)\right\rvert\leq\textstyle\sum_{\ell=1}^{L}\mathbbm{1}_{\\{{\bm{k}}^{\ell}\neq{\bm{k}}^{\prime\ell}\\}}\leq\lVert{\bm{k}}^{1:L}-{\bm{k}}^{\prime{1:L}}\rVert.$

That is, $f_{1}\in\text{Lip(1)}$. Applying (11) and Lemma 4 to (12), we get

$\mathbb{E}\left[\textstyle\sum_{{\bm{k}}\neq\bm{0}}X_{{\bm{k}}}\right]-L\cdot\mathbb{P}(\overline{{\bm{K}}}\neq
0)=O\left(r^{0.5}\right).$

Similarly, to prove (6), we write the difference in cost between the two
systems as follows:

$\displaystyle\mspace{23.0mu}\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}]-L\cdot\mathbb{E}\Big{[}h(\overline{{\bm{K}}})\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\sum_{\ell=1}^{L}h({\bm{K}}^{\ell})\Big{]}+\mathbb{E}\Big{[}\sum_{\ell=L+1}^{\infty}h({\bm{K}}^{\ell})\Big{]}-\mathbb{E}\Big{[}\sum_{\ell=1}^{L}h(\overline{{\bm{K}}}^{\ell})\Big{]}$
$\displaystyle=\Gamma\cdot f_{2}\big{(}{\bm{K}}^{1:L}\big{)}-\Gamma\cdot
f_{2}\big{(}\overline{{\bm{K}}}^{1:L}\big{)}+\mathbb{E}\Big{[}\sum_{\ell=L+1}^{\infty}h({\bm{K}}^{\ell})\Big{]}\cdot\max_{{\bm{k}}\in\mathcal{K}}h({\bm{k}}),$

where $\Gamma$ is the Lipschitz constant of $h(\cdot)$, and $f_{2}$ is given
by
$f_{2}({\bm{k}}^{1:L})=\frac{1}{\Gamma}\sum_{\ell=1}^{L}h({\bm{k}}^{\ell}).$
For any ${\bm{k}}^{1:L},{\bm{k}}^{\prime{1:L}}\in\mathcal{K}^{L}$, we have

$\left\lvert
f_{2}({\bm{k}}^{1:L})-f_{2}({\bm{k}}^{\prime{1:L}})\right\rvert=\left\lvert{\Gamma}^{-1}\textstyle\sum_{\ell=1}^{L}(h({\bm{k}}^{\ell})-h({\bm{k}}^{\prime\ell}))\right\rvert\leq\textstyle\sum_{\ell=1}^{L}\lVert{\bm{k}}^{\ell}-{\bm{k}}^{\prime\ell}\rVert=\lVert{\bm{k}}^{1:L}-{\bm{k}}^{\prime{1:L}}\rVert,$

i.e., $f_{2}\in\text{Lip(1)}$. Applying (11) and Lemma 4 to (12), we get

$\textstyle\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}]-L\cdot\mathbb{E}\left[h(\overline{{\bm{K}}})\right]=O\left(r^{0.5}\right).$

To show (7) and (8), noting that $\left\lceil
N\right\rceil=\Theta\left(r\right)$, we have

$N(\sigma)=\textstyle\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}]\leq\left\lceil\overline{N}\right\rceil+O\left(r^{0.5}\right)=\left(1+O\left(r^{-0.5}\right)\right)\cdot\overline{N},$

and

$\displaystyle C(\sigma)$
$\displaystyle=\frac{\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}]}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]}$
$\displaystyle=\frac{\sum_{j=1}^{J}\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\pi({\bm{k}})+O\left(r^{-0.5}\right)}{1-\pi(\bm{0})+O\left(r^{-0.5}\right)}$
$\displaystyle=\left(1+O\left(r^{-0.5}\right)\right)\cdot\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]$
$\displaystyle\leq\left(1+O\left(r^{-0.5}\right)\right)\cdot\epsilon.$

This completes the proof.

∎

##  5 Conclusion

In this paper, we study a new setting of stochastic bin-packing in service
systems that features time-varying item sizes. Since our formulation is
motivated by the problem of virtual-machine scheduling in computing systems,
we use the terminology of jobs and servers, where jobs correspond to items,
whose sizes are their resource requirements, and servers correspond to bins.
Then the time-varying item sizes capture the widespread phenomenon in practice
that jobs’ resource requirements vary over time. Our goal is to design a job
dispatch policy to minimize the expected number of active servers in steady
state, subject to a constraint on the cost of resource contentions. Our main
result is the design of a policy that achieves an optimality gap of
$O(\sqrt{r})$, where $r$ is the scaling factor of the arrival rate. When
specialized to the traditional setting where jobs’ resource requirements
remain fixed over time, this result improves upon the state-of-the-art $o(r)$
optimality gap. Our technical approach highlights a novel policy conversion
framework, Join-Requesting-Server, that reduces the policy design problem to
that in a single-server system.

There are several potential directions that may be worth further exploration.
One direction is to strengthen the optimality result within the current
setting. Specifically, it is interesting to investigate (i) whether it is
possible to achieve an optimality gap smaller than $\Theta(\sqrt{r})$; and
(ii) whether there exist asymptotically optimal policies whose average cost
rate of resource contention satisfies the budget strictly instead of
asymptotically.

We are also interested in extending our technique to the optimal control of
other systems with similar structures. Intuitively, this technique could be
applied to systems with many components that evolve independently but are
weakly coupled by certain constraints. Viewing each component as a server, we
can define a suitable single-server problem and then design a policy for the
original system to track the dynamics of the optimal single-server solution.
Below we list several variations of our model that can potentially be analyzed
using the proposed technique.

  * •

A model where jobs running on each server will be put into a local queue when
there are resource contentions. The goal thus becomes finding the optimal
trade-offs between the number of active servers and the waiting time of the
jobs.

  * •

A model that allows each server to have a Markovian state that affects the
dynamics of the jobs running on the server.

  * •

A model that allows jobs to migrate to different servers at the cost of
migration delays.

  * •

A closed-system model where jobs re-enter the system after completion.

A third possible direction is to tackle the problem when the arrival rates and
the parameters in the job model are unknown. A possible approach is to develop
an approximate version of the JRS framework, where the optimal single-server
policy and the simulator for the virtual jobs are both learned from data. It
is desirable to design such an approximate framework whose performance
degrades gracefully as the approximation error increases.

## References

  * Ayyadevara et al. (2022) Ayyadevara N, Dabas R, Khan A, Sreenivas KVN (2022) Near-Optimal Algorithms for Stochastic Online Bin Packing. _Proc. Int. Conf. Automata, Languages and Programming (ICALP)_ , volume 229, 12:1–12:20.
  * Bashir et al. (2021) Bashir N, Deng N, Rzadca K, Irwin D, Kodak S, Jnagal R (2021) Take it to the limit: Peak prediction-driven resource overcommitment in datacenters. _Proc. European Conf. Computer Systems (EuroSys)_ , 556–573 (Online Event, United Kingdom).
  * Braverman (2022) Braverman A (2022) The prelimit generator comparison approach of stein’s method. _Stoch. Syst._ 12(2):181–204.
  * Braverman and Dai (2017) Braverman A, Dai JG (2017) Stein’s method for steady-state diffusion approximations of $m/\mathit{Ph}/n+m$ systems. _Ann. Appl. Probab._ 27:550–581.
  * Braverman et al. (2017) Braverman A, Dai JG, Feng J (2017) Stein’s method for steady-state diffusion approximations: an introduction through the Erlang-A and Erlang-C models. _Stoch. Syst._ 6(2):301–366.
  * Buchbinder et al. (2021) Buchbinder N, Fairstein Y, Mellou K, Menache I, Naor JS (2021) Online virtual machine allocation with lifetime and load predictions. _ACM SIGMETRICS Perform. Eval. Rev._ 49(1):9–10.
  * Cloud (2023) Cloud G (2023) Overcommitting CPUs on sole-tenant VMs. https://cloud.google.com/compute/docs/nodes/overcommitting-cpus-sole-tenant-vms.
  * Coffman et al. (1983) Coffman EG Jr, Garey MR, Johnson DS (1983) Dynamic bin packing. _SIAM J. Comput._ 12(2):227–258.
  * Coffman et al. (1984) Coffman EG Jr, Garey MR, Johnson DS (1984) _Approximation Algorithms for Bin-Packing — An Updated Survey_ , 49–106 (Springer Vienna).
  * Courcobetis and Weber (1990) Courcobetis C, Weber R (1990) Stability of on-line bin packing with random arrivals and long-run-average constraints. _Probab. Eng. Inf. Sci._ 4(4):447–460.
  * Courcoubetis and Weber (1986) Courcoubetis C, Weber RR (1986) Necessary and sufficient conditions for stability of a bin-packing system. _J. Appl. Probab._ 23(4):989–999.
  * Csirik et al. (2006) Csirik J, Johnson DS, Kenyon C, Orlin JB, Shor PW, Weber RR (2006) On the sum-of-squares algorithm for bin packing. _J. ACM_ 53(1):1–65.
  * Delimitrou and Kozyrakis (2014) Delimitrou C, Kozyrakis C (2014) Quasar: Resource-efficient and QoS-aware cluster management. _Proc. Int. Conf. Architectural Support for Programming Languages and Operating Systems (ASPLOS)_ , 127–144 (Salt Lake City, UT).
  * Foundation (2023) Foundation AS (2023) Apache mesos: Oversubscription. https://mesos.apache.org/documentation/latest/oversubscription/.
  * Freund and Banerjee (2019) Freund D, Banerjee S (2019) Good prophets know when the end is near. _Available at SSRN: https://ssrn.com/abstract=3479189_ .
  * Ghaderi et al. (2014) Ghaderi J, Zhong Y, Srikant R (2014) Asymptotic optimality of bestfit for stochastic bin packing. _SIGMETRICS Perform. Eval. Rev._ 42(2):64–66.
  * Gupta and Radovanović (2020) Gupta V, Radovanović A (2020) Interior-point-based online stochastic bin packing. _Oper. Res._ 68(5):1474–1492.
  * Kleinrock (1975) Kleinrock L (1975) _Queueing Systems_ (John Wiley & Son).
  * Li et al. (2014) Li Y, Tang X, Cai W (2014) On dynamic bin packing for resource allocation in the cloud. _Proc. Ann. ACM Symp. Parallelism in Algorithms and Architectures (SPAA)_ , 2–11 (Prague, Czech Republic).
  * Lo et al. (2015) Lo D, Cheng L, Govindaraju R, Ranganathan P, Kozyrakis C (2015) Heracles: Improving resource efficiency at scale. _Proc. ACM/IEEE Ann. Int. Symp. Computer Architecture (ISCA)_ , 450–462 (Portland, OR).
  * Meyn (2007) Meyn S (2007) _Control Techniques for Complex Networks_ (USA: Cambridge University Press), 1st edition.
  * Reiss et al. (2012) Reiss C, Tumanov A, Ganger GR, Katz RH, Kozuch MA (2012) Heterogeneity and dynamicity of clouds at scale: Google trace analysis. _Proc. ACM Symp. Cloud Computing (SoCC)_ (San Jose, CA).
  * Stolyar (2013) Stolyar AL (2013) An infinite server system with general packing constraints. _Oper. Res._ 61(5):1200–1217.
  * Stolyar (2017) Stolyar AL (2017) Large-scale heterogeneous service systems with general packing constraints. _Adv. Appl. Probab._ 49:61–83.
  * Stolyar and Zhong (2013) Stolyar AL, Zhong Y (2013) A large-scale service system with packing constraints: Minimizing the number of occupied servers. _ACM SIGMETRICS Perform. Eval. Rev._ 41(1):41–52.
  * Stolyar and Zhong (2015) Stolyar AL, Zhong Y (2015) Asymptotic optimality of a greedy randomized algorithm in a large-scale service system with general packing constraints. _Queueing Syst._ 79:117–143.
  * Stolyar and Zhong (2021) Stolyar AL, Zhong Y (2021) A service system with packing constraints: Greedy randomized algorithm achieving sublinear in scale optimality gap. _Stoch. Syst._ 11:83–111.

##  Appendix A Proof of Theorem 2 (Lower Bound)

###### Proof.

It is sufficient to show that given an infinite-server policy $\sigma$ for
$\mathcal{P}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (2), we have
$N(\sigma)\geq\overline{N}^{*}$. To this end, we will construct a single-
server policy $\overline{\sigma}$ such that the resulting system configuration
$\overline{{\bm{K}}}(\infty)$ in steady state satisfies:

$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]=C(\sigma)\leq\epsilon,$
(13)
$\displaystyle\overline{\lambda}_{i}=\frac{\lambda_{i}r}{N(\sigma)},\quad\forall
i\in\mathcal{I}.$ (14)

Let $\pi$ be the distribution of $\overline{{\bm{K}}}(\infty)$, then
$(N(\sigma),\overline{\sigma},\pi)$ is a feasible solution to the problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (3). As
a result, we have $N(\sigma)\geq\overline{N}^{*}$.

The construction of the single-server policy $\overline{\sigma}$ involves
simulating an infinite-server system under $\sigma$ from the empty
configuration. At time $0$, the policy $\overline{\sigma}$ randomly chooses
the $\ell$-th server in the infinite-server system with probability
$p^{\ell}$, for $\ell=1,2,\cdots$. It then requests jobs for the single-server
system according to a policy $\overline{\sigma}^{\ell}$. The key to our policy
$\overline{\sigma}^{\ell}$ is to make the single-server system emulate the job
assignment at the $\ell$-th server of the simulated infinite-server system,
but without incurring idleness. We first construct the policy
$\overline{\sigma}^{\ell}$, and then specify the probabilities $p^{\ell}$.

Let us start by introducing some useful notation. Let
$\overline{{\bm{K}}}^{\ell}(t)$ be the single-server system configuration
under $\overline{\sigma}^{\ell}$ at time $t$ and ${\bm{K}}^{\ell}(t)$ be the
configuration of the $\ell$-th simulated server in the infinite-server system
under $\sigma$. We define a stochastic process $\\{s^{\ell}(t),t\geq 0\\}$ as
follows:
$s^{\ell}(t)=\max_{\tau}\big{\\{}\tau:\int_{0}^{\tau}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(x)\neq\bm{0}\\}}dx=t\big{\\}}$.
The “$\max$” is well-defined because the integral is continuous in $\tau$.
Intuitively, $s^{\ell}(t)$ gives the maximum time when the accumulative busy
time of the $\ell$-th server is $t$. Note that $\\{s^{\ell}(t),t\geq 0\\}$ is
only discontinuous when ${\bm{K}}^{\ell}(\tau)$ reaches $0$, thus it is right-
differentiable with derivative equal to $1$ at any point.

We construct $\overline{\sigma}^{\ell}$ and the simulation of the infinite-
server system under $\sigma$ in a way such that:

$\overline{{\bm{K}}}^{\ell}(t)={\bm{K}}^{\ell}(s^{\ell}(t))\quad\forall t.$
(15)

That is, we want that the single-server system has the same dynamic of the
simulated $\ell$-th server except skipping the idle period. To this end, we
couple the two systems as follows:

  1. 1.

When the $\ell$-th simulated server ${\bm{K}}^{\ell}(s^{\ell}(t))$ receives a
type $i$ job, we let the single-server system $\overline{{\bm{K}}}^{\ell}(t)$
request a type $i$ job at time $t$. For each such job, its phase transition
process in the $\ell$-th simulated server is the same as that in the single-
server system. That is, when we observe any internal transition or departure
event in $\overline{{\bm{K}}}^{\ell}(t)$, we produce a same event on the
$\ell$-th simulated server ${\bm{K}}^{\ell}(s^{\ell}(t))$.

  2. 2.

The simulations of the rest of the infinite-server system under policy
$\sigma$ are driven by independently generated random seeds.

It is not hard to see that the simulated infinite server-system has the same
stochastic behavior as an uncoupled system under $\sigma.$ Moreover, as we
couple all the events that happen in $\overline{{\bm{K}}}^{\ell}(t)$ and
${\bm{K}}^{\ell}(s^{\ell}(t))$, together with the facts that
$\overline{{\bm{K}}}^{\ell}(t)$ and ${\bm{K}}^{\ell}(s^{\ell}(t))$ are
piecewise constant and
$\overline{{\bm{K}}}^{\ell}(0-)={\bm{K}}^{\ell}(s^{\ell}(0-))=\bm{0}$, we get
(15).

Next we claim that (15) implies the following relationship between the steady-
state cost of the single-server system under $\overline{\sigma}^{\ell}$ and
the steady-state cost of the $\ell$-th simulated server in the infinite-server
system under $\sigma$:

$\mathbb{E}\left[h\left(\overline{{\bm{K}}}^{\ell}(\infty)\right)\right]=\frac{\mathbb{E}\left[h({\bm{K}}^{\ell}(\infty))\right]}{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}.$
(16)

This is because for all ${\bm{k}}\neq\bm{0}$, we have

$\displaystyle\frac{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)={\bm{k}}\right)}{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\lim_{S\to\infty}\frac{\int_{0}^{S}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s)={\bm{k}}\\}}ds}{\int_{0}^{S}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s)\neq\bm{0}\\}}ds}\stackrel{{\scriptstyle(b)}}{{=}}\lim_{T\to\infty}\frac{\int_{0}^{T}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s^{\ell}(t))={\bm{k}}\\}}dt}{\int_{0}^{T}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s^{\ell}(t))\neq\bm{0}\\}}dt}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s^{\ell}(t))={\bm{k}}\\}}dt\stackrel{{\scriptstyle(d)}}{{=}}\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathbbm{1}_{\\{\overline{{\bm{K}}}^{\ell}(t)={\bm{k}}\\}}dt$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{=}}\mathbb{P}\left(\overline{{\bm{K}}}^{\ell}(\infty)={\bm{k}}\right),$

where $(a)$ and $(e)$ hold because long-run averages converge to steady-state
expectations; $(b)$ is due to the fact that
$\int_{0}^{T}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s^{\ell}(t))={\bm{k}}\\}}dt=\int_{0}^{s^{\ell}(T)}\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s)={\bm{k}}\\}}ds,$
for any ${\bm{k}}\neq\bm{0}$; $(c)$ is due to the fact that
$\mathbbm{1}_{\\{{\bm{K}}^{\ell}(s^{\ell}(t))\neq\bm{0}\\}}=1;$ and $(d)$
follows from (15).

Let $\overline{\lambda}_{i}^{\ell}$ be the long-run request rates of type $i$
jobs in the single-server system under $\overline{\sigma}^{\ell}$, and
$\lambda_{i}^{\ell}$ be the throughput of type $i$ jobs of $\ell$-th simulated
server under $\sigma$. By the construction of $\overline{\sigma}^{\ell}$, the
single-server system requests jobs based on the arrival events of the
$\ell$-th simulated server, we have
$\overline{\lambda}_{i}^{\ell}=\frac{\lambda_{i}^{\ell}}{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)},\forall
i\in\mathcal{I}.$

With the constructed policies $\\{\overline{\sigma}^{\ell},\ell=1,2,\dots\\}$,
we are ready to define the policy $\overline{\sigma}$. We let
$\overline{\sigma}$ choose an index $\ell$ with probability $p^{\ell}$ at time
$0$, and then follow $\overline{\sigma}^{\ell}$. We set the probability
$p^{\ell}$ as

$\displaystyle
p^{\ell}=\frac{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}{\sum_{\ell^{\prime}=1}^{\infty}\mathbb{P}\left({\bm{K}}^{\ell^{\prime}}(\infty)\neq\bm{0}\right)}=\frac{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]},\quad\forall\ell=1,2,\dots$
(17)

where the second inequality uses the fact that
$\sum_{\ell^{\prime}=1}^{\infty}\mathbb{P}\left({\bm{K}}^{\ell^{\prime}}(\infty)\neq\bm{0}\right)=\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]$.
Then under $\overline{\sigma}$, we have

$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\right]$
$\displaystyle=\sum_{\ell=1}^{\infty}p^{\ell}\mathbb{E}\left[h\left(\overline{{\bm{K}}}^{\ell}(\infty)\right)\right]\stackrel{{\scriptstyle(a)}}{{=}}\sum_{\ell=1}^{\infty}\frac{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]}\cdot\frac{\mathbb{E}\left[h({\bm{K}}^{\ell}(\infty))\right]}{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}$
$\displaystyle=\frac{\sum_{\ell=1}^{\infty}\mathbb{E}\left[h({\bm{K}}^{\ell}(\infty))\right]}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]}=\frac{\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}(\infty)]}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]}=C(\sigma),$

where $(a)$ follows from (16) and (17). Observe that under
$\overline{\sigma}$, $\overline{{\bm{K}}}(\infty)\neq\bm{0}$ almost surely, we
thus have

$\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]=\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\right]=C(\sigma),$

which proves (13). Moreover, for each $i\in\mathcal{I}$ the request rate
$\overline{\lambda}_{i}$ is given by

$\displaystyle\overline{\lambda}_{i}$
$\displaystyle=\sum_{\ell=1}^{\infty}p^{\ell}\cdot\overline{\lambda}_{i}^{\ell}\stackrel{{\scriptstyle}}{{=}}\sum_{\ell=1}^{\infty}\frac{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]}\cdot\frac{\lambda_{i}^{\ell}}{\mathbb{P}\left({\bm{K}}^{\ell}(\infty)\neq\bm{0}\right)}=\frac{\sum_{\ell=1}^{\infty}\lambda_{i}^{\ell}}{\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}(\infty)]}\stackrel{{\scriptstyle}}{{=}}\frac{\lambda_{i}r}{N(\sigma)}.$

This proves (14). By the argument presented at the beginning of the proof, we
get $N(\sigma)\geq\overline{N}^{*}$. ∎

##  Appendix B Proofs Needed for Theorem 3 (Conversion Theorem)

###  B.1 Proof of Lemma 1

###### Proof.

We will show that under the JRS policy, the Markov chain for the system state
(represented as
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$)
has a unique stationary distribution by first arguing that it is
${\bm{k}}^{0}$-irreducible (here being ${\bm{k}}^{0}$-irreducible means the
Markov chain has a state that can be reached by all other states through
transitions, “${\bm{k}}^{0}$” in “${\bm{k}}^{0}$-irreducible” does not refer
to any specific states), and then use Foster-Lyapunov theorem to show the
positive recurrence (see, e.g., Meyn 2007). Combining
${\bm{k}}^{0}$-irreducibility with positive recurrence, we can conclude that
the Markov chain under study has a unique stationary distribution.

First, we show that the Markov chain
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$
is ${\bm{k}}^{0}$-irreducible. Specifically, observe that the Markov chain
starting from any state
$(({\bm{k}}^{\ell})_{\ell=1,2,\dots},\bm{\zeta}^{1:L},\bm{\eta}^{1:L})$ can
reach the state
$({\bm{k}}^{1:L},(\bm{0})_{\ell=L+1,\dots},\bm{0}^{1:L},\bm{0}^{1:L})$, after
experiencing a sequence of departures and arrivals that clears up all the
tokens, virtual jobs and jobs on backup servers. Further, letting
$\widetilde{{\bm{k}}}$ be the configuration reachable by all other
configuration in the single-server system under the policy
$\overline{\sigma}$, we argue that starting from any states of the form
$({\bm{k}}^{1:L},(\bm{0})_{\ell=L+1,\dots},\bm{0}^{1:L},\bm{0}^{1:L})$, the
Markov chain can reach the state
$(\widetilde{{\bm{k}}}^{1:L},(\bm{0})_{\ell=L+1,\dots},\bm{0}^{1:L},\bm{0}^{1:L})$.
Because for any $\ell\leq L$, there is a transition path from
${\bm{k}}^{\ell}$ to $\widetilde{{\bm{k}}}$, consider the sequence of events
where each ${\bm{K}}^{\ell}(t)$ transitions independently following the path,
and the jobs arrive right after ${\bm{K}}^{\ell}(t)$ making a request, so that
the tokens are checked out before ${\bm{K}}^{\ell}(t)$ has a further
transition. In this way, each ${\bm{K}}^{\ell}(t)$ with $\ell\leq L$ can
eventually reach $\widetilde{{\bm{k}}}$ from ${\bm{k}}^{\ell}$. This proves
the ${\bm{k}}^{0}$-irreducibility of
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$.

Next, we show that
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$
satisfies the Foster-Lyapunov criterion, i.e.,

$\widehat{G}g\leq-1+b\mathbbm{1}_{\\{S\\}},$ (18)

where $g$ is a non-negative function of the states, $S$ is a finite set, $b$
is a finite number, and $\widehat{G}$ is the infinitesimal generator of the
continuous time Markov chain. Let $t_{i}$ be the expected remaining time in
system when a job in in phase $i$ for each $i\in\mathcal{I}$. According to the
job model, we have the recurrence relation

$\Big{(}\mu_{i\perp}+\sum_{i^{\prime}\in\mathcal{I}\colon i^{\prime}\neq
i}\mu_{ii^{\prime}}\Big{)}t_{i}=\sum_{i^{\prime}\in\mathcal{I}\colon
i^{\prime}\neq i}\mu_{ii^{\prime}}t_{i^{\prime}}\quad\forall i\in\mathcal{I}.$
(19)

We construct a Lyapunov function $g$ as follows:

$g(({\bm{k}}^{\ell})_{\ell=1,2,\dots},\bm{\zeta}^{1:L},\bm{\eta}^{1:L})=\sum_{i\in\mathcal{I}}\sum_{\ell=1}^{\infty}t_{i}k_{i}^{\ell}+\sum_{i\in\mathcal{I}}\sum_{\ell=1}^{\infty}t_{i}\zeta_{i}^{\ell}.$
(20)

Using the relation (19), it can be verified that the drift of $g$ satisfies

$\displaystyle\mspace{23.0mu}\widehat{G}g(({\bm{k}}^{\ell})_{\ell=1,2,\dots},\bm{\zeta}^{1:L},\bm{\eta}^{1:L})$
$\displaystyle\leq\sum_{i\in\mathcal{I}}\Big{(}\lambda_{i}t_{i}r-\sum_{\ell=1}^{\infty}k_{i}^{\ell}\Big{)}+\sum_{i\in\mathcal{I}}\Big{(}\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}}^{\ell})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}t_{i}-\sum_{\ell=1}^{\infty}\zeta_{i}^{\ell}\Big{)}$
$\displaystyle\leq\sum_{i\in\mathcal{I}}\Big{(}\lambda_{i}t_{i}r+L\cdot\max_{{\bm{k}}\in\mathcal{K}}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}})}\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}a_{i}t_{i}\Big{)}-\Big{(}\sum_{\ell=1}^{\infty}\sum_{i\in\mathcal{I}}k_{i}^{\ell}+\sum_{\ell=1}^{\infty}\sum_{i\in\mathcal{I}}\zeta_{i}^{\ell}\Big{)},$

where the first inequality uses the fact that the virtual jobs are generated
at a rate no faster than the total rate of job requests. Then the Foster-
Lyapunov criterion in (18) is satisfied with $b$ and $S$ given by

$b=\sum_{i\in\mathcal{I}}\Big{(}\lambda_{i}t_{i}r+L\cdot\max_{{\bm{k}}\in\mathcal{K}}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}})}\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}a_{i}t_{i}\Big{)},$
$S=\left\\{(({\bm{k}}^{\ell})_{\ell=1,2,\dots},\bm{\zeta}^{1:L},\bm{\eta}^{1:L})\colon
g(({\bm{k}}^{\ell})_{\ell=1,2,\dots},\bm{\zeta}^{1:L},\bm{\eta}^{1:L})\leq
b+1\right\\}.$

By the Foster-Lyapunov theorem,
$(({\bm{K}}^{\ell}(t))_{\ell=1,2,\dots},\bm{\zeta}^{1:L}(t),\bm{\eta}^{1:L}(t))$
is positive recurrent. ∎

###  B.2 Proof of Lemma 2

The proof of Lemma 2 follows the framework of Stein’s method (see, for
example, Braverman and Dai 2017, Braverman et al. 2017, Braverman 2022). To
bound the distance between $\overline{{\bm{K}}}^{1:L}$ and
$\widehat{{\bm{K}}}^{1:L}$, observe that because $f\in\text{Lip}(1)$ and
$\sum_{\ell=1}^{L}\sum_{i\in\mathcal{I}}\eta_{i}^{\ell}=O\left(r^{0.5}\right)$,
it suffices to bound bound the Wasserstein distance between
$\overline{{\bm{K}}}^{1:L}$ and $\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}$,
i.e.,

$\sup_{f\in\text{Lip}(1)}\left\\{\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]\right\\}=O\left(r^{0.5}\right),$
(21)

where $f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}$ is a valid
expression because $\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\in\mathcal{K}$ as
discussed in (4.2). Stein’s method usually consists of three steps: generator
comparison, Stein factor bounds, and moment bounds. In our case, due to the
finiteness of the state space $\mathcal{K}$, we only need to do the generator
comparison and the Stein factor bounds. The generator comparison step
justifies the intuition that most normal servers behave like single-server
systems under the policy $\overline{\sigma}$, and the Stein’s facotr bounds
corresponds to the intuition that the small fraction of servers whose dynamics
deviate from single-server systems do not contribute much to the overall
distance
$\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]$.

####  B.2.1 System Dynamics and Generator

To prepare for the proof, we first look into the dynamics of the two systems
under study. In particular, we write out the generators of
$\overline{{\bm{K}}}^{1:L}(t)$ and
$(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t))$, which are used in the
Stein’s method arguments.

We first examine the dynamics of the single-server system under Markovian
policy $\overline{\sigma}$. There are four types of events that change a
single-server system’s configuration: internal transitions, departures,
reactive requests, and proactive requests (see Section 4.1). The change of
configuration due to any event can be represented by the diagram

${\bm{k}}\overset{}{\rightarrow}{\bm{k}}^{\prime}\overset{}{\rightarrow}{\bm{k}}^{\prime}+\bm{a},$

where the arrow ${\bm{k}}\to{\bm{k}}^{\prime}$ denotes an internal transition
or a departure from configuration ${\bm{k}}$ to ${\bm{k}}^{\prime}$ if
${\bm{k}}\neq{\bm{k}}^{\prime}$; the arrow
${\bm{k}}^{\prime}\to{\bm{k}}^{\prime}+\bm{a}$ denotes a reactive request that
adds $\bm{a}$ jobs to the system if ${\bm{k}}\neq{\bm{k}}^{\prime}$, and
denotes a proactive request if ${\bm{k}}={\bm{k}}^{\prime}$. We call the above
change of configuration a transition, and denote its rate as
$\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}$. Let the set of possible
$({\bm{k}}^{\prime},\bm{a})$ pairs in a transition be $E({\bm{k}})$.

We define the total transition rate at configuration ${\bm{k}}$ as
$\gamma_{{\bm{k}}}\triangleq\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}})}\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})},$ and define the
maximal transition rate
$\gamma_{\max}=\max_{{\bm{k}}\in\mathcal{K}}\gamma_{{\bm{k}}}$. Since
$\mathcal{K}$ is a finite set, we must have $\gamma_{\max}<\infty$. Also
observe that the request rate of type $i$ jobs is given by

$\overline{\lambda}_{i}\triangleq\sum_{{\bm{k}}}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}})}\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}a_{i}\cdot\pi({\bm{k}}),$
(22)

where $\pi$ denotes the stationary distribution of single-server configuration
under policy $\overline{\sigma}.$

Next we focus on the dynamics of $L$ i.i.d. copies of single-server systems.
Consider the generator $\overline{G}$ of the corresponding Markov chain
$\\{\overline{{\bm{K}}}^{1:L}(t)\\}$, which is a linear operator on functions
$g\colon\mathcal{K}^{L}\to\mathbb{R}$ defined as:

$\overline{G}g({\bm{k}}^{1:L})\triangleq\frac{d}{dt}\mathbb{E}\left[g\left(\overline{{\bm{K}}}^{1:L}(t)\right)\middle|\overline{{\bm{K}}}^{1:L}(0)={\bm{k}}^{1:L}\right]\Big{|}_{t=0},$
(23)

and we call the resulting function $\overline{G}g(\cdot)$ the drift of
$g(\cdot)$. Based on the transition rates defined above, we have

$\overline{G}g({\bm{k}}^{1:L})=\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}}^{\ell})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(\cdot,{\bm{k}}^{\prime}+\bm{a},\cdot)-g(\cdot,{\bm{k}}^{\ell},\cdot)\right),$
(24)

where $g(\cdot,{\bm{k}}^{\prime}+\bm{a},\cdot)-g(\cdot,{\bm{k}}^{\ell},\cdot)$
is a shorthand for
$g({\bm{k}}^{1},\dots,{\bm{k}}^{\ell-1},{\bm{k}}^{\prime}+\bm{a},{\bm{k}}^{\ell+1},\dots,{\bm{k}}^{L})-g({\bm{k}}^{1:L})$,
i.e., we use $\cdot$ to omit the entries that agree with ${\bm{k}}^{1:L}$.

Similarly, for the infinite-server system, consider the generator
$\widehat{G}$ of $(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t))$ defined as

$\widehat{G}\psi({\bm{k}}^{1:L},\bm{\eta}^{1:L})\triangleq\frac{d}{dt}\mathbb{E}\left[\psi\left(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t)\right)\middle|\widehat{{\bm{K}}}^{1:L}(0)={\bm{k}}^{1:L},{\bm{K}}^{1:L}(0)=\bm{\eta}^{1:L}\right]\Big{|}_{t=0},$
(25)

for any function $\psi\colon(\mathcal{K}\times\mathcal{K})^{L}\to\mathbb{R}$.
Observe that for each $\ell$, the transition of
$(\widehat{{\bm{K}}}^{\ell}(t),\bm{\eta}^{\ell}(t))$ from
$({\bm{k}},\bm{\eta})$ to
$({\bm{k}}^{\prime},\bm{\eta}+\bm{a}\mathbbm{1}_{\\{\bm{\eta}=\bm{0}\\}})$
occurs at the rate $\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}$ for each
$({\bm{k}}^{\prime},\bm{a})\in E({\bm{k}})$, and any arrivals or virtual
arrivals do not change the sum
$\widehat{{\bm{K}}}^{\ell}(t)+\bm{\eta}^{\ell}(t)$. Consider any function
$g\colon\mathcal{K}^{L}\to\mathbb{R}$ and the function
$\psi({\bm{k}}^{1:L},\bm{\eta}^{1:L})=g({\bm{k}}^{1:L}+\bm{\eta}^{1:L})$. The
drift of $\psi$ under $\widehat{G}$ is

$\displaystyle\widehat{G}\psi({\bm{k}}^{1:L},\bm{\eta}^{1:L})$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}}^{\ell})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(\cdot,{\bm{k}}^{\prime}+\bm{a},\cdot)-g(\cdot,{\bm{k}}^{\ell},\cdot)\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$
$\displaystyle\mspace{23.0mu}+\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}}^{\ell})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(\cdot,{\bm{k}}^{\prime}+\bm{\eta}^{\ell},\cdot)-g(\cdot,{\bm{k}}^{\ell}+\bm{\eta}^{\ell},\cdot)\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}.$
(26)

In this context
$g(\cdot,{\bm{k}}^{\prime}+\bm{\eta}^{\ell},\cdot)-g(\cdot,{\bm{k}}^{\ell}+\bm{\eta}^{\ell},\cdot)$
is a shorthand for
$g({\bm{k}}^{1}+\bm{\eta}^{1},\dots,{\bm{k}}^{\ell-1}+\bm{\eta}^{\ell-1},{\bm{k}}^{\prime}+\bm{\eta}^{\ell},{\bm{k}}^{\ell+1}+\bm{\eta}^{\ell+1},\dots,{\bm{k}}^{L}+\bm{\eta}^{L})-g({\bm{k}}^{1:L}+\bm{\eta}^{1:L})$.
In other words, we use $\cdot$ to omit the entries of $g$’s input that agree
with the corresponding entries of ${\bm{k}}^{1:L}+\bm{\eta}^{1:L}$.

######  Remark 1.

In (26), although $g$ is only defined on the domain $\mathcal{K}^{L}$, it is
valid to write ${\bm{k}}^{1:L}+\bm{\eta}^{1:L}$ as its input because we always
have $\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell}\in\mathcal{K}$, i.e., the
total number of real jobs, virtual jobs, and tokens on a normal server never
exceeds $K_{\max}$. To see why this is true, the single-server policy
$\overline{\sigma}$ requests jobs only when there are no tokens on the server,
and it will not request more than $K_{\max}-n$ jobs if there are already $n$
real and virtual jobs on the server.

####  B.2.2 Proof of Lemma 2.

###### Proof.

Step 1: Generator Comparison. For any $f\in\text{Lip}(1)$, consider the
Poisson equation (see, for example, Braverman 2022) that solves for
$g_{f}\colon\mathcal{K}\to\mathbb{R}$:

$\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-f\big{(}{\bm{k}}^{1:L}\big{)}=\overline{G}g_{f}\big{(}{\bm{k}}^{1:L}\big{)},$
(27)

We let ${\bm{k}}^{1:L}=\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}$ in (27) and
take the expectation. This results in

$\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]=\mathbb{E}\left[\overline{G}g_{f}\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right].$
(28)

On the other hand, because $(\widehat{{\bm{K}}}^{1:L}(t),\bm{\eta}^{1:L}(t))$
is a finite-state Markov chain, we have

$\mathbb{E}\left[\widehat{G}\psi_{f}\big{(}\widehat{{\bm{K}}}^{1:L},\bm{\eta}^{1:L}\big{)}\right]=0,$
(29)

where $\psi_{f}$ is given by

$\psi_{f}\big{(}\widehat{{\bm{K}}}^{1:L},\bm{\eta}^{1:L}\big{)}=g_{f}(\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}).$

Subtracting (29) from (28), we get

$\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]=\mathbb{E}\left[\overline{G}g_{f}\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}-\widehat{G}\psi_{f}\big{(}\widehat{{\bm{K}}}^{1:L},\bm{\eta}^{1:L}\big{)}\right].$
(30)

We want to show that $\overline{G}$ and $\widehat{G}$ are close so that we can
bound the RHS of (30).

Now we plug in the formula of the generators in (24) and (26) to the RHS of
(30) and get

$\displaystyle\mspace{23.0mu}\left\lvert\overline{G}g_{f}\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}-\widehat{G}\psi_{f}\big{(}\widehat{{\bm{K}}}^{1:L},\bm{\eta}^{1:L}\big{)}\right\rvert$
$\displaystyle=\Bigg{\lvert}\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E(\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell})}\gamma_{\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g_{f}\big{(}\cdot,{\bm{k}}^{\prime}+\bm{a},\cdot\big{)}-g_{f}\big{(}\cdot,\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},\cdot\big{)}\right)\cdot\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}$
$\displaystyle\mspace{23.0mu}-\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E(\widehat{{\bm{K}}}^{\ell})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g_{f}\big{(}\cdot,{\bm{k}}^{\prime}+\bm{\eta}^{\ell},\cdot\big{)}-g_{f}\big{(}\cdot,\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},\cdot\big{)}\right)\cdot\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}\Bigg{\rvert}$
(31)
$\displaystyle\leq\sum_{\ell=1}^{L}\gamma_{\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell}}\cdot\sup_{({\bm{k}}^{\prime},\bm{a})\in
E(\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell})}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime}+\bm{a},\cdot\big{)}-g_{f}\big{(}\cdot,\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},\cdot\big{)}\right\rvert\cdot\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}$
$\displaystyle\mspace{23.0mu}+\sum_{\ell=1}^{L}\gamma_{\widehat{{\bm{K}}}^{\ell}}\cdot\sup_{({\bm{k}}^{\prime},\bm{a})\in
E(\widehat{{\bm{K}}}^{\ell})}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime}+\bm{\eta}^{\ell},\cdot\big{)}-g_{f}\big{(}\cdot,\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},\cdot\big{)}\right\rvert\cdot\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}$
(32) $\displaystyle\leq
2\gamma_{\max}\cdot\sum_{\ell=1}^{L}\sup_{{\bm{k}}^{\prime}\in\mathcal{K}}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime},\cdot\big{)}-g_{f}\big{(}\cdot,\widehat{{\bm{K}}}^{\ell}+\bm{\eta}^{\ell},\cdot\big{)}\right\rvert\cdot\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}$
(33) $\displaystyle\leq
2\gamma_{\max}\cdot\sup_{{\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K}}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime},\cdot\big{)}-g_{f}\big{(}\cdot,{\bm{k}},\cdot\big{)}\right\rvert\cdot\sum_{\ell=1}^{L}\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}},$
(34)

where in
$g_{f}\big{(}\cdot,{\bm{k}}^{\prime},\cdot\big{)}-g_{f}\big{(}\cdot,{\bm{k}},\cdot\big{)}$
we have omit the entries that agree with ${\bm{k}}^{1:L}$. The equality is
true because each of the $\ell$-th terms in $\overline{G}$ and $\widehat{G}$
is equal if $\bm{\eta}^{\ell}=\bm{0}$ for all $i$. For the first and second
inequalities, recall that $\gamma_{{\bm{k}}}$ is the total transition rate
given by $\gamma_{{\bm{k}}}\triangleq\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}})}\gamma_{{\bm{k}},({\bm{k}}^{\prime},\bm{a})}$. and
$\gamma_{\max}=\max_{{\bm{k}}\in\mathcal{K}}\gamma_{{\bm{k}}}$. Observe that

$\sum_{\ell=1}^{L}\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}\leq\sum_{\ell=1}^{L}\sum_{i\in\mathcal{I}}\eta_{i}^{\ell}\leq\left\lvert\mathcal{I}\right\rvert\cdot\eta_{\max}=O\left(r^{0.5}\right).$

therefore (34) can be further bounded by

$\displaystyle\left\lvert\overline{G}g_{f}\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}-\widehat{G}\psi_{f}\big{(}\widehat{{\bm{K}}}^{1:L},\bm{\eta}^{1:L}\big{)}\right\rvert$
$\displaystyle\leq
2\gamma_{\max}\cdot\sup_{{\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K}}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime},\cdot\big{)}-g_{f}\big{(}\cdot,{\bm{k}},\cdot\big{)}\right\rvert\cdot\sum_{\ell=1}^{L}\sum_{i\in\mathcal{I}}\eta_{i}^{\ell}$
$\displaystyle\leq
2\gamma_{\max}\cdot\sup_{{\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K}}\left\lvert
g_{f}\big{(}\cdot,{\bm{k}}^{\prime},\cdot\big{)}-g_{f}\big{(}\cdot,{\bm{k}},\cdot\big{)}\right\rvert\cdot
O\left(r^{0.5}\right).$

To prove (21), it remains to show that

$\sup_{{\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K}}\left\lvert
g_{f}(\cdot,{\bm{k}}^{\prime},\cdot)-g_{f}(\cdot,{\bm{k}},\cdot)\right\rvert=O\left(1\right).$
(35)

Step 2: Stein Factor Bound. To prove (35), observe that the following
$g_{f}(\cdot)$ is a solution to the Poisson equation (27) (see, for example,
Braverman 2022).

$g_{f}\big{(}{\bm{k}}^{1:L}\big{)}=\int_{0}^{\infty}\mathbb{E}\left[\left(f\big{(}\overline{{\bm{K}}}^{1:L}(t)\big{)}-\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]\right)\middle|\overline{{\bm{K}}}^{1:L}(0)={\bm{k}}^{1:L}\right]dt.$
(36)

This allows us to bound the difference of $g_{f}$ using coupling.
Specifically, we define the coupling of two systems, each consisting of $L$
i.i.d. copies of the single-server system under $\overline{\sigma}$. The two
systems are initialized with configurations $(\cdot,{\bm{k}}^{\prime},\cdot)$
and $(\cdot,{\bm{k}},\cdot)$ that only differ at the $\ell$-th server, where
we omit the entries that agree with ${\bm{k}}^{1:L}$. Let
$\big{(}\overline{{\bm{K}}}^{{1:L},1}(t),\overline{{\bm{K}}}^{{1:L},2}(t)\big{)}$
be the joint configuration of the two systems, which is actually $2L$ i.i.d.
copies of the single-server system. As a result, we can specify the couplings
$(\overline{{\bm{K}}}^{\ell^{\prime},1}(t),\overline{{\bm{K}}}^{\ell^{\prime},2}(t))$
for different $\ell^{\prime}$ separately. For $\ell^{\prime}\neq\ell$, the
corresponding server in the two systems have the same initial configurations,
so we can always keep their configurations identical. For the $\ell$-th
servers, we let them evolve independently following their own dynamics until a
stopping time $\tau_{\text{mix}}$, when their configurations meet. After that,
we can use coupling to keep their configurations identical. Under this
coupling, it is not hard to see that

$\displaystyle\left\lvert
g_{f}(\cdot,{\bm{k}}^{\prime},\cdot)-g_{f}(\cdot,{\bm{k}},\cdot)\right\rvert$
$\displaystyle=\left\lvert\int_{0}^{\infty}\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{{1:L},1}(t)\big{)}-f\big{(}\overline{{\bm{K}}}^{{1:L},2}(t)\big{)}\right]dt\right\rvert$
$\displaystyle\leq\int_{0}^{\infty}\mathbb{E}\left[\left\lvert
f\big{(}\overline{{\bm{K}}}^{{1:L},1}(t)\big{)}-f\big{(}\overline{{\bm{K}}}^{{1:L},2}(t)\big{)}\right\rvert\right]dt$
$\displaystyle\leq\mathbb{E}\left[\int_{0}^{\infty}\left\lvert
f\big{(}\overline{{\bm{K}}}^{{1:L},1}(t)\big{)}-f\big{(}\overline{{\bm{K}}}^{{1:L},2}(t)\big{)}\right\rvert
dt\right]$
$\displaystyle\leq\mathbb{E}\left[\int_{0}^{\infty}\big{\lVert}\overline{{\bm{K}}}^{{1:L},1}(t)-\overline{{\bm{K}}}^{{1:L},2}(t)\big{\rVert}dt\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{\infty}\sum_{\ell^{\prime}=1}^{L}\big{\lVert}\overline{{\bm{K}}}^{\ell^{\prime},1}(t)-\overline{{\bm{K}}}^{\ell^{\prime},2}(t)\big{\rVert}dt\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{\tau_{\text{mix}}}\big{\lVert}\overline{{\bm{K}}}^{\ell,1}(t)-\overline{{\bm{K}}}^{\ell,2}(t)\big{\rVert}dt\right].$
(37)

For each pair of ${\bm{k}},{\bm{k}}^{\prime}$, observe that because
$\overline{\sigma}$ is a ${\bm{k}}^{0}$-irreducible policy,
$\mathbb{E}[\tau_{\text{mix}}]$ is finite; and because $\mathcal{K}$ is a
finite set,
$\big{\lVert}\overline{{\bm{K}}}^{\ell,1}(t)-\overline{{\bm{K}}}^{\ell,2}(t)\big{\rVert}$
is uniformly bounded. All these finite quantities depend on a single-server
system under a policy $\overline{\sigma}$ that is independent of $r$. As a
result, the last expression in (37) is of constant order. Moreover, because
there are finite pairs of $({\bm{k}},{\bm{k}}^{\prime})$, the supremum
$\sup_{{\bm{k}},{\bm{k}}^{\prime}}\mathbb{E}\left[\int_{0}^{\tau_{\text{mix}}}\lVert\overline{{\bm{K}}}^{\ell,1}(t)-\overline{{\bm{K}}}^{\ell,2}(t)\rVert
dt\right]$ is also of constant order, independent of $r$. This proves the
Stein factor bound in (35). Together with the generator comparison, we have
proved

$\sup_{f\in\text{Lip}(1)}\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]=O\left(r^{0.5}\right).$
(21)

The goal of the lemma is a straightforward consequence of (21), as shown
below. Because
$\sum_{\ell=1}^{L}\sum_{i\in\mathcal{I}}\eta_{i}^{\ell}\leq\left\lvert\mathcal{I}\right\rvert\eta_{\max}=O\left(r^{0.5}\right)$,
for any $f\in\text{Lip}(1)$, we have

$\left\lvert\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}+\bm{\eta}^{1:L}\big{)}\right]\right\rvert\leq\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{i\in\mathcal{I}}\eta_{i}^{\ell}(t)\right]=O\left(r^{0.5}\right).$
(38)

Therefore,

$\sup_{f\in\text{Lip}(1)}\mathbb{E}\left[f\big{(}\overline{{\bm{K}}}^{1:L}\big{)}\right]-\mathbb{E}\left[f\big{(}\widehat{{\bm{K}}}^{1:L}\big{)}\right]=O\left(r^{0.5}\right).$
(39)

This finishes the proof.

∎

###  B.3 Proofs of Lemma 3 and Lemma 4

In this subsection, we prove Lemma 3 and Lemma 4 together. We begin by
introducing some notations. We use
$\bm{U}^{\ell}\triangleq(\widehat{{\bm{K}}}^{\ell},\bm{\eta}^{\ell})$ to
represent the state of the $\ell$-th server, and use $\bm{U}^{1:L}$ to
represent the joint state of the first $L$ servers. We also use the lower case
$\bm{u}^{\ell}$, $\bm{u}^{1:L}$ to represent the realizations of the
corresponding random variables. We denote the total number of type $i$ virtual
jobs as $V_{i}\triangleq\sum_{\ell=1}^{L}\zeta_{i}^{\ell}$ for
$i\in\mathcal{I}$, and its realization as $v_{i}$. We denote the total number
of type $i$ jobs on backup servers as $Y_{i}$ for $\in\mathcal{I}$, and its
realizations as $y_{i}$. We also denote the total number of type $i$ tokens
throughout the system as $Z_{i}\triangleq\sum_{\ell=1}^{L}\eta_{i}^{\ell}$,
and its realization as $z_{i}$. Our goal can be rewritten as proving
$\mathbb{E}[V_{i}]=O\left(r^{0.5}\right)$ and
$\mathbb{E}[Y_{i}]=O\left(r^{0.5}\right)$ for each $i\in\mathcal{I}$.

We first give an overview of the proof. Observe that in our model, the
expected time that a jobs stay in the system is fixed. As a result, bounding
the number of virtual jobs or jobs on backup servers in the system is
equivalent to bounding the rate that they are generated, according to Little’s
Law. By our construction of the policy, the rate of generating those jobs are
closely related to the dynamics of the total number of type $i$ tokens
$Z_{i}(t)$.

To describe the dynamics of $Z_{i}(t)$, we first introduce two functions
$dv_{i}$ and $dy_{i}$:

$\displaystyle dv_{i}(a_{i},z_{i})$
$\displaystyle\triangleq(z_{i}+a_{i}-\eta_{\max})^{+},$ $\displaystyle
dy_{i}(z_{i})$ $\displaystyle\triangleq(1-z_{i})^{+}.$

The function $dv_{i}$ represents the increment in the number of type $i$
virtual jobs due to the event that the total number of type $i$ tokens on the
first $L$ servers exceeds the token limit $\eta_{\max}$. The function $dy_{i}$
corresponds to the increment in the total number of type $i$ jobs on backup
servers due to the event that a type $i$ job arrives to the system without
seeing a type $i$ token. For a function
$g\colon(\mathcal{K}\times\mathcal{K})^{L}\to\mathbb{R}$ that only depends on
the number of type $i$ tokens $z_{i}$, its drift can be written as

$\displaystyle\widehat{G}g(\bm{u}^{1:L})$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E({\bm{k}}^{\ell})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(z_{i}+a_{i}-dv_{i}(a_{i},z_{i}))-g(z_{i})\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$
$\displaystyle\mspace{15.0mu}+\lambda_{i}r\left(g(z_{i}-1+dy_{i}(z_{i}))-g(z_{i})\right).$
(40)

We abuse the notation of $g$ here. For ease of exposition, we will simply
write $dv_{i}$ and $dy_{i}$ to represent $dv_{i}(a_{i},z_{i})$ and
$dy_{i}(z_{i})$.

By construction, the total number of type $i$ tokens $\\{Z_{i}(t\\})$ is a
stochastic process constrained within $[0,\eta_{\max}]$. Note that $Z_{i}(t)$
increases when some servers request new tokens, and decreases when a real or
virtual arrival checks out the token or when some servers have the excessive
tokens removed. When $Z_{i}(t)$ is away from the boundaries, the average rate
that it increases is approximately equal to $\lambda_{i}r$, and the average
rate that $Z_{i}(t)$ decreases is given by

$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E(\widehat{{\bm{K}}}^{\ell})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]\approx\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}}^{\prime},\bm{a})\in
E(\overline{{\bm{K}}}^{\ell})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\right]=L\cdot\overline{\lambda}_{i}\approx\lambda_{i}r,$

where we have used the approximations that
$\widehat{{\bm{K}}}^{\ell}\overset{d}{\approx}\overline{{\bm{K}}}^{\ell}$,
$\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\approx 1$ and
$L=\lceil\overline{N}\rceil\approx\overline{N}$.

As $\\{Z_{i}(t)\\}$ randomly moves up and down with approximately the same
rate and reflects on the boundaries of $0$ and $\eta_{\max}$, it behaves as a
reflected simple symmetric random walk. Intuitively speaking, the steady-state
distribution of $Z_{i}$ is approximately an uniform distribution over
$[0,\eta_{\max}]$. Recall that $dv_{i}$ and $dy_{i}$ can only be non-zero when
$Z_{i}(t)$ is near the boundaries. Since the length of the interval
$\eta_{\max}=\Theta\left(r^{0.5}\right)$, we can expect that $dv_{i}$ and
$dy_{i}$ diminish as $r\to\infty$.

In the proof, we first establish the relationship between $\mathbb{E}[V_{i}]$,
$\mathbb{E}[Y_{i}]$ and $dv_{i}$, $dy_{i}$ using Little’s Law. Then we derive
bounds on $dv_{i}$ and $dy_{i}$ by analyzing the drift of several test
functions of $Z_{i}$. This step is implicitly based on the intuition that
$Z_{i}$ is approximately uniform distributed over $[0,\eta_{\max}]$, with the
tokens being generated and eliminated at similar rates. Finally, we invoke
Lemma 2 to show that the tokens are indeed generated and eliminated at similar
speeds, which leads to bounds on $dv_{i}$ and $dy_{i}$.

Finally, we make some additional remarks on the notations. First, $dv_{i}$ and
$dy_{i}$ depend on total number of type $i$ tokens $z_{i}$ and the number of
newly requested jobs $a_{i}$, although we omit the dependency expression for
ease of exposition. Second, we abuse the notation $dv_{i}$ and $dy_{i}$ to
denote the corresponding random variables. We also write
$\sum_{({\bm{k}},\bm{a}^{\prime})}$ as a shorthand for
$\sum_{({\bm{k}}^{\prime},\bm{a})\in E({\bm{k}}^{\ell})}$ when the context is
clear.

####  B.3.1 Proofs of Lemma 3 and Lemma 4.

We are now ready to prove Lemma 3 and Lemma 4.

###### Proof.

Step 1: Bounding Virtual Jobs and Jobs on Backup Servers using Little’s Law.
We first apply Little’s Law to $V_{i}$ and $Y_{i}$. Because the expected time
that a job stays in the system is no more than $\frac{1}{\mu_{\min}}$ with
$\mu_{\min}\triangleq\min_{i}\mu_{i\perp}$, we have

$\displaystyle\mathbb{E}[V_{i}]$
$\displaystyle\leq\frac{1}{\mu_{\min}}\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right],$
(41) $\displaystyle\mathbb{E}[Y_{i}]$
$\displaystyle\leq\frac{1}{\mu_{\min}}\mathbb{E}\left[\lambda_{i}r\cdot
dy_{i}\right].$ (42)

Step 2: Drift Analysis. The above two equations (41) and (42) suggest that we
can derive upper bounds on $\mathbb{E}[V_{i}]$ and $\mathbb{E}[Y_{i}]$ by
analyzing the following two terms:

  * •

$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]$,
interpreted as the average rate that $Z_{i}(t)$ reflects on the boundary at
$\eta_{\max}$;

  * •

$\mathbb{E}\left[\lambda_{i}r\cdot dy_{i}\right]$, interpreted as the average
rate that $Z_{i}(t)$ reflects on the boundary at $0$.

We establish the relationships of $dv_{i}$, $dy_{i}$ and $Z_{i}$ by analyzing
the drift of two test functions $g$.

Letting $g(z_{i})=z_{i}$ and taking steady-state expectation over its drift,
by (B.3) and the fact that the drift is zero in steady state, we get

$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}(a_{i}-dv_{i})\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r(-1+dy_{i})\right]=0.$
(43)

Similarly, letting $g(z_{i})=z_{i}^{2}$ and taking steady-state expectation
over its drift, one can verify that

$\displaystyle\mspace{23.0mu}\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]$
(44)
$\displaystyle=\frac{1}{\eta_{\max}}\cdot\mathbb{E}\left[\left(\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right)\cdot
Z_{i}\right]$ (45)
$\displaystyle\mspace{2.0mu}+\frac{1}{2\eta_{\max}}\cdot\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(a_{i}^{2}-(dv_{i})^{2}\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\cdot(1-(dy_{i})^{2})\right].$
(46)

Readers may refer to the complete calculation at the end of this subsection.

Step 3: Estimating the Terms Obtained from Drift Analysis. We will first focus
on bounding
$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]$
analyzing the two terms in (45) and (46) separately. Then we invoke (43) to
bound $\mathbb{E}\left[\lambda_{i}r\cdot dy_{i}\right]$.

The term in (46) is easy to deal with. Observe that the number of jobs
requested each time should be no more than the maximal number of jobs that a
server can hold, i.e., $a_{i}\leq K_{\max}$, so

$\displaystyle\eqref{eq:sat:ttk-second-order-term2}$
$\displaystyle\leq\frac{1}{2\eta_{\max}}\cdot\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}K_{\max}^{2}+\lambda_{i}r\right]$
$\displaystyle\leq\frac{1}{2\eta_{\max}}\cdot\mathbb{E}\left[\sum_{\ell=1}^{L}\gamma_{\max}K_{\max}^{2}+\lambda_{i}r\right]$
$\displaystyle=O\left(r^{0.5}\right),$ (47)

where in the second inequality we have used the fact that the total rate is
uniformly bounded by $\gamma_{\max}$, and the last step uses the facts that
$L=O\left(r\right)$ and $\eta_{\max}=\Theta\left(r^{0.5}\right)$.

To bound the term in (45), first observe that $Z_{i}\leq\eta_{\max}$, which
implies that

$\eqref{eq:sat:ttk-second-order-
term1}\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right\rvert\right].$

The term on the RHS of the above equation is the expected absolute difference
between the rates of generating and eliminating type $i$ tokens, which is can
be shown to be small relative to $r$. Specifically, we claim that

$\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right\rvert\right]=O\left(r^{0.5}\right).$
(48)

To show (48), first notice that we can remove the indicator
$\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$ without introducing much error:

$\displaystyle\mspace{23.0mu}\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right\rvert\right]$
$\displaystyle\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}\right\rvert\right]$
$\displaystyle\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\gamma_{\max}K_{\max}\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}\right\rvert\right]$
$\displaystyle\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+\gamma_{\max}K_{\max}\left\lvert\mathcal{I}\right\rvert\eta_{\max},$

where the first inequality is due to triangular inequality, the second
inequality is due to the definition of $\gamma_{\max}$, and the last
inequality is due to the fact that
$\sum_{\ell=1}^{L}\mathbbm{1}_{\\{\bm{\eta}^{\ell}\neq\bm{0}\\}}\leq\left\lvert\mathcal{I}\right\rvert\eta_{\max}$.
It remains to bound the term
$\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]$,
which can be seen as showing that the rate of generating type $i$ tokens
concentrates around the type $i$ jobs’ arrival rate $\lambda_{i}r$. It is
natural to think of using some Law of Large Numbers. Unfortunately,
$\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}$
is not a sum of i.i.d. random variables due to dependencies among
$\widehat{{\bm{K}}}^{\ell}$ for different $\ell$’s. As a result, we want to
invoke the Wasserstein distance bound in Lemma 2 to replace
$\widehat{{\bm{K}}}^{\ell}$ in the above expression with
$\overline{{\bm{K}}}^{\ell}$. We define the function $f({\bm{k}}^{1:L})$ as

$f({\bm{k}}^{1:L})=\frac{1}{2\gamma_{\max}K_{\max}}\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert.$
(49)

We claim that $f\in\text{Lip}(1)$. For any two ${\bm{k}}^{{1:L},1}$,
${\bm{k}}^{{1:L},2}$,

$\displaystyle\mspace{23.0mu}2\gamma_{\max}K_{\max}\cdot\left(f({\bm{k}}^{{1:L},1})-f({\bm{k}}^{{1:L},2})\right)$
$\displaystyle=\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell,1},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert-\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell,2},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert$
$\displaystyle\leq\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell,1},({\bm{k}}^{\prime},\bm{a})}a_{i}-\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell,2},({\bm{k}}^{\prime},\bm{a})}a_{i}\right\rvert$
$\displaystyle=\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\left(\gamma_{{\bm{k}}^{\ell,1},({\bm{k}}^{\prime},\bm{a})}-\gamma_{{\bm{k}}^{\ell,2},({\bm{k}}^{\prime},\bm{a})}\right)a_{i}\mathbbm{1}_{\\{({\bm{k}}^{\ell,1})\neq({\bm{k}}^{\ell,2})\\}}\right\rvert$
$\displaystyle\leq\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\left\lvert\gamma_{{\bm{k}}^{\ell,1},({\bm{k}}^{\prime},\bm{a})}-\gamma_{{\bm{k}}^{\ell,2},({\bm{k}}^{\prime},\bm{a})}\right\rvert\cdot
K_{\max}\cdot\lVert{\bm{k}}^{\ell,1}-{\bm{k}}^{\ell,2}\rVert$
$\displaystyle\leq\sum_{\ell=1}^{L}2\gamma_{\max}K_{\max}\cdot\lVert{\bm{k}}^{\ell,1}-{\bm{k}}^{\ell,2}\rVert$
$\displaystyle=2\gamma_{\max}K_{\max}\cdot\lVert{\bm{k}}^{{1:L},1}-{\bm{k}}^{{1:L},2}\rVert,$

where the first inequality is due to triangular inequality; the second
inequality uses the fact that $a_{i}\leq K_{\max}$ and
$\mathbbm{1}_{\\{{\bm{k}}^{1,\ell}\neq{\bm{k}}^{2,\ell}\\}}\leq\lVert{\bm{k}}^{1,\ell}-{\bm{k}}^{2,\ell}\rVert$;
the third inequality uses triangular inequality, the fact that the total rate
at a configuration ${\bm{k}}$ is bounded by $\gamma_{\max}$ and the property
of the $L^{1}$ norm $\lVert\cdot\rVert$. Therefore, $f\in\text{Lip}(1)$. The
Lipschitz continuity of $f$ allows us to invoke Lemma 2 and get

$\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]-\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]\leq
2\gamma_{\max}K_{\max}\cdot O\left(r^{0.5}\right).$

Therefore,

$\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+O\left(r^{0.5}\right).$
(50)

Observe that under a Markovian policy, the request rate of type $i$ jobs can
be written as
$\overline{\lambda}_{i}=\mathbb{E}[\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}]$,
so we have

$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right]=\overline{\lambda}_{i}\cdot\lceil\overline{N}\rceil-\lambda_{i}r=O\left(1\right).$
(51)

Moreover, because
$\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}$
are i.i.d. for $\ell=1,\dots,L$, we have

$\displaystyle\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left[\left(\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right)^{2}\right]}+O\left(1\right)$
$\displaystyle=\sqrt{\sum_{\ell=1}^{L}\text{Var}\left(\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\right)}+O\left(1\right)$
$\displaystyle=O\left(r^{0.5}\right).$

Therefore, combining the arguments above, we get

$\displaystyle\mspace{23.0mu}\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right\rvert\right]$
$\displaystyle\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+O\left(r^{0.5}\right)$
$\displaystyle\leq\mathbb{E}\left[\left\lvert\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\overline{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}-\lambda_{i}r\right\rvert\right]+O\left(r^{0.5}\right)$
$\displaystyle\leq O\left(r^{0.5}\right),$

which proves (48). This implies that the term in (45) is also in
$O\left(r^{0.5}\right)$.

Combining the bounds on the terms in (45) and (46), we get

$\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]=O\left(r^{0.5}\right).$

Finally, we bound $\mathbb{E}\left[\lambda_{i}r\cdot dy_{i}\right]$. We
rearrange the terms in (43) and get

$\displaystyle\mathbb{E}\left[\lambda_{i}r\cdot dy_{i}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}(-a_{i}+dv_{i})\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\right]$
$\displaystyle=\mathbb{E}\left[-\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\right]+O\left(r^{0.5}\right).$

By (48), we have
$\mathbb{E}\left[-\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\right]=O\left(r^{0.5}\right)$.
Therefore,

$\mathbb{E}\left[\lambda_{i}r\cdot dy_{i}\right]=O\left(r^{0.5}\right).$

We invoke the equations (41) and (42) that we get at the beginning of the
proof, and conclude that

$\mathbb{E}[V_{i}]\leq\frac{1}{\mu_{\min}}\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]=O\left(r^{0.5}\right).$
$\mathbb{E}[Y_{i}]\leq\frac{1}{\mu_{\min}}\mathbb{E}\left[\lambda_{i}r\cdot
dy_{i}\right]=O\left(r^{0.5}\right).$

This finishes the proof. ∎

####  B.3.2 Deriving the equality in (44)

We show the calculation detail of deriving the following equality.

$\displaystyle\mspace{23.0mu}\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]$
(44)
$\displaystyle=\frac{1}{\eta_{\max}}\cdot\mathbb{E}\left[\left(\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right)\cdot
Z_{i}\right]$ (45)
$\displaystyle\mspace{2.0mu}+\frac{1}{2\eta_{\max}}\cdot\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(a_{i}^{2}-(dv_{i})^{2}\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\cdot(1-(dy_{i})^{2})\right].$
(46)

The equality is obtained by considering the drift of the function
$g(z_{i})=z_{i}^{2}$, which is zero in steady state. Recall that the drift of
$g(z_{i})$ is given by

$\displaystyle\widehat{G}g(z_{i},\bm{u}^{1:L})$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(z_{i}+a_{i}-dv_{i})-g(z_{i})\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$
$\displaystyle\mspace{15.0mu}+\lambda_{i}r\left(g(z_{i}-1+dy_{i})-g(z_{i})\right).$
(B.3)

We will first calculate $g(z_{i}+a_{i}-dv_{i})-g(z_{i})$, then
$g(z_{i}-1+dy_{i})-g(z_{i})$, and finally plug them into the (B.3).

The calculation of $g(z_{i}+a_{i}-dv_{i})-g(z_{i})$ utilizes the following
property of $dv_{i}$:

$(z_{i}+a_{i}-dv_{i})\cdot dv_{i}=\eta_{\max}\cdot dv_{i}.$ (52)

This property follows from the definition
$dv_{i}=(z_{i}+a_{i}-\eta_{\max})^{+}$. Intuitively, this is because $dv_{i}$
is the “force” that pushes $z_{i}$ back when it hits the boundary at
$\eta_{\max}$. Using the property, we have

$\displaystyle\mspace{23.0mu}g(z_{i}+a_{i}-dv_{i})-g(z_{i})$
$\displaystyle=(z_{i}+a_{i}-dv_{i})^{2}-z_{i}^{2}$
$\displaystyle=(z_{i}+a_{i}-dv_{i})^{2}-(z_{i}+a_{i}-dv_{i}-a_{i}+dv_{i})^{2}$
$\displaystyle=(z_{i}+a_{i}-dv_{i})^{2}-\left((z_{i}+a_{i}-dv_{i})^{2}+2(-a_{i}+dv_{i})\cdot(z_{i}+a_{i}-dv_{i})+(-a_{i}+dv_{i})^{2}\right)$
$\displaystyle=2(a_{i}-dv_{i})\cdot(z_{i}+a_{i}-dv_{i})-(-a_{i}+dv_{i})^{2}$
$\displaystyle=2a_{i}\cdot(z_{i}+a_{i}-dv_{i})-(-a_{i}+dv_{i})^{2}-2dv_{i}\cdot\eta_{\max}$
$\displaystyle=2a_{i}\cdot
z_{i}+a_{i}^{2}-(dv_{i})^{2}-2dv_{i}\cdot\eta_{\max}.$

The second last equality is due to (52), and the rest are all algebraic
manipulations.

We carry out a similar calculation for $g(z_{i}-1+dy_{i})-g(z_{i})$:

$\displaystyle\mspace{23.0mu}g(z_{i}-1+dy_{i})-g(z_{i})$
$\displaystyle=2(-1+dy_{i})\cdot z_{i}+(-1+dy_{i})^{2}$
$\displaystyle=-2z_{i}+2z_{i}\cdot dy_{i}+1-2dy_{i}+(dy_{i})^{2}$
$\displaystyle=-2z_{i}+1+2(z_{i}-1+dy_{i})\cdot dy_{i}-(dy_{i})^{2}$
$\displaystyle=-2z_{i}+1-(dy_{i})^{2}.$