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

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 are 𝑛 users who wish to transmit data through a specific node, each of them obtains a share 1/(𝑛+𝑊) of the service capacity to feed traffic into the queue of the node, whereas the remaining fraction 𝑊/(𝑛+𝑊) is used to serve the queue; here 𝑊>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 𝑊=1 has been addressed before; the present paper focuses on the intrinsically harder case 𝑊>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.

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 as queues: 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 are 𝑛 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/(𝑛+1) of the available medium capacity, which is the same fraction as is allocated to each of the 𝑛 “sending nodes.” In other words, as soon as 𝑛>1, the node's input rate exceeds its output rate, and hence the excess traffic accumulates in the node's buffer; only when 𝑛=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 𝑊>0, such that each of the 𝑛 sending nodes obtains a fraction 1/(𝑛+𝑊) of the capacity, while the remaining 𝑤/(𝑛+𝑊) is allocated to serving the queue. Clearly, when 𝑊>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 𝑛<𝑊 (rather than just for 𝑛=0). The price paid is that the traffic remains longer at the sending nodes.

There is one important modelling feature, that applies to the case 𝑊>1, that needs to be mentioned here. When the queue is empty and the number of sending nodes 𝑛 is below 𝑊, it does not make sense to assign each of the sending notes just a share 1/(𝑛+𝑊): it would imply that the (available) service rate is strictly larger than the input rate, and that a fraction (𝑊−𝑛)/(𝑛+𝑊)>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 𝑛<𝑊, then half of the capacity is allocated to serving the queue of the relay node, whereas the other half is shared equally among the 𝑛 sending nodes (such that the input and output rates are equal, the queue remains empty, and all available capacity is used). Notice that when 𝑊≤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 𝑛<𝑊 jobs in the system.

The case 𝑊=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 case 𝑊=1 (and in fact also 𝑊<1) has nice features that are lost for 𝑊>1. As a consequence, the analysis of [4] for 𝑊=1 does not carry over to 𝑊>1. The main differences are the following.

(1) In the first place, as we explained, the case 𝑊=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 𝑊>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 for 𝑊=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 𝑊>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 case 𝑊=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 𝑊>1, the periods of net output are not exponentially distributed, and therefore we lose this nice feature.

The above explains that the analysis for 𝑊>1 is considerably more involved than for 𝑊=1.

For various values of 𝑊, 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 𝑊>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 Section 2—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 for 𝑊=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].

Section 6 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 a FIFO manner) and leaves when ready. When there are 𝑛 flows active and the queue is nonempty, each flow transmits traffic into the queue at rate 𝐶/(𝑛+𝑊), while a rate 𝑊𝐶/(𝑛+𝑊) is used to serve the queue; as a consequence, the queue drains when the number of flows present is below 𝑊. For ease, we assume that 𝑊 is noninteger; we come back to this issue later in Section 7. When there are 𝑛 flows active and the queue is empty, all flows transmit at rate 𝐶/(2𝑛), while the queue is served at rate 𝐶/2, so that the queue remains empty.

Suppose that we impose the admission control policy that the system accommodates maximally 𝑁∈ℕ flows simultaneously; in this way, each active flow is guaranteed at least a transmission rate 𝐶/(𝑁+𝑊), and the queue at least 𝑊𝐶/(𝑁+𝑊). We assume that 𝑁>𝑊 (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 time 𝑡 by 𝑊_𝑡. We let 𝑁_𝑡 denote the number of flows present (i.e., feeding traffic into the queue) at time 𝑡. A pictorial illustration is given in Figure 1.

(a)

(b)

(c)

Figure 1 

Illustration of the model. The top panel depicts the number of flows 𝑁𝑡 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 𝑊𝑡 as a function of time; we have chosen 𝑊=3/2. Note that when 𝑁𝑡=1 the rate per source equals 𝐶/2 when the buffer is empty, and 2𝐶/5 when the buffer is nonempty. Also note that the queue builds up when 𝑁𝑡≥2.

We introduce the following notation.

(i) 𝑊_𝑡>0: the “busy mode.” It is not hard to see, under the assumption of exponentially distributed flow sizes (with mean 𝜇⁻¹) and interarrival times (with mean 𝜆⁻¹), that, during periods that 𝑊_𝑡>0, the process 𝑁_𝑡 behaves as a Markov chain on {0,…,𝑁}, with generator matrix 𝑄_b∶=⎛⎜⎜⎜⎜⎜⎜⎜⎝𝜇−𝜆𝜆¹𝑐−𝜇𝜇₁𝑐−𝜆𝜆₂𝑐−𝜇⋱₂𝑐−𝜆⋱𝜆𝜇⋱_𝑁𝑐−𝜇_𝑁𝑐⎞⎟⎟⎟⎟⎟⎟⎟⎠,(2.1)where 𝜇_𝑛∶=𝜇𝑛/(𝑛+𝑊); the subscript “b” stands for “busy.” We define 𝜈_𝑛∶=𝜇_𝑛𝐶+𝜆1{𝑛<𝑁}.

When 𝑁_𝑡=𝑛, the aggregate traffic rate generated by the flows is 𝑟_𝐼,𝑛∶=𝐶𝑛/(𝑛+𝑊), while the queue's output rate is 𝑟_𝑂,𝑛∶=𝑊𝐶/(𝑛+𝑊), such that the net rate of change of the queue is𝑟_𝐴,𝑛∶=𝑟_𝐼,𝑛−𝑟_𝑂,𝑛=𝑐𝑛−𝑊𝑛+𝑊.(2.2)Define 𝑅_𝐼∶=diag{𝑟_𝐼}, 𝑅_𝑂∶=diag{𝑟_𝑂}, and 𝑅_𝐴∶=𝑅_𝐼−𝑅_𝑂.

Busy periods are periods in which 𝑊_𝑡 is positive all the time. With 𝑛₊∶=⌈𝑊⌉, it is evident that the number of active flows at the beginning of a busy period equals 𝑛₊. The number of active flows at the end of a busy period is in {0,…,𝑛₋}, with 𝑛₋∶=⌊𝑊⌋.

(ii) 𝑊_𝑡=0: the “idle mode.” Let idle periods be periods in which 𝑊_𝑡=0 all the time. An idle period ends as soon as 𝑁_𝑡=𝑛₊. During the idle period necessarily 𝑁_𝑡∈{0,…,𝑛₋}. One could say that 𝑁_𝑡 behaves as a Markov chain on {0,…,𝑛₋} until 𝑁_𝑡 jumps from 𝑛₋ to 𝑛₊ (i.e., the start of a new busy period). The corresponding rate matrix (which is not a bona fide generator matrix) is𝑄_i∶=⎛⎜⎜⎜⎜⎜⎜⎜⎝𝜇−𝑐𝜆/2−𝜇𝜇𝑐/𝑐𝜆2/2−𝜆−𝜇𝑐⋱/𝜆2−𝜆⋱𝜆𝜇⋱𝑐/2−𝜇𝑐/2−𝜆⎞⎟⎟⎟⎟⎟⎟⎟⎠;(2.3)the subscript “i” stands for “idle.”

2.2. Stability Condition

To make sure that the steady-state workload 𝑊^⋆ is finite a.s., the mean drift should be negative when 𝑊_𝑡 is large. Since 𝑁_𝑡 behaves essentially like a stationary Markov process with generator 𝑄_b when 𝑊_𝑡 is large, it follows that 𝑊_𝑡 can only escape to ∞ when ∑^𝑁𝑛=0𝜋_𝑛𝑟_𝐴,𝑛≥0, denoting by ^→𝜋≡(𝜋₀,…,𝜋_𝑁)^T the invariant measure of 𝑄_b. Hence, we should require that ∑^𝑁𝑛=0𝜋_𝑛𝑟_𝐴,𝑛<0. Elementary Markov chain analysis yields that, with 𝜚∶=𝜆/(𝜇𝐶),𝜋_𝑛=𝜌^𝑛𝔹(𝑛+𝑊,𝑛)∑^𝑁𝑚=0𝜚^𝑚𝔹(𝑚+𝑊,𝑚),(2.4)where, for 𝑘≥ℓ≥0,𝔹(⋅,⋅) can be regarded as a “generalised binomial coefficient”:𝔹(𝑘,ℓ)∶=Γ(𝑘+1)Γ(ℓ+1)Γ(𝑘−ℓ+1).(2.5)The stability condition becomes ∑^𝑁𝑛=0𝑟_𝐴,𝑛𝜋_𝑛<0.

It is instructive to show how this condition simplifies in the situation of 𝑁→∞. Due to (recognize the probability density function of the negative binomial distribution)^∞_𝑚=0𝜚^𝑚𝔹(𝑚+𝑤,𝑚)=1(1−𝜚)^𝑊+1,(2.6)we have to verify whether^∞_𝑛=0𝜋_𝑛⋅𝑛−𝑊𝑛+𝑊=^∞_𝑛=0𝜌^𝑛𝔹(𝑛+𝑤,𝑛)(1−𝜚)^𝑊+1⋅𝑛−𝑊𝑛+𝑊<0.(2.7)Using identity (2.6) again, writing𝑛−𝑊𝑛+𝑊=𝑛+2𝑛𝑊−1,(2.8)and observing that^∞_𝑛=0𝜚^𝑛𝔹(𝑛+𝑤,𝑛)2𝑛𝑛+𝑊=^∞_𝑛=1𝜚^𝑛𝔹(𝑛+𝑤,𝑛)2𝑛=2_∞𝑛+𝑊_𝑛=1𝜚^𝑛𝔹(𝑛+𝑤−1,𝑛−1)=2𝜚^∞_𝑛=0𝜚^𝑛𝔹(𝑛+𝑤,𝑛)=2𝜚(1−𝜚)^𝑊+1,(2.9)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 𝑊). 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 workload 𝑊^⋆ of the queue introduced in the previous section (jointly with the steady-state number 𝑁^⋆ 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 𝑄_b and the traffic rate matrices 𝑅_𝐼, 𝑅_𝑂, and 𝑅_𝐴 apply also when the buffer is empty (so Model II has no “idle mode”).

The procedure of relating a feedback system (Model I, 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 model without feedback, that is, Model II: the generator matrix 𝑄_b and the traffic rate matrices 𝑅_𝐼, 𝑅_𝑂, and 𝑅_𝐴 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 𝑁→∞) applies. Denote by 𝑉_𝑡 the buffer content of this system at time 𝑡 (where 𝑉^⋆ is its stationary version), and by 𝑀_𝑡 the number of flows present at time 𝑡 (where 𝑀^⋆ is its stationary version). Define also 𝐹^→(𝑥)∶=(𝐹₀(𝑥),…,𝐹_𝑁(𝑥))^T, where𝐹_𝑛(𝑥)∶=ℙ𝑉^⋆≤𝑥,𝑀^⋆=𝑛.(3.1)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 the 𝐹_𝑛(𝑥) can be determined; they obey the system of linear differential equations:𝑅𝐹_𝐴→^(𝑥)=𝑄𝐹^Tb→(𝑥).(3.2)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𝑧^→𝑣𝑅_𝐴=^→𝑣𝑄_b,(3.3)with eigenvectors ^→𝑣₀ up to ^→𝑣_𝑁 and corresponding eigenvalues 𝑧₀,…,𝑧_𝑁. Notice that 𝑟_𝐴,𝑛≠0, for all 𝑛=0,…,𝑁, so 𝑅_𝐴 is invertible. Then, [9, Theorem 1] says that all eigenvalues 𝑧 are real and simple.
Moreover, observe that the number of states of 𝑀_𝑡 in which 𝑉_𝑡 drains (or remains empty) is 𝑛₊; in the other 𝑁−𝑛₋ the buffer level increases. Provided that the stability condition is satisfied, [9, Theorem 1] entails that there are 𝑁−𝑛₋ negative eigenvalues, one eigenvalue that equals 0, and 𝑛₋ positive eigenvalues. Put the eigenvalues in increasing order; let (^→𝑣_𝑚)_𝑛 refer to the 𝑛th component of ^→𝑣_𝑚. Then the above, in conjunction with the fact that 𝐹_𝑛(∞) lies between 0 and 1 for all 𝑛, implies that in the representation𝐹_𝑛(𝑥)=^𝑁_𝑚₌₀𝑐_𝑚⋅^→𝑣_𝑚_𝑛⋅𝑒^𝑧𝑛𝑥,(3.4)the terms 𝑁−𝑛₊+1 up to 𝑁 are 0. As 𝑧_{𝑁−𝑛−}=0 and ^→𝑣_{𝑁−𝑛−}=𝜋, it follows that the requirement 𝐹^→(∞)=𝜋 implies 𝑐_{𝑁−𝑛−}=1. We obtain𝐹_𝑛(𝑥)=𝜋_𝑛+^𝑁−𝑛+_𝑚₌₀𝑐_𝑚⋅^→𝑣_𝑚_𝑛⋅𝑒^𝑧𝑚𝑥.(3.5)Now, only the 𝑐_𝑚 (for 𝑚=0 up to 𝑁−𝑛₊) need to be determined. These follow from the fact that 𝐹_𝑛(0)=0 for 𝑛=𝑛₊,…,𝑁. These are 𝑁−𝑛₋ equations in the same number of unknowns, and can be determined explicitly in terms of 𝑧₀,…,𝑧_𝑁−𝑛+, 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 𝑛₊; at the end of the busy period the number of flows present is at most 𝑛₋.
Denote by 𝑃^→ the distribution of the number of flows present at the end of the busy period. Let the matrices 𝑄^bDD, 𝑄^bDU, 𝑄^bUD, 𝑄^bUU be the submatrices that are obtained by partitioning 𝑄_b into down-states (i.e., states 𝑛 such that 𝑟_𝐴,𝑛<0) and up-states (𝑟_𝐴,𝑛>0); similarly, 𝐹^→(𝑥) is partitioned in 𝐹^→D(𝑥) and 𝐹^→U(𝑥).
Then, it is not hard to prove that𝑃^→=1𝐹^→D(0)𝑄^bDD,𝟏⋅^→𝐹^D(0)𝑄^bDD,(3.6)see [10, Equation (5.9)]; ⟨⋅,⋅⟩ denotes the inner product of two vectors. The mean idle period 𝔼𝐼 is given by𝔼𝐼=−∑^{𝑛−𝑛=0}𝐹_𝑛(0)𝐹^→D(0)𝑄^bDD,𝟏.(3.7)Finally, the mean busy period 𝔼𝐵 can be calculated. According to “renewal reward” ∑^𝑛−_𝑛=0𝐹_𝑛(0)=𝔼𝐼/(𝔼𝐼+𝔼𝐵), so that𝔼𝐵=𝔼𝐼⋅1−∑^{𝑛−𝑛=0}𝐹_𝑛(0)∑^𝑛−_𝑛=0𝐹_𝑛(0).(3.8)

3.2. Analysis of Buffer Content Distribution

Now, we turn back to Model I, as described in Section 2.1. Our goal is to show that the steady-state buffer content 𝑊^⋆ 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 𝑉^⋆ 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 Models I and II a busy period starts with 𝑛₊ flow present. Also, the distribution of the length of the busy period 𝐵 is the same for both models, as well as the distribution 𝑃^→ 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 𝔼𝐼 of Model II. Let us therefore consider the mean idle period 𝔼𝐽 of Model I.

As in [7, Lemma 2.3], we have that the mean idle period of Model I equals𝔼𝐽=−𝑃^→𝑄⁻¹_i,𝟏;(3.9)the expected amount of time during this idle time in which there are 𝑛 flows present, say 𝔼(𝐽_𝑛), is (−𝑃^→𝑄⁻¹_i)_𝑛, that is, the 𝑛th entry of −𝑃^→𝑄⁻¹_i. This follows from the fact that the mean time 𝔼_𝑚(𝐽_𝑛) spent in 𝑛 during the idle time, given that at the beginning of the idle time 𝑚 flows were present, satisfies the linear system, for 𝑚=1,…,𝑛₋,𝔼_𝑚(𝐽_𝑛)=1𝜆+𝜇𝐶/21𝑚=𝑛+𝜆+𝜆𝜇𝐶/2⋅𝔼^𝑚+1𝐽^𝑛+𝜆+𝜇𝐶𝜇/𝐶2/2⋅𝔼^𝑚−1𝐽^𝑛,(3.10)with 𝔼_𝑛+(𝐽_𝑛)=0, and𝔼₀𝐽^𝑛=1𝜆1𝑛=0+𝔼¹𝐽^𝑛.(3.11)We now have collected all the required elements to determine the distribution of the steady-state buffer content 𝑊^⋆. Analogously to [7, Theorem 2.4] we obtain the following result for 𝐺_𝑛(𝑥)∶=ℙ(𝑊^⋆≤𝑥,𝑁^⋆=𝑛).

Theorem 3.1. For all 𝑥≥0, 𝐺_𝑛(𝑥)=𝐹_𝑛(𝑥)−𝐹_𝑛(0)⋅𝔼𝐼+𝔼𝐵𝔼𝐽+𝔼𝐵+𝔼𝐽^𝑛𝔼𝐽+𝔼𝐵.(3.12)

Proof. This is proven as follows. We first condition on 𝑊^⋆ being positive or zero:𝐺_𝑛(𝑥)=ℙ𝑊^⋆≤𝑥,𝑁^⋆=𝑛∣𝑊^⋆>0ℙ𝑊^⋆>0+ℙ𝑊^⋆=0,𝑁^⋆=𝑛.(3.13)By applying the renewal-reward theorem, the latter probability can be rewritten asℙ𝑊^⋆=0,𝑁^⋆=𝑛=𝔼𝐽^𝑛𝔼𝐽+𝔼𝐵;(3.14)notice that these probabilities equal 0 for 𝑛=𝑛₊,…,𝑁. Also, from “renewal-reward,”ℙ𝑊^⋆>0ℙ𝑉^⋆>0=𝔼𝔼𝐽𝐼++𝔼𝔼𝐵𝐵.(3.15)Hence, we are left with determining the first probability in the right-hand side of (3.13). We first rewrite it asℙ𝑁^⋆=𝑛∣𝑊^⋆>0ℙ𝑊^⋆>0−ℙ𝑊^⋆>𝑥,𝑁^⋆=𝑛∣𝑊^⋆>0ℙ𝑊^⋆>0.(3.16)Now, recall that the distribution of (𝑊^⋆,𝑁^⋆), conditional on 𝑊^⋆>0, is the same as the distribution of (𝑉^⋆,𝑀^⋆), conditional on 𝑉^⋆>0. Combining this with (3.15), we obtainℙ𝑀=^⋆=ℙ𝑛𝑉∣𝑉^⋆>0ℙ𝑊^⋆>0−ℙ𝑉^⋆>𝑥,𝑀^⋆=𝑛∣𝑉^⋆>0ℙ𝑊^⋆>0_⋆>0,𝑀^⋆=𝑛−ℙ𝑉^⋆>𝑥,𝑀^⋆=𝑛×𝔼𝐼+𝔼𝐵=𝐹𝔼𝐽+𝔼𝐵_𝑛(𝑥)−𝐹_𝑛(0)⋅𝔼𝐼+𝔼𝐵𝔼𝐽+𝔼𝐵.(3.17)This proves the claim.

Upon combining the above theorem with representation (3.5), we find the following useful result.

Corollary 3.2. For all 𝑥≥0, Theorem 3.1, in conjunction with (3.5), defines numbers 𝛼_𝑛 and 𝛽_𝑛𝑚 (with 𝑛=0,…,𝑁 and 𝑚=0,…,𝑁−𝑛₊) such that 𝐺_𝑛(𝑥)=𝛼_𝑛+^𝑁−𝑛+_𝑚₌₀𝛽_𝑛𝑚𝑒^𝑧𝑚𝑥.(3.18)Here, 𝐺_𝑛(0)>0 if and only if 𝑛∈{0,…,𝑛₋}; 𝐺_𝑛(0) is given by 𝐺_𝑛(0)=ℙ𝑊^⋆=0,𝑁^⋆=𝑛=𝛼_𝑛+^𝑁−𝑛+_𝑚₌₀𝛽_𝑛𝑚.(3.19)The probability of 𝑛 flows in the system is given by ℙ(𝑁^⋆=𝑛)=𝐺_𝑛(∞)=𝛼_𝑛 (where ∑^𝑁𝑛=0𝛼_𝑛=1). The Laplace transform of 𝑊^⋆ reads, for 𝑠>0, 𝔼𝑒^{−𝑠𝑊⋆}1𝑁^⋆=𝑛=𝛼_𝑛+^𝑁−𝑛+_𝑚₌₀𝑠𝛽_𝑛𝑚𝑠−𝑧_𝑚.(3.20)

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 flows after 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

Let 𝐷^⋆ 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 𝑂(0,𝑡) denote the amount of output capacity available in the interval [0,𝑡). 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 𝑟_𝑂,𝑛 until the particle has been served (in fact even until the queue is empty). Define, for 𝑧≥0, the random variable 𝜏_𝑧 as the time until 𝑧 units of service have become available:𝜏_𝑧∶=inf𝑡≥0∶𝑂(0,𝑡)=𝑧=inf𝑡≥0∶^𝑡0𝑟_{𝑂,𝑁𝑠}d𝑠=𝑧;(4.1)notice that 𝑂(0,𝑡) is increasing in 𝑡. Then, analogously to [4, Section 4.1], with some abuse of notation,𝔼𝑒^{−𝑠𝐷⋆}=^𝑁_𝑛=0^∞0𝔼𝑒^{−𝑠𝜏𝑧}∣𝑁^⋆=𝑛ℙ𝑊^⋆=𝑧,𝑁^⋆=𝑛d𝑧.(4.2)Hence, to further compute this expression, we need to evaluate 𝔼(𝑒^{−𝑠𝜏𝑧}∣𝑁^⋆=𝑛). Here, we can use [4, Proposition 4.1]:𝔼𝑒^{−𝑠𝜏𝑧}∣𝑁^⋆=𝑛=exp𝑅^𝑂−1𝑄_b−𝑠𝑅^𝑂−1𝑧𝟏_𝑛,(4.3)where 𝟏 denotes an (𝑁+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 𝛾_𝑚𝑛, with 𝑚,𝑛=0,…,𝑁,𝔼𝑒^{−𝑠𝜏𝑧}∣𝑁^⋆=𝑛=^𝑁_𝑚₌₀𝛾_𝑚𝑛𝑒^{𝛿𝑚(𝑠)𝑧}.(4.4)We thus obtain the following result.

Theorem 4.1. For 𝑠>0, 𝔼𝑒^{−𝑠𝐷⋆}=^𝑁_𝑁𝑛=0_𝑚=0𝛾_𝑚𝑛𝔼𝑒^{𝛿𝑚(𝑠)𝑊⋆}1𝑁^⋆=𝑛,(4.5)where the 𝛾_𝑚𝑛 are as in (4.4). The 𝛿_𝑛(𝑠), for 𝑛=0,…,𝑁, are the eigenvalues of 𝑅^𝑂−1𝑄_b−𝑠𝑅^𝑂−1 (which are negative). An expression for 𝔼(𝑒^{−𝑠𝑊⋆}1{𝑁^⋆=𝑛}), with 𝑠>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 (𝑊₀=𝑥) that the fluid particle sees (say that it arrives at time 0), and the number of flows that are then present (𝑁₀=𝑛). Define, for given 𝑛=1,…,𝑁 and 𝑥≥0, the transform of the queueing delay:𝜉_𝑛(𝑠∣𝑥)∶=𝔼𝑒^−𝑠𝐷∣𝑊₀=𝑥,𝑁₀=𝑛.(4.6)Then, we also introduce the transform of 𝜉_𝑛(𝑠∣𝑥) with respect to the workload 𝑥:𝐾_𝑛(𝑠,𝑡)=^∞0𝑒^−𝑡𝑥𝜉_𝑛(𝑠∣𝑥)d𝑥;(4.7)we say that 𝐾_𝑛(𝑠,𝑡) is a “double transform.” Below, we show how to use these double transforms to derive 𝔼𝑒^{−𝑠𝐷⋆}.

Our first goal is to characterise the 𝐾_𝑛(𝑠,𝑡), 𝑛=1,…,𝑁, for fixed 𝑠≥0 and 𝑡>0. We do this by expressing 𝐾_𝑛(𝑠,𝑡) in terms of 𝐾_𝑚(𝑠,𝑡) (with 𝑚≠𝑛) as follows. Condition on the time until the service rate changes; this time has an exponential distribution with mean 𝜈^𝑛−1. Hence,𝐾_𝑛(𝑠,𝑡)=^∞0𝑒^−𝑡𝑥^{𝑥/𝑟𝑂,𝑛}₀𝜈_𝑛𝑒^{−𝜈𝑛𝑦}𝑒^−𝑠𝑦𝜆1𝑛<𝑛𝜈_𝑛𝜉_𝑛+1𝑠∣𝑥−𝑟^𝑂,𝑛𝑦+𝜇^𝑛𝐶𝜈_𝑛+𝜉_𝑛−1𝑠∣𝑥−𝑟^𝑂,𝑛𝑦d𝑦_{∞𝑥/𝑟𝑂,𝑛}𝜈_𝑛𝑒^{−𝜈𝑛𝑦}𝑒^{−𝑠𝑥/𝑟𝑂,𝑛}d𝑦d𝑥.(4.8)A straightforward change of variable (𝑥−𝑟_𝑂,𝑛𝑦=𝑧) then yields that𝐾_𝑛(𝑠,𝑡)=1𝜈_𝑛+𝑠+tr_𝑂,𝑛𝜆1𝑛<𝑛𝐾^𝑛+1(𝑠,𝑡)+𝜇_𝑛𝑐𝐾_𝑛−1(𝑠,𝑡)+𝑟_𝑂,𝑛.(4.9)For given 𝑠 and 𝑡, these are 𝑁+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 𝐾_𝑛(𝑠,𝑡).

Our second goal is to show how these 𝐾_𝑛(𝑠,𝑡) yield an expression for 𝔼𝑒^{−𝑠𝐷⋆}. At an arbitrary point in time, the distribution function of the workload (jointly with the number of flows present) is given by 𝐺_𝑛(⋅), as given by (3.18). But, as the corresponding density is the weighted sum of exponentials, it entails that knowledge of the 𝐾_𝑛(𝑠,𝑡) gives an expression for the Laplace transform of the virtual delay:𝔼𝑒^{−𝑠𝐷⋆}=^𝑁_𝑛=0_[0,∞)=𝜉_𝑛(𝑠∣𝑥)d𝐺_𝑛(𝑥)_𝑁_𝑛=0𝜉^𝑛(𝑠∣0)𝐺_𝑛(0)+_{(0,∞)𝑁−𝑛}⁻_𝑚₌₀=𝜉_𝑛(𝑠∣𝑥)𝑧_𝑚𝛽_𝑛𝑚𝑒^𝑧𝑚𝑥d𝑥_𝑛−_𝑛=0𝐺_𝑛(0)+^𝑁_{𝑛=0𝑁−𝑛}⁻_𝑚=0𝑧_𝑚𝛽_𝑛𝑚𝐾_𝑛𝑠,−𝑧^𝑚;(4.10)recall that 𝑧_𝑚<0 for 𝑚=0,…,𝑁−𝑛₋.

Theorem 4.2. For 𝑠>0, 𝔼𝑒^{−𝑠𝐷⋆}=^𝑛−_𝑛=0𝐺_𝑛(0)+^𝑁_{𝑛=0𝑁−𝑛}⁻_𝑚=0𝑧_𝑚𝛽_𝑛𝑚𝐾_𝑛𝑠,−𝑧^𝑚.(4.11)

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 delay 𝐷^⋆ “at an arbitrary point in time” and delay 𝐷^⋆ “seen by an arbitrary fluid particle.” The correction to be made is analogous to [4, Section 4.2] and rather straightforward:𝔼𝑒^{−𝑠𝐷⋆}=^𝑁_𝑛=0∑𝑟_𝑁𝑘=0𝜋^𝑖,𝑛_𝑘𝑟_𝐼,𝑘^𝑁_𝑚=0𝛾_𝑚𝑛𝔼(𝑒^{𝛿𝑚(𝑠)𝑊⋆}1𝑁^⋆=𝑛),(4.12)compared to [11, Proposition 7.2].

5. Flow Transfer Delay Distribution

In this section, we focus on the time 𝐹 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 𝑊=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 (𝑊₀=𝑥) and the number of flows (𝑁₀=𝑛). Define, for 𝑠≥0, 𝑡>0, the Laplace transform of the flow transfer time (conditional on 𝑊₀=𝑥 and 𝑁₀=𝑛) and its transform with respect to 𝑥:𝜂𝐿_𝑛(𝑠∣𝑥)∶=𝔼𝑒^−𝑠𝐹∣𝑊₀=𝑥,𝑁₀=𝑛_𝑛;(𝑠,𝑡)=^∞0𝑒^−𝑡𝑥𝜂_𝑛(𝑠∣𝑥)d𝑥.(5.1)For later reference, we also introduce, for 𝑛=1,…,𝑁 and 𝑚>𝑊,𝜅_𝑛(𝑠,𝑡)∶=𝜈_𝑛+𝑠−tr_𝐴,𝑛;𝑡_𝑛(𝑠)∶=𝜈^𝑛+𝑠𝑟_𝐴,𝑛;ℓ_𝑛,𝑚(𝑠)∶=𝐿_𝑛𝑠,𝑡^𝑚(𝑠).(5.2)Notice that, for 𝑛>𝑊, 𝑡_𝑛(𝑠) is positive.

In Section 5.1, we find the 𝐿_𝑛(𝑠,𝑡); 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 the 𝐿_𝑛(𝑠,𝑡), 𝑛=1,…,𝑁. We do so by distinguishing between “down-states” (𝑛 with 𝑟_𝐴,𝑛>0) and “up-states” (with 𝑟_𝐴,𝑛<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 that 𝑛 is an “up-state”: 𝑛∈{𝑛₊,…,𝑁}. It is elementary to see that, conditioning on the first event taking place after 𝑤 units of time,𝐿_𝑛(𝑠,𝑡)=^∞0𝑒^−𝑡𝑥^∞0𝜈_𝑛𝑒^{−𝜈𝑛𝑤}𝑒^−𝑠𝑤×𝜆1𝑛<𝑁𝜈_𝑛𝜂_𝑛+1𝑠∣𝑥+𝑟^𝐴,𝑛𝑤+𝑛−1𝑛𝜇^𝑛𝐶𝜈_𝑛𝜂_𝑛−1𝑠∣𝑥+𝑟^𝐴,𝑛𝑤+1𝑛𝜇^𝑛𝐶𝜈_𝑛d𝑤(d5𝑥.3.)This is the sum of three integrals. The third equals𝐿^𝑛(3)(𝑠,𝑡)=𝜇^𝑛𝐶𝑛⋅1𝑡⋅𝜈_𝑛1+𝑠=∶𝜇^+𝑛(𝑠,𝑡).(5.4) Consider the first, and perform the change of variable 𝑥+𝑟_𝐴,𝑛𝑤=𝑦:𝐿^𝑛(1)(𝑠,𝑡)==^∞0𝑒^−𝑡𝑥^∞0𝜆1𝑛<𝑛𝑒^{−𝜈𝑛𝑤}𝑒^−𝑠𝑤𝜂_𝑛+1𝑠∣𝑥+𝑟^𝐴,𝑛𝑤d𝑤d𝑥_∞0^{𝑦/𝑟𝐴,𝑛}₀=𝑒^{−𝑡(𝑦−𝑟𝐴,𝑛𝑤)}𝜆1𝑛<𝑛𝑒^{−𝜈𝑛𝑤}𝑒^−𝑠𝑤𝜂_𝑛+1(𝑠∣𝑦)d𝑤d𝑦_∞0𝑒^−𝑡𝑦𝜆1𝑛<𝑛𝑒^{(𝑡−𝜈𝑛/𝑟𝐴,𝑛−𝑠/𝑟𝐴,𝑛)𝑦}−1tr_𝐴,𝑛−𝜈_𝑛−𝑠𝜂^𝑛+1=𝜆_𝑛(𝑠∣𝑦)d𝑦(𝑠,𝑡)𝐿_𝑛+1(𝑠,𝑡)−𝐿_𝑛+1𝑠,𝜈^𝑛+𝑠𝑟_𝐴,𝑛,where𝜆^𝑛(𝑠,𝑡)∶=𝜆1𝑛<𝑁𝜅_𝑛(𝑠,𝑡).(5.5)Similarly,𝐿^𝑛(2)(𝑠,𝑡)=𝜇_𝑛(𝑠,𝑡)𝐿_𝑛−1(𝑠,𝑡)−𝐿_𝑛−1𝑠,𝜈^𝑛+𝑠𝑟_𝐴,𝑛,where𝜇^𝑛(𝑠,𝑡)∶=𝑛−1𝑛𝜇^𝑛𝐶𝜅_𝑛(𝑠,𝑡).(5.6)We arrive at𝐿_𝑛(𝑠,𝑡)=𝜆_𝑛(𝑠,𝑡)𝐿_𝑛+1𝑠,𝑡−ℓ_𝑛+1,𝑛(𝑠)+𝜇_𝑛(𝑠,𝑡)𝐿_𝑛−1(𝑠,𝑡)−ℓ_{𝑛−1,𝑛}(𝑠)+𝜇^+𝑛(𝑠,𝑡).(5.7)Later, it will turn out to be also useful to consider the representation𝜅_𝑛(𝑠,𝑡)𝐿_𝑛(𝑠,𝑡)=𝜆1𝑛<𝑛⋅𝐿_𝑛+1(𝑠,𝑡)+𝑛−1𝑛𝜇^𝑛𝑐⋅𝐿_𝑛−1(𝑠,𝑡)+𝜇^𝑛𝐶𝑛⋅𝜅^𝑛(𝑠,𝑡)−𝜆1𝑛<𝑛⋅ℓ_𝑛+1,𝑛𝑡⋅𝜈_𝑛1+𝑠(𝑠)−𝑛−1𝑛𝜇^𝑛𝑐⋅ℓ_{𝑛−1,𝑛}(𝑠).(5.8)

Down-State
Now, assume that 𝑛 is a down-state: 𝑛∈{0,…,𝑛₋}. In this case, we must distinguish between the cases that the process remains in state 𝑛 shorter, respectively, longer than −𝑥/𝑟_𝐴,𝑛 (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𝐿_𝑛(𝑠,𝑡)=^∞0𝑒^−𝑡𝑥^{−𝑥/𝑟𝐴,𝑛}₀𝜈_𝑛𝑒^{−𝜈𝑛𝑤}𝑒^−𝑠𝑤×𝜆𝜈_𝑛𝜂_𝑛+1𝑠∣𝑥+𝑟^𝐴,𝑛𝑤+𝑛−1𝑛𝜇^𝑛𝐶𝜈_𝑛𝜂_𝑛−1𝑠∣𝑥+𝑟^𝐴,𝑛𝑤+1𝑛𝜇^𝑛𝐶+𝜈_𝑛d𝑤_{∞−𝑥/𝑟𝐴,𝑛}𝜈_𝑛𝑒^{−𝜈𝑛𝑥}𝑒^{𝑠𝑥/𝑟𝐴,𝑛}𝜂_𝑛(𝑠∣0)d𝑤d𝑥.(5.9)With 𝜇^−𝑛(𝑠,𝑡)∶=𝜇_𝑛𝐶/(𝑛⋅𝜅_𝑛(𝑠,𝑡)⋅𝑡), this simplifies to𝐿_𝑛(𝑠,𝑡)=𝜆_𝑛(𝑠,𝑡)𝐿_𝑛+1(𝑠,𝑡)+𝜇_𝑛(𝑠,𝑡)𝐿_𝑛−1(𝑠,𝑡)+𝜇^−𝑛(𝑠,𝑡)−𝑟^𝐴,𝑛𝜅_𝑛(𝑠,𝑡)⋅𝜂^𝑛(𝑠∣0).(5.10)As indicated, our goal is to generate a system of equations for the 𝐿_𝑛(𝑠,𝑡) (with 𝑠 and 𝑡 fixed); we therefore wish to express 𝜂_𝑛(𝑠∣0) in terms of the 𝐿_𝑛(𝑠,𝑡). This can be done as follows.
(i) First, for 𝑛∈{0,…,𝑛₋},𝜂_𝑛(𝑠∣0)=1𝜈_𝑛+𝑠𝜆𝜂^𝑛+1(𝑠∣0)+𝑛−1𝑛𝜇^𝑛𝑐⋅𝜂_𝑛−1(𝑠∣0)+𝜇^𝑛𝐶𝑛,(5.11)whereas for 𝑛∈{𝑛₊,…,𝑁},𝜂_𝑛(𝑠∣0)=𝜆1𝑛<𝑁𝑟_𝐴,𝑛⋅𝐿_𝑛+1𝑠,𝜈^𝑛+𝑠+𝑛−1𝑟_𝐴,𝑛𝑛⋅𝜇^𝑛𝐶𝑟_𝐴,𝑛⋅𝐿_𝑛−1𝑠,𝜈^𝑛+𝑠𝑟_𝐴,𝑛+𝜇^𝑛𝐶=𝜆1𝑛<𝑁𝑛⋅𝜈_𝑛1+𝑠𝑟_𝐴,𝑛⋅ℓ_𝑛+1,𝑛(𝑠)+𝑛−1𝑛⋅𝜇^𝑛𝐶𝑟_𝐴,𝑛⋅ℓ_{𝑛−1,𝑛}(𝑠)+𝜇^𝑛𝐶𝑛⋅𝜈_𝑛1+𝑠.(5.12)(ii) Now, consider the vector ^→𝜂(𝑠)≡(𝜂₁(𝑠∣0),…,𝜂_𝑛−(𝑠∣0))^T. Define for 𝑛,𝑚=1,…,𝑛₋,𝑎_𝑛,𝑚(𝑠)∶=⎧⎪⎪⎨⎪⎪⎩−𝜈−(𝜆𝑛+,−𝑠,1)/𝑛⋅𝜇_𝑛iiff𝑚𝑚==𝑛𝑛;+1;_𝑛𝑐,if𝑚=𝑛+1.(5.13)and 0 else. The corresponding matrix is called 𝐴(𝑠), that is, 𝐴(𝑠)≡(𝑎_𝑛,𝑚(𝑠))^𝑛−_𝑛,𝑚=1; for 𝑠>0 we have that 𝐴(𝑠) is diagonally dominant and hence invertible. Also, ^→𝑢(𝑠)≡(𝑢₁(𝑠),…,𝑢_𝑛−(𝑠))^T, with𝑢_𝑛(𝑠)∶=⎧⎪⎨⎪⎩𝜇^𝑛𝜇𝑐/𝑛,if𝑛=1,…,𝑛₋−1,_𝑛−𝑐/𝑛₋+𝜂_𝑛+(𝑠∣0),if𝑛=𝑛₋.(5.14)Then, (5.11) implies that ^→𝜂(𝑠)=(𝐴(𝑠))^→−1𝑢(𝑠). In other words, once we know 𝜂_𝑛+(𝑠∣0), we can compute 𝜂₁(𝑠∣0),…,𝜂_𝑛−(𝑠∣0). (iii) Let 𝑎^{−1𝑛,𝑚}(𝑠)∶=(𝐴(𝑠))^{−1𝑛,𝑚}. Now, (5.12) entails that, for 𝑛=1,…,𝑛₋,𝜂_𝑛(𝑠∣0)=^𝑛−−1_𝑚₌₁𝑎^{−1𝑛,𝑚}(𝑠)𝜇^𝑚𝐶𝑚+𝑎^{−1𝑛,𝑛−}(𝑠)𝜇^𝑛−𝐶𝑛₋=𝛼_𝑛+𝜂_𝑛+(𝑠∣0)(𝑠)+𝛽_𝑛(𝑠)ℓ_{𝑛++1,𝑛+}(𝑠)+𝛾_𝑛(𝑠)ℓ_{𝑛−,𝑛+}(𝑠),(5.15) where𝛼_𝑛(𝑠)∶=^𝑛−_𝑚₌₁𝑎^{−1𝑛,𝑚}(𝑠)𝜇^𝑚𝐶𝑚+𝑎^{−1𝑛,𝑛−}(𝑠)𝜇^𝑛−𝐶𝑛₋⋅1𝛽𝜈_𝑛−+𝑠_𝑛;(𝑠)∶=𝑎^{−1𝑛,𝑛}₋(𝑠)⋅𝜆𝑟_𝐴,𝑛+;𝛾_𝑛(𝑠)∶=𝑎^{−1𝑛,𝑛}₋(𝑠)⋅𝑛⁻𝑛₊⋅𝜇^𝑛+𝐶𝑟_𝐴,𝑛+.(5.16) Inserting this into (5.10), we have found the following relation for down-states 𝑛:𝜅_𝑛(𝑠,𝑡)𝐿_𝑛(𝑠,𝑡)=𝜆⋅𝐿_𝑛+1(𝑠,𝑡)+𝑛−1𝑛𝜇^𝑛+𝜇𝑐⋅𝐿_𝑛−1(𝑠,𝑡)_𝑛𝐶𝑛𝑡−𝑟^𝐴,𝑛⋅𝛼_𝑛(𝑠)+𝛽_𝑛(𝑠)ℓ_{𝑛++1,𝑛+}(𝑠)+𝛾_𝑛(𝑠)ℓ_{𝑛−,𝑛+}(𝑠);(5.17)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 𝜙_𝑛(⋅,⋅) and 𝜓_𝑚,𝑘(⋅),𝐿^→_𝑛(𝑠,𝑡)=^𝑁_𝑚=1𝑄(𝑠,𝑡⁻¹_𝑛,𝑚(𝜙_𝑚(𝑠,𝑡)+^𝑁_𝑘=𝑛+𝜓_𝑚,𝑘(𝑠)ℓ_𝑚,𝑘(𝑠));(5.18)here, the 𝑁×𝑁 matrix 𝑄(𝑠,𝑡) is given through𝑄(𝑠,𝑡)_𝑛,𝑚∶=⎧⎪⎪⎨⎪⎪⎩𝜆−(𝑛1𝜅−_𝑛𝑛(1𝑠<,)/𝑡𝑛)𝑛,⋅𝜇,_𝑛iiff𝑚𝑚==𝑛𝑛;+1;𝑐,if𝑚=𝑛−1.(5.19)In other words, if the transforms ℓ_𝑛,𝑚(𝑠), for 𝑛=1,…,𝑁 and 𝑚=𝑛₊,…,𝑁, would be known, then, for fixed 𝑠 and 𝑡, the values of the 𝐿_𝑛(𝑠,𝑡) can be found directly from solving a system of linear equations. The rest of this subsection is devoted to explaining how to identify the ℓ_𝑛,𝑚(𝑠), for given 𝑠>0. We first prove a useful lemma.

Lemma 5.1. Consider, for fixed 𝑠>0, the 𝑡∈ℂ such that de𝑡(𝑄(𝑠,𝑡))=0. There are 𝑁−𝑛₋ such values such that Re𝑡>0.

Proof. First, rewrite 𝑄=𝑄^⌣−𝐷+𝑡𝑅_𝐴, with𝑄̆_𝑛,𝑚∶=⎧⎪⎪⎨⎪⎪⎩𝜆−(𝑛1𝜆−1𝑛1𝑛≤)/𝑛≤𝑛⋅𝑛𝜇_𝑛−(𝑛−1)/𝑛⋅𝜇_𝑛𝑐iiff𝑚𝑚==𝑛𝑛,+1,𝑐if𝑚=𝑛−1,(5.20)𝑑_𝑛∶=𝑛⁻¹⋅𝜇_𝑛𝐶+𝑠>0, 𝐷∶=diag{𝑑⃗}. Observe that 𝑄^⌣ is a generator matrix. Notice also that solutions of det(𝑄(𝑠,𝑡))=0 are the eigenvalues of −𝑅^𝐴−1(𝑄^⌣−𝐷).
(i) We first focus on properties of the eigenvalues of −𝑅^𝐴⌣−1𝑄. Recall that it follows from Theorem 2, part 3 of [8] that −𝑅^𝐴⌣−1𝑄 has as many eigenvalues in the right half plane as the number of up-states in 𝑅_𝐴, that is, 𝑀∶=𝑁−𝑛₋; there is also one eigenvalue of −𝑅^𝐴⌣−1𝑄 equal to 0 (note that 𝑄^⌣ is singular), and the remaining 𝑛₋−1 are in the left half plane. “Ger𝑧∈ℂ∶|||𝑧+^⌣𝑄_𝑛,𝑛𝑟_𝐴,𝑛|||≤|||𝑄^⌣_𝑛,𝑛𝑟_𝐴,𝑛|||;(5.21)gorin” states that all eigenvalues are in at least one of the disks𝑛the −𝑄^⌣_𝑛,𝑛/𝑟_𝐴,𝑛th disk is a circle in the complex plane around |𝑄^⌣_𝑛,𝑛/𝑟_𝐴,𝑛| of radius 𝑄^⌣_𝑛,𝑛<0 (which therefore goes through 0). Notice that −𝑅^𝐴−1(𝑄^⌣−𝜀𝐷) 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 𝜀>0. for small 𝜁(𝜀,𝑡)∶=det(𝑄^⌣−𝜀𝐷+𝑡𝑅_𝐴)=0. Observe that these solve the equation 𝜕𝜕𝑡𝜁(𝜀,𝑡)|||_𝜀,𝑡=0(5.22) As seen in the proof of [8, Theorem 2],𝜁has 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 −𝑅^𝐴⌣−1𝑄 are positive). Hence, replacing −𝑅^𝐴−1(𝑄^⌣−𝜀𝐷) by −𝑅^𝐴−1(𝑄^⌣−𝜀𝐷) moves the zero eigenvalue to the left.
“Ger𝑧∈ℂ∶|||𝑧+^⌣𝑄_𝑛,𝑛−𝜀𝑑_𝑛𝑟_𝐴,𝑛|||≤|||𝑄^⌣_𝑛,𝑛𝑟_𝐴,𝑛|||;(5.23)gorin” implies that all eigenvalues of 𝜀=0 are in at least one of the disks𝜁(𝜀,𝑡)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 𝑁−𝑛₋ in the coefficients of the characteristic polynomial.
Conclude that for small 𝑛₋, there are then 𝜀 in the right half plane, and the remaining 𝜀=1 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 −𝑅^𝐴⌣−1𝑄 further. The special case −𝑅^𝐴−1(𝑄^⌣−𝜀𝐷) proves our claim.

(a)

(b)

Figure 2 

Eigenvalues in the complex plane. (a) Eigenvalues of 𝑁=5; (b) eigenvalues of ℓ𝑛,𝑚(𝑠). In this example →ℓ(𝑠), and there are 3 up-states and 2 down-states.

Now, we are able to characterise the transforms 𝑛=𝑛₊,…,𝑁 as follows.

STEP 1. Determine linear equations for the entries of 𝜅_𝑛,𝑚(𝑠)∶=𝜅_𝑛(𝑠,𝑡_𝑚(𝑠)).
(a)First, focus on (5.8); these relate to 𝑛≠𝑚. We introduce the notation 𝑚=𝑛₊,…,𝑁. Then, for 𝑡=𝑡_𝑚(𝑠) and 𝜅_𝑛,𝑚(𝑠)ℓ_𝑛,𝑚(𝑠)=+𝜆1𝑛𝑛−≤1𝑛ℓ_𝑛+1,𝑚(𝑠)−ℓ_𝑛+1,𝑛(𝑠)𝑛𝜇^𝑛𝑐ℓ_{𝑛−1,𝑚}(𝑠)−ℓ_{𝑛−1,𝑛}(𝑠)+𝜅_𝑛,𝑚(𝑠)𝜇^+𝑛𝑠,𝑡^𝑚(𝑠).(5.24), we obtain by inserting 𝜅_𝑛,𝑛(𝑠)=𝜅_𝑛(𝑠,𝑡_𝑛(𝑠))=0:𝑛=𝑚Notice that 𝑀, so that for 𝑛=1,…,𝑛₋ both sides reduce to 0; hence, these 𝑡=𝑡_𝑚(𝑠) equations are meaningless.(b)Now, consider (5.17); these relate to 𝑚=𝑛₊,…,𝑁. Plugging in 𝜅_𝑛,𝑚(𝑠)ℓ_𝑛,𝑚(𝑠)=𝜆ℓ_𝑛+1,𝑚(𝑠)+𝑛−1𝑛𝜇^𝑛𝑐ℓ_{𝑛−1,𝑚}(𝑠)+𝜇^𝑛𝐶−𝑟_𝐴,𝑛𝑛𝑡_𝑚(𝑠)⋅𝛼_𝑛(𝑠)+𝛽_𝑛(𝑠)ℓ_{𝑛++1,𝑛+}(𝑠)+𝛾_𝑛(𝑠)ℓ_{𝑛−,𝑛+}(𝑠).(5.25) for ^→ℓ(𝑠), we obtainℓ_𝑛,𝑚(𝑠)

STEP 2. Reduce dimension of the vector 𝑛=1,…,𝑁. The sets of (5.24) and (5.25) enable us to express the 𝑚=𝑛₊,…,𝑁 for 𝑛≠𝑚 and ^→ℓ₊(𝑠)∶=(ℓ_𝑛+,𝑛+(𝑠),…,ℓ_𝑁,𝑁(𝑠))^T, but 𝜑_𝑚(⋅,⋅), in terms of 𝜗_𝑚,𝑘(⋅). We have thus identified functions 𝑛=1,…,𝑁 and 𝐿_𝑛(𝑠,𝑡)=^𝑁_𝑚₌₁𝑄(𝑠,𝑡)⁻¹_𝑛,𝑚(𝜑_𝑚(𝑠,𝑡)+^𝑁_𝑘_=𝑛+𝜗_𝑚,𝑘(𝑠)ℓ_𝑘,𝑘(𝑠)).(5.26) such that for 𝑁,ℓ_𝑛+,𝑛+(𝑠),…,ℓ_𝑁,𝑁(𝑠)In other words, when the 𝐿_𝑛(𝑠,𝑡) functions 𝐿_𝑛(𝑠,𝑡)=det𝑄_𝑛(𝑠,𝑡)det𝑄(𝑠,𝑡),(5.27) would be known, we would have found the 𝑄_𝑛(𝑠,𝑡).

STEP 3. Apply Lemma 5.1. By virtue of Cramer's rule, we obtain from (5.26) that𝑄(𝑠,𝑡)where 𝑛 is defined as 𝑚 but with the 𝜑_𝑚(𝑠,𝑡)+∑^{𝑁𝑘=𝑛}₊𝜗_𝑚,𝑘(𝑠)ℓ_𝑘,𝑘(𝑠)th column replaced by a vector of which the 𝑡th entry is 𝑠>0. For any 𝜏₁(𝑠),…,𝜏_{𝑁−𝑛−}(𝑠) in the right half plane, this should have a finite norm. Now fix 𝑁−𝑛₋, and use Lemma 5.1. Denote the zeroes of the denominator by ℓ_𝑛+,𝑛+(𝑠),…,ℓ_𝑁,𝑁(𝑠). Conclude that each zero of the denominator should correspond to a zero of the numerator. This yields 𝑠,𝑡>0 equations that determine 𝐿^→(𝑠,𝑡)=𝑄(𝑠,𝑡)⁻¹𝜑^→(𝑠,𝑡)+Θ(𝑠)^→ℓ₊(𝑠),(5.28).

The above results are summarised in the following theorem.

Theorem 5.2. For ^→𝜙(𝑠), Θ(𝑠)with 𝑠>0 and 𝑁−𝑛₋ defined by (5.26). For any 𝑡 there are det𝑄(𝑠,𝑡)=0 values of 𝜏₁(𝑠),…,𝜏_{𝑁−𝑛−}(𝑠) in the right half plane such that ^→ℓ₊(𝑠), say 𝑠>0. The vector 𝑡→𝜏_𝑚(𝑠) follows, for fixed 0, from letting 𝑚=1,…,𝑁−𝑛₋. in the numerator of (5.27) and equating this to 𝜂_𝑛(𝑠∣0), for 𝑛=0,…,𝑁

Remark 5.3. It is easily verified that, in passing, we have also found a procedure to compute the 𝑛=1,…,𝑁, for 𝐻_𝑛∶=𝐺^𝑛−1(∞)∑^{𝑁−1𝑚=0}𝐺_𝑚(∞).(5.29), compared to (5.11), (5.12), (5.15).

5.3. Flow Transfer Delay

It is clear that, due to PASTA, the number of customers present at (i.e., just after) arrival of an (accepted) flow has distribution (𝐻_𝑛(𝑥)∶=𝐺^𝑛−1(𝑥)∑^{𝑁−1𝑚=0}𝐺_𝑚(∞)=𝐴×(𝛼^𝑛−1+^{𝑁−𝑛−}_𝑚₌₀𝛽_{𝑛−1,𝑚}𝑒^𝑧𝑚𝑥),𝐴∶=(^𝑁−1_𝑚₌₀𝛼_𝑚)⁻¹;(5.30))∑^𝑁𝑛=1𝐻_𝑛(∞)=1.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𝔼𝑒^−𝑠𝐹=^𝑁_𝑛=1_[0,∞)𝜂_𝑛(𝑠∣𝑥)d𝐻_𝑛(𝑥).(5.31)observe that indeed 𝔼𝑒^{−𝑠𝐷⋆} Hence,𝑠>0Mimicking the derivation of 𝔼𝑒^−𝑠𝐹=^𝑛+_𝑛=1𝐴𝐺_𝑛−1(0)⋅𝜂_𝑛(𝑠∣0)+^𝑁_{𝑛=1𝑁−𝑛}⁻_𝑚=0𝐴𝑧^𝑚𝛽_{𝑛−1,𝑚}⋅𝐿^𝑛𝑠,−𝑧^𝑚.(5.32), we obtain the following result.

Theorem 5.4. For 𝐿_𝑛(⋅,⋅), 𝜂_𝑛(⋅∣0)Expressions for 𝑆 and 𝐻_𝑛(⋅) are available from the previous subsection.

6. Sojourn Time Distribution

In this section, we study the sojourn time 𝑆 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 𝐹, as defined through (5.30). To derive the Laplace transform of 𝑆, we first need to describe the workload increment (which can be positive or negative) during the flow transfer time 𝐹, see Section 6.1. Then, in Section 6.2 we put the components together, and derive the transform of 𝑊_𝐹+.

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 of 𝑁_𝐹+, but rather its joint distribution with the workload 𝜂_𝑛(𝑠∣𝑤) and the number of flows present 𝐿_𝑛(𝑠,𝑡) at the end of the transfer (not counting the flow that just left). To this end we introduce the counterparts of 𝜂𝐿_𝑛,𝑚(^→𝑠∣𝑥)∶=𝔼𝑒^{−𝑠1𝐹−𝑠2𝑊𝐹+}1𝑁_𝐹+=𝑚∣𝑊₀=𝑥,𝑁₀=𝑛_𝑛,𝑚;(^→𝑠,𝑡)=^∞0𝑒^−𝑡𝑥𝜂_𝑛,𝑚(^→𝑠∣𝑥)d𝑥,(6.1) and ^→𝑠≡(𝑠₁,𝑠₂).:ℓ_{𝑛,𝑚,𝑘}(^→𝑠)∶=𝐿_𝑛,𝑚(^→𝑠,𝑡_𝑘(𝑠₁))with 𝐿_𝑛,𝑚(^→𝑠,𝑡) Similarly to before, we also define 𝐿_𝑛(𝑠,𝑡).

𝑛>𝑊 can be derived essentially in the same fashion as 𝐿_𝑛,𝑚(^→𝑠,𝑡)=+𝜆_𝑛𝜇𝑠¹_𝑛,𝑡𝐿_𝑛+1,𝑚(^→𝑠,𝑡)−ℓ_{𝑛+1,𝑚,𝑛}(^→𝑠)+𝜇_𝑛𝑠¹,𝑡𝐿_{𝑛−1,𝑚}(^→𝑠,𝑡)−ℓ_{𝑛−1,𝑚,𝑛}(^→𝑠)𝐶𝑛⋅𝑡+1𝑠₁⋅1𝜈_𝑛+𝑠₁+𝑠₂𝑟_𝐴,𝑛⋅1𝑛−1=𝑚,(6.2) was found in Sections 5.1 and 5.2. We sketch this procedure. The counterpart of (5.7) is, for 𝑛<𝑊,𝐿_𝑛,𝑚(^→𝑠,𝑡)=+𝜆_𝑛𝜇𝑠¹,𝑡𝐿_𝑛+1,𝑚(^→𝑠,𝑡)+𝜇_𝑛𝑠¹,𝑡𝐿_{𝑛−1,𝑚}(^→𝑠,𝑡)_𝑛𝐶𝑛⋅𝜅_𝑛1𝑠¹,𝑡⋅1𝑡+𝑠₂⋅1𝑛−1=𝑚−𝑟^𝐴,𝑛𝜅_𝑛𝑠¹,𝑡⋅𝜂^𝑛,𝑚(^→𝑠∣0).(6.3)whereas (5.10) generalises to, for 𝜂_𝑛,𝑚(^→𝑠∣0),ℓ_{𝑛,𝑚,𝑘}(^→𝑠)In the last equation, as before, the ℓ_{𝑛,𝑚,𝑘}(^→𝑠) can be expressed in terms of the 𝔼𝑒^−𝑠𝑆==𝔼_∞0𝑒^{−𝑠𝐹−𝑠𝜏𝑊𝐹+}_𝑁[0,∞)_{𝑛=1𝑁−1𝑚=0}×𝔼𝑒_{−𝑠𝜏𝑦}𝔼𝑒^−𝑠𝐹1𝑊_𝐹+=𝑦,𝑁_𝐹+=𝑚∣𝑊₀=𝑥,𝑁₀=𝑛∣𝑁₀=𝑚d𝐻_𝑛(𝑥)d𝑦.(6.4). Then, it remains to find the transforms {𝑊₀=0}; 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{𝑊₀>0}Consider the cases {𝑊₀=0} and ^∞0𝑁_{𝑛=1𝑁−1𝑚=0}𝔼𝑒^−𝑠𝐹1𝑊_𝐹+=𝑦,𝑁_𝐹+=𝑚∣𝑊₀=0,𝑁₀=𝑛^𝑁_𝑘=0=𝛾_𝑘𝑚𝑒^{𝛿𝑘(𝑠)𝑦}𝐴𝐺_𝑛−1(0)d𝑦_𝑁_{𝑛=1𝑁−1𝑁𝑚=0}_𝑘=0=𝛾_𝑘𝑚𝔼𝑒^{−𝑠𝐹+𝛿𝑘(𝑠)𝑊𝐹+}1𝑁_𝐹+=𝑚∣𝑊₀=0,𝑁₀=𝑛𝐴𝐺_𝑛−1(0)_𝑁_{𝑛=1𝑁−1𝑁𝑚=0}_𝑘=0𝛾_𝑘𝑚⋅𝐴𝐺_𝑛−1(0)⋅𝜂_𝑛,𝑚𝑠,−𝛿^𝑘(𝑠)∣0.(6.5) separately. The contribution due to {𝑊₀>0} amounts to_𝑁(0,∞)_{𝑛=1𝑁−1𝑁𝑚=0}_𝑘=0𝛾_𝑘𝑚𝜂_𝑛,𝑚𝑠,−𝛿^𝑘(𝑠)∣𝑥^{𝑁−𝑛−}_𝑗=0𝐴𝑧^𝑗=𝛽_{𝑛−1,𝑗}𝑒^𝑧𝑗𝑥d𝑥_𝑁_{𝑛=1𝑁−1𝑁𝑚=0}_{𝑘=0𝑁−𝑛}⁻_𝑗=0𝛾_𝑘𝑚⋅𝐴𝑧_𝑗𝛽_{𝑛−1,𝑗}⋅𝐿^𝑛,𝑚𝑠,−𝛿^𝑘(𝑠),−𝑧_𝑗.(6.6)Similarly, we find that the contribution due to 𝑠>0 is𝔼𝑒^−𝑠𝑆

Theorem 6.1. For 𝑊, 𝑊∈ℕ 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 weight 𝑟_𝐴,𝑊=0 was noninteger. If 0,…,𝑊, 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 𝑊); if we enumerate the equations 0=^𝑁_𝑛=0𝑞_𝑛,𝑊𝐹_𝑛(𝑥),or𝐹_𝑊(𝑥)=−1𝑞_𝑊,𝑊_𝑛≠𝑊𝑞_𝑛,𝑊𝐹_𝑛(𝑥),(7.1), then the 𝑞_𝑛,𝑚th equation reads(𝑛,𝑚)where 𝑄_b is the 𝑛-entry of 𝑛≠𝑊. Then the 𝑟_𝐴,𝑛𝐹^𝑛(𝑥)=_𝑚≠𝑊𝑞^{−𝑤𝑚,𝑛}𝐹_𝑚(𝑥),where𝑞^{−𝑤𝑚,𝑛}∶=𝑞_𝑚,𝑛−𝑞^𝑊,𝑛𝑞_𝑊,𝑊⋅𝑞_𝑚,𝑊.(7.2)th differential equation (where 𝑞^{−𝑊𝑚,𝑛}≥0) becomes𝑚≠𝑛Interestingly, ∑_𝑛≠𝑊𝑞^{−𝑊𝑚,𝑛}=0 for 𝑄^−𝑊_b∶=(𝑞^{−𝑊𝑚,𝑛})_𝑚,𝑛, and 𝑚,𝑛=0,…,𝑁, and hence the 𝑚,𝑛≠𝑊 (with 𝑁, but 𝑅^𝑎−𝑤𝐹^→−𝑊^(𝑥)=𝑄^−𝑤_b𝐹^T→−𝑊(𝑥),(7.3)) correspond to an 𝑅^𝐴−𝑊-dimensional generator matrix. In self-evident notation, we have arrived at𝐹_𝑛(⋅)with all entries of 𝑛≠𝑊 not equal to 0. In this way, the steady-state buffer content distribution of Model II can be determined: we first find the 𝐹_𝑊(⋅) for 𝑊, and then we use (7.1) to derive 𝑊.
The buffer content distribution follows in the same fashion as in Section 3.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 𝑁=∞. 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 weight 𝑁^⋆(𝑊). For ease, we lift the assumption of performing admission control; in other words: we take 𝑊 Let 𝑊^⋆(𝑊) be the steady-state number of flows in the system (i.e., transmitting traffic into the queue), for a given weight 2𝑊^⋆𝑠(𝑊)+𝑊^⋆(𝑊); 𝑊^⋆𝑠(𝑊) 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 𝑊^⋆(𝑊), 𝑊 denoting the traffic at the sources, and 4𝜆⁻¹⋅𝜚²/(1−2𝜚) 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 𝑊.; 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 𝑊, independently of the choice of 𝑟_𝐼,𝑛 It can be seen that pathwise the amount of traffic that is at the sources increases in 𝑊 (as the 𝑊=∞ decrease), so that the amount of traffic in the queue decreases in 𝑛.

(i) In case (1−2𝜚)(2𝜚)^𝑛, 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 𝔼𝑁^⋆(∞)=2𝜚/(1−2𝜚). flows in the system is 𝑊=0.; the mean number of flows is 𝐶.

(ii) A second extreme case is 𝑛 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 (1−𝜚)𝜚^𝑛 In this case, the probability of 𝔼𝑁^⋆(0)=𝜚/(1−𝜚). flows in the system is 𝑊.; the mean number of flows is 2𝔼𝑁^⋆(𝑤)⋅1𝜇,(7.4)

Evidently, the mean amount of work in the queue does depend on 𝔼𝑊^⋆(∞)=𝑐1𝔼𝑊_⋆𝜆⋅14−𝜚2²𝜚−𝜇2𝐶⋅1−2𝜚2𝜚=0;(0)=𝑐1𝜆⋅14−𝜚2²𝜚−𝜇2𝐶⋅1−𝜚𝜚=𝜇1⋅(1−𝜚)2(𝜚1−2𝜚).(7.5) First, observe that the mean amount of traffic (measured in processing time) that the flows still need to inject into the queue is𝑊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𝑊These arguments can be used to quantify the tradeoff between the flow transfer delay (which increases in 𝑊=1) and the queueing delay of the last particle (which decreases in 𝔼𝑊^⋆(1)=1𝜇⋅(1−𝜚4)(𝜚1²−2𝜚),(7.6)). The formula for 𝜚<1/2 was already given in [3]𝑊^⋆which is indeed between (7.5) for 𝐷^⋆.

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 of 𝐹 and 𝑆 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 𝑊=1 and lim_𝑥→∞1𝑥logℙ𝑊^⋆>𝑥=−𝜗^⋆,(7.7) 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 𝜗^⋆>0, in the spirit of𝐷^⋆,𝐹where the constant 𝑆 follows from the solution of the corresponding eigensystem; likewise, one could consider the logarithmic tail asymptotics of ℙ(𝑊^⋆>𝑥)=𝜓(𝑥)𝑒^{−𝜗⋆𝑥}, and 𝜓(⋅). As this just provides us with decay rate (i.e., it says that 𝑥⁻¹⋅log𝜓(𝑥)→0 for an unknown function −𝛿 such that 1−𝛿), 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 𝑥, 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 lim_𝑥↓0𝔼𝑆(𝑥)>0., say 𝑊, is such that 𝑊 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 𝑊, one may wonder what value of 𝑀𝑥=𝜆(Γ−𝑐𝐼)𝑥 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 𝑀 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].

Acknowledgments

The authors would like to thank Hans van den Berg and Frank Roijers (TNO ICT) for useful remarks.