WAITING TIME VARIANCES IN CONTINUOUS-TIME POLLING SYSTEMS Mandyam M. Srinivasan Department of Industrial & Operations Engineering University of Michigan Ann Arbor, MI 48109-2117 Technical Report 91-37 December 1991

Waiting Time Variances in Continuous-Time Polling Systems Mandyam M. Srinivasan Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, MI 48109-2117 Technical Report 91-37 December 1991 Abstract We present an efficient technique to compute waiting time variances in polling systems that use the exhaustive and gated service disciplines. The polling systems can have a mix of service disciplines, wherein some of the nodes can use the exhaustive service discipline while the other nodes use the gated service discipline. The technique allows the user to obtain the waiting time variances at select nodes without having to obtain waiting times variances at all nodes simultaneously. Once the waiting time variances have been evaluated for a given configuration, then the IS technique allows the user to rapidly evaluate several alternate configurations having different switchover time parameters, with just a few elementary arithmetic operations. Keywords: Cyclic Queues; Polling Systems; Algorithms; Waiting Time Variances

1. Introduction A polling system is a queueing system in which a single server services M nodes. The server cycles around the system, visiting each node in turn, and attends to the waiting customers using some prespecified service discipline. Real-world applications which are modeled and analyzed using polling systems can be found, for example, in computer communication and material handling systems, and there is a vast body of literature devoted to the analysis of polling systems. An excellent survey of the research done on polling systems can be found in Takagi [10]. A number of service disciplines have been studied, which determine the service strategy adopted by the server at a node. Two of the most commonly studied service disciplines are i) the exhaustive service discipline, where the server continues to serve all customers at the node until there are no more customers at the node, ii) the gated service discipline, where the server only serves those customers present at the node at the instant the node is polled. A number of techniques that have been developed for computing the mean waiting times in polling systems that use the exhaustive and gated service disciplines. However, there is very little work on computing the variances of these waiting times. To our knowledge, there are two papers by Aminetzah [1] and Ferguson [5] for computing waiting time variances. These algorithms are based on solving a set of O(M3) equations. As observed in Takagi [10], these equations are quite complicated and require the use of a symbolic formula manipulation language. No computational results are presented. In a recent paper [8], we developed a powerful technique, termed the Individual Station (IS) technique, to compute mean waiting times for continuous-time polling systems using the exhaustive and/or gated service disciplines. In this paper, we further develop the IS technique and present an iterative algorithm to compute the waiting time variances for these systems. Although the iterative algorithm apparently requires the computation of an infinite sequence of terms, we can show that the number of iterations it requires is logarithmic in the accuracy required, using arguments similar to those presented by Levy [7]. The IS technique allows us to analyze systems with mixed service disciplines, wherein some nodes are served using the exhaustive service discipline while the other nodes are served using the gatedl service discipline. The technique also allows the user to determine the waiting time variances at one or more select nodes without having to obtain the waiting time variances at all the nodes, simultaneously. Furthermore, once waiting time variances are evaluated for a configuration, then the user can rapidly evaluate several alternate configurations that have different values for the switchover time parameters, with just a few elementary arithmetic operations. 1

The notation involved in the intermediate stages during the development of the algorithm may seem a little complex. However, the reader will find that very little notation in involved in the final expressions for the mean and variance of waiting times, and that the resulting algorithm is very simple and easily implemented. We provide computational results to demonstrate the efficiency of the algorithm. 2. The Waiting Time Moments The waiting time at node m is the time spent by an arbitrary customer, waiting for service, and is denoted by the random variable Wm. Customers are served according to the first-come-firstserved rule. In this section, we derive expressions for the first and second moments of Wm, denoted by E(Wm) and E(W2). This paper is concerned with computing the second moment/variance of the waiting time. Algorithms for computing the first moment of the waiting time have already been presented in an earlier paper [8]. We first present the notation. Unless otherwise specified, it is implicit that the index for any summation is over the range 1 through M. It is also implicit that the index is i) reset to 1 if it becomes M+1, and ii) reset to M if it drops down to 0. We adopt the convention that an empty product equals 1. A single server serves requests from customers at the M nodes. Customers arrive at each node according to independent Poisson processes with rate Xm at node m. The service time, Bm, of a customer at node m is an independent random variable, with its ith moment denoted by E(B ). The Laplace Stieltjes Transform (LST) of the service time distribution is denoted by 3m(.). The busy period, Hm, induced by a single customer at node m has its ith moment denoted by E(HI), and its LST denoted by Tim(.). The first three moments of Hm are (Kleinrock [6]): E(BE(B2m) E(B3 Xm(EB 2m))2 E(Hm) - E(Hm) = E E(Hm) = m 4 + 3 (3.1) 1-Pm ' (1-pm)3 ' (1-pm)4 (1-pm)5 The time taken by the server to switch from node m to node m+l is an independent random variable, Sm, with its ith moment denoted by E(S ). The variance and LST of Sm are denoted by Var(Sm), and am(.), respectively. The traffic intensity at node m is denoted by Pm = XmE(Bm), and p = Xm pm denotes the server utilization. For the polling system to be stable, we require p to be less than one. If the polling system is stable, then it can be shown that the cycle time, C, which is the expected time between two successive visits by the server to a node is the same for all nodes, and is given by C = m E(Sm)/(l-p). We are interested in obtaining the waiting time variance, Var(Wm), at node m, for m=l,...,M. Let ~E denote the set of nodes that is served using the 2

exhaustive service discipline, and let av denote the set of nodes that is served using the gated service discipline. Let W*(s) denote the LST of the waiting time distribution. Let Fm(zi,..., zM) denote the probability generating function (p.g.f.) for the number of customers at each node at the instant the serve-r polls node m, and for notational ease, let Fm(xm) = Fm(zl,..., Zm4, Xm, zm+n,1,~., zM). For polling systems that use the exhaustive or gated service discipline, Wm(s) is expressed in terms of Fm(xm) as follows (see, for instance, Takagi [9]): W *(S)1 lFm(1-s/Xm)ME32a W,(s) = E mE, (3.2a) C~ s - km + XMP,(s) W *(S) = 1 Fm(p(s))- Fm(1-s/Xm) C SXm+XmPm(S) M. (3.2b) Let fQ) denote the ith factorial moment of the number of customers present at node m when the server polls node m. These terms are obtained as: aiFm(zm) f~i> 1azmi (3.3) m = 'z=1 From equations (3.2a) and (3.2b), we obtain E(Wm) = -DW (S)/asI,=o as (this requires two applications of L'Hospital's rule) f(2)m kmE(B2) E(W = m + 2(1-Pm)' mEE, (3.4a) 2Xmfm f(2) E(Wi) = m (l+PM), inCa. (3.4b) 2Xmfm We obtain E(W,) = W (/ as (this requires three applications of L'Hospital's rule) ' E 2 f f3) ~( E(Wm) = E(Wm) + _+ meE, (3.5a) )(l-Pm) 3 X2 fm 3(1-pm) XmE(Bm) fIM) E(Wm2) = E(Wm) + - (l + PM + PM M EinCa. (3.5b) Ol+pM) 3 X2 fm T'herefore, if we can obtain the factorial moments, f~, f(, and f(gm then we can determine E(W,,,), and Var(Wm) = E(Wm2) - (E(Wm))2. However, obtaining these factorial moments is not straightforward, since the generating function Fm(.) is expressed, recursively, in terms of the 3

generating function FM-,1(). It is well known (see for instance, Takagi [9]) that the generating functions are related by the following expressions: Fm+i(zli,...,zM) = Fm('lM[Xk#M (Xk-XkkzlI) aYm(Xk (Lk-Xkzk))j, ME e E (3.6a) Fm+i(zl,....,zM,) = Fm(P3m[Xk (Xk-XLkZk)]) Om(Xk (Xk-kkzk)),) MEa. (3.6b) Equations (3.6a) and (3.6b) form the starting point for the IS technique, which obtains an expression for Fm(.) solely in terms of functions of input parameters. We can now differentiate the resulting expression for Fm(.) three times to obtain the three factorial moments, and then use equations (3.4a) through (3.5b) to obtain the first and second moments of the waiting times. The technique we develop below computes E(WI) and E(W21). The first and second moments of the waiting times at the other nodes are obtained in an identical manner. We first define a recursive (nested) function, ymo), as follows: 'ym(O) = Zm, m=l,.,M, (3.7a) 7m() TIrM( X~k-XkYk(}4)]+ [Xk-XkYkOj)]), j>O, MEtr, (3.7b) k<m k>m 'YmO) PMX [Xk-XkAU-l)] + [ Rk-XYkO)]), j >0, M E a. (3.7c) k__ k>m Let amo) = O(TM [Xk-kYkOj)]+ [Xk-XkYkOj+1)]) j~~~0, M=l,...,)M. (3.8) k~m k>m W e can now rewrite both equations (3.6a) and (3.6b) in the following general form: Fj(y1 (0),...,ym(O)) = FM(y1j(0),...,yNi1(O),yM(l)) caM(O). Recursively expressing FM(.) in terms of FM-,1(.), and FM-,(.) in terms of FM-2(., and so on, we obtain: Fj(y1(O)j...5,ym(0)) = Fi(yj(l)9...,,ym(l))flGm(O). We now continue to recursively express Fi(yfj),...,ym(j)) in terms of Fj(y(jO+1),....,YMj+ 1)),' for j ~ 1, and this results in the following expression for F, (y1 (O),... O n-i Fj(y1(O)9...,yM(0)) = Fj(y1(n),,...,,yM(n)) 17a(7 or0)) (3.9) 4

It can be shown (refer Cooper [2], and Eisenberg [4], for example), that limn....,yi(n) = 1. Hence, letting n ->00o in equation (3.9), we obtain: F 1Q(yj(0), I.. gym (O)) = 00 am ') - 0 j —om (3.10) We now have an equation for Fl(.) that is expressed only in terms of (complex) functions of the cYm() values. Let z = (zl,..., zM), and 1 = (,.,).Differentiating equation (3.10) once, twice, and thrice with respect to zi, and setting z = 1, we obtain the following expressions for the factorial moments of Fl(.): fi a= Iz=19 (3.l1la) f(2). =fi2+X j=O m ~2a(j) - _mj)_2 and (3.1llb) j=O (3arn - aim) az2mo + 2amo 3\".(3.ll 1c) -ZI 3akj az2amj __zl IZ In the above equations, for the sake of clarity, we let f1 denote the factorial moment, f('), and will continue to do so henceforth. Equations (3. 1 l a) through (3. 1 i c) hold regardless of whether node 1 is served according to the exhaustive service discipline, or the gated service discipline. 4. Expressions for the Factorial Moments of Fm(zm) F~rom equation (3.8), we note that the derivatives in equations (3.1 la) through (3. 1 lc) can be obtained in terms of aiym/azji1, for i = 1,...,) 3. Thus, if we define a1i and (4.1) Xwo() = IEX T)(j+ Nf j+I) ], j~>0,) k~m k>m (4.2) then. the derivatives in equations (3. 1 l a) through (3. 1 i c) are expressed, for j ~ 0, as: I z=i = E(Sm) X(M)), a20m I 1 E(Sm) X(2(j) + E(S 2) (X~(1j))2, (4.3a) (4.3b) az1. lz1= = E(S3) (X)(1)0)3 + 3 E(S2) X~(1j) X(2)(j) + E(Sm) X(3Cj). 5

Note that equations (4.3a) and (4.3b) are valid for both m GE ~, and m e: a-. We can obtain a recursive expression for V4(ji~), i = 1, 2, 3, for j > 0, by differentiating equations (3.7b) and (3.7c) with respect to zi the appropriate number of times, and setting z equal to 1. The expressions for these terms will depend on whether m e= t, or m = ax. We will first obtain the recursive expressions for m E E. Let YT() = I AV ~)() + V wP0+1) I= X(0) - y(M(j), j~O. (4.4) k<m k>m Analogous to equations (4.3a) and (4.3b) for the derivatives of Gmoj) evaluated at z = 1, recursive expressions for Vi44o(), m E E, and j > 0, are obtained as N~))= XmnE(Hm) Y(0(-1), VI(2)() = km(E(Hm)y(2)(0V) + E(H2j) (y(1)(-1))2) (4.5a) 14i(1)() = Xm(E(HI) (y(1)(j-1))3 + 3 E(H2j) y(1(j-1) y(2(j-1) + E(Hm) y~(3j-1)) (4.5b) For m E a, and j > 0, the corresponding expressions for xp(i(j) are: Om~)= XME(Bm) X(1)n?}-1)3, Oji(2)() = km(E(Bm)X(2)(i1) + E(B2) (X(1(j-1))2) (4.6a) Vm))= Xm(E(B3) (X(1)(-1))3 + 3 E(B 2) X(1)(-1) X(2)(-1) + E(Bm) X~(3jl)). (4.6b) Let Pm Pm ~~~~~~~~~XE(B2) -M 'Am, and 1 2 (4.7) l-Pm l+Pm PM We first account for all the unknowns, Vj,(1)(), j > 0, using the following observation. Observation 1: a) The term 01()(O) = Xkk~ 'z1 = XI when k=1; it is equal to zero otherwise. b) The term Vqi()(O) = 0, i > 0, for all k. (Refer to equation 4. 1.)U Based on Observation 1, we obtain (Srinivasan [8]): Nf())= iI ' I' V (1)0i\=, ier=a, and (4.8a) m j=O ipm j=Oi20 = = (2)(1)(j))2 (4.8b) m j=O m j=1 m j= 1p 6

We obtain an expression for m j~~=l V t3)(j) using an approach similar to the approach used to obtain the above equations, namely, (4.8a) and (4.8b). First, note that for m E E;, from equations (4.4) and (4.5a), Y(l)(j-1) = W)(j)/mE(Hm). We can use equations (3.1) and (4.7) to rewrite this as wjl)(j) Yl)(j-1) = m, me~, j>0. (4.9) Om Again, from equations (4.4) and (4.5a), for m e E, and j > 0, (2) lm2)() - mE (Hm2) m( -) (j) Xm2E(Hm) V) 2)0) m-1) u ="Y m --- - Tn 'l1c" m M. (4 10) Y(((1 I XmE(B-) 1 - am Om Om equation (3.1), we arrive at the following intermediate result after some elementary algebra:(11) (m)(i) = (4>((j)) - 3 1-r + )(j) )(j) +6mY()(j-). (4. Pm(-Pm) 1Pm -Pm (4.11) From Observation 1, (3m)(0) = 0, for all m, and hence,?J Y(3)(-1) = Xj Xk Vm N.)J) Therefore, if we add 5m (3)(j) to both sides of equation (4.11), and sum both sides of the equation over j=l,...,oo, we obtain: OC f (3)(j) (l+S1m) j=l ( (1)(J))3 (rmE(B3) +m ((1))() 1-pm ( p ~ m -302 +3 ~ Ev lit 11t(2)( J) + j=l 1-P m Pj= j=l k We now sum both sides of the above equation over all m, and collect terms, to get (note that ( 1+m) = 1/(1-m)): to yM XIE(B 3 15m2 )\00 1~m 00 ~j (3m)) = - 3 (-(l)(j))3 + 3 m V ()0)2)). (4.12) m j=l m (1p)p3 1-p)j= L m 1-p j=l Following an identical approach, we can show that equation (4.12) holds even for m E i. From equations (3.1 la), (4.2), and (4.3a), noting that (1)(0) = 0 for k > 1 we get: 7

00 00 fl= -] E(Sm) [ I())+ ()(j+l)] E(Sm) v)() m j=O k<m k>m m j=O k Therefore, from equation (4.8a), fi is obtained as: fl = xl(l-pl) C, 1 E, fl = C, 1 E. (4.13) In a similar manner, from equations (3.1 lb), (4.2), (4.3a), and (4.8b), we get: 00 00 f(2) = f2 + C 61m ((lm)(j)2 + + Var(Sm) (X()(j))2. (4.14) m j=l m j=0 From equations (3.1 1c), (4.2), (4.3b), and (4.12), the expression for f(3) is obtained as: 00 00 f3) = 3ff(2) - 2fi3 + I T(Sm) X (X ()(j))3 + E Var(Sm) I X (l(j) X(2)(j) m j=0 m j=0 _ m'm)4nE m _ V + CE 3 3 2 (Vl)(1i))3 + 3 C m ()(j) t (2)(j) (4.15) m Pm j= m j=1 where T(Sm) denotes the third central moment of the random variable Sm. (Note that Var(Sm) is the second central moment of Sm.) Equations (4.14) and (4.15) are valid regardless of whether node 1 is served using the exhaustive service discipline or using the gated service discipline. Once we compute the three factorial moments, the mean waiting time and the second moment of the waiting time at node 1 can be obtained using equations (3.4) and (3.5). While the first factorial moment is directly computed, computing f(2) and f(3) apparently require the computation of infinite sums. However, we can show, using arguments similar to those presented by Levy [7], that the number of iterations required to compute these moments is logarithmic in the accuracy required. Let pmin = minm(pm). Then the overall complexity of the algorithm to determine the variance of the waiting time at one station can be shown to be O(M log~e), where e is the accuracy required, and gt = [(p - pmin)/(l - Pmin)]2. In section 5, we present the algorithm to compute the three factorial moments, and thereby the variance of the waiting time. In section 6, we present a number of computational results, and report on the execution time of the algorithm. In section 6, we also present a closed form expression for f(2) and f(3), for the special case of a pure exhaustive service system with two nodes. 8

5. The Algorithm In order to compute the factorial moments using equations (4.13) through (4.15), we must evauat te trm 14(V(j),) Nf~(2)0) X(1)(j), and X(2~)(), M = 1,... M, iteratively, starting at j = 0. It is convenient, at this stage, to define the terms, XmoJ) and X(2), as follows: xi xm cymoj) Xl az and xmi = X12 (5.1) Inste-.ad of computing WN$Ok) and 'ij2)j), we compute Xmoi) and X(2kj), starting with Xi(O) = 1, Xm(O) = 0 for m ~ 1, and. X~(2j) = 0 for all m. From equations (4. 1) through (4.6) we get: 0m4i x4nmj) = iIXkl + XOI k<m k>m = PM [ XXk(J-l)+ X~kJ)],9 k~m k>m mec-t, j >0, inC a, j >O. (5.2a) (5.2b) From equations (4. 1) through (4.6), we also get k<m k>m PM xM. k~m Xk k>m xT +PM M), MEtr, j >O. inca, j>0. (5. 3a) (5.3b) Equatos(5.2) and (5.3) suggest that the computational effor tootiimJ nd X( j)fo given j requires 0(M) operations. However, we can express these terms as functions of neighboring terms with 0(1) arithmetic operations, independent of M. The computational formula we derive will, however, depend on the service discipline used for node m+1. We derive the computational formula for Xm) when m E and M+1 E v. The expressions for the other cases are obtained in a similar manner. We consider a node m ~ M. From equation (5.2a) applied at m and equation (5.2b) applied at m+1, we obtain the following relationship between Xmoj) and xm+PlJ) when mE tandm+l E a: Xm=i) = 6m II XkCJ-l) + XkO) + Xm+1(J) - XmOhl) - Xm+1i(}+] k~m+l k>m+l = m I xrn+ij)/Prn+i + Xm+iOj) - Xm(.V1) - Xm+1(hl)]. = 8m [Xm+ilj)/Am+l - Xmohl) - Xm+i(WD]1. 9

The expressions for all the possible combinations are obtained in a similar manner, and presented below as equations (5.5a) through (5.9). We use the convention Xm+l(k) = Xl(k-1) if m = M. Similarly, Xm2+)(k) = %(2)(k-1) if m = M. There is one special case to be considered; the case where m = M, and j = 1. For this special case, we have: XM(l) = 8M x^(1) = PM X2(1), PM Me E; (5.5a) XM(1) = PM, X(1l) = PM M2(1), PM Me a. (5.5b) In general, we have the following expressions. Let: Y (2) ) (im+ Am+(j) Pmi - Pn 2 m )+ Pm2m i Am+l(j Pm Pm+1 m+ PM Xm = x2+)l() m+l 2 m+ Bm+lo() - AI- PI 10)+ -m' Am~l Pm Case i) me E, m+l e L: Xm(i) = 6m [Xm )- Xm(-1)], Pm+1 (2)(j) = m [Am+l(j) -X(2)-1)] m m~m (5.6) Caseii) me E, m+le a: = Xm+i(i) Xno) = Sm[A - m(o4>) - Xm+i(}4)], ~Lmr+l (m (ji-) = B) - X (j-1) - -, (5.7) Case iii) mea G, m+ 1 E: Xm(j) = P [Xm+l() ] Pm+l X(2)j) = Pm[Am+l(j)], (5.8) Case iv) me G, m+l e G: Xm(j) = Pm [Xm+l(J) _ m+lj-1)], Am+l Once we compute Xm4(jl) and (2)(j+1), From equations (4.4), (4.9) and (5.1), (2)(j) = Pm Bm+l j) -X (2)l(-1) m m [Bm~l V) -m+ (5.9) we quickly compute X(l)(j), and X(2)(j) as follows.,,y\1lyuJ~vrrr~r -M M/ /-1 ---~~ X()(j) = X Xm(j) + Y()() = 1 [ Xm(j) + m+l ], mv m x [xo + rn1 me E, (5.10a) (5.10b) (2)(i+l) 2 X(2)(j) =2 x(2)(i)+ Y(2m)(j) =- 2 [ X(2mj) + - 2m(+l)] me. ~ m Pm 10

Similarly, from equations (4.2), (4.6a), and (5. 1), we obtain X(l)(j) = 1XI mO+l), m E. (5.l la) PM X(2kJ+1) ~ X(2)(j) = X2[xm 2 ~j+),Mr=. (5lb M 1 ~PM PM mo The number of arithmetic operations required to compute these terms is thus independent of the number of nodes in the system. Hence, for a given j, computing these terms for all m requires 0(M) operations. Let 00 00 00 (PM Om = XXn() COin= Xm(J)X(2k)0, (5.12a) j=1 j=1 j=1 00 X ~~~ ~~00 00 (X()Xm__2 j= K1 = M0 ' = Mi I Inn - XmuE(B~) (.3 P3~ The expressions for fp2) and f~3) are now rewritten as: f~2) f 2+ X2(c (P 'Var(Sm) ~m (5.14a) M ~~~~m) =~3 3 f, P2 -2f 3 + 3 (T(Sm)am +3Var(Sm) Km + (((7cm- 312)em + 3lmCom)). (5. 14b) The iterative algorithm to compute the waiting time variances is given below as Algorithm 1. Algorithm 1. 0. Initialize values: Set a) j = 0, b) X1O) = 1, and yX(j) = 0, m = 2,...,) M, c) X1~(2j) = 0, for allm, and d) (pm= Om = om=4m =cam = Km =Ofor allmr. Choose atolerance valuecF. 1. Setj =j+1, and for each ma) compute Xm() and X(2~) )using equations (5.5) through (5.10), b) compute XV(j) and X(2(j) using equations (5.l1Oa.) through (5.1 ib). 2. Update the values of (pm, Om, Com, ~m am.) and Kin, for m = 1,...,) M, using equations (5. 12a) and (5.12b). 11I

3. If Xm() < e, for all m, go to step 4. Otherwise, go to step 1. 4. Compute f(2) and f(3) using equations (5.14a) and (5.14b), and compute E(Wi) and E(W2) using equations (3.4a) through (3.5b). Obtain Var(W1) = E(W2) - (E(W1))2. Algorithm 1 is executed for each node whose waiting time variance is desired. Each time, a new set of terms, (Pm, 0m, C0, am, and Km, need to be computed, for m = 1,..., M. One way to compute these terms would be to renumber the nodes. Therefore, if performance measures are required for, say, node 2, then we can consider renumbering nodes so that node 2 is now node 1, node 3 is now node 2, and so on. However, note that in order to compute (Pm, em, (om, m, cm, and Km, the only data parameters involved are Pm, 8m, 1m and Am. We use a temporary array variable, say Pm, to store the permuted Pm values for every other node whose performance measures are desired. We use an offset to index the nodes and assign the appropriate Pm values to the array Dpm, m=l,...,M. We use the same offset to recompute new 8m, 1m and Am values and proceed with Algorithm 1. The Special Cases of Pure Exhaustive and Pure Gated Systems If all the nodes are served using the same service discipline, then we can reduce the computational effort almost by half. For the pure exhaustive service system, the terms X(1)(j) and X(2)(j) are obtained from equations (5.6), (5.10a), and (5. 10b) as: X( m(j) = X m+l(j+l) X()(j) = Xm+ 2 + me. (5.15) m(2ll)( 1 Om+_..__..~l (j+l) m E 1E; (5.15) Pm+l M Pm+1 Pm+l Note that these terms are expressed solely as functions of Xm+1(j+1) and X (2)+1)* Corresponding expressions for m a, om and Km, are easily seen to be: 4m = (Pm+l/Pm2+l m = 0m+1/Pm3+l, me E, m M, (5.16a) M = (pI/P2 + 1, a = 01/p3 +1, me E, (5.16b) Km = (Wm+l - m+m+l )m+l)/pm2 mE AE. (5.16c) Therefore, if we let m = C m + Var(Sm-i)/Pm, (5.17) then the equations for f2) and f(3) for the pure exhaustive service system can be written as: 1 1 12

f(2) = f22 + X E (Pm 2m + 2 Var(SM), 1 1 h.dm 1 me E, (5.18a) m f(3) -X 3 flf(2)-2f3 + E0m[C7m+ T(S )-3m~m] +3C0mtm} +3 T(SM), me E(5.18b) m Pm Unlike the case of the pure exhaustive service system, for the pure gated service system X(1)(j) and X(m)(j) are expressed in terms of the parameters for node m only. The expressions for m, (nm and Km in this case are: m/p2 m = q)n, Km = m/Pm m = (Om -m Om)/Pm, m Ga, (5.19) and so, setting 'Tm = C m + Var(Sm)/pm, (5.20) the equations for f(2) and f(3) are: f(2) = f2 + 2 m Tm me a, (5.21a) m f(3) 3 fl f(2)_ 2 f3l +X31 { Om [ m + T(Sm) 31mm] + 3, m E a. (5.2 1 b) = + 3m'm'}, m i n. (5.21b) m Pm Thus, for either the pure exhaustive service system or the pure gated service system, we only evaluate (Pm, Om and om. This reduces the number of computations almost by half, compared to that needed for a system with a mix of the two service disciplines. 6. Computational Experience We conducted a number of experiments to determine the efficiency and speed of convergence of the algorithm, for various stopping tolerances,. (Recall that the algorithm terminates when the largest value of Xm(J), for a given j, is less than or equal to ~.) We found that a stopping tolerance of 0.0001 gave very accurate results, with the mean waiting times and variances matching the exact values up to about seven decimal places for the examples we report on below. (We did not check for numerical accuracy beyond seven decimal places.) We first present several examples, to give the reader a feel for actual numerical values. The program to run the experiments was coded in Pascal. The experiments were conducted on an IBM 9021-720 running the MTS operating system. 13

Example 1. The number of nodes = 5, with arrival rates 0.7, 0.8, 1.2, 0.1, and 0.2. The service times for the 5 nodes are uniformly distributed over the following ranges: [0, 0.8], [0, 1], [0, 0.2], [0.5, 1.5], and [0.2, 0.6]. The switchover times are exponentially distributed with means 0.5, 0.5, 0.4, 1.0, and 1.2. The utilizations at the five nodes are 0.28, 0.40, 0.12, 0.10, and 0.08 respectively, and p = 0.98. The mean waiting times and the waiting time variances are presented below for the following cases: i) the pure exhaustive service system, ii) the pure gated service system, iii) a system with a mix of service disciplines, where nodes 2, and 5 use the gated service discipline, while nodes 1, 3 and 4 use the exhaustive service discipline. Note that the uniform distribution with range [a,b] has 2nd and 3rd moments (b3 - a3)/3(b-a), and (b4 - a4)/4(b-a), respectively. The exponential distribution with parameter g has third moment 6/43 and third central moment 2/W3. i) The Pure Exhaustive Service System: The mean waiting times are E(W1) = 46.9850, E(W2) = 39.0331, E(W3) = 57.3817, E(W4) = 58.6123, and E(W5) = 59.9624. The waiting time variances are Var(W1) = 1,640.6459, Var(W2) = 1,118.7125, Var(W3) = 2,458.9258, Var(W4) = 2,556.9872, Var(W5) = 2,683.2556. The algorithm took 208 iterations to converge, and required 93.03 millseconds of CPU time. ii) The Pure Gated Service System: The mean waiting times are E(W1) = 72.2632, E(W2) = 79.0607, E(W3) = 63.2353, E(W4) = 62.1746, and E(W5) = 60.9668. The waiting time variances are Var(Wi) = 1,728.0505, Var(W2) = 1,741.9405, Var(W3) = 1,800.3384, Var(W4) = 1,829.5007, Var(W5) = 1,835.4479. The algorithm took 352 iterations to converge, and required 157.97 millseconds of CPU time. iii) The Mixed Service System: The mean waiting times are E(W1) = 42.2697, E(W2) = 82.1887, E(W3) = 51.7842, E(W4) = 52.8793, and E(W5) = 63.4502. The waiting time variances are Var(W1) = 1,142.2205, Var(W2) = 2,100.5381, Var(W3) = 1,731.4286, Var(W4)= 1,796.6408, Var(Ws) = 2,143.8311. The algorithm took 299 iterations to converge, and required 192.66 millseconds of CPU time. Example 2. The number of nodes = 5. The arrival rates are 0.20, 0.10, 0.16, 0.04, and 0.04. The service times at the nodes are uniformly distributed over the following ranges: [0, 6], [0, 4], [0, 1], [1, 2], and [0, 2]. The switchover times are exponentially distributed with means 1.0, 1.0, 0.8, 2.0, and 2.0. The utilizations at the five nodes are 0.60, 0.20, 0.08, 0.06, and 0.04 respectively and, as in the previous example, p = 0.98. The mean waiting times and the waiting time variances are presented below for: i) the pure exhaustive service system, ii) the pure gated 14

service system, iii) a system with a mix of service disciplines, where nodes 2, and 5 use the gated service discipline, while nodes 1, 3 and 4 use the exhaustive service discipline. i) The Pure Exhaustive Service System: The mean waiting times are E(W1) = 122.1462, E(W2) = 245.0320, E(W3) = 282.8141, E(W4) = 289.7239, and E(W5) = 296.4669. The waiting time variances are Var(Wi) = 13,245.9090, Var(W2) = 56,148.0645, Var(W3) = 75,071.5516, Var(W4) = 78,832.1916, Var(Ws) = 82,618.0849. The algorithm took 176 iterations to converge, and required 78.49 millseconds of CPU time. ii) The Pure Gated Service System: The mean waiting times are E(W1) = 361.2300, E(W2) = 270.6877, E(W3) = 243.5746, E(W4) = 239.2250, and E(W5) = 234.7594. The waiting time variances are Var(W1) = 36,125.1632, Var(W2) = 31,690.2523, Var(W3) = 32,496.2060, Var(W4) = 32,824.0134, Var(Ws) = 33,115.3918. The algorithm took 403 iterations to converge, and required 179.99 millseconds of CPU time. iii) The Mixed Service System: The mean waiting times are E(W1) = 115.3107, E(W2)= 345.9230, E(W3) = 267.3364, E(W4) = 273.8346, and E(W5) = 302.9340. The waiting time variances are Var(Wi) = 11,145.7619, Var(W2) = 75,274.2534, Var(W3) = 63,820.7406, Var(W4) = 67,018.8124, Var(Ws) = 75,651.9396. The algorithm took 201 iterations to converge, and required 121.91 millseconds of CPU time. Example 3. The number of nodes = 10, with arrival rate equal to 0.8 at node 1, and equal to 0.1 at all the other nodes. The service times are exponentially distributed with means 0.5, 0.5, 0.4, 1.0, 0.4, 0.5, 0.5, 0.4, 1.0, and 0.4. The switchover times are uniformly distributed over the following ranges: [0, 1], [0, 1], [0, 0.8], [0.5, 1.5], [0, 0.8], [0, 1], [0, 1], [0, 0.8], [0.5, 1.5], and [0, 0.8]. The utilizations at the ten nodes are 0.40, 0.05, 0.04, 0.10, 0.04, 0.05, 0.05, 0.04., 0.10, and 0.04, respectively, and p = 0.91. The mean waiting times and the waiting time variances are presented below for: i) the pure exhaustive service system, ii) the pure gated service system, iii) a system with a mix of service disciplines, where nodes 1 through 5 use the gated service discipline, while nodes 6 through 10 use the exhaustive service discipline. i) T'he Pure Exhaustive Service System: The mean waiting times are E(Wi) = 23.2652, E(W2) = 36.8740, E(W3) = 37.3111, E(W4) = 34.8626, and E(W5) = 37.2287, E(W6) = 36.8793, E(W7) = 36.9212, E(W8) = 37.3629, E(W9) = 34.9234, and E(W1o) = 37.2999. The waiting time variances are Var(W1) = 313.1774, Var(W2) = 825.2952, Var(W3) = 847.6837, Var(W4) = 729.0181, Var(W5) = 840.5882, Var(W6) = 825.7207, Var(W7) = 829.3474, Var(W8) = 15

852.1895, Var(Wg) = 734.0381, Var(W1o) = 846.8985. The algorithm took 50 iterations to converge, and required 90.73 millseconds of CPU time. ii) The Pure Gated Service System: The mean waiting times are E(W1) = 50.4670, E(W2) = 37.9978, E(W3) = 37.6446, E(W4) = 39.8357, and E(W5) = 37.5621, E(W6) = 37.9356, E(W7) = 37.9392, E(W8) = 37.5820, E(W9) = 39.7653, and E(W10) = 37.4834. The waiting time variances are Var(W1) = 548.3302, Var(W2) = 658.8032, Var(W3) = 664.4823, Var(W4) = 640.3986, Var(Ws) = 658.1578, Var(W6) = 654.1311, Var(W7) = 654.4568, Var(W8) = 659.7992, Var(Wg) = 635.3322, Var(Wlo) = 652.2293. The algorithm took 75 iterations to converge, and required 135.62 millseconds of CPU time. iii) The Mixed Service System: The mean waiting times are E(W1) = 50.6987, E(W2) = 38.1981, E(W3) = 37.8473, E(W4) = 40.0546, and E(W5) = 37.7763, E(W6) = 34.5286, E(W7) = 34.5329, E(W8) = 34.9119, E(W9) = 32.5897, and E(Wo1) = 34.7659. The waiting time variances are Var(W) = 565.0151, Var(W2) = 674.7923, Var(W3) = 680.7920, Var(W4) = 657.3561, Var(W5) = 675.3881, Var(W6) = 629.7939, Var(W7) = 630.2239, Var(Wg) = 645.5318, Var(Wg) = 550.7752, Var(W10) = 634.0304. The algorithm took 73 iterations to converge, and required 187.40 millseconds of CPU time. The IS technique we have developed is based on non-zero switchover times. However, the technique can be used to obtain the mean and variance of the waiting time at a node for systems in which the switchover times are arbitrarily close to zero as we show in example 4 below. The data in example 4 is the same as in example 3, except that the mean switchover times are very close to zero, and the variances are set equal to zero. Computational Remark:1 Note that the terms (pm, Om m, m, am, and Km, depend only on n the arrival rates and the service time parameters. These terms are independent of the switchover time parameters. Hence, once these terms are computed for a given set of arrival rate and service time parameters, we can obtain the mean and variance of waiting time distributions at the various nodes, for different switchover time parameters, through just a few arithmetic operations. U Example 4. The data parameters are exactly as in example 3, except that the mean switchover time at each node are all equal to 0.000001, and the variance of the switchover time is equal to 0. The mean waiting times and the waiting time variances are presented below for: i) the pure exhaustive service system, ii) the pure gated service system, iii) a system with a mix of service 1 The author thanks Professor Robert Cooper for a very fruitful discussion on this feature of the IS Technique. 16

disciplines, where nodes 1 through 5 use the gated service discipline, while nodes 6 through 10 use the exhaustive service discipline. i) The Pure Exhaustive Service System: The mean waiting times are E(W1) = 4.5488, E(W2) = 7.2406, E(W3) = 7.3670, E(W4) = 6.7888, and E(W5) = 7.2846, E(W6) = 7.2473, E(W7) = 7.2904, E(W8) = 7.4216, E(W9) = 6.8530, and E(W1o) = 7.3598. The waiting time variances are Var(Wi) = 37.6933, Var(W2) = 102.7229, Var(W3) = 106.2822, Var(W4) = 89.3060, Var(W5) = 104.1149, Var(W6) = 102.8147, Var(W7) = 103.9472, Var(W) = 107.6527, Var(W9) = 90.8709, Var(Wlo) = 106.1094. ii) The Pure Gated Service System: The mean waiting times are E(W1) = 6.8397, E(W2) = 5.2727, E(W3) = 5.2319, E(W4) = 5.5535, and E(W5) = 5.1496, E(W6) = 5.2119, E(W7) = 5.2164, E(Wg) = 5.1717, E(W9) = 5.4857, and E(W1o) = 5.0740. The waiting time variances are Var(Wi) = 63.3489, Var(W2) = 49.2457, Var(W3) = 49.0801, Var(W4) = 52.0203, Var(Ws) = 47.7020, Var(W6) = 48.2603, Var(W7) = 48.3642, Var(Ws) = 48.1301, Var(W9) = 50.9621, Var(Wlo) = 46.4690. iii) The Mixed Service System: The mean waiting times are E(W1) = 7.0694, E(W2) = 5.4712, E(W3) = 5.4329, E(W4) = 5.7706, and E(W5) = 5.3621, E(W6) = 4.9190, E(W7) = 4.9242, E(Wg) = 4.9927, E(W9) = 4.5401, and E(W10) = 4.8476. The waiting time variances are Var(Wi) = 68.0742, Var(W2) = 53.3218, Var(W3) = 53.2248, Var(W4) = 56.4554, Var(Ws) = 52.0771, Var(W6) = 46.0762, Var(W7) = 46.2455, Var(W8) = 47.6637, Var(W9) = 38.5235, Var(W1o) = 44.8867. A few remarks are made based on the experiments. It may be observed that in all the examrples, the mean waiting times in the pure gated service system are higher at the nodes with higher utilizations. For the pure exhaustive service system, however, the reverse holds: the mean waiting times are higher at nodes with lower utilizations. (For the exhaustive service system, the same pattern of behavior appears to hold even for the waiting time variances.) It is observed that the waiting time variances at the nodes for the pure gated service system are somewhat close to one another; we observed this behavior in all the experiments conducted. (In fact, the closeness is even more pronounced if we compare the standard deviations of the waiting time distributions, instead of the variances.) It may also be observed that for the mixed service systems, the waiting time variances at the nodes served using the gated service discipline are fairly close to one another. This suggests that while the gated service system results in higher mean waiting times compared to the exhaustive service system, the waiting time variances are "smoothed 17

out." This is a positive factor which has to be given due consideration in the design and operation of a polling system. The relative closeness of the waiting time variances certainly does not appear to hold for the pure exhaustive service system. On the other hand, interestingly enough, we observe that the standard deviation of the waiting time distribution tracks the mean of the distribution somewhat closely, in the case of the pure exhaustive service system. Also, as remarked upon earlier, for pure exhaustive service systems the standard deviation (at a node) is high relative to that at other nodes, when the corresponding mean waiting time is high. Computational Efficiency In general, the computational results suggest that the algorithm requires roughly the same number of iterations for a given value of p, for any number of nodes. A primary advantage of the IS technique is that the user can obtain the waiting time variance at select nodes without having to evaluate waiting time variances at all nodes simultaneously. It would be of interest to compare the relative efficiency of the IS technique with the algorithm of Ferguson [5]. We note that the algorithm of Ferguson involves the solution of O(M3) equations (and, thereby, requires O(M9) operations). Since the IS technique is an iterative technique, it is therefore difficult to make any immediate conclusions about the relative efficiency of either algorithm. Also, since there are no computational results reported in the paper by Ferguson, we are unable to make any comparisons on execution times. However, we can indirectly estimate the efficiency of the IS technique as follows. First, we observe, based on the expressions for f(2) and f(3), that the IS technique for computing waiting time variances requires about three times the computational effort required to just compute mean waiting times using the same approach. In an earlier paper (Srinivasan [8]), we had published extensive results on the computational effort required to obtain the mean waiting times. The numerical comparisons indicated that for p values of about 0.90 or less, the iterative algorithm to compute mean waiting times compared very favorably with an O(M3) algorithm, even when the iterative algorithm was used to obtain mean waiting times at all the nodes. It is, therefore, not unreasonable to conclude that if the IS technique is used to compute waiting time variances at all the nodes, it would still compare very favorably with any algorithm involving O(M3) computations, at least for p values of about 0.90 or less. Note, in addition, that the IS technique has a major advantage in that one can compute waiting time variances at select nodes, without having to compute all the waiting time variances at all nodes simultaneously. 18

The special case of the system with two nodes n this case, we can obtain the waiting time variances at the nodes in closed form, for the case of the pure exhaustive service system. It is also possible to obtain the waiting time variances for either the pure gated service system or the mixed service system with two nodes, without having to perform any iterations at all. However, the expression for the waiting time variances in these latter cases is a little more complicated, and involves the solution of a system of upto three equations; we do not present these cases here. Define I = 1 82, and 8mim Pm Pm (6.1) For the pure exhaustive service siystem with two nodes, the terms l() and X2() are obtained very simply as: rJ xi(j) = rJ, x2() = -, The terms X(2)(j) and )(j) are obtained, for j > 1, as: 1v uu~ —/CL VLLLL~IIJ/1 3 j 1l. (6.2) X2)(j) = r X(2)(-1+ )2 (+ I ~j X(2)) = FX2)-1)+ (-+ (6.3) x = F(2)81 IF 51 The terms (pm, Om, and on, m = 1, 2, are thus evaluated in closed form: r2 -= lr2 r2 (92 = -(1-F2)?(1-r2) r3 01 = -, 1-r3 r3 02 = (-r3) 5 (1-3) (6.4a) (6.4b) (6.4c) 1 + C I2/61 )1l = - 01, ' l-r2 181 + C2 r_ (022 - 92 C1' r(l-r2) Substituting these values in equations (5.18a) and (5.18b), we obtain f2) and f(3). The first and second moments of the waiting time at node 1 are then obtained from equations (3.4a) and (3.5a). 19

7. Summary and Conclusions We have obtained the waiting time variances in a continuous-time polling system using a powerful technique, which we term the Individual Station (IS) technique. The polling system can adopt either the pure exhaustive service discipline, the pure gated service discipline, or a system with a mix of both disciplines. The IS technique allows the user to obtain the waiting time variances at one or more select nodes in the polling system. Furthermore, the IS technique has the following additional advantage. The resulting expression for the second moment (or the mean) of the waiting time at a node is composed of two parts, one of which is independent of the switchover time parameters. Therefore, once the means and second moments of the waiting times have been obtained for a given set of service time and switchover time parameters, the user can compute means and second moments of the waiting times for many alternate configurations involving different switchover time parameters, with just a few elementary arithmetic operations. The IS technique for computing waiting time variances requires about three times the computational effort required to obtain the mean waiting times using the same iterative algorithm. In a previous paper we had shown, empirically, that the IS technique using the iterative algorithm computes the mean waiting times at all the nodes in a time comparable to an O(M3) algorithm, as long as p is about 0.90 or less. Therefore, this suggests that the performance of the iterative algorithm for computing waiting time variances at all nodes would also be comparable to an O(M3) algorithm, as long as p is about 0.90 or less. There are a number of topics to be explored, on further applications of the IS technique. Although we can obtain the performance measures for systems in which the switchover times are arbitrarily close to zero, the technique is based on non-zero switchover times. It would be of interest to develop the IS technique for systems in which the swithcover times are all equal to zero. It would also be of interest to consider using the IS technique to analyze other polling systems. 20

References [1] Y. Aminetzah. 1975. An exact approach to the polling system. Ph.D. dissertation, McGill University, Montreal, P.Q., Canada. [2] R.B. Cooper and G. Murray. 1969. Queues served in cyclic order. Bell System Technical Journal, v 48, pp. 675-689. [3] R.B. Cooper. 1970. Queues served in cyclic order: waiting times. Bell System Technical Journal, v 49, pp. 399-413. [4] M. Eisenberg. 1972. Queues with periodic service and changeover time. Operations Research, v 20, pp. 440-451. [5] M.J. Ferguson. 1986. Computation of the variance of the waiting time for token rings. IEEE Journal on Selected Areas in Communications, v SAC-4, pp. 775-782. [6] L. Kleinrock. 1975. Queueing Systems - Volume I: Theory, John Wiley, New York. [7] H. Levy. 1991. A note on the complexity of Swartz's method for calculating the expected delay in non-symmetric cyclic polling systems. Operations Research Letters, v 10, pp. 363-368. [8] M.M. Srinivasan 1991. The Individual Station Technique for the Analysis of ComtinuousTime Polling Systems. Technical Report 91-36, Department of Industrial and Operations Engineering, The University of Michigan, Ann Arbor, MI 48109-2117. [9] H. Takagi. 1986. Analysis of polling systems, MIT Press, Cambridge, MA. [10] H. Takagi. 1990. Queueing analysis of polling models: an update. Stochastic Analysis of Computer and Communication Systems, Elsevier Science Publishers B.V. (NorthHolland), Amsterdam. 21