# A Fluid Model for a Relay Node in an Ad Hoc Network: Evaluation of Resource Sharing Policies

Authors: Michel Mandjes; Werner Scheinhardt
Journal: Journal of Applied Mathematics and Stochastic Analysis
 (2008)

---

## Abstract

Fluid queues offer a natural framework for analyzing waiting times in a relay node of an ad hoc network. Because of the resource sharing policy applied, the input and output of these queues are coupled. More specifically, when there aren users who wish to transmit data through a specific node, each of them obtains a share 1/(n+W) of the service capacity to feed traffic into the queue of the node, whereas the remaining fraction W/(n+W) is used to serve the queue; here W>0 is a free design parameter. Assume now that jobs arrive at the relay node according to a Poisson process, and that they bring along exponentially distributed amounts of data. The case W=1 has been addressed before; the present paper focuses on the intrinsically harder case W>1, that is, policies that give more weight to serving the queue. Four performance metrics are considered: (i) the stationary workload of the queue, (ii) the queueing delay, that is, the delay of a “packet”  (a fluid particle) that arrives at an arbitrary point in time, (iii) the flow transfer delay, (iv) the sojourn time, that is, the flow transfer time increased by the time it takes before the last fluid particle of the flow is served. We explicitly compute the Laplace transforms of these random variables.


---

## Body

## 1. Introduction

Ad hoc networks
are self-configuring networks of
mobile routers, connected by wireless links. They enable infrastructure-free
communication: no fixed equipment is needed, but instead each client acts as a
hub. When information needs to be transmitted across the network, it is sent
from the sender to the receiver by relaying the packets along the intermediate
hubs. An excellent survey on ad hoc networks, with special emphasis on
quality-of-service aspects, is [1].On a somewhat more abstract level the nodes in ad hoc
networks can be regarded asqueues:
information packets arrive and are relayed, and when the arrival rate
(temporarily) exceeds the departure rate, the buffer content of the queue
builds up. These queues, however, have the interesting modelling feature that
the available transmission capacity at any specific node is used both to (i)
“pull” information packets from the “predecessor hubs” into the queue and (ii)
“push” information packets from the queue towards “successor hubs” (and
eventually the destination client).Now, consider the situation that at some point in time
there aren stations that transmit traffic via the same
relay node. If the standard resource sharing policy is used, then the
relay-node is assigned a share 1/(n+1) of the available medium capacity, which is the
same fraction as is allocated to each of the n “sending nodes.” In other words, as soon as n>1,
the node's input rate exceeds its output rate, and hence the excess traffic
accumulates in the node's buffer; only when n=0 the queue drains. This explains why relay
nodes are prone to becoming bottlenecks.The above observations have motivated the development
of alternative resource sharing policies, that assign more “weight” to serving
the relay node; see for instance [2] and references therein. With the so-called enhanced distributed channel access
(edca) protocol, it will be possible to set a parameter W>0,
such that each of the n sending nodes obtains a fraction 1/(n+W) of the capacity, while the remaining w/(n+W) is allocated to serving the queue. Clearly,
when W>1 this has a benign effect on the buffer content
of the queue, compared to the standard resource sharing policy described above:
the queue drains for all n<W (rather than just for n=0). The price paid is that the traffic remains
longer at the sending nodes.There is one important modelling feature, that applies
to the caseW>1,
that needs to be mentioned here. When the queue is empty and the number of
sending nodes n is below W,
it does not make sense to assign each of the sending notes just a share 1/(n+W):
it would imply that the (available) service rate is strictly larger than the
input rate, and that a fraction (W−n)/(n+W)>0 is left unused. For that reason, the
edca protocol (ieee 802.11E) was
augmented with an “idle mode”: if the queue is idle and n<W,
then half of the capacity is allocated to serving the queue of the relay node,
whereas the other half is shared equally among the n sending nodes (such that the input and output
rates are equal, the queue remains empty, and all available capacity is used).
Notice that when W≤1 this special rate allocation (during periods
in which the buffer is empty) is not required, as it cannot be that both the
buffer is empty and that there are n<W jobs in the system.The caseW=1,
corresponding to the standard resource sharing policy, was proposed and
analysed in detail in [3, 4]. Seen from a queueing-theoretic perspective, this is
a model “with coupled input and output,” in that the capacity is shared between
the input and output process. It was assumed that flows (or jobs) arrive at the
relay node according to a Poisson process and initiate a data transfer. For the
special case of exponentially distributed flow sizes, [4] explicitly gave the Laplace
transforms (and tail asymptotics) of several performance measures.Importantly, the caseW=1 (and in fact also W<1) has nice features that are lost for W>1.
As a consequence, the analysis of [4] for W=1 does not carry over to W>1. The main differences are the
following.(1) In the first place, as we explained, the caseW=1 does not require the “idle mode” rate
allocation that was introduced above (during periods in which the buffer is
empty). As a consequence, the number of flows present evolves independently of
the buffer content process (and hence the distribution of the number of flows
present can be computed independently of the distribution of the buffer content
of the queue). This nice property is lost in the case W>1;
one could say that there is then some sort of “feedback” from the buffer content
to the flows, in the sense that the buffer content has impact on the flow
behaviour, and hence the distributions of the number of sources present and the
buffer content cannot be determined separately.(ii) In the second place, suppose we wish to analyse
the flow transfer delay, defined as the time between the arrival of the flow
and the epoch that its last fluid particle has been transmitted into the queue.
We know that forW=1 the queue cannot become empty during this flow
transfer. Therefore, all the “state information” that we have to keep track of
is the number of flows present when entering (and not the buffer content); for W>1 the queue can become empty, and therefore we have to
take into account the buffer content, as seen upon arrival, as well.(iii) In caseW=1,
the buffer content decreases only during periods that there are no flows in the
system, and these intervals are exponentially distributed; it turned out in
[4] that this
property entails that a direct translation is possible in terms of a related
classical M/G/1 queueing model. For W>1,
the periods of net output are not exponentially distributed, and therefore we
lose this nice feature.The above
explains that the analysis forW>1 is considerably more involved than for W=1.For various values ofW,
the model described above has been extensively validated in
[2], by ad hoc network
simulations that include all the details of the widely used
ieee 802.11 mac-protocol. As
indicated above, the alternative resource sharing policies can be enforced in
real systems by deploying the recently standardised ieee
802.11E protocol; [2]
also indicates how to map the parameter settings of ieee
802.11E on our model parameters.The goal of the present paper is to extend the results
of [4] to the case W>1. As in [4], four performance metrics are considered: (i) the
stationary workload of the queue, (ii) the queueing delay, that is, the delay of
a “packet” (a fluid particle) that arrives at the queue at an arbitrary point
in time, (iii) the flow transfer delay, that is, the time elapsed between
arrival of a flow and the epoch that all its traffic has been put into the
queue, and (iv) the sojourn time, that is, the flow transfer time increased by
the time it takes before the last fluid particle of the flow is served.We introduce the model (including a graphical
illustration) and notation in Section2—it is noted that in our model an admission control policy is in
force. We also present some preliminaries as well as the stability condition.
In Section 3  the steady-state
workload (i.e., buffer content) distribution of the queue is characterised.
This is done by relying on techniques from [5–7]. Unfortunately, we cannot exploit the resemblance a
related M/G/1 queueing model, as was possible in [4]. We find the distribution
function of the steady-state workload, in terms of the solution of an
eigensystem and a set of linear equations. In Section 4, we study the queueing delay of
a tagged fluid particle that arrived at time 0. A full characterisation of its
Laplace transform can be given; the computations turn out to be relatively
straightforward, based on the observation that during the queueing delay the
queue cannot enter the “idle mode.”Above we already indicated that the analysis of the
flow transfer delay is much more involved than forW=1,
mainly due to the fact that the buffer can become empty during the flow
transfer, so that the allocation gets into the “idle mode.” The derivation of
the Laplace transform requires the solution of various complex systems of
equations. The results can be found in Section 5. In order to prove that the
number of equations matches the number of unknowns, we need to show that a
certain eigensystem has sufficiently many eigenvalues in the right half-plane;
this we showed by using an elegant and powerful lemma of Sonneveld [8].Section6 concentrates on the so-called sojourn time, which is
defined as the flow transfer time increased by the time it takes before the
last fluid particle of the flow has left the queue; in other words, the sojourn
time is the time it takes for the flow to go through the relay node. Relying on
the results of Section 5, the sojourn time can be decomposed into a number of
known components. Section 7 concludes and identifies a number of topics for future research.


## 2. Model and Preliminaries

In this
section, we first give a detailed description of our model and introduce
notation. Then, we derive the stability condition.


### 2.1. Model

The following
model was verbally motivated in the Introduction. Consider a queueing system at
which flows arrive according to a Poisson process transmits traffic into a
queue (which is served in afifo manner) and leaves when
ready. When there are n flows active and the queue is nonempty, each
flow transmits traffic into the queue at rate C/(n+W),
while a rate WC/(n+W) is used to serve the queue; as a consequence,
the queue drains when the number of flows present is below W.
For ease, we assume that W is noninteger; we come back to this issue
later in Section 7. When there are n flows active and the queue is empty, all flows
transmit at rate C/(2n),
while the queue is served at rate C/2,
so that the queue remains empty.Suppose that we impose the admission control policy
that the system accommodates maximallyN∈ℕ flows simultaneously; in this way, each active
flow is guaranteed at least a transmission rate C/(N+W),
and the queue at least WC/(N+W). We assume that N>W (as otherwise the queue remains empty).The above dynamics define a queueing process, for any
given initial buffer level and initial number of flows present; we denote the
buffer content at timet by Wt.
We let Nt denote the number of flows present (i.e.,
feeding traffic into the queue) at time t.
A pictorial illustration is given in Figure 1.Illustration of the model. The top panel depicts the
number of flowsNt as a function of time, the middle panel the
rate allocated to each source as a function of time, and the bottom panel the
workload Wt as a function of time; we have chosen W=3/2.
Note that when Nt=1 the rate per source equals C/2 when the buffer is empty, and 2C/5 when the buffer is nonempty. Also note that
the queue builds up when Nt≥2.
(a)(b)(c)We introduce the following notation.(i)Wt>0:
the “busy mode.” It is not hard to see, under the assumption of exponentially
distributed flow sizes (with mean μ−1) and interarrival times (with mean λ−1), that, during periods thatWt>0,
the process Nt behaves as a Markov chain on {0,…,N},
with
generator matrix
(2.1)Qb:=(−λλμ1c−μ1c−λλμ2c−μ2c−λλ⋱⋱⋱μnc−μnc),where
μn:=μn/(n+W);
the subscript “b” stands for “busy.” We define νn:=μnC+λ1{n<N}.WhenNt=n,
the aggregate traffic rate generated by the flows is rI,n:=Cn/(n+W),
while the queue's output rate is rO,n:=WC/(n+W),
such that the net rate of change of the queue is(2.2)ra,n:=ri,n−ro,n=cn−wn+w.Define RI:=diag{rI}, RO:=diag{rO},
and RA:=RI−RO.Busy periods are periods in which Wt is positive all the time. With n+:=⌈W⌉, it is evident that the number of active flows at the beginning of a
busy period equals n+.
The number of active flows at the end of a busy period is in {0,…,n−},
with n−:=⌊W⌋.(ii)Wt=0:
the “idle mode.” Let idle periods be periods in which Wt=0 all the time. An idle period ends as soon as Nt=n+.
During the idle period necessarily Nt∈{0,…,n−}. One could say that Nt behaves as a Markov chain on {0,…,n−} until Nt jumps from n− to n+ (i.e., the start of a new busy period). The
corresponding rate matrix (which is not a bona fide generator matrix)
is(2.3)Qi:=(−λλμc/2−μc/2−λλμc/2−μc/2−λλ⋱⋱⋱μc/2−μc/2−λ);the subscript “i” stands for
“idle.”


### 2.2. Stability Condition

To make sure
that the steady-state workloadW⋆ is finite a.s., the mean drift should be negative
when Wt is large. Since Nt behaves essentially like a stationary Markov
process with generator Qb when Wt is large, it follows that Wt can only escape to ∞ when ∑n=0NπnrA,n≥0,
denoting by π→≡(π0,…,πN)T the invariant measure of Qb.
Hence, we should require that ∑n=0NπnrA,n<0. Elementary Markov chain analysis yields that,
with ϱ:=λ/(μC),(2.4)πn=ρnB(n+w,n)∑m=0nϱmB(m+w,m),where, for k≥ℓ≥0,𝔹(⋅,⋅) can be regarded as a “generalised binomial
coefficient”:(2.5)B(k,ℓ):=Γ(k+1)Γ(ℓ+1)Γ(k−ℓ+1).The stability condition becomes ∑n=0NrA,nπn<0.It is instructive to show how this condition
simplifies in the situation ofN→∞. Due to (recognize the probability density
function of the negative binomial distribution)(2.6)∑m=0∞ϱmB(m+w,m)=1(1−ϱ)w+1,we have to verify
whether(2.7)∑n=0∞πn⋅n−wn+w=∑n=0∞ρnB(n+w,n)(1−ϱ)w+1⋅n−wn+w<0.Using identity (2.6) again,
writing(2.8)n−wn+w=2nn+w−1,and observing that(2.9)∑n=0∞ϱnB(n+w,n)2nn+w=∑n=1∞ϱnB(n+w,n)2nn+w=2∑n=1∞ϱnB(n+w−1,n−1)=2ϱ∑n=0∞ϱnB(n+w,n)=2ϱ(1−ϱ)w+1,it is readily verified that the
stability requirement reduces to 2ϱ−1<0.
In other words, for the system to be stable, it is required that ϱ<1/2 (irrespective of the value of W). This result makes sense as essentially all
traffic has to be “served” twice: first it has to be transmitted from the
sources to the queue, and then it has to be served by the queue.


## 2.1. Model

The following
model was verbally motivated in the Introduction. Consider a queueing system at
which flows arrive according to a Poisson process transmits traffic into a
queue (which is served in afifo manner) and leaves when
ready. When there are n flows active and the queue is nonempty, each
flow transmits traffic into the queue at rate C/(n+W),
while a rate WC/(n+W) is used to serve the queue; as a consequence,
the queue drains when the number of flows present is below W.
For ease, we assume that W is noninteger; we come back to this issue
later in Section 7. When there are n flows active and the queue is empty, all flows
transmit at rate C/(2n),
while the queue is served at rate C/2,
so that the queue remains empty.Suppose that we impose the admission control policy
that the system accommodates maximallyN∈ℕ flows simultaneously; in this way, each active
flow is guaranteed at least a transmission rate C/(N+W),
and the queue at least WC/(N+W). We assume that N>W (as otherwise the queue remains empty).The above dynamics define a queueing process, for any
given initial buffer level and initial number of flows present; we denote the
buffer content at timet by Wt.
We let Nt denote the number of flows present (i.e.,
feeding traffic into the queue) at time t.
A pictorial illustration is given in Figure 1.Illustration of the model. The top panel depicts the
number of flowsNt as a function of time, the middle panel the
rate allocated to each source as a function of time, and the bottom panel the
workload Wt as a function of time; we have chosen W=3/2.
Note that when Nt=1 the rate per source equals C/2 when the buffer is empty, and 2C/5 when the buffer is nonempty. Also note that
the queue builds up when Nt≥2.
(a)(b)(c)We introduce the following notation.(i)Wt>0:
the “busy mode.” It is not hard to see, under the assumption of exponentially
distributed flow sizes (with mean μ−1) and interarrival times (with mean λ−1), that, during periods thatWt>0,
the process Nt behaves as a Markov chain on {0,…,N},
with
generator matrix
(2.1)Qb:=(−λλμ1c−μ1c−λλμ2c−μ2c−λλ⋱⋱⋱μnc−μnc),where
μn:=μn/(n+W);
the subscript “b” stands for “busy.” We define νn:=μnC+λ1{n<N}.WhenNt=n,
the aggregate traffic rate generated by the flows is rI,n:=Cn/(n+W),
while the queue's output rate is rO,n:=WC/(n+W),
such that the net rate of change of the queue is(2.2)ra,n:=ri,n−ro,n=cn−wn+w.Define RI:=diag{rI}, RO:=diag{rO},
and RA:=RI−RO.Busy periods are periods in which Wt is positive all the time. With n+:=⌈W⌉, it is evident that the number of active flows at the beginning of a
busy period equals n+.
The number of active flows at the end of a busy period is in {0,…,n−},
with n−:=⌊W⌋.(ii)Wt=0:
the “idle mode.” Let idle periods be periods in which Wt=0 all the time. An idle period ends as soon as Nt=n+.
During the idle period necessarily Nt∈{0,…,n−}. One could say that Nt behaves as a Markov chain on {0,…,n−} until Nt jumps from n− to n+ (i.e., the start of a new busy period). The
corresponding rate matrix (which is not a bona fide generator matrix)
is(2.3)Qi:=(−λλμc/2−μc/2−λλμc/2−μc/2−λλ⋱⋱⋱μc/2−μc/2−λ);the subscript “i” stands for
“idle.”


## 2.2. Stability Condition

To make sure
that the steady-state workloadW⋆ is finite a.s., the mean drift should be negative
when Wt is large. Since Nt behaves essentially like a stationary Markov
process with generator Qb when Wt is large, it follows that Wt can only escape to ∞ when ∑n=0NπnrA,n≥0,
denoting by π→≡(π0,…,πN)T the invariant measure of Qb.
Hence, we should require that ∑n=0NπnrA,n<0. Elementary Markov chain analysis yields that,
with ϱ:=λ/(μC),(2.4)πn=ρnB(n+w,n)∑m=0nϱmB(m+w,m),where, for k≥ℓ≥0,𝔹(⋅,⋅) can be regarded as a “generalised binomial
coefficient”:(2.5)B(k,ℓ):=Γ(k+1)Γ(ℓ+1)Γ(k−ℓ+1).The stability condition becomes ∑n=0NrA,nπn<0.It is instructive to show how this condition
simplifies in the situation ofN→∞. Due to (recognize the probability density
function of the negative binomial distribution)(2.6)∑m=0∞ϱmB(m+w,m)=1(1−ϱ)w+1,we have to verify
whether(2.7)∑n=0∞πn⋅n−wn+w=∑n=0∞ρnB(n+w,n)(1−ϱ)w+1⋅n−wn+w<0.Using identity (2.6) again,
writing(2.8)n−wn+w=2nn+w−1,and observing that(2.9)∑n=0∞ϱnB(n+w,n)2nn+w=∑n=1∞ϱnB(n+w,n)2nn+w=2∑n=1∞ϱnB(n+w−1,n−1)=2ϱ∑n=0∞ϱnB(n+w,n)=2ϱ(1−ϱ)w+1,it is readily verified that the
stability requirement reduces to 2ϱ−1<0.
In other words, for the system to be stable, it is required that ϱ<1/2 (irrespective of the value of W). This result makes sense as essentially all
traffic has to be “served” twice: first it has to be transmitted from the
sources to the queue, and then it has to be served by the queue.


## 3. Buffer Content Distribution

In this
section, we study the steady-state workloadW⋆ of the queue introduced in the previous
section (jointly with the steady-state number N⋆ of sources present). We do this by relating
the workload of our model (to which we refer as Model i)
to the workload in a slightly different system (Model ii):
a model in which the generator Qb and the traffic rate matrices RI, RO,
and RA apply also when the buffer is empty (so Model
ii has no “idle mode”).The procedure of relating a feedback system (Modeli,
in which the sources react to the buffer content) to an (easier) nonfeedback
system (Model ii, in which the flows behave independently
of the buffer content) resembles that of  [7, Section 2].The distribution of the steady-state workload is
characterised in terms of the solution of a certain eigensystem (and a number
of additional linear equations). It also enables us to compute the
corresponding Laplace transform, which we use several times in the next
sections.


### 3.1. Preliminary Results

In this
subsection, we consider the modelwithout feedback, that is, Model
ii: the generator matrix Qb and the traffic rate matrices RI, RO, and RA apply not only when the buffer content is
positive, but also when the buffer is empty. We assume that the stability
criterion derived above (which reduces to ϱ<1/2 when N→∞) applies. Denote by Vt the buffer content of this system at time t (where V⋆ is its stationary version), and by Mt the number of flows present at time t (where M⋆ is its stationary version). Define also F→(x):=(F0(x),…,FN(x))T, where(3.1)Fn(x):=ℙ(V⋆≤x,⁢M⋆=n).Model ii
has been studied extensively in the literature; we now recall a number of basic
properties, which turn out to become useful when analysing Model
i, see Section 3.2.Buffer Content Distribution
It is well known from
the literature how theFn(x) can be determined; they obey the system of
linear differential equations:(3.2)RaF→′(x)=QbTF→(x).Owing to the special birth-death
structure, we can use explicit results obtained by van Doorn et al. [9].
A central role in the analysis is played by the
eigensystem(3.3)zv→Ra=v→Qb,with eigenvectors v→0 up to v→N and corresponding eigenvalues z0,…,zN. Notice that rA,n≠0,
for all n=0,…,N,
so RA is invertible. Then, [9, Theorem 1] says that all
eigenvalues z are real and simple.
Moreover, observe that the number of states ofMt in which Vt drains (or remains empty) is n+;
in the other N−n− the buffer level increases. Provided that the
stability condition is satisfied, [9, Theorem 1] entails that there are N−n− negative eigenvalues, one eigenvalue that
equals 0, and n− positive eigenvalues. Put the eigenvalues in increasing
order; let (v→m)n refer to the nth component of v→m.
Then the above, in conjunction with the fact that Fn(∞) lies between 0 and 1 for all n,
implies that in the representation(3.4)Fn(x)=∑m=0ncm⋅(v→m)n⋅eznx,the terms N−n++1 up to N are 0. As zN−n−=0 and v→N−n−=π,
it follows that the requirement F→(∞)=π implies cN−n−=1. We obtain(3.5)Fn(x)=πn+∑m=0n−n+cm⋅(v→m)n⋅ezmx.Now, only the cm (for m=0 up to N−n+) need to be determined. These follow from the
fact that Fn(0)=0 for n=n+,…,N.
These are N−n− equations in the same number of unknowns, and
can be determined explicitly in terms of z0,…,zN−n+,
as described in [9],
Section 4.Busy and Idle Periods
Elwalid and Mitra [10] give explicit expressions
for a number of quantities that are related to the busy and idle periods of the
queue. A busy period is, as before, defined as a period in which the buffer content
is positive, whereas an idle period is a period in which the buffer is empty.
It is easily seen that at the beginning of a busy period the number of flows
present is equal to n+;
at the end of the busy period the number of flows present is at most n−.
Denote byP→ the distribution of the number of flows
present at the end of the busy period. Let the matrices QbDD, QbDU, QbUD, QbUU be the submatrices that are obtained by
partitioning Qb into down-states (i.e., states n such that rA,n<0) and up-states (rA,n>0); similarly, F→(x) is partitioned in F→D(x) and F→U(x).
Then, it is not hard to prove that(3.6)P→=1〈F→D(0)QbDD,1〉⋅F→D(0)QbDD,see [10, Equation (5.9)]; 〈⋅,⋅〉 denotes the inner product of two vectors. The
mean idle period 𝔼I is given by(3.7)EI=−∑n=0n−Fn(0)〈F→D(0)QbDD,1〉.Finally, the mean busy period 𝔼B can be calculated. According to “renewal
reward” ∑n=0n−Fn(0)=𝔼I/(𝔼I+𝔼B), so that(3.8)EB=EI⋅1−∑n=0n−Fn(0)∑n=0n−Fn(0).


### 3.2. Analysis of Buffer Content Distribution

Now, we turn
back to Modeli, as described in Section 2.1. Our goal is
to show that the steady-state buffer content W⋆ of Model i (in which
there are different queueing dynamics when the buffer is empty or nonempty) is
intimately related to the steady-state buffer V⋆ content of Model ii (in
which there is no distinction between an empty and nonempty buffer).Let us start by making a number observations. First,
observe that in both Modelsi and ii
a busy period starts with n+ flow present. Also, the distribution of the
length of the busy period B is the same for both models, as well as the
distribution P→ of the number of flows present at the end of
the busy period. In other words, the difference between the models lies just in
the duration of the idle periods. In Section 3.1, we already found the mean
idle period 𝔼I of Model ii. Let us
therefore consider the mean idle period 𝔼J of Model i.As in [7, Lemma 2.3], we have that the mean idle period of
Model i equals(3.9)EJ=〈−P→Qi−1,1〉;the expected amount of time
during this idle time in which there are n flows present, say 𝔼(Jn),
is (−P→Qi−1)n,
that is, the nth entry of −P→Qi−1.
This follows from the fact that the mean time 𝔼m(Jn) spent in n during the idle time, given that at the
beginning of the idle time m flows were present, satisfies the linear
system, for m=1,…,n−,(3.10)Em(Jn)=(1λ+μc/2)1{m=n}+λλ+μc/2⋅Em+1(Jn)+μc/2λ+μc/2⋅Em−1(Jn),with 𝔼n+(Jn)=0,
and(3.11)E0(Jn)=1λ1{n=0}+E1(Jn).We now have collected all the required elements to
determine the distribution of the steady-state buffer content W⋆.
Analogously to [7,
Theorem 2.4] we obtain the following result for Gn(x):=ℙ(W⋆≤x,N⋆=n).Theorem 3.1.
For allx≥0,(3.12)Gn(x)=(Fn(x)−Fn(0))⋅EI+EBEJ+EB+E(Jn)EJ+EB.Proof.
This is proven as follows. We first condition onW⋆ being positive or zero:(3.13)Gn(x)=ℙ(W⋆≤x,N⋆=n∣W⋆>0)ℙ(W⋆>0)+ℙ(W⋆=0,N⋆=n).By applying the renewal-reward
theorem, the latter probability can be rewritten as(3.14)ℙ(W⋆=0,N⋆=n)=E(Jn)EJ+EB;notice that these probabilities
equal 0 for n=n+,…,N. Also, from “renewal-reward,”(3.15)ℙ(W⋆>0)ℙ(V⋆>0)=EI+EBEJ+EB.Hence, we are left with determining
the first probability in the right-hand side of (3.13). We first rewrite it
as(3.16)ℙ(N⋆=n∣W⋆>0)ℙ(W⋆>0)−ℙ(W⋆>x,N⋆=n∣W⋆>0)ℙ(W⋆>0).Now, recall that the
distribution of (W⋆,N⋆),
conditional on W⋆>0,
is the same as the distribution of (V⋆,M⋆),
conditional on V⋆>0. Combining this with (3.15), we
obtain(3.17)ℙ(M⋆=n∣V⋆>0)ℙ(W⋆>0)−ℙ(V⋆>x,M⋆=n∣V⋆>0)ℙ(W⋆>0)=(ℙ(V⋆>0,M⋆=n)−ℙ(V⋆>x,M⋆=n))×EI+EBEJ+EB=(Fn(x)−Fn(0))⋅EI+EBEJ+EB.This proves the claim.Upon combining the above theorem with representation
(3.5), we find the following useful result.Corollary 3.2.
For
allx≥0,
Theorem 3.1,
in conjunction with (3.5),
defines numbers αn and βnm (with n=0,…,N and m=0,…,N−n+) such that(3.18)Gn(x)=αn+∑m=0n−n+βnmezmx.Here, Gn(0)>0 if and only if n∈{0,…,n−}; Gn(0) is given by(3.19)Gn(0)=ℙ(W⋆=0,N⋆=n)=αn+∑m=0n−n+βnm.The probability of n flows in the system is given by ℙ(N⋆=n)=Gn(∞)=αn (where ∑n=0Nαn=1).
The Laplace transform of W⋆ reads, for s>0,(3.20)E(e−sW⋆1{N⋆=n})=αn+∑m=0n−n+sβnms−zm.


## 3.1. Preliminary Results

In this
subsection, we consider the modelwithout feedback, that is, Model
ii: the generator matrix Qb and the traffic rate matrices RI, RO, and RA apply not only when the buffer content is
positive, but also when the buffer is empty. We assume that the stability
criterion derived above (which reduces to ϱ<1/2 when N→∞) applies. Denote by Vt the buffer content of this system at time t (where V⋆ is its stationary version), and by Mt the number of flows present at time t (where M⋆ is its stationary version). Define also F→(x):=(F0(x),…,FN(x))T, where(3.1)Fn(x):=ℙ(V⋆≤x,⁢M⋆=n).Model ii
has been studied extensively in the literature; we now recall a number of basic
properties, which turn out to become useful when analysing Model
i, see Section 3.2.Buffer Content Distribution
It is well known from
the literature how theFn(x) can be determined; they obey the system of
linear differential equations:(3.2)RaF→′(x)=QbTF→(x).Owing to the special birth-death
structure, we can use explicit results obtained by van Doorn et al. [9].
A central role in the analysis is played by the
eigensystem(3.3)zv→Ra=v→Qb,with eigenvectors v→0 up to v→N and corresponding eigenvalues z0,…,zN. Notice that rA,n≠0,
for all n=0,…,N,
so RA is invertible. Then, [9, Theorem 1] says that all
eigenvalues z are real and simple.
Moreover, observe that the number of states ofMt in which Vt drains (or remains empty) is n+;
in the other N−n− the buffer level increases. Provided that the
stability condition is satisfied, [9, Theorem 1] entails that there are N−n− negative eigenvalues, one eigenvalue that
equals 0, and n− positive eigenvalues. Put the eigenvalues in increasing
order; let (v→m)n refer to the nth component of v→m.
Then the above, in conjunction with the fact that Fn(∞) lies between 0 and 1 for all n,
implies that in the representation(3.4)Fn(x)=∑m=0ncm⋅(v→m)n⋅eznx,the terms N−n++1 up to N are 0. As zN−n−=0 and v→N−n−=π,
it follows that the requirement F→(∞)=π implies cN−n−=1. We obtain(3.5)Fn(x)=πn+∑m=0n−n+cm⋅(v→m)n⋅ezmx.Now, only the cm (for m=0 up to N−n+) need to be determined. These follow from the
fact that Fn(0)=0 for n=n+,…,N.
These are N−n− equations in the same number of unknowns, and
can be determined explicitly in terms of z0,…,zN−n+,
as described in [9],
Section 4.Busy and Idle Periods
Elwalid and Mitra [10] give explicit expressions
for a number of quantities that are related to the busy and idle periods of the
queue. A busy period is, as before, defined as a period in which the buffer content
is positive, whereas an idle period is a period in which the buffer is empty.
It is easily seen that at the beginning of a busy period the number of flows
present is equal to n+;
at the end of the busy period the number of flows present is at most n−.
Denote byP→ the distribution of the number of flows
present at the end of the busy period. Let the matrices QbDD, QbDU, QbUD, QbUU be the submatrices that are obtained by
partitioning Qb into down-states (i.e., states n such that rA,n<0) and up-states (rA,n>0); similarly, F→(x) is partitioned in F→D(x) and F→U(x).
Then, it is not hard to prove that(3.6)P→=1〈F→D(0)QbDD,1〉⋅F→D(0)QbDD,see [10, Equation (5.9)]; 〈⋅,⋅〉 denotes the inner product of two vectors. The
mean idle period 𝔼I is given by(3.7)EI=−∑n=0n−Fn(0)〈F→D(0)QbDD,1〉.Finally, the mean busy period 𝔼B can be calculated. According to “renewal
reward” ∑n=0n−Fn(0)=𝔼I/(𝔼I+𝔼B), so that(3.8)EB=EI⋅1−∑n=0n−Fn(0)∑n=0n−Fn(0).


## 3.2. Analysis of Buffer Content Distribution

Now, we turn
back to Modeli, as described in Section 2.1. Our goal is
to show that the steady-state buffer content W⋆ of Model i (in which
there are different queueing dynamics when the buffer is empty or nonempty) is
intimately related to the steady-state buffer V⋆ content of Model ii (in
which there is no distinction between an empty and nonempty buffer).Let us start by making a number observations. First,
observe that in both Modelsi and ii
a busy period starts with n+ flow present. Also, the distribution of the
length of the busy period B is the same for both models, as well as the
distribution P→ of the number of flows present at the end of
the busy period. In other words, the difference between the models lies just in
the duration of the idle periods. In Section 3.1, we already found the mean
idle period 𝔼I of Model ii. Let us
therefore consider the mean idle period 𝔼J of Model i.As in [7, Lemma 2.3], we have that the mean idle period of
Model i equals(3.9)EJ=〈−P→Qi−1,1〉;the expected amount of time
during this idle time in which there are n flows present, say 𝔼(Jn),
is (−P→Qi−1)n,
that is, the nth entry of −P→Qi−1.
This follows from the fact that the mean time 𝔼m(Jn) spent in n during the idle time, given that at the
beginning of the idle time m flows were present, satisfies the linear
system, for m=1,…,n−,(3.10)Em(Jn)=(1λ+μc/2)1{m=n}+λλ+μc/2⋅Em+1(Jn)+μc/2λ+μc/2⋅Em−1(Jn),with 𝔼n+(Jn)=0,
and(3.11)E0(Jn)=1λ1{n=0}+E1(Jn).We now have collected all the required elements to
determine the distribution of the steady-state buffer content W⋆.
Analogously to [7,
Theorem 2.4] we obtain the following result for Gn(x):=ℙ(W⋆≤x,N⋆=n).Theorem 3.1.
For allx≥0,(3.12)Gn(x)=(Fn(x)−Fn(0))⋅EI+EBEJ+EB+E(Jn)EJ+EB.Proof.
This is proven as follows. We first condition onW⋆ being positive or zero:(3.13)Gn(x)=ℙ(W⋆≤x,N⋆=n∣W⋆>0)ℙ(W⋆>0)+ℙ(W⋆=0,N⋆=n).By applying the renewal-reward
theorem, the latter probability can be rewritten as(3.14)ℙ(W⋆=0,N⋆=n)=E(Jn)EJ+EB;notice that these probabilities
equal 0 for n=n+,…,N. Also, from “renewal-reward,”(3.15)ℙ(W⋆>0)ℙ(V⋆>0)=EI+EBEJ+EB.Hence, we are left with determining
the first probability in the right-hand side of (3.13). We first rewrite it
as(3.16)ℙ(N⋆=n∣W⋆>0)ℙ(W⋆>0)−ℙ(W⋆>x,N⋆=n∣W⋆>0)ℙ(W⋆>0).Now, recall that the
distribution of (W⋆,N⋆),
conditional on W⋆>0,
is the same as the distribution of (V⋆,M⋆),
conditional on V⋆>0. Combining this with (3.15), we
obtain(3.17)ℙ(M⋆=n∣V⋆>0)ℙ(W⋆>0)−ℙ(V⋆>x,M⋆=n∣V⋆>0)ℙ(W⋆>0)=(ℙ(V⋆>0,M⋆=n)−ℙ(V⋆>x,M⋆=n))×EI+EBEJ+EB=(Fn(x)−Fn(0))⋅EI+EBEJ+EB.This proves the claim.Upon combining the above theorem with representation
(3.5), we find the following useful result.Corollary 3.2.
For
allx≥0,
Theorem 3.1,
in conjunction with (3.5),
defines numbers αn and βnm (with n=0,…,N and m=0,…,N−n+) such that(3.18)Gn(x)=αn+∑m=0n−n+βnmezmx.Here, Gn(0)>0 if and only if n∈{0,…,n−}; Gn(0) is given by(3.19)Gn(0)=ℙ(W⋆=0,N⋆=n)=αn+∑m=0n−n+βnm.The probability of n flows in the system is given by ℙ(N⋆=n)=Gn(∞)=αn (where ∑n=0Nαn=1).
The Laplace transform of W⋆ reads, for s>0,(3.20)E(e−sW⋆1{N⋆=n})=αn+∑m=0n−n+sβnms−zm.


## 4. Queueing Delay Distribution

It is clear that it is a nontrivial step to translate
the steady-state workload distribution into the queueing delay distribution.
Importantly, to study the delay of a fluid particle arriving at time, say, 0,
the arrivals and departures of flowsafter 0 have impact. In Sections
4.1 and
4.2, we analyse the so-called virtual
queueing delay, that is, the delay experienced by a fluid particle arriving
at a random point in time (i.e., a “time average”); this is done through a
direct approach in Section 4.1 and through so-called “double transforms” is
Section 4.2. Section 4.3 characterises the queueing delay of an arbitrary fluid
particle (i.e., a “traffic average”).


### 4.1. Virtual Queueing Delay

LetD⋆ denote the delay experienced by a fluid
particle arriving at the queue in steady state, say for ease at time 0; this
type of delay is sometimes referred to as virtual queueing delay. Let O(0,t) denote the amount of output capacity available
in the interval [0,t).
If the fluid particle arrives at an empty queue, then the virtual delay is
clearly zero; if the fluid particle arrives at a nonempty queue, then the queue
is drained according to the rates rO,n until the particle has been served (in fact
even until the queue is empty). Define, for z≥0,
the random variable τz as the time until z units of service have become
available:(4.1)τz:=inf⁡{t≥0:O(0,t)=z}=inf{⁡t≥0:∫0tro,Nsds=z};notice that O(0,t) is increasing in t.
Then, analogously to [4, Section 4.1],  with some abuse of
notation,(4.2)Ee−sD⋆=∑n=0n∫0∞E(e−sτz∣N⋆=n)ℙ(W⋆=z,N⋆=n)dz.Hence, to further compute this
expression, we need to evaluate 𝔼(e−sτz∣N⋆=n).
Here, we can use [4, Proposition 4.1]:(4.3)E(e−sτz∣N⋆=n)=(exp⁡((Ro−1Qb−sRo−1)z)1)n,where 1 denotes an (N+1)-dimensional vector with 1s. As the same
proposition entails that the eigenvalues are simple and negative (hence real
numbers), it allows us to write, for constants γmn,
with m,n=0,…,N,(4.4)E(e−sτz∣N⋆=n)=∑m=0nγmneδm(s)z.We thus obtain the following
result.Theorem 4.1.
Fors>0,(4.5)Ee−sD⋆=∑n=0n∑m=0nγmnE(eδm(s)W⋆1{N⋆=n}),where the γmn are as in (4.4).
The δn(s),
for n=0,…,N,
are the eigenvalues of RO−1Qb−sRO−1 (which are negative). An expression for 𝔼(e−sW⋆1{N⋆=n}),
with s>0,
is available from Corollary 3.2.


### 4.2. A Second Approach: Double Transforms

We now proceed
with demonstrating a second approach, which relies on the concept of “double
transforms.” We feel that this is instructive, as this approach is used
extensively in the remainder of the paper (when analysing the flow transfer
delay and the sojourn time).Let us first condition on the buffer content (W0=x) that the fluid particle sees (say that it
arrives at time 0), and the number of flows that are then present (N0=n). Define, for given n=1,…,N and x≥0,
the transform of the queueing delay:(4.6)ξn(s∣x):=E(e−sD∣W0=x,N0=n).Then, we also introduce the
transform of ξn(s∣x) with respect to the workload x:(4.7)Kn(s,t)=∫0∞e−txξn(s∣x)dx;we say that Kn(s,t) is a “double transform.” Below, we show how to
use these double transforms to derive 𝔼e−sD⋆.Our first goal is to characterise theKn(s,t), n=1,…,N, for fixed s≥0 and t>0.
We do this by expressing Kn(s,t) in terms of Km(s,t) (with m≠n) as follows. Condition on the time until the
service rate changes; this time has an exponential distribution with mean νn−1.
Hence,(4.8)Kn(s,t)=∫0∞e−tx(∫0x/ro,nνne−νnye−sy(λ1{n<n}νnξn+1(s∣x−ro,ny)+μncνnξn−1(s∣x−ro,ny))dy+∫x/ro,n∞νne−νnye−sx/ro,ndy)dx.A straightforward change of
variable (x−rO,ny=z) then yields that(4.9)Kn(s,t)=1νn+s+tro,n(λ1{n<n}Kn+1(s,t)+μncKn−1(s,t)+ro,n).For given s and t,
these are N+1 linear equations in the same number of
unknowns. It is easy to see that the corresponding linear system is diagonally
dominant, and hence there is a unique solution. This enables us to find the Kn(s,t).Our second goal is to show how theseKn(s,t) yield an expression for 𝔼e−sD⋆.
At an arbitrary point in time, the distribution function of the workload
(jointly with the number of flows present) is given by Gn(⋅),
as given by (3.18). But, as the corresponding density is the weighted sum of
exponentials, it entails that knowledge of the Kn(s,t) gives an expression for the Laplace transform
of the virtual delay:(4.10)Ee−sD⋆=∑n=0n∫[0,∞)ξn(s∣x)dGn(x)=∑n=0n(ξn(s∣0)Gn(0)+∫(0,∞)∑m=0n−n−ξn(s∣x)zmβnmezmxdx)=∑n=0n−Gn(0)+∑n=0n∑m=0n−n−zmβnmKn(s,−zm);recall that zm<0 for m=0,…,N−n−.Theorem 4.2.
Fors>0,(4.11)Ee−sD⋆=∑n=0n−Gn(0)+∑n=0n∑m=0n−n−zmβnmKn(s,−zm).


### 4.3. “Packet Average” Queueing Delay

The previous
subsections presented expressions for the Laplace transform of the queueing
delay “at an arbitrary point in time” (a “time average”). Clearly, there is a
bias between the delayD⋆ “at an arbitrary point in time” and delay D¯⋆ “seen by an arbitrary fluid particle.” The
correction to be made is analogous to  [4, Section 4.2] and rather
straightforward:(4.12)Ee−sD¯⋆=∑n=0n(ri,n∑k=0nπkri,k)∑m=0nγmnE(eδm(s)W⋆1{N⋆=n}),compared to [11, Proposition 7.2].


## 4.1. Virtual Queueing Delay

LetD⋆ denote the delay experienced by a fluid
particle arriving at the queue in steady state, say for ease at time 0; this
type of delay is sometimes referred to as virtual queueing delay. Let O(0,t) denote the amount of output capacity available
in the interval [0,t).
If the fluid particle arrives at an empty queue, then the virtual delay is
clearly zero; if the fluid particle arrives at a nonempty queue, then the queue
is drained according to the rates rO,n until the particle has been served (in fact
even until the queue is empty). Define, for z≥0,
the random variable τz as the time until z units of service have become
available:(4.1)τz:=inf⁡{t≥0:O(0,t)=z}=inf{⁡t≥0:∫0tro,Nsds=z};notice that O(0,t) is increasing in t.
Then, analogously to [4, Section 4.1],  with some abuse of
notation,(4.2)Ee−sD⋆=∑n=0n∫0∞E(e−sτz∣N⋆=n)ℙ(W⋆=z,N⋆=n)dz.Hence, to further compute this
expression, we need to evaluate 𝔼(e−sτz∣N⋆=n).
Here, we can use [4, Proposition 4.1]:(4.3)E(e−sτz∣N⋆=n)=(exp⁡((Ro−1Qb−sRo−1)z)1)n,where 1 denotes an (N+1)-dimensional vector with 1s. As the same
proposition entails that the eigenvalues are simple and negative (hence real
numbers), it allows us to write, for constants γmn,
with m,n=0,…,N,(4.4)E(e−sτz∣N⋆=n)=∑m=0nγmneδm(s)z.We thus obtain the following
result.Theorem 4.1.
Fors>0,(4.5)Ee−sD⋆=∑n=0n∑m=0nγmnE(eδm(s)W⋆1{N⋆=n}),where the γmn are as in (4.4).
The δn(s),
for n=0,…,N,
are the eigenvalues of RO−1Qb−sRO−1 (which are negative). An expression for 𝔼(e−sW⋆1{N⋆=n}),
with s>0,
is available from Corollary 3.2.


## 4.2. A Second Approach: Double Transforms

We now proceed
with demonstrating a second approach, which relies on the concept of “double
transforms.” We feel that this is instructive, as this approach is used
extensively in the remainder of the paper (when analysing the flow transfer
delay and the sojourn time).Let us first condition on the buffer content (W0=x) that the fluid particle sees (say that it
arrives at time 0), and the number of flows that are then present (N0=n). Define, for given n=1,…,N and x≥0,
the transform of the queueing delay:(4.6)ξn(s∣x):=E(e−sD∣W0=x,N0=n).Then, we also introduce the
transform of ξn(s∣x) with respect to the workload x:(4.7)Kn(s,t)=∫0∞e−txξn(s∣x)dx;we say that Kn(s,t) is a “double transform.” Below, we show how to
use these double transforms to derive 𝔼e−sD⋆.Our first goal is to characterise theKn(s,t), n=1,…,N, for fixed s≥0 and t>0.
We do this by expressing Kn(s,t) in terms of Km(s,t) (with m≠n) as follows. Condition on the time until the
service rate changes; this time has an exponential distribution with mean νn−1.
Hence,(4.8)Kn(s,t)=∫0∞e−tx(∫0x/ro,nνne−νnye−sy(λ1{n<n}νnξn+1(s∣x−ro,ny)+μncνnξn−1(s∣x−ro,ny))dy+∫x/ro,n∞νne−νnye−sx/ro,ndy)dx.A straightforward change of
variable (x−rO,ny=z) then yields that(4.9)Kn(s,t)=1νn+s+tro,n(λ1{n<n}Kn+1(s,t)+μncKn−1(s,t)+ro,n).For given s and t,
these are N+1 linear equations in the same number of
unknowns. It is easy to see that the corresponding linear system is diagonally
dominant, and hence there is a unique solution. This enables us to find the Kn(s,t).Our second goal is to show how theseKn(s,t) yield an expression for 𝔼e−sD⋆.
At an arbitrary point in time, the distribution function of the workload
(jointly with the number of flows present) is given by Gn(⋅),
as given by (3.18). But, as the corresponding density is the weighted sum of
exponentials, it entails that knowledge of the Kn(s,t) gives an expression for the Laplace transform
of the virtual delay:(4.10)Ee−sD⋆=∑n=0n∫[0,∞)ξn(s∣x)dGn(x)=∑n=0n(ξn(s∣0)Gn(0)+∫(0,∞)∑m=0n−n−ξn(s∣x)zmβnmezmxdx)=∑n=0n−Gn(0)+∑n=0n∑m=0n−n−zmβnmKn(s,−zm);recall that zm<0 for m=0,…,N−n−.Theorem 4.2.
Fors>0,(4.11)Ee−sD⋆=∑n=0n−Gn(0)+∑n=0n∑m=0n−n−zmβnmKn(s,−zm).


## 4.3. “Packet Average” Queueing Delay

The previous
subsections presented expressions for the Laplace transform of the queueing
delay “at an arbitrary point in time” (a “time average”). Clearly, there is a
bias between the delayD⋆ “at an arbitrary point in time” and delay D¯⋆ “seen by an arbitrary fluid particle.” The
correction to be made is analogous to  [4, Section 4.2] and rather
straightforward:(4.12)Ee−sD¯⋆=∑n=0n(ri,n∑k=0nπkri,k)∑m=0nγmnE(eδm(s)W⋆1{N⋆=n}),compared to [11, Proposition 7.2].


## 5. Flow Transfer Delay Distribution

In this
section, we focus on the timeF it takes for an arbitrary arriving flow to
transmit its traffic. We define the flow transfer delay as the time between
arrival and the epoch that its last fluid particle has been transmitted into
the queue. Realise that the flow transfer time depends on the buffer content
and number of flows that the tagged flow sees upon arrival. Due to the
pasta-property, these coincide with the corresponding
time-averages. Recall that the case W=1,
as addressed in [4]
is simpler, as the buffer content seen upon arrival does not play a role.Let us first condition on the buffer content (W0=x) and the number of flows (N0=n). Define, for s≥0, t>0,
the Laplace transform of the flow transfer time (conditional on W0=x and N0=n) and its transform with respect to x:(5.1)ηn(s∣x):=E(e−sF∣W0=x,N0=n);Ln(s,t)=∫0∞e−txηn(s∣x)dx.For later reference, we also
introduce, for n=1,…,N and m>W,(5.2)κn(s,t):=νn+s−tra,n;tn(s):=νn+sra,n;ℓn,m(s):=Ln(s,tm(s)).Notice that, for n>W, tn(s) is positive.In Section5.1, we find the Ln(s,t);
this is in terms of auxiliary transforms that are determined in Section 5.2. We
conclude this section by presenting the transform of the flow transfer delay;
see Section 5.3.


### 5.1. A System of Equations for the Double Transform

We now deduce
a system of equations for theLn(s,t), n=1,…,N. We do so by distinguishing between
“down-states” (n with rA,n>0) and “up-states” (with rA,n<0). The idea is that for up-states during the
time till the first event (new arrival or departure) the buffer content cannot
become 0, while for down-states this is possible. As a consequence, these cases
have to be dealt with differently.Up-State
First, assume thatn is an “up-state”: n∈{n+,…,N}.
It is elementary to see that, conditioning on the first event taking place
after w units of time,(5.3)Ln(s,t)=∫0∞e−tx(∫0∞νne−νnwe−sw×(λ1{n<n}νnηn+1(s∣x+ra,nw)+n−1nμncνnηn−1(s∣x+ra,nw)+1nμncνn))dw⁢dx.This is the sum of three
integrals. The third equals(5.4)Ln(3)(s,t)=μncn⋅1t⋅1νn+s=:μn+(s,t).
Consider the first, and perform the change of variable x+rA,nw=y:(5.5)Ln(1)(s,t)=∫0∞e−tx(∫0∞λ1{n<n}e−νnwe−swηn+1(s∣x+ra,nw))dw⁢dx=∫0∞∫0y/ra,ne−t(y−ra,nw)λ1{n<n}e−νnwe−swηn+1(s∣y)dw⁢dy=∫0∞e−tyλ1{n<n}(e(t−νn/ra,n−s/ra,n)y−1tra,n−νn−s)ηn+1(s∣y)dy=λn(s,t)(Ln+1(s,t)−Ln+1(s,νn+sra,n)),whereλn(s,t):=λ1{n<n}κn(s,t).Similarly,(5.6)Ln(2)(s,t)=μn(s,t)(Ln−1(s,t)−Ln−1(s,νn+sra,n)),whereμn(s,t):=n−1nμncκn(s,t).We arrive at(5.7)Ln(s,t)=λn(s,t)(Ln+1(s,t)−ℓn+1,n(s))+μn(s,t)(Ln−1(s,t)−ℓn−1,n(s))+μn+(s,t).Later, it will turn out to be
also useful to consider the representation(5.8)κn(s,t)Ln(s,t)=λ1{n<n}⋅Ln+1(s,t)+n−1nμnc⋅Ln−1(s,t)+μncn⋅κn(s,t)t⋅1νn+s−λ1{n<n}⋅ℓn+1,n(s)−n−1nμnc⋅ℓn−1,n(s).Down-State
Now, assume thatn is a down-state: n∈{0,…,n−}.
In this case, we must distinguish between the cases that the process remains in
state n shorter, respectively, longer than −x/rA,n (which is a positive number); in the former
case, the buffer does not become empty before the first event, whereas in the
latter case it does. In more detail, we have(5.9)Ln(s,t)=∫0∞e−tx(∫0−x/ra,nνne−νnwe−sw×(λνnηn+1(s∣x+ra,nw)+n−1nμncνnηn−1(s∣x+ra,nw)+1nμncνn)dw+∫−x/ra,n∞νne−νnxesx/ra,nηn(s∣0)dw)dx.With μn−(s,t):=μnC/(n⋅κn(s,t)⋅t),
this simplifies to(5.10)Ln(s,t)=λn(s,t)Ln+1(s,t)+μn(s,t)Ln−1(s,t)+μn−(s,t)−ra,nκn(s,t)⋅ηn(s∣0).As indicated, our goal is to
generate a system of equations for the Ln(s,t) (with s and t fixed); we therefore wish to express ηn(s∣0) in terms of the Ln(s,t).
This can be done as follows.
(i)
First, forn∈{0,…,n−},(5.11)ηn(s∣0)=1νn+s(ληn+1(s∣0)+n−1nμnc⋅ηn−1(s∣0)+μncn),whereas for n∈{n+,…,N},(5.12)ηn(s∣0)=λ1{n<n}ra,n⋅Ln+1(s,νn+sra,n)+n−1n⋅μncra,n⋅Ln−1(s,νn+sra,n)+μncn⋅1νn+s=λ1{n<n}ra,n⋅ℓn+1,n(s)+n−1n⋅μncra,n⋅ℓn−1,n(s)+μncn⋅1νn+s.(ii)
Now, consider the vectorη→(s)≡(η1(s∣0),…,ηn−(s∣0))T. Define for n,m=1,…,n−,(5.13)an,m(s):={−λ,ifm=n+1;νn+s,ifm=n;−(n−1)/n⋅μnc,ifm=n+1.and 0 else. The corresponding
matrix is called A(s),
that is, A(s)≡(an,m(s))n,m=1n−;
for s>0 we have that A(s) is diagonally dominant and hence invertible.
Also, u→(s)≡(u1(s),…,un−(s))T,
with(5.14)un(s):={μnc/n,ifn=1,…,n−−1,μn−c/n−+ηn+(s∣0),ifn=n−.Then, (5.11) implies that η→(s)=(A(s))−1u→(s). In other words, once we know ηn+(s∣0),
we can compute η1(s∣0),…,ηn−(s∣0).(iii)
Letan,m−1(s):=(A(s))n,m−1.
Now, (5.12) entails that, for n=1,…,n−,(5.15)ηn(s∣0)=∑m=1n−−1an,m−1(s)μmcm+an,n−−1(s)(μn−cn−+ηn+(s∣0))=αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s), where(5.16)αn(s):=∑m=1n−an,m−1(s)μmcm+an,n−−1(s)μn−cn−⋅1νn−+s;βn(s):=an,n−−1(s)⋅λra,n+;γn(s):=an,n−−1(s)⋅n−n+⋅μn+cra,n+.
Inserting this into (5.10), we have found the following
relation for down-states n:(5.17)κn(s,t)Ln(s,t)=λ⋅Ln+1(s,t)+n−1nμnc⋅Ln−1(s,t)+μncnt−ra,n⋅(αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s));notice the similarity with the
equation for the up-states (5.8).


### 5.2. Determining the Auxiliary Transforms

From (5.8) and (5.17), it follows that for known functions ϕn(⋅,⋅) and ψm,k(⋅),(5.18)L→n(s,t)=∑m=1n((Q¯(s,t))−1)n,m(ϕm(s,t)+∑k=n+nψm,k(s)ℓm,k(s));here, the N×N matrix Q¯(s,t) is given through(5.19)(Q¯(s,t))n,m:={λ1{n<n},ifm=n+1;−κn(s,t),ifm=n;(n−1)/n⋅μnc,ifm=n−1.In other words, if the
transforms ℓn,m(s),
for n=1,…,N and m=n+,…,N,
would be known, then, for fixed s and t,
the values of the Ln(s,t) can be found directly from solving a system of
linear equations. The rest of this subsection is devoted to explaining how to
identify the ℓn,m(s),
for given s>0. We first prove a useful lemma.Lemma 5.1.
Consider,
for fixeds>0,
the t∈ℂ such that det(Q¯(s,t))=0. There are N−n−
such values such
that Ret>0.Proof.
First, rewriteQ¯=Q⌣−D+tRA,
with(5.20)Q˘n,m:={λ1{n≤n}ifm=n+1,−λ1{n≤n}−(n−1)/n⋅μncifm=n,(n−1)/n⋅μncifm=n−1,dn:=n−1⋅μnC+s>0, D:=diag{d→}.
Observe that Q⌣ is a generator matrix. Notice also that
solutions of det(Q¯(s,t))=0 are the eigenvalues of −RA−1(Q⌣−D).
(i)  We first focus on properties of the
eigenvalues of−RA−1Q⌣.
Recall that it follows from Theorem 2, part 3 of
[8]  that −RA−1Q⌣ has as many eigenvalues in the right half
plane as the number of up-states in RA,
that is, M:=N−n−;
there is also one eigenvalue of −RA−1Q⌣ equal to 0 (note that Q⌣ is singular), and the remaining n−−1 are in the left half plane. “Geršgorin” states
that all eigenvalues are in at least one of the disks(5.21){z∈ℂ:|z+Q⌣n,nra,n|≤|Q⌣n,nra,n|};the nth disk is a circle in the complex plane
around −Q⌣n,n/rA,n of radius |Q⌣n,n/rA,n| (which therefore goes through 0). Notice that Q⌣n,n<0 implies that the number of disks in the right
(left) half plane equals the number of up-states (down-states, resp.).
These observations are illustrated in the left panel of Figure 2.
(ii)  Now, consider the eigenvalues of−RA−1(Q⌣−εD) for small ε>0. Observe that these solve the equation ζ(ε,t):=det(Q⌣−εD+tRA)=0. As seen in the proof of [8, Theorem 2],(5.22)∂∂tζ(ε,t)|ε,t=0has the same sign as the mean
drift, that is, negative. Likewise, the derivative of ζ with respect to ε is positive (use that all diagonal entries of D are positive). Hence, replacing −RA−1Q⌣ by −RA−1(Q⌣−εD) moves the zero eigenvalue to the left.
“Geršgorin” implies
that all eigenvalues of−RA−1(Q⌣−εD) are in at least one of the
disks(5.23){z∈ℂ:|z+Q⌣n,n−εdnra,n|≤|Q⌣n,nra,n|};compared to the situation of ε=0 this means that the disks corresponding to the
up-states (that were in the right half plane) move to the right (with the same
radius); likewise, the disks corresponding to the down-states move to the left.
This implies that all eigenvalues in the left (right) half plane remain in the left (right) half plane,
because of the continuity of the solutions of ζ(ε,t) in the coefficients of the characteristic
polynomial.
Conclude that for smallε>0,
there are then N−n− in the right half plane, and the remaining n− in the left half plane. This is illustrated in
the right panel of Figure 2.
(iii) Observe that the same arguments imply this
classification remains valid when increasingε further. The special case ε=1 proves our claim.Eigenvalues in the complex  plane. (a)
Eigenvalues of−RA−1Q⌣;
(b) eigenvalues of −RA−1(Q⌣−εD).
In this example N=5,
and there are 3 up-states and 2 down-states.
(a)(b)Now, we are able to characterise the transformsℓn,m(s) as follows.Step 1.
Determine linear equations for the entries ofℓ→(s).
(a)
First,
focus on (5.8); these relate to n=n+,…,N.
We introduce the notation κn,m(s):=κn(s,tm(s)).
Then, for n≠m and m=n+,…,N,
we obtain by inserting t=tm(s):(5.24)κn,m(s)ℓn,m(s)=λ1{n≤n}(ℓn+1,m(s)−ℓn+1,n(s))+n−1nμnc(ℓn−1,m(s)−ℓn−1,n(s))+κn,m(s)μn+(s,tm(s)).Notice that κn,n(s)=κn(s,tn(s))=0,
so that for n=m both sides reduce to 0; hence, these M equations are meaningless.(b)
Now,
consider (5.17); these relate to n=1,…,n−.
Plugging in t=tm(s) for m=n+,…,N,
we obtain(5.25)κn,m(s)ℓn,m(s)=λℓn+1,m(s)+n−1nμncℓn−1,m(s)+μncntm(s)−ra,n⋅(αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s)).Step 2.
Reduce dimension of the vectorℓ→(s). The sets of  (5.24) and
(5.25) enable us to express the ℓn,m(s) for n=1,…,N and m=n+,…,N,
but n≠m,
in terms of ℓ→+(s):=(ℓn+,n+(s),…,ℓN,N(s))T.  We have thus identified functions φm(⋅,⋅) and ϑm,k(⋅) such that for n=1,…,N,(5.26)Ln(s,t)=∑m=1n((Q¯(s,t))−1)n,m(φm(s,t)+∑k=n+nϑm,k(s)ℓk,k(s)).In other words, when the N functions ℓn+,n+(s),…,ℓN,N(s) would be known, we would have found the Ln(s,t).Step 3.
Apply Lemma5.1. By virtue of Cramer's
rule, we obtain from (5.26) that(5.27)Ln(s,t)=det⁡Q¯n(s,t)det⁡Q¯(s,t),where Q¯n(s,t) is defined as Q¯(s,t) but with the nth column replaced by a vector of which the mth entry is φm(s,t)+∑k=n+Nϑm,k(s)ℓk,k(s).
For any t in the right half plane, this should have a
finite norm. Now fix s>0,
and use Lemma 5.1. Denote the zeroes of the denominator by τ1(s),…,τN−n−(s).
Conclude that each zero of the denominator should correspond to a zero of the
numerator. This yields N−n− equations that determine ℓn+,n+(s),…,ℓN,N(s).The above results are summarised in the following
theorem.Theorem 5.2.
Fors,t>0,(5.28)L→(s,t)=(Q¯(s,t))−1(φ→(s,t)+Θ(s)ℓ→+(s)),with ϕ→(s) and Θ(s) defined by (5.26).
For any s>0 there are N−n− values of t in the right half plane such that detQ¯(s,t)=0,
say τ1(s),…,τN−n−(s).
The vector ℓ→+(s) follows, for fixed s>0,
from letting t→τm(s) in the numerator of (5.27) and equating this to 0,
for m=1,…,N−n−.Remark 5.3.
It is easily verified that, in passing,
we have also found a procedure to
compute  theηn(s∣0),
for n=0,…,N,
compared to  (5.11), (5.12),
(5.15).


### 5.3. Flow Transfer Delay

It is clear
that, due topasta, the number of customers present at
(i.e., just after) arrival of an
(accepted) flow has distribution (n=1,…,N)(5.29)Hn:=Gn−1(∞)∑m=0n−1Gm(∞).For determining the flow
transfer delay, however, it is also necessary to know the amount of work found
in the buffer. The joint distribution of the number of flows and the buffer
content is given by(5.30)Hn(x):=Gn−1(x)∑m=0n−1Gm(∞)=A×(αn−1+∑m=0n−n−βn−1,mezmx),A:=(∑m=0n−1αm)−1;observe that indeed ∑n=1NHn(∞)=1. Hence,(5.31)Ee−sF=∑n=1n∫[0,∞)ηn(s∣x)dHn(x).Mimicking the derivation of 𝔼e−sD⋆,
we obtain the following result.Theorem 5.4.
Fors>0,(5.32)Ee−sF=∑n=1n+AGn−1(0)⋅ηn(s∣0)+∑n=1n∑m=0n−n−(Azmβn−1,m)⋅Ln(s,−zm).Expressions for Ln(⋅,⋅) and ηn(⋅∣0) are available from the previous subsection.


## 5.1. A System of Equations for the Double Transform

We now deduce
a system of equations for theLn(s,t), n=1,…,N. We do so by distinguishing between
“down-states” (n with rA,n>0) and “up-states” (with rA,n<0). The idea is that for up-states during the
time till the first event (new arrival or departure) the buffer content cannot
become 0, while for down-states this is possible. As a consequence, these cases
have to be dealt with differently.Up-State
First, assume thatn is an “up-state”: n∈{n+,…,N}.
It is elementary to see that, conditioning on the first event taking place
after w units of time,(5.3)Ln(s,t)=∫0∞e−tx(∫0∞νne−νnwe−sw×(λ1{n<n}νnηn+1(s∣x+ra,nw)+n−1nμncνnηn−1(s∣x+ra,nw)+1nμncνn))dw⁢dx.This is the sum of three
integrals. The third equals(5.4)Ln(3)(s,t)=μncn⋅1t⋅1νn+s=:μn+(s,t).
Consider the first, and perform the change of variable x+rA,nw=y:(5.5)Ln(1)(s,t)=∫0∞e−tx(∫0∞λ1{n<n}e−νnwe−swηn+1(s∣x+ra,nw))dw⁢dx=∫0∞∫0y/ra,ne−t(y−ra,nw)λ1{n<n}e−νnwe−swηn+1(s∣y)dw⁢dy=∫0∞e−tyλ1{n<n}(e(t−νn/ra,n−s/ra,n)y−1tra,n−νn−s)ηn+1(s∣y)dy=λn(s,t)(Ln+1(s,t)−Ln+1(s,νn+sra,n)),whereλn(s,t):=λ1{n<n}κn(s,t).Similarly,(5.6)Ln(2)(s,t)=μn(s,t)(Ln−1(s,t)−Ln−1(s,νn+sra,n)),whereμn(s,t):=n−1nμncκn(s,t).We arrive at(5.7)Ln(s,t)=λn(s,t)(Ln+1(s,t)−ℓn+1,n(s))+μn(s,t)(Ln−1(s,t)−ℓn−1,n(s))+μn+(s,t).Later, it will turn out to be
also useful to consider the representation(5.8)κn(s,t)Ln(s,t)=λ1{n<n}⋅Ln+1(s,t)+n−1nμnc⋅Ln−1(s,t)+μncn⋅κn(s,t)t⋅1νn+s−λ1{n<n}⋅ℓn+1,n(s)−n−1nμnc⋅ℓn−1,n(s).Down-State
Now, assume thatn is a down-state: n∈{0,…,n−}.
In this case, we must distinguish between the cases that the process remains in
state n shorter, respectively, longer than −x/rA,n (which is a positive number); in the former
case, the buffer does not become empty before the first event, whereas in the
latter case it does. In more detail, we have(5.9)Ln(s,t)=∫0∞e−tx(∫0−x/ra,nνne−νnwe−sw×(λνnηn+1(s∣x+ra,nw)+n−1nμncνnηn−1(s∣x+ra,nw)+1nμncνn)dw+∫−x/ra,n∞νne−νnxesx/ra,nηn(s∣0)dw)dx.With μn−(s,t):=μnC/(n⋅κn(s,t)⋅t),
this simplifies to(5.10)Ln(s,t)=λn(s,t)Ln+1(s,t)+μn(s,t)Ln−1(s,t)+μn−(s,t)−ra,nκn(s,t)⋅ηn(s∣0).As indicated, our goal is to
generate a system of equations for the Ln(s,t) (with s and t fixed); we therefore wish to express ηn(s∣0) in terms of the Ln(s,t).
This can be done as follows.
(i)
First, forn∈{0,…,n−},(5.11)ηn(s∣0)=1νn+s(ληn+1(s∣0)+n−1nμnc⋅ηn−1(s∣0)+μncn),whereas for n∈{n+,…,N},(5.12)ηn(s∣0)=λ1{n<n}ra,n⋅Ln+1(s,νn+sra,n)+n−1n⋅μncra,n⋅Ln−1(s,νn+sra,n)+μncn⋅1νn+s=λ1{n<n}ra,n⋅ℓn+1,n(s)+n−1n⋅μncra,n⋅ℓn−1,n(s)+μncn⋅1νn+s.(ii)
Now, consider the vectorη→(s)≡(η1(s∣0),…,ηn−(s∣0))T. Define for n,m=1,…,n−,(5.13)an,m(s):={−λ,ifm=n+1;νn+s,ifm=n;−(n−1)/n⋅μnc,ifm=n+1.and 0 else. The corresponding
matrix is called A(s),
that is, A(s)≡(an,m(s))n,m=1n−;
for s>0 we have that A(s) is diagonally dominant and hence invertible.
Also, u→(s)≡(u1(s),…,un−(s))T,
with(5.14)un(s):={μnc/n,ifn=1,…,n−−1,μn−c/n−+ηn+(s∣0),ifn=n−.Then, (5.11) implies that η→(s)=(A(s))−1u→(s). In other words, once we know ηn+(s∣0),
we can compute η1(s∣0),…,ηn−(s∣0).(iii)
Letan,m−1(s):=(A(s))n,m−1.
Now, (5.12) entails that, for n=1,…,n−,(5.15)ηn(s∣0)=∑m=1n−−1an,m−1(s)μmcm+an,n−−1(s)(μn−cn−+ηn+(s∣0))=αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s), where(5.16)αn(s):=∑m=1n−an,m−1(s)μmcm+an,n−−1(s)μn−cn−⋅1νn−+s;βn(s):=an,n−−1(s)⋅λra,n+;γn(s):=an,n−−1(s)⋅n−n+⋅μn+cra,n+.
Inserting this into (5.10), we have found the following
relation for down-states n:(5.17)κn(s,t)Ln(s,t)=λ⋅Ln+1(s,t)+n−1nμnc⋅Ln−1(s,t)+μncnt−ra,n⋅(αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s));notice the similarity with the
equation for the up-states (5.8).


## 5.2. Determining the Auxiliary Transforms

From (5.8) and (5.17), it follows that for known functions ϕn(⋅,⋅) and ψm,k(⋅),(5.18)L→n(s,t)=∑m=1n((Q¯(s,t))−1)n,m(ϕm(s,t)+∑k=n+nψm,k(s)ℓm,k(s));here, the N×N matrix Q¯(s,t) is given through(5.19)(Q¯(s,t))n,m:={λ1{n<n},ifm=n+1;−κn(s,t),ifm=n;(n−1)/n⋅μnc,ifm=n−1.In other words, if the
transforms ℓn,m(s),
for n=1,…,N and m=n+,…,N,
would be known, then, for fixed s and t,
the values of the Ln(s,t) can be found directly from solving a system of
linear equations. The rest of this subsection is devoted to explaining how to
identify the ℓn,m(s),
for given s>0. We first prove a useful lemma.Lemma 5.1.
Consider,
for fixeds>0,
the t∈ℂ such that det(Q¯(s,t))=0. There are N−n−
such values such
that Ret>0.Proof.
First, rewriteQ¯=Q⌣−D+tRA,
with(5.20)Q˘n,m:={λ1{n≤n}ifm=n+1,−λ1{n≤n}−(n−1)/n⋅μncifm=n,(n−1)/n⋅μncifm=n−1,dn:=n−1⋅μnC+s>0, D:=diag{d→}.
Observe that Q⌣ is a generator matrix. Notice also that
solutions of det(Q¯(s,t))=0 are the eigenvalues of −RA−1(Q⌣−D).
(i)  We first focus on properties of the
eigenvalues of−RA−1Q⌣.
Recall that it follows from Theorem 2, part 3 of
[8]  that −RA−1Q⌣ has as many eigenvalues in the right half
plane as the number of up-states in RA,
that is, M:=N−n−;
there is also one eigenvalue of −RA−1Q⌣ equal to 0 (note that Q⌣ is singular), and the remaining n−−1 are in the left half plane. “Geršgorin” states
that all eigenvalues are in at least one of the disks(5.21){z∈ℂ:|z+Q⌣n,nra,n|≤|Q⌣n,nra,n|};the nth disk is a circle in the complex plane
around −Q⌣n,n/rA,n of radius |Q⌣n,n/rA,n| (which therefore goes through 0). Notice that Q⌣n,n<0 implies that the number of disks in the right
(left) half plane equals the number of up-states (down-states, resp.).
These observations are illustrated in the left panel of Figure 2.
(ii)  Now, consider the eigenvalues of−RA−1(Q⌣−εD) for small ε>0. Observe that these solve the equation ζ(ε,t):=det(Q⌣−εD+tRA)=0. As seen in the proof of [8, Theorem 2],(5.22)∂∂tζ(ε,t)|ε,t=0has the same sign as the mean
drift, that is, negative. Likewise, the derivative of ζ with respect to ε is positive (use that all diagonal entries of D are positive). Hence, replacing −RA−1Q⌣ by −RA−1(Q⌣−εD) moves the zero eigenvalue to the left.
“Geršgorin” implies
that all eigenvalues of−RA−1(Q⌣−εD) are in at least one of the
disks(5.23){z∈ℂ:|z+Q⌣n,n−εdnra,n|≤|Q⌣n,nra,n|};compared to the situation of ε=0 this means that the disks corresponding to the
up-states (that were in the right half plane) move to the right (with the same
radius); likewise, the disks corresponding to the down-states move to the left.
This implies that all eigenvalues in the left (right) half plane remain in the left (right) half plane,
because of the continuity of the solutions of ζ(ε,t) in the coefficients of the characteristic
polynomial.
Conclude that for smallε>0,
there are then N−n− in the right half plane, and the remaining n− in the left half plane. This is illustrated in
the right panel of Figure 2.
(iii) Observe that the same arguments imply this
classification remains valid when increasingε further. The special case ε=1 proves our claim.Eigenvalues in the complex  plane. (a)
Eigenvalues of−RA−1Q⌣;
(b) eigenvalues of −RA−1(Q⌣−εD).
In this example N=5,
and there are 3 up-states and 2 down-states.
(a)(b)Now, we are able to characterise the transformsℓn,m(s) as follows.Step 1.
Determine linear equations for the entries ofℓ→(s).
(a)
First,
focus on (5.8); these relate to n=n+,…,N.
We introduce the notation κn,m(s):=κn(s,tm(s)).
Then, for n≠m and m=n+,…,N,
we obtain by inserting t=tm(s):(5.24)κn,m(s)ℓn,m(s)=λ1{n≤n}(ℓn+1,m(s)−ℓn+1,n(s))+n−1nμnc(ℓn−1,m(s)−ℓn−1,n(s))+κn,m(s)μn+(s,tm(s)).Notice that κn,n(s)=κn(s,tn(s))=0,
so that for n=m both sides reduce to 0; hence, these M equations are meaningless.(b)
Now,
consider (5.17); these relate to n=1,…,n−.
Plugging in t=tm(s) for m=n+,…,N,
we obtain(5.25)κn,m(s)ℓn,m(s)=λℓn+1,m(s)+n−1nμncℓn−1,m(s)+μncntm(s)−ra,n⋅(αn(s)+βn(s)ℓn++1,n+(s)+γn(s)ℓn−,n+(s)).Step 2.
Reduce dimension of the vectorℓ→(s). The sets of  (5.24) and
(5.25) enable us to express the ℓn,m(s) for n=1,…,N and m=n+,…,N,
but n≠m,
in terms of ℓ→+(s):=(ℓn+,n+(s),…,ℓN,N(s))T.  We have thus identified functions φm(⋅,⋅) and ϑm,k(⋅) such that for n=1,…,N,(5.26)Ln(s,t)=∑m=1n((Q¯(s,t))−1)n,m(φm(s,t)+∑k=n+nϑm,k(s)ℓk,k(s)).In other words, when the N functions ℓn+,n+(s),…,ℓN,N(s) would be known, we would have found the Ln(s,t).Step 3.
Apply Lemma5.1. By virtue of Cramer's
rule, we obtain from (5.26) that(5.27)Ln(s,t)=det⁡Q¯n(s,t)det⁡Q¯(s,t),where Q¯n(s,t) is defined as Q¯(s,t) but with the nth column replaced by a vector of which the mth entry is φm(s,t)+∑k=n+Nϑm,k(s)ℓk,k(s).
For any t in the right half plane, this should have a
finite norm. Now fix s>0,
and use Lemma 5.1. Denote the zeroes of the denominator by τ1(s),…,τN−n−(s).
Conclude that each zero of the denominator should correspond to a zero of the
numerator. This yields N−n− equations that determine ℓn+,n+(s),…,ℓN,N(s).The above results are summarised in the following
theorem.Theorem 5.2.
Fors,t>0,(5.28)L→(s,t)=(Q¯(s,t))−1(φ→(s,t)+Θ(s)ℓ→+(s)),with ϕ→(s) and Θ(s) defined by (5.26).
For any s>0 there are N−n− values of t in the right half plane such that detQ¯(s,t)=0,
say τ1(s),…,τN−n−(s).
The vector ℓ→+(s) follows, for fixed s>0,
from letting t→τm(s) in the numerator of (5.27) and equating this to 0,
for m=1,…,N−n−.Remark 5.3.
It is easily verified that, in passing,
we have also found a procedure to
compute  theηn(s∣0),
for n=0,…,N,
compared to  (5.11), (5.12),
(5.15).


## 5.3. Flow Transfer Delay

It is clear
that, due topasta, the number of customers present at
(i.e., just after) arrival of an
(accepted) flow has distribution (n=1,…,N)(5.29)Hn:=Gn−1(∞)∑m=0n−1Gm(∞).For determining the flow
transfer delay, however, it is also necessary to know the amount of work found
in the buffer. The joint distribution of the number of flows and the buffer
content is given by(5.30)Hn(x):=Gn−1(x)∑m=0n−1Gm(∞)=A×(αn−1+∑m=0n−n−βn−1,mezmx),A:=(∑m=0n−1αm)−1;observe that indeed ∑n=1NHn(∞)=1. Hence,(5.31)Ee−sF=∑n=1n∫[0,∞)ηn(s∣x)dHn(x).Mimicking the derivation of 𝔼e−sD⋆,
we obtain the following result.Theorem 5.4.
Fors>0,(5.32)Ee−sF=∑n=1n+AGn−1(0)⋅ηn(s∣0)+∑n=1n∑m=0n−n−(Azmβn−1,m)⋅Ln(s,−zm).Expressions for Ln(⋅,⋅) and ηn(⋅∣0) are available from the previous subsection.


## 6. Sojourn Time Distribution

In this
section, we study the sojourn timeS of flows in the system, which is defined as
the flow transfer time, increased by the time it takes to serve the last
particle of the flow. These components are not independent. Due to
pasta, the joint distribution of the workload and the
number of flows that a new (accepted) flow sees upon arrival, is given by Hn(⋅),
as defined through (5.30). To derive the Laplace transform of S,
we first need to describe the workload increment (which can be positive or
negative) during the flow transfer time F,
see Section 6.1. Then, in Section 6.2 we put the components together, and
derive the transform of S.


### 6.1. Joint Transform of Flow Transfer Time and Workload Increment

In the sequel
it will turn out that, in order to characterise the distribution of the sojourn
time, we do not just need the distribution ofF,
but rather its joint distribution with the workload WF+ and the number of flows present NF+ at the end of the transfer (not counting the
flow that just left). To this end we introduce the counterparts of ηn(s∣w) and Ln(s,t):(6.1)ηn,m(s→∣x):=E(e−s1F−s2WF+1{NF+=m}∣W0=x,N0=n);Ln,m(s→,t)=∫0∞e−txηn,m(s→∣x)dx,with s→≡(s1,s2). Similarly to before, we also define ℓn,m,k(s→):=Ln,m(s→,tk(s1)).Ln,m(s→,t) can be derived essentially in the same fashion
as Ln(s,t) was found in Sections 5.1 and
5.2. We sketch
this procedure. The counterpart of (5.7) is, for n>W,(6.2)Ln,m(s→,t)=λn(s1,t)(Ln+1,m(s→,t)−ℓn+1,m,n(s→))+μn(s1,t)(Ln−1,m(s→,t)−ℓn−1,m,n(s→))+μncn⋅1t+s1⋅1νn+s1+s2ra,n⋅1{n−1=m},whereas (5.10) generalises to, for n<W,(6.3)Ln,m(s→,t)=λn(s1,t)Ln+1,m(s→,t)+μn(s1,t)Ln−1,m(s→,t)+μncn⋅1κn(s1,t)⋅1t+s2⋅1{n−1=m}−ra,nκn(s1,t)⋅ηn,m(s→∣0).In the last equation, as before,
the ηn,m(s→∣0) can be expressed in terms of the ℓn,m,k(s→).
Then, it remains to find the transforms ℓn,m,k(s→);
these can be determined as in Section 5.2.


### 6.2. Sojourn Time

The sojourn
time can be decomposed into (i) the flow transfer delay and (ii) the time it
takes to serve the traffic that is in the buffer at the end of the flow
transfer delay (i.e., the time it takes to serve the last particle of the
tagged flow). This allows us to write, as
in [4, Section 6.3], with the usual abuse of notation(6.4)Ee−sS=E(e−sF−sτWF+)=∫0∞∫[0,∞)∑n=1n∑m=0n−1E(e−sF1{WF+=y,NF+=m}∣W0=x,N0=n)×E(e−sτy∣N0=m)dHn(x)dy.Consider the cases {W0=0} and {W0>0} separately. The contribution due to {W0=0} amounts to(6.5)∫0∞∑n=1n∑m=0n−1E(e−sF1{WF+=y,NF+=m}∣W0=0,N0=n)∑k=0nγkmeδk(s)yAGn−1(0)dy=∑n=1n∑m=0n−1∑k=0nγkmE(e−sF+δk(s)WF+1{NF+=m}∣W0=0,N0=n)AGn−1(0)=∑n=1n∑m=0n−1∑k=0nγkm⋅AGn−1(0)⋅ηn,m(s,−δk(s)∣0).Similarly, we find that the
contribution due to {W0>0} is(6.6)∫(0,∞)∑n=1n∑m=0n−1∑k=0nγkmηn,m(s,−δk(s)∣x)∑j=0n−n−(Azjβn−1,j)ezjxdx=∑n=1n∑m=0n−1∑k=0n∑j=0n−n−γkm⋅(Azjβn−1,j)⋅Ln,m(s,−δk(s),−zj).Theorem 6.1.
Fors>0, 𝔼e−sS is given by the sum of Expressions (6.5)
and (6.6).


## 6.1. Joint Transform of Flow Transfer Time and Workload Increment

In the sequel
it will turn out that, in order to characterise the distribution of the sojourn
time, we do not just need the distribution ofF,
but rather its joint distribution with the workload WF+ and the number of flows present NF+ at the end of the transfer (not counting the
flow that just left). To this end we introduce the counterparts of ηn(s∣w) and Ln(s,t):(6.1)ηn,m(s→∣x):=E(e−s1F−s2WF+1{NF+=m}∣W0=x,N0=n);Ln,m(s→,t)=∫0∞e−txηn,m(s→∣x)dx,with s→≡(s1,s2). Similarly to before, we also define ℓn,m,k(s→):=Ln,m(s→,tk(s1)).Ln,m(s→,t) can be derived essentially in the same fashion
as Ln(s,t) was found in Sections 5.1 and
5.2. We sketch
this procedure. The counterpart of (5.7) is, for n>W,(6.2)Ln,m(s→,t)=λn(s1,t)(Ln+1,m(s→,t)−ℓn+1,m,n(s→))+μn(s1,t)(Ln−1,m(s→,t)−ℓn−1,m,n(s→))+μncn⋅1t+s1⋅1νn+s1+s2ra,n⋅1{n−1=m},whereas (5.10) generalises to, for n<W,(6.3)Ln,m(s→,t)=λn(s1,t)Ln+1,m(s→,t)+μn(s1,t)Ln−1,m(s→,t)+μncn⋅1κn(s1,t)⋅1t+s2⋅1{n−1=m}−ra,nκn(s1,t)⋅ηn,m(s→∣0).In the last equation, as before,
the ηn,m(s→∣0) can be expressed in terms of the ℓn,m,k(s→).
Then, it remains to find the transforms ℓn,m,k(s→);
these can be determined as in Section 5.2.


## 6.2. Sojourn Time

The sojourn
time can be decomposed into (i) the flow transfer delay and (ii) the time it
takes to serve the traffic that is in the buffer at the end of the flow
transfer delay (i.e., the time it takes to serve the last particle of the
tagged flow). This allows us to write, as
in [4, Section 6.3], with the usual abuse of notation(6.4)Ee−sS=E(e−sF−sτWF+)=∫0∞∫[0,∞)∑n=1n∑m=0n−1E(e−sF1{WF+=y,NF+=m}∣W0=x,N0=n)×E(e−sτy∣N0=m)dHn(x)dy.Consider the cases {W0=0} and {W0>0} separately. The contribution due to {W0=0} amounts to(6.5)∫0∞∑n=1n∑m=0n−1E(e−sF1{WF+=y,NF+=m}∣W0=0,N0=n)∑k=0nγkmeδk(s)yAGn−1(0)dy=∑n=1n∑m=0n−1∑k=0nγkmE(e−sF+δk(s)WF+1{NF+=m}∣W0=0,N0=n)AGn−1(0)=∑n=1n∑m=0n−1∑k=0nγkm⋅AGn−1(0)⋅ηn,m(s,−δk(s)∣0).Similarly, we find that the
contribution due to {W0>0} is(6.6)∫(0,∞)∑n=1n∑m=0n−1∑k=0nγkmηn,m(s,−δk(s)∣x)∑j=0n−n−(Azjβn−1,j)ezjxdx=∑n=1n∑m=0n−1∑k=0n∑j=0n−n−γkm⋅(Azjβn−1,j)⋅Ln,m(s,−δk(s),−zj).Theorem 6.1.
Fors>0, 𝔼e−sS is given by the sum of Expressions (6.5)
and (6.6).


## 7. Discussion and Concluding Remarks

In this paper,
we have considered a relay node in an ad hoc network, fed by a Poisson stream
of exponentially distributed jobs. We have characterised its performance in
terms of (the Laplace transforms of) the buffer content, the queueing delay,
the flow transfer delay, and the sojourn time.Integer Weights
In our analysis, we throughout assumed that the
weightW was noninteger. If W∈ℕ,
the analysis is slightly more involved. We now indicate how the analysis should
be adapted. In the first place, one of the coupled differential equations in
(3.2) has left-hand side 0 (because rA,W=0); if we enumerate the equations 0,…,W,
then the Wth equation reads(7.1)0=∑n=0nqn,wFn(x),orFw(x)=−1qw,w∑n≠wqn,wFn(x),where qn,m is the (n,m)-entry of Qb.
Then the nth differential equation (where n≠W) becomes(7.2)ra,nFn′(x)=∑m≠wqm,n−wFm(x),whereqm,n−w:=(qm,n−qw,nqw,w⋅qm,w).Interestingly, qm,n−W≥0 for m≠n,
and ∑n≠Wqm,n−W=0,
and hence the Qb−W:=(qm,n−W)m,n (with m,n=0,…,N,
but m,n≠W) correspond to an N-dimensional generator matrix. In self-evident
notation, we have arrived at(7.3)Ra−w(F→−w)′(x)=(Qb−w)TF→−w(x),with all entries of RA−W not equal to 0. In this way, the steady-state
buffer content distribution of Model ii can be determined:
we first find the Fn(⋅) for n≠W,
and then we use (7.1) to derive FW(⋅).
The buffer content distribution follows in the same
fashion as in Section3.2. It is of crucial importance to choose a definition
of busy and idle periods (i.e., one has to choose whether periods with an empty
buffer and W flows in the system belong to the busy or to
the idle periods), and then to consistently use this definition. It is readily
verified that for the other performance measures no specific problems arise.Limiting Cases
We now consider a number of interesting limiting
choices for the weightW. For ease, we lift the assumption of performing
admission control; in other words: we take N=∞. Let N⋆(W) be the steady-state number of flows in the
system (i.e., transmitting traffic into the queue), for a given weight W; W⋆(W) is defined analogously.As argued in [3], the total amount of traffic in the system has the same dynamics as an M/M/1 queue; in
this M/M/1 all job sizes are inflated by a factor 2 (as they have to be
processed twice). This total amount of traffic is to be understood as 2Ws⋆(W)+W⋆(W), Ws⋆(W) denoting the traffic at the sources, and W⋆(W) the traffic at the queue; the factor 2 is due
to the fact that traffic at the sources still needs to be processed twice, as
opposed to traffic at the queue. Importantly, the evolution of the total queue
is independent of W;
realise that the total queue is work-conserving. It follows from the
Pollaczek-Khinchine formula that the mean amount of work (measured in
processing time) in the system is 4λ−1⋅ϱ2/(1−2ϱ),
independently of the choice of W. It can be seen that pathwise the amount of traffic that is at
the sources increases in W (as the rI,n decrease), so that the amount of traffic in
the queue decreases in W.(i) In caseW=∞,
the queue has maximum weight, and never builds up. Always half of the capacity
is dedicated to the flows, and the other half to the queue. It can be verified
that the queue becomes a normal processor-sharing queue. The flow transfer
delays and the sojourn times coincide. Elementary computations reveal that the
steady-state distribution of the number of flows in the system is geometrically
distributed. The probability of n flows in the system is (1−2ϱ)(2ϱ)n;
the mean number of flows is 𝔼N⋆(∞)=2ϱ/(1−2ϱ).(ii) A second extreme case isW=0. Then, the queue is only served when the flows
have transmitted all their traffic into the queue; when the flows have
something to transmit, the queue grows at a rate C. In this case, the probability of n flows in the system is (1−ϱ)ϱn;
the mean number of flows is 𝔼N⋆(0)=ϱ/(1−ϱ).Evidently, the mean amount of work in the queuedoes depend on W. First, observe that the mean amount of traffic
(measured in processing time) that the flows still need to inject into the
queue is(7.4)2EN⋆(w)⋅1μ,where the factor 2 reflects the
fact that the traffic still needs to be processed twice. We conclude that the
mean amount of work (measured in processing time) at the queue
is(7.5)EW⋆(∞)=c(1λ⋅4ϱ21−2ϱ−2μc⋅2ϱ1−2ϱ)=0;EW⋆(0)=c(1λ⋅4ϱ21−2ϱ−2μc⋅ϱ1−ϱ)=1μ⋅2ϱ(1−ϱ)(1−2ϱ).These arguments can be used to
quantify the tradeoff between the flow transfer delay (which increases in W) and the queueing delay of the last particle
(which decreases in W). The formula for W=1 was already given in [3](7.6)EW⋆(1)=1μ⋅4ϱ2(1−ϱ)(1−2ϱ),which is indeed between (7.5) for ϱ<1/2.Numerical Aspects
The numerical techniques to
be used in the approach presented in this paper (namely, solving eigensystems,
solving linear systems, and numerically inversion of Laplace transforms) are
techniques that are well established, and for which efficient and reliable
computer code is available. The distributions ofW⋆ and D⋆ can be found without the need of performing
any Laplace inversion; as a result, their evaluation boils down to solving a
linear system of differential equations, comparable to those highlighted in,
for example, [5, 6]. Determining the distribution of F and S does require Laplace inversion. It is noted
that recently, substantial progress has been made with respect to this type of
inversion techniques. Besides the “classical” reference [12],
we wish to draw attention
to significant recent progress by den Iseger, as reported on in [13]; the latter reference also
provides a fairly complete literature overflow.When focusing on obtaining numerical output, there are
several alternatives to the approach described in this paper. A first
alternative is to rely on tail asymptotics, as done in [4] for W=1,
in the spirit of(7.7)lim⁡x→∞1xlog⁡ℙ(W⋆>x)=−ϑ⋆,where the constant ϑ⋆>0 follows from the solution of the corresponding
eigensystem; likewise, one could consider the logarithmic tail asymptotics of D⋆,F,
and S.
As this just provides us with decay rate (i.e., it says that ℙ(W⋆>x)=ψ(x)e−ϑ⋆x for an unknown function ψ(⋅) such that x−1⋅logψ(x)→0), one could use importance sampling-based
simulations to improve on this, where the twisted distribution can be computed
as in [14]. For
sojourn times, such importance sampling schemes can be set up as in [15].Subjects for Future Research
We mention the following directions for further research.
(i)Multiple bottlenecks. In some situations,
the scenario of a single bottleneck link may be an oversimplification of
reality, and in such cases one could study multiple bottleneck links that share capacity.
The complicating factor is that then the dynamics of the flows feeding into one
queue will be affected by the workload process in other queues; the queues
cannot be analysed separately. This gives the model the flavor of
coupled-processors systems as studied in, for example, [16].
(ii)Other flow-size distributions. Another
challenging extension is to consider nonexponential flow sizes; particularly
the impact heavy-tailed jobs is interesting to study. Suppose for instance that
the flow sizes have a regularly varying distribution of index −δ,
then it is an open question whether the sojourn times are regularly varying of
index 1−δ,
as is the case in the M/G/1 fifo queue [17, 18],
or regularly varying of
index −δ,
as is the case in the M/G/1 processor sharing (ps) queue
[19] (or perhaps
regularly varying of another index).
(iii)Other queueing disciplines. In this study,
as well as in [3, 4],
the queue was supposed to operate under the
fifo discipline. This introduces some “unfairness” in that
even small jobs can incur significant delay. Put differently the sojourn time
of a job of size x,
say S(x),
is such that limx↓0𝔼S(x)>0. If the scheduling discipline in the queue
would be ps (rather than fifo), then
this limit would be 0; in this sense, ps could be regarded
as a remedy for unfairness.
(iv)Weight selection. Now, that we are able to
evaluate the performance of the relay node for a given weight W,
one may wonder what value of W should be chosen. As argued above, there is a
tradeoff between the flow transfer delay and the queueing delay of the last
particle, imposing some cost structure, and optimal value for W can be selected.In a network setting, each node chooses its own
weight. A high weight may be beneficial for the node itself, but harmful for
other nodes. In view of this it may make sense to charge nodes for their
weight. Pricing schemes could provide incentives for users to act as transit
nodes on multihop paths, compared to [20, 21].



---

*Source: 518214-2008-07-13.xml*