Open networks of infinite server queues with non-homogeneous multivariate batch Poisson arrivals

In this paper, we consider the occupancy distribution for an open network of infinite server queues with multivariate batch arrivals following a non-homogeneous Poisson process, and general service time distributions. We derive a probability generating function for the transient occupancy distribution of the network and prove that it is necessary and sufficient for ergodicity that the expected occupancy time for each batch be finite. Further, we recover recurrence relations for the transient probability mass function formulated in terms of a distribution obtained by compounding the batch size with a multinomial distribution.


Introduction
Infinite server queues with Poisson batch arrivals -denoted M X t /G/∞ queues in Kendall's notation -have elicited significant research attention, from the canonical work of Shanbhag [22], Reynolds [20], and Brown and Ross [2] in the 1960s, to studies by Chatterjee and Mukherjee [3] and Liu and Templeton [13] and the recent treatise of Daw and Pender [6].The independence between customers within the queue -termed the 'non-interacting property' by [13] -allows for an analytic characterisation of the numbers of busy servers and departures, typically in the form of a probability generating function (PGF) [2,3,6].
A natural generalisation is an open network of infinite server queues, with multivariate batch arrivals following a non-homogeneous Poisson process and general service time distributions.With customers arriving one at a time, networks of M t /G/∞ queues have been studied by Harrison and Lemoine [7], Keilson and Servi [9], and Massey and Whitt [14], amongst others.In contrast, networks with multivariate batch arrivals have received comparatively little attention.Models have been proposed for specific applications such as private-line telecommunication services [15] and malarial relapses [16], while stochastic fluid networks have been shown to arise as scaling limits of appropriate batch-arrival infinite server queueing networks [10].To the best of our knowledge, results for more general networks of •/G/∞ queues with multivariate batch arrivals are yet to appear in the literature.
Here, we analyse the occupancy distribution in a network of infinite server queues, with multivariate batches arriving according to a non-homogeneous Poisson process.Our model assumptions are set out in Section 2. We derive the PGF for the transient queue length distribution in Theorem 3.1, extending the construction we adopted in Mehra et al. [16]; similar arguments have been presented previously by [2,7,3] and others.
While infinite server queues with single arrivals are necessarily stable when the expected service time is finite, the queue length Markov chain can be transient or null-recurrent when batch arrivals are permitted.Stability conditions for the M X /G/∞ queue, accommodating heterogenous customers within each batch, have been characterised by Cong [5].When service times are exponential, Yajima, Phung-Duc, and Masuyama [30] have shown that it is necessary and sufficient for an infinite server queue with a batch Markovian arrival process to be stable that the batch size distribution has a bounded logarithmic moment.For networks with multivariate batch arrivals, we extend these results in Section 4 to show that a necessary and sufficient condition for ergodicity is that the expected occupancy time for each batch be finite (Theorem 4).When the customer occupancy time distributions have exponential tails, this is equivalent to the logarithmic moment condition of [5,30] (Corollaries 4.1.2and 4.1.3).
We recover recurrence relations for the probability mass function (PMF) of the time-dependent queue length distribution, formulated with respect to the PMF of the batch size compounded with a multinomial distribution in Theorem 5.1.The utility of these formulae is constrained by the underlying compound distribution.In the context of the M X /G/∞ queue, Willmot and Drekic [27,28,29] found that for distributions within 'Sundt and Jewell's family' -encompassing the binomial, Poisson, logarithmic, geometric and negative binomial distributions, and zero-inflated analogues thereof [26] -the recurrence relations may be computationally tractable.
Here, we find that univariate batch sizes within 'Sundt and Jewell's family' [26] may likewise yield tractable formulae when the initial queue allocation for each customer is independent and identically distributed (i.i.d.) (Section 5.1.1).

Model structure
Consider an open network of J infinite server queues, indexed by j = 1, . . ., J, such that • Arrivals occur at points T 1 , T 2 , . . . of a non-homogeneous Poisson process with rate λ(t).
• At time point T i , for j = 1, . . ., J, Σ ij customers arrive at queue j.The vectors Σ i = (Σ i1 , . . ., Σ iJ ) are independent vector-valued random variables with non-negative integer components selected from a multivariate distribution with probability generating function defined on a subset D ⊆ C J that contains the J-dimensional unit polydisc with |z j | ≤ 1.
Note that the distribution of Σ i does not depend on i.We use the generic random variable S to denote a random variable with this distribution.
• The th customer that enters queue j at time point T i is marked with a stochastic process X j i, (t) where X j i, (t) = k if the customer will be in queue k at time T i + t.For each i and = 1, . . ., Σ ij , the X j i, (t) are selected independently with a law such that P (X j i, (t) = k) = q j k (t).
Under this model, the routes of customers through the network (that is, the sequence of visited nodes and associated service times) are mutually independent, conditional on the node of entry.
This 'non-interacting property' [13] enables explicit analysis.Here, we focus on the number of customers in the network at time t, denoted by the random vector N(t) = (N 1 (t), . . ., N J (t)).
For notational convenience, we introduce the random vector C i (t) = (C i,1 (t), . . ., C i,J (t)), with giving the number of customers from the ith batch that are in queue k at time T i + t.It follows from our assumptions above that the distribution of C i (t) is independent of i and we denote by C(t) a generic random variable with this distribution.

The transient queue length distibution
We first characterise the time-dependent PGF for the state of the network, generalising the construction of Mehra et al. [16] (which was used to analyse an open network of infinite server queues tailored to within-host superinfection and hypnozoite dynamics for Plasmodium vivax malaria).
Theorem 3.1.(Multivariate PGF) For z = (z 1 , . . ., z J ) ∈ D Proof.We use similar reasoning to Mehra et al. [16].In short, we condition on the multivariate batch size, and the sequence of arrival times.For related systems, similar proofs have been previously presented by others, including Brown and Ross [2], Harrison and Lemoine [7], and Chatterjee and Mukherjee [3].
We begin by deriving the PGF for the generic random vector C(t) described in Equation (2).
Since X j i, (t) are independent random variables taking values in {1, . . ., J}, the number of customers originating in queue j that are in queue k at time T i +t has a multinomial distribution with parameters q j k (t).Therefore, conditional on a multivariate batch of size .
By the law of total expectation, it follows that Under the assumption that arrivals follow a non-homogeneous Poisson process, the number of where Following Harrison and Lemoine [7], conditional on the event {M (t) = m}, the distribution of T 1 , . . ., T m is the same as the distribution of the order statistics of m i.i.d.variables with density Using the law of total expectation, we thus deduce Substituting Equation ( 5) into Equation ( 6) yields expression (4).

Necessary and sufficient conditions for ergodicity
Here, we characterise the stability, that is, ergodicity, of the network, assuming a homogeneous arrival rate λ(t) = λ.For M/G/∞ queues, a necessary and sufficient condition for ergodicity is that the occupation time distribution has finite expectation.The analogous constraint in this setting is that at least one customer from each batch will be present somewhere in the network for a time W that has finite expectation.
A customer who arrives to queue j at time T has left the network the network before time T + t with probability The expected time that at least one customer from a batch is present in the network can then be written where we have used Tonelli's theorem to swap the order of summation and integration; and noted that the resultant integrand can be written as a sum of two absolutely convergent series.
This expression is central to the ergodicity of the network.By the dominated convergence theorem, for all z ∈ Ω, Observe that the transient PGF φ t (z) for N(t) defines a family of functions {φ t } t≥0 that is both analytic and bounded on Ω, with |φ t (z)| ≤ 1 for all t ≥ 0. By Montel's theorem, there exists a subsequence φ t i which converges uniformly on all compact subsets K ⊂ Ω (Theorem 1.7.3 of [11]).It follows from a theorem of Weirstrass (Theorem 1.7.1 of [11]) that lim t→∞ φ t is analytic on Ω, with the uniform convergence of partial derivatives on all compact subsets K ⊂ Ω.In particular, for all v ∈ Z J , we note that It is therefore sufficient to consider the infinite time limit of the multivariate PGF (4) As in Cong [5], it is necessary and sufficient for the limiting function to be a PGF that the multivariate limit lim where (z 1 , . . ., z J ) → (1, . . ., 1) along any path through Ω.Since the generating function is analytic inside Ω, this is equivalent to the condition lim For notational convenience, we denote the integrand and observe from Equation ( 7) that Since G S is a multivariate PGF, for fixed τ , H(z, τ ) is a decreasing function of z in the domain [0, 1].In particular, for each τ ≥ 0 we have pointwise convergence H(z, τ ) → 0 as z → 1, with Then by the dominated convergence theorem, We thus see that E[W ] < ∞ is a sufficient condition for Equation (8) to hold.Now, suppose that the network is ergodic.Then the limiting probability of the queue being empty can be written In direct analogy to Lemma 1 of Cong [5], under the constraint that the time spent in the network by each customer has finite expectation -which is necessary for batch occupancy time to have finite expectation, that is, E[W ] < ∞ -we recover a sufficient condition for ergodicity.
Proof.Similarly to Cong [5], we observe that Therefore, from Equation ( 7), it follows that A finite expected network occupation time for each customer For M X /M/∞ queues, Cong [5] and Yajima, Phung-Duc, and Masuyama [30] have shown that a necessary and sufficient condition for stability is the batch size S has a finite logarithmic moment, that is, Given an exponential service time of mean duration µ per customer, the expected occupancy time for a batch can written where γ e denotes the Euler-Mascheroni constant [25], in the case of the M X /M/∞ queue We can generalise this observation to recover a sufficient condition for ergodicity when the network occupancy time distributions Q j (t) have exponentially-bounded tails, that is, there exists δ > 0, t 0 > 0 such that Q j (τ ) ≤ e −δτ for all τ ≥ t 0 and j ∈ {1, . . ., J}. Proof.By assumption, there exist δ > 0, t 0 > 0 such that since G S is a generating function.Equation ( 7) then yields the bound Therefore, a sufficient condition for E[W ] < ∞ -and, by Theorem 4.1, the ergodicity of the network -is that Using the binomial expansion and identity 0.155.4 of [8], we compute.
Therefore, E[log(S + 1)] < ∞ is a sufficient condition for the ergodicity of the network.
For Markovian queueing networks, where the time spent by a customer is each node is exponentiallydistributed, we can strengthen Corollary 4.1.2to recover a necessary and sufficient condition for ergodicity.
Proof.Sufficiency follows directly from Corollary 4.1.2.By assumption, there exists η > δ j and From Equation ( 9) [25], we observe that In a similar vein, we can recover a sufficient condition for ergodicity in the case of fat-tailed customer occupancy time distributions.
Corollary 4.1.4.Suppose there exists α > 1 and t 0 > 0 such that Q j (t) ≤ t −α for all j ∈ {1, . . ., J} and t ≥ t 0 .Then the network is ergodic if Proof.Using similar reasoning to Corollary 4.1.2,we obtain the lower bound Performing a change of variables u = 1 − t −α , we find that the condition is sufficient to ensure E[W ] < ∞, where the interchanging of the order of summation and integration is justified by Tonelli's theorem.Using the integral representation of the beta function and identity 8.384.1 of [8], we compute Substituting Equation ( 11) in (10) and interchanging the order of summation (which is justified by Tonelli's theorem) yields By identity 6.1.46 of [1] Γ(k + 1) where we have interchanged the order of summation.Noting that (identity 0.121 of [8]), we further deduce that

Recurrence relations for the PMF of the transient occupancy distribution
For a single batch arrival (M X /G/∞) queue with a homogeneous arrival process, Willmot and Drekic [27,28,29] recognised the PGF governing the number of customers at a fixed time t as that of a compound Poisson distribution.They observed that application of the Panjer-Adleson recursion scheme [19,24] yields recurrence relations for the PMF of the time-dependent queue length distribution, formulated with respect to the PMF of a distribution obtained by compounding the batch size with a binomial distribution.Using the transient PGF given by Equation ( 4), we can likewise recover a recurrence relation for the PMF of the time-dependent occupancy distribution, formulated in terms of the PMF of C(t) (Equation ( 2)), using an analogous argument to Panjer [19].
Theorem 5.1.(A recurrence relation for the multivariate PMF.) Proof.Following the proof of Theorem 3.1, given that there are M (t) = m arrivals in the interval [0, t), the PMF for the number of customers in each queue 1, . . ., J at time t attributable to each arrival event is i.i.d. with PMF Denote by f (m) the m-fold convolution of f .Then by the law of total probability, the PMF of N(t) can be written in the form Now, by symmetry, that is, given that the sum of m i.i.d.random vectors, each with PMF f , is (n 1 , . . ., n J ), the conditional mean of the v th element is n v /m (Relation II of Panjer [19]).
Using Equations ( 13) and ( 14), for (n 1 , . . ., n J ) = 0, we obtain the expression where interchanging the order of summation is justified since the series is absolutely convergent.
Sequential application of Equation ( 15) along the lower triangular integer lattice yields precisely the recurrence relation (12).The expression (12) can also be obtained by computing partial derivatives of the multivariate PGF (4), using the formulae provided in Miatto [18].
Through an application of Faa di Bruno's formula, we noted in Mehra et al. [16] that the transient distribution for the M X t /G/∞ queue can also be formulated in terms of Bell polynomials [4], which yield an analogous recurrence structure.For queueing networks with multivariate batch arrivals, alternative representations of the transient PMF in terms of multivariate Bell polynomials [21] can also be recovered.
Likewise, we can derive recurrence relations for the joint factorial moments of N(t), formulated with respect to those of C(t).

Some special cases
The recurrence relations derived in Theorem 3.1 are formulated in terms of the PMF for the distribution of C(t).The tractability of these formulae are therefore constrained by the multivariate PMF for C(t).Here, we provide explicit formulae for several special cases.

Univariate batch sizes
For M X /G/∞ queues, Willmot and Drekic [27,28] observed that batch sizes within 'Sundt and Jewell's family of discrete distributions' [23,26] may yield tractable formulae.This class of distributions has a probability mass function that is characterised by a recurrence relation of the form negative binomial and logarithmic distributions, and zero-inflated analogues thereof.
In the context of queueing networks, we restrict our attention to multvariate batch size distributions S that can be expressed in the form where S is a univariate random variable within Sundt and Jewell's family [23,26] and p j ≥ 0 with j p j = 1.This corresponds to a univariate batch, with each incoming customer independently assigned an entry point j ∈ {1, . . ., J}.For notational convenience, we set to be the probability that a customer is situated in queue k time t after their arrival into the network.In the corollaries below, we state explicit formulae for binomial, Poisson, negative binomial and logarithmic batch sizes.

Constant batch sizes
Suppose S = s is constant.Then C(t−τ ) is the sum of independent, but non-identical categorical random variables, that is, C(t − τ ) follows a Poisson multinomial distribution [12].Using the exact DFT-CT method proposed by Lin, Wang, and Hong [12], we can recover the PMF for C(t − τ ), which can then be plugged into the recurrence relation (12) to compute the transient occupancy distribution.

A model of superinfection and hypnozoite dynamics for vivax malaria
In Mehra et al. [16], we constructed an open network of infinite server queues with exponential service times to model within-host superinfection and hypnozoite dynamics for Plasmodium vivax malaria.A schematic of the queueing network, with nodes labelled {1, . . ., k, N L, A, D, C, P, P C}, is replicated from Figure 3 of Mehra et al. [16] below.
The batch takes the form of a single arrival into queue P with probability p prim (representing a primary blood-stage infection) and a geometrically-distributed arrival into queue 1 (representing the establishment of dormant hypnozoites).Each hypnozoite then undergoes a compartmental process, accounting for death (queue D), or progression through successive latency compartments (queues 2, . . ., k) and subsequent activation giving rise to a blood-stage relapse (queue A).We further account for the clearance of each primary infection (queue P to P C) and relapse (queue A to C). Analysis of the queueing network allows for the characterisation of quantities of epidemiological interest, with applications to population-level transmission modelling [17].

Conclusion
We

Corollary 4 . 1 . 1 .
Suppose the network occupation time distributions Q j (t) have finite expectation.Then the network is ergodic if and a finite mean batch size E[S j ] < ∞ for each j ∈ {1, . . ., J} is thus a sufficient condition for E[W ] < ∞ and, by Theorem 4.1, the ergodicity of the network.

Figure 1 :
Figure 1: Model schematic replicated from Figure 3 of Mehra et al.[16] have extended results for single infinite-server queues with batch arrivals and open networks of infinite-server queues with single arrivals to a general network of infinite server queues with multivariate batch arrivals arriving according to a non-homogeneous Poisson process.Theorem 3.1 gives an expression for the PGF of the transient queue length distribution, Theorem 4.1 a necessary and sufficient condition for ergodicity and Theorem 5.1 a recurrence relation for the multivariate probability mass function.