Heavy traffic analysis of roving server networks
aa r X i v : . [ m a t h . P R ] N ov Heavy traffic analysis of roving server networks
M.A.A. Boon ∗ R.D. van der Mei † E.M.M. Winands ‡ October 6, 2018
Abstract
This paper studies the heavy-traffic (HT) behaviour of queueing networks with a singleroving server. External customers arrive at the queues according to independent renewalprocesses and after completing service, a customer either leaves the system or is routed toanother queue. This type of customer routing in queueing networks arises very naturallyin many application areas (in production systems, computer- and communication networks,maintenance, etc.). In these networks, the single most important characteristic of the systemperformance is oftentimes the path time, i.e. the total time spent in the system by an arbitrarycustomer traversing a specific path. The current paper presents the first HT asymptotic forthe path-time distribution in queueing networks with a roving server under general renewalarrivals. In particular, we provide a strong conjecture for the system’s behaviour under HTextending the conjecture of Coffman et al. [
8, 9 ] to the roving server setting of the currentpaper. By combining this result with novel light-traffic asymptotics we derive an approxi-mation of the mean path-time for arbitrary values of the load and renewal arrivals. Thisapproximation is not only highly accurate for a wide range of parameter settings, but is alsoexact in various limiting cases. Keywords: queueing network, path times, waiting times, heavy traffic, approximation
Mathematics Subject Classification:
This paper considers heavy-traffic (HT) limits for queueing networks with a single roving serverthat visits the queues in a cyclic order according to the gated and exhaustive service. Customersfrom the outside arrive at the queues according to general renewal processes, and the servicetime and switch-over time distributions are general as well. After receiving service at queue i , acustomer is either routed to queue j with probability p i , j , or leaves the system with probability p i ,0 . This model can be seen as an extension of the classical polling model (in which customersalways leave the system upon completion of their service) by customer routing. ∗ Eurandom and Department of Mathematics and Computer Science, Eindhoven University of Technology, P.O. Box513, 5600MB Eindhoven, The Netherlands. Email: [email protected]. † Department of Mathematics, Section Stochastics, VU University, De Boelelaan 1081a, 1081HV Amsterdam, TheNetherlands and Centre for Mathematics and Computer Science (CWI), 1098 SJ Amsterdam, The Netherlands. Email:[email protected]. ‡ University of Amsterdam, Korteweg-de Vries Institute for Mathematics, Science Park 904, 1098 XH Amsterdam,The Netherlands. Email: [email protected]. [
3, 13, 17, 33 ] for overviews of polling systems andtheir applications). In many applications, however, the interarrival times are not exponentiallydistributed. Therefore, in this paper we study a network in which the arrival process at each ofthe queues is a general renewal process. For open networks an important characteristic of thesystem performance is the path time, defined as the total time spent in the system by an arbitrarycustomer traversing a specific path. That is, oftentimes service level agreements are made on thetotal time required for all the service requests of a customer in a system to be finished. Moreover,in many computer-communication and production-inventory systems the single most importantperformance measure is often not an aggregate measure like the mean path time, rather theprobability that this path time exceeds a pre-defined threshold. In view of dimensioning suchsystems, the importance of the path time distribution as a performance measure of interest isevident. Due to the routing of customers, which leads to non-renewal arrival processes at thequeues and to strong interdependence of the waiting times within a path time, simulation appearsto be the only practical recourse at the present time. In these circumstances, one naturally resortsto asymptotic estimates. In particular, in this paper we study the path-time distribution in aqueueing network with customer routing under HT conditions.The motivation for studying the HT regime - which is also the most challenging regime from ascheduling point of view - is twofold. First, an attractive feature of HT asymptotics is that in manycases they lead to strikingly simple expressions for the performance measures of interest. Thisremarkable simplicity of the HT asymptotics leads to structural insights into the dependence of theperformance measures on the system parameters and gives fundamental insights in the behaviourof the system in general. A second appealing feature of HT asymptotics is that they form anexcellent basis for developing simple accurate approximations for the performance measures forstable systems.The introduced queueing network is very general, which is illustrated by the fact that manyspecial cases have been studied in the past. Some special case configurations are standard pollingsystems [ ] , tandem queues [
20, 35 ] , multi-stage queueing models with parallel queues [ ] ,feedback vacation queues [
7, 34 ] , symmetric feedback polling systems [
32, 34 ] , systems with awaiting room [
1, 31 ] , and many others. Due to the intrinsic complexity of the model, previousstudies on the network in its full generality were restricted to queue lengths and waiting timedistributions for stable systems under the assumption of Poisson arrivals (see [
6, 28, 29, 30 ] ). Al-though the results in these papers are exact they lack an explicit analysis with simple expressions,leading to the need of using numerical techniques to determine performance measures of inter-est. Moreover, these results are limited to the waiting-time distribution instead of the practicallymost interesting and theoretically most challenging path-time distribution.Besides that we have a theoretical interest in the proposed queueing network, the presentwork is motivated by the fact that customer routing in polling systems arises naturally in a host ofapplication areas (in production systems, computer- and communication networks, maintenance,etc.). Some examples are a manufacturing system where products undergo service in a numberof stages or in the context of rework [ ] , a Ferry based Wireless Local Area Network (FWLAN)in which nodes can communicate with each other or with the outer world via a message ferry [ ] , a dynamic order picking system where the order picker drops off the picked items at thedepot where sorting of the items is performed [ ] , and an internal mail delivery system wherea clerk continuously makes rounds within the offices to pick up, sort and deliver mail [ ] .Motivated by the attractiveness of HT asymptotics, several approaches have been proposed toobtain HT limits for polling systems. HT limits have been rigorously proven for systems with Pois-2on arrivals (cf. [
21, 37 ] ) and renewal arrivals (cf. [
8, 9, 14, 38 ] . A central role in these papers isthe Heavy-Traffic Averaging Principle (HTAP) which means that the total scaled workload may beconsidered as a constant during a cycle, whereas the workload of the individual queues changemuch faster according to deterministic trajectories, or a fluid model. In this paper, we also followthis well established path of deriving HT limits for systems with general renewal arrivals, see also [
18, 19, 24, 25, 26 ] .One of the main contributions of the current paper is that we present the first HT asymptoticfor the path-time distribution in queueing networks with a roving server under general renewalarrivals. The main building blocks of this path-time distribution are inevitably the waiting-timedistributions at the individual queues, for which the current paper also presents the first exactHT asymptotics . In particular, we provide a strong conjecture for the systems behaviour underHT extending the polling conjecture of Coffman et al. [
8, 9 ] to the roving server setting of thecurrent paper. Our conjecture is validated in three ways. Firstly, as stated before we followa well established and accepted line of thinking. Secondly, our conjecture corresponds to theaforementioned rigorously proven distributional limits in special cases of our network. Thirdly,we give some numerical examples that illustrate that the correct limiting behaviour has beenderived. As an important by-product of our analytical framework, we obtain exact asymptotics inthe large switch-over time regime as well.The second main contribution concerns the derivation of a simple approximation of the meanpath-time for arbitrary values of the load and renewal arrivals by combining the HT results withnewly derived light-traffic (LT) asymptotics. This approximation technique is shown to be exact invarious limiting cases, and is known to be highly accurate for a wide range of parameter settings( [
5, 10 ] ). Moreover, it satisfies the Pseudo Conservation Law (PCL), and consequently it leads toexact closed-form results of the mean waiting time for symmetric systems with Poisson arrivals.The resulting expressions are very insightful, simple to implement, and suitable for optimizationpurposes. In particular, the approximation shows explicitly how the path times depend on thesystem parameters such as the routing probabilities p i , j .Our presentation of the analysis, and the analysis itself, may be considered to be somewhatinformal throughout. Providing a rigorous presentation of our results, however, would be anextremely interesting, but notoriously difficult area for further research. Furthermore, it wouldtake us far afield from our main goal, i.e. to obtain fundamental insights into customer routing insystems. Lastly, we would like to stress that the current paper concerns a continuous-time cyclicsystem with gated or exhaustive service in each queue, but that all results can be extended todiscrete time, to periodic polling, to batch arrivals, or to systems with different branching-typeservice disciplines such as globally gated service.The structure of the present paper is as follows. In Section 2 we introduce the model and therequired notation. In Section 3 we use the HTAP to derive exact queue-length, waiting-time, andpath-time asymptotics for networks with gated and / or exhaustive service. Based on these resultswe develop novel approximations in Section 4. In the penultimate section, we present somepractical cases that illustrate the versatility of the queueing network and the importance of thepath-time distribution in practice. Finally, we present some conclusions and points for discussionin the last section. Some preliminary results restricted to mean waiting times were derived in [ ] . Model and notation
In this paper we consider a queueing network consisting of N ≥ Q , . . . , Q N .External customers arrive at Q i according to a general renewal arrival process with rate λ i , andhave a generally distributed service requirement B i at Q i , with mean value b i : = E [ B i ] , and sec-ond moment b ( ) i : = E [ B i ] . Throughout this paper we assume that all random variables havefinite second moments. The queues are served by a single server in cyclic order. Whenever theserver switches from Q i to Q i + , a switch-over time R i is incurred, with mean r i . The cycle time C i is the time between successive moments when the server arrives at Q i . The total switch-over timein a cycle is denoted by R = P Ni = R i and its first two moments are r : = E [ R ] and r ( ) : = E [ R ] .Indices throughout the paper are modulo N , so Q N + actually refers to Q . All service timesand switch-over times are mutually independent. Each queue receives exhaustive or gated ser-vice. Exhaustive service means that each queue is served until no customers are present anymore,whereas gated service means that only those customers present at the server’s arrival at Q i will beserved before the server switches to the next queue. This queueing network can be modelled as a polling system with the specific feature that it allows for routing of the customers: upon comple-tion of service at Q i , a customer is either routed to Q j with probability p i , j , or leaves the systemwith probability p i ,0 . Note that P Nj = p i , j = i , and that the transition of a customerfrom Q i to Q j takes no time. The model under consideration has a branching structure, which isdiscussed in more detail by Resing [ ] . The total arrival rate at Q i is denoted by γ i , which isthe unique solution of the following set of linear equations: γ i = λ i + P Nj = γ j p j , i , i =
1, . . . , N .The offered load to Q i is ρ i : = γ i b i and the total utilisation is ρ : = P Ni = ρ i . We assume thatthe system is stable, which means that ρ should be less than one (see [ ] ). The total servicetime of a customer is the total amount of service given during the presence of the customer in thenetwork, denoted by e B i , and its first two moments by ˜ b i : = E [ e B i ] and ˜ b ( ) i : = E [ e B i ] . The firsttwo moments are uniquely determined by the following set of linear equations: For i =
1, . . . , N ,˜ b i = b i + N X j = ˜ b j p i , j , (2.1)˜ b ( ) i = b ( ) i + b i N X j = ˜ b j p i , j + N X j = ˜ b ( ) j p i , j . (2.2)We study this model under heavy-traffic conditions, i.e., we increase the load of the systemuntil it reaches the point of saturation, ρ ↑
1. As the total load of the system increases, the visittimes, cycle times, and waiting times become larger and will eventually grow to infinity. Forthis reason, we scale them appropriately and consider the scaled versions. We consider severalvariables as a function of the load ρ in the system. For each variable x that is a function of ρ , we denote its value evaluated at ρ = ˆ x . Scaling is done by varying the interarrivaltimes of the external customers. To be precise, the limit is taken such that the external arrivalrates λ , . . . , λ N are increased, while keeping the service and switch-over time distributions, therouting probabilities and the ratios between these arrival rates fixed. For ρ =
1, the genericinterarrival time of the stream in Q i is denoted by ˆ A i . Reducing the load ρ is done by scaling theinterarrival times, i.e., taking the random variable A i : = ˆ A i /ρ as generic interarrival time at Q i .Hence, the rate of the arrival stream at Q i satisfies λ i = / E [ A i ] . Furthermore, we define arrivalrates ˆ λ i = / E [ ˆ A i ] , and proportional load at Q i , ˆ ρ i = ρ i /ρ (“proportional” because P Ni = ˆ ρ i = HT asymptotics
To obtain HT-results for the waiting-time distributions, we use HT results for polling systemswithout customer routing, which are obtained by [
8, 9, 14 ] . The key observation in these papersis the occurrence of a so-called Heavy Traffic Averaging Principle (HTAP). When a polling systembecomes saturated, two limiting processes take place. Let V denote the total workload of thesystem, i.e., the total service requirement of all customers present in the system including thepossible residual service time of a customer being served. As the load offered to the system, ρ ,tends to 1, the scaled total workload ( − ρ ) V tends to a Bessel-type diffusion. However, the work in each queue is emptied and refilled at a faster rate than the rate at which the total workload ischanging. This implies that during the course of a cycle, the total workload can be considered asconstant, while the workloads of the individual queues fluctuate according to a fluid model. TheHTAP relates these two limiting processes.Although rigorous proofs have only been presented for standard polling models (see [
8, 9,14 ] ), results in [
18, 19, 24, 25, 26 ] support the conjecture that the HTAP holds for a muchwider class of systems. As in these papers, we make the crucial assumption that the HTAP holdswithout providing a rigorous proof of convergence. That is, the HTAP occurs due to a time scaledecomposition that is inherent in the heavy traffic scaling. Therefore, in HT the multi-dimensionalindividual workload processes move along a path in the constant workload hyperplane. Forsystems without routing and switch-over times, this path is described in detail by Jennings [ ] .Using the results from this section, this can easily be adapted to our model with customer routing.Furthermore, it is known that the scaled total workload tends to a Bessel-type limit as the systembecomes saturated. More colloquially, the routing of customers impacts the individual workloads,which impels us to considerably modify and extend the HT analysis of polling models, but doesnot affect the time scale decomposition which directly implies that the HTAP should also hold inthe current setting.We provide justification for this conjecture in four ways. Firstly, we follow a well-founded lineof thinking. Secondly, our conjecture corresponds to the rigorously proven distributional limit inthe special case of a two-queue polling model and branching-type polling models with Poissonarrivals. Thirdly, in the next section we present numerical examples supporting the conjecturethat the correct limiting behaviour has been derived. Finally, in the appendix, we provide atheorem and proof for the limiting queue-length distributions under the assumption of Poissonarrivals. In the next subsection we start by discussing the fluid model and subsequently discussthe limiting distribution of the scaled total workload. We use these results to obtain the HT limitof the scaled waiting time and path-time distributions. Workload.
In this section we consider the fluid version of the queueing network with a singleshared server, where the work travels as fluid from one station to another. This fluid modelis carefully constructed in such a way that its behaviour corresponds to the multi-dimensionalindividual workload processes, moving along a path in the constant workload hyperplane asdiscussed in the introduction of this section. An important aspect of this “corresponding fluidmodel” is the absence of switch-over times, as they become negligible under HT conditions. Westart by studying the fluid limit of the per-queue workload, which is obtained by multiplying by ( − ρ ) and letting ρ ↑
1. For our model, the fluid limit of the workload at Q i is a piecewiselinear function. During the visit time of Q k , denoted by V k , k =
1, . . . , N , external fluid particles,5orresponding to the customers in the original model, flow into Q i at rate ˆ λ i . Each of these fluidparticles brings along ˜ b i units of work into the system. Simultaneously, work is being processed in Q k at rate one. Since P Ni = ˆ λ i ˜ b i =
1, the total workload remains constant throughout the courseof a cycle. In steady state, the length of one cycle is fixed and denoted by c . In this paragraph wewill show that there is a simple linear relation between c and the total workload in the system,denoted by v . For this reason, we may regard c as a system parameter in the fluid model, insteadof v . PSfrag replacements V i V i + V i + V i + V i + N − δ i c ˆ γ i ˜ b i c c ˆ ρ i ˆ γ i , i ˜ b i c Figure 1: Mean amount of work in Q i in the fluid limit that arises when the system is in heavytraffic. The length of one cycle is c .Although work is processed at rate one, due to the internal routing work is flowing out of Q k at rate 1 + b k N X i = p k , i ˜ b i = ˜ b k b k ,which is greater than (or equal to) one. The reason for this anomaly is that work decreases in Q k either because of the service of fluid particles (customers) in this queue, or because work isshifted due to internal routing of fluid. Work including rerouted fluid particles is flowing into Q i ,during V k , at rate ˆ γ i , k ˜ b i , where ˆ γ i , k : = ˆ λ i + p k , i / b k , for i , k =
1, . . . , N . It is straightforward toverify that ˜ b k / b k = P Ni = ˆ γ i , k ˜ b i . Figure 1 depicts a graphical representation of the mean amountof work in Q i in the fluid limit throughout the course of a cycle, the length of which is a constant,denoted by c . Using the fact that fluid particles (external and internal) flow into Q i during V k atrate ˆ γ i , k and the fact that the length of V k in the fluid model is equal to ˆ ρ k c , one can show thatthe fluid limit of the mean amount of work in Q i at the beginning of a visit to Q j is j − X k = i ˆ ρ k ˆ γ i , k ˜ b i c for j = i +
1, . . . , i + N . (3.1)This reduces to ˆ γ i ˜ b i c for j = i + N . We have used that in the fluid limit the fraction of time thatthe server is visiting Q j is ˆ ρ j ( j =
1, . . . , N ) . Combining these observations, one can obtain thefollowing expression for δ i , defined as the ratio of the fluid limit of the average amount of workat Q i and the length of a cycle (see Figure 1). 6 efinition 3.1 For i =
1, . . . ,
N , δ i = ˆ ρ i ˜ b i (ˆ γ i + ˆ ρ i ˆ γ i , i ) + i + N − X j = i + ˆ ρ j ˆ ρ j ˜ b i ˆ γ i , j + j − X k = i ˆ ρ k ˜ b i ˆ γ i , k ! . (3.2)We introduce the notation v i for the average amount of work in Q i , v i = δ i c .As the total inflow in all queues is equal to the total outflow per time unit, the total amount ofwork during a cycle remains constant at level v = N X i = v i = N X i = δ i c = δ c ,where δ = P Ni = δ i . Amount of fluid.
We first study the number of fluid particles in each queue at various epochsduring the cycle. We denote by X fluidi , k the number of fluid particles in Q i at the beginning of a visitto Q k . Using (3.1) and the fact that the amount of fluid in Q i is equal to the amount of work in Q i divided by ˜ b i , we obtain X fluidi , k = k − X j = i ˆ ρ j ˆ γ i , j c for k = i +
1, . . . , i + N .Again, note that X fluidi , i + N can also be written as X fluidi , i = ˆ γ i c .As a consequence, the amount of fluid at an arbitrary moment during V k , denoted by L fluidi , k , isuniformly distributed on the interval hP k − j = i ˆ ρ j ˆ γ i , j c , P kj = i ˆ ρ j ˆ γ i , j c i for k = i +
1, . . . , i + N − ” ˆ ρ i ˆ γ i , i c , ˆ γ i c — for k = i . (3.3)We can now obtain an expression for the distribution of the amount of fluid in Q i at an arbitraryepoch, by conditioning on the visit period: L fluidi d = L fluidi , k w.p. ˆ ρ k , (3.4)where L fluidi , k is uniformly distributed as in (3.3). For notational reasons, we introduce the notation L fluidi : = L fluidi / c and L fluidi , k : = L fluidi , k / c . One can consider L fluidi as a standardised version of L fluidi ,not depending on the cycle length c . 7 aiting times. For the fluid model under consideration we are interested in the waiting timedistribution of an arbitrary fluid particle, internal or external. We define the waiting time as thetime between the arrival in a queue, and the moment of departure from this queue (even if theparticle is routed to another, or even the same queue). If we condition on the event that anarbitrary fluid particle arrives in Q i during V k , its conditional waiting time consists of the residualpart of V k , the visit periods V k + , . . . , V i − , and the processing of the amount of fluid that hasarrived in Q i during the elapsed part of the cycle, i.e., V i , . . . , V k − plus the elapsed part of V k . Let u k be the fraction of V k that has elapsed at the arrival epoch of a fluid particle in Q i . We denoteby W fluidi , k ( u k ) , for 0 ≤ u k ≤
1, the conditional waiting time of an arbitrary fluid particle arrivingin Q i during V k , at the moment that a fraction u k of this visit period has elapsed. We have W fluidi , k ( u k ) = k − X j = i − N ˆ ρ j c ˆ γ i , j b i + u k ˆ ρ k c ˆ γ i , k b i | {z } P i , k ( u k ) + ( − u k ) ˆ ρ k c + i − X j = k + ˆ ρ j c | {z } R i , k ( u k ) , (3.5)for i =
1, . . . , N , k = i − N , . . . , i −
1, and 0 ≤ u k ≤
1. As (3.5) indicates, we split W fluidi , k ( u k ) into two parts, namely P i , k ( u k ) representing the waiting time due to the customers that arrivedduring the elapsed (Past) part of the cycle c , and R i , k ( u k ) , which is the time until the next visitof the server to Q i (Residual cycle). This division into two parts turns out to be useful in the nextparagraph, for computing the path-time distributions.The waiting-time distribution follows after unconditioning. During V k fluid flows into Q i atrate ˆ γ i , k . Hence, the probability that an arbitrary fluid particle arrives during V k , given that it ar-rives in Q i , is π i , k : = ˆ γ i , k ˆ ρ k / ˆ γ i . The elapsed fraction of the visit period V k is uniformly distributedon [
0, 1 ] . These two results yield the following expression for the waiting-time distribution of anarbitrary fluid particle in Q i : W fluidi d = P i , k ( U k ) + R i , k ( U k ) w.p. π i , k = c (cid:16) + k − X j = i − N ˆ ρ j (ˆ γ i , j b i − ) + U k ˆ ρ k (ˆ γ i , k b i − ) (cid:17) w.p. π i , k , (3.6)for i =
1, . . . , N and k = i − N , . . . , i −
1. The random variables U , . . . , U N are independent andUniform [
0, 1 ] distributed. As before, we introduce the notation W fluidi : = W fluidi / c to represent astandardised version of W fluidi , not depending on the cycle length c . Path times.
In this paragraph we derive the path-time (the total time spent in the system)distribution of customers - or, in this case, fluid particles - traversing a specific path throughthe network. We denote the time spent in the system by a fluid particle traversing the path Q i , Q i , . . . , Q i M by W fluidi , i ,..., i M , where i k ∈ {
1, 2, . . . , N } for k =
1, 2, . . . , M and M ≥
1. Whenconsidering path times, each fluid particle enters the system as an external fluid particle in Q i .Just like in the previous paragraph, we condition on the visit period and the length of the elapsedpart of this visit period at its arrival epoch. Assume that a tagged fluid particle arrives in Q i atthe moment that the server is visiting Q k and a fraction u k of V k has elapsed ( ≤ u k ≤ ) . Wedenote its conditional path time as W fluidi , i ,..., i M ; k ( u k ) . Its time spent in Q i is exactly W fluidi , k ( u k ) . Atthe moment that the tagged fluid particle leaves Q i and is routed to Q i , the elapsed time of V i P i , k ( u k ) . As a consequence, the elapsed fraction of V i is P i , k ( u k ) / ( ˆ ρ i c ) . This result leads tothe following recursive equation for the conditional path time, W fluidi , i ,..., i M ; k ( u k ) = W fluidi , k ( u k ) + W fluidi ,..., i M ; i ‚ P i , k ( u k )ˆ ρ i c Œ .We define W fluidi M ; i M − ( · ) : = W fluidi M , i M − ( · ) to ensure that the recursion always ends.The path-time distribution follows after unconditioning, noting that the probability that an external fluid particle enters the system during V k is ˆ ρ k : W fluidi , i ,..., i M d = W fluidi , i ,..., i M ; k ( U k ) w.p. ˆ ρ k , (3.7)for i , i , . . . , i M ∈ {
1, 2, . . . , N } , and k =
1, 2, . . . , N . The random variables U , . . . , U N are in-dependent and Uniform [
0, 1 ] distributed. Note that W fluidi , i ,..., i M / c does not depend on the cy-cle time c anymore. Similar to the standardised waiting time, we now introduce the notation W fluidi , i ,..., i M : = W fluidi , i ,..., i M / c , which turns out to be useful later. In this subsection we briefly discuss the fluid model when some of the queues are served exhaus-tively. Rather than repeating all of the analysis of the previous subsection, we mainly focus onthe differences compared to a model with gated service. Note that changing the service disciplineof a particular queue has no effects on any of the other queues.
Workload.
Assume that the service discipline at a certain queue, say Q e with e =
1, . . . , N , isexhaustive service. When comparing the fluid trajectory of the workload of Q e during the courseof a cycle to the case with gated service (as depicted in Figure 1), the only change is that thetrajectory needs to be moved downwards such that the queue is empty at the end of the visitperiod. More precisely, if Q e receives exhaustive service, then we define δ e = δ gatede − ˆ ρ e ˆ γ e , e ˜ b e , (3.8)where δ gatede is the value of δ e given by (3.2) for the case that Q e would have received gatedservice. Amount of fluid.
Since the amount of fluid in Q i is equal to the amount of work in Q i divided by˜ b i , and since we have just established that the amount of work for an exhaustively served queue Q e is equal to the amount of work this queue would have under gated service, minus ˆ ρ e ˆ γ e , e ˜ b e ,we directly obtain that the amount of fluid at an arbitrary moment during V k , denoted by L fluide , k ,is uniformly distributed on the interval hP k − j = e + ˆ ρ j ˆ γ i , j c , P kj = e + ˆ ρ j ˆ γ i , j c i for k = e +
1, . . . , e + N − ” (ˆ γ e − ˆ ρ e ˆ γ e , e ) c — for k = e . (3.9)9 aiting times. Determining the waiting-time distribution of an arbitrary fluid particle also in-volves the same steps as in the gated case, except for fluid particles arriving in Q e during V e ,obviously. W fluide , k ( u k ) = ( − u k ) ˆ ρ k c + e − X j = k + ˆ ρ j c | {z } R e , k ( u k ) + k − X j = e − N + ˆ ρ j c ˆ γ e , j b e + u k ˆ ρ k c ˆ γ e , k b e | {z } P e , k ( u k ) ( k = e ) , ( − u e )(ˆ γ e c − ˆ γ e , e ˆ ρ e c ) b e | {z } P e , e ( u e ) , ( k = e ) .Note that R e , e ( u e ) = Q e is exhaustive. Path times.
We consider the conditional path time W fluide , i ,..., i M ; k ( u k ) , where Q e receives exhaustiveservice. As in the previous paragraph, we need to treat the case e = k separately. A fluid particlearriving in Q e during V e , given that a fraction u e of this visit period has elapsed, will wait for W fluide , e ( u e ) before leaving this queue. This waiting time, which can also be denoted as P e , e ( u e ) ,is a fraction P e , e ( u e )ˆ ρ e c of the visit period V e . This gives the following expression for the conditionalpath time: W fluide , i ,..., i M ; k ( u k ) = W fluide , k ( u k ) + W fluidi ,..., i M ; e (cid:16) P e , k ( u k )ˆ ρ e c (cid:17) , ( e = k ) , W fluide , e ( u e ) + W fluidi ,..., i M ; e (cid:16) u e + P e , e ( u e )ˆ ρ e c (cid:17) , ( e = k ) . In this section we expand upon the HTAP for polling models by relating the limit processes ofthe total workload process and the workload of the individual queues. The main difference withpolling models is that in the roving server network, (the direction of) the shifting of individualworkloads is not only determined by which queue is being served but also by the internal routingof the customers. To this end, we return to the original model under HT conditions.
Workload.
We denote by V the total amount of work in the system at the start of a cycle starting,without loss of generality, at the beginning of visit period V . As far as the total amount of workis concerned, the system behaves like a polling system in heavy traffic with external customersbringing in an amount of work e B i in Q i , but with work shifting from one queue to another uponthe service completion of a customer. For polling systems with general renewal arrivals the HTlimit of the scaled total amount of work at the beginning of a cycle is conjectured by Olsen andVan der Mei [ ] . An adaptation of the conjecture in [ ] , in accordance with the proof for thecase of Poisson arrivals in [ ] , to our model leads to the following result. Conjecture 3.2
Define σ = N X i = ˆ λ i € Var [ e B i ] + (ˆ λ i ˜ b i ) Var [ ˆ A i ] Š , α = r δ/σ , µ = /σ ,10 ith δ as defined in Definition 3.1. Then, for ρ ↑ , ( − ρ ) V has a Gamma distribution with shapeparameter α and rate parameter µ . For more details we refer to [ ] (who, in turn, refer to a result from [ ] ). Cycle time.
Subsequently, the diffusion limit of the total workload process and the workloadin the individual queues can be related using the HTAP. To this end, we start with the cycle-timedistribution under HT scalings, which follows from Conjecture 3.2 and the fluid analysis carriedout in the first part of this section. The length of a cycle depends on the amount of work atthe beginning of that cycle (which may be any arbitrarily chosen moment). Denote by C ( x ) thelength of a cycle, given that a total amount of x work is present at its beginning. In steady state,we have the following relation δ C ( x ) = x , (3.10)where δ is defined as in Definition 3.1. Hence, given an amount of work x , the cycle time is C ( x ) = x /δ . In the fluid model, all cycles have the same length. However, in the stochasticmodel, the cycle lengths are random. Note that, although the cycle-time distribution depends onthe service disciplines and the chosen starting point of the cycle, the mean cycle time is always E [ C i ] = E [ C ] = r / ( − ρ ) (cf. [ ] ). We are now ready to formulate the second conjecture,concerning the limiting distribution of the scaled length-biased cycle time. Conjecture 3.3
For ρ ↑ , we find that ( − ρ ) C i converges in distribution to a random variablehaving a Gamma distribution with shape parameter α and rate parameter δµ . Queue lengths.
Given the cycle-time distribution, we can finally find the scaled queue-lengthdistributions under HT conditions. We use the fluid analysis, in combination with the conjecturedcycle time distribution, to find the limiting distribution of the scaled queue lengths. In the fluidanalysis the cycle time had a fixed length c . Due to the HTAP we can replace the constantcycle time from the fluid analysis by the random variable C i , the scaled length-biased cycle time.If a random variable X has a Gamma distribution with parameters a and b , its length-biasedequivalent X has a Gamma distribution with parameters a + b . Obviously, the replacementof c by C i can only be carried out because of the independence between the length of the cycletime and the uniformly distributed random variables appearing in (3.4). The following conjecturesummarises this result. A theorem and proof for the scaled queue length distribution under theassumption of Poisson arrivals can be found in the Appendix.
Conjecture 3.4 As ρ ↑ , the scaled queue length ( − ρ ) L i converges in distribution to the productof two independent random variables. The first has the same distribution as L fluidi , and the secondrandom variable Γ has the same distribution as the limiting distribution of the scaled length-biasedcycle time, ( − ρ ) C i . For i =
1, . . . , N ; k = i , . . . , i + N − , and ρ ↑ , ( − ρ ) L i d → Γ × L fluidi , k w.p. ˆ ρ k , (3.11) where Γ is a random variable having a Gamma distribution with parameters α + and δµ , and L fluidi , k is the “standardised” number of a type-i particles in the fluid model during V k , introducedbelow (3.4) . The random variables Γ and L fluidi , k are independent. aiting times. The distributions of the scaled waiting times are obtained in a similar manner,exploiting the HTAP to combine the results from the fluid analysis and the scaled, length-biasedcycle-time distribution, which leads to the following conjecture.
Conjecture 3.5 As ρ ↑ , the scaled waiting time ( − ρ ) W i converges in distribution to the productof two independent random variables. The first has the same distribution as W fluidi , and the secondrandom variable Γ has the same distribution as the limiting distribution of the scaled length-biasedcycle time, ( − ρ ) C i . For i =
1, . . . , N ; k = i − N , . . . , i − , and ρ ↑ , ( − ρ ) W i d → Γ × W fluidi , (3.12) where Γ is a random variable having a Gamma distribution with parameters α + and δµ , and W fluidi is the “standardised” waiting time of a type-i particle in the fluid model, introduced below (3.6) . The two random variables Γ and W fluidi are independent. The (HT limit of the) mean waiting time of an arbitrary customer in Q i obviously follows from(3.12), but an easier way to find it, is by application of Little’s Law to the mean queue length at Q i , which is simply the mean amount of work in Q i divided by the mean total service time. Corollary 3.6
For i =
1, . . . ,
N , ( − ρ ) E [ W i ] → ‚ r + σ δ Œ δ i ˆ γ i ˜ b i , ( ρ ↑ ) . (3.13) Path times.
The derivation of the limiting distribution of the scaled path times proceeds alongthe exact same lines, starting with the fluid model. Due to the HTAP, we may again replace theconstant cycle time c from the fluid analysis by the random variable C i to obtain the limitingdistribution of the scaled path times. The conjecture below summarises this result, which isconsistent with the heavy-traffic snapshot principle [ ] . Conjecture 3.7 As ρ ↑ , the scaled path time ( − ρ ) W i , i ,..., i M converges in distribution to theproduct of a random variable having the same distribution as W fluidi , i ,..., i M and a random variable Γ having the same distribution as the limiting distribution of the scaled length-biased cycle time, ( − ρ ) C i . For i , k =
1, . . . ,
N , and ρ ↑ , ( − ρ ) W i , i ,..., i M d → Γ × W fluidi , i ,..., i M , (3.14) where Γ is a random variable having a Gamma distribution with parameters α + and δµ , and thedistribution of W fluidi , i ,..., i M is given by (3.7) . In Subsection 3.3 we have used the HTAP to easily derive the limiting distributions of the scaledworkload at cycle beginnings, the cycle times, the waiting times, and path times. This principlecan be used in the exact same manner when some of the queues receive exhaustive service.Therefore, we will not repeat all these steps to summarise the results for systems with exhaustiveservice. Basically, one only has to take the expressions from the fluid model, and replace theconstant cycle time c by a random variable with the same limiting distribution as ( − ρ ) C i , i.e.,a Gamma distribution with parameters α + δµ (as defined earlier in this section, taking δ : = δ e as defined in (3.8) for exhaustive service).12 .5 Poisson arrivals In the current paper we have derived the system behaviour under heavy traffic for systems withgeneral renewal arrival processes based on the partially conjectured HTAP. Van der Mei [ ] has developed a unifying framework to derive rigorous proofs of the heavy-traffic behaviourof branching-type polling models with (compound) Poisson arrivals. By applying this stepwiseapproach in conjunction with the results of the previous section to the model under consideration,one can rigorously prove the HT asymptotics in queueing networks served by a single sharedserver under the assumption of Poisson arrivals. In the Appendix we provide such a proof for thequeue-length distributions.
In HT the system reaches saturation due to an increase in the total utilisation ρ . However, thesystem might also get saturated due to an increase of the total switch-over time r . These twoasymptotic regimes show, however, significantly different behaviour. In [
40, 41 ] it was shown forpolling systems that the scaled cycle and intervisit times converge in probability to deterministicquantities in the case that the (deterministic) switch-over times tend to infinity. One has tocompare this with the Gamma distribution which is prevalent in the scaled cycle time in thediffusion limit of the present section. The results for polling systems with increasing switch-overtimes of [
40, 41 ] can be extended to the setting of the current paper. That is, as a consequence ofthe scaled cycle time converging to a constant, a fluid limit is obtained implying that the scaleddelay converges in distribution to a mixture of uniform distributions (cf. Formula (3.6)). The HT distributions derived in the preceding section may be used directly as an approximationfor the waiting-time and path-time distributions in non-heavy-traffic systems. However, they tendto perform poorly under low or moderate traffic. Therefore, in this section we combine these HTasymptotics with newly developed LT results leading to highly accurate approximations for themean performance measures for the whole range of load values. We will assess the accuracy ofthe resulting approximations in Section 5.
In order to derive an approximation for the mean waiting time, we study the LT limit of W i whichcan be found by conditioning on the customer type (external or internally routed). Theorem 4.1
For i =
1, . . . ,
N ,W i d → ( R res w.p. λ i /γ i , R j , i w.p. γ j p j , i /γ i , ( ρ ↓ ) , (4.1) where R res is a residual total switch-over time, with probability density functionf R res ( t ) = − P ( R ≤ t ) E [ R ] ,13 nd R j , i = R j + R j + + · · · + R i − is the sum of the switch-over times between Q j and Q i , which is acyclic sum if i ≤ j. Only if i = j and service in Q i is exhaustive, we have R j , i = . Proof:
In LT we ignore all O ( ρ ) terms, which implies that we can consider a customer as being alonein the system. Equation (4.1) can be interpreted as follows. An arbitrary customer in Q i has arrivedfrom outside the network with probability λ i /γ i . In this case he has to wait for a residual totalswitch-over time R res . If a customer in Q i arrives after being served in another queue, say Q j (withprobability γ j p j , i /γ i ), he has to wait for the mean switch-over times R j , . . . , R i − . ƒ The LT limit of the mean waiting time directly follows from (4.1): E [ W i ] → λ i γ i r ( ) r + i − X j = i ∗ γ j p j , i γ i i − X k = j r k , ( ρ ↓ ) , (4.2)where i ∗ = i − N if Q i has gated service, and i ∗ = i − N + Q i has exhaustive service. Sub-sequently, we construct an interpolation between the LT and HT limits that can be used as anapproximation for the mean waiting time for arbitrary ρ . For i =
1, . . . , N , E [ W approxi ] = w LTi + ( w HTi − w LTi ) ρ − ρ , (4.3)where w LTi and w HTi are the LT and HT limits of the mean waiting time respectively, as given in(4.2) and (3.13). Because of the way E [ W approxi ] is constructed, it has the nice properties thatit is exact as ρ ↓ ρ ↑
1. Furthermore, if we have Poisson arrivals, it satisfies a so-calledpseudo-conservation law for the mean waiting times, which is derived in [ ] . This implies thatthe E [ W approxi ] yields exact results for symmetric (and, hence, single-queue) systems. Finally, itcan be shown that this approximation is exact in the limiting case of deterministic set up timesthat tend to infinity (see [
40, 41 ] ).The astute reader has already noticed that the LT result (4.2) is a first-order Taylor expansionof the mean waiting time at ρ =
0, which can be naturally extended with the m th derivatives ofthe mean waiting time with respect to ρ at ρ =
0. Together with the HT limit one has m + ( m + ) th degree polynomial interpolation (cf. [ ] ). Therefore, it is not inconceivable that the approximation can be refined further, but sincethe primary goal of this paper has been the derivation of the path-time time distributions underHT conditions such refinements are beyond the scope of the paper. Moreover, the presentedfirst-order polynomial interpolation is already quite accurate as can be seen in the numericalevaluation. The principle used to determine waiting-time approximations, can be applied to path times inthe exact same manner. For reasons of compactness, we will not present the details in this paper.Almost all of the required ingredients to develop a path-time approximation have been discussed.The only missing piece of information is the LT limit of the mean path-time. However, it caneasily be found given the LT limit of the mean waiting time (4.2): E [ W i , i ,..., i M ] → r ( ) r + b i + M X j = (cid:16) r i j − , i j + b i j (cid:17) , ( ρ ↓ ) , (4.4)14here r j , i = E [ R j , i ] , with R j , i as defined in (4.1).In the current section we have focussed on an interpolation of the mean waiting time and pathtimes for the ease of presentation. It goes without saying that following the same recipe one couldderive an identical interpolation approximation for higher moments as well and, subsequently,fit a phase-type distribution for the complete distribution. Alternatively, one could follow theidea of [ ] , in which a redined HT distribution is derived as an approximation of the completedistribution for arbitrary loads. However, due to the strong correlation between the waiting timesof a customer within a specific path, we have observed in numerical tests that the accuracy of thisapproximation technique is not as high as for standard polling models. We would like to endwith noting that for the mean path times these correlations obviously do no play a role, whencethe developed approximation for the mean figures gives excellent results as illustrated in the nextsection. In the current section we present four practical cases from completely different application areas,indicating the versatility of the studied roving server network and showing the practical usage ofthe developed asymptotics and approximations. Moreover, these cases clearly illsutrate the factthat the path time is the most important performance measures in most practical applications.Special cases of the network are, for example, standard polling systems [ ] , tandem queues [ ] , multi-stage queueing models with parallel queues [ ] , feedback vacation queues [
7, 34 ] ,symmetric feedback polling systems [
32, 34 ] , systems with a waiting room [
1, 31 ] . We wouldlike to remind the reader that the only practical alternative to our expressions for a completeperformance analysis of the studied network is simulation. The first application is a multiproduct system with random yields introduced by [ ] , which is astylised practical case of a producer of plastic bumpers for automobiles. Every item is producedin a production queue and is defect with probability p ; a defect item has to be reproduced in thesame production queue. However, defect items are first routed to a temporary storage queue,which is served immediately after the current queue, before they are routed back to the originalqueue. The combination of exhaustively served production and storage queues implies that newlyarriving items are served during the current cycle, whereas defect items will be reproduced in thenext cycle. There are no switch-over time and service time required for these storage queues,implying that the storage queues do not contribute to the utilisation of the system.We consider a system consisting of five production queues with relative loads ( ) ,Poisson arrivals, exponentially distributed service times with mean 1, and switch-over times withconstant value 5. A typical feature of these type of production systems is that the switch-overtimes are relatively large compared to the processing times, in contrast to the next numericalexample that we discuss. We have additional five storage queues to which the defect items arerouted. Each item is defect with probability 0.25. We run simulations for various values of ρ , andwe compare the results with the limiting HT distributions. For each value of ρ we run 100 simu-lations, each of length 10 . These settings yield extremely accurate estimates for the probabilitydistribution functions. 15 onvergence to the HT limit. In this example, mostly due to the relatively large switch-overtimes, the (scaled) path-time distributions converge relatively fast to their HT limits. This isillustrated in Table 1 and in Figure 2. Table 1 depicts the means, standard deviations and tailprobabilities of the (scaled) path-time distributions for nine paths. Only production queues areincluded in the path names. For example, path 2 → → ρ < ρ = ρ = ρ → → → → → → → → → ρ → → → → → → → → → ρ → → → → → → → → → ρ
20 50 80 20 50 80 20 50 800.8 0.191 0.154 0.159 0.152 0.120 0.130 0.113 0.089 0.1030.9 0.187 0.161 0.175 0.144 0.123 0.142 0.103 0.090 0.1130.95 0.186 0.165 0.184 0.141 0.125 0.149 0.098 0.090 0.1181.00 0.184 0.168 0.192 0.138 0.126 0.156 0.093 0.091 0.123Table 1: Numerical results for Example 1.16 ccuracy of the approximation.
Table 2 shows the approximated means using the path-timeapproximation discussed in Section 4. For nine paths, the means of the (scaled) path-time ap-proximations are given. Only production queues are included in the path names. For example,path 2 → → scaled path times ( − ρ ) W i , i ,..., i M , becausenow they can be compared directly to the simulated values in Table 1. The results show that theaccuracy is very high for all values of the load and for all paths.Approximated scaled path times: Means ρ → → → → → → → → → Secondly, we use an example that was introduced by Katayama [ ] , who studies call processingin packet switching systems composed of multi-processors. In this application, each processorhas two kinds of buffers in parallel for receiving data packets from switching links and from linesand / or subscribers. During the first stage, called input processing, the processor polls both buffersfor data packets and moves these packets to a queue in the second stage. During this secondstage, called packet processing, the jobs are actually executed. This system can be modelled as anetwork consisting of three queues. Customers arrive at Q and Q , and are routed to Q afterbeing served (see Figure 3). This tandem queueing model with parallel queues in the first stageis a special case of a roving server network. We simply put p = p = p = p i , j are zero. We use the same values as in [ ] : Poisson arrivals with λ = λ /
10, service timesare deterministic with b = b =
1, and b =
5. The server serves the queues exhaustively, incyclic order: 1, 2, 3, 1, . . . . The only difference with the model discussed in [ ] is that weintroduce (deterministic) switch-over times r = r =
2. We assume that no time is required toswitch between the two queues in the first stage, so r = Convergence to the HT limit.
For the two possible paths, the means, standard deviations andtail probabilities of the (scaled) path-time distributions are given in Table 3 and in Figure 4. Thetail probabilities in Table 3 are given for path lengths of 20. The results for ρ < ρ = ccuracy of the approximation. Table 4 shows the approximated means using the path-timeapproximation discussed in Section 4 for both paths. Again the quality of the approximation isexcellent for all possible loads and all paths: The relative approximation error is less than 2% forall values of ρ .We have also tested the accuracy of the approximation for different interarrival-time distribu-tions, with squared coefficient of variation (SCV) equal to respectively and 2. In the first casewe have fitted a mixed Erlang distribution, and in the second case a hyperexponential distribu-tion. In these cases, the approximation still gives good results but slightly less accurate than forthe case with Poisson arrivals: when the SCV equals 2, the relative approximation error is lessthan 4% for all values of ρ . When the SCV is equal to , the relative error is less than 9%. Asexpected, the relative difference between the simulated and approximated values is biggest when ρ lies between 0.3 and 0.7. This application is taken from Sidi et al. [ ] , who consider a token ring network with onefile-server and K workstations, that transmit file requests to the file-server, which in turn repliesto the different stations by sending back the files they requested. The performance measure ofinterest is the response time of a file request, which is the time from the request generationby the workstation (being queued at that time at the output queue of the station, awaiting itsturn to be transmitted) until the file arrives back at the station. We can model this system as aroving server network with K + Q , . . . , Q K and are routedto Q K + after their service completion. Once served at Q K + the customer leaves the system.We are interested in the path time W i , K + for i =
1, . . . , K . In our numerical example we take N = K + =
11, identical arrival rates at the K workstations and no external arrivals at the fileserver, i.e. λ K + =
0. Since arrival processes in this application area tend to exhibit high variation,we assume that the interarrival times follow a hyperexponential distribution with balanced means(see, e.g., [ ] ) and squared coefficient of variation equal to 4. The routing probabilities are all0 except p i , K + = i =
1, . . . , K . The service times B , . . . , B K are all deterministic with value0.1, and the service times at Q K + are exponentially distributed with mean 1. The switch-overtimes in this kind of system are typically very small compared to the service times, so we takeexponentially distributed switch-over times with mean 0.01. The service discipline is gated at allqueues.We have chosen to highlight only a few typical results for this example, in Table 5 and Figure5, where the means and standard deviations of the scaled path times are shown. The values for ρ =
0, 9, 0.95, 0.98 have been obtained by simulation, whereas the values for ρ = ρ (see Figure5). Most interestingly, this is not caused by simulation inaccuracy, but by the fact that the pathtimes densities are (almost completely) concentrated on a discrete set of values for smaller loadsdue to the combination of small switch-over times, deterministic service times and deterministic routing. This effect disappears when the loads increase, but as a consequence the convergence tothe limiting HT distribution is not as fast as in the previous example. Nevertheless, despite theanomalous behaviour of the system’s performance for small loads, the derived asymptotic turnsout to be exact again. 18 .4 Example 4: A production system with a joint packaging queue and randomyield In the fourth example we look at a production system with two types of products and a jointpackaging queue and random yield. A set-up which is typically observed in many practical ap-plications. Product type A requires three operations, at queue 1, queue 3 and queue 5, beforejoining the packaging queue 7, whereas the other product type B visits queue 2, queue 4 andqueue 6 before it is sent to the joint packaging queue. Every production step fails with proba-bility p = A and queue 2 for product B . Packaging is, however, always successful. Weassume that all queues are served exhaustively.We study a system with external arrivals at queue 1 and queue 2 with intensity 0.12 and0.04, respectively. Furthermore we assume Erlang distributed arrivals with SCV 0.25, Erlangdistributed service times with means 1 (for the regular queues) and 2 (for the packaging queue)and SCV 0.5, and switch-over times with constant value 5.For two possible paths, the simulated scaled path-time distributions are depicted in Figure 6.The analytically computed limiting distribution at ρ = p . Therefore, we vary p while keeping all the other parameters fixed. Note thatincreasing p will also increase the total workload of the system. When considering the load ofthe system as a function of p , the stability condition can be rewritten to p < Q isextremely accurate. This is as expected, since p = ρ = ρ > E [ W ] simplifies to(5.1). As a consequence, the approximation, is very suitable for optimisation purposes. E [ W approx ] = € p − p + p − p + p − p + Š − × € p − p + p − p + p − p + p − p + p + p − p + Š . (5.1)We can therefore conclude that, although the four analyzed cases come from completely dif-ferent fields of application, have dissimilar parameter settings and differ significantly in system’sbehaviour, the derived HT asymptotics are shown to accurately describe the performance in allof these cases. Moreover, in all of these models the key performance metric is obviously the totaltime spent in the system by an arbitrary customer, i.e., the path time. We hope that this observa-tion illustrates the general nature of our framework, which extends and unifies the HT analysisof many queueing models. 19 Conclusions and suggestions for further research
In the current paper, we have analyzed a queueing network with customer routing, where a singleshared server serves the queues in a cyclic order. We have not only studied the waiting time ofa customer at a certain queue, but also the path time (the total time spent in the system byan arbitrary customer traversing a specific path). The main complicating factor in this analysisis the routing of customers, which leads to non-renewal arrival processes at the queues andto strong interdependence of the waiting times at the queues. These factors prohibit an exactexplicit analysis with closed-form expressions and, therefore, it is natural to resort to asymptoticestimates. That is, we have obtained easily computable expressions for both the waiting-time andpath-time distribution in HT. Combining these HT asymptotics with newly developed LT limitsleads to highly accurate approximations for the mean performance measures for the whole rangeof load values. The strength of this refined heavy-traffic approximation lies in its simplicity, whichopens up interesting possibilities for optimization of the system performance with respect to theroutes of the customers.
Acknowledgements
The authors are grateful to the anonymous referees, who have provided valuable comments re-sulting in improved readability of the current manuscript.
References [ ] O. M. E. Ali and M. F. Neuts. A service system with two stages of waiting and feedback ofcustomers.
Journal of Applied Probability , 21:404–413, 1984. [ ] K. B. Athreya and P. E. Ney.
Branching Processes . Springer-Verlag, Berlin, 1972. [ ] M. A. A. Boon, R. D. van der Mei, and E. M. M. Winands. Applications of polling systems.
Surveys in Operations Research and Management Science , 16:67–82, 2011. [ ] M. A. A. Boon, R. D. van der Mei, and E. M. M. Winands. Queueing networks with a singleshared server: light and heavy traffic.
SIGMETRICS Performance Evaluation Review , 39(2):44–46, 2011. [ ] M. A. A. Boon, E. M. M. Winands, I. J. B. F. Adan, and A. C. C. van Wijk. Closed-form waitingtime approximations for polling systems.
Performance Evaluation , 68:290–306, 2011. [ ] M. A. A. Boon, R. D. van der Mei, and E. M. M. Winands. Waiting times in queueingnetworks with a single shared server.
Queueing Systems , 74(4):403–429, 2013. [ ] O. J. Boxma and U. Yechiali. An M / G / Journal of Applied Probability , 34:773–784, 1997. [ ] E. G. Coffman, Jr., A. A. Puhalskii, and M. I. Reiman. Polling systems with zero switchovertimes: A heavy-traffic averaging principle.
The Annals of Applied Probability , 5(3):681–719,1995. 20 ] E. G. Coffman, Jr., A. A. Puhalskii, and M. I. Reiman. Polling systems in heavy-traffic: ABessel process limit.
Mathematics of Operations Research , 23:257–304, 1998. [ ] J. L. Dorsman, R. D. van der Mei, and E. M. M. Winands. A new method for derivingwaiting-time approximations in polling systems with renewal arrivals.
Stochastic Models , 27(2):318–332, 2011. [ ] Y. Gong and R. de Koster. A polling-based dynamic order picking system for online retailers.
IIE Transactions , 40:1070–1082, 2008. [ ] S. E. Grasman, T. L. Olsen, and J. R. Birge. Setting basestock levels in multiproduct systemswith setups and random yield.
IIE Transactions , 40(12):1158–1170, 2008. [ ] D. Grillo. Polling mechanism models in communication systems – some application ex-amples. In H. Takagi, editor,
Stochastic Analysis of Computer and Communication Systems ,pages 659–699. North-Holland, Amsterdam, 1990. [ ] O. B. Jennings. Averaging principles for a diffusion-scaled, heavy-traffic polling station with K job classes. Mathematics of Operations Research , 35(3):669–703, 2010. [ ] T. Katayama. A cyclic service tandem queueing model with parallel queues in the first stage.
Stochastic Models , 4:421–443, 1988. [ ] V. Kavitha and E. Altman. Queueing in space: design of message ferry routes in static adhocnetworks. In
Proceedings ITC21 , 2009. [ ] H. Levy and M. Sidi. Polling systems: applications, modeling, and optimization.
IEEETransactions on Communications , 38:1750–1760, 1990. [ ] D. Markowitz and L. Wein. Heavy traffic analysis of dynamic cyclic policies: A unifiedtreatment of the single machine scheduling problem.
Operations Research , 49(2):246–270,2001. [ ] D. Markowitz, M. Reiman, and L. Wein. The stochastic economic lot scheduling problem:Heavy traffic analysis of dynamic cyclic policies.
Operations Research , 48(2):136–154, 2000. [ ] S. S. Nair. A single server tandem queue.
Journal of Applied Probability , 8(1):95–109, 1971. [ ] T. L. Olsen and R. D. van der Mei. Polling systems with periodic server routeing in heavytraffic: distribution of the delay.
Journal of Applied Probability , 40:305–326, 2003. [ ] T. L. Olsen and R. D. van der Mei. Polling systems with periodic server routing in heavytraffic: renewal arrivals.
Operations Research Letters , 33:17–25, 2005. [ ] M. P. Quine. The multi-type Galton-Watson process with immigration.
Journal of AppliedProbability , 7(2):411–422, 1970. [ ] M. Reiman and L. Wein. Dynamic scheduling of a two-class queue with setups.
OperationsResearch , 46(4):532–547, 1998. [ ] M. Reiman and L. Wein. Heavy traffic analysis of polling systems in tandem.
OperationsResearch , 47(4):524–534, 1999. 21 ] M. Reiman, R. Rubio, and L. Wein. Heavy traffic analysis of the dynamic stochasticinventory-routing problem.
Transportation Science , 33(4):361–380, 1999. [ ] J. A. C. Resing. Polling systems and multitype branching processes.
Queueing Systems , 13:409–426, 1993. [ ] D. Sarkar and W. I. Zangwill. File and work transfers in cyclic queue systems.
ManagementScience , 38(10):1510–1523, 1992. [ ] M. Sidi and H. Levy. Customer routing in polling systems. In P. King, I. Mitrani, andR. Pooley, editors,
Proceedings Performance ’90 , pages 319–331. North-Holland, Amsterdam,1990. [ ] M. Sidi, H. Levy, and S. W. Fuhrmann. A queueing network with a single cyclically rovingserver.
Queueing Systems , 11:121–144, 1992. [ ] L. Takács. A queuing model with feedback.
Revue française d’automatique, d’informatique etde recherche opérationnelle. Recherche opérationnelle , 11(4):345–354, 1977. [ ] H. Takagi. Analysis and applications of a multiqueue cyclic service system with feedback.
IEEE Transactions on Communications - TCOM , 35(2):248–250, 1987. [ ] H. Takagi. Analysis and application of polling models. In G. Haring, C. Lindemann, andM. Reiser, editors,
Performance Evaluation: Origins and Directions , volume 1769 of
LectureNotes in Computer Science , pages 424–442. Springer Verlag, Berlin, 2000. [ ] T. Takine, H. Takagi, and T. Hasegawa. Sojourn times in vacation and polling systems withBernoulli feedback.
Journal of Applied Probability , 28(2):422–432, 1991. [ ] M. Taube-Netto. Two queues in tandem attended by a single server.
Operations Research ,25(1):140–147, 1977. [ ] H. C. Tijms.
Stochastic models: an algorithmic approach . Wiley, Chichester, 1994. [ ] R. D. van der Mei. Towards a unifying theory on branching-type polling models in heavytraffic.
Queueing Systems , 57:29–46, 2007. [ ] R. D. van der Mei and E. M. M. Winands. A note on polling models with renewal arrivalsand nonzero switch-over times.
Operations Research Letters , 36:500–505, 2008. [ ] W. Whitt.
Stochastic-Process Limits: An Introduction to Stochastic-Process Limits and TheirApplication to Queues . Springer Series in Operations Research and Financial Engineering.Springer, Berlin, 2002. [ ] E. M. M. Winands. On polling systems with large setups.
Operations Research Letters , 35:584–590, 2007. [ ] E. M. M. Winands. Branching-type polling systems with large setups.
OR Spectrum , 33(1):77–97, 2011. 22 ppendixA Proof for Poisson arrivals
In this appendix we present a proof for Conjecture 3.4 for the special case of Poisson arrivals,reformulated in the theorem below. We conclude the appendix with a short discussion.
Theorem A.1
Assume that external customers arrive at Q i according to a Poisson process with rate λ i . As ρ ↑ , the scaled queue length ( − ρ ) L i converges in distribution to the product of twoindependent random variables. The first has the same distribution as L fluidi , and the second randomvariable Γ has the same distribution as the limiting distribution of the scaled length-biased cycle time, ( − ρ ) C i . For i =
1, . . . , N ; k = i , . . . , i + N − , and ρ ↑ , ( − ρ ) L i d → Γ × L fluidi , k w.p. ˆ ρ k , (A.1) where Γ is a random variable having a Gamma distribution with parameters α + and δµ , and L fluidi , k is the “standardised” number of a type-i particles in the fluid model during V k , introducedbelow (3.4) , independent of Γ . The parameters of the Gamma distribution, in the case of Poissonarrivals, are defined as α = r δ/ ˜ b res , µ = ˜ b res , with δ as defined in Definition 3.1. The proof of Theorem A.1 is structured in the form of 9 small steps. Steps 1 – 7 involve findingthe limiting scaled joint queue-length distribution at visit beginnings and completions. We relyheavily on existing limiting results for polling systems and branching processes, exploiting thesimilarities between these models and roving server networks, closely following the frameworkdeveloped in [ ] . The remaining steps follow the approach in [ ] and [ ] , expressing queuelengths at arbitrary moments in terms of the queue-length distributions at visit beginnings, inheavy traffic. We start by giving some preliminary results, required for the remainder of theproof. Preliminaries.
As shown in [
6, 30 ] the joint queue-length process at visit beginnings constitutesan N -dimensional multi-type branching process (MTBP) with immigration. As preliminaries, wegive a brief overview of existing limiting results on critical MTBPs, required for the remainderof this proof. We start by introducing some additional notation. An N -dimensional vector v hascomponents ( v , . . . , v N ) , and we denote | v | : = P Ni = v i . Finally, I E is the indicator function on theevent E .Let Z = { Z n , n =
0, 1, . . . } be an N -dimensional multi-type branching process, where Z n =( Z ( ) n , . . . , Z ( N ) n ) is an N -dimensional vector denoting the state of the process in the n -th gen-eration. The MTBP is defined by its offspring generating functions and its immigration func-tions. The one-step offspring generating function is denoted by f ( z ) = ( f ( ) ( z ) , . . . , f ( N ) ( z )) ,with z = ( z , . . . , z N ) , and f ( i ) ( z ) = X j ,..., j N ≥ p ( i ) ( j , . . . , j N ) z j · · · z j N N , (A.2)23here p ( i ) ( j , . . . , j N ) is the probability that a type- i particle produces j k particles of type k ( i , k =
1, . . . , N ) . The immigration function is denoted as follows: g ( z ) = X j ,..., j N ≥ q ( j , . . . , j N ) z j · · · z j N N , (A.3)where q ( j , . . . , j N ) is the probability that a group of immigrant consists of j k particles of type k ( i , k =
1, . . . , N ) . Denote the mean immigration vector by g : = ( g , . . . , g N ) , where g i : = ∂ g ( z ) ∂ z i | z = ( i =
1, . . . , N ) , (A.4)and the mean offspring matrix , or simply mean matrix , by M = € m i , j Š , with m i , j : = ∂ f ( i ) ( z ) ∂ z j | z = ( i , j =
1, . . . , N ) . (A.5)Thus, for a given type- i particle, m i , j is the mean number of type- j “children” it has in the nextgeneration. Similarly, for a type- i particle, the second-order derivatives are denoted by the matrix K ( i ) = (cid:16) k ( i ) j , k (cid:17) , with k ( i ) j , k : = ∂ f ( i ) ( z ) ∂ z j ∂ z k | z = , ( i , j , k =
1, . . . , N ) . (A.6)Denote by v = ( v , . . . , v N ) and w = ( w , . . . , w N ) the left and right eigenvectors correspondingto the largest real-valued, positive eigenvalue ξ of M , commonly referred to as the maximumeigenvalue (cf., e.g., [ ] ), normalized such that v ⊤ = v ⊤ w =
1. (A.7)The following conditions are necessary and sufficient conditions for the ergodicity of the process Z (cf. [ ] ): ξ < X j + ··· + j N > q ( j , . . . , j N ) log ( j + · · · + j N ) < ∞ . (A.8)Note that ξ plays a role similar to ρ in the roving server network. The hat-notation introduced inSection 2 is also used here to indicate that x is evaluated at ξ =
1. Moreover, for ξ ≥ π ( ξ ) : =
0, and π n ( ξ ) : = n X r = ξ r − , n =
1, 2, . . . . (A.9)Quine [
23, Theorem 4 ] derives the following property for critical MTBPs. Property A.2
Assume that all derivatives of f ( z ) through order two exist at z = and that < g i < ∞ ( i =
1, . . . , N ) . Then π n ( ξ ) Z ( ) n ...Z ( N ) n d → A ˆ v ... ˆ v N Γ( α , 1 ) as ( ξ , n ) → ( ∞ ) , (A.10) where ˆ v = (ˆ v , . . . , ˆ v N ) is the normalized the left eigenvector of ˆ M , and where Γ( α , 1 ) is a gamma-distributed random variable with scale parameter 1 and shape parameter α : = A ˆ g ⊤ ˆ w = A N X i = ˆ g i ˆ w i , wi th A : = N X i = ˆ v i € ˆ w ⊤ ˆ K ( i ) ˆ w Š >
0. (A.11)24 tep 1: Characterise the embedded MTBP.
We now return to the roving server network. Inthis step we show how the joint queue-length process at successive polling instants at Q can bedescribed as an MTBP with immigration in each state. To this end, let X : = € X ,1 . . . , X N Š be the N -dimensional vector that describes the joint queue length at an arbitrary polling instant at Q .Let X i , n be the number of type- i customers in the system at the n -th polling instant of the serverat Q , for i =
1, . . . , N and n =
0, 1, . . ., and let X n : = € X n , . . . , X N , n Š (A.12)be the joint queue-length vector at the n -th polling instant at Q . In this appendix we will providethe results for gated service only. The results for exhaustive service can be obtained along theexact same lines, but require significantly more complex notations (see, for example, [ ] ) andare omitted here.The following result describes the process { X n , n =
0, 1, . . . } as an MTBP. Theorem A.3
The discrete-time process { X n , n =
0, 1, . . . } constitutes a N -dimensional MTBP withimmigration in each state. Denote by B ∗ i ( s ) and R ∗ i ( s ) the Laplace-Stieltjes transforms (LSTs) ofrespectively the service times B i and the switch-over times R i . The probability generating function(PGF) of the offspring function is given by the following expression: For | z i | ≤ ( i =
1, . . . , N ) ,f ( z ) : = € f ( ) ( z ) , . . . , f ( N ) ( z ) Š , (A.13) where for i =
1, . . . ,
N the function f ( i ) ( z ) is the unique solution off ( i ) ( z ) : = B ∗ i (cid:0) N X j = λ j € − z j Š (cid:1) · N X j = p i , j f ( j ) ( z ) , (A.14) and where the PGF of the immigration function is given byg ( z ) : = N Y i = R ∗ i (cid:16) i X j = λ j € − z j Š + N X j = i + λ j ( − f ( j ) ( z )) (cid:17) . (A.15) Proof:
Relations (A.13)-(A.15) can be obtained along the lines of Resing [ ] for the case ofclassical gated service, using simple generating-function manipulations. The only difference isthat the offspring function f ( i ) ( z ) has changed in the sense that each customer served at Q i iseffectively replaced not only by all customers that arrive at the different queues during its servicestime (leading to PGF B ∗ i ( P Nj = λ j ( − z j )) in Equation (A.14)), but in addition by a fresh customerat Q j (which creates an additional offspring f ( j ) ( z ) ), with probability p i , j . Next, Equation (A.15)stems from the fact that the immigration consists of the contributions of newly arriving customersthat arrive during the switch-over times R i , i =
1, . . . , N . ƒ Step 2: Compute the mean offspring matrix for the roving server network.
The followingresult gives a characterization of the mean offspring matrix M defined in (A.5). Lemma A.4
The mean offspring matrix M = € m i , j Š can be expressed by M = M · · · M N , (A.16)25 here for k =
1, . . . ,
N , the elements of the matrix M k = (cid:16) m ( k ) i , j (cid:17) are given by: For i , j =
1, . . . ,
N ,i = k, m ( k ) i , j = I { i = j } , (A.17) and for i = k, m ( k ) i , j = λ j b i + p i , j for j =
1, . . . , N . (A.18) Proof:
The result can be obtained directly from Theorem A.3 by taking the partial derivatives ofthe offspring function defined in (A.13) and (A.14). ƒ Step 3: Determine the left and right eigenvectors of M at ρ = The following result givesthe left and right eigenvectors (normalized according to (A.7)) of the mean offspring matrix M ,defined in (A.16)-(A.18), evaluated at ρ = Lemma A.5
The right eigenvector ˆ w of the mean matrix ˆ M , normalized such that ˆ w ⊤ ˆ w = ,corresponding with maximum eigenvalue 1, is given by ˆ w = ˆ w ... ˆ w N : = | ˜ b | − ˜ b , wi th ˜ b : = ˜ b ... ˜ b N . (A.19) The corresponding left eigenvector ˆ v, normalized such that ˆ v ⊤ ˆ w = , corresponding with maximumeigenvalue 1, is given by ˆ v = ˆ v ... ˆ v N : = | ˜ b | δ ˆ u , wi th u : = u ...u N , wher e u i : = λ i N X j = i ρ j + N X j = i γ j p j , i , (A.20) and where δ : = ˆ u ⊤ ˜ b = N X i = N X j = i + ˆ ρ i ˆ ρ j + N X i = ˜ b i N X j = i ˆ γ j p j , i . (A.21) Proof:
First, it is readily seen by using equations (A.17)-(A.18) and (2.1) that for k =
1, . . . , N ,we have N X j = ˆ m ( k ) k , j ˜ b j = N X j = € ˆ λ j b k + p k , j Š ˜ b j = b k N X j = ˆ λ j ˜ b j + N X j = p k , j ˜ b j = b k + N X j = p k , j ˜ b j = ˜ b k . (A.22)This immediately implies that ˆ M k ˆ w = ˆ w for k =
1, . . . , N , and hence from (A.16) that ˆ M ˆ w =ˆ M · · · ˆ M N ˆ w = ˆ w , which shows that ˆ w indeed is a right eigenvector of ˆ M with eigenvalue 1.Similar arguments can be used to show that ˆ M ⊤ ˆ v = ˆ v . ƒ Remark A.6
Although δ defined in (A.21) may appear different, at first sight, than δ defined inDefinition 3.1, it can be shown that they are in fact identical. tep 4: Mean immigration function at ρ = We now proceed to specify the mean immigra-tion vector g , defined in (A.7), for the model under consideration. Considering the evolution ofthe N -dimensional state vector as a discrete-time Markov chain { X n , n =
0, 1, . . . } at successivepolling instants at Q , the “immigrants” in the n -th generation are the customers present a time n that are not children of any of the customers present at time n −
1. Denote the mean immigrationvector by g = ( g , . . . , g N ) , where g i stands for the mean number of type- i immigrants. Lemma A.7
The mean immigration function is given by g = ( g , . . . , g N ) , where for j =
1, . . . ,
N ,g j = N X i = r i λ j I { j ≤ i } + N X k = i + λ k m k , j ! . (A.23) Moreover, ˆ g ⊤ ˆ w = | b | − r . (A.24) Proof:
Equation (A.23) follows directly from Theorem A.3 by differentiating once with respectto s j and subsituting s = (
1, . . . , 1 ) . Next, to prove (A.24), assume ρ =
1. We first observe that itfollows from (A.23) that the mean number of type- j customers that immigrate during a cycle isgiven by ˆ g j = N X i = r i ˆ λ j I { j ≤ i } + N X k = i + ˆ λ k ˆ m k , j ! . (A.25)This implies ˆ g ⊤ ˆ w : = | ˜ b | − N X j = ˆ g j ˜ b j = | ˜ b | − N X i = j r i i X j = ˜ b j ˆ λ j + N X k = i + ˆ λ k N X j = ˆ m k , j ˜ b j (A.26) = | ˜ b | − N X i = r i i X j = ˜ b j ˆ λ j + N X k = i + ˆ λ k ˜ b k = | ˜ b | − r N X i = ˆ ρ i = | ˜ b | − r , (A.27)by using the definition in (A.23) and the results in Lemma 2. ƒ Step 5: Determine the parameter A . The following result gives an expression for the scalingparameter A , defined in (A.11). Lemma A.8 A = | ˜ b | − δ − · ˜ b res , (A.28)where ˜ b res = E [ e B ] E [ e B ] denotes the expected residual extended service time e B of an arbitrary customerin the system. More precisely, E [ e B k ] = N X i = λ i E [ e B ki ] / N X j = λ j . Proof:
This result can be proven following the same lines as in [ ] (Remark 5) for the case ofbranching-type service policies for the case without customer routing (i.e., p i , j = i , j ). ƒ tep 6: Asymptotic properties for the maximum eigenvalue of M . The following resultdescribes the limiting behaviour of the maximum eigenvalue ξ ( ρ ) of the matrix M defined inLemma 1, considered as a function of ρ , as ρ goes to 1. Lemma A.9
The maximum eigenvalue ξ = ξ ( ρ ) satisfies the following properties: (1) ξ < if and only if ρ < , ξ = if and only if ρ = and ξ > if and only if ρ > ; (2) ξ ( ρ ) is a continuous function of ρ ; (3) lim ρ ↑ ξ ( ρ ) = ξ ( ) = ; (4) the derivative of ξ ( ρ ) at ρ = is given by ξ ′ ( ) = lim ρ ↑ − ξ ( ρ ) − ρ = δ , (A.29) where δ is defined in (A.21). Proof:
The proof can be obtained by following the approach in [ ] (Section 3.2), which provesfor the case without customer routing (i.e., p i , j = i , j ). ƒ Step 7: The joint queue-length vector at visit beginnings and completions.
We are nowready to present the HT result for the state vector at polling instants. Without loss of generality,we focus on the evolution of the state vector at embedded polling instants at Q . Theorem A.10
The joint queue-length vector at polling instants at Q has the following asymptoticbehaviour: ( − ρ ) X ...X N d → ˜ b res δ ˆ u ... ˆ u N Γ( α , 1 ) ( ρ ↑ ) , (A.30) where α = r δ/ ˜ b res , (A.31) and where ˆ u i ( i =
1, . . . , N ) and δ are defined in (A.20) and (A.21). Proof:
First, note that the process that describes the evolution of the state vector { X n , n =
0, 1, . . . } at successive polling instants at Q constitutes an N -dimensional MTBP with offspring function f ( s ) and immigration function g ( z ) defined in Theorem A.3, and with mean matrix M definedin Lemma A.4. Moreover, from Theorem A.3 it is readily verified that the assumptions of Property1 on the finiteness of the second-order derivatives of f ( z ) and the mean immigration function g are satisfied. Then using Property A.2 it follows that1 π n ( ξ ) · X ⊤ n d → A · ˆ v · Γ( α , 1 ) as ( ξ , n ) → ( ∞ ) , (A.32)where A , ˆ v and α are given in (A.11). Hence, translating this to the polling model and usingLemmas 1-5, it readily follows from (A.32) that ( − ρ ) X ⊤ n d → δ · A · ˆ v · Γ( α , 1 ) as ( ρ , n ) → ( ∞ ) , (A.33)where expressions for δ , A , ˆ v and α are given in (A.21), (A.28), (A.20) and (A.31). Combiningthese expressions leads to the result. ƒ orollary A.11 Let X scaledi , k : = lim ρ ↑ ( − ρ ) X i , k denote the scaled number of customers in Q i at thebeginning of a visit to Q k . Its LST is equal to E h e − ω X scaledi , k i = (cid:18) δ/ ˜ b res δ/ ˜ b res + ω P k − j = i ˆ ρ j ˆ γ i , j (cid:19) r δ/ ˜ b res k = i +
1, . . . , i + N − (cid:16) δ/ ˜ b res δ/ ˜ b res + ω ˆ γ i (cid:17) r δ/ ˜ b res k = i . (A.34) Step 8: The marginal queue-length distribution.
Let L scaledi , k : = lim ρ ↑ ( − ρ ) L i , k denote thescaled number of customers in Q i at an arbitrary moment during a visit to Q k , and let L scaledi : = lim ρ ↑ ( − ρ ) L i denote the scaled number of customers in Q i at an arbitrary moment. Theorem A.12 E h e − ω L scaledi i = N X k = ˆ ρ k E h e − ω L scaledi , k i , (A.35) where E h e − ω L scaledi , i i = (cid:129) δ/ ˜ b res δ/ ˜ b res + ω ˆ ρ i ˆ γ i , i ‹ r δ/ ˜ b res − (cid:16) δ/ ˜ b res δ/ ˜ b res + ω ˆ γ i (cid:17) r δ/ ˜ b res (ˆ γ i − ˆ γ i , i ˆ ρ i ) r ω , (A.36) E h e − ω L scaledi , k i = (cid:18) δ/ ˜ b res δ/ ˜ b res + ω P k − j = i ˆ ρ j ˆ γ i , j (cid:19) r δ/ ˜ b res − (cid:18) δ/ ˜ b res δ/ ˜ b res + ω P kj = i ˆ ρ j ˆ γ i , j (cid:19) r δ/ ˜ b res ˆ ρ k ˆ γ i , k r ω , i = k . (A.37) Proof:
To prove Theorem A.12 we need the following results, summarised in (A.38)–(A.41),which directly follow from Equations (3.6)–(3.10) in [ ] . E [ z L i ] = N X k = (cid:129) ˆ ρ k E [ z L ( Vk ) i ] + ( − ρ ) r j r E [ z L ( Rk ) i ] ‹ , (A.38)where L ( V k ) i and L ( R k ) i denote the number of customers in Q i at an arbitrary moment during V k and R k respectively.Denote by X ∗ i , k ( z ) and Y ∗ i , k ( z ) the PGF of the number of customers in Q i at the beginning andend of V k , respectively. L ( V i ) i ( z ) = (cid:0) X ∗ i , i ( z ) − Y ∗ i , i ( z ) (cid:1) ( − ρ )( − B ∗ i (cid:0) λ i ( − z ) (cid:1) ) λ i ( − z ) ρ i r (cid:0) − B ∗ i ( λ i ( − z ))( − p i , i + p i , i z ) / z ) (A.39) L ( V k ) i ( z ) = (cid:0) X ∗ i , k ( z ) − Y ∗ i , k ( z ) (cid:1) ( − ρ )( − B ∗ k (cid:0) λ i ( − z ) (cid:1) ) λ i ( − z ) ρ k r (cid:0) − B ∗ k ( λ i ( − z ))( − p k , i + p k , i z )) , for k = i .(A.40)Note that Y ∗ i , k ( z ) = X ∗ i , k + ( z ) R ∗ k ( λ i ( − z )) . (A.41)Our goal is to find the limiting distribution of ( − ρ ) L i as ρ ↑
1. In the limit, when substituting z = e − ω ( − ρ ) in (A.38) and letting ρ ↑
1, the terms that correspond to the queue lengths during29witch-over times vanish, caused by the ( − ρ ) term. Intuitively this is exactly what one wouldexpect, as switch-over times become negligible in heavy traffic. In order to prove (A.36) and(A.37), we substitute z = e − ω ( − ρ ) , λ i = ˆ λ i ρ , and ρ i = ˆ ρ i ρ in (A.39) and (A.40) respectively,and evaluate the Taylor series of the resulting functions near ρ = k = i this results inlim ρ ↑ E [ e − ω ( − ρ ) L ( Vi ) i ] = E h e − ω X scaledi , i + i − E h e − ω X scaledi , i i ω r ˆ ρ i ( − ˆ λ i b i − p i , i ) / b i , (A.42)and for the case k = i we obtainlim ρ ↑ E [ e − ω ( − ρ ) L ( Vk ) i ] = E h e − ω X scaledi , k i − E h e − ω X scaledi , k + i ω r ˆ ρ k (ˆ λ i b i + p k , i ) / b k . (A.43)After substitution of (A.34) this leads to the following result:lim ρ ↑ E [ e − ω ( − ρ ) L ( Vi ) i ] = (cid:129) δ/ ˜ b res δ/ ˜ b res + ˆ ρ i ˆ γ i , i ω ‹ r δ/ ˜ b res − (cid:16) δ/ ˜ b res δ/ ˜ b res +ˆ γ i ω (cid:17) r δ/ ˜ b res r (ˆ γ i − ˆ ρ i ˆ γ i , i ) ω , (A.44)lim ρ ↑ E [ e − ω ( − ρ ) L ( Vk ) i ] = (cid:18) δ/ ˜ b res δ/ ˜ b res + ω P k − j = i − N ˆ ρ j ˆ γ i , j (cid:19) r δ/ ˜ b res − (cid:18) δ/ ˜ b res δ/ ˜ b res + ω P kj = i − N ˆ ρ j ˆ γ i , j (cid:19) r δ/ ˜ b res r ˆ ρ k ˆ γ i , k ω . (A.45) Step 9: Proving Theorem A.1.
To prove that Theorem A.1 agrees with the results obtained inStep 8, we need the following, well-known property that is frequently used in HT limits of pollingsystems.
Property A.13
Let Γ be a random variable with a Gamma distribution with shape parameter α ∗ (to avoid confusion since we already have defined a parameter α ) and rate parameter β . Let U be auniform random variable on the interval [ a , b ] , independent of Γ . The LST of the product U × I isequal to E [ e − ω U I ] = (cid:16) β/ a β/ a + ω (cid:17) α ∗ − − (cid:16) β/ b β/ b + ω (cid:17) α ∗ − ( b − a )( α ∗ − ) ω/β . (A.46)Theorem A.1 states that the scaled queue length L scaledi is distributed as L scaledi , k with probability ˆ ρ k . This is in accordance with what one would obtain when substituting z = e − ω ( − ρ ) in (A.38)and letting ρ ↑ L scaledi , k is the product of two random variables. The first is a Gamma( α + δµ ) distribution, where α + δµ are the parameters of the Gamma distribution of the scaled,length-biased cycle time, discussed in Theorem A.1. When the external arrival process is a Poissonprocess, the parameter σ is equal to ˜ b ( ) / ˜ b , which means that α = r δ/ ˜ b res , µ = ˜ b res , (A.47)with δ as defined in Definition 3.1.The second random variable in the product is L fluidi , k , which is uniformly distributed. Remem-ber that L fluidi , k = L fluidi , k / c , so taking c = L fluidi , k . Substituting the following values in (A.46) leads to the LST of the limitingdistribution of the scaled number of customers in Q i at an arbitrary epoch during V k , α ∗ = r δ/ ˜ b res + β = δ ˜ b res , a = ˆ ρ j ˆ γ i , j , b = ˆ γ i , ( i = k ) a = k − X j = i ˆ ρ j ˆ γ i , j , b = k X j = i ˆ ρ j ˆ γ i , j , ( i = k ) It is quickly verified that these substitutions result in an expression completely equivalent to(A.44) for the case i = k , and to (A.45) for i = k . This concludes the proof of Theorem A.1. ƒ Discussion
In this appendix we have provided a proof for the scaled queue-length distributions in heavytraffic, relying on the framework developed in [ ] . In his paper, Van der Mei uses the distribu-tional form of Little’s Law to obtain the distributions of the scaled waiting times in polling modelswith Poisson arrivals. We note that the distributional form of Little’s Law cannot be applied toour model because of the internal routing, as discussed extensively in [ ] . In [ ] neverthelessa mathematical framework has been developed to derive the LSTs of the steady-state waiting-time distributions. As such, this framework provides an excellent basis to prove the HT limits forwaiting times without resorting to the distributional form of Little’s Law. The derived LSTs forsteady-state waiting-time distributions are given in the form of recursive expressions, which willsimplify to the elegant closed-form expressions presented in this paper, after taking the HT limit.For path times , there are currently no steady-state results available at all, due to the complexdependencies between successive visit times. In order to prove the HT limit of the scaled pathtimes, one would first have to apply the techniques in [ ] to find path-time LSTs in steady state,which is a separate study by itself.We conclude this discussion by noting that the mean waiting times can easily be obtained byapplying Little’s Law, which does not require the assumption of Poisson arrivals.31 Ρ =
Ρ =
Ρ =
Ρ =
20 40 60 800.0000.0050.0100.0150.0200.0250.0300.035
Ρ =
Ρ =
Ρ =
Ρ = →
20 40 60 80 1000.0000.0050.0100.0150.0200.025
Ρ =
Ρ =
Ρ =
Ρ = → → ρ < ρ = Sfrag replacements Q Q Q λ λ Server
Figure 3: Tandem queues with parallel queues in the first stage, as discussed in Example 2.Scaled path times: Means ρ → → ρ → → ρ → → = Ρ =
Ρ =
Ρ =
110 20 30 400.000.050.100.15 → Ρ =
Ρ =
Ρ =
Ρ =
110 20 30 400.000.020.040.060.080.100.12 → ρ < ρ = ρ → → ρ →
11 2 →
11 3 →
11 4 →
11 5 →
11 6 →
11 7 →
11 8 →
11 9 →
11 10 → ρ →
11 2 →
11 3 →
11 4 →
11 5 →
11 6 →
11 7 →
11 8 →
11 9 →
11 10 → = Ρ =
Ρ =
Ρ = → Ρ =
Ρ =
Ρ =
Ρ = → →
11 and 10 → = r= r= r= r= r=
120 40 60 80 1000.0000.0050.0100.0150.0200.0250.0300.035 → → → r= r= r= r= r= r=
120 40 60 80 1000.0000.0050.0100.0150.0200.0250.030 → → → → → → → → → Simulation Approximation0.02 0.04 0.06 0.08 0.10 0.12 0.14 p0100200300400500 E @ W D Figure 7: The mean waiting time in the first queue as a function of the rework probability pp