Two queues with random time-limited polling

TL;DR

This paper analyzes a two-queue polling system with a single server that spends random, exponentially distributed times at each queue, deriving workload distributions and queue length results, including heavy traffic and tail asymptotics.

Contribution

It introduces a novel analysis of a polling model with random, time-limited server visits, providing explicit workload and queue length distributions.

Findings

Derived steady-state workload distribution.

Obtained heavy traffic and tail asymptotics.

Calculated joint queue length distribution for exponential service times.

Abstract

In this paper, we analyse a single server polling model with two queues. Customers arrive at the two queues according to two independent Poisson processes. There is a single server that serves both queues with generally distributed service times. The server spends an exponentially distributed amount of time in each queue. After the completion of this residing time, the server instantaneously switches to the other queue, i.e., there is no switch-over time. For this polling model we derive the steady-state marginal workload distribution, as well as heavy traffic and heavy tail asymptotic results. Furthermore, we also calculate the joint queue length distribution for the special case of exponentially distributed service times using singular perturbation analysis.

Equations

$${\rho }_1<\frac{c_2}{c_1+c_2}\text{leavevmode}\text{ and }\text{leavevmode}{\rho }_2<\frac{c_1}{c_1+c_2},$$

$$\mathbb{E}(e^{-s{V}_1})=\frac{s\left[{\lambda }_1\mathbb{E}(B_1)\left(c_1+c_2\right)-c_2\right]\left[c_1+c_2+{\lambda }_1\left(1-{\tilde{b}}_1(s)\right)\right]}{\left[\left(c_2+{\lambda }_1\left(1-{\tilde{b}}_1(s)\right)\right)\left(c_1+{\lambda }_1\left(1-{\tilde{b}}_1(s)\right)-s\right)-c_1{c}_2\right]\left(c_1+c_2\right)}.$$

$$\begin{array}{rl}\mathbb{E}(e^{-s{V}_1(T_1+T_2)}|V_1(T_1)=y)& =e^{-sy}{\tilde{f}}_{T_2}({\lambda }_1(1-{\tilde{b}}_1(s)))\\ & =e^{-sy}\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}.\end{array}$$

$$\mathbb{E}(e^{-s{V}_1(T_1)}|V_1(0)=v)=$$

$${\int }_{\sigma =0}^{\infty }$$

$$-s{U}_1(t-v){\int }_{u=0}^{t-v}{e}^{(s-{\lambda }_1(1-{\tilde{b}}_1(s)))(t-u-v)}\mathbb{P}(V_1(u+v)=0|V_1(0)=v)\mathrm{d}u,$$

$$\begin{array}{rl}& \mathbb{E}(e^{-s{V}_1(T_1)}|V_1(0)=v)\\ & =\frac{c_1{e}^{-sv}}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}\\ & -{\int }_{t=v}^{\infty }s{c}_1{e}^{-c_{1}t}{\int }_{u=0}^{t-v}{e}^{(s-{\lambda }_1(1-{\tilde{b}}_1(s)))(t-u-v)}\mathbb{P}(V_1(u+v)=0|V_1(0)=v)\mathrm{d}u\mathrm{d}t.\end{array}$$

$${\int }_{t=0}^{\infty }{e}^{-st}\mathbb{P}(V_1(t)=0|V_1(0)=v)\mathrm{d}t=\frac{e^{-(s+(1-\mu (s,1)){\lambda }_1)v}}{s+(1-\mu (s,1)){\lambda }_1},$$

$${\int }_{t=v}^{\infty }s{c}_1{e}^{-c_{1}t}{\int }_{u=0}^{t-v}$$

$$=\frac{s{c}_1}{c_1-s+{\lambda }_1(1-{\tilde{b}}_1(s))}\frac{e^{-(c_1+(1-\mu (c_1,1)){\lambda }_1)v}}{c_1+(1-\mu (c_1,1)){\lambda }_1}.$$

$$\begin{array}{rl}& \mathbb{E}(e^{-s{V}_1(T_1)}|V_1(0)=v)\\ & =\frac{c_1{e}^{-sv}}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}-\frac{s{c}_1}{c_1-s+{\lambda }_1(1-{\tilde{b}}_1(s))}\frac{e^{-(c_1+(1-\mu (c_1,1)){\lambda }_1)v}}{c_1+(1-\mu (c_1,1)){\lambda }_1}.\end{array}$$

$$\begin{array}{rl}& \mathbb{E}(e^{-s{V}_1(T_1+T_2)}|V_1(0)=v)\\ & ={\int }_{y=0}^{\infty }\mathbb{E}(e^{-s{V}_1(T_1+T_2)}|V_1(T_1)=y)f_{V_1}(V_1(T_1)=y|V_1(0)=v)\mathrm{d}y\\ & =\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}{\int }_{y=0}^{\infty }{e}^{-sy}{f}_{V_1}(V_1(T_1)=y|V_1(0)=v)\mathrm{d}y\\ & =\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}[\frac{e^{-sv}{c}_1}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}\\ & -\frac{c_1{e}^{-sv}}{c_1-s+{\lambda }_1(1-{\tilde{b}}_1(s))}\frac{e^{-(c_1+(1-\mu (c_1,1)){\lambda }_1)v}}{c_1+(1-\mu (c_1,1)){\lambda }_1}],\end{array}$$

$$\begin{array}{rl}& \mathbb{E}(e^{-s{V}_1(T_1+T_2)})\\ & {\int }_{v=0}^{\infty }\mathbb{E}(e^{-s{V}_1(T_1+T_2)}|V_1(0)=v)f_{V_1(0)}(v)\mathrm{d}v\\ & ={\int }_{v=0}^{\infty }[\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}[e^{-sv}\frac{c_1}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}\\ & -s\frac{c_1}{c_1-s+{\lambda }_1(1-{\tilde{b}}_1(s))}\frac{e^{-(c_1+(1-\mu (c_1,1)){\lambda }_1)v}}{c_1+(1-\mu (c_1,1)){\lambda }_1}]]f_{V_1(0)}(v)\mathrm{d}v,\end{array}$$

$$\mathbb{E}(e^{-s{V}_1(T_1+T_2)})$$

$$={\int }_{v=0}^{\infty }[\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}[e^{-sv}\frac{c_1}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}$$

$$-s\frac{c_1}{c_1-s+{\lambda }_1(1-{\tilde{b}}_1(s))}\frac{e^{-(c_1+(1-\mu (c_1,1)){\lambda }_1)v}}{c_1+(1-\mu (c_1,1)){\lambda }_1}]]f_{V_1(T_1+T_2)}(v)\mathrm{d}v.$$

$$\begin{array}{rl}\mathbb{E}(& e^{-s{V}_1(T_1+T_2)})\left[\frac{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}{c_2}-\frac{c_1}{c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s}\right]\\ =& -\frac{s{c}_1}{\left[c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s\right]\left(c_1+(1-\mu (c_1,1)){\lambda }_1\right)}\\ & \mathbb{E}(e^{-\left(c_1+\left(1-\mu (c_1,1)\right){\lambda }_1\right)V_1(T_1+T_2)}).\end{array}$$

$$\mathbb{E}(e^{-\left(c_1+\left(1-\mu (c_1,1)\right){\lambda }_1\right)V_1(T_1+T_2)})$$

$$=-\frac{\left[{\lambda }_1\mathbb{E}(B_1)c_1+{\lambda }_1\mathbb{E}(B_1)c_2-c_2\right]\left(c_1+(1-\mu (c_1,1)){\lambda }_1\right)}{c_1{c}_2}.$$

$$\mathbb{E}(e^{-s{V}_1(T_1+T_2)})=\frac{s\left[{\lambda }_1\mathbb{E}(B_1)(c_1+c_2)-c_2\right]}{\left[c_2+{\lambda }_1(1-{\tilde{b}}_1(s))\right]\left[c_1+{\lambda }_1(1-{\tilde{b}}_1(s))-s\right]-c_1{c}_2}.$$

$$\mathbb{E}(e^{-s{V}_1})$$

$$=\mathbb{E}(e^{-s{V}_1}|S=1)\frac{c_2}{c_1+c_2}+\mathbb{E}(e^{-s{V}_1}|S=2)\frac{c_1}{c_1+c_2}.$$

$$\mathbb{E}(e^{-s{V}_1}|S=1)=\mathbb{E}(e^{-s{V}_1(T_1)}),\mathbb{E}(e^{-s{V}_1}|S=2)=\mathbb{E}(e^{-s{V}_1(T_1+T_2)}).$$

$$\mathbb{E}(e^{-s{V}_1(T_1+T_2)})=\mathbb{E}(e^{-s{V}_1(T_1)})\frac{c_2}{c_2+{\lambda }_1(1-{\tilde{b}}_1(s))}.$$

$$V_1\stackrel{\text{d}}{=}{V}_{\text{M/G/1}}+Y,$$

$$\begin{array}{rl}& \mathbb{E}(e^{-sY})\\ & =\frac{c_2-{\rho }_1(c_1+c_2)}{(1-{\rho }_1)(c_1+c_2)}\left[1-\frac{s{c}_1}{(c_2+c_1-s){\lambda }_1(1-{\tilde{b}}_1(s))+{\left({\lambda }_1(1-{\tilde{b}}_1(s))\right)}^2-s{c}_2}\right].\end{array}$$

$$\mathbb{E}(V_1)=\frac{{\rho }_1(c_1+c_2)}{c_2-{\rho }_1(c_1+c_2)}\left[\frac{1}{2}\frac{\mathbb{E}(B_1^2)}{\mathbb{E}(B_1)}+\frac{c_1}{(c_1+c_2{)}^2}\right],$$

$$\begin{array}{rl}\mathbb{V}\mathrm{ar}(V_1)=& \frac{{\rho }_1(c_1+c_2)}{c_2-{\rho }_1(c_1+c_2)}\\ & \left[\frac{1}{3}\frac{\mathbb{E}(B_1^3)}{\mathbb{E}(B_1)}+\frac{1}{4}\frac{{\rho }_1}{1-{\rho }_1}\frac{{\left(\mathbb{E}(B_1^2)\right)}^2}{{\left(\mathbb{E}(B_1)\right)}^2}+\frac{c_1}{(c_1+c_2{)}^2}\frac{\mathbb{E}(B_1^2)}{\mathbb{E}(B_1)}+\frac{c_1}{(c_1+c_2{)}^3}\right].\end{array}$$

$$\mathbb{E}(e^{-sY})=\frac{A_0}{1-{\rho }_1}\left[1-\frac{c_1}{c_1+c_2}\left(\frac{1}{A_0+s{A}_1-\frac{s^2}{2}{A}_2+o(s^2)}\right)\right],$$

$$A_0=$$

$$A_2=$$

$$\left(\frac{c_2}{c_1+c_2}-{\rho }_1\right)V_1\stackrel{d}{⟶}Z,$$

$$\mathbb{E}(e^{-s{V}_1})=\mathbb{E}(e^{-s{V}_{\text{M/G/1}}})\mathbb{E}(e^{-sY}).$$

$$\begin{array}{rl}& \underset{{\rho }_1\uparrow \frac{c_2}{c_1+c_2}}{\lim }\mathbb{E}\left(e^{-s\left(\frac{c_2}{c_1+c_2}-{\rho }_1\right)V_1}\right)\\ & =\underset{{\rho }_1\uparrow \frac{c_2}{c_1+c_2}}{\lim }\mathbb{E}\left(e^{-s\left(\frac{c_2}{c_1+c_2}-{\rho }_1\right)V_{\text{M/G/1}}}\right)\mathbb{E}\left(e^{-s\left(\frac{c_2}{c_1+c_2}-{\rho }_1\right)Y}\right).\end{array}$$

Full text

This paper is dedicated to Tomasz Rolski, in friendship, respect and admiration. His love of applied probability and never-ending curiosity are a blessing for our field.

Keywords: polling model; workload decomposition; heavy traffic; heavy tail asymptotics; singular perturbation analysis; time-scale separation; geometric ergodicity

1 Introduction

In this paper, we are interested in the performance analysis of a single server polling model with a special service discipline (i.e., the criterion which determines how many customers are served during a visit of the server to a queue). A typical polling model consists of multiple queues, attended by a single server who visits the queues in some order to render service to the customers waiting at the queues. Moving from one queue to another, the server incurs a (possibly zero) switch-over time. Once the server is at one of the queues, the server serves the customers of that queue based on a service discipline and according to some service time distribution.

Polling models were initially introduced in the 1950’s but mostly gained their popularity during the 1990’s. This popularity rise was due to the wide range of applicability of polling models, especially for the modelling of computer-communica-tion systems and protocols, traffic signal management, and manufacturing, see, e.g., [42, 43, 46] for a series of comprehensive surveys and [8, 33, 41] for extensive overviews of the applicability of polling systems.

The performance analysis of polling models has received considerable attention, see, e.g., [40]. In particular, in the polling literature much attention has been given to determining the probability generating function (PGF) of the joint queue length distribution under stationarity and at various epochs. A wide range of service disciplines has been considered, including exhaustive service (per visit to a queue, the server continues to serve all customers until it empties) and gated service (per visit to a queue, the server serves only those customers which are already present at the start of the visit). In [37], Resing shows that the joint queue length PGF of polling models in which the service discipline satisfies the so-called branching property equals the (known) PGF of a multi-type branching process with immigration. Service disciplines which satisfy the branching property include the exhaustive and gated disciplines. Polling systems with disciplines which do not satisfy the branching property usually defy an exact analysis. In our paper, we assume that the server spends an exponentially distributed amount of time at each queue. Upon the completion of this residing time at each queue, the server instantaneously switches to another queue according to a cyclic order. Such a service protocol does not exhibit the branching property, which complicates the analysis significantly. We concentrate on the two-queue model and, whenever possible, suggest extensions to the multi-queue model. A similar service discipline has been considered in [23], [1], [48], [21], and the references therein.

Related literature

In [23], Eliazar and Yechiali consider a multi-queue polling system under the Randomly Timed Gated (RTG) service discipline. The RTG discipline operates as follows: whenever the server enters a station, a timer is activated. If the server empties the queue before the timer’s expiration, the server moves on to the next queue. Otherwise (i.e., if there is still work in the station when the timer expires), the server obeys one of the following rules, each leading to a different model: (1) The server completes all the work accumulated up to the timer’s expiration and then moves on to the next node. (2) The server completes only the service of the job currently being served, and moves on. (3) The server stops working immediately and moves on. The model suggested in this manuscript bears resemblance to rule (3), however, in our case if a queue becomes empty, the server does not switch, and only does so when the timer expires. Eliazar and Yechiali, in [23], produce a recursive expression for the PGF of the number of customers in the queues of the polling model, while the case of two queues is sketched in [19], by a transformation to a boundary value problem. In [30], Katayama using a level-crossing approach, obtains the Laplace-Stieltjes transforms (LSTs) and the moment formulas for the waiting times and the sojourn times, and based on these expressions, he also proves a decomposition property.

In [1], the authors consider a polling model with Poisson batch arrivals and phase-type service times, and an exponential service timer. The authors establish a relation for the PGF of the number of customers in the queue at the beginning and at the end of the server’s visit to a queue. This is used as an input for a numerical scheme that is used to approximate the joint queue-length distribution at the server departure instants from the queues.

In [48], Xie et al. consider a single server multi-queue system, in which the server visits the individual queues for a fixed amount of time in a deterministic, cyclic order. They refer to the fixed residing time as the orientation time. They argue that such a service discipline comes with two operational advantages: it enables to i) keep the frequency of switching at a predetermined level (thus controlling the total cost, if there is a switching cost), ii) balance the time that the server spends in each queue (since, contrary to exhaustive or gated service disciplines, this discipline does not depend on the number of customers present in the various queues).

In [21], the authors assume a random visit (residing) time for each queue, which is independent of the number of customers present at each queue, and a preemptive-repeat with resampling service strategy. This autonomous service discipline is motivated from wireless ad hoc networks with movable communication hops. Another application is in single upstream tree-based ethernet passive optical networks, in which the central optical line terminal dedicates the channel to a specific user (e.g., the user with the highest priority) for a random amount of time, see [31] and the references therein. For more applications on this type of autonomous service disciplines, the interested reader is referred to [2]. For all aforementioned applications, we consider it natural to assume that the service strategy is preemptive-resume and that the switch-over time is negligible in comparison to the service time and the residing time.

Paper originality

In this paper, we initially devote attention to the individual queues. When focusing on a single queue, the model can be interpreted as a service system with vacations: we interpret the time that the server visits the other queue as a vacation period. Vacation queues - and priority queues for which the mathematical analysis is similar - are well studied in the queueing literature starting with the work of White and Christie [47] (exponentially distributed service times and vacations), Gaver[28], Thiruvengadam [45] and Avi-Itzhak and Naor [5] (the latter three assuming generally distributed service times and vacations). All these works assume that the service periods have an exponential distribution, but vary for example in the assumptions regarding whether interrupted services are resumed or repeated and in the metrics of interest. Takagi [42] provides an excellent overview of vacation and priority systems. The interested reader is also referred to Federgruen and Green [25] for phase-type distributed service periods, to Takine and Sengupta [44] for Markovian arrival processes and to Fiems et al. [26] for a more recent publication with various sorts of service disruptions. For a more extensive overview of the literature, we refer to the recent survey of Krishnamoorthy et al. [32].

A particular feature of a large class of vacation queues is that the stationary workload and queue length distributions obey a stochastic decomposition property, as first observed by Gaver [28] and Miller [35]. Fuhrmann and Cooper [27] give conditions for such a queue length decomposition to hold. Our model does not satisfy these conditions, but we show that it does allow a stochastic decomposition of the stationary workload. It is notable that the workload can be decomposed into two independent terms: the amount of work of an M/G/1 queue, and the amount of work when the server is not serving the first queue, due to either an idle period or due to a visit at the second queue. The second term in the decomposition is at first sight surprising, as one would typically expect that the decomposition is identical with that of the standard M/G/1 queue with exhaustive service and multiple vacations, as is the case for the RTG, see [30, Remarks 2 and 4]. In that respect, the exact analysis of the marginal workload distribution was quite helpful to us. The use of the decomposition property further facilitates our heavy traffic and heavy tail analyses, as we can use known results for the M/G/1 queue and restrict our attention to the second term appearing in the stochastic decomposition.

In this paper, we also devote attention to the joint workload distribution. Using similar probabilistic arguments as in the analysis of the marginal workload, we show that the joint LST of the workloads at the queues satisfies a functional equation (6.4), which is, in the case of identical queues, then reduced to a Dirichlet boundary value problem. Thus, one can numerically evaluate the mapping from the contour defined by the kernel of the functional equation to the unit circle, and obtain a solution to the joint distribution. However, we have to note that, depending on the service time distribution, the kernel of the functional equation does not have the typical (quadratic) polynomial form, which complicates tremendously the analysis and differentiates it from the known results of the literature, see, e.g., [18, 19, 24]. In the case of exponentially distributed service times, we propose perturbation analysis for the calculation of the joint queue length distribution, as a methodological alternative to the boundary value problem approach. By appropriately scaling the arrival and service rates by a factor $\varepsilon$, the invariant probability measure of the perturbed Markov chain is written as a power series expression in terms of $\varepsilon$, whose coefficients form a geometric sequence, that can be used for both exact and numerical calculations. Furthermore, we show that there exists a computationally stable updating formula for the calculation of the perturbed invariant measure. To approximate the joint distribution numerically, one needs to solve a large system of equations, for which we indicate two possible approaches, but we do not pursue these in this paper.

Paper overview

The paper is organized as follows: In Section 2, we describe the two-dimensional polling model under consideration. In Section 3, we present the LSTs of the model’s marginal workload distributions in steady-state at an arbitrary epoch. In Section 4, we show that a single queue’s marginal workload satisfies a decomposition property and then by using the decomposition property in the light tailed case we obtain the heavy traffic limit of the marginal workload distributions in steady-state. In Section 5, we discuss the heavy tail asymptotics of the marginal workload distributions in steady-state, and in that case we also discuss the heavy traffic behavior. We then discuss, in Section 6, open problems arising in the calculation of the joint workload distribution. Assuming exponentially distributed service times, we calculate in Section 7 the joint queue length distribution in steady-state at an arbitrary epoch using singular perturbation analysis. Several possible future research directions are discussed in Section 8.

2 Model description and notation

In this paper, we consider a two-queue polling model. Customers arrive to queue $i$ according to a Poisson process at rate $\lambda_{i}$, $i=1,2$. There is a single server, that serves both queues according to the first come first serve (FCFS) discipline. The service times of customers in queue $i$ are independent and identically generally distributed positive random variables, say $B_{i}$, $i=1,2$. We denote the LST of the service time $B_{i}$ by $\tilde{b}_{i}(s)=\mathbb{E}(e^{-sB_{i}})$, with $\mathrm{Re}\leavevmode\nobreak\ s\geq 0$, $i=1,2$.

A special feature of the polling model under consideration is that the server spends an exponentially distributed amount of time at queue $i$ with rate $c_{i}$, $i=1,2$. Upon completion of the residing time at queue $i$, the server instantaneously switches to the other queue, i.e., there is no switch-over time. Furthermore, if upon completion of the residing time, the server is providing service to a customer, this service is interrupted and resumed at the next visit of the server to the queue. More explicitly, we assume that if a server resumes the service after being interrupted, the server continues from where the service stopped instead of starting from the beginning, i.e., the service is preemptive–resume. We denote the LST of the residing time $T_{i}$ of the server in queue $i$ by $\tilde{f}_{T_{i}}(s)=\mathbb{E}(e^{-sT_{i}})$, with $\mathrm{Re}\leavevmode\nobreak\ s\geq 0$, $i=1,2$, with $T_{i}$ exponentially distributed with rate $c_{i}$ and probability density function $f_{T_{i}}(t)=c_{i}e^{-c_{i}t}$, $t\geq 0$, $i=1,2$.

Stability condition

For the two-queue polling model under consideration the stability condition (sufficient and necessary) is

[TABLE]

with $\rho_{i}=\lambda_{i}\,\mathbb{E}(B_{i})$, $i=1,2$.

The stability condition can be proven by appropriately adapting and extending the proof of Altman et al. in [4]. To this purpose, one would need to calculate the expected increase of the workload during a cycle (i.e., the time between two successive arrivals of the server at the first queue) and use an extension of Foster’s criterion known as the positivity/regularity criterion (V2), cf. [34]. The derivations of the expected increase of the workload during a cycle are similar to the analysis performed for the proof of Theorem 3.1. This guarantees that the stability condition is sufficient. In order to show that it is also necessary, one may use the expression for the expectation of the steady-state workload of each queue, cf. Equation (4.3).

Remark 2.1**.**

Equivalently, one can prove the stability condition by adapting the steps presented in [15]. More concretely, the two-queue polling system under consideration is said to be stable if the workloads at each queue, at the polling instants (i.e., the instant when the server arrives at a queue), have a proper limiting distribution, and the mean cycle time is finite, as time tends to infinity. By definition, the latter always holds, as the cycle time for the polling model under consideration is given by the sum of the two exponentially distributed residing times at the two queues. Following a similar argumentation as in [15], one would need to ensure that, under the above stability condition, each queue is stable. Note that each queue in isolation behaves like an M/G/1 queue with a service speed governed by a two-state Markov chain. The stability condition (in the sense that the workload has a proper limiting distribution) of such a queueing system is studied in [6].

Remark 2.2**.**

Intuitively, the stability condition for the first queue can be interpreted as follows: the long-run proportion of time the server spends in the first queue is equal to $c_{2}/(c_{1}+c_{2})$, thus the long-run rate of service in the first queue is $c_{2}/(c_{1}+c_{2})\mathbb{E}(B_{1})$. Hence, for the first queue to be stable it is needed that the arrival rate is strictly smaller than the long-run rate of service, which corresponds to the left hand side of (2.1). The stability condition for the second queue can be interpreted in an analogous manner.

3 Marginal workload analysis

In this section, we derive the distribution of the marginal workload in steady-state at an arbitrary epoch. As discussed in the introductory section, the individual queues behave as vacation systems: from the perspective of one queue, the server is on vacation when it resides at the other queue. In this section, we give a direct derivation of the stationary marginal workload distributions.

We let $V_{i}(t)$ denote the workload at time $t$, $t\geq 0$, of queue $i$, $i=1,2$, and $V_{i}$ denote the steady-state workload of queue $i$ at an arbitrary epoch, $i=1,2$.

Theorem 3.1**.**

The LST of the workload of the first queue in steady-state under the stability condition (2.1) is given by

[TABLE]

A symmetric formula holds for the LST of $V_{2}$ under the stability condition (2.1).

Proof.

The derivation of the LST of the steady-state workload for the first queue is performed by considering the renewal process at the instances the server arrives at the first queue, i.e., the inter-renewal times are identical in distribution to $T_{1}+T_{2}$, with $T_{i}\sim\mathrm{Exp}(c_{i})$, $i=1,2$.

To structure the exposition, the proof of the theorem is split into five steps. A key point of the proof is the derivation of $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})})$; this is achieved in Step 4, after we derive $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})}|V_{1}(T_{1})=y)$ in Step 1, $\mathbb{E}(e^{-sV_{1}(T_{1})}|V_{1}(0)=v)$ in Step 2, and subsequently $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})}|V_{1}(0)=v)$ in Step 3. Finally, in Step 5, we calculate $\mathbb{E}(e^{-sV_{1}})$ using the PASTA property and the result of Step 4.

Step 1: Calculation of $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})}|V_{1}(T_{1})=y)$.

During $(T_{1},T_{1}+T_{2}]$ the server serves only customers in the second queue, so the workload in the first queue only increases by the sum of the service times of all the customers that arrived within this interval. The increments occur according to a compound Poisson process. So,

[TABLE]

Step 2: Calculation of $\mathbb{E}(e^{-sV_{1}(T_{1})}|V_{1}(0)=v)$.

Note that

[TABLE]

In order to calculate the right hand side of (3.3), we use [17, p. 262, Equation (4.99)]

[TABLE]

with $\mathrm{Re}\leavevmode\nobreak\ s\geq 0$, $t\geq 0$, and $U_{1}(x)=0$, if $x<0$, and $U_{1}(x)=1$, otherwise. Hence, Equation (3.3) in light of [17, p. 262, Equation (4.99)] yields

[TABLE]

For the calculation of the integrals in the right hand side of Equation (3.4) we use [17, p. 260,  Equation (4.92)], for $\mathrm{Re}\leavevmode\nobreak\ s\geq 0$, $t\geq 0$,

[TABLE]

with $\mu(s,1)$ denoting the LST of the busy period distribution of the M/G/1 queue with arrival rate $\lambda_{1}$ and service time LST $\tilde{b}_{1}(s)$; $\mu(s,1)$ is the root of $z=\tilde{b}_{1}\left(s+(1-z)\lambda_{1}\right)$ with the smallest absolute value, cf. [17, p. 250]). A lengthy but straightforward calculation, that involves interchanging the integrations, yields, for $\mathrm{Re}\leavevmode\nobreak\ s\geq 0$,

[TABLE]

Combining (3.4) and (3.5) yields

[TABLE]

Step 3: Calculation of $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})}|V_{1}(0)=v)$.

[TABLE]

where the second equation comes from Equation (3.2) and the third from Equation (LABEL:Step2).

Step 4: Calculation of $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})})$ in steady-state.

Observe that

[TABLE]

with $f_{V_{1}(0)}(v)$ the probability density function of $V_{1}(0)$. Now observe that in steady-state $V_{1}(T_{1}+T_{2})$ has the same distribution as $V_{1}(0)$. So we can rewrite (3.8) as follows:

[TABLE]

So,

[TABLE]

Taking the limit as $s\rightarrow 0$ in (3.9) and using L’Hôpital’s rule yields

[TABLE]

Hence,

[TABLE]

Step 5: Calculation of $\mathbb{E}(e^{-sV_{1}})$ in steady-state.

Firstly, let us denote by $S=1$ (respectively by $S=2$) the event of the server residing in the first (respectively second) queue. Then,

[TABLE]

Because of the memoryless property of the exponential distribution it is obvious that

[TABLE]

The latter term is given by (3.10), while the former term is calculated using the same argument as in the derivation of Equation (3.2):

[TABLE]

Substituting (3.12), for $\mathbb{E}(e^{-sV_{1}(T_{1})})$, and (3.10) in Equation (3.11) yields (3.1), which concludes the proof.

Similarly, we can also calculate the LST of the workload of the second queue. ∎

Remark 3.1**.**

It is not difficult to extend the above results to the case that the $T_{2}$ periods are non-exponential, see, e.g., [28, 45, 5], and to the case that the arrival process during those periods is a different compound Poisson process than during the $T_{1}$ periods, see, e.g., [44] and [26]. One could even allow a more general non-decreasing Lévy process (subordinator) during those $T_{2}$ periods. During $T_{1}$ periods, one could also allow the input process to be a subordinator. However, we do note that it is considerably more complicated to consider non-exponential $T_{1}$ periods, see, [25].

4 Workload decomposition and Heavy Traffic analysis

In this section, we show that the steady-state workload $V_{1}$ (similarly for $V_{2}$) can be decomposed into two independent terms, one corresponding to the steady-state workload of the first queue in isolation, i.e., the M/G/1 queue with arrival rate $\lambda_{1}$ and service times $B_{1}$ (to be called: corresponding M/G/1 queue), and the second corresponding to the amount of work when the server is not serving the first queue, due to either an idle period or due to a visit at the second queue.

Assuming that the first three moments of $B_{1}$ are finite and then using the decomposition of $V_{1}$, we determine the mean, the variance, and the heavy traffic limit of the workload $V_{1}$. Furthermore, in this and in the next section, we use the decomposition to obtain various asymptotic (heavy traffic and/or heavy tail) results.

Corollary 4.0.1**.**

The steady-state amount of work of the first queue, $V_{1}$, is distributed as a sum of two independent random variables $V_{\text{\tiny M/G/1}}$ and $Y$, i.e.,

[TABLE]

where $V_{\text{\tiny M/G/1}}$ is the steady-state amount of work in the corresponding M/G/1 queue, and $Y$ is the steady-state amount of work when the server is not serving at the first queue. The LST of the random variable $Y$ is given as

[TABLE]

Proof.

The workload decomposition result follows from [10, Theorem 2.1]; it is readily verified that all conditions of that theorem are satisfied. And the LST of $Y$ can be directly obtained by dividing the LST of $V_{1}$ (which is given by Eqn. (3.1)) by the known LST of the M/G/1 queue (cf. [17, p. 257, Equation (4.90)]). ∎

Remark 4.1**.**

$\mathbb{E}(e^{-sY})$* could also have been obtained by writing it as a weighted sum of two known terms: (i) 1, which is the LST of the (zero) workload in the first queue during an idle part of the visit period, and (ii) $\mathbb{E}(e^{-sV_{1}(T_{1}+T_{2})})$, which is given in (3.10). PASTA implies that the latter term is also the LST of the workload in the first queue at an arbitrary epoch of a visit period of the other queue.*

We now use the decomposition result (4.1) to determine the mean and the variance of $V_{1}$.

Theorem 4.1**.**

The expectation of the steady-state workload of the first queue, $\mathbb{E}(V_{1})$, is

[TABLE]

and the corresponding variance, $\mathbb{V}\mathrm{ar}(V_{1})$, is

[TABLE]

Proof.

The mean and variance can be obtained by using the decomposition result (4.1). For this purpose, we can separately calculate the mean and the variance of the M/G/1 queue, cf. [17, p. 256], as well as the mean and the variance corresponding to the random variable $Y$. For the latter we use Equation (4.2) (after dividing its numerator and denominator by $s$). Combining these results yields Equations (4.3) and (4.4). ∎

Remark 4.2**.**

Equation (4.3) and Equation (4.4) for $c_{2}\to\infty$ (or equivalently for $c_{1}\to 0$) yield the corresponding expressions for the mean and the variance of the M/G/1 queue, cf. [17, p. 256].

Now, we study the behavior of the workload $V_{1}$ in heavy traffic, i.e., when $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$. In Corollary 4.0.1 we have shown that $V_{1}$ can be written as the sum of the independent random variables $V_{\text{\tiny M/G/1}}$ and $Y$. Since most of the results related to the M/G/1 queue are already known, we take a closer look at $\mathbb{E}(e^{-sY})$, with the assumption that the first three moments of $B_{1}$ are finite. Substituting $\tilde{b}_{1}(s)=1-s\mathbb{E}(B_{1})+\frac{s^{2}}{2}\mathbb{E}(B_{1}^{2})-\frac{s^{3}}{3!}\mathbb{E}(B_{1}^{3})+o(s^{3})$ in (4.2) and rearranging the terms yields, for $s\downarrow 0$,

[TABLE]

with

[TABLE]

Equation (4.5) will play a very important role in the proof of Theorem 4.2 and also in the next section where we study the tail behavior of the workload $V_{1}$.

Theorem 4.2**.**

Assume that $\mathbb{E}(B_{1}^{2})<\infty$. For $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$,

[TABLE]

with $Z$ an exponentially distributed random variable with mean $\frac{c_{1}c_{2}}{(c_{1}+c_{2})^{3}}+\frac{c_{2}}{c_{1}+c_{2}}\frac{1}{2}\frac{\mathbb{E}(B_{1}^{2})}{\mathbb{E}(B_{1})}$.

Proof.

To obtain the heavy traffic limit of $V_{1}$ one can use the workload decomposition. Corollary 4.0.1 implies that

[TABLE]

Replacing $s$ by $sA_{0}=s(\frac{c_{2}}{c_{1}+c_{2}}-\rho_{1})$, cf. (4.6), in the above equation and taking the limit $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$ yields

[TABLE]

The first term in the right hand side obviously tends to one for $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$, as the corresponding M/G/1 queue is in heavy traffic only when $\rho_{1}\uparrow 1$. In order to calculate the limit for the second term in (4.10), we replace $s$ by $sA_{0}=s(\frac{c_{2}}{c_{1}+c_{2}}-\rho_{1})$, cf. (4.6), in (4.5), which yields

[TABLE]

Taking the limit $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$ in (4.11) yields

[TABLE]

with $A_{1}$ given in (4.6). From (4.12), (4.6), and (4.10) the statement of the theorem follows by noticing that the right hand side of (4.12) corresponds to the LST of an exponentially distributed random variable with mean $A_{1}$. ∎

Remark 4.3**.**

Letting $c_{2}\to\infty$, Theorem 4.2 indicates that the heavy traffic result reduces to that of an ordinary M/G/1 queue.

5 Heavy tail asymptotics

In this section, we discuss the tail behavior of the workload in the case of heavy tailed service time distributions. We also study the heavy traffic behavior of the workload $V_{1}$ when the service time distribution $B_{1}$ is regularly varying. To do this analysis, we now introduce the definition of a regularly varying random variable/distribution.

Definition 5.1**.**

The distribution function of a random variable $B_{1}$ on $[0,\infty)$ is called regularly varying of index $-\nu$, with $\nu\in\mathbb{R}$, if

[TABLE]

with $L(x)$ a slowly varying function at infinity, i.e., $\lim\limits_{x\to\infty}\frac{L(\alpha x)}{L(x)}=1$, for all $\alpha>1.$

Theorem 5.1**.**

If the service time distribution of the random variable $B_{1}$ is regularly varying of index $-\nu$, with $\nu\in(1,2)$, then the workload of the first queue under the stability condition (2.1) is regularly varying at infinity of index $1-\nu$. More precisely,

[TABLE]

Proof.

To prove that $V_{1}$ is regularly varying at infinity, one can again use the decomposition property of the workload $V_{1}$. From Corollary 4.0.1, we get

[TABLE]

In the M/G/1 queue it follows from [16] that $\mathbb{P}(V_{\text{\tiny M/G/1}}>x)$ is regularly varying of index $1-\nu$ at infinity if and only if the tail of the service time distribution $\mathbb{P}(B_{1}>x)$ is regularly varying of index $-\nu$ at infinity, and one has

[TABLE]

Now we have to compute $\mathbb{P}(Y>x)$ for $x\uparrow\infty$. Our main tool is the Tauberian theorem of [7, Theorem 8.1.6], which relates the behavior of a regularly varying function at infinity and the behavior of its LST near [math]. Applying this theorem to (5.1) gives

[TABLE]

and hence

[TABLE]

Substituting Equation (5.5) in (4.2) yields, for $s\downarrow 0$:

[TABLE]

Simplifying, we get

[TABLE]

Applying the Tauberian theorem of [7, Theorem 8.1.6] once again, now in the reverse direction, yields

[TABLE]

From Equation (5.4) and (5.6), we see that both $V_{\text{\tiny M/G/1}}$ and $Y$ are regularly varying random variables of index $1-\nu$. Using the workload decomposition property (4.1) and a well known result regarding the tail behavior of the sum of two independent regularly varying random variables of the same index, see [38], yields

[TABLE]

with $C_{1}$ and $C_{2}$ the coefficients of the tail $x^{1-\nu}$ for $V_{\text{\tiny M/G/1}}$ and $Y$ in (5.4) and (5.6), respectively. Substituting the coefficients from (5.4) and (5.6) concludes the proof of the theorem. ∎

Remark 5.1**.**

Letting $c_{2}\to\infty$ in Equation (5.7) yields

[TABLE]

which is, as expected, the result for an ordinary M/G/1 queue.

Remark 5.2**.**

Theorem 5.1 is closely related to [14, Theorem 4.1] for a single server queue with alternating high and low service speeds. In [14] both the service time distribution and the distribution of the periods of low service speed are regularly varying. If the latter tail is less heavy than the tail of the service time distribution, then our formula (5.2) displays exactly the same tail behavior as [14, Formula (4.1)].

In the next theorem we discuss how Theorem 5.1 can be generalized to the case of subexponential (residual) service times.

Definition 5.2**.**

A distribution function $\mathbb{P}(B_{1}\leq x)$, $x\geq 0$, is called subexponential if

[TABLE]

for any $n\geq 2$, with $B_{11},\dots,B_{1n}$ independent and identical copies of $B_{1}$.

Theorem 5.2**.**

If the distribution of the residual service time requirement, say $B_{1}^{r}$, is subexponential, then $V_{1}$ is also subexponential and

[TABLE]

Heuristic proof. The asymptotic relation in (5.9) can be proved formally using sample-path techniques along the following lines. We assume the system is in stationarity and focus on the workload at time $t=0$. If the workload level $V_{1}$ at this time is very large, then that is most likely due to the prior arrival of a customer with a large service requirement $B_{1}$, at some time $t=-y$. We can observe that from time $t=-y$ onward, the workload decreases nearly linearly with rate $\frac{c_{2}}{c_{1}+c_{2}}-\rho_{1}$. So in order for the workload at time $t=0$ to exceed the level $x$, the service requirement $B_{1}$ must be larger than $x+y\left(\frac{c_{2}}{c_{1}+c_{2}}-\rho_{1}\right)$. Since customers arrive according to a Poisson process with rate $\lambda_{1}$, the distribution of the workload $V_{1}$ for large $x$ can be computed as

[TABLE]

A change of variable $z:=x+y\left(\frac{c_{2}}{c_{1}+c_{2}}-\rho_{1}\right)$ in Equation (5.10) yields

[TABLE]

which leads to Relation (5.9).

This proof can be made rigorous by providing lower and upper bounds for $\mathbb{P}(V_{1}>x)$ that in the limit coincide. The lower bound is easily obtained by using the law of large numbers. The upper bound is more difficult; one has to give a formal version of the statement “exceedance of a high level $x$ happens as a consequence of a single big jump”, and one has to show that other exceedance scenarios (like two rather big jumps) do not contribute to the asymptotics of the exceedance probability. We refer to [49, Section 2.4] for a detailed exposition of this technique.

Remark 5.3**.**

Note that, indeed, Relation (5.9) contains the result of Theorem 5.1 as a special case, since $B_{1}$ being regularly varying at infinity of index $-\nu$, with $\nu\in(1,2)$, has a subexponential distribution. In this regularly varying case, we have

[TABLE]

In the above equation by applying the regular varying function property from [7, p. 26, Proposition 1.5.8], we get

[TABLE]

Combining (5.11) and (5.12), we obtain Theorem 5.1.

Now we concentrate on a heavy traffic limit theorem for $V_{1}$ in the heavy tailed case. To do this analysis, we first scale $V_{1}$ by the coefficient of contraction $\Delta(\rho_{1})$. Similarly to [11, p. 188, Equation (4.24)], we define the coefficient of contraction $\Delta(\rho_{1})$ as the unique root of the following equation in $x$

[TABLE]

such that $\Delta(\rho_{1})\downarrow 0$ for $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$.

Theorem 5.3**.**

If the service time distribution of the random variable $B_{1}$ is regularly varying of index $-\nu$, with $\nu\in(1,2)$, then the heavy traffic limiting distribution of workload $V_{1}$ of the first queue in the heavy tailed case is given by the Mittag-Leffler distribution:

[TABLE]

Proof.

We can obtain the heavy traffic limit in the heavy tailed case by using the workload decomposition property (4.0.1) and its LST version (4.9). The heavy traffic limit can be computed by replacing $s$ by $s\Delta(\rho_{1})$ in Equation (4.9) and taking the limit $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$, which yields

[TABLE]

Just as in the light tailed case (cf. Theorem 4.2), the contribution of $V_{\text{\tiny M/G/1}}$ becomes negligible compared to the contribution of $Y$. To calculate the limit for the second factor in (5.15), we use (4.2).

[TABLE]

with $f(s\Delta(\rho_{1}))=\frac{\rho_{1}(1-\tilde{b}_{1}(s\Delta(\rho_{1})))}{\mathbb{E}(B_{1})}$. Taking the limit $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$ in (5.16) yields

[TABLE]

since $\frac{s\Delta(\rho_{1})}{c_{1}+c_{2}}\left(\frac{f(s\Delta(\rho_{1}))}{s\Delta(\rho_{1})}\right)^{2}\to 0$, $f(s\Delta(\rho_{1}))\to 0$ and $\Delta(\rho_{1})\to 0$ when $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$. After rearranging the terms of (5.17) we get

[TABLE]

Since $B_{1}$ is regularly varying, we get by using [11, Lemma 5.1 (iv)],

[TABLE]

Using [11, p. 188, Equation (4.22)], we know that

[TABLE]

From the definition of the coefficient of contraction $\Delta(\rho_{1})$ as the unique root of Equation (5.13) such that $\Delta(\rho_{1})\downarrow 0$ for $\rho_{1}\uparrow\frac{c_{2}}{c_{1}+c_{2}}$, we have

[TABLE]

Furthermore, from the definition of a slowly varying function $L(\cdot)$, we get $\frac{L\left(\frac{1}{s\mathbb{E}(B_{1})\Delta(\rho_{1})}\right)}{L\left(\frac{1}{\Delta(\rho_{1})}\right)}\to 1$, when $\Delta(\rho_{1})\downarrow 0$. Now by combining (5.18) - (5.21) we get

[TABLE]

Substituting (5.22) in Equation (5.15) concludes the proof of the theorem. ∎

Remark 5.4**.**

In [11] a class of service time distributions is considered that is slightly larger than the class of regularly varying distributions. Theorem 5.3 can be seen to hold under these conditions as well.

Remark 5.5**.**

At this stage it may be appropriate to discuss some of the advantages and disadvantages of working with timers instead of traditional polling disciplines like exhaustive, gated or $k$-limited. An advantage of working with timers appears to be the following. It is proven in [12] for $N$-queue cyclic polling models with exhaustive or gated service, that if the heaviest service time distribution (say, at $Q_{M}$) is regularly varying of index $-\nu$, then all waiting time distributions are regularly varying of index $1-\nu$. This is intuitively clear: There is a positive probability of a customer in $Q_{1}$ arriving during a service time of $Q_{M}$, and then its waiting time includes a residual type-$M$ service time – which is regularly varying of index $1-\nu$. That intuition also indicates that a similar tail behavior will occur for disciplines like $k$-limited.

On the contrary, in our model with exponential timers, the service time distribution at one queue has no effect at all on the waiting time distribution or workload distribution at the other queue (for the workload, this is also seen from Formula (3.1), which does not involve $\tilde{b}_{2}(s)$). In particular, the waiting time distribution at $Q_{1}$ will asymptotically behave exponentially if the service time tail is exponential; and if the service time distribution at $Q_{1}$ is regularly varying of index $-\zeta$ while the service time distribution at $Q_{2}$ is regularly varying of index $-\nu$ with $\zeta>\nu$, then the waiting time distribution will be regularly varying of index $1-\zeta<1-\nu$.

Next to protection against heavy tails at other queues, there is also some protection against long mean service times at other queues. If $\operatorname{\mathbb{E}}B_{2}$ is much larger than $\operatorname{\mathbb{E}}B_{1}$, then a type-$1$ customer in an ordinary polling model may experience a long mean delay because of the presence of type-$2$ customers. However, when timers are used, a customer of type $1$ with a short service time will not suffer much from the presence of type-$2$ customers with long mean service times. This is a similar protection phenomenon as round robin or processor sharing protecting short customers against having to wait long for a customer with a very long service time.

A disadvantage of using timers is that the system is not work conserving, and hence operates in a sense less efficient than one would wish. This is revealed in the fully symmetric case. It follows from (4.3) that the mean workload at $Q_{1}$, in the fully symmetric case (equal arrival rates, service time distributions and mean visit periods), equals

[TABLE]

so that the mean total workload equals

[TABLE]

On the other hand, in a symmetric two-queue polling model with, e.g., exhaustive or gated service, one has work conservation, and hence the mean total workload equals

[TABLE]

This is less than half the value of $\operatorname{\mathbb{E}}V_{1}+\operatorname{\mathbb{E}}V_{2}$. The main reason for this is the fact that the server is sometimes idle although there is work at the other queue. This even holds when $c_{1}=c_{2}\rightarrow\infty$; $Q_{1}$ operates effectively as if the server works at half speed for it, or as if the service times are twice as long.

6 Joint workload distribution

So far we have focused on the marginal workload distribution at the first queue. A much harder problem is to determine the steady-state joint workload distribution. In this section, we begin the exploration of this problem, outlining a possible approach as well as the mathematical complications arising.

Let $\tilde{v}(s_{1},s_{2}):=\mathbb{E}({e}^{-s_{1}V_{1}(T_{1}+T_{2})-s_{2}V_{2}(T_{1}+T_{2})})$ denote the steady-state joint workload LST at endings of visit periods at the second queue. Reiterating Steps 1 - 4 of Section 3, but now taking both workloads into account, leads after lengthy calculations to the following functional equation

[TABLE]

with $\omega_{1}(s_{2}):=c_{1}+\lambda_{2}(1-\tilde{b}_{2}(s_{2}))+\lambda_{1}(1-\mu(\zeta,1))$; as before, $\mu(s,1)$ is the busy period LST of the M/G/1 queue in isolation corresponding to the first queue.

Let us now restrict ourselves to the fully symmetric case $c_{1}=c_{2}=c$, $\lambda_{1}=\lambda_{2}=\lambda$, $\tilde{b}_{1}(s)=\tilde{b}_{2}(s)=\tilde{b}(s)$. Formula (6.1) then becomes

[TABLE]

Taking $s_{1}=s_{2}$ in (6.2) allows us to express $\tilde{v}(\omega_{1}(s_{2}),s_{2})$ in terms of $\tilde{v}(s_{2},s_{2})$, thus reducing (6.2) to

[TABLE]

Interchanging all indices, one obtains a mirrored equation of (6.3), and the two equations combined yield

[TABLE]

with $K(s_{1},s_{2})=c^{2}-\left(c-s_{1}+\lambda(1-\tilde{b}(s_{1}))+\lambda(1-\tilde{b}(s_{2}))\right)\left(c-s_{2}+\lambda(1-\tilde{b}(s_{1}))\right.\\ \left.+\lambda(1-\tilde{b}(s_{2}))\right)$. This is a so-called boundary value problem equation. Equations of this type have been studied in the monograph [18]. There an approach is outlined that, for the present problem, amounts to the following global steps:

Step 1:

Consider the zeros of the kernel equation $K(s_{1},s_{2})$, that have $\mathrm{Re}\leavevmode\nobreak\ s_{1},\ \mathrm{Re}\leavevmode\nobreak\ s_{2}\geq 0$. For such pairs $(s_{1},s_{2})$, $\tilde{v}(s_{1},s_{2})$ is analytic, and hence, for those pairs, the right hand side of (6.4) is equal to zero.

Step 2:

For the pairs $(s_{1},s_{2})$ satisfying Step 1, one needs to translate the fact that the right hand side of Equation (6.4) is zero into a Riemann or Riemann-Hilbert boundary value problem. The solution of such a problem yields $\tilde{v}(s_{1},s_{1})$ and $\tilde{v}(s_{2},s_{2})$. Then $\tilde{v}(s_{1},s_{2})$ follows via (6.4).

Unfortunately, the above steps do not constitute a simple, straightforward recipe. For example, several choices of zero pairs are possible in the present problem, and it is not a priori clear what is the best choice. A natural choice, due to the symmetry of the underlying problem, seems to be to restrict oneself to complex conjugate points, i.e., choose $(s_{1},s_{2})=(z,\bar{z})$. The kernel then becomes

[TABLE]

Taking

[TABLE]

indeed yields $K(z,\overline{z})=K\left(z(\theta),\overline{z}(\theta)\right)=0$, for all $\theta\in[0,2\pi]$, while it is readily checked that for each such $\theta$ there is a unique $z(\theta)$ with $\mathrm{Re}\leavevmode\nobreak\ z(\theta)\geq 0$.

Turning to Step 2, one sees that the $z(\theta)$ satisfying (6.5) describe a closed contour, say $L$, in the right half plane, for $\theta:0\rightarrow 2\pi$, while the fact that the right hand side of (6.4), after a division by $s_{1}s_{2}=z\bar{z}$, is zero for all these $(s_{1},s_{2})=(z(\theta),\overline{z}(\theta))$ translates into the following relation

[TABLE]

with $\tilde{v}(z,z)$ and $\tilde{b}^{r}(z)=\frac{1-\tilde{b}(z)}{z\mathbb{E}B}$ analytic inside $L$. The expression inside the square brackets of (6.6) is analytic inside $L$ apart from the pole at $z=0$. This results in a similar boundary value problem as has been treated in [13]. The solution of such a problem is known when $L$ is the unit circle. For other closed contours, one needs a conformal mapping of that contour onto the unit circle; several procedures are available for obtaining such conformal mappings.

Remark 6.1**.**

In a future study we aim to handle all the technicalities which arise in treating this boundary value problem with a pole. When we manage to solve the present symmetric problem, we are still faced with the more general asymmetric two-queue problem. Subsequently, one could turn to the joint queue length distribution. However, a complication there is that a switch between queues might occur during a service time, forcing one to keep track of the length of the residual service time. From that perspective, workload seems to be an easier quantity than queue length.

In the next section, we turn our focus to the joint queue length distribution, but restricting ourselves to exponential service time distributions, so we do not need to keep track of the residual service time. Instead of pursuing a boundary value approach, we explore a perturbation approach, which allows us to derive an analytic expansion for the joint queue length distribution.

7 Joint queue length distribution

In this section, we turn our attention to the steady-state joint queue length distribution, restricting ourselves to exponential service requirement distributions in both queues, with rates $\mu_{i}=1/\mathbb{E}(B_{i})$, $i=1,2$, respectively. Under this assumption, we do not need to keep track of the residual service times, which simplifies the analysis. However, a direct analytic derivation of the joint queue length distribution (or its PGF) turns out to be as challenging as the analysis presented in Section 6. To address this issue, in this section, we explore the use of parametric perturbation for the derivation of the joint queue length distribution. In what follows, we use the framework developed in Altman et al. [3]; we perturb the service and arrival rates by a common parameter, denoted by $\varepsilon\geq 0$, i.e., the perturbed service rate of the customers in queue $i$ is $\varepsilon\mu_{i}$, $i=1,2$, and arrivals occur according to two independent Poisson processes with perturbed rates $\varepsilon\lambda_{i}$, $i=1,2$. The parameters that are not perturbed are $c_{i},i=1,2$ i.e., the rates of the exponentially distributed durations that the server spends in each queue. Note that the stability condition (2.1) is not affected by this scaling.

The perturbed process is a continuous time Markov chain defined on the state space

[TABLE]

in which $n_{i}$ denotes the queue length in queue $i$, $i=1,2$, and the third element in the state space description reports the queue in which the server is active. Furthermore, let $\bm{G}(\varepsilon)$ denote the generator of the perturbed Markov process. We decompose this perturbed generator into the unperturbed generator $\bm{G}^{(0)}$ and the perturbation matrix $\bm{G}^{(1)}$,

[TABLE]

so as to investigate the dependence of the stationary joint queue length distribution on the parameter $\varepsilon$. The unperturbed generator $\bm{G}^{(0)}$ corresponds to the states depicting a change of the state of the server from one queue to the other; it is given by

[TABLE]

with $\bm{C}=\begin{bmatrix}-c_{1}&c_{1}\\ c_{2}&-c_{2}\\ \end{bmatrix}$, and $\bm{0}_{2\times 2}$ a $2\times 2$ matrix of zeros. Throughout the remainder of the paper we use this notation with subscripts to indicate the dimension when needed. When the dimension is clear from the context, the index is omitted; note that the dimension can be infinite.

The perturbation matrix $\bm{G}^{(1)}$ is defined in terms of its elements, with $n_{1}\geq 0$, $n_{2}\geq 0$, $k=1,2$,

[TABLE]

with $\mbox{1{1}}_{\{n_{k}\geq 1\}}$ an indicator function taking value 1, if $n_{k}\geq 1$, and 0, otherwise.

In order to implement the framework of Altman et al. [3], it is convenient to first define the transition probability matrix $\bm{P}(\varepsilon)=\bm{I}+\Delta\bm{G}(\varepsilon)$ of the corresponding (uniformized) discrete time perturbed Markov chain ($\bm{I}$ being the identity matrix). In order to simplify notation, in what follows, we assume without loss of generality that

[TABLE]

Note that indeed, this assumption simply entails a scaling of time. Still, it allows us to take $\Delta=1$ and it ensures that

[TABLE]

is a probability matrix for all $\varepsilon\in[0,1]$, which is convenient. We remind the reader that our ultimate goal is to find (or approximate) the stationary measure belonging to $\bm{G}(1)$ (and, equivalently, of the discrete time counter part $\bm{P}(1)$). In order to achieve that, we first establish the analyticity of the stationary distribution for $\varepsilon$ in a punctured neighborhood of 0, cf. Theorem 7.1 below. We emphasize that it is not guaranteed that the stationary distribution will be analytic up to $\varepsilon=1$. The analysis in [3] gives a lower bound for the radius of convergence, which in general turns out to be rather conservative.

Note that the perturbed transition probability matrix $\bm{P}(\varepsilon)$ can also be decomposed into the unperturbed probability matrix $\bm{P}^{(0)}$ and the perturbation matrix $\bm{G}^{(1)}$, with $\bm{P}^{(0)}=\bm{I}+\bm{G}^{(0)}$, i.e.,

[TABLE]

It is evident that the unperturbed process consists of several ergodic classes, making our setting fit the singular perturbation approach in [3].

7.1 Singular perturbation analysis: outline

Following the analysis performed in [3], we now formulate four conditions based on which the invariant probability measure of the perturbed Markov chain, denoted by $\pi(\varepsilon)$, is derived. These four conditions guarantee the analyticity of $\pi(\varepsilon)$ in $\varepsilon$ in a punctured neighborhood of zero. Furthermore, they guarantee that the coefficients of the power series $\pi(\varepsilon)=\sum_{m=0}^{\infty}\varepsilon^{m}\pi^{(m)}$ form a geometric sequence and, hence, that there exists a computationally stable updating formula for $\pi(\varepsilon)$, see [3].

In this subsection we only formulate the four conditions and give the main result of the section. The detailed mathematical proofs follow in the next subsection.

Condition 7.1**.**

The unperturbed Markov chain consists of several (denumerable) ergodic classes and there are no transient states.

There is an ergodic class for each \bm{i}\in\big{\{}(n_{1},n_{2}),\ n_{1},n_{2}\in\mathbb{N}\big{\}}, i.e., in an ergodic class, the numbers of customers in both queues are fixed. All ergodic classes are identical, and consist of two states, $k\in\{1,2\}$, indicating the queue being served.

Condition 7.2**.**

The Markov chains corresponding to the ergodic classes of the unperturbed Markov chain are uniformly Lyapunov stable i.e., for each ergodic class there exist a strongly aperiodic state $\alpha\in\{1,2\}$ (with a strictly positive probability on the corresponding diagonal element of the transition matrix $\bm{I}+\bm{C}$, with the matrix $\bm{C}$ given in (7.2)), constants $0<\delta<1$ and $b<\infty$, and a Lyapunov function \bm{u}=(\begin{array}[]{c c}u_{1}&u_{2}\end{array})^{\prime}, with $u_{i}\geq 1$, $i=1,2$, such that

[TABLE]

with $\bm{e}_{\alpha}$ a vector with 1 in the entry belonging to state $\alpha$ and zero in the other entry.

For the next condition, we first introduce the aggregated Markov chain [20, 22, 36], with generator $\bm{\Gamma}$, given in matrix form as follows

[TABLE]

with $\bm{V}$ and $\bm{W}$ defined as in [3, p. 844]; $\bm{V}$ (resp., $\bm{W}$) is a matrix whose rows (columns) correspond to the ergodic classes and its columns (rows) to the states in $\mathcal{S}$. The $\bm{i}$-th row of $\bm{V}$ is the invariant measure of the unperturbed Markov chain, given that the process starts in the $\bm{i}$-th ergodic class, with \bm{i}\in\big{\{}(n_{1},n_{2}),\ n_{1},n_{2}\in\mathbb{N}\big{\}}, i.e.,

[TABLE]

with $\bm{\tilde{C}}=\begin{bmatrix}c_{2}/(c_{1}+c_{2})&c_{1}/(c_{1}+c_{2})\\ \end{bmatrix}$. The $\bm{j}$-th column of $\bm{W}$ has ones in the components corresponding to the $\bm{j}$-th ergodic class and zeros in the other components, with \bm{j}\in\big{\{}(n_{1},n_{2}),\ n_{1},n_{2}\in\mathbb{N}\big{\}}, i.e.,

[TABLE]

with $\bm{1_{2\times 1}}=\begin{bmatrix}1\\ 1\end{bmatrix}$.

Hence, for $n_{1}\geq 0$, $n_{2}\geq 0$, the elements of the generator matrix $\bm{\Gamma}$ are:

[TABLE]

It is convenient to think of the aggregated Markov chain as the limiting joint queue length process as $\varepsilon\to 0$. In this limit, the server moves infinitely fast between the two queues, making them two independent M/M/1 queues with arrival rates $\lambda_{i}$ and service rates ${\mu_{i}}\frac{c_{1}c_{2}/{c_{i}}}{c_{1}+c_{2}}$, $i=1,2$. Based on this remark, one can immediately deduce that the invariant probability measure of the aggregated Markov chain is

[TABLE]

with $\tilde{\rho}_{i}=\frac{\lambda_{i}c_{i}(c_{1}+c_{2})}{\mu_{i}c_{1}c_{2}}$, $i=1,2$.

We are now ready to state the third condition.

Condition 7.3**.**

The aggregated Markov chain is irreducible and Lyapunov stable, i.e., there exist a strongly aperiodic state $\bar{\alpha}=(n_{1},n_{2})$ (with a strictly positive probability on the diagonal of the transition matrix $\bm{I}+\bm{\Gamma}$, with the matrix $\bm{\Gamma}$ given in (7.11)), constants $0<\overline{\delta}<1$, $\overline{b}<\infty$ and a Lyapunov function $\overline{\bm{u}}=\left(\overline{u}_{(n_{1},n_{2})}\right)_{(n_{1},n_{2})\in\mathbb{N}^{2}}$, with elements $\overline{u}_{(n_{1},n_{2})}\geq 1$, for all $n_{1},n_{2}\geq 0$, such that

[TABLE]

Condition 7.4**.**

The perturbation matrix $\bm{G}^{(1)}$ is $\tilde{\bm{u}}$-bounded (for $\tilde{u}_{\bm{i}k}=\overline{u}_{\bm{i}}u_{k}$, with \bm{i}\in\big{\{}(n_{1},n_{2}),\ n_{1},n_{2}\in\mathbb{N}\big{\}} and $k=1,2$) or, equivalently,

[TABLE]

is bounded by some constant $g>0$, cf. [3, p. 841].

Note that, because of the repetitive structure of $\bm{G}^{(0)}$, this assumption implies that $\bm{P}(\varepsilon)$ is $\tilde{\bm{u}}$-bounded for all $\varepsilon\geq 0$.

We can now state the main theorem of the section, which is based on [3, p. 845, Theorem 4.1].

Theorem 7.1**.**

Under Conditions 7.1–7.4, the perturbed Markov chain has a unique invariant probability measure, $\pi(\varepsilon)$, which is an analytic function of $\varepsilon$ in a neighborhood of 0,

[TABLE]

where $\bar{\pi}$ is the invariant probability measure of the aggregated Markov chain, cf. (7.12), and

[TABLE]

$\bm{V}$* and $\bm{W}$ are given in (7.9) and (7.10), respectively, and $\bm{H}$ and $\bm{\Phi}$ the deviation matrices of the unperturbed and aggregated Markov chains, respectively, are given by*

[TABLE]

and

[TABLE]

Here $\bm{\gamma}$ is the ergodic projection of the aggregated Markov chain, with generator $\bm{\Gamma}$ given in (7.11), i.e.,

[TABLE]

Remark 7.1**.**

We do not discuss the radius of convergence of the series in (7.15). Theorem 4.1 of [3] gives a (rather conservative) bound for the analyticity region.

Remark 7.2**.**

The invariant probability measure of the perturbed Markov chain can be calculated by the updating formula

[TABLE]

with $\varepsilon$ in a neighborhood of 0, cf. [3, p. 845, Theorem 4.1].

Remark 7.3**.**

In order to calculate the deviation matrix $\bm{\Phi}$, one may use the following equations

[TABLE]

We briefly describe two approaches to obtain the deviation matrix $\bm{\Phi}$: an analytic one involving PGFs and a numerical one. Both approaches require some additional work. The analytic approach, which involves the consideration of generating functions, leads to a boundary value problem for which we can employ Steps 1 and 2 discussed in Section 6. Performing these steps reveals a problem similar to the combinatorial random walk in the quadrant with transitions to the West, North, and South-East, cf. [9, Section 5.2]. In order to obtain the expression for $\bm{\Phi}$, we need to invert the obtained PGF. A numerical approach is to truncate the state space and solve numerically the corresponding finite system of equations above. We do remark that truncating the state space is a delicate task, since the entries of $\bm{\Phi}$ corresponding to states far from the origin are unbounded. We do not further investigate this in this paper.

7.2 Singular perturbation analysis: verification of the conditions

It remains to prove that Conditions 7.1 - 7.4 are satisfied and also to indicate how the deviation matrix of the unperturbed Markov chain, $\bm{H}$, is calculated.

Verification of Condition 7.1

As explained in the previous section, this condition follows directly from Equation (7.6).

Verification of Condition 7.2

Obviously, all ergodic classes are identical and contain two states ($k=1,2$), thus this condition is trivially met, but for the construction in the remainder it is useful to specify the Lyapunov function used.

First note that the strong aperiodicity follows from the uniformization condition (7.4). We can choose any of the two states as the strongly aperiodic state; in the following we choose $\alpha:=1$. To construct the Lyapunov function first we choose the constants $\delta$ and $b$ as $\delta\in\left(1-c_{2},1-\frac{c_{1}c_{2}}{c_{1}+c_{2}}\right)$, $b=1-\delta+\frac{c_{1}^{2}}{c_{2}}$. Then we can verify that the Lyapunov function

[TABLE]

satisfies (7.7). It also follows that, indeed, $\delta\in(0,1)$, $0<b<\infty$ and $u_{k}\geq 1$, $k=1,2$.

Verification of Condition 7.3

From the definition of the generator of the aggregated Markov chain, cf. (7.11), and the stability condition (2.1), it is immediately evident that the aggregated Markov chain is ergodic, since it behaves as two independent ergodic M/M/1 queues with arrival rate $\lambda_{i}$ and service rate $\frac{\mu_{i}}{c_{i}}\frac{c_{1}c_{2}}{c_{1}+c_{2}}$, $i=1,2$.

Now by using the uniformization condition (7.4), state $(0,0)$ is strongly aperiodic i.e., we may choose $\bar{\alpha}=(0,0)$. We proceed to describe the Lyapunov function $\overline{\bm{u}}$ and the constants $\bar{\delta}\in(0,1)$ and $\bar{b}$ which satisfy Condition 7.3. Note that Equation (7.13) is written as follows, for $n_{1},n_{2}\geq 0$,

[TABLE]

Solving the above equations with equality, after choosing

[TABLE]

yields the solution for $\bar{\delta}$ and $\bar{b}$. We choose

[TABLE]

and

[TABLE]

Note that due to the uniformization condition (7.4), indeed $\bar{\delta}\in(0,1)$, $0<\bar{b}<\infty$ and $\overline{u}_{(n_{1},n_{2})}\geq 1,$ for all $n_{1},n_{2}\geq 0.$

Verification of condition 7.4

To verify this assumption, we apply the definition, cf. (7.14), and show that

[TABLE]

with $g_{1}=\left(\frac{\mu_{1}c_{2}}{\lambda_{1}(c_{1}+c_{2})}\right)^{-\frac{1}{2}}\left(\mu_{1}+\frac{\mu_{1}c_{2}+\mu_{2}c_{1}}{c_{1}+c_{2}}\right)$ and $g_{2}=\left(\frac{\mu_{2}c_{1}}{\lambda_{2}(c_{1}+c_{2})}\right)^{-\frac{1}{2}}\left(\mu_{2}+\frac{\mu_{1}c_{2}+\mu_{2}c_{1}}{c_{1}+c_{2}}\right)$.

In order to do so, we use the following $\tilde{\bm{u}}$-norm

[TABLE]

with $\overline{u}_{(n_{1},n_{2})}$ given in (7.21) and $u_{k}$ given in (7.20).

Derivation of the deviation matrix of the unperturbed Markov chain

It follows from Condition 7.1, that the deviation matrix of the unperturbed Markov chain, $\bm{H}$, has the following block diagonal structure

[TABLE]

with $\bm{H}_{2\times 2}$ the deviation matrix of each ergodic class of the unperturbed Markov chain, i.e.,

[TABLE]

with $\bm{C}$ given in (7.2) and $\bm{c}$ the ergodic projection of the unperturbed Markov chain given as

[TABLE]

cf. [39, p. 64, Equation 4.1].

We evaluate (7.23) using the spectral decomposition (eigen-decomposition) of matrices $\bm{I}+\bm{C}$ and $\bm{c}$; the diagonal matrices containing the eigenvalues are $\bm{D}_{\bm{I}+\bm{C}}=\mathrm{diag}\{1,1-(c_{1}+c_{2})\}=\begin{bmatrix}1&0\\ 0&1-(c_{1}+c_{2})\end{bmatrix}$ and $\bm{D}_{\bm{c}}=\mathrm{diag}\{1,0\}=\begin{bmatrix}1&0\\ 0&0\end{bmatrix}$, respectively, and the corresponding matrix of eigenvectors is $\bm{M}=\begin{bmatrix}1&-c_{1}\\ 1&c_{2}\end{bmatrix}$. Naturally, in dimension 2, both matrices produce the same eigenvectors because $\bm{c}$ is the ergodic projection of $\bm{I}+\bm{C}$. Therefore, Equation (7.23) can be written as

[TABLE]

Combining (7.24) and (7.22) yields Equation (7.17).

8 Possible future directions

We have studied a single server two-queue polling model with a random residing time service discipline. More concretely, we considered that customers arrive at the two queues according to two independent Poisson processes. There is a single server that serves both queues with generally distributed service times. The server spends an exponentially distributed amount of time in each queue. After the completion of this residing time, the server instantaneously switches to the other queue, i.e., there is no switch-over time. A service discipline with a random residing time does not satisfy the so-called branching property [37], which significantly complicates the underlying analysis.

For this polling model, we derived the steady-state marginal workload distribution and used it to obtain several asymptotic results. We also discussed the complications arising in the calculation of the joint workload distribution. Furthermore, restricting ourselves to the case of exponential service times, we have calculated the joint queue length distribution using (singular) perturbation analysis. It is a topic for further research to determine how to efficiently truncate the system without inducing too large errors. The insights gained for the two-queue polling model, specifically for the derivation of the marginal workload, cf. Section 3, can be also used in the case of $N$ queues, $N>2$. In addition, one may generalize the compound Poisson input to a Lévy subordinator input process.

Another interesting topic for future research is to develop the framework for the derivation of the bivariate LST of the joint workload distribution, in particular in the asymmetric case, cf. Section 6, and for the derivation of the bivariate PGF of the joint queue length distribution in the case of exponential service requirements. In particular, the objective in such a setting is to develop an approach for the transformation of Equation (6.4), and its version for the asymmetric case, into a Riemann or Riemann-Hilbert boundary value problem. This requires, that we first choose the zeros of the kernel equation $K(s_{1},s_{2})$, so as to define a closed smooth contour. Thereafter, we need to show that Equation (6.4) on the contour reduces to the study of an analytic function (probably with the exception of one pole) with a known boundary condition. An interesting alternative direction would be to extend the framework developed by Fayolle et al. [24], of the systematic use of the kernel method using the group of birational transformations that leave the kernel equation unchanged. The challenge in our case is that the kernel $K(s_{1},s_{2})$ does not have the regular structure indicated in [24], but this does not seem to impose an insuperable obstacle, see, e.g., [29].

Acknowledgments

The authors gratefully acknowledge useful discussions with Ahmad Al Hanbali, Offer Kella, Socrates Olympios, Shelly Zacks and Bert Zwart. The research of Mayank Saxena was funded by the NWO TOP-C1 project of the Netherlands Organisation for Scientific Research. The research of Onno Boxma was partly done in the framework of the IAP BESTCOM project, funded by the Belgian government; it was also funded by the NWO Gravitation Program NETWORKS of the Dutch government. The work of Stella Kapodistria is supported by the NWO Gravitation Program NETWORKS of the Dutch government.