Infinite-server systems with Coxian arrivals

We consider a network of infinite-server queues where the input process is a Cox process of the following form: The arrival rate is a vector-valued linear transform of a multivariate generalized (i.e., being driven by a subordinator rather than a compound Poisson process) shot-noise process. We first derive some distributional properties of the multivariate generalized shot-noise process. Then these are exploited to obtain the joint transform of the numbers of customers, at various time epochs, in a single infinite-server queue fed by the above-mentioned Cox process. We also obtain transforms pertaining to the joint stationary arrival rate and queue length processes (thus facilitating the analysis of the corresponding departure process), as well as their means and covariance structure. Finally, we extend to the setting of a network of infinite-server queues.


Introduction
In queueing theory it is commonly assumed that the customer arrival process is a Poisson process. Some recent empirical studies (see, for example, [4] for references) suggest, however, that arrival processes may exhibit overdispersion, i.e., the variance of the number of arrivals in an interval is larger than the corresponding mean. This has triggered research on queueing systems with overdispersed input.
The focus of the present paper is on infinite-server queues with, as input process, a doubly stochastic process, also known as a Cox process. Such a process can be seen as a Poisson process in which the rate is not a constant; instead, the rate process { (t), t ∈ R} itself is a (nonnegative) stochastic process. As an immediate consequence of the law of total variance, Cox processes indeed are overdispersed. Before describing the main contributions of the present paper, we first provide a brief account of the existing literature in this area. Literature In [7], the arrival process of the infinite-server queues is a Cox process in which the arrival rate is a shot-noise process. More specifically, the jumps of the shot-noise process occur according to a homogeneous Poisson process and are i.i.d. (independent and identically distributed) with a general distribution, and between jumps, the shot-noise process decays exponentially fast. The main object of study in [7] is a feedforward network of infinite-server queues, and the main result is the joint transform of the shot-noise-driven arrival rates and the numbers of customers in the queues.
Infinite-server queues are also studied in [2] and [8], but the arrival process there is a Hawkes process, the so-called self-exciting process-jumps in the shot-noise process in turn increase the intensity of the occurrence of the jumps. Daw and Pender [2] present several interesting motivating examples. They consider deterministic jump sizes in the shot-noise process and study in particular the Hawkes/Ph/∞ and Hawkes/D/∞ queues, obtaining detailed expressions for moments and autocovariances. Koops et al. [8] allow generally distributed jump sizes. For the case of exponentially distributed service times, the joint Laplace and z-transform of the Hawkes intensity and the number of customers is characterized via a partial differential equation, and that PDE is exploited to obtain recursive expressions for (joint) moments of the Hawkes intensity and the number of customers. For the case of generally distributed service times, the Hawkes process is viewed as a branching process, which allows expressing the z-transform of the number of customers in terms of the solution of a fixed-point equation.

Main goals and results
In this paper we aim to develop a general framework for the study of infinite-server queues with, as input, a quite general Cox process and to provide an exact analysis of networks of such queues. The modeling framework of [7] is extended in several ways; in particular, the shot-noise process is multivariate and not driven by a compound Poisson process, but by a Lévy subordinator. The class of Lévy subordinators not only contains compound Poisson processes, but also, for example, linear drifts and Gamma processes (where the latter process has the notable property that it exhibits infinitely many jumps in any time interval of positive length).
The main results of the paper are compact, elegant, transform expressions for joint distributions of arrival rates and numbers of customers, and properties of the departure processes, which reveal an interesting generalization of the classical results for the number-of-customers process and the departure process in the M/G/∞ queue. While large parts of the paper are quite technical, we also try to make clear why a detailed analysis of this complicated network is possible. In a nutshell, crucial for the analysis are the following convenient properties: (i) Lévy processes have stationary and independent increments, (ii) in a shot-noise process with exponential decay, all shots simultaneously decay at the same exponential rate, independently of each other, and (iii) in (a network of) infinite-server queues, all customers move through the queues independently, without interfering with each other. Motivation Our motivation for this study is threefold. (i) We started studying queues with Coxian input as a framework naturally allowing the incorporation of overdispersion. (ii) When we realized that a detailed analysis of an isolated infinite-server queue with Coxian input is possible, we were motivated to develop a very general framework for studying networks of infinite-server queues with a multivariate Coxian input process that is of great generality but still allows a detailed exact analysis. Together with product-form networks and linear stochastic fluid networks (cf. [6]) these are some of the rare examples of queueing networks for which one can obtain the joint distribution of key performance measures. (iii) Finally, while infinite-server queues with overdispersed input have various applications, we were in particular inspired by Web server applications. Consider the number of simultaneous visitors to a popular Web site or online video [11,12]. While the arrival process of visitors to the Web site may well be a Poisson process, its rate may jump up due to some (external) event, then decay gradually, only to jump up again because of another event. Such an example formed one of the motivations for [2,7,8], which all study infinite-server queues with an overdispersed arrival process.
Cox input processes driven by a multivariate Lévy subordinator are quite relevant for the above-mentioned example of a Web site with a visit rate that jumps up because of a certain event and subsequently decays again. They allow one to take into account multiple, possibly related, events that affect the visit rate of the Web site. Those events could occur instantaneously (the Lévy subordinator being a compound Poisson process), but there might also, for example, be a linear (fluid) component.
Since we also consider networks of infinite-server queues, we can also model the situation in which a visitor to a Web site, after a while, clicks in order to move to another Web site. Organization of the paper In Sect. 2 we introduce the multivariate shot-noise process and we describe distributional properties of this process. Section 3 studies a single infinite-server queue with, as input process, the above-described Cox process. The general case of a network of infinite-server queues is dealt with in Sect. 4. Each of these three sections starts with a brief overview of related known results: the onedimensional shot-noise process (Sect. 2.1), the classical M/G/∞ queue (Sect. 3.1), and a tandem network of M/G/∞ queues (Sect. 4.1). Section 5 contains some conclusions and suggestions for further research.

Multivariate shot-noise
In this preliminary section, we first briefly introduce shot-noise (Sect. 2.1) and define multivariate shot-noise and its stationary version (Sect. 2.2), then we recall some basic facts on Poisson random measures (Sect. 2.3), and we conclude in Sect. 2.4 with describing distributional properties of the multivariate shot-noise process. In between jumps, the process decreases with a state-dependent rate r X(t). Possible representations are hence, with T 1 , T 2 , . . . denoting the jump epochs of the compound Poisson process,

Shot-noise processes
alternatively, X (·) is the unique solution of the stochastic integral equation (2) Fig. 1 The rate X (t) as a function of time t

Definition of multivariate shot-noise
After having defined the classical shot-noise process, we now describe the generalized version of the shot-noise process that we work with in this paper. Let X(·) be a (generalized) d-dimensional multivariate shot-noise process, which is defined as follows: Let J(·) be a d-dimensional subordinator, i.e., a d-dimensional Lévy process which is non-decreasing in all components. Examples of such nondecreasing Lévy processes include linear drifts, compound Poisson processes and Gamma processes (but, obviously, not Brownian motion). We define the exponent of J(·) by −η(·); it satisfies, for α ∈ R d + and with denoting transposition, where c ∈ R d + and ν is an associated Lévy measure satisfying cf. [9,13]. Also, let Q = (q i j ) be a (d ×d)-matrix with nonnegative diagonal and non-positive off-diagonal entries, and with all eigenvalues having strictly positive real parts. An example of such a matrix is Q = (I − P )D, where D is a positive diagonal matrix and P is a substochastic matrix satisfying P n → 0 as n → ∞. A detailed motivation for studying this setting is provided in [6,Section 4].
We now introduce, for X(0) componentwise strictly positive and independent of J(·), the multivariate shot-noise process through cf. (1) and (2) above. Alternatively, X(·) is the unique (strong) solution to the stochastic integral equation The cumulative external input to station i is X i (0) + J i (t), the internal input rate from station i to station j is j =i q i j X j (t), and the output rate (some of which goes to other stations, while the rest leaves the system) is q i X i (t), where q i = −q ii ≥ 0. The reason for the above assumption that the eigenvalues of Q have a strictly positive real part is that it ensures that all the coordinates of the matrix e −Qt vanish as t → ∞ and thus, from (3), it follows that the effect of each arriving particle vanishes as well.
The process X(·) is Markovian and strictly positive (i.e., never hits or crosses any of the axes). Moreover, it has a unique stationary/ergodic/limiting distribution. This limiting distribution is identical to that of So far we have considered J(t) for positive values of t, but we can extend J(·) to the whole real line (with J(0) = 0). It follows directly that is a stationary process satisfying the same conditions that X does, that is, relations (3) and (4). From here on we assume that X = X * . In Sect. 3 we shall consider an infinite-server queue with as arrival process a nonhomogeneous Poisson process with arrival rate function (t) = a X(t) at time t. But first, in the next subsection, we review some well-known properties and results for Poisson random measures. We need those properties and results in order to provide an elegant treatment of the transforms of the arrival process and numbers of customers in networks of infinite-server queues with Coxian arrivals. In Sect. 2.4 we already demonstrate that by deriving an expression for the d-dimensional Laplace-Stieltjes transform (LST) of X.

Properties of Poisson random measures
Readers who are not familiar with Poisson random measures are referred to, for example, Chapter 6 of Çinlar [1]. To proceed, we first recall that a Poisson random measure N on some measurable space (X, G ) is a random measure having the property that, for every pairwise disjoint A 1 , . . . , A n ∈ G , the random variables N ( is defined to be almost surely 0 or ∞, respectively. We note that it is well known that N is a Poisson random measure with a (σ -finite) mean measure μ if and only if for any G -measurable f : Observe that E f (s) N (ds) = f (s) dμ(s) (which holds even if N is not Poisson). Finally, it is also known, and actually easy to check (first for indicators, then simple functions, etc.), that if f 1 (·) and f 2 (·) are nonnegative G -measurable, then Clearly, when It is also well known (see, for example, Chapter 4 of [13]) that there exists a Poisson random measure and

g(s) h(s)ds,
with provided that the corresponding variances R g(s) g(s)ds and R h(s) h(s)ds are both finite. With σ i j denoting the (i, j)th coordinate of , it also holds that where ν i j is the (i, j)th marginal measure associated with ν (i.e., the Lévy measure associated with the (i, j)th coordinate of J(·)).

Distributional properties of multivariate shot-noise
Combining (5) with the powerful identity (8), we immediately obtain the following In the usual manner, moments can be found from the LST. The first moment takes the following form: With ∇η(0) the vector of first derivatives, We now identify the covariance matrix of X. From (9) it follows, as was also shown in [6,Thm. 5.2], that the covariance matrix of X is given by and is the unique solution of The next objective is to find an expression for the covariance between X(t) and X(t + h), again bearing in mind that we started the process at −∞. Now, clearly is independent of X(t), due to the independent increments property of J(·). As a consequence, with It now follows from (9) that, by taking g(s) := e −Q (t−s) x and h(s) := e −Q (t+h−s) y, for every x, y ∈ R d . It thus follows, using (11), that Combining the above, we find the following result.
Since Q and e −Q h commute, it also follows from (12) that In fact, employing the method of proof of [6, Thm. 5.2], h is the unique solution of this equation.

Remark 2.3
In the one-dimensional case, letting σ 2 and q denote the one-dimensional versions of the matrices and Q, respectively, we conclude that the covariance between X (t) and X (t + h) is given by 1 2 e −qh σ 2 /q.

Single infinite-server queue with a Cox input process
We consider an infinite-server queue in which, conditioned on J(·), the arrival process is a non-homogeneous Poisson process with rate function (t) = a X(t) at time t, where a ∈ R d + . It is throughout assumed that service times are i.i.d. and are independent of J(·) (hence also of X(·)) and the arrival process. The service times have distribution function F(·) and complementary distribution functionF(·), and in addition we define F(s, t) := F(t) − F(s) for s < t.
In this section, we start by reminding the reader of some well-known results concerning the classical M/G/∞ queue (Sect. 3.1), subsequently we derive, for the infinite-server queue with a Cox input process, the transform of queue lengths at various time epochs jointly with the numbers of arrivals in various time intervals (Sect. 3.2), and we conclude, in Sect. 3.3, with several calculations yielding moments as well as the joint steady-state transform of the number of customers in the infinite-server queue and the arrival rate.

The M/G/∞ queue
In this subsection we consider the classical M/G/∞ queue with fixed arrival rate λ and generic service time B. It is well known [3] that the number of customers L(t) at time t, starting with L(0) = 0, is Poisson distributed with mean This is easily seen by observing that the number of arrivals in [0, t] is Poisson(λt), that different customers do not interact and that, given that there were n arrivals in [0, t], the arrival epochs were independent and uniformly distributed on [0, t]. So the Poisson arrival process is thinned, and the number of customers still present at t follows a Poisson distribution. Thus, the generating function of L(t) immediately follows: The pleasing properties of non-interfering customer behavior and of the uniformly distributed arrival epochs also quickly lead to elegant results for joint distributions of customer numbers at various epochs. Another well-known result about the M/G/∞ queue, going back to Mirasol [10], is that in stationarity the departure process of the queue also is Poisson. We now study the question of how these results change in the setting of the Cox input process.

Transform of queue lengths and numbers of arrivals
Our first objective is to derive the transform of queue lengths at various points in time, jointly with the number of arrivals in corresponding intervals.
To this end, we consider n ∈ N time intervals, say (t 0 , t 1 ] up to (t n−1 , t n ] (where it is assumed that t i−1 < t i for i = 1, . . . , n). Our goal is to establish the joint transform of the queue lengths L i at each of the t i and the numbers of arrivals A i in each of the intervals (t i−1 , t i ], i.e., we shall compute the transform, for w ∈ R n+1 and z ∈ R n , Notice that a job arriving in the interval (t i−1 , t i ] contributes to A i (and does not contribute to any of the other A j ) and potentially contributes to L i up to L n . More precisely, supposing that the job arrives at s ∈ (t i−1 , t i ], with probability F(t j − s, t j+1 − s) it contributes to L i up to L j , for j ∈ {i, . . . , n} and defining t n+1 := ∞.

Conditional on (·) = λ(·) this concerns a Poisson contribution with parameter
conditional on (·) = λ(·) all these contributions are independent. In addition, we recall that the probability-generating function (evaluated in z) of a Poisson random variable with parameter λ equals exp(λ(z − 1)). Combining the above elements, we obtain where, with t −1 := −∞, Representation (15) generally holds, i.e., for any nonnegative arrival rate process (·). For the case of multivariate shot-noise, however, the expression can be made more explicit. To this end, we substitute (cf. By applying Eq. (10), we obtain that (15) equals | w, z)

ds a dr .
We have thus established the following result. Theorem 3.1 For any w ∈ R n+1 and z ∈ R n , Interestingly, the above result directly enables us to describe the distribution of the number of departures in all intervals. Let D i be the number of departures in (t i−1 , t i ]. Then we would like to evaluate , . . . , n − 1} and w n (v) := v −1 n , we thus find the following result.

Explicit calculations
Let L(t) be the number of customers present at time t. In this subsection we compute (i) the mean and variance of L(0), (ii) the joint transform of (0) and L(0), and (iii) the covariance between L(0) and L(t) for some t > 0. As before, we assume that the process started at −∞, so that it displays stationary behavior at time 0 (and consequently at t as well). In principle, (factorial) moments can be derived by differentiating (w, z) suitably often and plugging in w = 1 and z = 1, but elegant direct arguments can be given, as we show now. We first introduce some notation. Define β as the mean service time: The density of the residual service time B e is given by f e (s) :=F(s)/β.

• Mean and variance of 3(0) and L(0)
It is easily checked that E (0) = λ := a Q −1 and Var (0) = a 0 a. We therefore now concentrate on the mean and variance of L(0). Let P(μ) denote a Poisson random variable with parameter μ. It is straightforward that We thus find, with J shorthand notation for the entire path of the process J(·), For later use, we rewrite this expression as where the last step is due to time-reversibility of J(·). Applying (7), we obtain Using the law of total variance, (16), (18) and (9), Expression (19) can be substantially simplified, which we do later in this section.

• Joint transform of 3(0) and L(0)
Here the goal is to determine E e −v (0) w L(0) . By arguments similar to those used in Sect. 3, Applying (8), we thus obtain the following result.
It is possible to apply this joint transform to the computation of moments of L(0) and (0), as well as their mixed moments, but lower moments can typically be computed by direct arguments. The mean and variance of (0) and L(0) have been identified above. We therefore now focus on the covariance between (0) and L(0). Since as X(·) is a functional of J(·)), and since B e is independent of J , we find where we employ the fact that the covariance matrix between X(t) and X(t + h), which is the same as that between X(t − h) and X(t), is given by h = 0 e −Q h .

• Covariance between L(0) and L(t)
The starting point is the law of total covariance: We evaluate the two terms separately. The first term can be rewritten as follows: Let C 1 (t) denote the number of customers that arrive before time 0 and depart in (0, t]; C 2 (t) is the number of customers that arrive before time 0 and depart after t; finally, C 3 (t) is the number of customers that arrive in (0, t] and depart after t. Evidently, due to the conditional independence of these three quantities, with, mimicking the above arguments, where the last equality is due to the fact that the conditional distribution of C 2 (t) given J is Poisson. Thus, with F e (·) denoting the distribution function of B e andF e (·) the corresponding complementary distribution function, Next, we move to the second term. To this end, we first recall (cf. (16) and (17)) that As J has independent increments, when considering the covariance between where the last step follows by using (9). As it turns out, the last formula (as well as (19)) may be simplified, as follows: Let B e,1 and B e,2 be two i.i.d. copies of B e . Recalling (12) and denoting x + := x ∨ 0 and where, in the last equality, (11) has been used. For any s ∈ R we have that, recalling Equation (13), a e −Qs − 0 e −Q s + a = a 0 e −Q |s| a = a |s| a.
Thus, denoting by R(t) the autocorrelation function (with R(0) = 1), then adding the two terms yields In particular, when t = 0, (21) provides us with a simplified expression for Var L(0). The density of B e,1 − B e,2 is clearly symmetric around zero and is given by Since f e (·) = β −1F (·) is non-increasing on [0, ∞), this implies that g(·) is unimodal. We end this subsection by considering the one-dimensional case. Since both h(x) := e −q|x| (q > 0) and g(x) are symmetric and unimodal (around zero), it follows by [15] that so is their convolution. Alternatively, this follows also from [5] as h(·) is also log-concave. This implies that is unimodal with a mode at zero, hence non-increasing on [0, ∞), and thus so is R(·). A similar result holds if a is an eigenvector associated with a real-valued (positive) eigenvalue q of Q since then a 0 e −Q |x| a = a 0 a e −q|x| . However, in general, even though a 0 e −Q |x| a is a symmetric function, it is not clear if it is decreasing on [0, ∞) or if it is log-concave. However, since it vanishes as |x| → ∞, it is clear that R(t) vanishes as t → ∞.

The network case
In this section we consider networks of infinite-server queues with a Cox input process. There are n queues, with the arrival rate of queue k being a non-homogeneous Poisson process with rate function k (t) = a k X(t). This is the same as defining the vector of arrival processes to be (t) = AX(t) for some matrix A with rows a k . The X(t) process is the same multivariate shot-noise process as before. Notice that in this construction the arrival processes at the various queues are potentially dependent.
Define by p km (t) the probability that a job entering at queue k at time 0 is at queue m at time t; likewise, p k0 (t) is the probability that a job entering at queue k at time 0 has left the network by time t. In the case of exponentially distributed service times and probabilistic routing, these p km (t)'s can be computed more explicitly (relying on the machinery developed for phase-type distributions).
Before studying the joint queue length distribution in the network (Sect. 4.2) and presenting an example (Sect. 4.3), we briefly remind the reader of the structure of the joint queue length distribution in the case of two M/G/∞ queues in tandem (Sect. 4.1). This helps sharpen our intuition and will make it easier to digest and understand the results of Sect. 4.2.

Two M/G/∞ queues in tandem
Consider a model of two M/G/∞ queues 1, 2 in series, with arrival rate λ at queue 1, generic service time B i at queue i, i = 1, 2, and all customers moving from queue 1 to 2 before leaving the system. Queue 2 has no external arrivals. The same concepts that were underlying the results mentioned in Sect. 3.1 (see, in particular, (14)) quickly lead to the following expression for the generating function of the joint distribution of the two queue lengths L 1 (t), L 2 (t): Furthermore, the generating function of the steady-state joint queue length distribution is given by

Joint queue length
We focus on analyzing the stationary joint queue length distribution. In principle, virtually all quantities studied in the previous section can be derived again, at the expense of introducing rather heavy notation.
In this section, we let K m denote the stationary queue length at node m ∈ {1, . . . , n}. Our objective is to compute the joint probability-generating function Using the same arguments as in the previous section, we conclude that K m has a mixed Poisson distribution. More precisely, K m has a Poisson distribution with (random) parameter 0 −∞ n k=1 k (s) p km (−s) ds.
It follows that Substituting (cf. (6)) in (23) gives This leads to the following result.

Example
In this illustrative example we consider a two-node tandem system in which the service times at both nodes are exponential with parameter κ > 0. We assume for ease that there is only input at the first queue and that all customers move from the first queue to the second queue. It is immediate that, for t 0, use that the sojourn time in the system is Erlang (2). We find that Assume for ease that d = 1 (i.e., one-dimensional shot-noise). We choose a 1 = 1 (whereas a 2 = 0, as we assumed no external input to the second queue). As a consequence, we have that 1 (t) = X (t) and 2 (t) = 0, where X (t) can be represented as for some q > 0 and a (scalar) Lévy subordinator J (·) (with exponent −η(·)). Appealing to Thm. 4.1, To evaluate this expression, observe that It is easy to obtain moments of K 1 and K 2 from (24) and (25) by differentiation with respect to w 1 and w 2 , respectively. In particular, straightforward algebra yields Remark 4.2 Similar calculations can be performed for considerably more general networks. Application of Thm. 4.1 requires (i) expressions for the matrix exponential of Q and (ii) expressions for π k (s | w) (and therefore for p km (s)) for s 0. To this end, evaluate, using common techniques from linear algebra, the matrix exponential of Q ; this can be done relying on standard computer algebra. Assuming the eigenvalues of Q to have multiplicity 1 (where we note that higher multiplicities can be dealt with too, but require some additional care), we thus obtain a matrix whose entries are weighted sums of exponentials. Likewise, if the service times are of phase type, then the probabilities p km (s) can be written in terms of (specific entries of) a matrix exponential, and therefore (again under mild conditions) also as a weighted sum of exponentials. It thus follows that the integrand e −Q (s−r ) π k (s | w) is a mixture of exponentials (in s) and can be integrated by elementary calculus. This allows the computation of the (mixed) moments pertaining to the stationary queue lengths.
The probability-generating function of the queue length in the first queue is therefore We thus end up with .
A similar procedure yields the probability-generating function of the second queue, but the resulting expressions are slightly more complicated.

Conclusions and suggestions for further research
In this paper we have considered a network of infinite-server queues where the input process is a Cox process that allows much modeling flexibility; the arrival rate is represented as a linear combination of the components of a multivariate generalized shot-noise process. We have derived various distributional properties of the multivariate shot-noise process, subsequently exploiting them to obtain the joint transform of the numbers of customers, at consecutive time epochs, in an infinite-server queue with, as input process, such a Cox process. We have also derived the joint steadystate transform of the vectors of arrival rate and queue length, as well as their means and covariance structure, and we have studied the departure process from the queue. Finally, we extended our analysis to the setting of a network of infinite-server queues, allowing the arrival processes at the various queues to be dependent Cox processes.
In a follow-up study [14] we investigate several related aspects. Specifically, (i) we develop a recursive scheme that allows us to obtain higher-order moments of ( (t), L(t)), (ii) we derive asymptotics of the queue length process, under assumptions regarding the tail behavior of the shot-noise process, (iii) we study the heavy-traffic behavior of the queue length process, and (iv) we obtain several numerical results for means, variances and correlations; this includes a numerical exploration of the sensitivity of various performance metrics with respect to the service time distribution. It could also be investigated to what extent the assumption of exponential decay in the shot-noise process can be relaxed.