#### 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 (http://creativecommons.org/licenses/by/4.0/), 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 ). https://arxiv.org/abs/1607.06624
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 )