Networks of \(\cdot /G/\infty \) queues with shot-noise-driven arrival intensities

Queueing Systems, Mar 2017

We study infinite-server queues in which the arrival process is a Cox process (or doubly stochastic Poisson process), of which the arrival rate is given by a shot-noise process. A shot-noise rate emerges naturally in cases where the arrival rate tends to exhibit sudden increases (or shots) at random epochs, after which the rate is inclined to revert to lower values. Exponential decay of the shot noise is assumed, so that the queueing systems are amenable to analysis. In particular, we perform transient analysis on the number of jobs in the queue jointly with the value of the driving shot-noise process. Additionally, we derive heavy-traffic asymptotics for the number of jobs in the system by using a linear scaling of the shot intensity. First we focus on a one-dimensional setting in which there is a single infinite-server queue, which we then extend to a network setting.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Networks of \(\cdot /G/\infty \) queues with shot-noise-driven arrival intensities

Networks of ·/ G/∞ queues with shot-noise-driven arrival intensities D. T. Koops 0 1 O. J. Boxma 0 1 M. R. H. Mandjes 0 1 0 Eurandom and Department of Mathematics and Computer Science, Eindhoven University of Technology , Eindhoven , The Netherlands 1 Korteweg-de Vries Institute, University of Amsterdam , Amsterdam , The Netherlands We study infinite-server queues in which the arrival process is a Cox process (or doubly stochastic Poisson process), of which the arrival rate is given by a shot-noise process. A shot-noise rate emerges naturally in cases where the arrival rate tends to exhibit sudden increases (or shots) at random epochs, after which the rate is inclined to revert to lower values. Exponential decay of the shot noise is assumed, so that the queueing systems are amenable to analysis. In particular, we perform transient analysis on the number of jobs in the queue jointly with the value of the driving shot-noise process. Additionally, we derive heavy-traffic asymptotics for the number of jobs in the system by using a linear scaling of the shot intensity. First we focus on a onedimensional setting in which there is a single infinite-server queue, which we then extend to a network setting. Infinite-server queue; Stochastic arrival rate; Cox; Doubly stochastic Poisson; Shot noise; Network; Functional central limit theorem; Heavy traffic 1 Introduction In the queueing literature, one has traditionally studied queues with Poisson input. The Poisson assumption typically facilitates explicit analysis, but it does not always align well with actual data; see, for example, [11] and references therein. More specifically, statistical studies show that in many practical situations, Poisson processes underestimate the variability of the queue’s input stream. This observation has motivated research on queues fed by arrival processes that better capture the burstiness observed in practice. The extent to which burstiness takes place can be measured by the dispersion index, i.e., the ratio of the variance of the number of arrivals in a given interval, and the corresponding expected value. In arrival streams that display burstiness, the dispersion index is larger than unity (as opposed to Poisson processes, for which it is equal to unity), a phenomenon that is usually referred to as overdispersion. It is desirable that the arrival process of the queueing model takes the observed overdispersion into account. One way to achieve this is to make use of Cox processes, which are Poisson processes, conditional on the stochastic time-dependent intensity. It is an immediate consequence of the law of total variance that Cox processes do have a dispersion index larger than unity. Therefore, this class of processes makes for a good candidate to model overdispersed input processes. In this paper we contribute to the development of queueing models fed by input streams that exhibit overdispersion. We analyze infinite-server queues driven by a particular Cox process, in which the rate is a (stochastic) shot-noise process. The shot-noise process we use is one in which there are only upward jumps (or shots) that arrive according to a homogeneous Poisson process. Furthermore, we employ an exponential ‘response’ or ‘decay’ function, which encodes how quickly the process will decline after a jump. In this case, the shot-noise process is a Markov process; see [16, p. 393]. There are several variations on shot-noise processes; see, for example, [10] for a comprehensive overview. It is not a novel idea to use a shot-noise process as stochastic intensity. For instance, in insurance mathematics, the authors of [5] use a shot-noise-driven Cox process to model the claim count. They assume that disasters happen according to a Poisson process, and each disaster can induce a cluster of arriving claims. The disaster corresponds to a shot upwards in the claim intensity. As time passes, the claim intensity process decreases, as more and more claims are settled. Another example of shot-noise arrival processes is found in the famous paper [13], where it is used to model the occurrences of earthquakes. The arrival process considered in [13] has one crucial difference with the one used in this paper: it makes use of Hawkes processes [9], which do have a shot-noise structure, but have the special feature that they are self-exciting. More specifically, in Hawkes processes, an arrival induces a shot in the arrival rate, whereas in our shot-noise-driven Cox model these shots are merely exogenous. The Hawkes process is less tractable than the shot-noisedriven Cox process. A very recent effort to analyze ·/G/∞ queues that are driven by a Hawkes process has been made in [8], where a functional central limit theorem is derived for the number of jobs in the system. In this model, obtaining explicit results (in a non-asymptotic setting), as we are able to do in the shot-noise-driven Cox variant, is still an open problem. In order to successfully implement a theoretical model, it is crucial to have methods to estimate its parameters from data. The shot-noise-driven Cox process is attractive since it has this property. Statistical methods that filter the unobservable intensity process, based on Markov Chain Monte Carlo (MCMC) techniques, have been developed; see [3] and references therein. By filtering, they refer to the estimation of the intensity process in a given time interval, given a realized arrival process. Subsequently, given this ‘filtered path’ of the intensity process, the parameters of the shot-noise process can be estimated by a Monte Carlo version of the expectation maximization (EM) method. Furthermore, the shot-noise-driven Cox process can also be easily simulated; see, for example, the thinning procedure described in [12]. In this paper we study networks of infinite-server queues with shot-noise-driven Cox input. We assume that the service times at a given node are i.i.d. samples from a general distribution. The output of a queue is routed to a next queue, or leaves the network. Infinite-server queues have the inherent advantage that jobs do not interfere with one another, which considerably simplifies the analysis. Furthermore, infinite-server systems are frequently used to produce approximations for corresponding finite-server systems. In the network setting, we can model queueing systems that are driven by correlated shot-noise arrival processes. With regards to applications, such a system could, for example, represent the call centers of a fire department and police department in the same town. The contributions and organization of this paper are as follows. In this paper we derive exact and asymptotic results. The main result of the exact analysis is Theorem 4.6, where we find the joint Laplace transform of the numbers of jobs in the queues of a feedforward network, jointly with the shot-noise-driven arrival rates. We build up toward this result as follows. In Sect. 2 we introduce notation and we state the important Lemma 2.1 that we repeatedly rely on. Then we derive exact results for the single infinite-server queue with a shot-noise arrival rate in Sect. 3.1. Subsequently, in Sect. 3.2, we show that after an appropriate scaling the number of jobs in the system satisfies a functional central limit theorem (Theorem 3.4); the limiting process is an Ornstein–Uhlenbeck (OU) process driven by a superposition of a Brownian motion and an integrated OU process. We then extend the theory to a network setting in Sect. 4. Before we consider full-blown networks, we first consider a tandem system consisting of an arbitrary number of infinite-server queues in Sect. 4.1. Then it is argued in Sect. 4.2 that a feedforward network can be seen as a number of tandem queues in parallel. We analyze two different ways in which dependency can enter the system through the arrival process. Firstly, in Model (M1), parallel service facilities are driven by a multidimensional shot-noise process in which the shots are simultaneous (which includes the possibility that all shot-noise processes are equal). Secondly, in Model (M2), we assume that there is one shot-noise arrival intensity that generates simultaneous arrivals in all tandems. In Sect. 5 we finish with some concluding remarks. 2 Notation and preliminaries Let ( , F , {Ft }t≥0, P) be a probability space, in which the filtration {Ft }t≥0 is such that (·) is adapted to it. A shot-noise process is a process that has random jumps at Poisson epochs, and a deterministic ‘response’ or ‘decay’ function, which governs the behavior of the process. See [16, Section 8.7] for a brief account of shot-noise processes. The shot noise that we use in this paper has the following representation: Fig. 1 Sample path of shot-noise process Fig. 2 A realization of arrival process corresponding to the sample path of the arrival rate process in Fig. 1 (t ) = (0)e−r t + i=1 Bi e−r (t−ti ), where the Bi ≥ 0 are i.i.d. shots from a general distribution, the decay function is exponential with rate r > 0, PB is a homogeneous Poisson process with rate ν, and the epochs of the shots, that arrived before time t , are labeled t1, t2, . . . , tPB (t). As explained in the introduction, the shot-noise process serves as a stochastic arrival rate to a queueing system. It is straightforward to simulate a shot-noise process; for an illustration of a sample path, consider Fig. 1. Using the thinning method for nonhomogeneous Poisson processes [12], and using the sample path of Fig. 1 as the arrival rate, one can generate a corresponding sample path for the arrival process, as is displayed in Fig. 2. Typically, most arrivals occur shortly after peaks in the shot-noise process in Fig. 1, as expected. We write (i.e., without argument) for a random variable with distribution equal to that of limt→∞ (t ). We now present well-known transient and stationary moments of the shot-noise process, see Appendix 2 and, for example [16]: with B distributed as B1, (t ) = (0)e−r t + (t ) = (1 − e−r t ), (1 − e−2r t ), We remark that, for convenience, we throughout assume (0) = 0. The results can be readily extended to the case in which (0) is a nonnegative random variable, at the cost of a somewhat more cumbersome notation. In the one-dimensional case, we denote β(s) = E e−s B , and in the multidimensional case, where s = (s1, s2, . . . , sd ), for some integer d ≥ 2, now denotes a vector, we write The following lemma will be important for the derivation of the joint transform of (t ) and the number of jobs in system, in both single- and multi-node cases. Lemma 2.1 Let (·) be a shot-noise process. Let f : R × Rd → R be a function which is piecewise continuous in its first argument, with at most a countable number of discontinuities. Then it holds that = exp ⎜ ν ⎝ t ⎛ f (u, z) (u) du − s (t ) f (u, z)e−ru du Proof See Appendix 1. 3 A single infinite-server queue In this section we study the MS/G/∞ queue. This is a single infinite-server queue, of which the arrival process is a Cox process driven by the shot-noise process (·), as defined in Sect. 2. First we derive exact results in Sect. 3.1, where we find the joint transform of the number of jobs in the system and the shot-noise rate, and derive expressions for the expected value and variance. Subsequently, in Sect. 3.2, we derive a functional central limit theorem for this model. 3.1 Exact analysis We let Ji be the service requirement of the i -th job, where J1, J2, . . . are assumed to be i.i.d.; in the sequel J denotes a random variable that is equal in distribution to J1. Our first objective is to find the distribution of the number of jobs in the system at time t , in the sequel denoted by N (t ). This can be found in several ways; because of the appealing underlying intuition, we here provide an argument in which we approximate the arrival rate on intervals of length by a constant, and then let ↓ 0. This procedure works as follows. We let (t ) = (ω, t ) be an arbitrary sample path of the driving shot-noise process. Given (t ), the number of jobs that arrived in the interval [k , (k + 1) ) and are still in the system at time t has a Poisson distribution with parameter P( J > t − (k + Uk )) · (k ) + o( ), where U1, U2, . . . are i.i.d. standard uniform random variables. Summing over k yields that the number of jobs in the system at time t has a Poisson distribution with parameter t/ −1 k=0 P( J > t − (k which converges, as P( J > t − u) (u) du. The argument above is not new: a similar observation was mentioned in, for example [6], for deterministic rate functions. Since (·) is actually a stochastic process, we conclude that the number of jobs has a mixed Poisson distribution, with the expression in Eq. (3) as random parameter. As a consequence, we find by conditioning on Ft , ξ(t, z, s) := E z N (t)e−s (t) = E e−s (t) E z N (t) | Ft In addition, by the law of total variance we find Var N (t ) = Var We have found the following result. Theorem 3.1 Let (·) be a shot-noise process. Then t ⎛ P( J > t − u)e−ru du + se−r(t−v) Proof The result follows directly from Lemma 2.1 and Eq. (4). E N (t ) = E (u) P( J > t − u) du. = E exp (z − 1)P( J > t − u) (u) du − s (t ) . = 2 t/ −1 t/ −1 The latter expression we can further evaluate, using an approximation argument that resembles the one we used above. Using a Riemann sum approximation, we find (u) P( J > t − u) du = lim↓0 Var ⎝ i=0 0 v t E (u) P( J > t − u) du. We can make this more explicit using the corresponding formulas in (2). Example 3.2 (Exponential case) Consider the case in which J is exponentially distributed with mean 1/μ and (0) = 0. Then we can calculate the mean and variance explicitly. For μ = r , E N (t ) = ⎧ t → ⎪⎨ μ(1 − e−rt ) − r (1 − e−μt ) Var N (t ) = 1 − e−2rt − 2r t (1 + r t )e−2rt 3.2 Asymptotic analysis This subsection focuses on deriving a functional central limit theorem (FCLT) for the model under study, after having appropriately scaled the shot rate of the shot-noise process. In the following, we assume that the service requirements are exponentially distributed with rate μ, and we point out how it can be generalized to a general distribution in Remark 3.6 below. We follow the standard approach to derive the FCLT for infinite-server queueing systems; we mimic the argumentation used in, for example [1,14]. As the proof has a relatively large number of standard elements, we restrict ourselves to the most important steps. We apply a linear scaling to the shot rate of the shot-noise process, i.e., ν → nν. It is readily checked that under this scaling, the steady-state level of the shot-noise process and the steady-state number of jobs in the queue blow up by a factor n. It is our objective to prove that, after appropriate centering and normalization, the process recording the number of jobs in the system converges to a Gaussian process. In the n-th scaled model, the number of jobs in the system at time t , denoted by Nn(t ), has the following (obvious) representation: with An(t ) denoting the number of arrivals in [0, t ], and Dn(t ) the number of departures, Nn(t ) = Nn(0) + An(t ) − Dn(t ). Here, An(t ) corresponds to a Cox process with a shot-noise-driven rate, and therefore we have, with n(s) the shot-noise in the scaled model at time s and SA(·) a unit-rate Poisson process, An(t ) = SA Dn(t ) = SD in line with our previous assumptions, we put n(0) = 0. For our infinite-server model the departures Dn(t ) can be written as, with SD(·) a unit-rate Poisson process (independent of SA(·)), We start by identifying the average behavior of the process Nn(t ). Following the reasoning of [1], assuming that Nn(0)/n ⇒ ρ(0) (where ‘⇒’ denotes weak convergence), Nn(t )/n converges almost surely to the solution of E (u) du − Now we move to the derivation of the FCLT. Following the approach used in [1], we proceed by studying an FCLT for the input rate process. To this end, we first define ˆ n(t ) := √n n(t ) − E Kˆ n(t ) := ˆ (u) du. The following lemma states that Kˆ n(·) converges to an integrated Ornstein– Uhlenbeck (OU) process, corresponding to an OU process ˆ (·) with a speed of mean reversion equal to r , long-run equilibrium level 0, and variance σ 2 := ν E B2/(2r ). Lemma 3.3 Assume that for the shot sizes, distributed as B, it holds that E B, E B2 < ∞. Then Kˆ n(·) ⇒ Kˆ (·) as n → ∞, where Proof This proof is standard; for instance from [2, Prop. 3], by putting the λd in that paper to zero, it follows that ˆ n(·) ⇒ ˆ (·). This implies Kˆ n(·) ⇒ Kˆ (·), using (10) together with the continuous mapping theorem. Interestingly, the above result entails that the arrival rate process displays meanreverting behavior. This also holds for the job count process in standard infinite-server queues. In other words, the job count process in the queueing system we are studying can be considered as the composition of two mean-reverting processes. We make this more precise in the following. From now on we consider the following centered and normalized version of the number of jobs in the system: Nˆ n(t ) := √n in which ˆ satisfies, with W1(·) a standard Brownian motion, Kˆ (t ) = ˆ (u) du, ˆ (u) du. We assume that Nˆ n(0) ⇒ Nˆ (0) as n → ∞. To prove the FCLT, we rewrite Nˆ n(t ) in a convenient form. Mimicking the steps performed in [1] or [14], with S¯ A(t ) := SA(t ) − t , S¯D(t ) := SD(t ) − t , Rn(t ) := S¯ A − S¯D Nn(u)du , and using the relation (9), we eventually obtain Nˆ n(t ) = Nˆ n(0) + Our next goal is to apply the martingale FCLT to the martingales Rn(t )/√n; see, for background on the martingale FCLT, for instance [7] and [17]. The quadratic variation equals which converges to 0t E (u)du + μ 0t ρ(u)du. Appealing to the martingale FCLT, the following FCLT is obtained. Theorem 3.4 The centered and normalized version of the number of jobs in the queue satisfies an FCLT: Nˆ n(·) ⇒ Nˆ (·) as n → ∞, where Nˆ (t ) solves the stochastic integral equation Nˆ (t ) = Nˆ (0) + E (u) + μρ(u) dW2(u) + Kˆ (t ) − μ with W2(·) a standard Brownian motion that is independent of the Brownian motion W1(·) we introduced in the definition of Kˆ (·). Remark 3.5 In passing, we have proven that the arrival process as such obeys an FCLT. With we find that Aˆn(t ) ⇒ Aˆ(t ) as n → ∞, where Aˆn(t ) := √n An(t ) − Aˆ(t ) := Remark 3.6 The FCLT can be extended to non-exponential service requirements, by making use of [15, Thm. 3.2]. Their approach relies on two assumptions: • The arrival process should satisfy an FCLT; • The service times are i.i.d. nonnegative random variables with a general c.d.f. independent of the arrival process. As noted in Remark 3.5, the first assumption is satisfied for the model in this paper. The second assumption holds as well. In the non-exponential case, the results are less clean; in general, the limiting process can be expressed in terms of a Kiefer process, cf. for example [4]. 4 Networks Now that the reader is familiar with the one-dimensional setting, we extend this to networks. In this section, we focus on feedforward networks in which each node corresponds to an infinite-server queue. Feedforward networks are defined as follows. Definition 4.1 (feedforward network) Let G = (V , E ) be a directed graph with nodes V and edges E . The nodes represent infinite-server queues and the directed edges between the facilities demonstrate how jobs move through the system. We suppose that there are no cycles in G, i.e., there is no sequence of nodes, starting and ending at the same node, with each two consecutive nodes adjacent to each other in the graph, consistent with the orientation of the edges. We focus on feedforward networks to keep the notation manageable. In Theorem 4.6, we derive the transform of the numbers of jobs in all nodes, jointly with the shot-noise process(es) for feedforward networks. Nonetheless, we provide Example 4.5, to show that analysis is in fact possible if there is a loop, but at the expense of more involved calculations. Since all nodes represent infinite-server queues, one can see that whenever a node has multiple input streams, it is equivalent to multiple infinite-server queues that work independently from each other, but have the same service speed and induce the same service requirement for arriving jobs. Consider Fig. 3 for an illustration. The reason why this holds is that different job streams move independently through the system, without creating waiting times for others. Therefore, merging streams do not increase the complexity of our network. The same holds for ‘splits’ in job streams. By this we mean that after jobs finished their service at a server, they move to server i with probability qi (with i qi = 1). Then, one can simply sample the entire path that the job will take through the system, at the arrival instance at its first server. If one recognizes the above, then all feedforward networks reduce to parallel tandem systems in which the first node in each tandem system is fed by external input. The procedure to decompose a network into parallel tandems consists of finding all paths between nodes in which jobs either enter or leave the system. Each of these paths will subsequently be considered as a tandem queue, which are then set in parallel. Fig. 3 Since the jobs are not interfering with each other, the network on the left is equivalent to the graph on the right. Node 3’ is a copy of node 3: It works at the same speed and induces the same service requirements To build up to the main result, we first study tandem systems in Sect. 4.1. Subsequently, we put the tandem systems in parallel in Sect. 4.2 and finally we present the main theorem and some implications in Sect. 4.3. 4.1 Tandem systems As announced, we proceed by studying tandem systems. In Sect. 4.2 below, we study d parallel tandem systems. In this subsection, we consider the i -th of these tandem systems, where i = 1, . . . , d. Suppose that tandem i has Si service facilities and the input process at the first node is Poisson, with a shot-noise arrival rate i (·). We assume that jobs enter node i 1. When they finish service, they enter node i 2, etc., until they enter node i Si after which they leave the system. We use i j as a subscript referring to node j in tandem system i and we refer to the node as node i j . Hence Ni j (t ) and Ji j denote the number of jobs in node i j at time t , and a copy of a service requirement, respectively, where j = 1, . . . , Si . Fix some time t > 0. Again we derive results by splitting time into intervals of length . Denote by Mi j (k, ) the number of jobs present in node i j at time t that have entered node i 1 between time k and (k + 1) ; as we keep t fixed we suppress it in our notation. Because jobs are not interfering with each other in the infinite-server realm, we can decompose the transform of interest: j=1 j=1 Supposing that the arrival rate is a deterministic function of time λi (·), by conditioning on the number of arrivals in the k-th interval, j=1 m=0 ∞ e−λi (k ) (λi (k ) )m = exp λi (k )( fi (k , z) − 1) , j=1 in which pi (u) ( pi j (u), respectively) denotes the probability that the job that entered tandem i at time u has already left the tandem (is in node j , respectively) at time t . pi (u) = P ⎝ pi j (u) = P ⎝ =1 =1 Ji < t − u, =1 Recognizing a Riemann sum and letting following form: ↓ 0, we conclude that Eq. (12) takes the j=1 = exp and we consequently find j=1 = E exp 4.2 Parallel (tandem) systems Now that the tandem case has been analyzed, the next step is to put the tandem systems as described in Sect. 4.1 in parallel. We assume that there are d parallel tandems. There are different ways in which dependence between the parallel systems can be created. Two relevant models are listed below, and illustrated in Fig. 4. In case of a stochastic rate process i (·), we obtain j=1 Therefore, it holds that j=1 i (u)( fi (u, z) − 1) du . i (u)( fi (u, z) − 1) du − s i (t ) . (15) Fig. 4 Model (M1) is illustrated on the left, and Model (M2) is illustrated on the right. The rectangles represent tandem systems, which consist of an arbitrary number of nodes in series (M1) Let ≡ (·) be a d-dimensional shot-noise process ( 1, . . . , d ) where the shots in all i occur simultaneously (the shot distributions and decay rates may be different). The process i , for i = 1, . . . , d, corresponds to the arrival rate of tandem system i . Each tandem system has an arrival process, in which the Cox processes are independent given their shot-noise arrival rates. (M2) Let ≡ (·) be the shot-noise rate of a Cox process. The corresponding Poisson process generates simultaneous arrivals in all tandems. Remark 4.2 The model in which there is essentially one shot-noise process that generates arrivals for all queues independently is a special case of Model (M1). This can be seen by setting all components of = ( 1, . . . , d ) equal, by letting the shots and decay rate be identical. In Model (M1), correlation between the shot-noise arrival rates induces correlation between the numbers of jobs in the different queues. In Model (M2), correlation clearly appears because all tandem systems have the same input process. Of course, the tandem systems will not behave identically because the jobs may have different service requirements. In short, correlation across different tandems in Model (M1) is due to linked arrival rates, and correlation in Model (M2) is due to simultaneous arrival epochs. We feel that both versions are relevant, depending on the application, and hence we analyze both. Analysis of (M1)— Suppose that the dependency is of the type as in Model (M1). This means that the shots, in each component of , occur simultaneously. Recall the definition of fi as stated in Eq. (13). It holds that e−si i (t) j=1 = E exp ⎝ i=1 i (u)( fi (u, z) − 1) du − si i (t) ⎠ , where the last equality holds due to (15). Analysis of (M2)— Now suppose that the dependency in this model is of type (M2), i.e., there is one shot-noise process that generates simultaneous arrivals in the parallel tandem systems. First we assume a deterministic arrival rate function λ(·). Let Mi j (k, ) be the number of jobs present in tandem system i at node j at time t that have arrived in the system between k and (k + 1) . Note that i=1 j=1 ⎞ z Ni j (t) i j i=1 j=1 To further evaluate the right-hand side of the previous display, we observe that we can write m=0 j=1 j =1 i=1 f (u, z) := in this definition p 1,..., d ≡ p 1,..., d (u) equals the probability that a job that arrived at time u in tandem i is in node i at time t [cf. Eq. (14)]. The situation that i = Si + 1 means that the job left the tandem system; we define zi,Si +1 = 1. In a similar fashion as before, we conclude that i=1 j=1 ⎞ z Ni j (t)e−s (t) i j with f defined in Eq. (17). = E exp (u)( f (u, z) − 1) du − s (t ) , (18) Example 4.3 (Two-node parallel system) In the case of a parallel system of two infinite-server queues, f (u, z) simplifies to 1=1 2=1 f (u, z11, z21) = Remark 4.4 (Routing) Consider a feedforward network with routing. As argued in the beginning of this section, the network can be decomposed as a parallel tandem system. In case there is splitting at some point, then one decomposes the network as a parallel system, in which each tandem i receives the job with probability qi , such that qi = 1. This can be incorporated simply by adjusting the probabilities contained in fi in Eq. (16), which are given in Eq. (14), so that they include the event that the job joined the tandem under consideration. For instance, the expression for pi (u) in the left equation in (14) would become Q = i, Ji < t − u , where Q is a random variable with a generalized Bernoulli (also called ‘categorical’) distribution, where with qi = 1; the right equation in (14) is adjusted similarly. Other than that the analysis is the same for the case of splits. Remark 4.5 (Networks with loops) So far we have only considered feedforward networks. Networks with loops can be analyzed as well, but the notation becomes quite cumbersome. To show the method in by which networks with loops and routing can be analyzed, we consider a specific example. Suppose that arrivals enter node one, after which they enter node two. After they have been served in node two, they go back to node one with probability η, or leave the system with probability 1 − η. In this case, with similar techniques as before, we can find E z1N1(t)z2N2(t) = exp (u)( f (u, z1, z2) − 1) du , =1 f (u, z1, z2) = P( job(u) left system) + zi P( job(u) is in node i ), i=1 in which job(u) is the job that arrived at time u and we are examining the system at time t . Now, if we denote service times in the j -th node by J ( j), then, at a specific time t , P( job(u) left system) = k=0 i=1 Analogously, P( job(u)is in node 1) equals, by conditioning on the job having taken k loops, likewise, P( job(u)is in node 2) equals > t − u, Jk(+1)1 + For example, in the case where all Ji( j) are independent and exponentially distributed with mean 1/μ, we can calculate those probabilities explicitly. Indeed, if we denote by Y a Poisson process with rate μ, then, for example > t − u, Jk(+1)1 + k=0 k=0 i=1 i=1 = P(Y (t − u) = 2k + 1) P( job(u)is in node 2) = f (u, z1, z2) = z1 i=1 i=1 i=1 m=0 ηm e−μ(t−u) (μ(t − u))2m+1 m=0 ∞ ηm e−μ(t−u) (μ(t − u))2m+1 (2m + 1)! m=0 ηm (1 − η)F (2m+2,μ)(t − u) m=0 = z1e−μ(t−u) cosh μ√η(t − u) + z2 A similar calculation can be done for the probability that the job is in node one. Recalling that a sum of exponentials has a Gamma distribution, we can write where F (2m+2,μ) denotes the distribution function of a -distributed random variable with rate μ and shape parameter 2m + 2. 4.3 Main result In this subsection, we summarize and conclude with the following main result. Recall Definition 4.1 of a feedforward network. In the beginning of Sect. 4 we argued that we can decompose a feedforward network into parallel tandems. In Sect. 4.2 we studied exactly those systems, leading up to the following results. Theorem 4.6 Suppose we have a feedforward network of infinite-server queues, where the input process is a Poisson process with shot-noise arrival rate. Then the network can be decomposed into parallel tandem systems. In Model (M1), it holds that i (u)( fi (u, z) − 1) du − si i (t) ⎠ with fi (·, ·) as defined in Eq. (13) and where g(s, v) is a vector-valued function in which component i is given by si e−ri (t−v) − eri v ( fi (u, z) − 1)e−ri u du. i=1 Furthermore, in Model (M2), i=1 j=1 = E exp (u)( f (u, z) − 1) du − s (t ) ( f (u, z) − 1)e−ru du with f (·, ·) as defined in Eq. (17). Proof These are Eqs. (16) and (18) to which we applied Lemma 2.1. Next we calculate covariances between nodes in tandem and parallel thereafter. Covariance in Tandem System— Consider a tandem system consisting of two nodes and we want to analyze the covariance between the number of jobs in the nodes. Dropping the index of the tandem system, denote by N1(·) and N2(·) the number of jobs in nodes 1 and 2, respectively. Using Eq. (16), differentiation yields E N2(t ) = P( J1 < t − u, J1 + J2 > t − u) E E N1(t )N2(t ) = E = Cov = 2 = 2 P( J1 < t − u, J1 + J2 > t − u) (u) du P( J1 > t − v) (v) dv P( J1 < t − u, J1 + J2 > t − u) (u) du, P( J1 > t − v) (v) dv cf. Eq. (2) for the last equality. Covariance parallel (M1)— Consider a parallel system consisting of two nodes only. In order to study covariance in the parallel (M1) case, we need a result about the covariance of the corresponding shot-noise process. Lemma 4.7 Let 1(·), 2(·) be shot-noise processes of which the jumps occur simultaneously according to a Poisson arrival process with rate ν. Let the decay be exponential with rate r1, r2, respectively. Then it holds that, for δ > 0, Cov( 1(t ), 2(t + δ)) = e−r2δ Cov( 1(t ), 2(t )) = e−r2δ ν E B11 B12 (1 − e−(r1+r2)t ), r1 + r2 which, in the case 1 = 2, reduces to corresponding to [16, p. 394] and Eq. (2). Proof See Appendix 2. By making use of Eq. (16), we find Cov( i (t ), i (t + δ)) = e−ri δ Var i (t ), for i = 1, 2, E N1(t )N2(t ) = E 1(u) P( J1 > t − u) du 2(v) P( J2 > t − v) dv . Cov(N1(t), N2(t)) = Cov 1(u) P(J1 > t − u) du, 2(v) P(J2 > t − v) dv t t ν E B11 B12 1 − e−(r1+r2)v e−r2(u−v) P(J1 > t − u) P(J2 > t − v) du dv, 0 v r1 + r2 where we made use of the fact that, for u ≥ v, Cov( 1(u), 2(v)) = 1 − e−(r1+r2)v e−r2(u−v), = 2 = 2 cf. Lemma 4.7. The following proposition compares the correlations present in Model (M1) and (M2). In the proposition, we refer to the number of jobs in queue j in Model (Mi ) at time t as N (ji)(t ), for i = 1, 2. We find the anticipated result that the correlation in Model (M2) is stronger than in Model (M1). Proposition 4.8 Let (·) be the shot-noise process that generates simultaneous arrivals in both queues and let 1(·), 2(·) be processes that have simultaneous jumps and generate arrivals in both queues independently. Suppose that 1(t ) =d 2(t ) =d (t ), for t ≥ 0. Then, for any t ≥ 0, Cov(N1(t ), N2(t )) = Cov Covariance parallel (M2)— Extracting the mixed moment from the transform in Eq. (18), we derive directly that E N1(t )N2(t ) = E (u) P( J1 > t − u, J2 > t − u) du (u) P( J1 > t − u) du (u) P( J2 > t − u) du . (u) P( J1 > t − u) du, (u) P( J2 > t − u) du E (u) P( J1 > t − u, J2 > t − u) du. Proof Because of the assumption 1(t ) =d 2(t ) =d (t ), we have that, for all combinations i, j ∈ {1, 2}, the N ( j)(t ) are equal in distribution. Therefore, it is sufficient i to show that ≤ Cov The expressions for the covariances, which are derived earlier in this section, imply that − Cov = E (u) P( J1 > t − u, J2 > t − u) du, which is nonnegative, as desired. 5 Concluding remarks We have considered networks of infinite-server queues with shot-noise-driven Coxian input processes. For the single queue, we found explicit expressions for the Laplace transform of the joint distribution of the number of jobs and the driving shot-noise arrival rate, as well as a functional central limit theorem of the number of jobs in the system under a particular scaling. The results were then extended to a network context: we derived an expression for the joint transform of the numbers of jobs in the individual queues, jointly with the values of the driving shot-noise processes. We included the functional central limit theorem for the single queue, but it is anticipated that a similar setup carries over to the network context, albeit at the expense of considerably more involved notation. Our future research will include the study of the departure process of a single queue; the output stream should remain Coxian, but of another type than the input process. Acknowledgements The authors thank Marijn Jansen for insightful discussions about the functional central limit theorem in this paper. The research for this paper is partly funded by the NWO Gravitation Project NETWORKS, Grant Number 024.002.003. The research of Onno Boxma was also partly funded by the Belgian Government, via the IAP Bestcom Project. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix 1: Proof of Lemma 2.1 There are various ways to prove this result; we here include a procedure that intensively relies on the probabilistic properties of the shot-noise process involved. Observe that, recognizing a Riemann sum, k=1 With PB (t ) a Poisson process with rate ν and the Ui i.i.d. samples from a uniform distribution on [0, 1], it holds that (k ) = =1 i=PB (( −1) )+1 Bi e−r Ui e−r(k− ) . We thus obtain that the expression in (20) equals (where the equality follows by interchanging the order of the summations) k=1 =1 i=PB (( −1) )+1 k= Bi e−r Ui e−r(k− ) Bi e−r Ui f (k , z)e−r(k− ) , which behaves as Furthermore, we have the representation We conclude that E z N (t)e−s (t) equals =1 i=PB(( −1) )+1 Bi e−r Ui er =1 f (u, z)e−ru du i=PB (( −1) )+1 Bi e−r Ui . =1 i=P(( −1) )+1 Bi e−r(t− Conditioning on the values of PB ( ) − PB (( − 1) ), for = 1, . . . , t / , and using that the Bi are i.i.d., we find that the expression in Eq. (21)) is equal to the limit, as ↓ 0, of the following expression =1 =1 k =0 which, in the limit ↓ 0, equals ⎛ e−ν exp ⎜ ν ⎝ which can be written as =1 The lemma now follows from continuity of the exponent and the definition of the Riemann integral. Appendix 2: Proof of Lemma 4.7 Let PB (·) be the Poisson process with rate ν, corresponding to the occurrences of shots, and let Et,δ(n) be the event that PB (t + δ) − PB (t ) = n. By conditioning on the number of shots in the interval (t, t + δ], we find We proceed by rewriting the conditional expectation as i=1 1 E( 1(t) 2(t + δ) | Et,δ(n)) = δn denoting by Ft1,...,tn,δ(n) the event Et,δ(n) and the arrival epochs are t1, . . . , tn. Note that we have, due to Eq. (1), conditional on Ft1,...,tn,δ(n), the distributional equality f (u, z)e−ru du − se−r(t− + Ui ) f (u, z)e−ru du − se−r(t− + Ui ) − e−r Ui er f (u, z)e−ru du n=0 and consequently = E E Bi2e−r2(t+δ−ti ). i=1 Note that for all i = 1, . . . , n we have After unconditioning Eq. (23) with respect to the arrival epochs by integrating over all ti from t to t + δ and dividing by δn, we thus obtain E( 1(t) 2(t + δ) | Et,δ(n)) = E and hence, denoting i := limt→∞ i (t ) for i = 1, 2, = E 1(t) E 2 + e−r2δ E 1(t) 2(t) − E 1(t) E 2 , where we made use of E i = ν E B1i /ri . It follows that Cov( 1(t ), 2(t + δ)) = E 1(t ) 2(t + δ) − E 1(t ) E 2(t + δ) 1(t ) 2(t )e−r2δ + (1 − e−r2δ) E 1(t ) E 1(t ) E 2(t )e−r2δ + (1 − e−r2δ) E where the equality E 2(t + δ) = E 2(t )e−r2δ + (1 − e−r2δ) E 2 is used, which can be directly checked using the expressions for the mean in Eq. (2). This proves the first equality in Eq. (19). The proof of the second equality follows from Cov( 1(t ), 2(t )) = E 1(t ) 2(t ) − E 1(t ) E 2(t ), in which n=0 e−r1(t−u) du e−r2(t−v) dv ν2 E B11 E B12(1 − e−r1t )(1 − e−r2t ) (1 − e−(r1+r2)t ), i (t ), for i = 1, 2, is given in Eq. (2). 1. Anderson , D. , Blom , J. , Mandjes , M. , Thorsdottir , H. , de Turck , K. : A functional central limit theorem for a Markov-modulated infinite-server queue . Method. Comput. Appl. Probab . 18 , 153 - 168 ( 2016 ) 2. Bar-Lev , S. , Boxma , O. , Mathijsen , B. , Perry , D. : A blood bank model with perishable blood and impatience . Eurandom reports ( 2015 ) 3. Centanni , S. , Minozzo , M. : A Monte Carlo approach to filtering for a class of marked doubly stochastic Poisson processes . J. Am. Stat. Assoc . 101 (476), 1582 - 1597 ( 2006 ) 4. Csörgo˝ , M. , Révész , P. : Strong Approximations in Probability and Statistics. Academic Press, New York ( 1981 ) 5. Dassios , A. , Jang , J.-W.: Pricing of catastrophe reinsurance and derivatives using the Cox process with shot noise intensity . J. Financ. Stoch . 7 , 73 - 95 ( 2003 ) 6. Eick , S. , Massey , W. , Whitt , W. : The physics of the Mt /G/∞ queue. Manag. Sci. 39 ( 2 ), 241 - 252 ( 1993 ) 7. Kurtz , T. , Ethier , S. : Markov Processes: Characterization and Convergence . Wiley, New York ( 1986 ) 8. Gao , X. , Zhu , L.: A functional central limit theorem for stationary Hawkes processes and its application to infinite-server queues . ( 2016 ). 9. Hawkes , A.G. : Point spectra of some mutually exciting point processes . J. R. Stat. Soc . 33 , 438 - 443 ( 1971 ) 10. Iksanov , A. , Jurek , Z. : Shot noise distributions and selfdecomposability . Stoch. Anal. Appl . 21 ( 3 ), 593 - 609 ( 2003 ) 11. Kim , S.-H. , Whitt , W.: Are call center and hospital arrivals well modeled by nonhomogenous Poisson processes? Manuf . Serv. Oper. Manag. 16 ( 3 ), 464 - 480 ( 2014 ) 12. Lewis , P.A. : Simulation of non-homogenous Poisson processes by thinning . Naval Res . Logist . Q. 26 ( 3 ), 403 - 413 ( 1979 ) 13. Ogata , Y. : Statistical models for earthquake occurences and residual analysis for point processes . J. Am. Stat. Assoc . 83 , 9 - 27 ( 1988 ) 14. Pang , G. , Talreja , R. , Whitt , W.: Martingale proofs of many-server heavy-traffic limits for Markovian queues . Probab. Surv . 4 , 193 - 267 ( 2007 ) 15. Pang , G. , Whitt , W.: Two-parameter heavy-traffic limits for infinite-server queues . Queueing Syst . 65 , 325 - 364 ( 2010 ) 16. Ross , S.M. : Stochastic Processes . Wiley, New York ( 1996 ) 17. Whitt , W.: Proofs of the Martingale FCLT. Probab. Surv . 4 , 268 - 302 ( 2007 )

This is a preview of a remote PDF:

D. T. Koops, O. J. Boxma, M. R. H. Mandjes. Networks of \(\cdot /G/\infty \) queues with shot-noise-driven arrival intensities, Queueing Systems, 2017, 1-25, DOI: 10.1007/s11134-017-9520-7