Expansions for Steady-State Characteristics of (max,+)–Linear Systems

François Baccelli^{1}  Sven Hasenfuss^{2}  Volker Schmidt^{2}

1 INRIA Sophia Antipolis, 2004, Route des Lucioles, B.P. 93, F-06902 Sophia Antipolis Cedex, France

            E-mail: baccelli@sophia.inria.fr

2 Universität Ulm, Abteilung Stochastik, Helmholtzstr. 18, D-89069 Ulm, Germany

            E-mail: hasenfus@mathematik.uni-ulm.de
                  schmidt@mathematik.uni-ulm.de

This paper on computational probability was written in honour to Marcel Neuts, who contributed so centrally to this field.

Abstract

An expansion formula is given for the expected value of a general class of functions of steady state variables in open, stochastic (max,+)–linear systems with Poisson input. Expansions for Laplace transforms, moments, distribution functions and tail probabilities of these steady state variables are considered as specific instances of our general expansion formula. Such (max,+)–linear systems are known to allow one to represent a class of discrete event networks called stochastic event graphs. Examples of such event graphs pertaining to queueing theory and Petri net theory are given in the paper in order to illustrate the proposed expansion method, including at the numerical level.

Keywords: Queueing networks, stochastic Petri nets, Poisson input, vectorial recurrence equation, stationary state variables, marked point processes, factorial moment measures, perturbation analysis, admissibility.

AMS classification: 60K25, 60G55, 90B22

1 Introduction

The present paper bears on a general expansion technique which can be applied to the stationary characteristics of two different though related classes of stochastic discrete event systems with Poisson input:

1. A subclass of queueing networks. Within the framework of queueing theory, the expansion technique can be applied to the subclass of open queueing networks which admit a representation in terms of an evolution equation involving only the max and + operations. This class contains of course the elementary M/G/1/∞ queue, which is known to be reducible to the Lindley equation indeed. Queues in tandem, acyclic Fork-Join networks, acyclic networks with finite capacity and blocking, synchronized queueing networks provide other higher dimension (and less trivial) examples of such systems. Several examples of such queueing networks will be considered in Section 4 below. It is unfortunately impossible to give an extensive list of the queueing networks with such a property, which explains the interest of the second class below.

2. A subclass of stochastic Petri networks. Within this setting, a more compact definition can be given: the expansion technique applies to all open stochastic Petri networks which fall in the class of stochastic event graphs. All the above queueing networks can be seen as special incarnations of such stochastic event graphs. Other and more intricate systems such as the Kanban example considered in Section 4 below fall in this class too. Roughly speaking, such a system can be reduced to a sequence of random matrices in the (max, +) algebra. There are as many incarnations of such systems as there are ways of choosing the dimension and the entries of these matrices.

Expansions for this type of systems were first derived for the stationary first moment of waiting times in [4]. The main innovation of the present paper lies in the extension of the method to more complex stationary characteristics than the first moment, like for instance stationary Laplace transforms, distribution functions, tail probabilities, higher moments, etc.

The new mathematical result of the paper is the expansion formula of Theorem 1, which gives the general expression of the coefficients in the expansion of an arbitrary function of the steady state variables of a (max,+)–linear system. Because of space constraints, we could not give a complete proof of Theorem 1. The following compromise was adopted: all new steps required for this proof are given in the present paper; all steps which are similar to those used in the proof of the corresponding theorem in [4] are only summarized. We nevertheless give a reference to the internal report [2] which is available on our web pages, and which contains the complete proof of all steps.

In place of a complete proof, we give several examples pertaining to both queueing theory and Petri net theory, and propose in Section 4 a detailled numerical exploitation of the expansion formula for various types of stationary characteristics. We do not know yet under what conditions such expansions are analytic expansions, and we would like to stress that this is quite an intriguing open question at this stage. Nevertheless, we hope that the generality of the formula and the quality of the numerical results will be sufficient to trigger the interest of both practically and theoretically oriented readers.

It is worthwhile noticing that the expansion method can be applied to non-Markovian settings, and that it is not restricted to problems of small dimensions. In particular, the case with several unbounded coordinates poses no specific difficulties.

2 General Setting

This section gives a brief summary of the mathematical setting for (max, +)-linear systems. For a more complete exposition on this or for more potential examples, the reader should refer to [1] and [4] and the literature therein.

2.1 Summary of Setting

The present paper focuses on the unique solution W (·) ∈ IR^{α} of a vectorial (max, +)–linear functional equation of the form

W (θω) = A(ω) ⊗ C(ω) ⊗ W (ω) ⊕ B(θω),  (1)

which arises in the theory of open multidimensional discrete event systems. In this equation,

• θ is a shift operator defined on some probability space (Ω, F, IP), with generic element ω, and such that θ is IP-invariant and IP-ergodic;

• A(·), B(·) and C(·) are random matrices which characterize the topological and the stochastic structure of the system: A(·) and B(·) are associated with the system itself, and C(·) with its input point process;

• The reference algebra is the so-called (max, +) algebra: when operating on scalars, ⊕ stands for max and ⊗ for +, whereas the ⊗-product of two matrices, say a of size p × q and b of size q × r is the p × r matrix a ⊗ b with entries

(a ⊗ b)_{i,j} = ⨁_{k=1}^{q} a_{i,k} ⊗ b_{k,j}

and the ⊕-sum of two matrices, say a and b both of size p × q is the p × q matrix a ⊕ b with entries

(a ⊕ b)_{i,j} = a_{i,j} ⊕ b_{i,j}.

We will concentrate on the stable case, i.e. when the unique solution of (1) is finite (the conditions for stability are given in the above references).

As it was mentioned in the introduction, vectorial equations of the form given in (1) typically arise when looking at the temporal evolution of open, single input, stochastic event graphs. In such systems, the epochs X_{n}^{i} at which events occur on node i (for example, certain arrivals or departures in a queueing network) satisfy a vectorial recurrence equation of the form

X_{n+1} = A_{n} ⊗ X_{n} ⊕ B_{n+1} ⊗ T_{n+1}  (2)

with X_{n} = (X_{n}^{1}, …, X_{n}^{α})^{⊤}, where {T_{n}} is the sequence of IR-valued arrival times and where the random matrices A_{n} ∈ IR^{α,α} and B_{n} ∈ IR^{α,1} characterize the topological and the stochastic structure of the system (for example, service times and routing structure in a queueing network).

It is easy to check that the vector W_{n} = (W_{n}^{1}, …, W_{n}^{α})^{⊤} with W_{n}^{i} = X_{n}^{i} − T_{n} (which can be interpreted as a waiting time vector in the queueing context) then satisfies a vectorial recurrence equation of the form

W_{n+1} = A_{n} ⊗ C_{n} ⊗ W_{n} ⊕ B_{n+1},  (3)

where C_{n} = Diag(−τ_{n}), is the diagonal α × α (max, +) matrix with all its diagonal entries equal to −τ_{n} = T_{n} − T_{n+1}, and all its non-diagonal entries equal to ε = −∞.

Under stationarity and ergodicity conditions on the firing times in the underlying event graph (service times and inter-arrival times in the queueing network), which will contain the case with independent and i.i.d. firing (service) times, the sequences A_{n}, B_{n} and C_{n} are stationary and ergodic, which means that they can be represented as A(θ^{n}ω), B(θ^{n}ω) and C(θ^{n}ω) respectively, where θ is a mere translation on a space of sequences. As it was shown in [1], under the above assumptions, in the stable case, {W_{n}} couples with a unique stationary sequence {W (θ^{n})}, where W (·) is the unique finite random vector solution of (1), which is given by the series

W = B ⊕ ⨁_{n≥1} A(θ^{−1}) ⊗ C(θ^{−1}) ⊗ … ⊗ A(θ^{−n}) ⊗ C(θ^{−n}) ⊗ B(θ^{−n}).  (4)

Let D_{0}, D_{1}, … be the following α-dimensional vectors: D_{0} = B(θ^{−1}), and for k ≥ 1

D_{k} = (⨂_{n=1}^{k} A(θ^{−n})) ⊗ B(θ^{−k}).  (5)

Since C commutes with any matrix, W admits the following equivalent representation:

W = B_{0} ⊕ ⨁_{n≥1} Diag(T_{−n}) ⊗ D_{n}.  (6)

Various incarnations of such equations will be considered in Section 4.

Let IR^{+} = [0, ∞). The aim of the present paper is to show that for all f : IR^{+} → IR^{+} in a general class of functions, under certain assumptions, whenever the input process is Poisson, the expectation IEf(W^{i}), where W^{i} is the i-th component of the stationary state vector W, can be expanded into a finite power series with respect to the arrival intensity λ (Theorem 1). Formulas for the coefficients and algorithms to compute them are the main practical output of the paper. This result can be specialized to derive series expansions for

1. Laplace transforms IE (e^{−sW^{i}}) , s ≥ 0;

2. Higher-order moments IE (W^{i})^{ν}, ν ∈ IN = {1, 2, …};

3. Tail probabilities IP (W^{i} > x), x ≥ 0.

2.2 Technical Conditions for the Expansions

The support and stochastic assumptions of the forthcoming paragraphs will in particular be satisfied for any single input stochastic event graph with i.i.d. and integrable firing times, provided it is input-connected (see the above references) and stable.

Support Assumptions We assume that

1. each entry of the α × α-dimensional random matrix A is either a.s. non-negative or a.s. equal to ε = −∞, i.e.

A_{i,j} ≥ 0 or A_{i,j} = ε,  a.s.  (7)

and all entries on the diagonal of A are a.s. non-negative.

2. there exists an integer 0 < α′ ≤ α such that the first α′ coordinates of the α-dimensional random vectors B are a.s. non-negative, i.e. B^{i} ≥ 0, for all 1 ≤ i ≤ α′.

3. the α′ first coordinates of D_{k} are non-decreasing in k:

0 ≤ D_{0}^{i} ≤ D_{1}^{i} ≤ …,  (8)

for all i = 1, …, α′.

Stochastic Assumptions Throughout the paper, we assume that T = {T_{n}} is a homogeneous Poisson process with intensity λ, and that Z = {Z_{n}} where Z_{n} = (A_{n}, B_{n}) = (A(θ^{n}), B(θ^{n})) is integrable and independent of T. Since the sequences T and Z are independent, we can also assume that T_{n+1}(ω) − T_{n}(ω) = τ(θ^{n}ω) for some random variable τ(·) > 0. For all n ∈ IN_{0} = IN ∪ {0} and ℓ ∈ IN, let

G_{n} = ⨁_{i=1}^{α} {(A_{−1} ⊗ A_{−2} ⊗ … ⊗ A_{−n} ⊗ (B_{−n} ⊕ O))^{i}} ,  (9)

and

H_{n,ℓ} = ⨁_{i=1}^{α′} {(A_{−(nℓ+1)} ⊗ A_{−(nℓ+2)} ⊗ … ⊗ A_{−(n+1)ℓ} ⊗ (B_{−(n+1)ℓ} ⊕ O))^{i}} ,  (10)

where O is the α-dimensional column vector with all its components equal to zero. We assume that with probability one

lim_{ℓ→∞} H_{0,ℓ} = ∞.  (11)

Let a denote the maximal Lyapunov exponent of evolution equation (3), i.e.

a = lim_{r→∞} \frac{1}{r} \{(⨂_{k=0}^{r} A_{−k})_{i,j}\} .  (12)

Assume that λ < a^{−1}. Then the system is stable. Furthermore, for some r′ ∈ IN large enough, for all s ≥ r′,

λ < s (IE H_{0,s})^{−1}.  (13)

A key assumption that we make is that {H_{n,r′}}_{n≥0} is a sequence of 1-dependent random variables. In the event graph case, this is satisfied whenever the sequences of firing times are i.i.d., be it in the independent case or in the dependent case (like for instance in the case when the n-th service time on all nodes of a queueing network is the very same sample of an i.i.d. sequence, as it is the case for the tandem queue example considered in Section 4.3).

An immediate consequence of the above stochastic assumptions is that there exists an r ∈ IN, r ≥ r′ such that

λ < (r − 1) (IE H_{0,r})^{−1}  (14)

and {H_{n,r}}_{n≥0} is a sequence of 1-dependent random variables. In what follows, we will most often drop the index r and will only write H_{n} instead of H_{n,r}, but will always implicitly assume that r is chosen as above.

2.3 A Family of Polynomials

We now recall a class of polynomials which were introduced for the expansion of IEW in [4]. These polynomials are given by the following formula

p_{k}(x_{0}, x_{1}, …, x_{k−1}) = ∑_{(i_{0},i_{1},...,i_{k−1})∈S_{k}}(−1)^{γ_{k}(i_{0},i_{1},...,i_{k−1})} \frac{x_{0}^{i_{0}}}{i_{0}!} \frac{x_{1}^{i_{1}}}{i_{1}!} … \frac{x_{k−1}^{i_{k−1}}}{i_{k−1}!} ,  (15)

where

S_{k} = {(i_{0}, i_{1}, …, i_{k−1}) ∈ {0, 1, …, k}^{k} : i_{0} + i_{1} + … + i_{k−1} = k and if
i_{s} = l > 1, then i_{s−1 mod k} = i_{s−2 mod k} = … = i_{s−l+1 mod k} = 0}, (16)

and γ_{k}(i_{0}, i_{1}, …, i_{k−1}) = 1 + ∑_{s=0}^{k−1} 1(i_{s} > 0), for all k ≥ 1, with 1(i_{s} > 0) = 1 whenever i_{s} > 0 and 1(i_{s} > 0) = 0 otherwise.

The polynomials p_{k} satisfy the following properties (see [4]):

• Translation Invariance:

p_{k}(x_{0} + u, x_{1} + u, …, x_{k−1} + u) = p_{k}(x_{0}, x_{1}, …, x_{k−1})  (17)

for all non-decreasing sequences of real numbers 0 ≤ x_{0} ≤ x_{1} ≤ … ≤ x_{k−1}, for all u ≥ 0 and k = 2, 3, . . ..

• Integral Representation:

p_{k}(x_{0}, …, x_{k−1}) = ∑_{n=0}^{k−2} ∫_{x_{n}}^{x_{n+1}} \{p_{k−1}(\underbrace{v, …, v}_{n}, x_{n+1}, …, x_{k−1})

− p_{k−1}(\underbrace{v, …, v}_{n+1}, x_{n+1}, …, x_{k−2})\} dv,  (18)

for all non-decreasing sequences of real numbers 0 ≤ x_{0} ≤ x_{1} ≤ … ≤ x_{k−1}, k ≥ 2,

with p_{1}(x_{0}) = x_{0}.

For further interesting properties of these polynomials, see also [3].

3 Expansion Formula

In this section, we suppose that all the assumptions made on {(A_{n}, B_{n})} in Section 2 hold. Let f : IR^{+} → IR^{+} be some Riemann-integrable function (e.g. f(x) monotone or continuous). Associated with f, we define the following class of functions: q_{1}(x_{0}) = f(x_{0}), and

q_{k+1}(x_{0}, x_{1}, …, x_{k})
= ∑_{n=0}^{k} \binom{k}{n}(−1)^{k−n}F^{[k]}(x_{n}) − ∑_{n=0}^{k−1} ∑_{j=n}^{k−1} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{n})
{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})},  (19)

for 1 ≤ k ≤ m, with F^{[0]}(x) = f(x) and F^{[n]}(x), n ≥ 1, is recursively defined by a suitably chosen version of the indefinite Riemann-integral ∫ F^{[n−1]}(x) dx, x ∈ IR^{+}. The functions p_{k} in (19) are the polynomials defined in (15).

Theorem 1 Let f : IR^{+} → IR^{+} be a Riemann-integrable, non-negative function which satisfies f(x) ≤ c x^{ν} for all x ≥ 0, where ν ∈ IN_{0} is a non-negative integer and c ∈ IR^{+} is a finite constant.

1. If the (ν + m + 2)-th moment of H_{n} is finite, (where H_{n} was defined in (10)), then, for 1 ≤ i ≤ α′,

IE[f(W^{i})] = ∑_{k=0}^{m} λ^{k} IE[ q_{k+1}(D_{0}^{i}, D_{1}^{i}, …, D_{k}^{i}) ] + O(λ^{m+1}).  (20)

2. Assume that there exists a θ^{∗} > 0, such that for all θ ∈ [0, θ^{∗}),

IE[ e^{θG_{n}}] ≤ L_{θ} (ϕ(θ))^{n} ,  (21)

(where G_{n} was defined in (9)), for some finite functions ϕ(θ) > 1 and L_{θ}. Then, for 1 ≤ i ≤ α′ the expectation IE[f(W^{i})] is infinitely differentiable in λ for all λ in a right neighborhood of 0, expansion (20) holds for all m, and

lim_{λ↓0} \frac{d^{n}}{dλ^{n}} IE[f(W^{i})](λ) = IE[ q_{n+1}(D_{0}^{i}, D_{1}^{i}, …, D_{n}^{i}) ] .  (22)

For instance, in order to derive a finite series expansion of the Laplace transform of the components W^{i} of the stationary state vector W, we choose f(x) = e^{−sx}, for all x ∈ IR^{+} and s ≥ 0. Since f(x) is bounded by 1 for all x ∈ IR^{+}, Theorem 1 can be applied to obtain an expansion of IE e^{−sW^{i}}. To explicitly compute the coefficients of this expansion, we choose F^{[n]}(x) = ∫_{x}^{∞} F^{[n−1]}(u) du. In principle, it is possible to consider any version of the indefinite integral ∫ F^{[n−1]}(x) dx for F^{[n]}(x) in Theorem 1. However, depending on the particular function f(x), a suitably chosen version of the indefinite integral may help to reduce computational efforts significantly.

To give another example, taking f(x) = x gives a function which satisfies the boundedness requirement on f(x) with c = 1 and ν = 1. For this function, it is advantageous to take F^{[n]}(x) = ∫_{0}^{x} F^{[n−1]}(u) du, i.e. F^{[n]}(x) = \frac{x^{n+1}}{(n+1)!}, for n ∈ IN_{0}. Thus,

q_{k+1}(x_{0}, x_{1}, …, x_{k}) = ∑_{n=0}^{k} \binom{k}{n}(−1)^{k−n} \frac{(x_{n})^{k+1}}{(k + 1)!} − ∑_{n=0}^{k−1} ∑_{j=n}^{k−1} \binom{j}{n}(−1)^{j−n}

·\frac{(x_{n})^{j+1}}{(j + 1)!} \{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})\}.  (23)

Elementary but tedious calculations show that the polynomials obtained in (23) are equal to those given by (15), see Section 6.2 of [4], i.e. Theorem 1 is actually a generalization of the result obtained in [4] and (23) gives another recurrence equation for the polynomials p_{k}.

4 Examples

In essence, the examples below have already been considered in [4]. It was shown there that the technical conditions of Section 2 are satisfied, and this was used to derive expansions for the first moment of W^{i} in some particular cases. Whenever the integrability conditions of Theorem 1 of the present paper are satisfied, we can use general expansion formula (20) to derive expansions for further characteristics.

4.1 Queues in Tandem

Consider a network of κ single-server FIFO queues with infinite capacity in tandem and all queues initially empty. Let σ_{n}^{i} denote the service time of customer n in station i. The i-th component W^{i} of the random vector W describes the stationary waiting time of a randomly chosen customer until the beginning of his service on station i. It was shown in [4] that in this case

D_{n}^{i} = max_{1≤l_{n}≤...≤l_{1}≤i} { ∑_{k=l_{1}}^{i−1} σ_{0}^{k} + ∑_{k=l_{2}}^{l_{1}} σ_{−1}^{k} + … + ∑_{k=l_{n}}^{l_{n−1}} σ_{−n+1}^{k} + ∑_{k=1}^{l_{n}} σ_{−n}^{k}} .  (24)

In this section let us consider the distribution function of W^{i}. In other words, we want to find the distribution of the stationary waiting time from the arrival of a randomly chosen customer until the beginning of his service at station i, or, more formally, we want to find IP(W^{i} ≤ ξ) = 1 − IP(W^{i} > ξ), for ξ ≥ 0. We apply Theorem 1 to derive a finite series expansion of IP (W^{i} > ξ). Before we do so, we notice that IP (W^{i} > ξ) = IE 1(W^{i} > ξ) and thus choose f(x) = 1(x > ξ) for all x ≥ 0. Clearly, 0 ≤ f(x) ≤ 1, for all x ∈ IR^{+}. For this functional f(x) we use the version F^{[n]}(x) = \frac{1}{n!}(x − ξ)^{n} 1(x > ξ), for all n ∈ IN_{0}, of the indefinite integral ∫ F^{[n−1]}(x) dx, when applying Theorem 1.

Expansions — Stochastic Service Times Case To give a first numerical example, we take i = 2. Due to the acyclic nature of this queueing network whatever happens in queue 3 or any other queue placed behind queue 2 has no influence on W^{2}. Thus, without loss of generality, we can simply put κ = 2. From (24) we see that in this case

D_{0}^{2} = σ_{0}^{1},
D_{1}^{2} = max {σ_{−1}^{1} + σ_{0}^{1}, σ_{−1}^{1} + σ_{−1}^{2}},
D_{2}^{2} = max {σ_{−2}^{1} + σ_{−1}^{1} + σ_{0}^{1}, σ_{−2}^{1} + σ_{−1}^{1} + σ_{−1}^{2}, σ_{−2}^{1} + σ_{−2}^{2} + σ_{−1}^{2}}.

Let us assume that for each i = 1, 2, the service times σ_{n}^{i} of server i are given by a sequence of independent and identically distributed nonnegative random variables. In particular, let service times σ_{n}^{1} of server 1 follow a 2-Erlang(\frac{3}{4}) distribution, whereas σ_{n}^{2} is distributed according to a 3-Erlang(\frac{2}{3}) distribution. Thus, the mean service times of server 1 and 2 are 1.5 and 2, respectively. Then Theorem 1 yields

c_{0} = IE q_{1}(D_{0}^{2}) = IE 1(σ_{0}^{1} > ξ) = \frac{1}{3} (3 + 4 ξ)e^{− \frac{4}{3} ξ}

for the first coefficient of expansion (20). The coefficient of λ is found to be

c_{1} = IE q_{2}(D_{0}^{2}, D_{1}^{2}) = IE[(D_{1}^{2} − ξ)(1(D_{1}^{2} > ξ) − 1(D_{0}^{2} > ξ))]
= e^{− \frac{3}{2} ξ}(48 ξ^{2} + 1280 ξ + 12032) + e^{− \frac{4}{3} ξ}(\frac{8}{27} ξ^{3} + \frac{4}{3} ξ^{2} + \frac{19 456}{27} ξ − \frac{2 921 216}{243})
− e^{− \frac{17}{6} ξ}(\frac{8}{9} ξ^{3} + \frac{128}{27} ξ^{2} + \frac{896}{81} ξ + \frac{2 560}{243}).

Similar computations lead to higher-order coefficients. Note that these coefficients are all functions of ξ. Thus, for any fixed arrival intensity λ < 1/a, some m ∈ IN_{0} and for ξ ≥ 0, we get that

Θ_{m}(ξ) = 1 − (c_{0}(ξ) + c_{1}(ξ) λ + … + c_{m}(ξ) λ^{m})

is an approximation for IP(W^{2} ≤ ξ), the distribution function of W^{2} under the assumption of a fixed arrival intensity λ. Figure 1 shows the plot of Θ_{0}(ξ) as well as that of Θ_{1}(ξ) and Θ_{2}(ξ) together with an empirical distribution function of the stationary waiting time W^{2} which was obtained by a simulation study. The arrival intensity λ was chosen to be λ = 0.12 in this example. Since a = 2.5 for the given queueing network, this means that the traffic intensity ρ = λ a is equal to 0.3. For more details on how to compute a and thus ρ see [1]. It can be seen from Figure 1 that the second order approximation Θ_{2}(ξ) already matches the simulated curve very well.

Figure 1: Approximations Θ_{m}(ξ) of IP(W^{2} ≤ ξ) for ρ = 0.3.

Expansions — Deterministic Service Times Case In the next example we focus on the distribution of W^{3}, the stationary waiting time from arrival until the beginning of service at station 3 in a three station tandem queueing network. Another difference is that the service times σ^{i} of server i ∈ {1, 2, 3} are assumed to be deterministic. Without loss of generality we can and will assume that σ^{1} ≤ σ^{2} ≤ σ^{3}. In any other case, say σ^{i} > σ^{i+1} for some i < 3, we can consider the i-th queue and the (i + 1)-th queue as one single-server queue with service time σ^{i} + σ^{i+1} because in front of the (i + 1)-th server the waiting room is always empty. From (24) it can be seen that if we have deterministic service at all stations, D_{n}^{3} boils down to D_{n}^{3} = σ^{1} + σ^{2} + n σ^{3} for all n ≥ 0. To give a numerical example, we choose σ^{1} = 1.5, σ^{2} = 2 and σ^{3} = 2.5. This yields D_{n}^{3} = 3.5 + 2.5 n, for all n ≥ 0. Since all D_{n}^{3} are deterministic, taking the expectation of q_{k+1}(D_{0}^{3}, …, D_{k}^{3}) reduces to plugging in corresponding values of D_{n}^{3}. The first coefficients of the expansion are readily obtained:

c_{0} = q_{1}(3.5) = 1(ξ < 3.5),
c_{1} = q_{2}(3.5, 6) = (6 − ξ)(1(ξ < 6) − 1(ξ < 3.5)),
c_{2} = q_{2}(3.5, 6, 8.5) = \frac{1}{2}(8.5 − ξ)^{2}(1(ξ < 8.5) − 1(ξ < 3.5))
− (8.5 − ξ)(6 − ξ)(1(ξ < 6) − 1(ξ < 3.5)).

In Figure 2 we show how an empirical distribution function of W^{3} is approximated by those functions derived from expansion formula (20) of Theorem 1. For this plot the arrival intensity λ has been fixed at λ = 0.2, which corresponds to a traffic intensity ρ = 0.5.

Figure 2: Approximations Θ_{m}(ξ) of IP(W^{3} ≤ ξ) for ρ = 0.5.

4.2 A Two-Stage Kanban System

In this section we apply the obtained results to analyze a two-stage Kanban system. A Petri net representation of this system is given in Figure 3. Anybody not familiar with the Petri net modeling formalism can find an excellent overview of the theory of Petri nets in [9]. This two-stage Kanban example has already been studied in Section 4.2.4 of [4], where polynomial approximations for mean stationary waiting times are calculated. In case all service times are assumed deterministic an explicit representation for IEW^{i} is given in Section 3.1 of [3]. We want to apply Theorem 1 to investigate W^{4}, i.e. the stationary waiting time starting from the arrival of a customer to the system (or, equivalently, the firing of transition u), until this customer leaves Station 1 (firing of transition t_{4}). It was shown in [4] that

D_{0}^{4} = σ_{0}^{1},
D_{1}^{4} = σ_{−1}^{1} + σ_{0}^{1},
D_{2}^{4} = max {σ_{−2}^{1} + σ_{−1}^{1} + σ_{0}^{1}, σ_{−2}^{1} + σ_{−2}^{2}},

for this Kanban system.

Expectation We first recall the results concerning the expectation of W^{4} for a specific situation. As a numerical example we assume that service time σ_{n}^{1} of the machine at Station 1 is constant and equal to 3 units of time. Similar, service time σ_{n}^{2} of machine 2 at the second station is assumed constant and equal to 5 units of time. Then we learn that D_{0}^{4} = 3,

[figure: σ_{n}^{1}; σ_{n}^{2}; ; W^{4}; u; t_{1}; t_{2}; t_{3}; t_{4}; t_{5}; t_{6}; t_{7}; Station 1; Station 2]

Figure 3: Petri net representation of a two-stage Kanban system

D_{1}^{4} = 6, D_{2}^{4} = 9 and D_{n}^{4} = 13 + 5 (n − 3), for n ≥ 3. From [3], we know that for λ ∈ [0, 0.2),

IEW^{4} = ∑_{k=0}^{∞} λ^{k} p_{k+1}(D_{0}^{4}, D_{1}^{4}, …, D_{k}^{4})  (25)

= \frac{25 λ}{2 − 10 λ} + 13 − \frac{3}{λ}

+ \frac{(1 − 5 λ) (e^{5 λ} + e^{3 λ} + e^{λ} − 3 λ − λ e^{−2 λ} − \frac{5}{2} λ^{2} e^{−5 λ})}{λ}

= 3 + \frac{9}{2} λ + \frac{27}{2} λ^{2} + \frac{3 169}{24} λ^{3} + \frac{98 657}{60} λ^{4} + \frac{2 654 251}{360} λ^{5} + O(λ^{6}).

Figure 4 shows the exact solution of IEW^{4} as a function of λ as well as some polynomial approximations, which can also be derived from Theorem 1 when choosing f(x) = x. In this example the maximal Lyapunov exponent a can be shown to be a = 5, and thus, Figure 4 describes the situation for all traffic intensities 0 ≤ ρ ≤ 0.75. We can derive similar expansions for higher-order moments of W^{4}.

Variance Choosing f(x) = x^{2} and setting c = 1 allows us to apply Theorem 1 to get the following expansion for the second moment of W^{4}:

IE(W^{4})^{2} = ∑_{k=0}^{m} λ^{k} q_{k+1}(D_{0}^{4}, D_{1}^{4}, …, D_{k}^{4}) + O(λ^{m+1}),  (26)

for any fixed m ∈ IN, where the coefficient polynomials are given by

q_{1}(x_{0}) = x_{0}^{2}
q_{2}(x_{0}, x_{1}) = \frac{1}{3} (2x_{0}^{3} + x_{1}^{3} − 3x_{0}^{2}x_{1})
q_{3}(x_{0}, x_{1}, x2) = \frac{1}{12} (3x_{0}^{4} + 2x_{1}^{4} + x_{2}^{4} − 8x_{0}^{3}x_{1} − 4x_{1}^{3}x_{2} − 6x_{0}^{2}x_{2}^{2} + 12x_{0}^{2}x_{1}x_{2}) .

Figure 4: Approximations of the mean stationary waiting time IEW^{4}.

In this representation the upper indices denote powers. Additional coefficient polynomials can be obtained analogously. Using (25) and (26) we get an expansion of the variance of W^{4}:

Var (W^{4}) = IE (W^{4})^{2} − (IEW^{4})^{2}
= ∑_{k=0}^{m} λ^{k} q̃_{k+1}(D_{0}^{4}, D_{1}^{4}, …, D_{k}^{4}) + O(λ^{m+1}),  (27)

with

q̃_{k+1}(…) = q_{k+1}(…) − ∑_{j=0}^{k} p_{j+1}(…) p_{k+1−j}(…).  (28)

The arguments of the polynomials p_{j+1}(…), q_{j+1}(…) or q̃_{j+1}(…) are always given by the sequence D_{0}^{4}, …, D_{j}^{4}. Of course, these consideration are not restricted to the 4-th component W^{4} of the stationary state vector W. Furthermore, these considerations can easily be adapted to the stochastic case. Namely, if the moment conditions of Theorem 1 are satisfied such that IEW^{i} admits an expansion of order m, then it is also possible to give an expansion of IE(W^{i})^{2} of order m − 1. Thus, by the same reasoning as above we get

Var (W^{i}) = ∑_{k=0}^{m−1} λ^{k} q_{k+1}^{⋆}(D_{0}^{i}, D_{1}^{i}, …, D_{k}^{i}) + O(λ^{m}),  (29)

with

q_{k+1}^{⋆}(…) = IE q_{k+1}(…) − ∑_{j=0}^{k} IE p_{j+1}(…) IE p_{k+1−j}(…),  (30)

for k = 0, …, m − 1, whenever the service (firing) times σ_{n}^{i} are stochastic. We come back

Figure 5: Approximations of Var(W^{4}).

to our example where σ^{1} = 3 and σ^{2} = 5. Computing the coefficients in (27) according to formula (28) yields

Var (W^{4}) = 9 λ + \frac{189}{4} λ^{2} + \frac{44 731}{60} λ^{3} + \frac{223 291}{18} λ^{4} + \frac{15 335 719}{180} λ^{5} + O(λ^{6}).

Figure 5 gives the plots of various such approximations. The reference curve has been generated by simulation, where the range of the arrival intensity again corresponds to all traffic intensities ρ ∈ [0, 0.75).

4.3 Queues in Tandem — Revisited

We have seen above how the variance of a stationary state variable W^{i} can be expanded to a finite series. In this section we also give an example with stochastic service. Consider the tandem queue example from Section 4.1 with κ = 2. Assume now that σ_{n}^{i} = σ(n) for i = 1, 2, i.e. an arriving customer experiences identical service times at both servers. This model has for instance been closely investigated in [7]. For this network we get

D_{0}^{2} = σ(0),
D_{1}^{2} = max {σ(0) + σ(1), 2 σ(1)},
D_{2}^{2} = max {σ(0) + σ(1) + σ(2), 2 σ(1) + σ(2), σ(1) + 2 σ(2)},
D_{3}^{2} = max {σ(0) + σ(1) + σ(2) + σ(3), 2 σ(1) + σ(2) + σ(3),
σ(1) + 2 σ(2) + σ(3), σ(1) + σ(2) + 2 σ(3)}.

Additionally, we make the assumption that the service time sequence {σ(n)} is a sequence of i.i.d. random variables and σ(n) ∼ 2-Erlang(2) for all n. Using (29) and (30) we get the following approximation of the variance of W^{2}:

Var(W^{2}) = 2 + 28 λ + 43.81 λ^{2} + 104.56 λ^{3} + 242.30 λ^{4} + O(λ^{5}).

Figure 6 shows these approximations of the variance of the stationary variable W^{2} together with a plot generated by simulation. The range of λ given corresponds to all traffic intensities

Figure 6: Approximations of Var(W^{2}).

ρ ∈ [0, \frac{2}{3}). As long as we do not investigate a heavy-traffic situation, say for ρ < \frac{1}{2}, already polynomial approximations of lower order seem to provide satisfactory results.

5 Proof of the Main Theorem

The proof of Theorem 1 is based on representation (6), which allows one to express f(W) as an explicit functional of the Poisson point process T, and proceeds in two steps:

• the application of general techniques to derive expansions of functionals of a Poisson process in function of the intensity parameter λ; this step is considered in Section 5.1;

• the computation of the coefficients of the expansion for the specific functional considered in (6); this step, which is mostly based on combinatorics and real analysis, will be considered in Section 5.2.

5.1 Step 1: Applicability of Two General Expansion Techniques

5.1.1 Factorial Moment Expansions

This first method, which was introduced in [5], was shown to be quite useful for approximating the characteristics of several stochastic models (see e.g. [3], [4] or [6]). Although it can be applied to functionals of general point processes, we will only use it here in the Poisson case.

Expansion Kernels Let ψ be an IR^{α}-valued functional of the marked point process (T, Z) = {(T_{n}, Z_{n})}, i.e. a measurable mapping ψ : M × K^{∞} → IR^{α}, where M is the space of all realizations of the point process T = {T_{n}} and K^{∞} is the space of all sequences Z = {Z_{n}} of potential marks. Note that the sequence T of points may be infinite, finite or empty, whereas the sequence Z of potential marks is always two-sided infinite. Let Z_{n} denote the mark of point T_{n}. We represent a realization {t_{n}} of the point process T by the counting measure μ = ∑_{n} δ_{t_{n}}. By o we denote the null measure, representing an input with no arrivals (i.e. o(IR) = 0). For every x ∈ IR, let the restriction μ|^{x} of μ ∈ M be defined by

μ|^{x}(D) = μ(D ∩ (x, ∞)).

Furthermore, for any x ∈ IR and z ∈ K^{∞}, let

ψ_{x}(μ, z) = ψ(μ|^{x} + δ_{x}, z) − ψ(μ|^{x}, z).  (31)

Let k ≥ 1 be an arbitrary, but fixed integer. For any x_{1}, …, x_{k} ∈ IR, let ψ_{x_{1},...,x_{k}} be defined by iteration of the mapping ψ → ψ_{x}, i.e.

ψ_{x_{1},...,x_{k}}(μ, z) = ( … (ψ_{x_{1}})_{x_{2}} …)_{x_{k}}(μ, z).

Note that the functional ψ_{x_{1},...,x_{k}} can be written in the form

ψ_{x_{1},...,x_{k}}(μ, z) = \{ \begin{matrix}∑_{j=0}^{k}(−1)^{k−j} ∑_{π∈K_{k,j}} ψ(μ|^{x_{k}} + ∑_{i∈π} δ_{x_{i}}, z), & for x_{1} < … < x_{k}, \\ 0 & otherwise, \end{matrix} (32)

where K_{k,j} denotes the collection of all subsets of {1, …, k} containing j elements. Following [5], we call the functional ψ continuous at infinity if

lim_{x→∞} ψ(μ|^{x} + ν, z) = ψ(ν, z),  lim_{x→−∞} ψ(μ|^{x}, z) = ψ(μ, z)  (33)

for all μ, ν ∈ M, z ∈ K^{∞} with ν(IR) < ∞.

Representation Formula For the stationary Poisson process T with intensity λ, and for the stationary sequence Z of K-valued random variables which is independent of T, let P_{λ} denote the distribution of T, and Q the distribution of Z. A slight variant of the following result is given in [8], where the general expansion derived in [5] is specified to the Poisson case. With ψ^{i} we denote the i-th component of the IR^{α}-valued functional ψ. The differential of the Lebesgue measure on IR^{k} is denoted by d(x_{1}, …, x_{k}).

Theorem 2 Let m ≥ 1 be a fixed integer. If the functional ψ is continuous at infinity, if

∫_{IR_{k}} ∫_{K^{∞}} ∫_{M} | ψ_{x_{1},...,x_{k}}^{i}(μ, z) | P_{λ}(dμ) Q(dz) d(x_{1}, …, x_{k}) < ∞,  (34)

for all k = 1, …, m, and if

lim sup_{λ→0} ∫_{IR_{m+1}} ∫_{K^{∞}} ∫_{M} | ψ_{x_{1},...,x_{m+1}}^{i}(μ, z) |

P_{λ}(dμ) Q(dz) d(x_{1}, …, x_{m+1}) < ∞,  (35)

then

IE ψ^{i}(T, Z) = IE ψ^{i}(o, Z) + ∑_{k=1}^{m} λ^{k} ∫_{IR_{k}} IE ψ_{x_{1},...,x_{k}}^{i}(o, Z) d(x_{1}, …, x_{k}) + O(λ^{m+1}). (36)

Applicability of Factorial Moment Expansion Here we consider the following functionals ψ^{i}:

ψ^{i}(T, Z) = f (W^{i}) = f(D_{0}^{i} ⊕ ⨁_{n≥1}(D_{n}^{i} ⊗ T_{−n})),  (37)

where Z = {(A_{n}, B_{n})} and D_{n}^{i} denotes the i-th component of the random vector D_{n}, n ≥ 0, defined in (5). To prove that ψ^{i}(T, Z) satisfies all prerequisites of Theorem 2 and, in particular, integrability conditions (34) and (35), one can proceed in a way similar to that used in [4] for expanding the first moment IE W^{i}. The interested reader should refer to [2], where a detailed checking of these conditions can be found.

This result together with the explicit computation of the coefficients which is done in Section 5.2.1 gives the proof of the first assertion of Theorem 1.

5.1.2 Admissibility

In general, the conditions of Theorem 2 are not sufficient to ensure any differentiability property of IE f(W^{i}) in λ. However, it was shown in [4] that under Cramér type assumptions on the random variables G_{n} considered in (21), the function IE W^{i} is infinitely differentiable in λ in a right neighborhood of 0 and such that (22) holds. The proof of this property is based on the notion of admissibility defined in [10]. This infinite differentiability result can also be extended to the expectations of more general functions of W^{i}. The proof, which can again be found in [2], is similar to that for the first moment case, and we will skip it.

This result together with those of Section 5.2.1 complete the proof of the second assertion of Theorem 1. So, at least in this particular Cramér case, the expansion given in Theorem 1 is a Taylor type expansion.

5.2 Step 2: Calculation of the Coefficients

Assuming now that Theorem 2 can be applied to derive the expansion of IE f(W^{i}), we will calculate the coefficients of λ^{k}. We will first give a recursive representation of these coefficients and will use this result to reveal their explicit form later on.

5.2.1 Recursion Formula

According to Theorem 2 the coefficients of λ^{k} in the expansions are given by

IE ψ^{i}(o, Z),  ∫_{IR_{k}} IE ψ_{x_{1},...,x_{k}}^{i}(o, Z) d(x_{1}, …, x_{k}),  (38)

for k = 0, k ≥ 1, respectively, where ψ^{i} is the functional defined in (37). In what follows, {d_{n}^{i}} will denote a realization of the sequence {D_{n}^{i}} with 1 ≤ i ≤ α′. Introducing the notation q_{1}(d_{0}^{i}) = ψ^{i}(o, z) and

q_{k+1}(d_{0}^{i}, d_{1}^{i}, …, d_{k}^{i}) = ∫_{IR_{k}} ψ_{x_{1},...,x_{k}}^{i}(o, z) d(x_{1}, …, x_{k}),  (39)

for k ≥ 1, we obtain the following result.

Theorem 3 Under monotonicity assumption (8) it holds for each i ∈ {1, …, α′} that q_{1}(d_{0}^{i}) = f(d_{0}^{i}) and

q_{k+1}(d_{0}^{i}, d_{1}^{i}, …, d_{k}^{i}) = ∑_{p=0}^{k−1} ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} \{q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p}, d_{p+1}^{i} − u, …, d_{k}^{i} − u)

− q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p+1}, d_{p+1}^{i} − u, …, d_{k−1}^{i} − u)\} du,  (40)

for all k ≥ 1.

Proof From (37) it follows that

q_{1}(d_{0}^{i}) = ψ^{i}(o, Z) = f(d_{0}^{i}).

Using (32), (37) and (39) for k ≥ 1 we obtain

q_{k+1}(d_{0}^{i}, d_{1}^{i}, …, d_{k}^{i}) = ∫_{0}^{∞} ∫_{s_{1}}^{∞} … ∫_{s_{k−1}}^{∞} ∑_{j=0}^{k}(−1)^{k−j}

∑_{π∈K_{k,j}} f(d_{0}^{i} ⊕ ⨁_{n=1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{2} ds_{1},  (41)

where π_{(n)} denotes the n-th order statistic of π, i.e., the n-th smallest element of the j-tupel π ∈ K_{k,j}. As a next step, we decompose the outer integral in the following way:

∫_{0}^{∞} … = ∑_{p=0}^{k−1} ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} … + ∫_{d_{k}^{i}−d_{0}^{i}}^{∞} … = ∑_{p=0}^{k−1} ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} …  (42)

because from s_{1} ≤ … ≤ s_{k}, d_{0}^{i} ≤ … ≤ d_{k}^{i} and s_{1} ≥ d_{k}^{i} − d_{0}^{i} it follows that s_{m} ≥ s_{1} ≥ d_{k}^{i} − d_{0}^{i} ≥ d_{n}^{i} − d_{0}^{i} and therefore d_{n}^{i} − s_{m} ≤ d_{0}^{i} for all 1 ≤ m, n ≤ k. But this means that the integrand of the integral ranging from d_{k}^{i} − d_{0}^{i} to ∞ is equal to

∑_{j=0}^{k}(−1)^{k−j} ∑_{π∈K_{k,j}} f(d_{0}^{i}) = f(d_{0}^{i}) ∑_{j=0}^{k} \binom{k}{j}(−1)^{k−j} = 0.

For each of the remaining summands in (42), we decompose the summation according to

∑_{j=0}^{k}(−1)^{k−j} ∑_{π∈K_{k,j}} …

= ∑_{j=1}^{k}(−1)^{k−j} ∑_{π∈K_{k,j},π∋1} … + ∑_{j=0}^{k−1}(−1)^{k−j} ∑_{π∈K_{k,j},π∌1} …,  (43)

where the dots are placeholders for the function f(·) given in (41). We first consider that part of the integration involving only the first of these two sums, i.e., for 0 ≤ p ≤ k − 1,

∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} … ∫_{s_{k−1}}^{∞} ∑_{j=1}^{k}(−1)^{k−j} ∑_{π∈K_{k,j},π∋1} f(d_{0}^{i} ⊕ ⨁_{n=1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{1}

= ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} ∫_{s_{1}}^{∞} … ∫_{s_{k−1}}^{∞} ∑_{j=1}^{k}(−1)^{k−j}

∑_{π∈K_{k,j},π∋1} f(d_{0}^{i} ⊕ ⨁_{n=min{j,p}+1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{2} ds_{1},  (44)

because from s_{1} ≥ d_{p}^{i} − d_{0}^{i}, it follows that d_{n}^{i} − s_{π_{(n)}} ≤ d_{n}^{i} − s_{1} ≤ d_{p}^{i} − s_{1} ≤ d_{0}^{i}, for all 1 ≤ n ≤ p, and therefore, in particular, for all 1 ≤ n ≤ min{j, p}. We first consider the case p = 0. Because π ∋ 1, we see that s_{π_{(1)}} = s_{1}, and since in this case d_{1}^{i} − s_{1} ≥ d_{0}^{i}, we get that (44) is equal to

∫_{0}^{d_{1}^{i}−d_{0}^{i}} … ∫_{s_{k−1}}^{∞} ∑_{j=1}^{k}(−1)^{k−j} ∑_{π∈K_{k,j},π∋1} f((d_{1}^{i} − s_{1}) ⊕ ⨁_{n=2}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{1}

= ∫_{0}^{d_{1}^{i}−d_{0}^{i}} ∫_{0}^{∞} ∫_{s_{1}}^{∞} … ∫_{s_{k−2}}^{∞} ∑_{j=0}^{k−1}(−1)^{k−1−j}

∑_{π∈K_{k−1,j}} f((d_{1}^{i} − u) ⊕ ⨁_{n=1}^{j}((d_{n+1}^{i} − u) − s_{π_{(n)}})) ds_{k−1} … ds_{2} ds_{1} du

= ∫_{0}^{d_{1}^{i}−d_{0}^{i}} q_{k}(d_{1}^{i} − u, d_{2}^{i} − u, …, d_{k}^{i} − u) du,  (45)

where we made substitutions s_{1} → u, s_{n} → s_{n}−u for all 2 ≤ n ≤ k, j → j +1 and n → n+1, and finally substituted s_{n} → s_{n−1} for all n = 2, …, k. Comming back to the remaining cases 1 ≤ p ≤ k − 1 we apply the same transformations we just used for p = 0. This yields

∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} … ∫_{s_{k−1}}^{∞} ∑_{j=1}^{k}(−1)^{k−j} ∑_{π∈K_{k,j},π∋1} f(d_{0}^{i} ⊕ ⨁_{n=min{j,p}+1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{1}

= ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} ∫_{0}^{∞} ∫_{s_{1}}^{∞} … ∫_{s_{k−2}}^{∞} ∑_{j=0}^{k−1}(−1)^{k−1−j} ∑_{π∈K_{k−1,j}}

f(d_{0}^{i} ⊕ ⨁_{n=1}^{min{j,p−1}}(d_{0}^{i} − s_{π_{(n)}}) ⊕ ⨁_{n=min{j,p−1}+1}^{j}(d_{n+1}^{i} − u − s_{π_{(n)}}))

ds_{k−1} … ds_{2} ds_{1} du
= ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p}, d_{p+1}^{i} − u, …, d_{k}^{i} − u) du,  (46)

for p = 1, …, k − 1. Aggregating (45) and (46), we obtain for (44)

∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} … ∫_{s_{k−1}}^{∞} ∑_{j=1}^{k}(−1)^{k−j} ∑_{π∈K_{k,j},π∋1} f(d_{0}^{i} ⊕ ⨁_{n=1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{1}

= ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p}, d_{p+1}^{i} − u, …, d_{k}^{i} − u) du,  (47)

for all p = 0, …, k − 1. Using analogous reasoning for the remaining second part of the integration, we get

∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} … ∫_{s_{k−1}}^{∞} ∑_{j=0}^{k−1}(−1)^{k−j} ∑_{π∈K_{k,j},π∌1} f(d_{0}^{i} ⊕ ⨁_{n=1}^{j}(d_{n}^{i} − s_{π_{(n)}})) ds_{k} … ds_{1}

= − ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p+1}, d_{p+1}^{i} − u, …, d_{k−1}^{i} − u) du,  (48)

for all 0 ≤ p ≤ k − 1, and thus

q_{k+1}(d_{0}^{i}, d_{1}^{i}, …, d_{k}^{i}) = ∑_{p=0}^{k−1} ∫_{d_{p}^{i}−d_{0}^{i}}^{d_{p+1}^{i}−d_{0}^{i}} \{q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p}, d_{p+1}^{i} − u, …, d_{k}^{i} − u)

− q_{k}(\underbrace{d_{0}^{i}, …, d_{0}^{i}}_{p+1}, d_{p+1}^{i} − u, …, d_{k−1}^{i} − u)\} du,

which is the statement of Theorem 3.

5.2.2 Explicit Solution

In this section we want to derive the explicit representation (19) of the functions q_{k}(…), k ≥ 1, that show up in the calculation of the coefficients IE q_{k+1}(D_{0}^{i}, …, D_{k}^{i}) of λ^{k} in the series expansion of IE f(W^{i}) considered in Theorem 1.

Theorem 4 For k ≥ 1 and 0 ≤ x_{0} ≤ x_{1} ≤ … ≤ x_{k}, a solution to the integral recursion equation (40) stated in Theorem 3 is given by the functions q_{k} defined in Theorem 1, particularly in (19).

Proof We show (19) by induction with respect to k ≥ 1. Using (40) for k = 1 yields

q_{2}(x_{0}, x_{1}) = ∫_{0}^{x_{1}−x_{0}} {q_{1}(x_{1} − u) − q_{1}(x_{0})} du
= [ − F (x_{1} − u) − f(x_{0})u ]_{0}^{x_{1}−x_{0}}
= (F (x_{1}) − F (x_{0})) − f(x_{0}) {p_{1}(x_{1}) − p_{1}(x_{0})},

where we wrote F (x) for F^{[1]}(x) and thus, (19) holds for k = 1. In the following k will always be assumed to be greater than one, i.e. k ≥ 2. Equation (40) can be rewritten in the form:

q_{k+1}(x_{0}, x_{1}, …, x_{k})
= ∑_{p=0}^{k−1} ∫_{x_{p}−x_{0}}^{x_{p+1}−x_{0}} q_{k}(\underbrace{x_{0}, …, x_{0}}_{p}, x_{p+1} − u, …, x_{k} − u) du  (49)

− ∑_{p=0}^{k−1} ∫_{x_{p}−x_{0}}^{x_{p+1}−x_{0}} q_{k}(\underbrace{x_{0}, …, x_{0}}_{p+1}, x_{p+1} − u, …, x_{k−1} − u) du.  (50)

In what follows, (49) and (50) will be evaluated separately. Now assume (19) to be true for some natural k ≥ 2, i.e.

q_{k}(x_{1}, x_{2}, …, x_{k})
= ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{n+1}) − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{n+1})
{p_{k−1−j}(x_{n+2}, …, x_{k−j+n}) − p_{k−1−j}(x_{n+1}, …, x_{k−j+n−1})}.  (51)

Since the polynomials p_{n}(x_{0}, …, x_{n−1}) are translation-invariant with respect to a constant translation of all arguments, for n ≥ 2, and since p_{1}(x) = x implies p_{1}(x_{1} − u) − p_{1}(x_{0} − u) = p_{1}(x_{1}) − p_{1}(x_{0}), it follows from (51) that for p = 0

∫_{0}^{x_{1}−x_{0}} q_{k}(x_{1} − u, …, x_{k} − u) du

= ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}(F^{[k]}(x_{n+1}) − F^{[k]}(x_{n+1} + x_{0} − x_{1}))  (52)

− ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}\{p_{k−1−j}(x_{n+2}, …, x_{k−j+n}) − p_{k−1−j}(x_{n+1}, …, x_{k−j+n−1})\}

(−1)^{j−n}(F^{[j+1]}(x_{n+1}) − F^{[j+1]}(x_{n+1} + x_{0} − x_{1})).  (53)

For 1 ≤ p ≤ k − 1 the integrand is given by:

q_{k}(\underbrace{x_{0}, …, x_{0}}_{p}, x_{p+1} − u, …, x_{k} − u)

= ∑_{n=0}^{p−1} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{0}) + ∑_{n=p}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{n+1} − u)

− ∑_{n=0}^{p−1} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0})\{p_{k−1−j}(\underbrace{x_{0}, …, x_{0}}_{p−n−1}, x_{p+1} − u, …, x_{k−j+n} − u)

− p_{k−1−j}(\underbrace{x_{0}, …, x_{0}}_{p−n}, x_{p+1} − u, …, x_{k−j+n−1} − u)\}

− ∑_{n=p}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{n+1} − u)\{p_{k−1−j}(x_{n+2} − u, …, x_{k−j+n} − u)

− p_{k−1−j}(x_{n+1} − u, …, x_{k−j+n−1} − u)},

and therefore, if we use the translation-invariance property (17) of the polynomials p_{n} to shift all arguments by u and then make the substitution x_{0} + u → v, we see that

∫_{x_{p}−x_{0}}^{x_{p+1}−x_{0}} q_{k}(\underbrace{x_{0}, …, x_{0}}_{p}, x_{p+1} − u, …, x_{k} − u) du

= ∑_{n=0}^{p−1} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{0})\{p_{1}(x_{p+1}) − p_{1}(x_{p})\}  (54)

+ ∑_{n=p}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}\{F^{[k]}(x_{n+1} + x_{0} − x_{p})

− F^{[k]}(x_{n+1} + x_{0} − x_{p+1})}  (55)

− ∑_{n=0}^{p−1} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0}) ∫_{x_{p}}^{x_{p+1}} \{p_{k−1−j}(\underbrace{v, …, v}_{p−n−1}, x_{p+1}, …, x_{k−j+n})

− p_{k−1−j}(\underbrace{v, …, v}_{p−n}, x_{p+1}, …, x_{k−j+n−1})\} dv  (56)

− ∑_{n=p}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}\{F^{[j+1]}(x_{n+1} + x_{0} − x_{p}) − F^{[j+1]}(x_{n+1} + x_{0} − x_{p+1})\}

{p_{k−1−j}(x_{n+2}, …, x_{k−j+n}) − p_{k−1−j}(x_{n+1}, …, x_{k−j+n−1})}  (57)
= I_{p} + II_{p} + III_{p} + IV_{p},

for 1 ≤ p ≤ k − 1, where I_{p}, II_{p}, III_{p} and IV_{p} are compact notations for the expressions given by (54), (55), (56) and (57), respectively. So we see that (49) can be represented as the sum of the terms given by (52) and (53) as well as the expressions just introduced for 1 ≤ p ≤ k −1 above. We will now continue by summing up these terms in the following way:

∑_{p=1}^{k−1} I_{p} = F^{[k−1]}(x_{0}) ∑_{n=0}^{k−2} \binom{k − 1}{n}(−1)^{k−1−n} ∑_{p=n}^{k−2} \{p_{1}(x_{p+2}) − p_{1}(x_{p+1})\}

= −F^{[k−1]}(x_{0}) ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}p_{1}(x_{n+1}),  (58)

where we substituted p → p + 1 and then interchanged the order of summation to get the first equality. The last line was obtained after calculating the telescope-sum, splitting up the result into two summations and using the fact that

∑_{n=0}^{k−2} \binom{k − 1}{n}(−1)^{k−1−n} = −1.

If we compare (52) to (55) we notice that (52) can just as well be interpreted as a II_{0}. Thus, a substitution n → n − 1 and similar calculations as before yield

∑_{p=0}^{k−1} II_{p} = ∑_{n=1}^{k} \binom{k − 1}{n − 1}(−1)^{k−n} ∑_{p=0}^{n−1} \{F^{[k]}(x_{n} + x_{0} − x_{p}) − F^{[k]}(x_{n} + x_{0} − x_{p+1})\}

= ∑_{n=1}^{k} \binom{k − 1}{n − 1}(−1)^{k−n}F^{[k]}(x_{n}),  (59)

as well as

∑_{p=1}^{k−1} III_{p} = − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0}) ∑_{p=n}^{k−2} ∫_{x_{p+1}}^{x_{p+2}}

\{p_{k−1−j}(\underbrace{v, …, v}_{p−n}, x_{p+2}, …, x_{k−j+n})

− p_{k−1−j}(\underbrace{v, …, v}_{p−n+1}, x_{p+2}, …, x_{k−j+n−1})\} dv

= − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0})p_{k−j}(x_{n+1}, …, x_{k−j+n}).  (60)

The second equality holds because the polynomials p_{k−1−j} cancel out whenever p − n ≥ k − 1 − j or, equivalently, if p ≥ k − 1 − j + n, and therefore

∑_{p=n}^{k−2} ∫_{x_{p+1}}^{x_{p+2}} \{p_{k−1−j}(\underbrace{v, …, v}_{p−n}, x_{p+2}, …, x_{k−j+n})

− p_{k−1−j}(\underbrace{v, …, v}_{p−n+1}, x_{p+2}, …, x_{k−j+n−1})\} dv

= ∑_{p=0}^{k−2−j} ∫_{x_{p+n+1}}^{x_{p+n+2}} \{p_{k−1−j}(\underbrace{v, …, v}_{p}, x_{p+n+2}, …, x_{k−j+n})

− p_{k−1−j}(\underbrace{v, …, v}_{p+1}, x_{p+n+2}, …, x_{k−j+n−1})\} dv

= p_{k−j}(x_{n+1}, …, x_{k−j+n}),

where the last equality follows from (18). It can be seen that the IV_{p}’s show up only for 1 ≤ p ≤ k −2. Similarly as we did for the II_{p}’s above, we realize that (53) can be interpreted as a IV_{0}. So we include this term in the summation of the IV_{p}’s and obtain

∑_{p=0}^{k−1} IV_{p} = − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n} ∑_{p=0}^{n} \{F^{[j+1]}(x_{n+1} + x_{0} − x_{p})

− F^{[j+1]}(x_{n+1} + x_{0} − x_{p+1})}
{p_{k−1−j}(x_{n+2}, …, x_{k−j+n}) − p_{k−1−j}(x_{n+1}, …, x_{k−j+n−1})}
= − ∑_{n=1}^{k−1} ∑_{j=n}^{k−1} \binom{j − 1}{n − 1}(−1)^{j−n}\{F^{[j]}(x_{n}) − F^{[j]}(x_{0})\}
{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})}.  (61)

We can approach expression (50) in an almost identical way. For p = k − 1 we obtain q_{k}(x_{0}, …, x_{0}) = 0 and for p = 0, …, k − 2

∫_{x_{p}−x_{0}}^{x_{p+1}−x_{0}} q_{k}(\underbrace{x_{0}, …, x_{0}}_{p+1}, x_{p+1} − u, …, x_{k−1} − u) du

= ∑_{n=0}^{p} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{0})\{p_{1}(x_{p+1}) − p_{1}(x_{p})\}  (62)

+ ∑_{n=p+1}^{k−1}\binom{k − 1}{n}(−1)^{k−1−n}\{F^{[k]}(x_{n} + x_{0} − x_{p}) − F^{[k]}(x_{n} + x_{0} − x_{p+1})\}  (63)

− ∑_{n=0}^{p} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0}) ∫_{x_{p}}^{x_{p+1}} \{p_{k−1−j}(\underbrace{v, …, v}_{p−n}, x_{p+1}, …, x_{k−j+n−1})

− p_{k−1−j}(\underbrace{v, …, v}_{p−n+1}, x_{p+1}, …, x_{k−j+n−2})\} du  (64)

− ∑_{n=p+1}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}\{F^{[j+1]}(x_{n} + x_{0} − x_{p}) − F^{[j+1]}(x_{n} + x_{0} − x_{p+1})\}

{p_{k−1−j}(x_{n+1}, …, x_{k−j+n−1}) − p_{k−1−j}(x_{n}, …, x_{k−j+n−2})}  (65)
= i_{p} + ii_{p} + iii_{p} + iv_{p},

with i_{p}, ii_{p}, iii_{p} and iv_{p} representing expression (62), (63), (64) and (65), respectively. Using similar calculations as before we get the following expressions for the summations of the different terms:

∑_{p=0}^{k−2} i_{p} = −F^{[k−1]}(x_{0}) ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}p_{1}(x_{n}),  (66)

∑_{p=0}^{k−2} ii_{p} = − ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−n}F^{[k]}(x_{n}),  (67)

∑_{p=0}^{k−2} iii_{p} = − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0}) p_{k−j}(x_{n}, …, x_{k−j+n−1}),  (68)

∑_{p=0}^{k−2} iv_{p} = ∑_{n=1}^{k−2} ∑_{j=n+1}^{k−1} \binom{j − 1}{n}(−1)^{j−n}\{F^{[j]}(x_{n}) − F^{[j]}(x_{0})\}
{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})}  . (69)

So we know that q_{k+1}(x_{0}, x_{1}, …, x_{k}) for k ≥ 1 is given by the sum of (58), (59), (60) and (61), subtracted by the sum of (66), (67), (68) and (69). To continue we evaluate the following differences:

∑_{p=1}^{k−1} I_{p} − ∑_{p=0}^{k−2} i_{p} = − ∑_{n=0}^{k−1} \binom{k − 1}{n}(−1)^{k−1−n}F^{[k−1]}(x_{0})\{p_{1}(x_{n+1}) − p_{1}(x_{n})\},  (70)

whereas

∑_{p=0}^{k−1} II_{p} − ∑_{p=0}^{k−2} ii_{p} = F^{[k]}(x_{k}) + ∑_{n=1}^{k−1} \{\binom{k − 1}{n − 1} + \binom{k − 1}{n}\} (−1)^{k−n}F^{[k]}(x_{n})
+ (−1)^{k}F^{[k]}(x_{0})
= ∑_{n=0}^{k} \binom{k}{n}(−1)^{k−n}F^{[k]}(x_{n}),  (71)

since \binom{k−1}{n−1} + \binom{k−1}{n} = \binom{k}{n}. Furthermore,

∑_{p=1}^{k−1} III_{p} − ∑_{p=0}^{k−2} iii_{p} = − ∑_{n=0}^{k−2} ∑_{j=n}^{k−2} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{0}) \{ p_{k−j}(x_{n+1}, …, x_{k−j+n})
− p_{k−j}(x_{n}, …, x_{k−j+n−1}) }  (72)

and

∑_{p=0}^{k−2} IV_{p} − ∑_{p=0}^{k−2} iv_{p} = − ∑_{n=0}^{k−1} ∑_{j=n}^{k−1} \binom{j}{n}(−1)^{j−n}\{F^{[j]}(x_{n}) − F^{[j]}(x_{0})\}
{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})},  (73)

where we extended the summation to start from n = 0 instead of n = 1. This does not change the overall value of (73) as all these additional terms cancel out. Looking at the sum of (70) and (72), we realize that the result cancela out with the part of (73) involving F^{[j]}(x_{0}). Thus, we are left with (71) and the remaining part of (73), and we eventually obtain that

q_{k+1}(x_{0}, x_{1}, …, x_{k}) = ∑_{n=0}^{k} \binom{k}{n}(−1)^{k−n}F^{[k]}(x_{n}) − ∑_{n=0}^{k−1} ∑_{j=n}^{k−1} \binom{j}{n}(−1)^{j−n}F^{[j]}(x_{n})
{p_{k−j}(x_{n+1}, …, x_{k−j+n}) − p_{k−j}(x_{n}, …, x_{k−j+n−1})},

which is the statement of Theorem 4.

References

[1] Baccelli, F., G. Cohen, G.J. Olsder and J.P. Quadrat (1992) Synchronization and Linearity. J. Wiley & Sons, Chichester.

[2] Baccelli, F., S. Hasenfuss and V. Schmidt (1996) Expansions for steady-state characteristics in (max,+)–linear systems. Preprint. Ulm University / INRIA Report N^{∘} 2785.

[3] Baccelli, F., S. Hasenfuss and V. Schmidt (1997) Transient and stationary waiting times in (max,+)–linear systems with Poisson input. Queueing Systems. (to appear)

[4] Baccelli, F. and V. Schmidt (1996) Taylor series expansions for Poisson driven (max,+)– linear systems. Ann. Appl. Probab. 6, 138–185.

[5] Błaszczyszyn, B. (1995) Factorial moment expansion for stochastic systems. Stoch. Proc. Appl. 56, 321–335.

[6] Błaszczyszyn, B., Rolski, T. and V. Schmidt (1995) Light-traffic approximations in queues and related stochastic models. In: J.H. Dshalalow (ed.) Advances in Queueing: Theory, Methods and Open Problems. CRC Press, Boca Raton, 379–406.

[7] Boxma, O.J. (1979) On a tandem queueing model with identical service times at both counters. Adv. Appl. Probab. 11, 616–659.

[8] Kroese, D.P. and V. Schmidt (1996) Light-traffic analysis for queues with spatially distributed arrivals. Math. Oper. Res. 21, 135–157.

[9] Murata T. (1989) Petri nets: properties, analysis and applications. Proc. of the IEEE 77(4), 541–580.

[10] Reiman, M.I. and B. Simon (1989) Open queueing systems in light traffic. Math. Oper. Res. 14, 26–59.