Sharp Waiting-Time Bounds for Multiserver Jobs

1. Introduction

In today’s large-scale computing clusters behind cloud platforms, multiserver jobs have become increasingly prevalent, where a multiserver job is a job that demands to occupy multiple “servers” (which can be multiple physical servers, multiple Central Processing Unit (CPU) cores, etc.) simultaneously during its run time (Verma et al. 2015, Abadi et al. 2016, Lin et al. 2018, Tirmazi et al. 2020). For example, cloud platforms allow users to specify the number of CPU cores in their virtual machines (VMs) or containers, and this information can be utilized by centralized schedulers to make scheduling decisions (see, e.g., Verma et al. 2015 and the Google Kubernetes Engine (Google 2022)). Moreover, the number of “servers” that a multiserver job requests, which we refer to as the server need, is becoming increasingly large. This trend is driven by machine learning jobs from applications like TensorFlow in Abadi et al. (2016), where the jobs are highly parallel and require synchronization. According to the statistics from Google’s Borg Scheduler in Verma et al. (2015), the server needs in Borg can vary across six orders of magnitudes.

In this paper, we study the impact of multiserver jobs on the delay performance of large-scale computing systems using queueing models. Queueing models with multiserver jobs have been studied in the literature, but quantifying the delay performance is notoriously hard. Exact steady-state distributions can only be derived in highly simplified settings with two servers (Brill and Green 1984, Filippopoulos and Karatza 2008), whereas the majority of prior work has focused on characterizing stability conditions (Morozov and Rumyantsev 2016, Rumyantsev and Morozov 2017, Afanaseva et al. 2019, Grosof et al. 2020, Olliaro et al. 2023). However, even for stability, exact conditions are known only for the special cases where all jobs have the same service rate or where there are two job classes. We comment that concurrent to the conference version of our work (Hong and Wang 2022), Grosof et al. (2022a, b) study the delay performance of multiserver jobs in the traditional heavy traffic regime. A more detailed review of the related work is provided in Section 2.

A recent advance in understanding the delay of multiserver jobs is a characterization of the queueing probability in a large system by Wang et al. (2021), where the queueing probability is the probability that an arriving job has to queue rather than entering service immediately. Specifically, Wang et al. (2021) consider a multiserver job system with n servers and study the asymptotic scaling regimes where n becomes large. The scaling regimes allow different job types to have different arrival rates, server needs, and service rates. Among those parameters, server needs and arrival rates can scale up with n. Such scaling regimes capture the trend that different multiserver jobs can be highly heterogeneous, especially in terms of server needs. They establish an upper bound on the queueing probability, based on which they give a sufficient condition for the queueing probability to diminish as n goes to infinity.

Although the work of Wang et al. (2021) identifies when the queueing probability diminishes in large systems, which is a much desirable operating scenario, it does not provide much insight for differentiating between scheduling policies. In particular, their queueing probability upper bound holds for any scheduling policy that is reasonably work conserving (although the bound is presented only for the first come, first serve (FCFS) policy). Moreover, queueing probability does not directly translate to delay of jobs.

In this paper, we focus on the waiting time of jobs, which is the total time a job spends waiting in the queue (not receiving any service), under various scheduling policies. The waiting time is a performance metric that is directly related to job delay. Our goal is to establish bounds on the mean waiting time that are order-wise tight as the number of servers, n, scales. Such tight bounds will enable us to differentiate between policies based on their delay performance. We comment that there has been a line of work in the literature (Liu 2019; Liu and Ying 2020, 2022; van der Boor et al. 2020; Weng and Wang 2020; Weng et al. 2020; Liu et al. 2022) that focuses on quantifying when the mean waiting time diminishes in large systems for various queueing models. However, little is known on how fast the mean waiting time diminishes because of the lack of lower bounds. Our results provide the rate of diminishing when the mean waiting time does diminish, but our tight bounds on the mean waiting time are not limited to the “diminishing” scenario.

Because the first come, first serve policy is widely used as a default policy in practice and also receives the most attention from theoretical studies of multiserver jobs (Brill and Green 1984, Filippopoulos and Karatza 2008, Morozov and Rumyantsev 2016, Rumyantsev and Morozov 2017, Afanaseva et al. 2019, Grosof et al. 2020), in this paper, we will first examine FCFS and understand the exact order of the mean waiting time under it. Then, a natural question that arises is as follows. Can any policy outperform FCFS in terms of the mean waiting time? More generally, we aim to answer the following fundamental questions.

• What is the optimal order of the mean waiting time as the system scales?

• Which policy achieves the optimal order?

1.1. Model and Performance Metric

We consider a system that consists of n servers and I types of jobs. An example is illustrated in Figure 1. Suppose type i jobs need the simultaneous service of ${\mathcal{\ell }}_i$ servers. We sort the job types such that their server needs ${\mathcal{\ell }}_i$’s satisfy ${\mathcal{\ell }}_1\le {\mathcal{\ell }}_2\le \cdots \le {\mathcal{\ell }}_I$. Let the maximal server need ${\mathcal{\ell }}_{\text{max}}$ be ${\mathcal{\ell }}_{\text{max}}={\max }_{i\in \{1,2,\dots ,I\}}{\mathcal{\ell }}_i={\mathcal{\ell }}_I$, and we call type I jobs the maximal-need jobs.

The dynamics of the system are as follows. For each $i=1,2,\dots ,I$, type i jobs arrive to the system following a Poisson process with arrival rate λ_i. Upon arrival, a job either starts service immediately or waits in a centralized queue. When a type i job starts service, it leaves the queue and makes exclusive use of ${\mathcal{\ell }}_i$ servers. The job leaves the system after receiving enough service. The service time of a type i job follows an exponential distribution with service rate μ_i. The service times and arrival events are independent. During the operation of the system, a scheduling policy is used to determine which set of jobs to serve at any time. The scheduling policy is allowed to be preemptive (i.e., we can put a job in service back to the queue and resume its service later).

We measure the performance of our scheduling policy based on mean waiting time as defined below. Let $T_i^w(\infty )$ denote the waiting time of type i jobs in the steady state; then, the mean waiting time is defined as the steady-state expected waiting time averaged over all job types: that is,

where $\lambda \triangleq {\sum }_{i=1}^I{\lambda }_i$ is the total arrival rate.

1.2. Scaling Regimes

We study job delay in scaling regimes where the number of servers, n, goes to infinity. Specifically, we consider a sequence of systems with parameters scaling up jointly with n, and we analyze the growth/decrease rate of the mean waiting time. In the considered scaling regimes, the arrival rates λ_i and server needs ${\mathcal{\ell }}_i$ are allowed to scale with n, whereas the service rate μ_i and the number of job types I stay constant. One key parameter for specifying a scaling regime is the slack capacity δ, defined as $\delta =n-{\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}$, which is the expected number of idle servers in the steady state. Slack capacity is used to specify the heaviness of traffic, which is alternatively specified by load ρ given by $\rho ={\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{n{\mu }_i}$ in the literature.

For expositional purposes, we now parameterize the scaling regimes in a specific way below to demonstrate our results. Our general model is presented in Section 3. Suppose ${\mathcal{\ell }}_{\text{max}}=n^{\gamma },\hspace{0.17em}\delta =n^{\alpha }$ for some exponents $0\le \alpha ,\gamma <1$, and the total arrival rate $\lambda =\mathrm{\Theta }(n)$. We refer to such scaling regimes as polynomial scaling regimes. To classify and visualize different polynomial scaling regimes, we aggregate all scaling regimes with the same $(\alpha ,\gamma )$ pair into one point and plot all such points, as shown in Figure 2. We partition the set of exponent pairs $(\alpha ,\gamma )\in {[0,1)}^2$ into three triangles using the lines $\alpha =\gamma$ and $\alpha =\frac{1+\gamma }{2}$. The corresponding scaling regimes in the upper left in Figure 2 are in general “heavier” than the scaling regimes in the lower right because the former regimes have larger work variability and smaller slack capacity. We comment that the point $(\alpha ,\gamma )=(\frac{1}{2},0)$ is analogous to the celebrated Halfin–Whitt regime introduced in Halfin and Whitt (1981) and that $(\alpha ,\gamma )=(0,0)$ is analogous to the nondegenerate slowdown regime in Atar (2012) in traditional multiclass M/M/n models.

In the polynomial scaling regimes where $\frac{1+\gamma }{2}<\alpha <1$ (indicated in white in Figure 2), Wang et al. (2021) proved that the queueing probability diminishes under FCFS when n gets large. In this paper, we further show that in these regimes, both the queueing probability and the mean waiting time diminish at a fast rate as n gets large under a large class of policies.

Although all reasonable policies perform well when $\frac{1+\gamma }{2}<\alpha <1$, as we move on to polynomial scaling regimes with $\gamma <\alpha <\frac{1+\gamma }{2}$ (marked in blue in Figure 2), different policies begin to show separations in their mean waiting times, as presented in our main results below.

1.3. Results

Here, we present our main results specialized to the polynomial scaling regimes. These results hold under several assumptions: $\gamma <\alpha <\frac{\gamma +1}{2},\hspace{0.17em}\lambda =\mathrm{\Theta }(n)$, and ${\lambda }_I{\mathcal{\ell }}_I=\mathrm{\Theta }(n)$.¹ In Section 4, we will present the general results and comment on how they imply the results here.

• Mean waiting time under FCFS. The exact order of the mean waiting time under FCFS is given by

  $$\mathbb{E}{[T^w(\infty )]}^{\text{FCFS}}=\mathrm{\Theta }\left(n^{\gamma -\alpha }\right).$$(1)

• Mean waiting time lower bound. Under any policy, the mean waiting time satisfies

  $$\mathbb{E}[T^w(\infty )]=\mathrm{\Omega }\left(n^{-\alpha }\right).$$(2)

• Order-wise optimal policy. Consider a static priority policy that we call the smallest-need-first (SNF) policy, which preemptively prioritizes the jobs with smaller server needs. Then, the mean waiting time under SNF achieves the lower bound in (2): that is,

  $$\mathbb{E}{[T^w(\infty )]}^{\text{SNF}}=\mathrm{\Theta }\left(n^{-\alpha }\right).$$(3)

  Therefore, the SNF policy is order-wise optimal in the mean waiting time.

Comparing the mean waiting times under FCFS and under SNF, we can see that FCFS is strictly suboptimal when $\gamma >0$, and SNF improves upon FCFS by a factor of $\mathrm{\Theta }(n^{\gamma })$.

A key to proving the mean waiting time results above is the order-wise tight bounds on the expected workload we establish (Lemmas 1 and 2). In addition, although we consider the system under the traffic regime where $\gamma <\alpha <\frac{1+\gamma }{2}$, we still need to analyze “subsystems” that are in the lighter traffic regime. In this lighter regime, we show that the total server need decays faster than any polynomial (Lemma 8). All these lemmas hold under a very general class of policies, so they could be potentially relevant when we study policies other than FCFS and SNF.

1.3.1. Results on Queueing Probability.

As a by-product to our analysis, we further derive an upper bound on the queueing probability under any work-conserving policy, presented in Corollary 1, which significantly improves upon the queueing probability upper bound in Wang et al. (2021) in a slightly more constrained traffic regime.

1.3.2. Simulation Experiments.

The SNF policy we consider in our analysis is a preemptive priority policy. Preemption is sometimes undesirable in practice. Therefore, we use simulation experiments to also explore a nonpreemptive version of the priority policy, which we call the nonpreemptive smallest-need-first (SNF-NP) policy. SNF-NP serves a job with the smallest server need in the queue when enough numbers of servers free up. Our simulation experiments compare the mean waiting times under FCFS, SNF, and SNF-NP, along with a few other policies. The simulation results, presented in Section 10, show that SNF-NP has comparable performance with SNF, and both of them perform well relative to FCFS and other policies. We also investigate the weighted mean waiting times where jobs with larger server needs have larger weights. We find that SNF and SNF-NP remain competitive when the weights moderately depend on the server needs.

1.4. Technical Challenges

The main technical challenges in analyzing the considered multiserver-job system are rooted in the heterogeneity among job types in both their service rates and their server needs. Such heterogeneity makes the system dynamics multidimensional; neither the total number of jobs in service nor the total number of busy servers determine the current job departure rate. We comment that even for a classical multiclass M/M/n system, where there are multiple job types with different service rates but all job types have a server need of one, finding an optimal scheduling policy is known to be a hard problem, and solutions are available mostly in the so-called Halfin–Whitt heavy traffic regime through the diffusion control problem (Atar et al. 2004, Harrison and Zeevi 2004, Ata and Gurvich 2012). Compared with the classical multiclass M/M/n system, our multiserver-job system has an additional layer of intricacy because of the heterogeneous server needs, which makes it possible for the system to have servers idling while there are jobs waiting in the queue.

To address the challenges because of heterogeneity, our analysis relies on various state-space concentration results. State-space concentration is a phenomenon where the state concentrates around a subset of the state space in steady state, observed in queueing systems in heavy traffic or large-system regimes (Wang et al. 2018; Liu 2019; Liu and Ying 2020, 2022; Weng et al. 2020; Liu et al. 2022). In the multiserver-job system we consider, state-space concentration results are crucial for analyzing the system dynamics when the queue is nonempty. The scenario when the queue is nonempty is especially important to our scaling regimes because the queueing probability may not be diminishing even when the mean waiting time is diminishing. This contrasts with the analysis in the prior work of Wang et al. (2021), which focuses on diminishing queueing probability. Furthermore, our performance goal is to achieve the optimal order of the mean waiting time in large systems, deviating from the traditional performance goal of minimizing delay or certain long-run cost.

1.5. Organization of the Paper and the Relationship with the Conference Version

This paper is organized as follows. In Section 2, we review some additional related work that has not been discussed in Section 1. We present our model and assumptions in Section 3 and then formally state our three main theorems in Section 4. In Section 5, we give an overview of the proof structure and preliminaries of our main proof technique, the drift method. In Section 6, we state and prove Lemmas 1 and 2, which will be used to prove the three theorems in Sections 7, 8, and 9. In Section 10, we present the simulation results. Finally, we conclude the paper in Section 11.

This paper has the following differences from our previous conference version (Hong and Wang 2022). First, we have included a more comprehensive related work section (Section 2). Second, we have included proofs of some important lemmas and a theorem that were omitted because of the space limit in the conference paper. These are the proofs of Lemma 1, Lemma 2, and Theorem 1, which can be found in Section 6 and Section 7. Third, we have changed the name of the order-wise optimal policy that we propose from “P-priority” to “smallest need first” to better reflect the feature of the policy. Fourth, we have included more simulation results.

2. Related Work

In this section, we give a more detailed review of the prior work on multiserver-job models as well as some related models that are not covered in Section 1.

2.1. Multiserver-Job Model

As mentioned in Section 1, the majority of prior work on the multiserver-job model has either focused on characterizing stability conditions (Morozov and Rumyantsev 2016, Rumyantsev and Morozov 2017, Afanaseva et al. 2019, Grosof et al. 2020, Olliaro et al. 2023) or has been restricted to the highly specialized settings with two servers (Brill and Green 1984, Filippopoulos and Karatza 2008). There are relatively fewer papers that study the delay of the multiserver-job model, with different performance metrics. Wang et al. (2021) characterizes the queueing probability in a large multiserver-job model. Zychlinski et al. (2023) studies the optimal cumulative holding cost in a finite-horizon setting. The papers whose performance metrics are closest to our paper are Grosof et al. (2022a, b). Grosof et al. (2022a) characterizes the mean response time in a multiserver-job model under a policy called ServerFilling. Grosof et al. (2022b) proposes a variant of ServerFilling called ServerFilling shortest remaining service time (SRPT), which optimizes the mean response time in the traditional heavy traffic regime. The biggest distinction between their work and our work lies in the scaling regimes. In their work, the analysis of mean response time is asymptotically tight when the load of the system goes to one and the number of servers remains fixed; in contrast, we consider the scaling regimes where the load, number of servers, and server needs scale jointly. Another distinction is the assumptions on the server needs; their work assumes that the server needs can divide the total number of servers, whereas our work assumes that the maximal server need is small compared with the slack capacity.

2.2. Virtual Machine Scheduling

A problem related to the multiserver-job scheduling problem studied in this paper is the VM scheduling problem (see, e.g., Maguluri et al. 2012, 2014; Maguluri and Srikant 2014; Xie et al. 2015; Psychas and Ghaderi 2018, 2019; Stolyar and Zhong 2021). For the VM scheduling problem, typically the system consists of multiple servers, where each server has certain units of each type of resource (e.g., CPU, memory, storage). A VM job demands to occupy multiple units of each type of resource. Each VM job will be served on a single server. Some results for the VM scheduling problem in the traditional heavy traffic regime can be specialized to the multiserver job scheduling problem. To see this, consider a VM scheduling problem where the system consists of a single server and where there is a single resource type. Then, each unit of resource can be viewed as a server in the multiserver-job scheduling problem. With this specialization, the results in Maguluri et al. (2014) provide bounds on a linear combination of the queue lengths of different types of jobs. However, these bounds do not directly translate into bounds on mean job response time.

2.3. Multitask-Job Model

A multitask job is a job that consists of a batch of tasks that can run on servers in parallel, which is similar to a multiserver job in that both can occupy multiple servers at the same time. However, unlike a multiserver job, the tasks of a multitask job can have different run times and do not need to be executed simultaneously. The multitask job model has been considered under a wide variety of settings, and it is sometimes referred to as the batch arrival model (see, e.g., Miller 1959; Daw and Pender 2019; Daw et al. 2019, 2020). Recently, the multitask job model has also been extensively studied under the setting of parallel computing because of the popularity of large-scale data processing systems, such as MapReduce, Apache Hadoop, and Apache Spark (see, e.g., Weng and Wang 2020, Zubeldia 2020). The work closest to our work is Weng and Wang (2020), which shows diminishing queueing time for multiserver jobs in a load-balancing system where tasks of a job need to be dispatched to the queues at the servers upon arrival.

2.4. Dropping Model

When the multiserver-job system does not have any queueing space and allows incoming jobs to be dropped, it becomes a model that has been extensively studied in the literature, and we refer to it as the dropping model. In this model, one can design a dropping policy that decides whether to drop an incoming job or not based on the types of the incoming job and of the jobs currently in service. Under the policy that drops an incoming job only when it cannot fit into the servers (i.e., when its server need is larger than the number of available servers), the stationary distribution has a product form under exponentially distributed service times, as observed by Arthurs and Kaufman (1979). The results have been generalized by Whitt (1985) to allow jobs to demand multiple resource types (e.g., both CPU and input/output) and by van Dijk (1989) to allow general service time distributions. Tikhonenko (2005) further combined aspects of Whitt (1985) and van Dijk (1989). Different dropping policies, which mostly fall within the class of trunk reservation policies, have been designed to minimize the cost associated with dropping (Hunt and Kurtz 1994, Bean et al. 1995, Hunt and Laws 1997).

2.5. Streaming Model

The streaming model for a communication network resembles the multiserver-job model in many aspects. In a streaming model, the “servers” correspond to the bandwidth in the network, and the “jobs” are data flows, such as audio or video flows. Then, flows that require a fixed amount of bandwidth (Melikov 1996, Dasylva and Srikant 1999, Benameur et al. 2001, Ponomarenko et al. 2010), sometimes referred to as streaming flows, can be viewed as multiserver jobs. However, a communication network also features a network structure that the multiserver-job model does not have. A communication network usually has both streaming flows and flows that are flexible in their bandwidth needs, and streaming flows again operate in the dropping model. The performance metric in such a system typically combines the cost associated with dropping for streaming flows and the cost associated with delay for other flows.

3. Model

A basic description of the system parameters and dynamics has been given in Section 1. In this section, we provide formal descriptions of the scheduling policies, the system states, the scaling regime, and the concept of subsystems used in our analysis.

3.1. Scheduling Policies

A scheduling policy decides which jobs to put into service at any moment of time. We are interested in the following two policies.

• FCFS. Jobs are placed onto servers in a first come, first serve fashion until either the next job in queue does not fit or all the jobs are in service.

• SNF. Recall that the job types are indexed in a way such that ${\mathcal{\ell }}_1\le {\mathcal{\ell }}_2\le \cdots \le {\mathcal{\ell }}_I$. We assign priorities to job types such that a smaller index has a higher priority. Whenever there is a job arrival or departure, SNF preempts all the jobs in service and determines a new schedule from scratch. SNF starts from job type 1 and places as many type 1 jobs as possible onto servers. After this, if there are still servers available, SNF goes to the next priority level, type 2, and places as many type 2 jobs as possible onto servers. This procedure continues until no more jobs in the queue can fit into the servers.

3.2. System State

Under FCFS or SNF, a Markovian representation of the system state can be described as follows. The state u of the Markov chain is an ordered list of the jobs in the system, sorted in their order of arrival, and each entry of u describes the type of the corresponding job and whether the job is in service or not. Let the state space be denoted as $\mathcal{U}$. Although the state space is infinite dimensional, in our analysis, we typically only need to focus on three I-dimensional vectors defined below.

For any time t and each job type i, let $X_i(t)$ denote the number of type i jobs in the system, $Z_i(t)$ denote the number of type i jobs in service, and $Q_i(t)\triangleq X_i(t)-Z_i(t)$ denote the number of type i jobs waiting in the queue. Note that because the total number of servers in use cannot exceed n and we cannot serve more jobs than there are in the system, we have the following constraints:

where $[I]$ denotes the index set $\{1,2,\dots ,I\}$.

Let $X_i(\infty ),\hspace{0.17em}{Z}_i(\infty )$, and $Q_i(\infty )$ be random variables that follow the corresponding steady-state distributions when they exist. We sometimes use vector representations of these quantities for convenience. For example, we write $\mathbf{X}(t)=(X_1(t),X_2(t),\dots ,X_I(t))$. We define the vectors $\mathbf{Z}(t),\hspace{0.17em}\mathbf{Q}(t),\hspace{0.17em}\mathbf{X}(\infty ),\hspace{0.17em}\mathbf{Z}(\infty )$, and $\mathbf{Q}(\infty )$ in a similar way. Note that these random elements correspond to the n server system, and thus, their distributions depend on n. Throughout this paper, for conciseness, we often omit the $(\infty )$ in the steady-state random elements except in theorem or lemma statements.

Recall that the waiting time for a job is the total time that the job spends in the queue. Our performance metric is the mean waiting time $\mathbb{E}[T^w(\infty )]$, given by

where $T_i^w(\infty )$ is the waiting time of a type i job in the steady state. Note that the waiting time of a job is the summation of multiple periods if the job is preempted. In any case, we can still apply Little’s law (see, e.g., Harchol-Balter 2013) and rewrite the mean waiting time as

Bounding the mean waiting time is thus equivalent to bounding the expected total queue length.

3.3. Scaling Regimes

Recall that we consider scaling regimes where the number of servers, n, goes to infinity and the arrival rates λ_i and server needs ${\mathcal{\ell }}_i$ are allowed to scale with n, whereas the service rate μ_i and the number of job types I stay constant. The scaling regimes are specified by the slack capacity $\delta \triangleq n-{\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}$, the maximal server need ${\mathcal{\ell }}_{\text{max}}\triangleq {\max }_{i\in [I]}{\mathcal{\ell }}_i={\mathcal{\ell }}_I$, and another parameter called the work variability: ${\sigma }^2\triangleq {\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}$. Work variability reflects the variability of the “work” caused by job arrivals in terms of server-time product, which is $\frac{{\mathcal{\ell }}_i}{{\mu }_i}$ in expectation for each type i job. To help later presentation, we also define the load brought by type i jobs ρ_i as ${\rho }_i=\frac{{\lambda }_i{\mathcal{\ell }}_i}{n{\mu }_i}$.

We state our assumptions below. Throughout the paper, $\log \hspace{0.17em}n$ denotes natural logarithm.

Assumption 1

(Heavy Traffic Assumption). The slack capacity δ is small relative to $\sqrt{{\sigma }^2}$:

$$\delta =o\left(\frac{\sqrt{{\sigma }^2}}{\log n}\right).$$(5)

Assumption 2

(Maximal Server Need Assumption). There exists a constant ${ϵ}_0\in (0,1)$ such that

$${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta .$$(6)

Assumption 3

(Commonness Assumption). The load brought by the maximal-need jobs is not too small (note that ${\mathcal{\ell }}_1\le {\mathcal{\ell }}_2\le \cdots \le {\mathcal{\ell }}_I={\mathcal{\ell }}_{\text{max}}$):

$${\rho }_I\triangleq \frac{{\lambda }_I{\mathcal{\ell }}_I}{n{\mu }_I}=\omega \left(\sqrt{\frac{\delta \hspace{0.17em}\log \hspace{0.17em}n}{\sqrt{{\sigma }^2}}·\frac{{\mathcal{\ell }}_{\text{max}}}{n}}\log \hspace{0.17em}n\right).$$(7)

Assumption 1 guarantees that the traffic is not too light, whereas Assumption 2 guarantees that the system is stable under FCFS and SNF. In the simplified setting of Section 1 where ${\mathcal{\ell }}_{\text{max}}=n^{\gamma }$ and $\delta =n^{\alpha }$, the first two assumptions correspond to $\alpha <\frac{1+\gamma }{2}$ and $\alpha >\gamma$, which exclude the white and gray parts in Figure 2, respectively. Assumption 3 states that the load brought by the maximal-need jobs is not too small. To understand the right-hand side expression in Assumption 3, note that it is automatically satisfied when ${\rho }_I=\omega (\sqrt{\frac{{\mathcal{\ell }}_{\text{max}}}{n}}\log n)$. For example, when ${\mathcal{\ell }}_{\text{max}}=\mathrm{\Theta }(\sqrt{n})$, then it suffices to have ${\rho }_I=\omega (n^{-1/4}\hspace{0.17em}\log \hspace{0.17em}n)$. However, when the traffic becomes heavier (i.e., when $\frac{\delta \hspace{0.17em}\log \hspace{0.17em}n}{\sqrt{{\sigma }^2}}$ becomes smaller), Assumption 3 in (7) can be much weaker than ${\rho }_I=\omega (\sqrt{\frac{{\mathcal{\ell }}_{\text{max}}}{n}}\log n)$.

To have an intuitive view of the magnitudes of the parameters, we give the following asymptotics: ${\sigma }^2=O\left(n{\mathcal{\ell }}_{\text{max}}\right),\hspace{0.17em}\delta =o(\frac{n}{{(\log n)}^2})$, and ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta =o(\frac{n}{{(\log n)}^2})$; they can be verified from the assumptions.

3.4. Subsystems

When we analyze SNF and prove the lower bound for all policies, we frequently use the concept of the ith subsystem, which is the system that has all type j jobs in the original system with $j\le i$ and removes all type k jobs with $k\ge i$. In the ith subsystem, the slack capacity becomes ${\delta }_i=n-{\sum }_{j=1}^i\frac{{\lambda }_j{\mathcal{\ell }}_j}{{\mu }_j}$, and the work variability becomes ${\sigma }_i^2={\sum }_{j=1}^i\frac{{\lambda }_j{\mathcal{\ell }}_j^2}{{\mu }_j^2}$. Note that $\delta ={\delta }_I$ and ${\sigma }_I^2={\sigma }^2$. The maximal server need in the ith system is ${\mathcal{\ell }}_i$ because ${\mathcal{\ell }}_1\le {\mathcal{\ell }}_2\le \cdots \le {\mathcal{\ell }}_i$.

As i increases, the load of the ith subsystem gets heavier because δ_i becomes smaller. There is a critical index $i^{*}$ such that

(i.e., the $i^{*}$th subsystem is the smallest subsystem whose traffic regime is as heavy as that of the original system). Recall that we have assumed $\delta =o(\frac{\sqrt{{\sigma }^2}}{\log n})$, so the set in (8) contains at least the index I; thus, $i^{*}$ is well defined. Note that δ_i is monotonically decreasing, whereas ${\sigma }_i^2$ is monotonically increasing. Thus, the index $i^{*}$ serves as a division point; for any i with $i^{*}\le i\le I$, we have ${\delta }_i=o(\frac{\sqrt{{\sigma }_i^2}}{\log n})$, resulting in a heavier traffic regime, and for any i with $1\le i<i^{*}$, we have ${\delta }_i=\mathrm{\Omega }(\frac{\sqrt{{\sigma }_i^2}}{\log n})$, resulting in a lighter traffic regime.

4. Main Results

In this section, we first present our main results under the scaling regimes we specify in Section 3 as Theorems 1, 2, and 3. Then, we comment on how the theorems subsume the results in Section 1.3.

Theorem 1

(Mean Waiting Time Under FCFS). Consider the multiserver-job system with n servers satisfying Assumptions 1 and 2. Under the FCFS policy, for each $i\in [I]$, the expected waiting time of type i jobs satisfies

$$\mathbb{E}{[T_i^w(\infty )]}^{\text{FCFS}}\ge \frac{{\sigma }^2}{n(\delta +{\mathcal{\ell }}_{\text{max}})}·(1-o(1)),$$(9)

$$\mathbb{E}{[T_i^w(\infty )]}^{\text{FCFS}}\le \frac{{\sigma }^2}{n(\delta -{\mathcal{\ell }}_{\text{max}})}·(1+o(1)).$$(10)

Consequently,

$$\mathbb{E}{[T^w(\infty )]}^{\text{FCFS}}=\mathrm{\Theta }\left(\frac{{\sigma }^2}{n\delta }\right),\hspace{1em}\mathbb{E}{[T_i^w(\infty )]}^{\text{FCFS}}=\mathrm{\Theta }\left(\frac{{\sigma }^2}{n\delta }\right).$$(11)

Theorem 2

(Mean Waiting Time Lower Bound). Consider the multiserver-job system with n servers satisfying Assumptions 1 and 2. Under any policy, the mean waiting time satisfies

$$\mathbb{E}[T^w(\infty )]\ge \underset{i^{*}\le i\le I}{\max }\frac{1}{\lambda }\frac{{\mu }_{\text{min}}{\sigma }_i^2}{{\mathcal{\ell }}_i{\delta }_i}·(1-o(1))=\mathrm{\Omega }\left(\underset{i^{*}\le i\le I}{\max }\frac{1}{\lambda }\frac{{\sigma }_i^2}{{\mathcal{\ell }}_i{\delta }_i}\right),$$(12)

where $i^{*}$ is the critical index defined in (8) and ${\mu }_{\text{min}}={\min }_{i\in [I]}{\mu }_i$, and the expression represented by o(1) is independent of the policies.

Theorem 3

(Mean Waiting Time Under SNF). Consider the multiserver-job system with n servers satisfying Assumptions 1, 2, and 3. Under the SNF policy, the mean waiting time satisfies

$$\mathbb{E}{[T^w(\infty )]}^{\text{SNF}}\le \frac{1}{\lambda }{\displaystyle \sum _{i=i^{*}}^I\frac{{\mu }_{\text{max}}{\sigma }_i^2}{{\mathcal{\ell }}_i({\delta }_i-{\mathcal{\ell }}_i)}}·(1+o(1))=O\left(\underset{i^{*}\le i\le I}{\max }\frac{1}{\lambda }\frac{{\sigma }_i^2}{{\mathcal{\ell }}_i{\delta }_i}\right),$$(13)

where $i^{*}$ is the critical index defined in (8) and ${\mu }_{\text{max}}={\max }_{i\in [I]}{\mu }_i$. Consequently, the SNF policy achieves the optimal order of the mean waiting time.

We have a more general bound for the SNF policy that holds without Assumption 3. Interested readers can refer to Online Appendix A.

The results appearing in Section 1.3 are direct consequences of the above theorems. Recall that in the polynomial scaling regimes, we have ${\mathcal{\ell }}_{\text{max}}=\mathrm{\Theta }(n^{\gamma })$, $\delta =\mathrm{\Theta }(n^{\alpha })$, and $\lambda =\mathrm{\Theta }(n)$. The requirement $\gamma <\alpha <\frac{\gamma +1}{2}$ implies Assumptions 1 and 2; the requirement ${\lambda }_I{\mathcal{\ell }}_I=\mathrm{\Theta }(n)$ implies Assumption 3. Under these assumptions, the bounds in Theorems 1, 2, and 3 hold. It can be verified that ${\sigma }^2=\mathrm{\Theta }(n^{1+\gamma })$, so the expression $\frac{{\sigma }^2}{n\delta }$ in (11) has the order $\mathrm{\Theta }(n^{\gamma -\alpha })$; the critical index $i^{*}=I$, so the expression ${\max }_{i^{*}\le i\le I}\frac{1}{\lambda }\frac{{\sigma }_i^2}{{\mathcal{\ell }}_i{\delta }_i}$ in (12) and (13) equals $\frac{1}{\lambda }\frac{{\sigma }^2}{{\mathcal{\ell }}_{\text{max}}\delta }=\mathrm{\Theta }(n^{-\alpha })$.

5. Proof Road Map and Drift Method Preliminaries

We organize our proofs of the main results as follows; we first prove two important bounds for a quantity called workload given by ${\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i$ in Lemmas 1 and 2, respectively. Then, we convert the workload bounds to the waiting time bounds in Theorems 1 and 2 using properties of FCFS and a linear programming relaxation. For Theorem 3, we analyze SNF by considering each ith subsystems for $i\in [I]$. Some subsystems only need Lemma 1 and Lemma 2, whereas others require the addition of Lemma 8.

Our proof approach is closely related to the recently developed drift method (Eryilmaz and Srikant 2012, Maguluri and Srikant 2016). The drift method allows us to extract information from a continuous-time Markov chain ${\{\mathbf{S}(t)\}}_{t\ge 0}$ on state-space $\mathcal{S}$ by computing the drift of different test functions. Because S(t) is a Markov chain with countable state space and bounded transition rates, we can define the drift of any function f as

We call the operator G the generator of the Markov chain.

For a multiserver-job system, let I-dimensional real vectors $\mathbf{x},\mathbf{z}\in {\mathbb{R}}_{+}^I$ be possible realizations of state descriptors $\mathbf{X}(t)$ and $\mathbf{Z}(t)$, where recall that $\mathbf{X}(t)$ is the vector of the number of jobs in the system at time t and $\mathbf{Z}(t)$ is the vector of the number of jobs in service at time t. We focus on f that only depends on x (i.e., $f:{\mathbb{R}}_{+}^I\to \mathbb{R}$). The drift of f is of the form

where ${\mathbf{e}}_i\in {\mathbb{R}}_{+}^I$ is the vector whose ith entry is one and all the other entries are zero. Note that although f and $Gf$ are functions of the system state u, we write $f(\mathbf{x})$ and $Gf(\mathbf{x},\mathbf{z})$ to highlight the variables that affect their values.

We frequently use the following relation regarding the drift

Heuristically, this is because when $\mathbf{X}(0)$ and $\mathbf{Z}(0)$ follow the stationary distribution, $\mathbf{X}(t)$ and $\mathbf{Z}(t)$ also follow the stationary distribution; so, $f(\mathbf{X}(t))$ and $f(\mathbf{X}(0))$ have the same expectation. Rigorously speaking, this relation only holds for well-behaved functions and Markov processes. The conditions under which the relation holds are discussed in detail in Online Appendix B. Throughout the paper, we assume (16) holds for all f that we consider.

6. Workload Bounds

In this section, we prove two bounds for a quantity called workload given by ${\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i$. These bounds are fundamental to the proofs of the main theorems. In Lemma 1, we give a lower bound on the expected workload applicable to any policy. In Lemma 2, we give upper bounds on the expected workload under any $\delta \prime$-work-conserving policy, a class of policies defined in Definition 1. After stating these two lemmas, we go through preliminaries for the proofs. Finally, we show the full proof of the two lemmas at the end of the section.

Lemma 1

(Workload Lower Bound). Consider the multiserver-job system with n servers satisfying Assumptions 1 and 2. Under any policy, the expected workload is lower bounded as

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}{Q}_i(\infty )]\ge \frac{{\sigma }^2}{\delta }·(1-o(1)),$$(17)

where the expression represented by o(1) is independent of the policies.

Definition 1.

We call a policy $\delta \prime$-work conserving if the following equation holds

$${\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}{Z}_i(t)\ge \min ({\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}{X}_i(t),n-\delta \prime )\hspace{1em}\forall t\ge 0.$$(18)

Here, ${\sum }_{i=1}^I{\mathcal{\ell }}_i{Z}_i(t)$ is equal to the number of busy servers at time t, whereas ${\sum }_{i=1}^I{\mathcal{\ell }}_i{X}_i(t)$, which we call the total server need, is the potential number of busy servers if we can put all jobs at time t into service. Therefore, under a $\delta \prime$-work-conserving policy, either all jobs are in service, or there are at most $\delta \prime$ idling servers. Under any ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy, one can show that the system is stable when ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta$ (Assumption 2) holds. In particular, the system is stable under both FCFS and SNF.

Lemma 2

(Workload Upper Bound). Consider the multiserver-job system with n servers under a $\delta \prime$-work-conserving policy with $\delta \prime \le {ϵ}_0\delta$, where ${ϵ}_0\in (0,1)$ is the parameter in Assumption 2. Then, when $\delta =o(\frac{\sqrt{{\sigma }^2}}{\log n})$,

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}{Q}_i(\infty )]\le \frac{{\sigma }^2}{\delta -\delta \prime }·(1+o(1))=O\left(\frac{{\sigma }^2}{\delta }\right);$$(19)

when $\delta =\mathrm{\Omega }(\frac{\sqrt{{\sigma }^2}}{\log n})$,

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}{Q}_i(\infty )]=O(\sqrt{{\sigma }^2} \log n).$$(20)

Remark 1.

When Assumption 1 is satisfied (i.e., when $\delta =o(\frac{\sqrt{{\sigma }^2}}{\log n})$), the workload upper bound in Lemma 2 coincides with the workload lower bound in Lemma 1 order wise, which implies that the expected workload $\mathbb{E}[{\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i(\infty )]=\mathrm{\Theta }(\frac{{\sigma }^2}{\delta })$. Note that in this case, although the expected workload under all $\delta \prime$-work-conserving policies has the same order, the mean waiting time can vary among policies, as shown for FCFS and SNF in Theorems 1 and 3.

6.1. Preliminaries for Lemma 1 and Lemma 2

Our proofs focus on bounding the normalized work, defined as

where we write ${\overline{x}}_i\triangleq \frac{{\lambda }_i}{{\mu }_i}$ for notational simplicity. We claim that normalized work has the same expectation as the workload (i.e., $\mathbb{E}[\overline{W}]=\mathbb{E}[{\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i]$). To see this, recall that $Q_i=X_i-Z_i$. Now, consider the drift of X_i given by $G{X}_i={\lambda }_i-{\mu }_i{Z}_i$. One can verify that X_i satisfies $\mathbb{E}[G{X}_i]=0$, and thus, $\mathbb{E}[Z_i]=\frac{{\lambda }_i}{{\mu }_i}$. Therefore, the expected workload can be written as

Therefore, bounding the expected workload is equivalent to bounding the steady-state expectation of the normalized work $\mathbb{E}[\overline{W}]$.

We break $\mathbb{E}[\overline{W}]$ into three terms:

where $\overline{r}\in \mathbb{R}$ is up to our choice; ${(\overline{W}-\overline{r})}^{+}\triangleq \max \{\overline{W}-\overline{r},0\}$ denotes the positive part, and ${(\overline{W}-\overline{r})}^{-}\triangleq -\min \{\overline{W}-\overline{r},0\}$ denotes the negative part.

The major difficulty during the proofs is bounding the expectation of the positive part $\mathbb{E}[{(\overline{W}-\overline{r})}^{+}]$. This relies on the relation that $\mathbb{E}[Gf(\mathbf{X},\mathbf{Z})]=0$ as introduced in Section 5. In our proofs, we choose f to be piecewise quadratic functions to bound on the term

Because ${\sum }_{i=1}^I{\mathcal{\ell }}_i(Z_i-{\overline{x}}_i)={\sum }_{i=1}^I{\mathcal{\ell }}_i{Z}_i-n+\delta$, we will be able to bound $\mathbb{E}[{(\overline{W}-\overline{r})}^{+}]$ if we are able to give an accurate estimate of the number of busy servers ${\sum }_{i=1}^I{\mathcal{\ell }}_i{Z}_i$ when the normalized work $\overline{W}\ge \overline{r}$. To get a precise estimate, we exploit the state-space concentration result that says for each $i\in [I]$, X_i cannot be much smaller than ${\overline{x}}_i$ (i.e., ${(X_i-{\overline{x}}_i)}^{-}$ is small with high probability). Formally, this state-space concentration is established by Lemma 3, whose proof uses a sample-path coupling argument and is given in Online Appendix C.

Lemma 3.

Consider the multiserver-job system with n servers. For any nonnegative vector $c=(c_1,\dots ,c_I)\in {\mathbb{R}}_{+}^I$ independent of n, let $c_{\text{max}}={\max }_{i\in [I]}{c}_i,\hspace{0.17em}{\mu }_{\text{max}}={\max }_{i\in [I]}{\mu }_i$, and

$$\mathrm{\Phi }={\displaystyle \sum _{i=1}^I{c}_i}{\mathcal{\ell }}_i(X_i(\infty )-{\overline{x}}_i),$$

where ${\overline{x}}_i=\frac{{\lambda }_i}{{\mu }_i}$. Then, we have the three bounds below.

1. For any $K\ge 0$,

   $$\mathbb{P}\left(\mathrm{\Phi }\le -K\right)\le \exp (-\frac{K^2}{2c_{\text{max}}^2{\mu }_{\text{max}}{\sigma }^2}).$$(23)

2. For any $\alpha \ge 0$ and $\beta \ge 0$ such that $\alpha \beta \ge c_{\text{max}}^2{\mu }_{\text{max}}{\sigma }^2$ and any $j\ge 0$,

   $$\mathbb{P}\left(\mathrm{\Phi }\le -\alpha -\beta j\right)\le e^{-j}.$$(24)

3. Let ${\mathrm{\Phi }}^{-}=\max \{-\mathrm{\Phi },0\}$ to be the negative part of $\mathrm{\Phi }$. Then,

   $$\mathbb{E}\left[{\mathrm{\Phi }}^{-}\right]\le \sqrt{c_{\text{max}}^2{\mu }_{\text{max}}{\sigma }^2}.$$(25)

6.2. Proofs of Lemma 1 and Lemma 2

In this section, we show the full proofs of Lemma 1 and Lemma 2.

Proof of Lemma 1.

Recall that in Section 6.1, we have shown that $\mathbb{E}[{\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i]=\mathbb{E}[\overline{W}]$, where $\overline{W}$ is the normalized work given by $\overline{W}\triangleq {\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}\left(X_i-{\overline{x}}_i\right)$ and ${\overline{x}}_i\triangleq \frac{{\lambda }_i}{{\mu }_i}$. Therefore, our goal is equivalent to lower bounding $\mathbb{E}[\overline{W}]$. To do this, we first perform the following decomposition:

$$\mathbb{E}[\overline{W}]={\overline{r}}_1+\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{+}]-\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{-}].$$

Here, ${\overline{r}}_1$ is a properly chosen small number to be specified later. We bound the positive part $\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{+}]$ and the negative part $\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{-}]$ separately using different techniques.

We bound $\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{+}]$ by analyzing the Lyapunov drift of a function $f:{\mathbb{R}}_{+}^I\to \mathbb{R}$ defined below. Let $\mathbf{x},\mathbf{z}\in {\mathbb{R}}_{+}^I$ denote possible realizations of the state descriptors $\mathbf{X}(t),\mathbf{Z}(t)$. Then, f is defined as

$$f(\mathbf{x})={\phi }_M(w(\mathbf{x})-{\overline{r}}_1),$$

where $w(\mathbf{x})\triangleq {\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}(x_i-{\overline{x}}_i)$ is a possible realization of $\overline{W},\hspace{0.17em}{\phi }_M(s):\mathbb{R}\to \mathbb{R}$ is defined as

$${\phi }_M(s)=\{\begin{array}{ll}0\hfill & \text{if} s\le 0\hfill \\ s^2\hfill & \text{if} 0<s\le M\hfill \\ 2Ms-M^2\hfill & \text{if} s>M\hfill \end{array},$$

and M is a positive number preventing ${\phi }_M(s)$ from growing too fast as s gets large. The derivative of ${\phi }_M(s)$ is ${\phi }_M^{\prime }(s)=2\hspace{0.17em}\min \{s^{+},M\}$.

We will utilize the relation $\mathbb{E}[Gf(\mathbf{X},\mathbf{Z})]=0$, which is implied by Lemma 9 in Online Appendix B if we have $\mathbb{E}[|f(\mathbf{X})|]<\infty$. Because $f(\mathbf{x})$ grows linearly fast as x gets large, $\mathbb{E}[f(\mathbf{x})]<\infty$ follows if we have $\mathbb{E}[X_i]<\infty$ for all $i\in [I]$. On the other hand, the lemma statement holds trivially if $\mathbb{E}[X_i]=\infty$ for some i.

To calculate $Gf(\mathbf{x},\mathbf{z})$, we first decompose the drift $Gf(\mathbf{x},\mathbf{z})$ Equation (15) in the following way:

$$\begin{array}{l}Gf(\mathbf{x},\mathbf{z})={\displaystyle \sum _{i=1}^I\left({\lambda }_i-{\mu }_i{z}_i\right)}\frac{\partial f}{\partial x_i}(\mathbf{x})\hfill \\ +{\displaystyle \sum _{i=1}^I{\lambda }_i}(f(\mathbf{x}+{\mathbf{e}}_i)-f(\mathbf{x})-\frac{\partial f}{\partial x_i}(\mathbf{x}))+{\displaystyle \sum _{i=1}^I{\mu }_i}{z}_i(f(\mathbf{x}-{\mathbf{e}}_i)-f(\mathbf{x})+\frac{\partial f}{\partial x_i}(\mathbf{x})).\hfill \end{array}$$(26)

It is easy to see that the partial derivatives of f appearing in the first term are given by

$$\frac{\partial f}{\partial x_i}(\mathbf{x})=\frac{\partial w}{\partial x_i}(\mathbf{x}){\phi }_M^{\prime }(w(\mathbf{x})-{\overline{r}}_1)=\frac{2{\mathcal{\ell }}_i}{{\mu }_i}\min ({(w(\mathbf{x})-{\overline{r}}_1)}^{+},M).$$

To bound the remaining two terms, observe that

$$\begin{array}{l}f(\mathbf{x}+{\mathbf{e}}_i)-f(\mathbf{x})-\frac{\partial f}{\partial x_i}(\mathbf{x})={\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}(}\frac{\partial f}{\partial x_i}(\xi )-\frac{\partial f}{\partial x_i}(\mathbf{x}))d\xi \\ =\frac{{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}({\phi }_M^{\prime }(w(\xi )-{\overline{r}}_1)-{\phi }_M^{\prime }(w(\mathbf{x})-{\overline{r}}_1))}d\xi \\ \ge \frac{2{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}(w(\xi )-w(\mathbf{x}))}d\xi ·{\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<w(\mathbf{x})<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}\\ =\frac{{\mathcal{\ell }}_i^2}{{\mu }_i^2}·{\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<w(\mathbf{x})<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}},\end{array}$$

where the inequality is because of the fact that for any $s_1\ge s_2$,

$${\phi }_M^{\prime }(s_1)-{\phi }_M^{\prime }(s_2)\ge 2(s_1-s_2)·{\mathbf{1}}_{\{0\le s_2\le s_1\le M\}},$$(27)

which can be verified by brute force calculation. Similarly,

$$\begin{array}{l}f(\mathbf{x}-{\mathbf{e}}_i)-f(\mathbf{x})+\frac{\partial f}{\partial x_i}(\mathbf{x})={\displaystyle {\int }_{\xi \in [\mathbf{x}-{\mathbf{e}}_i,\mathbf{x}]}(}-\frac{\partial f}{\partial x_i}(\xi )+\frac{\partial f}{\partial x_i}(\mathbf{x}))d\xi \\ =\frac{{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x}-{\mathbf{e}}_i,\mathbf{x}]}(}-{\phi }_M^{\prime }(w(\xi )-{\overline{r}}_1)+{\phi }_M^{\prime }(w(\mathbf{x})-{\overline{r}}_1))d\xi \\ \ge \frac{2{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x}-{\mathbf{e}}_i,\mathbf{x}]}(-w(\xi )+w(\mathbf{x}))}d\xi ·{\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<w(\mathbf{x})<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}\\ =\frac{{\mathcal{\ell }}_i^2}{{\mu }_i^2}·{\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<w(\mathbf{x})<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}\\ \ge \frac{{\mathcal{\ell }}_i^2}{{\mu }_i^2}{\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<w(\mathbf{x})<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}.\end{array}$$

Plugging the above inequalities into the decomposition of the drift (26) and taking expectation on both sides because $\mathbb{E}[Gf(\mathbf{X},\mathbf{Z})]=0$, we have

$$\begin{array}{ll}0\hfill & \ge \mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left({\overline{x}}_i-Z_i\right)·2\min ({(\overline{W}-{\overline{r}}_1)}^{+},M)+{\displaystyle \sum _{i=1}^I(}\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i}){\mathbf{1}}_{\{{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}<\overline{W}<{\overline{r}}_1+M-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}].\hfill \end{array}$$(28)

The inequality still holds when we let $M\to \infty$ inside the expectation:

$$0\ge \mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left({\overline{x}}_i-Z_i\right)·2{(\overline{W}-{\overline{r}}_1)}^{+}+{\displaystyle \sum _{i=1}^I(}\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i}){\mathbf{1}}_{\{\overline{W}>{\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}\}}].$$(29)

Here, we are implicitly exchanging ${\lim }_{M\to \infty }$ and $\mathbb{E}$, which is legal by the dominated convergence theorem given that the random variable inside the expectation is dominated by another random variable with a finite expectation, $2n{\mathcal{\ell }}_{\text{max}}{(\overline{W}-{\overline{r}}_1)}^{+}+2n{\mathcal{\ell }}_{\text{max}}/{\mu }_{\text{min}}$. Using the facts that $\mathbb{E}[{\sum }_{i=1}^I(\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i})]=2{\sigma }^2$ and ${\sum }_{i=1}^I(\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i})\le \frac{2n{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}$, we get

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_1)}^{+}]\ge {\sigma }^2-\frac{n{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}·\mathbb{P}(\overline{W}\le {\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}).$$(30)

By Lemma 3 with $\mathrm{\Phi }=\overline{W},\hspace{0.17em}{c}_i=\frac{1}{{\mu }_i},\hspace{0.17em}K=2\sqrt{{\mu }_{\text{max}}{\sigma }^2\log n}/{\mu }_{\text{min}}$, we have $\mathbb{P}\left(\overline{W}\le -K\right)\le \frac{1}{n^2}.$ Note that this K satisfies $K=O(\sqrt{{\sigma }^2}\log n)$. We take ${\overline{r}}_1=-K-\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}}$; then, $\mathbb{P}(\overline{W}\le {\overline{r}}_1+\frac{{\mathcal{\ell }}_{\text{max}}}{{\mu }_{\text{min}}})\le \frac{1}{n^2}$. Moreover, observe that ${\sum }_{i=1}^I{\mathcal{\ell }}_i(Z_i-{\overline{x}}_i)\le n-(n-\delta )=\delta$, so

$$\mathbb{E}\left[{(\overline{W}-{\overline{r}}_1)}^{+}\right]\ge \frac{{\sigma }^2}{\delta }·(1-o(1)).$$(31)

The lower bound for negative part $\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{-}]$ follows from Lemma 3(c), which shows that $\mathbb{E}[{(\overline{W})}^{-}]=O(\sqrt{{\sigma }^2})$. Because $s\mapsto s^{-}$ is a monotonically decreasing function and ${\overline{r}}_1\le 0$,

$$\mathbb{E}[{(\overline{W}-{\overline{r}}_1)}^{-}]\le \mathbb{E}[{(\overline{W})}^{-}]=O(\sqrt{{\sigma }^2}).$$(32)

Therefore, combining (31), (32), and the fact that ${\overline{r}}_1=-O(\sqrt{{\sigma }^2}\log n+{\mathcal{\ell }}_{\text{max}})$, we have

$$\mathbb{E}[\overline{W}]=\frac{{\sigma }^2}{\delta }·(1-o(1))-O(\sqrt{{\sigma }^2}\log n+{\mathcal{\ell }}_{\text{max}})=\frac{{\sigma }^2}{\delta }·(1-o(1)),$$

and this finishes the proof of Lemma 1. □

Proof of Lemma 2.

Recall that in Section 6.1, we have shown that $\mathbb{E}[{\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}{Q}_i]=\mathbb{E}[\overline{W}]$, where $\overline{W}$ is the normalized work given by $\overline{W}\triangleq {\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}\left(X_i-{\overline{x}}_i\right)$ and ${\overline{x}}_i\triangleq \frac{{\lambda }_i}{{\mu }_i}$. Moreover, we have

$$\mathbb{E}[\overline{W}]\le {\overline{r}}_2+\mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}]$$(33)

for any number ${\overline{r}}_2$. Next, we will bound the term $\mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}]$ for a suitably chosen ${\overline{r}}_2$.

We bound $\mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}]$ by analyzing the Lyapunov drift of a function $f:{\mathbb{R}}_{+}^I\to \mathbb{R}$ defined below. Let $\mathbf{x},\mathbf{z}\in {\mathbb{R}}_{+}^I$ denote possible realizations of the state descriptors $\mathbf{X}(t),\mathbf{Z}(t)$. We define f as

$$f(\mathbf{x})=\phi \left(w(\mathbf{x})-{\overline{r}}_2\right),$$

where $\phi (s)={(s^{+})}^2$, and $w(\mathbf{x})\triangleq {\sum }_{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}(x_i-{\overline{x}}_i)$.

This proof relies on the relation that $\mathbb{E}[Gf(\mathbf{X},\mathbf{Z})]=0$, which is justified by Lemma 11 in Online Appendix B. To calculate $Gf(\mathbf{x},\mathbf{z})$, we decompose the drift Equation (15) in the following way:

$$\begin{array}{l}Gf(\mathbf{x},\mathbf{z})={\displaystyle \sum _{i=1}^I\left({\lambda }_i-{\mu }_i{z}_i\right)}\frac{\partial f}{\partial x_i}(\mathbf{x})\hfill \\ +{\displaystyle \sum _{i=1}^I{\lambda }_i}(f(\mathbf{x}+{\mathbf{e}}_i)-f(\mathbf{x})-\frac{\partial f}{\partial x_i}(\mathbf{x}))+{\displaystyle \sum _{i=1}^I{\mu }_i}{z}_i(f(\mathbf{x}-{\mathbf{e}}_i)-f(\mathbf{x})+\frac{\partial f}{\partial x_i}(\mathbf{x})).\hfill \end{array}$$(34)

It is easy to see that $\frac{\partial f}{\partial x_i}(\mathbf{x})$ in the first term of (34) is given by

$$\frac{\partial f}{\partial x_i}(\mathbf{x})=\frac{\partial r}{\partial x_i}(\mathbf{x})\phi \prime (w(\mathbf{x})-{\overline{r}}_2)=\frac{2{\mathcal{\ell }}_i}{{\mu }_i}{(w(\mathbf{x})-{\overline{r}}_2)}^{+}.$$

The other two terms in (34) can be bounded by constants independent of x and z as below:

$$\begin{array}{l}f(\mathbf{x}+{\mathbf{e}}_i)-f(\mathbf{x})-\frac{\partial f}{\partial x_i}(\mathbf{x})={\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}(}\frac{\partial f}{\partial x_i}(\xi )-\frac{\partial f}{\partial x_i}(\mathbf{x}))d\xi \\ =\frac{{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}(}\phi \prime (w(\xi )-{\overline{r}}_2)-\phi \prime (w(\mathbf{x})-{\overline{r}}_2))d\xi \\ \le \frac{2{\mathcal{\ell }}_i}{{\mu }_i}{\displaystyle {\int }_{\xi \in [\mathbf{x},\mathbf{x}+{\mathbf{e}}_i]}(w(\xi )-w(\mathbf{x}))}d\xi \\ =\frac{{\mathcal{\ell }}_i^2}{{\mu }_i^2},\end{array}$$

where the inequality is because $\phi \prime (s)=2s^{+}$ is 2-Lipschitz continuous. Similarly,

$$f(\mathbf{x}-{\mathbf{e}}_i)-f(\mathbf{x})+\frac{\partial f}{\partial x_i}(\mathbf{x})\le \frac{{\mathcal{\ell }}_i^2}{{\mu }_i^2}.$$

Because $\mathbb{E}[Gf(\mathbf{X},\mathbf{Z})]=0$, plugging the last inequality into the decomposition of the drift, (34), and taking expectation on both sides of (34), we get

$$0\le \mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left({\overline{x}}_i-Z_i\right)·2{(\overline{W}-{\overline{r}}_2)}^{+}]+\mathbb{E}[{\displaystyle \sum _{i=1}^I(}\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i})].$$(35)

Observe that because $\mathbb{E}[Z_i]=\frac{{\lambda }_i}{{\mu }_i}$, the second term in (35) can be straightforwardly computed as

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I(}\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}+\frac{Z_i{\mathcal{\ell }}_i^2}{{\mu }_i})]=2{\displaystyle \sum _{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i^2}{{\mu }_i^2}}=2{\sigma }^2.$$(36)

Rearranging the terms, we get the following key equation; for any real number ${\overline{r}}_2$,

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}]\le {\sigma }^2.$$(37)

Now, suppose we are able to show that

$$\gamma ·\mathbb{E}\left[{(\overline{W}-{\overline{r}}_2)}^{+}\right]\le \mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}]+o(1),$$(38)

for some $\gamma >0$; then, by (37), $\mathbb{E}\left[{(\overline{W}-{\overline{r}}_2)}^{+}\right]\le \frac{{\sigma }^2+o(1)}{\gamma }$, so $\mathbb{E}[\hspace{0.17em}\overline{W}\hspace{0.17em}]\le {\overline{r}}_2+\frac{{\sigma }^2+o(1)}{\gamma }$.

We devote the remainder of this proof to proving (38). The idea here is to use the following two events to further partition the probability space:

$${\mathcal{E}}_1=\{{\displaystyle \sum _{i=1}^I(}\frac{1}{{\mu }_{\text{min}}}-\frac{1}{{\mu }_i}){\mathcal{\ell }}_i(X_i-{\overline{x}}_i)>-K_2\},\hspace{1em}{\mathcal{E}}_2=\{\overline{W}\le \frac{n}{{\mu }_{\text{min}}}\},$$

where K₂ is a suitable number such that ${\mathcal{E}}_1$ happens with high probability. We break the term on the left of (37) based on the three cases ${\mathcal{E}}_1,\hspace{0.17em}{\mathcal{E}}_1^c\cap {\mathcal{E}}_2$, and ${\mathcal{E}}_1^c\cap {\mathcal{E}}_2^c$ to analyze them separately:

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}]=\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1\}}]$$(39)

$$+\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2\}}]$$(40)

$$+\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2^c\right\}}].$$(41)

• Case 1. Event ${\mathcal{E}}_1$ happens. Observe that the term in (39) is nonzero only when $\overline{W}>{\overline{r}}_2$ and event ${\mathcal{E}}_1$ happens, which implies that

  $$\begin{array}{l}{\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}(X_i-{\overline{x}}_i)>{\overline{r}}_2,\\ {\displaystyle \sum _{i=1}^I}(\frac{1}{{\mu }_{\text{min}}}-\frac{1}{{\mu }_i}){\mathcal{\ell }}_i(X_i-{\overline{x}}_i)>-K_2.\end{array}$$

  Adding up the above two inequalities and rearranging the terms yield

  $${\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(X_i-{\overline{x}}_i)>{\mu }_{\text{min}}({\overline{r}}_2-K_2).$$(42)

  When the above inequality holds, we can invoke the definition of the $\delta \prime$-work-conserving policy and the fact that ${\sum }_{i=1}^I{\mathcal{\ell }}_i{\overline{x}}_i=n-\delta$ to get

  $$\begin{array}{l}{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)\ge \min ({\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}{X}_i,n-\delta \prime )-{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}{\overline{x}}_i\hfill \\ =\min ({\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(X_i-{\overline{x}}_i),\delta -\delta \prime )\ge \min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime ).\hfill \end{array}$$(43)
  $$\begin{array}{l}\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1\right\}}]\hfill \\ \ge \min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )·\mathbb{E}[{\left(\overline{W}-{\overline{r}}_2\right)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1\right\}}].\hfill \end{array}$$(44)

• Case 2. Event ${\mathcal{E}}_1^c\cap {\mathcal{E}}_2$ happens. To bound the term in (40), we need to analyze the probability of event ${\mathcal{E}}_1^c$. Consider Lemma 3 with $\mathrm{\Phi }={\sum }_{i=1}^I(\frac{1}{{\mu }_{\text{min}}}-\frac{1}{{\mu }_i}){\mathcal{\ell }}_i(X_i-{\overline{x}}_i),\hspace{0.17em}{c}_i=\frac{1}{{\mu }_{\text{min}}}-\frac{1}{{\mu }_i},\hspace{0.17em}\alpha =\beta =\sqrt{\frac{{\mu }_{\text{max}}}{{\mu }_{\text{min}}^2}{\sigma }^2},\hspace{0.17em}j=3\hspace{0.17em}\log \hspace{0.17em}n$. It can be verified that $\alpha \beta \ge c_{\text{max}}^2{\mu }_{\text{max}}{\sigma }^2$. Let $K_2=\alpha +\beta j$. Then, we have

  $$\mathbb{P}\left({\mathcal{E}}_1^c\right)=\mathbb{P}({\displaystyle \sum _{i=1}^I(}\frac{1}{{\mu }_{\text{min}}}-\frac{1}{{\mu }_i}){\mathcal{\ell }}_i(X_i-{\overline{x}}_i)\le -K_2)\le \frac{1}{n^3}.$$

  Note that this choice of K₂ satisfies $K_2=O(\sqrt{{\sigma }^2}\log n)$.

  With the upper bound $\mathbb{P}\left({\mathcal{E}}_1^c\right)\le \frac{1}{n^3}$, we bound the term in (40) in the following way. Observe that ${\sum }_{i=1}^I{\mathcal{\ell }}_i(Z_i-{\overline{x}}_i)\ge -n$. Furthermore, when the event ${\mathcal{E}}_2$ occurs, $\overline{W}\le \frac{n}{{\mu }_{\text{min}}}$ and $0\le {(\overline{W}-{\overline{r}}_2)}^{+}\le \frac{n}{{\mu }_{\text{min}}}$. Consequently, we have the following two inequalities:

  $$\begin{array}{l}\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(Z_i-{\overline{x}}_i)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2\right\}}]\ge -n·\frac{n}{{\mu }_{\text{min}}}·\mathbb{P}\left({\mathcal{E}}_1^c\right)\ge -\frac{1}{{\mu }_{\text{min}}n},\\ \mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2\right\}}]\le \frac{n}{{\mu }_{\text{min}}}·\mathbb{P}\left({\mathcal{E}}_1^c\right)\le \frac{1}{{\mu }_{\text{min}}{n}^2}.\end{array}$$

  We can rearrange the above inequalities into a similar form as (44):

  $$\begin{array}{ll}\hfill & \mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(Z_i-{\overline{x}}_i)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2\right\}}]\hfill \\ \hfill & \hspace{1em}\hspace{1em}\ge \min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )·\mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2\right\}}]-\frac{n+\delta }{{\mu }_{\text{min}}{n}^2},\hfill \end{array}$$(45)
  where we have used the fact that $\min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )\le \delta$.

• Case 3. Event ${\mathcal{E}}_1^c\cap {\mathcal{E}}_2^c$ happens. Lastly, we bound the term in (41). Observe that ${\mathcal{E}}_2^c$ implies that ${\sum }_{i=1}^I{\mathcal{\ell }}_i{X}_i\ge n$, so a $\delta \prime$-work-conserving policy will make sure that ${\sum }_{i=1}^I{\mathcal{\ell }}_i{Z}_i\ge n-\delta \prime$. Therefore,

  $$\begin{array}{l}\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(Z_i-{\overline{x}}_i)·{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2^c\right\}}]\ge (\delta -\delta \prime )·\mathbb{E}[{(\overline{W}-{\overline{r}}_2)}^{+}{\mathbf{1}}_{\left\{\overline{W}>{\overline{r}}_2,{\mathcal{E}}_1^c,{\mathcal{E}}_2^c\right\}}].\hfill \end{array}$$(46)

Combining the results in three cases (44), (45), and (46), we have

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}\left(Z_i-{\overline{x}}_i\right)·{(\overline{W}-{\overline{r}}_2)}^{+}]\ge \min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )·\mathbb{E}\left[{(\overline{W}-{\overline{r}}_2)}^{+}\right]-\frac{n+\delta }{{\mu }_{\text{min}}{n}^2}.$$(47)

Therefore, we have shown (38) with $\gamma =\min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )$. By (37), we conclude that

$$\mathbb{E}[\overline{W}]\le {\overline{r}}_2+\mathbb{E}\left[{(\overline{W}-{\overline{r}}_2)}^{+}\right]\le {\overline{r}}_2+\frac{{\sigma }^2}{\min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )}+O\left(\frac{1}{n}\right).$$(48)

Next, we simplify (48) based on the order of δ. When $\delta =o\left(\frac{\sqrt{{\sigma }^2}}{\log n}\right)$, we take ${\overline{r}}_2=K_2+\frac{\delta -\delta \prime }{{\mu }_{\text{min}}}$ and get

$$\mathbb{E}[\overline{W}]\le K_2+\frac{\delta -\delta \prime }{{\mu }_{\text{min}}}+\frac{{\sigma }^2}{\delta -\delta \prime }+O\left(\frac{1}{n}\right)=\frac{{\sigma }^2}{\delta -\delta \prime }·(1+o(1)),$$

where we have used the fact that $K_2=O(\sqrt{{\sigma }^2}\log n)$ and $\delta =o\left(\frac{\sqrt{{\sigma }^2}}{\log n}\right)$.

When $\delta =\mathrm{\Omega }(\frac{\sqrt{{\sigma }^2}}{\log n})$, we have $\delta \ge C\frac{\sqrt{{\sigma }^2}}{\log n}$ for some C > 0 and n large enough. Because $\delta -\delta \prime \ge (1-{ϵ}_0)\delta$, we get $\delta -\delta \prime \ge (1-{ϵ}_0)C\frac{\sqrt{{\sigma }^2}}{\log n}$ for n large enough. Taking ${\overline{r}}_2=K_2+\frac{(1-{ϵ}_0)C}{{\mu }_{\text{min}}}\frac{\sqrt{{\sigma }^2}}{\log n}$ yields $\min ({\mu }_{\text{min}}({\overline{r}}_2-K_2),\delta -\delta \prime )\ge (1-{ϵ}_0)C\frac{\sqrt{{\sigma }^2}}{\log n}$, so

$$\mathbb{E}[\overline{W}]\le K_2+\frac{(1-{ϵ}_0)C\sqrt{{\sigma }^2}}{{\mu }_{\text{min}}\log n}+\frac{\sqrt{{\sigma }^2}\log n}{(1-{ϵ}_0)C}+O\left(\frac{1}{n}\right)=O(\sqrt{{\sigma }^2}\log n).$$

This finishes the proof. □

7. Proof of Theorem 1 (Waiting Times Under FCFS)

7.1. Proof Overview and Lemmas

We prove Theorem 1 in this section. We first present some intuition and state two lemmas. Then, we give a proof of the theorem using the lemmas. We leave the proofs of the lemmas to Online Appendix D.

The analysis of FCFS is based on the intuition that, under FCFS, the number of type i jobs in the queue is approximately proportional to its arrival rate λ_i: that is,

where $Q_{\mathrm{\Sigma }}\triangleq {\sum }_{i=1}^I{Q}_i$ is the total queue length. Using this relation, we can easily convert from the expected workload to mean waiting time:where the second equality is because of Little’s law and the fact that ${\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}=n-\delta$.

This intuition is formalized by considering a modified-FCFS policy, which serves the jobs in an FCFS order but will not serve the next job until there are at least ${\mathcal{\ell }}_{\text{max}}$ idle servers. In Lemma 4, we prove that under modified FCFS, (49) holds exactly, so the expected queue length of each type of jobs can be computed using the above argument. Moreover, when n is large, we expect modified FCFS to be a good approximation of FCFS. In fact, we can define an upper-bounding system as the multiserver-job system with n servers under modified FCFS and define a lower-bounding system as the multiserver-job system with $n+{\mathcal{\ell }}_{\text{max}}$ servers under modified FCFS. In Lemma 5, we prove that the queue length of each type of jobs in the original system is sandwiched between those in the two modified systems. Because the two modified systems themselves are very close to each other, we get a tight characterization of FCFS.

Lemma 4.

For both the upper-bounding system and the lower-bounding system, we have

$$\mathbb{E}[Q_i^{(U)}(\infty )]=\frac{{\lambda }_i}{\lambda }\mathbb{E}[Q_{\mathrm{\Sigma }}^{(U)}(\infty )]\hspace{1em}\forall i\in [I],$$(51)

$$\mathbb{E}[Q_i^{(L)}(\infty )]=\frac{{\lambda }_i}{\lambda }\mathbb{E}[Q_{\mathrm{\Sigma }}^{(L)}(\infty )]\hspace{1em}\forall i\in [I],$$(52)

where $Q_i^{(U)}(\infty )$ and $Q_i^{(L)}(\infty )$ denote the steady-state queue lengths of type i jobs in the upper- and lower-bounding systems, respectively; $Q_{\mathrm{\Sigma }}^{(U)}(\infty )\triangleq {\sum }_{i=1}^I{Q}_i^{(U)}(\infty )$ and $Q_{\mathrm{\Sigma }}^{(L)}(\infty )\triangleq {\sum }_{i=1}^I{Q}_i^{(L)}(t)$ denote the total queue lengths in the steady state.

Lemma 5.

$$\mathbb{E}[Q_i^{(L)}(\infty )]\le \mathbb{E}[Q_i(\infty )]\le \mathbb{E}[Q_i^{(U)}(\infty )]\hspace{1em}\forall i\in [I],$$(53)

where $Q_i^{(U)}(\infty )$ and $Q_i^{(L)}(\infty )$ denote the steady-state queue lengths of type i jobs in the upper- and lower-bounding systems, respectively.

7.2. Proof of Theorem 1

Proof of Theorem 1.

Recall that by Little’s law, $\mathbb{E}[T_i^w]=\frac{1}{{\lambda }_i}\mathbb{E}[Q_i]$ and $\mathbb{E}[T^w]=\frac{1}{\lambda }{\sum }_{i=1}^I\mathbb{E}[Q_i]$. To characterize $\mathbb{E}[Q_i]$, Lemma 5 suggests that we only need to characterize $\mathbb{E}[Q_i^{(L)}]$ and $\mathbb{E}[Q_i^{(U)}]$, the queue lengths in the lower-bounding system and the upper-bounding system, for each $i\in [I]$.

Observe that the lower-bounding system has $n+{\mathcal{\ell }}_{\text{max}}$ servers, with slack capacity $\delta +{\mathcal{\ell }}_{\text{max}}$. Because $\delta +{\mathcal{\ell }}_{\text{max}}\le (1+{ϵ}_0)\delta =o(\frac{\sqrt{{\sigma }^2}}{\log n})$, we can apply Lemmas 1 and 4 to get

$$\frac{1}{\lambda }\mathbb{E}[Q_{\mathrm{\Sigma }}^{(L)}]·{\displaystyle \sum _{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}}={\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}\mathbb{E}[Q_i^{(L)}]\ge \frac{{\sigma }^2}{\delta +{\mathcal{\ell }}_{\text{max}}}·(1-o(1)).$$

Recall the facts introduced in Section 3 that $\delta =o(n)$ and $\delta =n-{\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}$, so we have ${\sum }_{i=1}^I\frac{{\lambda }_i{\mathcal{\ell }}_i}{{\mu }_i}=n·(1-o(1))$. Therefore, the above equation implies that $\mathbb{E}[Q_{\mathrm{\Sigma }}^{(L)}]\ge \frac{\lambda {\sigma }^2}{n(\delta +{\mathcal{\ell }}_{\text{max}})}·(1-o(1))$. By Lemma 4,

$$\mathbb{E}[Q_i^{(L)}]\ge \frac{{\lambda }_i{\sigma }^2}{n(\delta +{\mathcal{\ell }}_{\text{max}})}·(1-o(1))\hspace{1em}\forall i\in [I].$$(54)

For the upper-bounding system, observe that it has n server and slack capacity δ and that it operates under an ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy. Applying the workload upper bound in Lemma 2 yields

$$\frac{1}{\lambda }{\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}\mathbb{E}[Q_i^{(U)}]\le \frac{{\sigma }^2}{\delta -{\mathcal{\ell }}_{\text{max}}}·(1+o(1)).$$

Following a similar argument as in the lower-bounding system, we have

$$\mathbb{E}[Q_i^{(U)}]\le \frac{{\lambda }_i{\sigma }^2}{n(\delta -{\mathcal{\ell }}_{\text{max}})}·(1+o(1))\hspace{1em}\forall i\in [I].$$(55)

Therefore, by Lemma 5, we have

$$\frac{{\lambda }_i{\sigma }^2}{n(\delta +{\mathcal{\ell }}_{\text{max}})}·(1-o(1))\le \mathbb{E}[Q_i]\le \frac{{\lambda }_i{\sigma }^2}{n(\delta -{\mathcal{\ell }}_{\text{max}})}·(1+o(1))\hspace{1em}\forall i\in [I],$$(56)

which implies the waiting time bounds in Theorem 1. □

8. Proof of Theorem 2 (Mean Waiting Time Lower Bound)

Proof of Theorem 2.

Recall that by Little’s law, $\mathbb{E}[T^w]=\frac{1}{\lambda }{\sum }_{i=1}^I\mathbb{E}[Q_i]$. Therefore, it suffices to show a lower bound on the total queue length. We fix a policy in the original system. For any i with $i^{*}\le i\le I$, we consider the ith subsystem by ignoring all job types with index greater than i. In the ith subsystem, it is always possible to achieve the same $\mathbb{E}[Q_j]$’s for $j\le i$ by imitating the service decisions taken by the original system. Therefore, we have ${\sum }_{j=1}^i\frac{{\mathcal{\ell }}_j}{{\mu }_j}\mathbb{E}[Q_j]\ge \frac{{\sigma }_i^2}{{\delta }_i}·(1-o(1))$, where the right-hand side expression is the workload lower bound of the ith subsystem according to Lemma 1. Then, the expected waiting time $\mathbb{E}[T^w]$ is lower bounded by the optimal value of the following linear programming problem:

$$\underset{\{q_j:j\in [I]\}}{\min }\hspace{1em}\frac{1}{\lambda }{\displaystyle \sum _{j=1}^I{q}_j}$$

$$\text{subject} \text{to}\hspace{1em}{\displaystyle \sum _{j=1}^i\frac{{\mathcal{\ell }}_j}{{\mu }_j}}{q}_j\ge \frac{{\sigma }_i^2}{{\delta }_i}·(1-o(1))\hspace{1em}{i}^{*}\le i\le I$$

$$q_j\ge 0\hspace{1em}\forall j\in [I],$$

where q_j corresponds to $\mathbb{E}[Q_j]$. It is easy to see that the objective value satisfies

$$\frac{1}{\lambda }{\displaystyle \sum _{j=1}^I{q}_j}\ge \frac{1}{\lambda }\frac{{\mu }_{\text{min}}}{{\mathcal{\ell }}_i}{\displaystyle \sum _{j=1}^i\frac{{\mathcal{\ell }}_j}{{\mu }_j}}{q}_j\ge \frac{{\mu }_{\text{min}}{\sigma }_i^2}{\lambda {\mathcal{\ell }}_i{\delta }_i}·(1-o(1))$$(57)

for any $i^{*}\le i\le I$, where in the first inequality, we have used the fact that ${\mu }_j\ge {\mu }_{\text{min}}$ and ${\mathcal{\ell }}_j\le {\mathcal{\ell }}_i$ for any $j\le i$. Note that when q_j is zero for each j with $j\ne i$, the first inequality becomes an equality up to a constant order factor in terms of μ_i and ${\mu }_{\text{min}}$. Because i in (57) can be an arbitrary integer satisfying $i^{*}\le i\le I$, the optimal objective value is no less than

$$\underset{i^{*}\le i\le I}{\max }\frac{{\mu }_{\text{min}}{\sigma }_i^2}{\lambda {\mathcal{\ell }}_i{\delta }_i}·(1-o(1)).$$(58)

Therefore,

$$\mathbb{E}[T^w]\ge \underset{i^{*}\le i\le I}{\max }\frac{{\mu }_{\text{min}}{\sigma }_i^2}{\lambda {\mathcal{\ell }}_i{\delta }_i}·(1-o(1)).$$

This completes the proof. □

Remark 2.

The proof of the lower bound provides some intuitions for choosing the SNF policy. We consider the simple case where $i^{*}=I$. By (57), the total queue length order wise achieves the lower bound when the queue consists of jobs with the largest server needs, suggesting that we should give low priorities to those jobs. This is in a similar spirit to the famous SRPT policy for single-server jobs, which leaves the jobs with large remaining service times in the queue (see, e.g., Harchol-Balter 2013).

9. Proof of Theorem 3 (Mean Waiting Time Under SNF)

To understand the behavior under the SNF policy, one key observation is that for each $i\in [I]$, the type i jobs are unaffected by type j jobs with j > i. As a result, we can learn about the original system by analyzing each ith subsystem, which is obtained by removing all jobs of type j with j > i. Some subsystems are under relatively heavier traffic or more precisely, subject to ${\delta }_i=o(\frac{\sqrt{{\sigma }_i^2}}{\log n})$, whereas some subsystems are under lighter traffic. For those subsystems under relatively heavier traffic, Lemmas 1 and 2 are enough for use; for those subsystems under lighter traffic, we sometimes use Lemma 8 below, which is a more refined bound on the expected total server need $\mathbb{E}[{\sum }_{i=1}^I{\mathcal{\ell }}_i{X}_i]$, proven based on the two tail bounds in Lemmas 6 and 7. The proofs of Lemma 6, Lemma 7, and Lemma 8 are in Online Appendix E.

Lemma 6.

Consider the multiserver-job system with n servers satisfying ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta$ (Assumption 2). Letting ${\overline{x}}_i=\frac{{\lambda }_i}{{\mu }_i}$, under any ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy, the normalized work has the following tail bound; for any ϵ such that $0<ϵ<{ϵ}_0$, there exists ${\alpha }_1=\frac{2ϵ\delta }{{\mu }_{\text{min}}}+\mathrm{\Theta }\left(\frac{n{\mathcal{\ell }}_{\text{max}}}{\delta }\log n\right)$ and ${\beta }_1=\mathrm{\Theta }\left(\frac{n{\mathcal{\ell }}_{\text{max}}}{\delta }\right)$ such that for any $j\ge 0$,

$$\mathbb{P}({\displaystyle \sum _{i=1}^I\frac{{\mathcal{\ell }}_i}{{\mu }_i}}(X_i(\infty )-{\overline{x}}_i)\ge {\alpha }_1+{\beta }_1·j)\le e^{-j}.$$(59)

Lemma 7.

Consider the multiserver-job system with n servers satisfying ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta$ (Assumption 2). Letting ${\overline{x}}_i=\frac{{\lambda }_i}{{\mu }_i}$, under any ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy, the total server need has the following tail bound; there exists α₂ and β₂ with ${\alpha }_2=\frac{\delta }{2}+\mathrm{\Theta }\left(\frac{n{\mathcal{\ell }}_{\text{max}}}{\delta }\log n\right)$ and ${\beta }_2=\mathrm{\Theta }\left(\frac{n{\mathcal{\ell }}_{\text{max}}}{\delta }\right)$ such that for any $j\ge 0$,

$$\mathbb{P}({\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(X_i(\infty )-{\overline{x}}_i)\ge {\alpha }_2+{\beta }_2·j)\le e^{-j}.$$(60)

The proof technique of the two lemmas is using state-space concentration successively. Lemma 6 relies on the state-space concentration implied by Lemma 3, whereas Lemma 7 relies on the state-space concentration implied by Lemma 3 and Lemma 6.

As a consequence of the first two lemmas, we can give a bound on the expectation of the total server need $\mathbb{E}[{\sum }_{i=1}^I{\mathcal{\ell }}_i{X}_i]$ in a different traffic regime than what is assumed in Assumption 1. This is useful for analyzing the dynamics of subsystems under SNF.

Lemma 8

(Total Server Need Upper Bound Under Lighter Traffic). Consider the multiserver-job system with n servers satisfying ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta$ (Assumption 2), and $\delta =\omega (\sqrt{n{\mathcal{\ell }}_{\text{max}}}\log n)$. Letting ${\overline{x}}_i=\frac{{\lambda }_i}{{\mu }_i}$, under any ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy, the expected total server need has the following upper bound:

$$\mathbb{E}[{\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}(X_i(\infty )-{\overline{x}}_i)]=\exp (-\mathrm{\Omega }(\frac{{\delta }^2}{n{\mathcal{\ell }}_{\text{max}}}\left)\right).$$(61)

As a quick digression, with Lemma 7, we can prove an upper bound on the queueing probability in Corollary 1. We comment that this queueing probability bound significantly improves on the bound in Wang et al. (2021) in a slightly more constrained traffic regime.

Corollary 1

(Queueing Probability). Consider the multiserver-job system with n servers satisfying ${\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta$ (Assumption 2), and $\delta =\omega (\sqrt{n{\mathcal{\ell }}_{\text{max}}}\log n)$. Then, the steady-state probability that a job experiences queueing has the following upper bound:

$$\mathbb{P}({\displaystyle \sum _{i=1}^I{\mathcal{\ell }}_i}{X}_i(\infty )\ge n)=\exp (-\mathrm{\Omega }(\frac{{\delta }^2}{n{\mathcal{\ell }}_{\text{max}}}\left)\right).$$(62)

Now, we are ready to prove Theorem 3.

Proof of Theorem 3.

Recall that by Little’s law, we have that $\mathbb{E}\left[T^w\right]=\frac{1}{\lambda }{\sum }_{i=1}^I\mathbb{E}[Q_i]$, and thus, it suffices to bound $\mathbb{E}[Q_i]$’s. We bound $\mathbb{E}[Q_i]$ separately for each $i\in [I]$ as follows:

$$\mathbb{E}\left[Q_i\right]\le \frac{1}{{\mathcal{\ell }}_i}\mathbb{E}\left[{\displaystyle \sum _{j=1}^i{\mathcal{\ell }}_j}{Q}_j\right]\le \frac{1}{{\mathcal{\ell }}_i}\mathbb{E}\left[{\displaystyle \sum _{j=1}^i{\mathcal{\ell }}_j}\left(X_j-{\overline{x}}_j\right)\right],$$(63)

where we have used the fact that Q_i’s are nonnegative, $Q_j=X_j-Z_j$, and $\mathbb{E}[Z_j]=\frac{{\lambda }_j}{{\mu }_j}={\overline{x}}_j$. Observe that under SNF, the dynamics of the first i types of jobs are unaffected by the rest of the jobs. Therefore, the expectation $\mathbb{E}[{\sum }_{j=1}^i{\mathcal{\ell }}_j(X_j-{\overline{x}}_j)]$ can be viewed as the expected total server need in the ith subsystem, which has slack capacity δ_i, maximal server need ${\mathcal{\ell }}_i$, and work variability ${\sigma }_i^2$. We discuss the bound on $\mathbb{E}[Q_i]$ in three cases based on different relationships of δ_i, ${\mathcal{\ell }}_i$, and ${\sigma }_i^2$.

• Case 1. ${\delta }_i=o(\frac{\sqrt{{\sigma }_i^2}}{\log n})$. Applying Lemma 2 to the ith subsystem, we have

  $$\mathbb{E}\left[{\displaystyle \sum _{j=1}^i{\mathcal{\ell }}_j}\left(X_j-{\overline{x}}_j\right)\right]\le {\mu }_{\text{max}}\mathbb{E}\left[{\displaystyle \sum _{j=1}^i\frac{{\mathcal{\ell }}_j}{{\mu }_j}}\left(X_j-{\overline{x}}_j\right)\right]\le \frac{{\mu }_{\text{max}}{\sigma }_i^2}{{\delta }_i-{\mathcal{\ell }}_i}·(1+o(1)).$$

  Therefore,

  $$\mathbb{E}[Q_i]\le \frac{{\mu }_{\text{max}}{\sigma }_i^2}{{\mathcal{\ell }}_i({\delta }_i-{\mathcal{\ell }}_i)}·(1+o(1)).$$

• Case 2. ${\delta }_i=\mathrm{\Omega }\left(\frac{\sqrt{{\sigma }_i^2}}{\log n}\right)$ and ${\mathcal{\ell }}_i=\omega (\frac{\delta \log n}{\sqrt{{\sigma }^2}}{\mathcal{\ell }}_{\text{max}})$. Applying Lemma 2 to the ith subsystem, we have

  $$\mathbb{E}\left[{\displaystyle \sum _{j=1}^i{\mathcal{\ell }}_j}\left(X_j-{\overline{x}}_j\right)\right]\le {\mu }_{\text{max}}\mathbb{E}\left[{\displaystyle \sum _{j=1}^i\frac{{\mathcal{\ell }}_j}{{\mu }_j}}\left(X_j-{\overline{x}}_j\right)\right]=O\left(\sqrt{{\sigma }_i^2}\log n\right).$$

  Combining the above equation with (63) and the facts that ${\mathcal{\ell }}_i=\omega \left(\frac{\delta \log n}{\sqrt{{\sigma }^2}{\mathcal{\ell }}_{\text{max}}}\right)$ and ${\sigma }_i^2\le {\sigma }^2$, we have

  $$\mathbb{E}[Q_i]\le O\left(\frac{\sqrt{{\sigma }_i^2}\log n}{{\mathcal{\ell }}_i}\right)=o(\frac{{\sigma }^2}{\delta }·\frac{1}{{\mathcal{\ell }}_{\text{max}}}).$$

• Case 3. ${\delta }_i=\mathrm{\Omega }(\frac{\sqrt{{\sigma }_i^2}}{\log n})$ and ${\mathcal{\ell }}_i=O\left(\frac{\delta \log n}{\sqrt{{\sigma }^2}}{\mathcal{\ell }}_{\text{max}}\right)$. We first show that ${\delta }_i=\omega (\sqrt{n{\mathcal{\ell }}_i}\log n)$. To see this, observe that ${\delta }_i=\mathrm{\Omega }(\frac{\sqrt{{\sigma }_i^2}}{\log n})$ implies i < I, so ${\delta }_i\ge \frac{{\lambda }_I{\mathcal{\ell }}_I}{{\mu }_I}=n{\rho }_I$. By Assumption 3, we have

  $${\delta }_i\ge n{\rho }_I=\omega \left(\sqrt{\frac{\delta \log n}{\sqrt{{\sigma }^2}}·n{\mathcal{\ell }}_{\text{max}}}\log n\right)=\omega \left(\sqrt{n{\mathcal{\ell }}_i}\log n\right).$$(64)

Therefore, we can apply Lemma 8 to the ith subsystem to get

$$\mathbb{E}\left[{\displaystyle \sum _{j=1}^i{\mathcal{\ell }}_j}\left(X_j-{\overline{x}}_j\right)\right]=\exp (-\mathrm{\Omega }(\frac{{\delta }_i^2}{n{\mathcal{\ell }}_i}\left)\right)\le \exp (-\mathrm{\Omega }({(\log n)}^2)),$$

where we have used ${\delta }_i=\omega (\sqrt{n{\mathcal{\ell }}_i}\log n)$ in the last inequality. Therefore, $\mathbb{E}[Q_i]\le \exp (-\mathrm{\Omega }({(\log n)}^2))$, which decays faster than any polynomial in n.

Combining the three cases, we get

$$\mathbb{E}\left[T^w\right]=\frac{1}{\lambda }{\displaystyle \sum _{i=1}^I\mathbb{E}}[Q_i]\le \frac{1}{\lambda }{\displaystyle \sum _{i=i^{*}}^I\frac{{\mu }_{\text{max}}{\sigma }_i^2}{{\mathcal{\ell }}_i({\delta }_i-{\mathcal{\ell }}_i)}}·(1+o(1)),$$

where we have used the fact that $\mathbb{E}\left[T^w\right]=\frac{1}{\lambda }{\sum }_{i=1}^I\mathbb{E}[Q_i]$. Here, the summation is taken from $i^{*}$ to I because Case 1 is equivalent to $i^{*}\le i\le I$. The second inequality is because the other two terms are of lower order compared with $\frac{{\mu }_{\text{max}}{\sigma }_I^2}{{\mathcal{\ell }}_I({\delta }_I-{\mathcal{\ell }}_I)}$.

Finally, note that I and ${\mu }_{\text{max}}$ are independent of n during scaling, and ${\mathcal{\ell }}_i\le {\mathcal{\ell }}_{\text{max}}\le {ϵ}_0\delta \le {ϵ}_0{\delta }_i$ for some ${ϵ}_0<1$. Therefore,

$$\frac{1}{\lambda }{\displaystyle \sum _{i=i^{*}}^I\frac{{\mu }_{\text{max}}{\sigma }_i^2}{{\mathcal{\ell }}_i({\delta }_i-{\mathcal{\ell }}_i)}}=\mathrm{\Theta }\left(\underset{i^{*}\le i\le I}{\max }\frac{1}{\lambda }\frac{{\sigma }_i^2}{{\mathcal{\ell }}_i{\delta }_i}\right).$$

This finishes the proof. □

10. Simulation Results

10.1. Experiment 1: Comparing SNF and SNF-NP with FCFS

Our theoretical analysis shows that there is an order gap in mean waiting times under FCFS and SNF for the scaling regimes satisfying Assumptions 1, 2, and 3. Here, we complement our asymptotic results by comparing the simulated performances of SNF and FCFS in finite systems. We also simulate the nonpreemptive variant of SNF, which we call SNF-NP. SNF-NP identifies jobs in the queue that can fit into the idle servers and prioritizes those with the smallest server needs for service.

We run the simulation experiments under two sets of parameters.² Parameter Set One satisfies all three assumptions, whereas Parameter Set Two does not satisfy Assumption 3. The parameters are specified in Figure 3. We can see that both SNF and SNF-NP outperform FCFS by a large margin. In particular, the ratio between the mean waiting times of FCFS and SNF is always greater than 3 when $n\ge 2^{10}$, and it gets to as large as 40.0 under Parameter Set One and 54.4 under Parameter Set Two.

A secondary finding is that SNF-NP performs slightly better than SNF. This phenomenon is also observed in the next few experiments, the reason for which is unclear from our current theory. We hypothesize that under SNF, the jobs with large server needs sometimes do not fit into the remaining servers after they are preempted, leading to more idle servers. SNF-NP alleviates this issue by prohibiting preemption. However, this advantage of SNF-NP is moderate because the maximal server need is smaller than the slack capacity in our scaling regimes.

10.2. Experiment 2: Comparing with More Policies

One limitation of our policies, SNF and SNF-NP, is that they do not utilize the service rate information. Although we can show the order-wise optimality of SNF, our analysis is restricted to the scaling regimes where the service rates μ_i’s do not scale with n (i.e., the values of μ_i’s cannot be too large or small). When the values of μ_i’s vary significantly across different job types, it becomes questionable whether scheduling jobs solely based on server needs is the most effective way. For example, it could be more natural to prioritize jobs with smaller expected sizes, defined as $\frac{{\mathcal{\ell }}_i}{{\mu }_i}$ for jobs of type i.

Another aspect not considered in SNF and SNF-NP is the packing effect; SNF and SNF-NP do not aim to minimize the number of idle servers, $n-{\sum }_{i=1}^I{\mathcal{\ell }}_i{Z}_i$, when choosing the subset of jobs to serve, which might harm the mean waiting times. Furthermore, the approach of SNF and SNF-NP is somewhat contrary to common bin-packing algorithms. For instance, the best-fit algorithm prioritizes packing jobs with larger server needs rather than those with smaller server needs.

Given the motivation discussed, we compare the mean waiting times of SNF and SNF-NP with the following policies that utilize the expected sizes or consider the packing effect.

• Smallest size first (SSF). SSF prioritizes jobs with smaller values of $\frac{{\mathcal{\ell }}_i}{{\mu }_i}$ preemptively; it is the special case of the $c\mu /m$ rule in Zychlinski et al. (2023), with c being the all-one vector.

• Largest need first (LNF). LNF prioritizes jobs with larger server needs preemptively; it can be seen as following the best-fit algorithm to improve the packing.

• ServerFilling Gittins (SFG). SFG is a policy designed for multiserver job models with general service time distributions under the assumption that the number of servers n and the server needs are powers of two (Grosof et al. 2022b). Adapted to our setting, SFG schedules jobs based on both the expected sizes and the server needs, aiming to prioritize jobs with smaller expected sizes while optimizing packing.

We simulate the three policies described above along with SNF, SNF-NP, and FCFS. The parameters are similar to those used in Figure 3(a), with two modifications; we set ${\mathcal{\ell }}_2$ to be a power of two and increase the service rates for job types 2 and 3. The simulations are conducted under two sets of parameters with different service rates, namely Parameter Set Three and Parameter Set Four, as illustrated in Figure 4, (a) and (b), respectively.

From Figure 4, we can see that LNF is always worse than SNF and SNF-NP for $n\ge 2^8$. The performance of SSF depends on its priority order. When it gives type 3 jobs the lowest priority (e.g., when $n\ge 2^8$ in Figure 4(a)), it performs similarly to SNF and SNF-NP; otherwise, it performs similarly to LNF (for all n in Figure 4(b)). The performance of SFG lies between LNF and SSF. None of the three policies outperform SNF and SNF-NP for $n\ge 2^8$.

The comparison between LNF and SFG with SNF and SNF-NP indicates that considering the packing effect does not bring much benefit to the mean waiting time in the scaling regimes of the two examples in Figure 4. This outcome aligns with our expectations because the slack capacity here is larger than the maximal server needs, which ensures that any ${\mathcal{\ell }}_{\text{max}}$-work-conserving policy maintains enough busy servers to quickly reduce the workload even if the packing is not optimal.

On the other hand, the relationships between SSF, LNF, and SNF are more interesting; they suggest that in these examples, a policy is likely to have small mean waiting times when the maximal-need jobs (i.e., type 3 jobs) are served with the lowest priority, even if those jobs may have large service rates and the smallest expected sizes, as we can see from the data points in Figure 4(b) with $n\ge 2^8$. Fully understanding this phenomenon is outside the scope of our theory because we assume that the service rates do not scale with n, whereas in Figure 4(b), the service rates can be even larger than the server needs. In Online Appendix F, we give further details of Experiment 2 and provide some preliminary hypotheses about this phenomenon. There we make an interesting observation that in the setting of Experiment 2, the workload is considerably smaller under SNF than under LNF. Intuitively, SNF make the jobs with large service rates and large server needs stay longer in the system by giving them low priorities; those jobs will help keep more servers busy during periods when there are fewer jobs in the system.

10.3. Experiment 3: Mean Waiting Times Weighted by Server Needs

We have been focusing on the unweighted mean waiting time. However, sometimes we may want to give the jobs with larger server needs larger weights—for example, when each multiserver job consists of many unit-server-need components and one cares about the mean waiting time averaged over the components.

In Figure 5, we consider mean waiting times weighted by $\sqrt{{\mathcal{\ell }}_i}$ (Figure 5(a)) and weighted by ${\mathcal{\ell }}_i$ (Figure 5(b)) under Parameter Set Three (the same model parameters as Figure 4(a)). We can see that when the mean waiting times are weighted by $\sqrt{{\mathcal{\ell }}_i}$, for $n\ge 2^{10}$, SFG outperforms other policies followed by SNF-NP, SNF, SSF, and lastly, LNF and FCFS. When the mean waiting times are weighted by ${\mathcal{\ell }}_i$, LNF shows the best performance for all n followed by SFG and other policies.

These two simulations indicate that generally, when the weights of the mean waiting times depend on the server needs, other policies could outperform SNF and SNF-NP. Nevertheless, as shown in Figure 5(a), SNF and SNF-NP still perform reasonably well when the weights moderately depend on the server needs.

11. Conclusion and Future Work

In this paper, we establish order-wise sharp bounds on the mean waiting times of multiserver jobs under the FCFS policy and the SNF policy. We also prove a lower bound for the mean waiting times under all policies. These bounds imply the order-wise optimality of SNF and the strict suboptimality of FCFS.

Our analysis has two steps; we first prove the workload bounds stated in Lemmas 1 and 2 for a general class of policies using the drift method, along with some novel state-space concentration arguments. Then, we use the workload bounds to analyze the mean waiting times under FCFS and SNF, and we prove a lower bound for the mean waiting times under all policies.

Apart from the theoretical analysis, we also perform simulation experiments. We first demonstrate through simulations the performance improvement of SNF compared with FCFS in finite systems and the fact that SNF-NP, the nonpreemptive variant of SNF, has comparable performance with SNF. We also compare SNF and SNF-NP with other policies that utilize the service rate information or consider packing effects. We find that in our examples, when the number of servers is large, SNF and SNF-NP outperform those policies. Finally, we investigate the weighted mean waiting times under these policies, where jobs with larger server needs have larger weights. We find that in this case, other policies could outperform SNF and SNF-NP, but SNF and SNF-NP still perform reasonably well when the weights moderately depend on the server needs.

There are several interesting directions for future work.

1. Derive a tighter bound on the mean waiting time under SNF when Assumption 3 is violated.

2. Design a policy that achieves order-wise optimality in the many-server scaling limit without the maximal server need (assumption Assumption 2).

3. Understand why the performance of SNF-NP appears to be slightly better than SNF.

4. Extend the analysis to the settings where the service rates scale with n; this could help us understand why SNF and SNF-NP outperform other policies in Experiment 2.