Queue-length balance equations in multiclass multiserver queues and their generalizations

A classical result for the steady-state queue-length distribution of single-class queueing systems is the following: the distribution of the queue length just before an arrival epoch equals the distribution of the queue length just after a departure epoch. The constraint for this result to be valid is that arrivals, and also service completions, with probability one occur individually, i.e., not in batches. We show that it is easy to write down somewhat similar balance equations for {\em multidimensional} queue-length processes for a quite general network of multiclass multiserver queues. We formally derive those balance equations under a general framework. They are called distributional relationships, and are obtained for any external arrival process and state dependent routing as long as certain stationarity conditions are satisfied and external arrivals and service completions do not simultaneously occur. We demonstrate the use of these balance equations, in combination with PASTA, by (i) providing very simple derivations of some known results for polling systems, and (ii) obtaining new results for some queueing systems with priorities. We also extend the distributional relationships for a non-stationary framework.


Introduction
A classical result for the steady-state queue-length distribution of single-class queueing systems is the following: the distribution of the queue length just before an arrival epoch equals the distribution of the queue length just after a departure epoch. The constraint for this result to be valid is that, with probability one, arrivals, and also service completions, occur individually, i.e., not in batches. The result then follows by a simple level-crossing argument: in steady state, the event that a customer arrives to find j customers present occurs just as often as the event that a customer leaves j customers behind, for all j = 0, 1, . . . . See [7], pp. 154-156, for a formal statement and proof (due to P.J. Burke, unpublished) of this result.
At first sight this level-crossing argument breaks down in higher dimensions, for example in the case of multiple customer classes. Indeed, with x ≥ 0 and e k being a unit vector with 1 in the kth coordinate and zero elsewhere, an m-dimensional process can leave state x because of an arrival of a customer of type i, and enter that state from state x + e k because of a departure of a customer of another type k. However, we shall argue that it is easy to write down a more global balance equation for multidimensional queue length processes for a large class of queues and queueing networks -also when service times are not exponentially distributed, and even when arrivals may occur in batches. We shall explore that fact to obtain a simple relation between the steady-state joint queue-length distribution at arrival epochs (which under various circumstances is equal to the time average distribution) and at service completion epochs. Once one has a relation between the probability generating function (PGF) at arbitrary epochs and at service completion epochs, one can find the former when one has the latter. The latter results are indeed known in an M/G/1 setting, where it is natural to look at departure epochs. This will yield both new results (for multiclass queueing models with fixed priorities and for the longer-queue model), as well as new and simple derivations of known results for, e.g., polling models.
The research for the present paper was initially motivated by the desire to provide an intuitive explanation of a result in [3] regarding the steady-state joint queue-length distribution in a large class of polling models. That distribution turned out to have a remarkably simple relation with a weighted sum of the joint queue length distributions at departure epochs of customers from each of the queues. In Section 2 we provide such an explanation. Although balance equations are intuitively appealing, their mathematical verification may require a large amount of work. This motivates us to derive distributional relationships for queue lengths in a unified way using a general tool. The so called rate conservation law is such a tool as demonstrated in [15] (also see [1,14]). This method is applicable to a general model, but requires Palm distributions, which may not be easy to understand. In Section 3 of this paper we take another approach, based on a time evolution of a sample path. This approach is parallel to the rate conservation law, but does not require Palm distributions, which are replaced by sample averages. We apply it to a general model, and derive a distributional relationship among different embedded epochs. In Section 4 a non-stationary version of the distributional relationship is derived with some error term, which vanishes as time goes to infinity. Our main result, viz. Theorem 1, as well as the non-stationary results, are novel to the best of our knowledge.
Literature review. Hébuterne [11] provides a generalization of the above-mentioned classical result of Burke in two directions: he allows (i) batch arrivals, with batches of random size, and (ii) batch services, with batches of fixed size. He also points out that emptying the queue up to N customers is beyond the scope of the analysis, because then the batch sizes are not independent of the system state. Fakinos [9] manages to treat a quite general grouparrival group-departure queue. He treats the batch size problem by assuming that customers within a departing group are randomly ordered, and that they leave the system according to their order. Papaconstantinou and Bertsimas [16] generalize Burke's result to the multi-server E k /G/s queue. Kim [13] combines the features of batch arrivals, batch services and multiple servers, also allowing multiple customer classes. he does not explicitly address the issue of customers in a departing group being randomly ordered. Hébuterne and Rosenberg [12] focus on the G/G/1 queue with batch services and finite capacity. Takine has obtained several relations between queue lengths at random instants and at departure instants; see in particular the very general Theorem 1 in [19], for a single server queue with multiple Markovian arrival streams -an extension of Markovian arrival processes to (possibly correlated) multiple arrival streams.
Organization of the paper. Section 2 provides a short proof of a result in [3] by using a multidimensional queue-length balance argument. Section 3 derives the distributional relationship for an open queueing network under a very general setting in Theorem 1. Extensions to the non-stationary case are discussed in Section 4. Some applications are presented in Section 5. Section 6 contains concluding remarks.

A balance equation for a class of polling models
In this section we provide a simple relation between the steady-state joint queue-length distribution at arbitrary epochs and at departure epochs for polling models. This relation, which is derived by introducing a multi-dimensional queue-length balance argument, is used to provide a short, but somewhat intuitive derivation of Theorem 1 of [3]. In the next section we shall extend that balance equation in a very general setting, and give a rigorous derivation. Let us first describe the polling model studied in [3].
Consider a system of m ≥ 1 infinite-buffer queues Q 1 , . . . , Q m and a single server S. Queues are indexed by J = {1, 2, . . . , m}. The service times of customers in Q i are i.i.d. (independent, identically distributed) positive random variables generically denoted by B i , with means b i := EB i . Denote the Laplace-Stieltjes transform (LST) of B i byB i (·). The server moves among the queues in a cyclic order. When S moves from Q i to Q i+1 , it incurs a switchover period. The durations of successive switchover times are i.i.d. non-negative random variables, which we generically denote by S i . Denote the LST of S i byS i (·) and assume that s i := ES i < ∞; let s := m i=1 s i . Customers arrive at Q i according to a Poisson process with rate λ i ; let λ := m i=1 λ i . We do not assume anything about the service disciplines at Q i . Define ρ i := λ i b i as the traffic intensity at Q i ; let ρ := m i=1 ρ i . We assume that ρ < 1, which is a necessary condition for the system to be stable. In what follows we shall write z for an m-dimensional vector in R m , z = (z 1 , . . . , z m ), and we assume that |z i | ≤ 1 for every i ∈ J. We implicitly use the convention that any index summation is modulo m, for example Assume that all the usual independence assumptions hold between the service times, the switchover times and the interarrival times. We assume that the ergodicity conditions are fulfilled and we restrict ourselves to results for the stationary situation. Now introduce the PGF of various joint queue-length distributions: V b i (z) and V c i (z) denote the PGFs of the joint queue-length distribution at visit beginnings and visit completions at Q i , while S b i (z) and S c i (z) denote the PGFs of the joint queue-length distribution at service beginnings and service completions at Q i ; L(z) denotes the PGF of the joint queue-length distribution at an arbitrary time in steady-state. Theorem 1 of [3] states that, with mean cycle time EC = s 1−ρ : with Σ(z) := m j=1 λ j (1 − z j ). Its proof in [3] is based on the following relations: (i) a balance relation for polling systems, which is due to Eisenberg [8] and which was generalized in [2]: Here γ i := 1/λ i EC represents the reciprocal of the mean number of customers served at Q i per visit, i.e., the long-term ratio of visit beginnings to service beginnings.
(ii) an obvious relation between queue lengths at the beginning and end of a service time: (iii) an obvious relation between queue lengths at the beginning and end of a switchover time: (iv) a stochastic mean value theorem, expressing L(z) as an average over the PGFs of the joint queue-length distribution at an arbitrary moment during a visit to Q i (X i (z)) and during a switchover period between Q i and Q i+1 (Y i (z)): where, for i ∈ J, whereB past i (·) andS past i (·) are the LST's of the past (elapsed) parts of B i and S i , respectively, that is, they are defined as Starting from (5), substituting (6) and (7), and using (2) and (3) to eliminate all S c i (z) and S b i (z), yields (1).

Remark 1
In [3] also zero switchover times are allowed; the same result (1) is shown to hold.
In Theorem 1 of [3] it was subsequently observed that one may simplify (1) as follows, by using (2) and (3): This formula is remarkably simple; please notice that it does not involve the service time distributions, and that the service disciplines at the various queues also do not play a role, which suggests that (1) is based on very general principles. This is the formula for which we would like to provide a short proof (see below). In combination with (2) -(4), it also gives a short proof of (1). In other words, one can obtain an expression for the PGF of the joint steady-state queue-length distribution in a large class of polling systems by just using the elementary balance equations (2) and (9) (see below), combined with the obvious relations (3) and (4).
Secondly, observe that, because of the Poisson arrival processes, L(z) is also the PGF of the joint queue-length distribution just before an arrival at Q i , i ∈ J by PASTA (Poisson Arrival See Time Averages, e.g., see [1,14]). Thirdly, invert the transform expressions on both sides of (9), yielding for x ≥ 0 and e i being the unit vector with 1 in the ith coordinate and zero elsewhere: where π d i (·) indicates that we consider the joint queue-length distribution right after a departure from Q i , and π e i (·) denotes that we view the system just before an external arrival at Q i . Fourthly, we reshuffle the terms: Finally, observe that the lefthand side of (11) represents the rate out of state x, and the righthand side represents the rate into that state. Indeed, the first term in the lefthand side corresponds to arrivals which find x customers in the system. The second term in the lefthand side is slightly less obvious. It corresponds to departures that take place in state x. Notice that the rate at which customers depart from Q i equals λ i (although the departure process will not be a Poisson process), and that π d i (x − e i ) is the fraction of departures from Q i which take the system out of state x. Similarly interpret the terms in the righthand side. We conclude that (8) amounts to a simple flow balance formula.

Remark 2
A similar flow balance argument was used in [5] to derive a queue-length expression in an M/G/1 FCFS queue with multiple customer classes.

Remark 3
Observe that (8) immediately gives the formula for the marginal distributions. Indeed, for a vector z m,i = (1, . . . , 1, z i , 1, . . . , 1), L(z m,i ) = S c i (z m,i ). From the well-known 'step' (level-crossing) argument it follows that S c i (z m,i ) is also the PGF of the queue-length distribution in Q i at an arrival epoch at Q i . By PASTA it is also the PGF of the steady-state distribution of Q i .
Next take z T = (z, . . . , z). (8) now states that the PGF of the distribution of the total queue length (in terms of z) equals m i=1 λ i S c i (z T )/ m j=1 λ j . This formula may be interpreted as follows. By PASTA, L(z T ) is also the PGF of the distribution of the total queue length at an arrival epoch. By a level-crossing argument, it follows that this equals the PGF of the distribution of the total queue length just after a departure epoch. The result now follows from the observation that a fraction λ i / m j=1 λ j of the departure epochs refers to a departure from Q i . (8) may be viewed as an m-dimensional version of the above-mentioned one-dimensional 'step' (level-crossing) relation that holds for queues with single arrivals and single departures.

Formal derivations under a general framework
In this section, we aim to derive distributional relationships at arrival and departure instants for various queues and their network models in a unified way, under general settings. Roughly speaking, these settings allow simultaneous external arrivals, simultaneous departures and routing at different stations; however, we do not allow an external arrival to coincide with a departure. We use their time evolutions in sample paths for deriving the relationships rather than using flow balance.
We describe a queueing network system under a fairly general framework. We consider an open queueing network system with m queues, where queues uniquely belong to service facilities, which are called stations. Queues in the same station may be distinguished by customer classes. Each station may have multiple servers, which may change in time. External arrivals at queues are general as long as they satisfy certain stationarity conditions. Customers completing service may be routed among queues depending on the state of the whole system. Thus, this model is quite general and very flexible.
To describe this model, we introduce a stochastic process. Queues are still indexed by J = {1, 2, . . . , m}. Let where X i (t) represents the length of queue i at time t, which includes customers in service. Here, each queue belongs to a single station. There is a mapping from queues to stations, which will be given when needed.
In addition to X(t), the following counting processes count the number of specified events until time t ≥ 0 for i ∈ J, • N e i (t) -external arrivals at queue i, • N r i (t) -internal arrivals at queue i (transition from some queue).
With N u (t) = (N u 1 (t), . . . , N u m (t)) for u = e, d, r, we consider the process All processes are assumed right-continuous with left limits. Let is similarly defined and is in Z m + for u = e, d, r, where Z + is the set of nonnegative integers.
For u = e, d, r, denote and assume that That is, external arrivals and service completions can not occur simultaneously.
We also need to define the intermediate state It differs from X(t) only at departure epochs and it describes the state "after" a departure and "before" an internal arrival at a different queue. Clearly, the following dynamics hold.
Because of i), X(t) and X d (t) are also finite. It may be natural to assume that |N r |(t) ≤ |N d |(t) for t ≥ 0, but we do not require it in this section.
Thus, X(t) is the state of the system at time t of an input-output system driven by counting processes N e , N d , N r . The dynamics of (12) and (13) indicates that we adopt the departure first framework. We have used queueing terminologies, but our results are valid as long as the above mathematical assumptions and (13) are satisfied.
In general, |N e |(t), |N d |(t) and N e i (t) and N d i (t) may have jumps greater than one, which is not convenient to describe the time evolution of Z(t). Thus, for u = e, d, we introduce Another basic assumption on the counting processes is iii) There exist finite and positive numbers λ u , u = e, d such that a.s. (almost surely) w.r.t. the underlying probability measure P.
We further assume the following ergodic type conditions. iv) There exist probability distributions π e and π d such that y, z), a.s., x, y, z ∈ Z m + .
From the definitions in iv), π e and π d are considered as the embedded stationary distributions just before arrival epochs and just after departure epochs but before internal arrivals, respectively. They correspond to Palm distributions concerning their counting processes in the time stationary framework (e.g., see [1]).
Since the process X(t) is vector valued, it is not so convenient for manipulations. So, we introduce a test function f : Z m + → R. Under the setting i)-iv), we will derive distributional relationships among characteristics at different embedded instants using the test function f . For this, we need the following lemma.

Lemma 1 If (15) holds, then, for any bounded function
Similarly, if (16) holds, then, for any bounded function h : Z 3m + → R, we have This lemma may look obvious, but its proof is not immediate because we need to verify the exchange of limits. We prove it in Appendix A.
We are now ready to prove distributional relationships. First, we denote the expectations under π e and π d by E e and E d , respectively. That is, Note that Y in E e represents sizes of externally arriving batches, while Y in E d represents sizes of departing batches.

Theorem 1 Under the setting i)-iv), for any bounded function
Proof Since f (X(t)) changes in time only at the counting instants t e n or t d n , we have (with ∆ being defined as earlier in this section) ∆f (X(t d ℓ )).
Recalling (12), we have , and this and (13) yield Substituting these X(t d ℓ ) and It follows from (22) and (23) that Dividing both sides of this equation by t and letting t → ∞ yields (21) by (14)- (16) and Lemma 1 because f is bounded.
The assumptions of Theorem 1 exclude arrivals and departures to occur simultaneously, but allow them to occur separately as multiple simultaneous external arrivals or multiple simultaneous departures and routing. The model as well as the distributional relationship may be too general for queueing networks. To make them more specific, we make the following assumption. v) There exist finite and nonnegative numbers λ d A for nonempty A ⊂ J, that is, A ∈ 2 J \ {∅}, such that where, with notation Note thatÑ d A counts instants when departures occur simultaneously from queues i ∈ A but there is no departure from queue j ∈ J \ A, while ∆Ñ d Thus, the setting i)-v) still allows batch arrivals and batch departures and simultaneous transfer of customers in a departing batch.
We will use the following notation. For each A ∈ 2 J \ {∅}, let, for x, y, z ∈ Z m + , SinceÑ d A exclusively counts the increasing epochs of |Ñ d | for different A's, we have which implies that λ d = A∈2 J \{∅} λ d A , and, for A ∈ {B ∈ 2 J |λ d B > 0}, π d A is a probability distribution on Z 3m + , which can be restricted to Z m + × S A × Z m + . Let t d A,n be the n th jump epoch ofÑ d A . Just as Lemma 1 does, the following lemma plays a key role; it is proved in Appendix B.

Corollary 1 Under the setting i)-v), for any bounded function
In this case, the summations over A in (29) can be reduced to those over i ∈ J, replacing A by i.
Until now, our distributional relationship may still be too general because no assumption is made on how the counting processes are generated from X(t) and other information. To describe this, a filtration is convenient. Let F t be the σ-field generated by all events up to time t, and let F t− = σ(∪ u<t F u ), that is, F t− is a σ-field generated by all events before time t. For a stopping time τ , let F τ − = σ(F 0 , {A ∩ {t < τ } ∈ F t }), where σ(A) is the σ-field generated by a family of events A. Using the filtration, the following assumptions are typically used under the setting i)-v).
(a1) t e n , t d i,n are stopping times with respect to {F t ; t ≥ 0}. This can always be realized by choosing a sufficiently large F t .
(a2) ∆N e (t e n ) is independent of F t e n − . That is, the sizes of batch arrivals are independent of the state of the system just before their arrival epochs.
That is, departures singly occur from one queue at a time. .
Under the setting i)-v) and the assumptions (a1)-(a4), ∆N r j (t d i,ℓ ) ≤ 1, and therefore Lemma 2 yields which is denoted by π d ij (x). We here recall that e i ∈ Z m + is the unit vector whose i-th entry is one and the other entries are zero. Thus, applying Corollary 1 for f (x) = z x , where we recall that z x = i∈J z x i i , we have the following relationship.

Remark 6
Under the assumptions of this corollary, the routing of departing customers may depend on all queue lengths in the network.
Corollary 2 is specialized to Corollary 3 if external arrivals to queues occur one at a time. Namely, vi) No simultaneous arrivals occur, and there exist finite numbers (some, but not all, possibly zero) λ e k for k ∈ J such that Corollary 3 Under the assumptions of Corollary 2, assume that vi) also holds. Define the π e k as π e k (x, y k ) = λ e λ e k π e (x, y k ), λ e k > 0, 0, λ e k = 0, then, for k ∈ {i ∈ J|λ e i > 0}, the π e k is a probability distribution on Z m+1 + , and (30) becomes where ϕ e k is the generating function of X under the conditional distribution π e k .
Corollary 3 immediately implies the following corollary.

Remark 7
In Section 5 we shall present several applications of the above theorem and corollaries. In particular, the polling result (8) of Section 2 is there shown to be a special case of Corollary 1.
Notice that the setup of this section includes the finite buffer case. This is done by having no arrivals to a queue during times in which it is saturated. This type of dependence is allowed by our setup. Some results for the single server queue with finite capacity are contained in [12].

Distributional relationship up to a given time
The purpose of this section is to derive a non-stationary version of Theorem 1, a distributional relationship up to a given time. We adopt the settings i)-iv) of Section 3, and consider the process Z(t) introduced in the beginning of that section. We first define the expected relative frequencies for bounded test functions g, h from Z 2m + , Z 3m + to R up to time t as g(X(t e n −), ∆N e (t e n ))1(|Ñ e |(t) > 0), For each bounded function f : Z m + → R, we define the following test functions. Let Then, (24) yields the following lemma.

Lemma 3 Under the setting i)-iv), for any bounded function
We may interpret Lemma 3 as a transient version of Theorem 1. It is notable that (34) holds without any stability condition, and its right-hand side vanishes as t → ∞ at most in linear order of t −1 because f is bounded. If there exists a unique probability measure such that (X(t), ∆N e (t), ∆N d (t), ∆N r (t)) is stationary, then R e t g, R d t h converge to the corresponding expectations under the Palm distributions involving |Ñ e |, |Ñ d |, respectively.

Some special cases and applications
In this section we consider several applications of the theorem and corollaries of Section 3.
We first note that, if nonzero N e k for k ∈ J are independent compound Poisson processes, then by PASTA the embedded stationary distributions π e and π e k are identical with the time stationary distributions.

Case 1: An m-class queue with batch arrivals
We consider an m-class single-node service facility, with m ≥ 1. We allow multiple servers. Customers arrive according to a Poisson process, possibly in batches. Customers of class i require service at the service facility according to service time distribution B i (·), i ∈ J. These distributions are assumed to be continuous, but not otherwise specified. No customers are lost; there is an infinite waiting room. After completion of their service, customers immediately leave. We assume that the steady-state joint queue-length distribution (numbers of customers of all classes in the system) exists. Its PGF is denoted by L(z). We also again (as in Section 2) denote the PGF of the steady-state joint queue-length distributions immediately after departure epochs of a class i-customer by S c i (z), i ∈ J. We do not specify according to which service discipline the customers are served; polling with FCFS within each class is just one of many options.
Theorem 2 Consider the above-described m-class single-node service facility. Assume that customers arrive according to a batch Poisson process with rate λ and that customers are served individually, in some non-specified order. Let an arbitrary batch arrival have size . Then the following relation holds between the PGF L(z) and the PGFs S c i (z), i ∈ J: Proof After using PASTA, Theorem 2 is a special case of (30) of Corollary 2 in which there is no routing.
Remark 8 Special cases of the above theorem are obtained by assuming that batches always contain only customers of one type. For the special case that batches have just one customer of class i with probability λ i λ , i ∈ J, (35) reduces to (8) that was obtained for the polling system that provided the initial motivation for the present study (but (8) obviously holds for a much more general class of service disciplines).
Case 2: Generalization of Theorem 2 to the case of batch services The following theorem generalizes the main result in [11], but is a special case of Theorem 2 of [19] which allows a more general arrival process (but in that theorem batch service is not considered).
Theorem 3 Consider the m-class single-node service facility of Theorem 2, with the additional assumption that customers of class-i are always served in batches of fixed size K i , i ∈ J; the start of a service of class-i customers is delayed until K i customers are present. Then the following relation holds between the PGF L(z) and the PGFs S c i (z), i ∈ J: Proof In view of Remark 5, and after using PASTA, Theorem 3 is a special case of (29) of Corollary 1 in which there is no routing, the external arrival batch Y (t e n ) is independent of F t e n − and the departing batch size Y i from queue i is some constant K i . In this case, it is easy to see that λ d i E[Y i ] = λ e /K i , and we obtain (36) from (29).
Case 3: Non-preemptive priority queues In this example we consider a non-preemptive priority queue with P customer classes. We first verify the equality between the PGFs as given by Theorem 2 for P = 2, and subsequently point out how one may use the theorem to obtain the steady-state joint queue-length distribution in that example for a P -class queue. Consider the M/G/1 queue with P classes of customers, with nonpreemptive priority in descending order 1, 2, . . . , P (so Class 1 has the highest priority). Let λ i denote the arrival rate of customers of class i, i = 1, 2. Takagi ( [18], Formula (2.87) on p. 311) presents the PGF Π(z 1 , z 2 , . . . , z P ) of the steady-state joint queue-length distribution immediately after an arbitrary customer departure epoch. For P = 2 he also obtains the PGF P (z 1 , z 2 , . . . , z P ) of the steady-state joint queue-length distribution at an arbitrary epoch ( [18], Formula (5.82b) on p. 397). We have verified that, indeed, for P = 2 classes one has (cf. Theorem 2 with single arrivals), The starting point for this verification was the obvious set of relations, with β i (z 1 , z 2 ) the PGF of the numbers of arrivals at both queues during one service of a class-i customer, i = 1, 2: cesses, and with continuous service time distributions. We have Markovian routing, a customer moving from Q i to Q k with probability p ik and leaving the system after its service completion in Q i with probability p i0 , i, k ∈ J. Define Λ i as the total flow through Q i per time unit, i ∈ J; these Λ i are the unique solution of the set of equations Let A i indicate that the system is viewed just before an arrival at Q i , D i that the system is viewed just after a departure from Q i , and I ik that the system is viewed just after a departure from Q i and just before the arrival of the departing customer at Q k . Letting j = (j 1 , j 2 , . . . , j m ), one can write down the following balance equations for the queue length vector X = (X 1 , X 2 , . . . , X m ): The (PGF of the) probabilities, given that we observe just after a real departure from Q i or that we observe just after a departure from Q i that will in an instant result in an arrival at Q k , are obviously the same. If one takes PGFs, one quickly sees that a special case of (33) is obtained.
We now use (42) to provide an alternative proof for the joint queue-length distribution in a queueing network with a single roving server as studied in [4,17]. Again consider a network of m queues with Markovian customer routing, as described above. In this particular example, we assume that a single server visits the queues in a fixed, cyclic order, requiring a switch-over time S i to move from Q i to Q i+1 . We do not make any assumptions regarding the service disciplines at each queue. This model, which can be regarded as a polling model with customer routing, has been studied by Sidi, Levy and Fuhrmann [17] who refer to this model as a queueing network with a roving server. Sidi et al. obtain the joint queue-length distribution at arbitrary moments, as well as the joint queue-length distribution at departure epochs. The waiting-time distributions are obtained in a different paper [4]. For us, it is slightly more convenient to refer to this latter paper in the analysis below, because the authors in [4] use the same definition of V c i (z), the PGF of the joint queue length at departure epochs, just after a departure from Q i and just before the arrival of the departing customer at the next queue.
Take the formulas (3.2)-(3.6) of [4]. From (3.2), which is the counterpart of our (2), one can express (in the notation of the present paper) the differences of PGFs at visit beginning and visit completion epochs into those at service beginning and service completion epochs: Here P i (z) := p i0 + m k=1 p ik z k , and EC = s/(1 − ρ) with ρ := m i=1 Λ i b i . Next use our relation (3) to express S b i (z) into S c i (z). Subsequently express L(z), in (3.4) of [4], which is the counterpart of (1) above, in differences V b i (z) − V c i (z), as was also done in [3]. This gives This is indeed in agreement with (42): the LHS of (44) gives the first and the fifth term in (42). The last term in the RHS gives the second plus the third term in (42), once we realize that p i0 + m k=1 p ik = 1, and that the conditional probabilities both refer to a service completion in Q i , no matter whether the condition is D i or I ik . The first term in the RHS gives the fourth plus fifth term in (42). One could argue that some results in [4] and [17] could have been derived faster by starting from (44).

Concluding remarks
This paper derives a distributional relationship, at different embedded epochs, for analyzing queues and their networks. As shown in Section 3, it has different forms according to the abstraction level of the model. This may both lead to new results and easier derivations of some known results. In Section 5 this is demonstrated for a few examples.
The relationship in Section 4 has a different nature than the rest of this paper because it does not require any stationarity of the processes of interest. Namely, it suggests that such an asymptotic relationship may enable us to obtain queueing characteristics with some error bounds, not assuming any stationarity condition. This is completely different from the standard analysis in queueing theory. Thus, it would be interesting to see whether it can yield useful results for the performance evaluation of queueing models. We leave this for future studies.

B Proof of Lemma 2
In view of Lemma 1, to prove (28) it suffices to prove that, for A ∈ 2 J \ {∅}, It follows from v) that, for each i ∈ J, ℓ ≥ 1, y ∈ S A , z ∈ Z m + , there is a unique k ≥ 1 such that ℓ ≤ k and and (14) and (25) imply Hence, for y ∈ S A , π d (x, y, z) = lim This proves (46) by the definition π d A , and therefore (28) holds. The fact that π d A is a probability distribution is immediate from (28) with h(x, y, z) ≡ 1.