key: cord-0144332-vrvg59yl
authors: Daw, Andrew; Fralix, Brian; Pender, Jamol
title: Non-Stationary Queues with Batch Arrivals
date: 2020-08-03
journal: nan
DOI: nan
sha: f82829db2cd70890b8b882d4bd9ff4f3e4fb4257
doc_id: 144332
cord_uid: vrvg59yl

Motivated by applications that involve setting proper staffing levels for multi-server queueing systems with batch arrivals, we present a thorough study of the queue-length process ${Q(t); t geq 0}$, departure process ${D(t); t geq 0}$, and the workload process ${W(t); t geq 0}$ associated with the M$_{t}^{B_{t}}$/G$_{t}$/$infty$ queueing system, where arrivals occur in batches, with the batch size distribution varying with time. Notably, we first show that both $Q(t)$ and $D(t)$ are equal in distribution to an infinite sum of independent, scaled Poisson random variables. When the batch size distribution has finite support, this sum becomes finite as well. We then derive the finite-dimensional distributions of both the queue-length process and the departure process, and we use these results to show that these finite-dimensional distributions converge weakly under a certain scaling regime, where the finite-dimensional distributions of the queue-length process converge weakly to a shot-noise process driven by a non-homogeneous Poisson process. Next, we derive an expression for the joint Laplace-Stieltjes transform of $W(t)$, $Q(t)$, and $D(t)$, and we show that these three random variables, under the same scaling regime, converge weakly, where the limit associated with the workload process corresponds to another Poisson-driven shot-noise process.

The M t /G t /∞ queueing system is arguably the most tractable time-varying queue studied in the literature, and is described as follows. Customers arrive to an area, consisting of infinitely many servers, in accordance to a non-homogeneous Poisson process with points {T n } n≥1 and arrival rate function λ : [0, ∞) → [0, ∞), and if a customer arrives to the system at time t, it brings with it a random amount of work having cumulative distribution function (CDF) F t for processing. We assume that the mth arrival to the system occurs at time T m , and it brings an amount of work S m for processing: hence, conditional on T m , the CDF of S m is F Tm . Finally, we let Λ : [0, ∞) → [0, ∞) denote the mean measure associated with the arrival process, where for each t ≥ 0,

For each real t ≥ 0, let Q(t) denote the number of customers present in the system at time t. Generally {Q(t); t ≥ 0} is not a Markov process, yet it is well-known that when Q(0) = 0 (or when the law of Q(0) is Poisson) the marginal distributions of {Q(t); t ≥ 0} are Poisson distributed: more particularly, assuming Q(0) = 0 with probability one, it can be shown that for each t > 0,

for each integer k ≥ 0, where F s (u) := 1 − F s (u)) for each u ≥ 0. Formula (1.1) can be proven in at least two different ways. One approach involves making use of a time-dependent thinning property of non-homogeneous Poisson processes: given our fixed t > 0, we say that if an arrival occurs at time s ∈ (0, t], we 'count' it with probability p t (s) := F s (t − s), independently of all other points in (0, t]. Then Q(t) is simply the number of counted points in (0, t], which is Poisson distributed with mean Another way to prove (1.1) is to simply note that {(T n , S n )} n≥1 correspond to the points of a spatial Poisson process on R 2 + , whose mean measure µ satisfies µ((a, b] × C) = b a C dF s (u)λ(s)ds for each a, b ∈ [0, ∞) satisfying a < b, and for each Borel measurable subset C of [0, ∞), where dF s (u) denotes Lebesgue-Stieltjes integration with respect to the CDF F s . Our primary objective is to illustrate how these ideas apply to non-stationary, infinite-server queueing systems with batch arrivals. Throughout we consider the infinite-server queueing system M Bt t /G t /∞, where batches of customers arrive in accordance to a non-homogeneous Poisson process {A(t); t ≥ 0} with rate function λ : [0, ∞) → [0, ∞). We denote the size (meaning number of customers) of the batch arriving at time t as B t , which is a random variable whose CDF depends on t, and we assume the amounts of work brought by customers within the batch has a joint distribution that depends on t. In general, we allow amounts of work within a given batch to be dependent, and no assumptions are placed on the distributions of work within a batch, but all batches are independent of each other. We associate with this infinite-server system the stochastic processes {Q(t); t ≥ 0} and {D(t); t ≥ 0}, where Q(t) denotes the number of customers present in the system at time t, and D(t) denotes the number of departures from the system that occur over the interval (0, t].

The research literature contains a large body of work addressing batch/bulk queueing systems: to the best of our knowledge, the first study featuring queues with batch arrivals is that of Miller Jr [21] . Since then, many other papers have been written that feature a study of queues with batch arrivals that operate under various different conditions: see for example Foster [10] , Shanbhag [28] , Brown and Ross [1] , Holman et al. [13] , Fakinos [9] , Chatterjee and Mukherjee [2] , Lucantoni [18] , Takagi and Takahashi [29] , Economou and Fakinos [7] , Masuyama and Takine [19] , Liu and Templeton [16] , Lee et al. [15] , Daw and Pender [4] . Later work has expanded the concept to a variety of related models, including priority queues and queues with server vacations. There are other papers in the literature that establish heavy traffic limit theorems for queues with batch arrivals: examples include Chiamsiri and Leonard [3] , Pang and Whitt [25, 26] . These papers show that under certain conditions, one can approximate a properly-scaled queue length process with a diffusion process-such as Brownian motion and Ornstein-Uhlenbeck processes-and also show that these approximations can be applied to even multi-server and non-Markovian queues.

Another recent application of batch queueing models is in the space of cloud-based data processing. In this case, the batches arriving to the system are collections of jobs submitted simultaneously. These jobs are then served by each being processed individually and returned. For more discussion, detailed models, and specific analysis for this setting, see works such as Lu et al. [17] , Pender and Phung-Duc [27] , Xie et al. [30] , Yekkehkhany et al. [31] and references therein. Another more relevant application is in infectious disease modeling such as COVID-19, see for example Kaplan [14] , Morozova et al. [23] , Palomo et al. [24] . In this setting, the results for patients who potentially have COVID-19 arrive in a large batch to be processed at a facility. Moreover, the data that we observe from COVID-19 is also of batch form as counts are made daily. Finally, an emerging application of batch queues is in context of autonomous vehicles moving in platoons (batches) down highways and roads, see for example Mirzaeian et al. [22] , Hampshire et al. [12] , Daw et al. [5] .

This paper contributes to the literature on infinite-server queues with batch arrivals in multiple ways. First, it was recently discovered in Daw and Pender [4] that the stationary distribution of a M B /G/∞ queueing system is equal in distribution to an infinite convolution of scaled Poisson random variables, where the means of these random variables are expressed in terms of the order statistics associated with the amounts of work found in an arriving batch. This discovery was made by first realizing that, for the case where all batch sizes are deterministic-more particularly, of size n for some integer n ≥ 1-the moment generating function (MGF) of the M B /M/∞ queue can be rewritten in a manner that allows for the MGF to be inverted analytically. Next, the authors extend this result to the M B /G/∞ case, again where each arriving batch consists of n customers, by reinterpreting the queueing system as a collection of n sub-queues. The authors then continue in Daw and Pender [4] by showing that the same results still apply when the batch sizes are random. Our first primary goal in this work is to illustrate how these observations carry over to the M Bt t /G t /∞ system, where we show that the marginal distributions of both the queue-length process, as well as the departure process of this queueing system are equal in distribution to an infinite sum of independent, scaled Poisson random variables. While this could also be carried out using the sub-queue construction featured in Daw and Pender [4] , we choose to instead use point process reasoning to find these marginal distributions. Next, we build on this point process approach even further, by illustrating how it can be used to calculate the finite-dimensional distributions of both the queue-length process and the departure process of the M Bt t /G t /∞ system. Right before submitting this paper, we discovered that in the works of Fakinos [9] , Economou and Fakinos [7] , the authors derive generating functions corresponding to the marginal distributions of both the queue-length process and the departure process for two different types of infinite-server queues with batch arrivals, and from their expressions (which also feature the use of order statistics within batches) one could theoretically deduce that these marginal distributions are scaled Poisson distributions, but to the best of our knowledge, Daw and Pender [4] is the first to observe this independent sum of scaled Poisson structure. This is an important observation, as once this fact is known, it provides a way of simulating the marginal distributions exactly. Not only that, the point process approach we use allows us to study, via a more elaborate thinning procedure, finitedimensional distributions with little difficulty, which to the best of our knowledge has not been done previously. We should also emphasize that our method of viewing the sum of scaled Poisson random variables through the lens of order statistics via our sub-queue perspective provides an explicit construction of the queue length and departure processes that gives a visual construction that is easy to understand.

Second, we show that for the M Bt t /G t /∞ queue, it is actually possible to derive the joint Laplace-Stieltjes transform (LST) of W (t), Q(t), and D(t), which corresponds to the workload at time t, the number of customers present in the system at time t, and the number of departures from the system in the interval (0, t], respectively. While some of the works cited above do feature a study of various auto-covariance functions associated with these processes, in particular the covariance of Q(t) and D(t), to the best of our knowledge, no one has yet to study the joint LST of W (t), Q(t), and D(t).

Finally, our third major contribution involves using the above-mentioned results to establish various scaling-limit theorems. Recent work by De Graaf et al. [6] , Daw and Pender [4] has shown that a novel batch scaling of M B /G/∞ queues converge to a shot noise process in steady state. Our work generalizes the work by De Graaf et al. [6] , Daw and Pender [4] in several ways. First, we consider infinite server queues featuring non-stationary (and general) service time distributions, non-stationary batch sizes, and non-stationary arrival rates. Second, we prove that the batch scaling extends to this general setting and we extend the results of De Graaf et al. [6] by establishing convergence of the finite dimensional distributions directly. In order to establish weak convergence of the finite-dimensional distributions of the rescaled queue-length processes, the authors of De Graaf et al. [6] use a sophisticated "convergence of the generator" approach, which is arguably less direct than the approach we provide.

This paper is organized as follows. In Section 2, we use the order statistics associated with the amounts of work found in each arriving batch to show that, for each t ≥ 0, both Q(t) and D(t) can be expressed as an infinite sum of independent, scaled Poisson random variables. We then continue by applying the same ideas to study the finite-dimensional distributions of both {Q(t); t ≥ 0} and {D(t); t ≥ 0}, through the calculation of joint LST and auto-covariance functions, and we use these expressions to establish a weak-convergence result that builds on a scaling-limit result found in De Graaf et al. [6] . In Section 3, we analyze the joint LST of W (t), Q(t), and D(t), and we use this expression to again derive a scaling limit result that builds on the scaling-limit result of De Graaf et al. [6] in a different way.

The thinning approach involves using the order statistics within each arriving batch to perform a thinning procedure to the non-homogeneous Poisson arrival process governing batch arrivals. The specifics of the approach will depend on the random variable we wish to study, so for now we focus on explaining how the procedure can be used to study the law of Q(t), when Q(0) = 0.

For each integer n ≥ 1, and each j ∈ {1, 2, . . . , n}, let S j,n (s) denote the amount of work brought by the j th customer contained in the batch of size n that arrives at times s, and let S j:n (s) denote the j th smallest amount of work found in the same batch. Then for each integer j ∈ {1, 2, . . . , n}, and each real a, b satisfying 0 ≤ a < b ≤ t, the random variable

counts the number of batches arriving in the interval (a, b] that are of size n, and are such that precisely j − 1 customers within the batch have departed by time t, and precisely n − (j − 1) customers from this batch are still present in the system at time t: we assign each such batch with the label (j; n). Using the fact that the both the size of each batch, as well as the amounts of work present in each batch are independent of all other batches, it follows from a standard thinning procedure of non-homogeneous Poisson processes that the random variables {Y j;n (a, b]} 1≤j≤n are independent, Poisson random variables, where the mean of Y j;n (a, b] is given by

where P s is the probability measure associated with the batch of customers that arrive at time s. Throughout we follow the convention that S 0;n (s) = 0 with probability one, and S n+1;n (s) = ∞ with probability one. t /G/∞ queueing system using Poisson random measures in the style of Eick et al. [8] , which allows us to represent the thinning perspective.

As a visualization of this idea, let us adapt the elegant Poisson random measure perspective shown in Figure 1 of Eick et al. [8] . In this diagram, the solid vertical lines mark the times of arrivals in the Poisson process. The dots along these lines then denote the lengths of the service durations within each arriving batch. Of course, by comparison to the M t /G/∞ queue considered by Eick et al. [8] , the batch arrivals mean that there are multiple service durations for each arrival epoch in the Poisson process. Because a customer is still in the system if her arrival time plus her service duration is greater than the current time, the total queue length is the number of points above the 45 • line. In this way, we can classify all of the arrival epochs up to t by the number of jobs within a batch that remain in the system at time t, meaning the number of points above this line. This classification of the arrival times yields our thinning of the Poisson process.

This style of reasoning leads us to our first result, which shows how to express Q(t) as a sum of independent, scaled Poisson random variables. The same can also be said of D(t), which is defined as the number of departures from the system in the interval (0, t].

Theorem 2.1. Assume Q(0) = 0 with probability one. For each t > 0,

Proof. If a batch arriving in the interval (0, t] is assigned label (j; n), precisely (n − j + 1) customers within that batch are still present at time t, and precisely j − 1 of those customers have departed by time t. Hence, the number of customers present in the system at time t from a batch with label (j; n) is (n − j + 1)Y j;n (0, t], and the number of departures over (0, t] of customers from a batch with label (j; n) is (j − 1)Y j;n (0, t]; adding over all possible labels completes the proof.

Remark Readers may observe that we include the term (0)Y n+1;n (0, t], as well as the term (0)Y 1;n (0, t] in Q(t) and D(t), respectively, which seems unnecessary since both terms are clearly zero with probability one. However, following this convention will make it easier later to express various joint Laplace-Stieltjes transforms associated with both {Q(t); t ≥ 0} and {D(t); t ≥ 0}. Now that we have shown the Poisson decomposition of the queue length and departure processes, we can use the representation to analyze the covariance between the two processes. The next result provides the covariance between Q(t) and D(t).

Proposition 2.2. For each t > 0, the covariance of Q(t) and D(t) is given by

Proof. The proof of this result exploits properties of Poisson processes and the decomposition of the queue length and departure process given in Theorem 2.1.

Proposition 2.2 shows that in general, Q(t) and D(t) are positively correlated, but when all batches are of size one with probability one, the result shows that

which is a well-known result that is addressed in e.g. Eick et al. [8] . The same Poisson random measure visualizations can reveal this dependence as well and demonstrate the difference between the batch and solitary arrival settings, as shown in Figure 2 . Let time t and offset u be fixed. Following the decomposition used in Eick et al. [8] , let us introduce the quantities A(t, u), B(t, u), and C(t, u) defined such that A(t, u) is the number of entities that arrive by time t and depart in the interval [t, t + u], B(t, u) is the number of entities arriving in [t, t + u] that remain in the system at time t + u, and C(t, u) is the number of entities arriving by time t that remain in the system at time t + u. Then, by definition we have that Q (n) (t) = A(t, u) + C(t, u) and Q (n) (t + u) = B(t, u) + C(t, u). By the independent increments of the Poisson process, we can note that B(t, u) is independent from A(t, u) and C(t, u). However, unlike the M t /G/∞ model studied in Eick et al. [8] , A(t, u) is not independent from C(t, u). This is a consequence of the batch arrivals, as there is dependency between the ordered service times within one batch. Using these definitions, we can express the auto-covariance in terms of these regions as

This thinning technique can also be used to derive the joint finite-dimensional distributions of both {Q(t); t ≥ 0} and {D(t); t ≥ 0}, but doing so requires a more elaborate thinning procedure. Given a collection of real numbers {t ℓ } m ℓ=1 satisfying 0 < t 1 < t 2 < . . . < t m and an integer n ≥ 1, we define the random variable Y j k ,j k+1 ,...,jm;n as

where Y j k ,j k+1 ,...,jm;n (a, b] should be interpreted as the number of batches of size n that arrive in the interval (a, b] satisfying the property that for each ℓ ∈ {k, k + 1, . . . , m}, exactly j ℓ − 1 customers from the batch have departed from the system at time t ℓ (meaning also that exactly n − j ℓ + 1 customers from the batch are still present in the system at time t ℓ ). Using well-known thinning properties of non-homogeneous Poisson processes, we can say that Y j k ,j k+1 ,...,jm (a, b] is a Poisson random variable that satisfies

These random variables contribute value to each Q(t ℓ ), as well as to each D(t ℓ ) value, for 1 ≤ ℓ ≤ m.

Our next result provides an expression for the joint LST of the finite-dimensional distributions of both {Q(t); t ≥ 0} and {D(t); t ≥ 0}, as well as the auto-covariance functions of both {Q(t); t ≥ 0} and {D(t); t ≥ 0}. is as follows: for α :

where

Furthermore, the auto-covariance functions of {Q(t); t ≥ 0} and {D(t); t ≥ 0} are as follows: for each t 1 , t 2 satisfying 0 < t 1 < t 2 ,

Proof. We begin by deriving both (2.4) and (2.5) . Considering first the random vector (Q(t 1 ), Q(t 2 )), we see that for (α 1 , α 2 ) ∈ R 2 + ,

These representations for Q(t), Q(t + u), D(t), and D(t + u) can be used to derive the autocovariance functions. Indeed,

{S j ℓ −1:n (s) ≤ t ℓ − s, S j ℓ :n (s) > t 1 − s} λ(s)ds and a similar argument can be used to establish (2.5) .

It remains to prove (2.2). Given any collection of real numbers 0 = t 0 < t 1 < t 2 < t 3 < · · · < t m−1 < t m , we have that for α = (α 1 , α 2 

This new representation shows that m k=1 (α k Q(t k ) + β k Q(t k )) can be expressed as a finite sum of independent, scaled Poisson random variables. Further exploitation of this observation gives Even though Theorem 2.3 holds with minimal assumptions placed on the arrival rates of batches and the services found within each batch, it is fairly clear that the above auto-covariance functions, as well as the LST of the finite-dimensional distributions, are extremely complicated, but this is the price we pay for relaxing our original assumptions as much as possible. If we consider specific settings, we are able to derive simplified expressions. For example, in the case of stationary exponential service we can cleanly relate the auto-covariance and the variance.

Proposition 2.4. If the service is exponentially distributed at rate µ > 0, the auto-covariance of the queue length is such that

for t, δ ≥ 0.

Proof. Since the queue length at time t + δ can be written as the queue length at t plus the number of arrivals in [t, t + δ) and less the number of departures in [t, t + δ), i.e.

we can decompose the auto-covariance accordingly. That is, by the definition of covariance we have that

Since both the future of the arrival process and the sequence of batch sizes are independent from the history of queue, these terms cancel one another. With the linearity of expectation, this then simplifies to

Given the queue length at time t, the number of departures on the interval [t, t + δ) can be written as a sum over all services that were completed. Using the memoryless-ness of exponential service, this means that

where S j ∼ Exp(µ) are mutually independent and also independent of Q(t). Through conditional expectation, we can also observe that

Using this observation and the analogous result for the mean number of departures, we can further simplify the auto-covariance to

which completes the proof.

In addition to the exponential case, we can also explicitly analyze the auto-covariance in the case where the service distribution is deterministic. In this case, we leverage the fact that when the service is deterministic, we can write the queue length as the difference of the arrival process at the current time and the arrival process at the current time but delayed by the constant service time.

Proposition 2.5. If the service is deterministic with constant ∆ > 0 and the batch size is fixed as size n, then auto-covariance of the queue length is such that Cov[Q (n) )(t), Q (n) )(t + δ)] = n 2 t t−δ+∆ λ(s)ds.

if ∆ > δ. Otherwise the auto-covariance is zero.

Proof. This result follows immediately after one makes the observation that in the deterministic service setting that Q (n) )(t) = n · (A(t) − A(t − ∆)) and

Thus, the covariance can be written as

we obtain the final result.

In our next result, we find that under the assumptions where, for a batch that arrives at time s, all services within that batch are i.i.d. with CDF F s , as well as independent of the batch size B s , the finite-dimensional distributions of both {Q(t); t ≥ 0} and {D(t); t ≥ 0} simplify considerably. is as follows: for α := (α 1 , α 2 , . . . , α m ) ∈ R m + , β := (β 1 , β 2 , . . . , β m ) ∈ R m + , we have

where for each k ∈ {1, 2, . . . , m}, each integer b ≥ 1, and each s ∈ [t k−1 , t k ), we have

Proof. Our objective now is to simplify the Laplace-Stieltjes transform found in (2.3). Using standard properties of order statistics associated with i.i.d. random variables, we find that for each s ∈ (t k−1 , t k ],

where F s (s, t] := P s (s < S 1 (s) ≤ t), and F s (t) = P s (S 1 (s) > t).

Plugging this probability into (2.2), after plugging (2.3) into (2.2), and combining all exponential terms yields, within the exponent, the summation

but this sum is simply

k,n,α,β (s).

Hence,

which proves the claim.

We conclude this section with a simple convergence result that builds on the recent work of [6] . Suppose {Z(t); t ≥ 0} is a fixed Poisson-driven shot-noise process defined as

for each t ≥ 0, with random (and nonnegative) jump size B s when {A(t); t ≥ 0} has a point at time s, and {F s } s≥0 is again a collection of cumulative distribution functions. Again, we assume {A(t); t ≥ 0} is a nonhomogeneous Poisson process, having rate function λ and mean measure Λ.

From this shot-noise process, we construct a sequence of infinite-server queueing systems {Q n (t); t ≥ 0} where for each fixed integer n ≥ 1, the system associated with {Q n (t); t ≥ 0} is such that the size of a batch arriving at time s is B (n) s := ⌈nB s ⌉, but all other random elements of the nth infinite-server queue are equal in distribution to the corresponding random elements from our fixed M Bt t /G t /∞ queue. Letting also {D n (t); t ≥ 0} represent the departure process associated with the nth system, we define for each t ≥ 0 the quantities

Theorem 2.7. For each integer m ≥ 1, and each collection of real numbers 0 < t 1 < t 2 < . . . < t m , we get Proof. Observe that for each α ≥ 0, β ≥ 0,

k,b,α/n,β/n (s)Ps(B (n) s =b)λ(s)ds meaning that establishing this limit amounts to finding

for each k ∈ {1, 2, . . . , m}. Fix an integer k, along with a real number s ∈ [t k−1 , t k ]. Then

Fixing now y > 0, define b n (y) := inf{b ≥ 1 : b/n ≥ y}. 

In this thinning decomposition of the queue length into a sum of scaled Poisson random variables, we are taking a retrospective approach. Given the current time, we classify all preceding arrival epochs based on their relation to the fixed present moment. However, one can actually reach the same Poisson sum decompositions through considering the queue as a collection of evolving and inter-related sub-systems. For the sake of example consider a fixed batch size, which may also be thought of as thinning the Poisson process by conditioning that the randomly drawn batch sizes are of a particular cardinality. Then, rather than retro-actively classifying the arrival epochs, we can classify the jobs within each batch based on the relative length, again leading us to order statistics. If the fixed batch size is n, one can suppose that there are n sub-queues, each with infinitely many servers. That is, let Q 1 , ..., Q n be infinite server queues. Suppose that the customer with the earliest service completion is sent to Q 1 , the customer with the second earliest is sent to Q 2 , and so on. When viewing each sub-system on its own, we see that Q j is an infinite server queue with single arrivals according to a Poisson process with rate λ(t) and service distribution matching that of S (j) , the j th order statistics of G. Thus, we can see that

through the literature for M/G/∞ queues, such as in Eick et al. [8] .

Now, we can also note that there is inherent correlation between these sub-queues, but we can also explicitly identify it. Because each order statistic can be written as a telescoping sum of the lower order statistics, i.e.

S k:n − S k−1:n , the service durations in a given sub-system are built out of pieces that are repeated in every higher indexed system. The first piece S 1:n is repeated n times, the second piece S 2:n is repeated n − 1 times, and so on. In each sub-queue we can use the thinning of Poisson random variables and the union of disjoint events to write the distribution of Q j as a sum of Poisson random variables, as given by

Poisson t 0 λ(s)P s (S k−1:n ≤ t − s, S k:n > t − s) ds .

Then, by using our observation of the repeated pieces between the sub-systems, we can re-assemble the full queue as the same decomposition that we have identified through the thinning perspective. This alternate splitting perspective may be useful in understanding various sub-systems of the queue, as originally discussed in Daw and Pender [4] . As a summary of this splitting perspective and the role of the repeated pieces, in Figure 3 we visualize how the different sub-queues are formed and how this yields the sum of scaled Poisson's identity.

We can also contrast the thinning and splitting ideas through the Poisson random measure diagrams. In Figure 4 , we distinguish the two concepts as vertical and horizontal classifications, respectively. As we discussed for Figure 1 , the thinning perspective classifies the vertical arrival lines based on their relation to the current time, specifically the number of arrivals left at time t. Hence, in the left hand side of Figure 4 the vertical lines are color coded for the number of points above the 45 • line. Because the splitting approach assigns jobs to sub-queues based on their relative length, the right hand side figure colors the points instead of the lines, grouping times together across the batches horizontally.

The splitting decomposition and its visualization in Figure 3 also provide a natural context for considering the workload process. Because there are infinitely many servers, all jobs begin being processing immediately. While a single server queue will process the workload at a unit rate whenever there is at least one customer in the system, an infinite server queue's workload decreases at rate equal to the number of jobs in the system. Thus, much like how the sub-systems feature repeated pieces of the same service lengths, the work brought in by a batch of size n will be processed at rate n through the first order statistic, then at rate n − 1 until the second, and so on.

Parallel Queues Independent Scaled Poissons In the following section, we will use this observation to further analyze and discuss the M Bt t /G/∞ workload process in detail.

Results on the workload process of time-varying infinite-server queueing systems with batch arrivals, even when each batch is of size one, appear to be scarce. The most relevant reference we found that even remotely addresses the workload process of time-varying infinite-server queues with Poisson arrivals is Goldberg and Whitt [11] , which is primarily concerned with the study of the last departure time from a M t /G/∞ queueing system, when the arrival process stops at some fixed, deterministic time τ . While Goldberg and Whitt [11] do not study the workload process in itself, Theorem 2.1 of Goldberg and Whitt [11] can be used to derive the LST of the workload process of the M t /G/∞ queue, as this result provides the conditional joint distribution, given Q(t) = n, of the remaining service times of the n customers present in the system at time t.

The next proposition is a slight generalization of Theorem 2.1 of Goldberg and Whitt [11] , in that it applies to the M t /G t /∞ system, and it can be proven in precisely the same manner, which involves conditioning on the order statistics associated with the thinned Poisson process associated with customers that are still present in the system at time t, then simplifying: we omit the details. 

Once Proposition 3.1 is known, it can be used to calculate the LST of W (t) for the M t /G t /∞ queue. Again, we omit the proof as it follows from conditioning on Q(t), then applying Proposition 3.1. Figure 4 : A comparison of the "thinning" and "splitting" perspectives. On the left, the "thinning" approach distinguishes the arrival epochs based on the number of entities from the batch that remain at time t. On the right, the "splitting" approach groups each entity by its order within its batch.

Proposition 3.2. The LST of W (t) is as follows: for each α ≥ 0,

where φ t is the LST associated with the CDF H t .

As an interesting aside, it is worth noting that Theorem 2.2 of Goldberg and Whitt [11] (namely, Identity (2.5) of Goldberg and Whitt [11] ) can be derived with our thinning approach from Section 2, without applying Proposition 3.1. Proof. Fix x ≥ 0, and let Y 0 (0, t] denote the number of jobs that arrive in the interval (0, t] that are still present in the system at time t + x. This random variable is a Poisson random variable with mean

proving the claim.

Just as we have used order statistics to understand the queue length and departure processes, we can also use these quantities to contextualize the M Bt t /G t /∞ workload process. In Figure 5 , we illustrate how the workload process moves over time. If there are no arrivals that occur, then the workload decreases linearly at rate 3 until the first customer leaves, then it proceeds down at rate S (1) +S (2) +S (3) S (2) +S (3) S (3) S (3) S (1) S (2) Figure 5 : A hypothetical sample path of the workload process, where the change in the rate of decrease can be seen as each order statistic is passed.

2 until the second customer leaves, and then finally it proceeds downward at a unit rate until the final customer leaves the system. From this perspective, the workload process can be described as a filtered point process, see e.g. Michel [20] . We present an alternative approach towards studying, for the M Bt t /G t /∞ queueing system, the joint LST of W (t), Q(t), and D(t), where our approach also makes use of conditioning on the number of arrivals that occur in the interval (0, t], but instead of making use of a result analogous to Proposition 3.1, we keep track of all arrivals in (0, t], then use indicator functions to describe W (t), Q(t), and D(t) once we know when all arrivals occur in (0, t]. Doing the calculations in this way allow us to more easily work with random batch sizes. where E s denotes expectation, conditional on having a batch arrival at time s. where E s 1 ,...,sm represents conditional expectation, given batches arrive at times s 1 , s 2 , . . . , s m . Furthermore, since batches are independent, (3.11) and this proves that the integrand of the multiple integral found in (3.10) is a symmetric function on [0, t] n . Hence, (3.11) simplifies to

and after plugging (3.12) into (3.9) and simplifying, we get

proving Theorem 3.4.

In light of Theorem 3.4, it is not difficult to see that the joint LST of W (t), Q(t), and D(t) simplifies significantly under the additional assumption that within a batch arriving at time s, the amounts of work are all i.i.d. with CDF F s . Corollary 3.5. Suppose that when a batch arrives at time s, each customer within that batch brings a generally distributed amount of work with CDF F s , independently of everyone else. Next, for each s, t ≥ 0, let X s be a random variable whose CDF is F s , and define the Laplace-Stieltjes transform

Then for α ≥ 0, β ≥ 0, and γ ≥ 0, Even though each LST φ s,t (α) typically does not simplify much further, it is noteworthy to realize that φ s,t (α) can be expressed reasonably well for the special case where F s is the CDF of a phase-type random variable. In particular, when F s is the CDF of an exponentially distributed random variable with rate µ s , we get φ s,t (α) = µ s µ s + α .

In the next corollary, we use this simple fact to show that the joint LST of Q(t) and W (t) simplifies considerably when all amounts of work are exponentially distributed.

Corollary 3.6. Suppose that when a batch arrives at time s, each customer within that batch brings an exponentially distributed amount of work with rate µ s , independently of everyone else (i.e. the amounts of work arriving at time s are i.i.d. exponentially distributed with rate µ s ). Then for α ≥ 0, β ≥ 0, and γ ≥ 0,

Bs λ(s)ds .

We close this section by reconsidering the scaling-limit regime examined in De Graaf et al.. We consider a shot-noise process fed by a non-homogeneous Poisson process {A(t); t ≥ 0} with rate function λ and mean measure Λ, and we assume that if {A(t); t ≥ 0} has a point at location s, the jump size of the shot-noise process at that point is B s , and the decay pattern associated with that point is F s (t − s), where F s is a CDF.

Just as was done in Section 2, we consider a sequence of queueing systems indexed by m, where {W m (t); t ≥ 0}, {Q m (t); t ≥ 0}, and {D m (t); t ≥ 0} correspond to the workload, queueing process and departure process associated with the mth M B (m) t t /G t /∞ queue, which is fed by the Poisson arrival process {A(t); t ≥ 0} associated with the above-mentioned shot-noise process. We further assume that if a batch arrives at time s, the amounts of work found within that batch are i.i.d. with CDF F s , and the number of jobs/customers found in the batch is B

Next, for each integer m ≥ 1 we define the rescaled processes Our next result, Theorem 3.7, shows that the random vectors {(W m (t), Q m (t), D m (t))} m≥1 converge weakly as m → ∞.

s (t−s)+βF s(t−s)+γFs(t−s)] λ(s)ds .

In particular, What is especially notable about this result is that even though both the marginal distributions of both {Q m (t); t ≥ 0} and {W m (t); t ≥ 0} converge weakly to the marginal distributions of shot-noise processes, the two shot-noise processes have different decay patters: the decay pattern associated with the scaling-limit of the queue-lengths is F s (t), while the decay pattern associated with the scaling-limit of the workload processes is E s [S 1 (s)]F e s (t). We can also see that, when F s corresponds to the CDF of an exponential random variable with rate µ (independent of s), the two scaling limits actually coincide if we further multiply W m (t) by µ. 

In this paper, we generalize many results about the infinite server queue to the case of non-stationary batch sizes. Using two different methods, the thinning approach and the splitting approach, we prove a novel decomposition of the M Bt t /G t /∞ queue length distribution and the departure process in terms of sums of scaled Poisson random variables. In fact our work generalizes Eick et al. [8] by giving a Poisson representation for batch queues. We show that the number of Poisson random variables needed to describe the queue length process is intimately connected to the support of the batch distribution. Thus, the infinite sum representation collapses into a finite sum when the batch distribution is supported on a finite set. We also discover that the independence property of the queue length and departure process is only true in the case when the jumps are of size one, thus showing that the independence result of Eick et al. [8] is a special case and not true for batch arrivals greater than or equal to two.

We also prove a batch scaling result, which shows that as the batch size gets larger, the queue length process converges to a non-stationary shot noise process. We show this result for the finite dimensional distributions and prove it by leveraging our Poisson decomposition. This method of proof not only allows us to generalize previous results, but also prove it in an elegant and insightful way. As a side result, we derive the joint Laplace-Stieltjes transform of the queue length and workload processes, which is the first of its kind in the queueing literature to our knowledge. Finally, from a simulation perspective, our work presents a new way of simulating from the M Bt t /G t /∞ queue length distribution by just simulating and adding Poisson random variables and without keeping track of the queue length process. We hope to apply our work to the setting of multiserver queues with batch arrivals and exploring how we can use our Poisson representations for staffing multiserver queues to achieve stable performance.

Another interesting direction of future research is to extend this method of decomposition to batch arrival queues with arrival processes that are not Poisson. While the thinning technique we've used is inherently based on the Poisson process, the sub-queue construction visualized in Figure 3 may be readily applied to other contexts. In particular, this notion of distilling a batch arrival infinite server queue down into a collection of solitary arrival queues may be quite useful in many different settings. The remaining question for future research is, however, identifying how these sub-systems are distributed. In this case, we have been able to capitalize on the literature for M t /G/∞ queues to recognize that the sub-systems are themselves Poisson distributed, but this is of course dependent on the Poisson arrival process. For arrival processes that are not Poisson processes but may be closely related, such as e.g. Cox processes, it may be promising to use this near-Poisson-ness to uncover the resulting distributions of these solitary arrival queues, and thus gain insight into the distribution of the batch arrival queues overall.

Some results for infinite server poisson queues

On the non-homogeneous service system

A diffusion approximation for bulk queues

On the distributions of infinite server queues with batch arrivals

Can teleoperations efficiently support autonomous vehicles? a critcal staffing question

Shot-noise fluid queues and infinite-server systems with batch arrivals. Performance evaluation

The infinite server queue with arrivals generated by a non-homogeneous compound poisson process and heterogeneous customers

The infinite server queue with arrivals generated by a non-homogeneous compound poisson process

Batched queuing processes

The last departure time from an M t /G/∞ queue with a terminating arrival process

Beyond safety drivers: Applying air traffic control principles to support the deployment of driverless vehicles

On the service system

Covid-19 scratch models to support local decisions

Batch arrival queue with N-policy and single vacation

Autocorrelations in infinite server batch arrival queues

Joinidle-queue: A novel load balancing algorithm for dynamically scalable web services. Performance Evaluation

New results on the single server queue with a batch Markovian arrival process

Analysis of an infinite-server queue with batch Markovian arrival streams

A point process approach to filtered processes

A contribution to the theory of bulk queues

A queueing model and analysis for autonomous vehicles on highways

A model for COVID-19 transmission in connecticut. medRxiv

Flattening the curve: Insights from queueing theory

Two-parameter heavy-traffic limits for infinite-server queues

Infinite-server queues with batch arrivals and dependent service times

/delayoff-setup queues with nonstationary arrivals

On infinite server queues with batch arrivals

Priority queues with batch Poisson arrivals

Pandas: Robust locality-aware scheduling with stochastic delay optimality

GB-PANDAS:: Throughput and heavy-traffic optimality analysis for affinity scheduling

This paper was partially composed when Andrew Daw was a doctoral student at Cornell, at which time he was supported by the National Science Foundation through a Graduate Research Fellowship under grant DGE-1650441.