THE INDIVIDUAL STATION TECHNIQUE FOR THE ANALYSIS OF CONTINUOUS-TIME POLLING SYSTEMS Mandyam M. Srinivasan Department of Industrial & Operations Engineering University of Michigan Ann Arbor, MI 48109-2117 Technical Report 91-36 December 1991

The Individual Station Technique for the Analysis of 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-36 December 1991 Abstract Polling systems find application in a wide variety of real-world systems, for example, in telecommunication systems and material handling systems, and there is continued interest in the analysis of polling systems for their performance. In this paper, we present a very promising technique to obtain mean waiting times for continuous-time polling systems which use the exhaustive and gated service disciplines. The technique enables the analysis of polling systems which have a mix of these service disciplines. Thus, some nodes could adopt the exhaustive service discipline while the other nodes use the gated service discipline. The technique also allows the user to obtain the mean waiting time at select nodes without having to obtain all the mean waiting times simultaneously.

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 telecommunication and material handling systems, and there is continued interest in developing new algorithms and applications for polling systems. An excellent survey of the research done on polling systems can be found in Takagi (1990). A number of service disciplines have been studied, which determine the service strategy adopted by the server at a node. 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, and iii) the k-limited service discipline in which the server serves at most k customers during a visit to that node. Several variants of these strategies have been proposed and studied in the past. In this paper we present a new and powerful technique for the analysis of continuous-time polling systems, which we term the Individual Station (IS) technique. The IS technique allows us to determine the mean waiting time at one or more select nodes without having to obtain mean waiting times at all the nodes, simultaneously. We develop the following new results: a) We provide an iterative algorithm to determine the exact mean waiting times for the asymmetric, continuous-time, exhaustive and gated service polling systems. b) We next present an algorithm which computes mean waiting times very rapidly without recourse to iteration. The mean waiting time at a node is determined, in an explicitform (through the solution to a system of six equations or less), and requires O(M2) operations per node. The algorithm is numerically exact for polling systems with upto seven nodes. We demonstrate that for polling systems with an arbitrarily larger number of nodes, the solution of this system of equations results in mean waiting time estimates very close to the exact values, with errors less than 0.0025%, even when utilizations approach 1. For the exhaustive service polling system with more than seven nodes, we also show that the mean waiting time obtained from the algorithm bounds the exact numerical value from below. If mean waiting times need to be computed even more accurately, they can be obtained with very little additional effort, simply by solving a slightly larger system of equations (say, a system of seven equations). c) We obtain the mean waiting times for a polling system which allows mixed service strategies, wherein some of the nodes are served using the exhaustive service discipline, while the other nodes are served using the gated service discipline. 1

1. Literature Review There is a considerable amount of literature devoted to polling systems analysis, and the reader is referred to the paper by Takagi (1990) for a comprehensive survey of the work done in this area. We will briefly review some of the papers most pertinent to our work. A number of numerical techniques have been proposed for computing the mean waiting times for the exhaustive and gated polling systems. Levy (1991) classifies these techniques as a) The buffer-occupancy equation method (Cooper and Murray 1969, Cooper and Murray 1970, Eisenberg 1972, Hashida 1972, Konheim and Meister 1974, Rubin and Demoraesl983, Kleinrock and Levy 1988, and Takagi 1986), b) The station-time equation method (Aminetzah 1975, Humblet 1978, Ferguson and Aminetzah 1985, and Baker and Rubin 1987), c) Carsten, Newhall and Posner's method (1977), d) Swartz's method (1980), and e) Sarkar and Zangwill's method (1989). Most of the above methods require the solution of large sets of linear equations (either M2 or M3 equations). The algorithm of Sarkar and Zangwill, however, obtains the mean waiting times by solving a system of M linear equations, thereby requiring only O(M3) operations. Swartz (1980) develops an iterative procedure for analyzing the discrete-time exhaustive service system, in which time is slotted in fixed intervals. Although this procedure does not guarantee an upper bound on the number of iterations required, it converges rather quickly when the utilization of the server is not very close to 1. Levy (1989, 1991) discusses the complexity of this algorithm, and the rate at which it converges. Levy (1991) observes that Swartz's method can be used to calculate the mean waiting time at a station independently of the mean waiting times at other stations. This is in contrast to the other methods which solve for all M mean waiting times, simultaneously. Levy shows that the method forms a contraction mapping and, therefore, the number of iterations it requires is logarithmic in the accuracy required. Levy concludes that "for a wide range of parameters, the method is the most efficient method known today for computing the expected delay in polling systems," and that it is desirable to apply this approach to polling systems other than the discrete-time exhaustive service system. The iterative algorithm we present in this paper (for the continuous-time exhaustive and gated service polling systems) is similar to Swartz's approach for the discrete-time exhaustive service system. It can, therefore, be shown to share the above property of logarithmic convergence. 2

We first consider the exhaustive service system, and derive the expression for the mean waiting time in Section 2. In Section 3, we present the iterative algorithm, using the IS technique. We next develop the algorithm which requires O(M2) operations to obtain the mean waiting time at a node. Section 4 presents a number of computational results for the two algorithms. These results serve to illustrate both the speed of convergence of the iterative algorithm, and demonstrate the accuracy of the second algorithm. In section 5, we extend the results of sections 3 and 4 to the gated service system. In section 6, we show how the IS technique easily extends to handle polling systems in which some of the nodes adopt the exhaustive service discipline while the others adopt the gated service discipline. Section 7 presents a summary and discusses ways in which the IS technique could be further extended. 2. The Model for the Exhaustive Service Polling System We present most of the notation that will be used in the paper. 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) res, et 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, according to the exhaustive service discipline. Customers arrive at each node according to independent Poisson processes with rate?m at node m. The service time of a customer at node m is an independent random variable, having finite mean bm and second moment b(). The Laplace Stieltjes Transform (LST) of the service time is denoted by 3m(.), and the LST of the busy period (Cooper 1981) induced by a single customer at node m is denoted by rlm(.). The time taken by the server to switch between nodes m and m+l is a random variable with mean sm, second moment s(2), variance Var(Sm), and LST am(.). The traffic intensity at node m is denoted by Pm = Xmbm, and p = Em Pm denotes the server utilization. For the polling system to be stable, we require p to be less than one. The sum of the mean switchover times is denoted by s = Xm Sm. We are interested in obtaining the mean waiting time, W,, at node m, for m=1,...,M. Let z = (z1,...,ZM), and let Fm(z) denote the probability generating function for the number of customers present at each node at the instant the server polls node m. The first and second factorial moments of the number of customers present at node m when the server polls the node are given by fm and f(m, respectively, where aF ( 1,z) fD2) 2Fm(z) fm = I z=1, and f() iZm2 Iz=1. (2.1) Mzm aZ,2 3

If we can obtain the first and second factorial moments of the number of customers present at node m when the server polls that node, then the mean waiting time is obtained as (Takagi 1986): ft2) Xmb(2) Wm = m + 1bm (2.2) 2XnfM 2(1-pm) Obtaining these moments is not straightforward, since the generating function Fm(.) is expressed, recursively, in terms of the generating function Fm-_(.). It is well known (see for instance, Takagi 1986) that the generating functions are related by the following expression: Fm+l(Zi,...,ZM) = Fm(Zl,..,m-l_,lm(Ykem (Xk-XkZk)),Zm+l,..,ZM) om(Xk (Xk-XkZk))- (2.3) This equation forms the starting point for the IS technique. We initially obtain an expression which expresses Fm(.) solely in terms of functions of input parameters. We then differentiate the resulting expression twice and use equation (2.2) to obtain the mean waiting times. The technique we develop below, computes W1. The mean waiting times at the other nodes are obtained in an identical manner. We first define a recursive (nested) function, ym(j) as follows: Ym(O) = Zm, m=1,...,M, (2.4a) Ym(i) = Tm( X [k-xkYk(i-1)] + [k-XkYk(j)]), j > 0, m = 1,.....M. (2.4b) k<m k>m Let om(j) = (m( [Xk-kYk(j)] + [Xk-kYk(j+l)]), j > 0, m = 1,...M. (2.5) k<m k>m We can now cast equation (2.3) in the following form: F,(y1(0),...,YM(O)) = FM(y1(O),...,'M-1(0),YM(1)) OM(O). Recursively expressing FM(.) in terms of FM-1(.), and FM-I(.) in terms of FM-2(.), and so on, we obtain: F1(yl(0),...,yM(0)) = Fl(yl(l),...,YM(l)) 1 am(0). m We now continue to recursively express F(yi(j),...,M(j)) in terms of F(yI(j+1),...,yM(j+l)), for j > 1, and this results in the following expression for F(y71(0),...,yM(0)): n-l Fi(y1(0),...,yM(0)) = F,(yl(n),...,yM(n)) (T m(j)) (2.6) j=0 m ) It can be shown (refer Cooper 1969, and Eisenberg 1972, for example), that limn->ooyl(n) = 1. Hence, letting n -- oo in equation (2.6), we obtain: 4

00 Fl(y(O),...,M(O)) = n m() (2.7) j=o {m ) We now have an equation for F1 that is only in terms of (complex) functions of the am(j) values. Differentiating equation (2.7) once and twice with respect to z1, and setting z = 1, we obtain: f I= m 'z=1, and (2.8a) j=O m 00 D2(ymo) oo {D}ym9 0 f(2) (= + i) 1 - ( j) I=1. (2.8b) j=O m j=O m Let Vm() = Xm a z=1, Ix1(j) = Xm a ) z=l, and (2.9) 6m = m (2.10) 1-Pm We obtain Vmo() and V(m)(j, for j > 0, by differentiating equation (2.4b) once and twice with respect to zl, and then setting z equal to 1, to get: VNm() = m[ Wk(j-l) + Wk(J)], j>0, and (2.11a) k<m k>m -2(j) = Sm[ X V(2)0-1) + A ()(j)] + -mb [ W Vk(j-1) + Z Vk(j)]2, j > 0. (2.11lb) k<m k>m ( 3 k<m k>m In a similar manner, from equation (2.5) we obtain, for j > 0: ar Iz=l = Sm[ Nk() + Nfk(j+l) ] and (2.12a) k<m k>m azmy2 z= = Sm [ g)(j) + V(2)(j+1)] + S(2)[ k(j) + Vk(j+l) ]2. (2.12b) k<m k>m k<m k>m Hence, from equations (2.8a) and (2.12a), interchanging the order of summations, we obtain 00 fi = Sm [ vk(j) + Vk(J+l) ] (2.13) m j=0 k<m k>m Observation 1: a) The term Wk(0) = k-Z I!z=l = Xi when k=l; it is equal to zero otherwise. b) The term V2(2)(0) = 0 for all k. U Based on Observation la), we can rewrite equation (2.13) as 5

00 00 00 fi = X Sm, X~ j, V X k(j) s= S kk(j) (2.14) m j=O k j=O k k j=O Let Xk =,j=O k(J). (2.15) To evaluate Xk, we use Observation la) once more. From equation (2.1 la), we obtain: X i + 81 ~ kel Xk, m = 1, Xm lm Ckzm Xk, otherwise. We add 6m Xm to both sides of the above equation, and divide both sides of the resulting identity by (1+5m). Next, we use the following facts: a) 6m/(l+6m) = Pm, and b) l/(l+8m) = (l-pm). This results in the following expression for xm: { (l-p)iAL + pi k Xk m = 1, Pm =k Xk otherwise. We now sum xm over all m, to get: k Xk = i 1-p (2.16) 1-p From equations (2.15) and (2.16), we thus obtain: 1-pl f, = 1P s (2.17) 1-p We can now obtain f(2) using an approach similar to that used to obtain fl. From equations (2.8b), (2.12a), and (2.12b), we obtain the following expression for f(2): 00 f(2) = f +, (Sm [ V(2)(J) + (2k(+1)] + Var(Sm)[ Vtk(O) + Nfk(j+l) ]2) j=O m k<m k>m k<m k>m Let Yk = jowl V)(O) Interchanging the order of summation wherever possible in the above equation, and using Observation lb) to note that V(2)(O) = 0 for all k, we get: 00 f2) = 2 + Sm Z Yk + Z Var(Sm) i [ xvk(j) + ~ k(j+l) ]2. (2.18) m k m j=O k<m k>m Similar to the approach used to determine ~m xm, we obtain Xm ym, after some algebra, as Xmbm() 00 X Ym -= X'1'iS^Qm)ul2 (2.19) m m (-P)Pm j=1 6

From equation (2.11a), for m=l,...,M-l, [ N k(J) k<m + vk(+l)]2 k>m [ X k(J) + Vm+l(i+l) k<m + VkO+l) ]2 k>m+l = [W1 ) +,m+l(j+l)]2 Bm+l 2m+l Pm+i (2.20) For m=M, from equation (2.1 la): [C vk()]2 = k M X2, j=O, 1' and [ k(i)]2 k<M i2(j) 2' P1 j>0. (2.21) Hence, from equations (2.18) through (2.21), we obtain: f(2) = f2 + s m Xnb(2) ~o ( -P)P~ 2 t2(J) + Var(Sm) m Pm oo j=l + X2 Var(SM). I ""M/ (2.22) We can rewrite equation (2.22) as f) = f + 1 1 Xmb(1-) 0 m (1-p)pm j=l Vm(i) + m Var(Sm,,,-) 00 2 Pm j-~ + X2 Var(SM). (2.23) It is now convenient to define two terms, Xm(J) and (pm(n), as follows. Xm(J) m aYm(O) azl ' and (2.24a) X1 (pm(n) 00 Z Xm(j)Xm(j-n)j=n+l (2.24b) From equations (2.23) and (2.24), we get an expression for f2) given below as Lemma 2.1. Lemma 2.1: f(2) 1 2mb(2) _m m 1-p Pm2 m Var(Sml) ) + 2 Var(SM). Pm2 1 Frcxm Lemma 2.1 and equation (2.17), we obtain the main result of this section, stated as: Theorem 2.2: The mean waiting time at node 1 for the exhaustive service system is: WI = + 2 1-p m (Pm(O) _ +( Vb(Sm() 1-p Var(SM)).(2.25) m -Var(Sm))+( (l-p,)p2m 2 +2s ) (l-pO 2s(l-p,) 7

The expressions for the mean waiting times at the other nodes can be obtained in a similar manner, simply by renumbering the indices. Thus, obtaining the mean waiting times reduces to computing the term (Pm(O), since the other terms are just input data. In the following section, we develop the two algorithms to compute pm(0). 3. Computing (pm(O) for the exhaustive service system The two algorithms we present in this section differ in computing (pm(O) as follows. The first algorithm iteratively computes X%(j), starting with j = 1 until the terms approach 0. In the second algorithm, we only compute x(ji) for j=l,...,M-l, and then solve a system of equations as explained subsequently. We first present a simple expression for recursively computing Xm(j). From equations (2.1 la) and (2.24a), after some elementary algebra, we obtain: %M(1) = sM; Xm(1) = (Xm+l(1)/Pm+l)6m, m=l,...,M-1, (3.1a) and for j > 1, XM(j) = (%l(j-l)/p - XM(j-1))6M; Xm(J) = (Xm+l(j)/Pm+l-Xm(j-l))6m, m=l,...,M-l. (3.1b) Algorithm 1: 0. Setj = 1, andcompute m(j), m=l,...,M, usingequation (3.1a). Initialize (m(O) = X(j) form = 1,...,M. Choose a tolerance value e. 1. Setj =j+l. 2. Compute Xmi), using equation (3.1b), and set (pm() = (Pm(O) + x2(j), m 1,..,M. If Xm(J) < e for all m, stop. Otherwise go to step 1. 3. Obtain the mean waiting time at node 1, using equation (2.25). Algorithm 1 is executed for each node whose mean waiting time is desired. Since computation of (pm(O) only needs the terms Pm and 6m = Pm/(l-pm), we can compute the mean waiting time for, say, node 2 simply by using a temporary variable, P,, to store permuted Pm values. To compute W2, we just set Pm = Pm+l, m=l,...,M, and (re)compute new 5m values using the pm values. For Algorithm 2, we need to define some additional terms, Fm, m=l,...,M. Let M FM = 1 8m, and Fm = m+ + K),1 < m<M, where (3.2) m k=m m-1 K(m) = 8 and c() = 6m K(i'1), m = 1,...,M. (3.3) j=i-1 8

Intuitively, K(i) denotes the sum of all possible combinations of i unique 8k terms, such that 5m is in every one of these combinations. The term Fm denotes the sum of all possible combinations of m unique 8k terms, plus the sum of all possible combinations of m+1 unique 6k terms +... + the sum of all possible combinations of M unique 8k terms. (There is, of course, only one combination of M unique 6k terms.) Note that all the Fm terms, m=l,...,M, are strictly positive. Lemma 3.1 shows that the sum of Fm, m = 2,...,M is less than one. The proof of Lemma 3.1 easily follows by induction on M, and it is omitted here. Lemma 3.1: M (1- n) II(1-Pm) = 1-p. U n=2 m The next lemma, Lemma 3.2, expresses Xm(J) in terms of Xm(j-n), for n>O. The proof of Lemma 3.2 follows by induction on the number of nodes and is omitted. Although Lemma 3.2 requires the definition of an additional variable, 6(, we remark that this variable is needed only to establish some initial conditions. It is not used subsequently. The variables 0( are defined recursively as follows: (m) = m, and for i > 1, ( = K-iLemma 3.2: M-1 Xm(j) = SXm(j-n) n+, j > M, m = 1,...,M, with (3.4) n=l ^ ~~M Xm(j) = Xm(j-n) Fn+ + e0) l (l+6n), j < M, m = 2,...,M, (3.5a) n=l n=m+l i ~M M i(J ) = ZXi(-)Fn+l - +) (l+sn), j < M. (3.5b) n=l m=j+2 n=m+l We now show how Lemma 3.2 enables us to obtain the mean waiting time at a node. To clarify the technique that we subsequently develop for the general case, we first consider the special case for M=3 nodes. It will be convenient to define: 00 m(n) = Xm(j)Xm(j-n). (3.6) j=M Note that 9

M-1 (Pm(n) = (pm(n) + Z Xm(j)Xm(j-n). (3.7) j=l We remind the reader that we are interested in computing (pm(O). 3.1. Algorithm 2: The special case of M=3 nodes For this case, we have from Lemma 3.2, 00 00 M0) = X2 (j) = [2 Xm(i-1) + r3 Xm(j-2)]2 j=3 j=3 00 = i [-r2 x2m(-1) + 2i X2(j-2) + 2r2r3m —1)Xm(ji-2)] j=3 r22 2(2) + F2 pm(O) + r2 (m2 (1)+2m(2)) + rF2 Pm(0) + 2r2r3Xm(l)Xm(2) + 2r2r3(m(1). The term m(1) is, in turn, obtained from another application of Lemma 3.2 as: 00 00 m(l) = X Xm(i)Xm(j-1) = E [F2 Xm(-1) + F3 Xm(j-2)] Xm(-1) j=3 j=3 = F2 2(2) + r2 (m(O) + r3Xm(1)Xm(2) + F3 (pm(l) This results in two equations in the two unknowns, (Pm(O) and (pm(l). Solving for qm(O), we get: X2m(2) (r2 + r32 + 2r2 r3/(1 -r3)) + X 2() r + 2m(1) Xm(2) r2r3/(1-r3) mW (O) = X 2 + 23 ) 3 XmlXm22. (3.8a) -"2- r]2- 2r2 F3/(1-F3) To compute (pm(O), from equation (3.7) we add x2 (1)+X2(2) to (Pm(O), and simplify the result to get X2(2)+ %m(i1) + 2Xm (1)Zm(2)r2r3/(1-r3) pm(O) = 2 (1) + 2 m(l)rXm(2)F2/iF3/ (3.8b) 1 - rF2 - _ F 22 F2 3/( 1-3) Note that for three nodes, F2 and F3 are computed very simply as F3 = 618263, and F2 = F3 + 6182 + 8283 + 8631. Note, too, that we only compute F2 and r3 once, since these terms are permutation-invariant. However, for each node for which mean waiting times are desired, we must compute the Xm(l) and Xm(2) terms, m=l,...M, after permuting the node indices as mentioned earlier. We can directly substitute the expression for (p(O), obtained from equation (3.8b), in equation (2.25) to obtain the mean waiting time. 10

3.2. Algorithm 2: The general case For an arbitrary number of nodes, using a similar approach as in section 3.1, we can show that (P,,(O) is obtained by solving the following system of equations: M M-2 M-i q5m(0) = bm(O) + Sm(0) r 2 + 2 A fm(i) Fn rn+i, (3.9a) n=2 i=l n=2 M-k-1 (pm(k) = bm(k) + E qm(n)rn+k+i, M-2 2 k > 0. (3.9b) n=l-k In equation (3.9b), we adopt the convention that ~m(-n) = PSm(n), n 2 0. In general, to obtain the Om(O) values for a system with M nodes, we first need to evaluate rn from n=2 upto n=M. The terms bm(k) in the above equations are data and are obtained as follows: M-1 k bm(0) = 2 g2(n), bm(k) = Xm(M-l-k+n) gm(n), k = 1,...,M-2, (3.10) n=1 n=l where gm(n) is computed, recursively, from M-1 gm(n) = E Xm(M-l-k+n) Fk+l. (3.11) k=n We remind the reader that in order to compute the mean waiting time at a node, the bm(n) terms have to be computed for each m. Computing the bm(n) terms from equations (3.10) and (3.11) would require O(M2) effort for each m, and hence it would require O(M3) operations to determine all the coefficients. However, as we now show, it is not necessary to evaluate all the bm(n) terms. This is based on the fact that the Fn terms converge to 0 very rapidly. For example, suppose we make the reasonable assumption that the value of rn is maximized, for any n>2, when the 6m values are all the same. Under this assumption, the largest possible for 8m would result when p approaches 1 and, since all Pm values are equal to p/M, 5m would approach 1/(M-l). Since Fn is the sum of all possible combinations of n unique 8k terms plus the sum of all possible combinations of n+1 unique 8k terms +... + the sum of all possible combinations of M unique 8k terms, this, in turn, implies that the largest possible value for Fn would be AM ( ")(MlJ. We know that the sum i=nM i =1Mi im_ + Mo ( i M = 1 + MlJ^, and that limM-(l + Ml )M = e. This, in conjunction with Lemma 3.1, leads us to the following observation: 11

Observation 2. The terms Fn decrease at an exponential rate with increasing n. For n > 2, the magnitude of Fn is on the order of n. Based on the above observation, instead of solving the system of equations specified by (3.9), we consider solving a reduced system of equations, ignoring terms involving rk, for all k greater than some n. We ran a number of experiments, with values of p ranging from 0.3 to 0.996, and with the number of nodes ranging from 3 to 50. In all cases, the value of 18 was found to be less than 0.000015. Based on these results, we concluded that it was reasonable to consider solving a reduced system of equations in which all terms involving Fn, n ~ 8, were disregarded. Thus, we set L = min(M,7), and content ourselves with solving the following system of equations: L L-2 L-i (m(O) = bm(0) + pm(0) r" 2 + 2 A pm(i) X Fn n+i, (3.12a) n=2 i=l n=2 L-k-1 pm(k) = bm(k) + X 'm(n)rn+k+l, L-2 k > 0, with (3.12b) n=l-k L-1 k bm(0) = g2(n), bm(k) = Xm(M-l-k+n) gm(n), k=1,...,L-2, and (3.13) n=l n=l L-1 gm(n) = S Xm(M-l-k+n) Fk+1. (3.14) k=n We will refer to Algorithm 2 with L = min(M,7), as Algorithm 2 at "level 7." Note that Algorithm 2 at level 7 returns the exact numerical solution if the number of nodes is 7 (or less). We can rewrite the system of equations given by (3.12) as Bx = b, where, for Algorithm 2 at level 7, x = (pm(0),..., pm(5)), and b = (bm(0),...,bm(5)). The matrix B is presented in Figure 3.1. n 0 1 2 3 4 5 1 - 7 r2 -2 6 nrn+l -25 2Fnn+2 -24 2rnn+3 -23 2rnrn+4 -22 Fnrn+5 _n=2 n n=2 n=2 =2 n=2 n=2 1 -r2 1-r3 -r4 -F5 -r6 -r7 2 -r3 -(r2+r4) l-F5 -r6 -F7 0 3 -F4 -(F3+r5) -(r2+r6) 1-r7 0 0 4 -rs -(F4+r6) -(r3+r7) -r2 1 0 5 -F6 -(F5+r7) -F4 -F3 -r2 1 Figure 3.1. The Matrix, B, for the system of equations at level 7 12

Lemma 3.3. The solution to the system of equations specified by (3.12) through (3.14) will result in a lower bound on the exact value of 5m(0), when M > L. Proof: The original system of equations (equations 3.9a and 3.9b) is in the form A(pm = bm, where A is an M-1 x M-1 matrix consists of terms involving Ins. It follows easily, from Lemma 3.1, that the matrix A is an M-matrix, i.e., a matrix that is diagonal-dominant, with it's offdiagonal elements non-positive. It is well known (Fiedler and Ptak 1962) that the inverse of such a matrix is positive. Let the inverse be denoted by A-1. We also know that the vector bm is positive, and so we know that the true solution pmn = A- bm > 0. Note that the matrix is an (M-1) x (M-1) square matrix and that it's row indices are 0,...,M-2. Let us denote the true solution by the vector (x ) where x denotes the first L-1 elements m(O) through im(L-2), and y denotes the elements qm(L-1) through qm(M-2). Analgously, let b denote the vector bm, obtained from equation (3.10), where b1 denotes elements bm(O) through bm(L-2), and b2 denotes elements bm(L-2) through bm(M-2). Denote by B, the square submatrix consisting of the first L-1 columns and rows of A, and let C denote the matrix consisting of all the remaining elements in the first L-1 rows of A. Thus, Bx + Cy = b. Since matrix A is an Mmatrix, we know that the matrix C is nonpositive, and so we must have Bx = b1 + (-Cy) > b1. Since we know that B-1 is positive (because A-1 is positive) the solution to x = B-bl represents a lower bound on the true solution. Furthermore, equation (3.13) only computes a lower bound (say, b,) on the true value of b. Thus, the resulting x is a lower bound on the exact solution. U For L = min(M,7), the expression for 0m(0) can be obtained in a more explicit form as follows. We first exploit the structure of matrix B, to express jm(i), i=2,...,5, only in terms ijm(O) and (m(1) as follows: (Pm(i) = um,i + vi m(O) + Wi m(l), i = 2,...,5, (3.15) where ) 6+I26 7 bm(3) + Um,2(F2 + T6) um,2 = (bm(2) + bm(34)7)/A, Um,3 l-r7 1-F7 Um,4 = bm(4) + Um,2(r3+r7) + Um,3r2, um,5 = bm(5) + um,2F4 + um,3r3 + Um,4r2, (3.16a) + 6+1r27 r4 + v2(r2 + F6) v2 = (F3 + r4 --- + 5r7)/A, v3 = 1-F7 1-r7 V4 == F5 + v2(r3+F7) + V3F2, v5 = F6 + v2F4 + v313 + v4r2, (3.16b) 13

W2 = (1'2+F14+(17'3+I75) 6+2]7+ r(4I6)A W4 =174+176 + W2(173+177) + w31F2, 173+175 + w2(17'2+176) W3 = iv W5 =175+1F77+ W2174 + W31F3 + W4172, (3. 16c) and A = '~~~~~I6 + r2l'7_ 1-177 (3.17) We now express 'm(l) in terms of ~m(O), (refer to the second row of Figure 3. 1), and this gives: =5Ml cm +.~ a ~+~mYLJ a, where (3.18) Cm = ao = bm(l) + F4Um,2 + IF5Um,3 + IF6Um,4 + IF7Um,5,) 172 + 174v2 + 175v3 + 176v4 + F77v5, (3. 19a) (3.19b) and al = 1 - (173 + 174w2 + J75W3 + 176w4 + 177w5). (3. 19c) From equations (3.15) through (3. 19c), we obtain the solution for ~m(0) presented below, after some simplification, as Lemma 3.4. Lemma 3.4. For Algorithm 2 at level 7, (pm(O) is obtained as M-1 (Pm(0) = X Xmi j=1 bm(0) - 2( bM(1)2 + bm(2)173 -bm(3) F4r- um,2 v2 A- cm -o 1-17l a,1 (3.20) 1-X17 2+ 2 ][-2+172 4 - V22A a1 n=2 n 2+ 3-1r We can, of course, simplify the denominator of equation (3.20) somewhat further. We have left it in the current form to indicate how the numerator and denominator are related. (If we replace bm(n) by 1+ for all the terms within the paranthesis in the numerator we obtain the terms within the paranthesis in the denominator of the equation.) If M < 7, then some of these terms are not computed. For example, if M = 4, then A = 1, and we only compute um,,i, vi and wi for i=2, and set the terms involving 175, 176 and 177 = 0, in equation (3.20). Similarly, for M=3, um,j, vi, and wi = 0 for all i, cm = bm(l), ao = I72, al = 1 -173, and bm(2) = 0 (note that k ranges only from 1 through M-2 in equation 3. 10). Thus, for M=3: bm,(0) + 2bm(1)172173/(l1173) %(O) =1 _ - 27 _]- 12 - 2 [217-13I( 1-13)' (3.21) 14

From equations (3.10) and (3.11), bm(0) = g2(1) + g2(2) = (Xm(2)F2 + Xm(l)F3)2 + (Xm(2)r3)2, and bm(l) = Xm(2) (Xm(2)F2 + Xm(l)r3). If we substitute these values in equation (3.21), we obtain the same expression for ~(m(O) as in equation (3.8a). Algorithm 2: 0. Initialization: Compute a) 6m, m=l,...,M (equation 2.10); b) Fn, n = 2,...,L, (equation (3.3); c) A (equation 3.17); d) vi and wi, i = 2,..., L-2, (equations 3.16b and 3.16c); and e) ao and al (equations 3.19b and 3.19c). Evaluate the denominator of equation (3.20). 1. For every node, k, whose mean waiting time is desired, renumber the nodes so that node k is node "1", node k+l is node "2",...., and execute steps 2 and 3. 2. For each m=l,...,M, do steps a) through d). a) Compute m(j), j = 1,...,M-1, (equation 3.1). b) Compute bm(j), j = 1,...,L-2, using (equation 3.13). c) Compute um,i, i=2,...,L-2, (equation 3.16a), and Cm (equation 3.19a). d) Compute (pm(O) (equation 3.20). 3. Obtain the mean waiting time using equation (2.25). The computational complexity of Algorithm 2 is obtained as follows. Evaluating rn requires 0(M2) operations, regardless of the level L. The rest of the computations in Step 0 requires 0(1) time. These computations are incurred regardless of the number of nodes being evaluated. Step 2 requires O(M2) operations each time it is carried out. The major computational effort in this step arises from steps 2a) and 2d). Hence, the overall computational effort is O(M2) for each node that is evaluated for its mean waiting time. In terms of storage requirements, the algorithm only needs one 2-dimensional array of size M1 x M-1 to store the Xm(J) terms, and several one-dimensional arrays of size M or less. Computational Remark: If M is large, then there can be significant savings in computation if we make use of equations (3.15) and (3.18) as follows. From equation (3.1), we obtain 21 ) P 2?O )) ~m(+i) =,jM Xm((i) = Pm+?j-1XM (5 + Xm(i-)Y = Pm+l ^2 + m() + m(M0-) + 2 jM Xm(J)Xm(j-1) = ~ m +m J)+2A m( = P Pm(0)( + ) + 2m(l) + X(M-l) (3.22) 15

The term,m(l) is obtained, from equation (3.18), as pm(l) = (Cm + aopm(O))/al. Thus, once we compute {m(O), we directly obtain?m+1(0) with just a few arithmetic operations, from equations (3.18) and (3.22). We can carry this approach further as follows. The term im+i(l) is obtained, similar to equation (3.22), as ml) = m+m1 +m(O)+~m(2) +l(l) = Pmljm()( +^-) + + Xm(M-)[Xm(M-l)/m + Xm(M-2)] (3.23) We compute (m(2) from equation (3.15), and use it in equation (3.23) to evaluate 5m+1(l). Having computed (m+i(0) and Fm+i(l), we obtain m+2(0) using equation (3.22) once again. It is, of course, possible to carry this idea further. Thus, we can compute 5m+3(0) if we know i5m+2(0) and (m+2(1). Note that, in order to compute Sm+2(1), we must compute (m+i(2), and so on. The savings in computation become insignificant after some time. In our implementation of Algorithm 2, we compute om(O), for m = 1, 4, 7,..., from equation (3.20), and then compute 5im(O) for m = 2, 3, 5, 6, 8, 9,... using the above approach. We can show that this approach still results in a lower bound on the exact mean waiting times. U The next section on computational experience demonstrates the accuracy of Algorithm 2 at level 7. We remind the reader that Algorithm 2 at level 7 obtains the exact mean waiting times if M<7. 4. Computational Experience We conducted a number of experiments with randomly generated networks, to test the performance of the two algorithms. We were interested in determining a) the execution times for Algorithm 1, and b) the execution times and accuracy of Algorithm 2. We varied the number of nodes from 3 to 50, and the server utilization, p, from 0.30 to 0.996. We report on some experimental results in this section, which appear to best represent the performance of these algorithms. The code was implemented in Pascal and the experiments were conducted on an IBM 9021-720 running the MTS operating system. The iterative algorithm was terminated as soon as the largest Xm(J) value obtained during an iteration dropped below 0.0001. Based on extensive comparisons with exact results, we found that this termination criterion gave accurate results. (The results were accurate upto 4 decimal places, and resulted in a percentage error of less than 0.0001%.) We observed that Algorithm 1 typically ran very quickly when the value of p was less than 0.90. We also observed that Algorithm 2 was uniformly extremely accurate and obtained mean waiting times very quickly. We ran Algorithm 2 at level 6 for p < 0.90, and this always gave the exact results (accurate upto 6 decimal places). For p = 0.99, we ran Algorithm 2 at level 7, and the mean waiting times obtained 16

from Algorithm 2, when compared to the corresponding exact mean waiting times, had a maximum percentage deviation of less than 0.0025% for all the experiments we conducted. At this point, we must emphasize that if even greater accuracy is required, then we only have to set L = min(M,8), and solve the "level 8" algorithm (this will involve the following additional computations: a onetime computation of the terms v6 and w6, and subsequent computations of bm(6) and um,6). In Table 1, we report on the maximum number of iterations and the running time of Algorithm 1, and the accuracy and running times for Algorithm 2. The results are reported for the following values of p: 0.30, 0.50, 0.75, 0.90, and 0.99, for M = 12, 24, 48 and 96 nodes. For each example considered, the algorithms were executed to obtain the mean waiting time for: a) 25% of the nodes, b) 50% of the nodes, c) 75% of the nodes, and d) 100% of the nodes. As expected, the execution time for Algorithm 1 increased in a linear manner with the number of nodes that were evaluated for mean waiting times. The execution time for Algorithm 2 was measured as the sum of two components: the time taken to complete the initialization phase, and the time taken for the second phase, namely to compute the mean waiting times. As expected, the initialization time, for a given p and M, was the same regardless of the number of nodes evaluated while the execution time for the second phase grew in a linear manner with the number of nodes evaluated (subject to insignificant variations in CPU times). Table 1 reports the execution times for these two phases. Table 1 indicates the growth of the execution times for the two algorithms grow with increasing values for p and M. Interestingly enough, Algorithm 1 uniformly appears to take roughly the same number of iterations for a given p, regardless of the number of nodes. In general, Algorithm 1 performs very well, so long as p is about 0.75 or less. For higher values of p, the performance of Algorithm 1 steadily degrades. Since it consistently takes about 28 to 35 iterations for p = 0.90, and since each iteration requires O(M) operations, this suggests that if p < 0.90 and M is larger than around 40, Algorithm 1 would require less than O(M2) operations to compute pm(O), m=1l,...,M, and, therefore, require less than O(M2) operations to find the mean waiting time at a node. This also suggests that if Algorithm 1 is used to compute all the mean waiting times, it would probably outperform an O(M3) algorithm that computes all the mean waiting times, for large values of M (say, for M > 40), and p < 0.90. Algorithm 2, of course, does not depend on the value of p, but on the other hand, it requires more time for increasing values of M. As far as the accuracy of Algorithm 2 is concerned, it never gave an error greater than 0.0022%. We emphasize that this is a percentage difference, namely that the relative error is multiplied by 100. We ran Algorithm at level 6 for p < 0.99, and it gave the exact values for the mean waiting times (correct to 4 decimal places) every time. (This explains the marginal increase in the execution time observed for p = 0.99.) As expected, from Table 1, 17

the observed growth in CPU time required by Algorithm 2 is of O(M3)to compute the mean waiting times at all nodes, and of 0(M2) to compute the mean waiting time at a single node. Table 1. Performance of Algorithms 1 and 2. Algorithm 1 Algorithm 2 l M Max. # Tot. cpu cpu time/ Max. Tot. cpu Initializn. cpu time/ iterations time (ms) node (ms) error (%)time(ms) time(ms)2 node(ms)3 0.30 5 10.04 0.84 0.0000 15.62 2.55 1.09 0.50 7 11.68 0.97 0.0000 15.67 2.61 1.09 12 0.75 14 17.60 1.47 0.0000 15.64 2.54 1.09 0.90 35 35.71 2.98 0.0000 15.54 2.64 1.08 0.99 368 316.55 41.66 0.0005 16.43 2.60 1.15 0.30 5 48.51 2.02 0.0000 95.53 10.17 3.56 0.50 7 55.54 2.31 0.0000 95.20 10.33 3.54 24 0.75 13 76.00 3.17 0.0000 94.79 10.14 3.53 0.90 34 146.87 6.12 0.0000 95.48 10.28 3.55 0.99 356 1,242.42 51.77 0.0022 99.76 10.37 3.72 0.30 4 275.48 5.74 0.0000 675.94 60.71 12.82 0.50 6 305.16 6.36 0.0000 675.41 60.74 12.81 48 0.75 13 392.82 8.18 0.0000 673.31 60.62 12.76 0.90 31 649.13 13.52 0.0000 676.33 60.74 12.82 0.99 331 4,738.71 98.72 0.0018 691.15 61.28 13.12 0.30 4 1,894.48 19.73 0.0000 5,081.83 437.11 48.38 0.50 6 2,002.55 20.86 0.0000 5,106.79 438.07 48.63 96 0.75 11 2,293.96 23.90 0.0000 5,080.99 440.98 48.33 0.90 28 3,235.45 33.70 0.0000 5,127.91 439.43 48.84 0.99 302 18,212.95 189.72 0.0000 5,158.67 440.15 49.15 We present some examples to give the reader some feel for actual numerical values. These examples also serve to illustrate the accuracy of Algorithm 2 at level 7, since all of them have very large p values. In one example with 12 nodes, we set p = 0.996. The percentage error for the largest deviation found in this example, was equal to 0.0016%. 1 Algorithm 2 was always run at level 6 for p < 0.99. 2 Initialization Time is incurred regardless of the number of nodes evaluated for their mean waiting times. 3 Excluding the initialization time 18

Example 1. The number of nodes = 5, with arrival rates 0.2, 0.8, 0.4, 0.2, and 0.1. The mean service times are 0.5, 0.5, 0.4, 1.0 and 1.2, and the second moments of the service times are 0.5, 0.5, 0.3, 2.0, and 2.9. The mean switchover times are 0.5, 0.5, 0.4, 1.0, and 1.2. The variances of the switchover times are equal to their means. For this problem, p = 0.98. The mean waiting times are W1 = 121.0880, W2 = 80.7446, W3 = 113.3191, W4 = 107.7545, W5 = 118.3033. Example 2. The number of nodes = 7, with arrival rates 0.2, 0.6, 0.4, 0.2, 0.1, 0.1, and 0.1. The service times have means 0.5, 0.5, 0.4, 1.0, 0.6, 0.6, and 1.1, and second moments 0.3, 0.3, 0.2, 1.5, 0.8, 0.8, and 1.8. The mean switchover times are 0.5, 0.5, 0.8, 1.0, 0.6, 0.6, and 1.1. The variances of the switchover times are equal to their means. For this problem, p = 0.99. The mean waiting times are W1 = 283.0562, W2 = 220.2503, W3 = 264.4854, W4 = 251.7250, W5 = 295.7516, W6 = 295.7251, W7 = 279.8502. Example 3. The number of nodes = 10, with arrival rate equal to 0.75 at node 1, and 0.2 at all the other nodes. The mean service times are bl = 0.50, b2 =... = b5 = 0.40, b6 = b7 = b8 = 0.3, and b9 = blo = 0.2. The second moments of the service times are as follows: 0.5 at node 1, 0.3 at nodes 2 through 5, 0.2 at nodes 6 through 8, and 0.05 at nodes 9 and 10. The mean switchover times are 0.5 at node 1, 0.2 at nodes 2 through 5, 0.3 at nodes 6 through 8, and 0.5 at nodes 9 and 10. The variances of the switchover times are 0.25 at node 1, 0.15 at nodes 2 through 5, 0.20 at nodes 6 through 8, and 0.25 at nodes 9 and 10. The value of p = 0.955. The mean waiting times obtained by Algorithm 2 matched the exact mean waiting times (correct to 4 decimal places). The mean waiting times are W1 = 28.8749, W2 = 42.5150, W3 = 42.5305, W4 = 42.5481, W5 = 42.5684, W6 = 43.5228, W7 = 43.5544, W8 = 43.5905, W9 = 44.5781, and W10 = 44.6425. Example 4. The number of nodes = 12, with Xi = 0.01 i, i = 1,...,12. The mean (second moment) of the service time is equal to 1.60 (2.00) at odd-numbered nodes. The mean (second moment) of the service time is equal to 1.00 (1.50) at even-numbered nodes. The switchover times have identical means and identical variances at all nodes, and they are, respectively, 0.50 and 0.25. The p value for this example is 0.996. (Algorithm 1 required over 1200 iterations to converge to the exact values.) Comparing the exact mean waiting times with the mean waiting times obtained by Algorithm 2, the largest deviation from the exact mean waiting time occurred at node 1. At node 1, the exact mean waiting time is 924.8319. The corresponding mean waiting time obtained by Algorithm 2 was 924.8174. This resulted in an absolute deviation of 0.0145, and a percentage deviation of 0.0016%. 19

Example 5. The number of nodes = 24. ki = 0.4 for nodes 1 and 13, Xi = 0.20 for nodes 2 and 14. For all other nodes, Xi = 0.01. The mean (second moment) of the service time is equal to 0.80 (2.00) at odd-numbered nodes, and equals 0.50 (1.50) at even-numbered nodes. The mean switchover times are all equal to 0.50, and the variances of the switchover times are all equal to 0.25. The value of p = 0.97. The mean waiting times are 173.8729 at nodes 1 and 13, and 230.1439 at nodes 2 and 14. The mean waiting times at nodes 3 through 12 range from 253.6835 to 254.5273. The mean waiting times at nodes i and i+14 are identical (upto 4 decimal places) for i=l,...,10. The mean waiting times obtained by Algorithm 2 matched the exact mean waiting times upto 4 decimal places. Example 6. The number of nodes = 48. Xi = 1.50 for nodes 1 and 2, 1.00 for nodes 3 and 4, 0.50 for nodes 5 through 10, 0.45 for nodes 11 through 22, 0.30 for nodes 23 through 34, and 0.20 for all other nodes. The mean (second moment) of the service time is 0.05 (0.05) at all nodes. The mean and variance of the switchover times are equal to 0.10 and 0.08, respectively, at all nodes. The value of p = 0.99. The exact mean waiting times at nodes 1 through 4 were W1 = 269.6124, W2 = 269.6075, W3 = 276.8928, and W4 = 276.8912. At nodes 5 through 10, Wi ranged from 284.1779 to 284.1786. At nodes 11 through 22, Wi ranged from 284.9075 to 284.9103. At nodes 23 through 34, Wi ranged from 287.0966 to 287.1007. At nodes 35 through 48, the mean waiting times ranged from 288.5583 to 288.5627. The mean waiting times obtained by Algorithm 2 were compared with the exact mean waiting times. The largest absolute deviation from the exact mean waiting time occurred at node 34. The exact mean waiting time at node 34 was 287.1007, while the mean waiting time obtained by Algorithm 2 was 287.0962, resulting in an absolute error of 0.0045, and a percentage deviation of about 0.0016%. In all these examples, the mean waiting times for nodes with identical characteristics are very nearly equal. This indicates that for systems where many nodes have similar characteristics (as is usually the case for large systems), it may be more appropriate, at least during the design stage of the network when many alternate configurations may need to be considered, to evaluate mean waiting times only for some select few nodes. For instance, in example 7, it may be adequate to evaluate the mean waiting times at nodes 1, 3, 5, 23, and 35. It is especially in such situations that the IS technique has a distinct advantage over other techniques where all nodes have to be evaluated simultaneously. 20

4.1. Comparison with the algorithm of Sarkar and Zangwill We compared the execution times for Algorithms 1 and 2 with the execution time for Sarkar and Zangwill's algorithm (the S&Z algorithm). Before we present the results, we make the disclaimer that while we tried to implement the S&Z algorithm as efficiently as possible, there is always a possibility that one could improve its implementation even further. The same disclaimer holds for the implementations of Algorithms 1 and 2. (Note that if Algorithm 2 is used to compute the mean waiting times at all nodes, it is an O(M3) algorithm, as is the S&Z algorithm.) We ran the S&Z algorithm using the same data that was used to generate Table 1. The algorithm was run on the same IBM 9021-720. The S&Z algorithm was coded in Fortran, which made it easier to call LNPACK (a library of subroutines for solving linear systems of equations, written in Fortran). The execution times for the three algorithms are presented graphically in Figures 4.1 through 4.4, for M = 12, 24, 48, and 96, respectively. These figures show the increase in execution times of the Algorithms 1 and 2, with the number of nodes that are evaluated for their mean waiting times. Since the execution time of Algorithm 1 increases with p, for each M, we present the execution times of Algorithm 1 for p = 0.50, 0.75 and 0.90. We do not present the execution time of Algorithm 1 at p = 0.99, since it is relatively very large. Unlike Algorithm 1, Algorithm 2 and the S&Z algorithm do not depend on p. Therefore we do not present the execution times of Algorithm 2 and the S&Z algorithm for different p values; the execution times presented are the averages of the execution times for the different p values. Based on these figures, Algorithm 2 would be preferred over the S&Z algorithm if mean waiting times are desired at less than about 60% of the nodes. Algorithm 2 takes about 60% more time than the S&Z algorithm, to evaluate all the mean waiting times, using a single computer. However, it is noted that Algorithm 2, following the initialization phase, computes the mean waiting times for individual nodes independently. Therefore, using Algorithm 2, one could obtain mean waiting times at all nodes faster than the S&Z algorithm if another process was spawned on a similar computer, once the initialization phase ends. The two processes could work in parallel, each computing mean waiting times at, say, half the nodes, reducing the execution time for the second phase of Algorithm 2 by half. The developments in distributed computing makes such an approach very feasible. The performance of Algorithm 1 steadily improves (for p < 0.90) as the number of nodes increase, relative to the performance of the other two algorithms. For M = 96 nodes, Algorithm 1 computes all the mean waiting times faster than the S&Z algorithm even at p = 0.90. This agrees with an earlier observation we made about the relative efficiency of Algorithm 1 as compared to an O(M3) algorithm for large M and p less than 0.90. 21

40 m E u (A I Wr *o E f.3 r. u 30 20 10 100% Fraction of nodes evaluated Figure 4.1: Execution times for M = 12 nodes 100% Fraction of nodes evaluated Figure 4.2: Execution times for M = 24 nodes 22

O13, Algorithm 1, Rho=0.90 O 600 -0- Algorithm 2 0 Sarkar&Zangwill * 400 Q) P1200 0%/ 25% 50%1c 75% 100%/ Fraction of nodes evaluated Figure 4.3: Execution times for M = 48 nodes 6000 Algorithm 1, Rho=0.50 5000 --- Algorithm 1, Rho=0.75 U0- Algorithm 1, Rho=0.90 0- Algorithm 2 0 Q) 4000 N Sarkar&Zangwill E 3000 p2000 Q 1000 0 0%/0 25% 50% 75% 100%0 Fraction of nodes evaluated Figure 4.4: Execution times for M = 96 nodes 23

5. The Gated Service Polling System The analysis for the gated service polling system is very similar to the analysis for the exhaustive service polling system. As before, we obtain the mean waiting times at the nodes, after computing the first and second factorial moments of the number of customers present at the node when the server polls it. The mean waiting time is given in terms of fm and f(2) as (Takagi 1986): f(2) l+pm Wm-2fm n (5.1) As before, we first obtain fi and f(). To compute these terms, we require a redefinition of the recursive funtion y as follows: Ym(O) = Zm, m= 1,...,M, (5.2a) m(j) = Pm( [k- XkYk(-1)] + [k-XkYk(J)]), j>O, m=l,...,M. (5.2b) k<m k>m The term Cmo(j) is still defined by equation (2.5), with the redefined Ym terms. Proceeding exactly as before, we obtain: 00 F1(y1(0),...,yM(O)) = n (m am(j)) (5.3) Equations (2.8), (2.9), and (2.12) still hold. Equation (2.11), however, is modified as follows: m(j) = Pm [ k(j-l) + W Wk(j)], j>O, and (5.4a) k<m k>m m)(J= Pm[ v (2)(j-1) + P N-(2k)(J)] + nb(2) [ W k(j-l) + VIk(j)]2 j > 0. (5.4b) k<m k>m k<m k>m Following a similar analysis as used to obtain fl for the exhaustive service case, we get fi = 1. (5.5) 1-p The term f(2) is given by the same expression as for the exhaustive service case: 00 f(2) = f + Sm y Yk + Var(Sm) [ X Vk(j) + Vk(j+l) ]2 (5.6) m k m j=O k<m k>m where y m = m.) (mb(2) 00 m m (P)Pm j(5.7) m ( 1 -p)Pm2 j=m 24

However, unlike the exhaustive service case, where [Xk~m 'Ik(J) + ~k>m lVk(+l) i2 depended on node m+1 (refer to equation 2.20), for the gated service case we observe, from equation (5.4a), that this term simply evaluates to ~Vmoj+l). Hence, equation (5.6) is rewritten as: 52' f2 ~ Xmb(m2) 200 + Var(Sm) 20,= + s~ V +. (5.8) m (1_p)p2~ j=1 m Pm2 j Define, as before, the terms Xmoj), (pm(n), and om(n) as xmOi) = VA4m~) _ 2OmIYmA) 1(5.9) XI XI az< 00 00 (Pm(n) = XmOi)XmOj-n), and m(n)= XmO)Xm(J-n). (5.10) j=n+l j=M The expression for f~2) for the gated service polling system is given below as Lemma 5. 1. Lemma 5.1: f2) = f2+ XM(m0( ~~ Var(Sm)) 1 1m 1-p PM 2 PM2 From Lemma 5. 1, and equations (5. 1) and (5.5), we obtain: Theorem 5.2: The mean waiting time at node 1 for the gated service polling system is: =, 2 1-pi +(l+p )X 2P(O (2 2s Vari(SM)) (5.11) 5.1. Computing (pm(O) for the gated service system As before, we present two algorithms for computing (Pm(O). Let 6m = PM/(1+Pm). (5.12) Similar to the exhaustive service case, the Xmoi) terms are computed using the following expression, which is obtained from equations (5.4a) and (5.9): XMOl) = PM; Xm(1) = (Xm+i(1)I 8n+D Pm' m=1,.. - 1 M-, (5.13a) and. for j> 1, XM(j) = (x1(-1)161 - X1(i-2))pM; Xmoj) = (Xm+loj)16m+l - Xm+1(.hl))pm.) m=1,....,M-1.(5.13b) 25

The iterative algorithm to compute pm(O) is identical to the one we developed for the exhaustive service case, and so we do not repeat it here. The algorithm which computes Tm(O) by solving a system of equations is also very similar to the corresponding one for the exhaustive service case. Some differences, however, exist, which we clarify here. Define the terms Ki(m and Qn as: m-l K( = Pm, m=1,...,M; K( = Pm E K(i1); (5.14) j=i-1 M QM = I Pm, and n = 2m+1 +, K, 1 < m < M. (5.15) m k=m Intuitively, Q1, for example, denotes the sum of all possible combinations of unique Pm terms taken singly (i.e., pl+.. +pM) + the sum of all possible combinations of 2 unique Pm terms +... + the sum of all possible combinations of M unique Pm terms. The terms, rn are defined as: M rl = ai, and Fn = (-l)n+1 (m22) Q n=2,... M. (5.16) m=n Defining the terms Om(j) as follows: OM(1) = PM, OM(j) =0, j>, (5.17a) Om(j) = Pm(Om+l(j)/6m+l - Om+(j-1)), m > 1, j = 1,...,M-1, (5.17b) (j) = M=2m, j = 1,...,M-l, (5.17c) / ~m=2 we can obtain an expression for Xm(j) in terms of Xm(j-n) stated below as Lemma 5.3. The proof of Lemma 5.3 is by induction on M, and is omitted. Lemma 5.3 M Xm(j) = E Xm(J-n) Fn, j 2 M, (5.18) n=l with Xm(i) = jXm(-n)rn + Om(j), j<M-l (5.19) n=l Note that the rn terms in equation (5.18) for the gated service polling system range from n=l, through M, as opposed to the corresponding equation (3.4) for the exhaustive service polling system, where the rn terms range from 2 through M. Based on equation (5.18), we can derive the following system of equations for obtaining q5(O): 26

M M-1 M-i in(0) = bm(O) + ~m(0) r 2 + 2 E xm(i) E rn r+i n=l i=l n=l (5.20a) M-k pm(k) = bm(k) + ~m(n)rn+k, n=l-k M bm(O) = 2 g2(n), bm(k n=l M-1 >k>0, (5.20b) k:) = Xm(M-l-k+n) gm(n), k = 1,...,M-1, n=l (5.21) where gm(n) is computed, recursively, from M gm(n) = E Xm(M-l-k+n) Ik. (5.22) k=n As before, in equation (5.20b) we use the convention that im(-k) = (m(k). This system of equations is, once again, of the form A(pm = bm, and can be solved for the mean waiting times. For the gated service polling system with M nodes, we have M equations in M unknowns, namely, the Sm(k)s, k=O,...,M-l. As before, the terms Fn decrease at an exponential rate with n. However, we were unable to prove that by ignoring terms involving rn, say, for n > 7, we would obtain a lower bound on the mean waiting time. For the gated service polling system, we observe that if we only consider terms upto '6, ignoring terms involving Fn, for n > 7 (the resulting solution is termed the solution to Algorithm 2 at "level 6"), excellent results are obtained if we also disregard all Q terms for n 2 7, in computing the F's. Lemma 5.4 presents the level 6 system of equations. It may be observed that they are strikingly similar to the level 7 equations for the exhaustive service polling system. Also, note that Algorithm 2 executed at level 6 returns the exact numerical results for a system with upto 6 nodes. Lemma 5.4. For Algorithm 2 at level 6, pm(O) is obtained as Er3 6 ao'\ M-1 bm(O) - 2(bm(l)1F + bm(2)r2 - bm(3) 6 - Um,2v2A - cma) (Pm(O) = m2(j) +, r(5.23) j=l 62r6 1j i 2 + 2 2+3: -v 2A- 2 n=l 1-F6 al J where Cm = bm(l) + r3um,2 + r4um,3 + r5Um,4 + r6Um,5, ao = rl + r3v2 + r4v3 + r5V4 + r6v5, al = 1 - (F2 + r3w2 + r4w3 + r5w4 + 6ws5), and (5.24a) (5.24b) (5.24c) 27

A = 1 - r4 - (FI + r)5f --- — 6 - r6(r2 + r6). 1-16r The terms ui, vi, and wi, for i = 2,...,5 are given below: (5.24d) Um,2 = (bm(2) + bm(3) 6 + b(4)r6)/, 1-F6 Um,4 = bm(4) + Um,2(r2+r6) + Um,3rl, V2 = (r2 + r3-5+l6 + r4r6)/A, I-F6 V4 = r4 + v2(r2+r6) + V3rl, w2 = (r1l+3+(12+r4)I5+I1F6 + r6(r3+r5))/A, 1-W6 W4 = 3+1r5 + W2(12+1r6) + W311, bm(3) + Um,2(r1 + Fs) Um,3 - 1- 6 1-F6 Um,5 = bm(5) + Um,2F3 + Um,3r2 + Um,41l. (5.25) r3 + v2(rl + r5) V3 = 1 -V5 = F5 + v2r3 + V41rl. (5.26) r2+r4 + w2(rl+r5) W3 = 1-F6 ws5 = r4 + r6 + w2r3 + w3F2 + W4F1. (5.27) 5.2. Computational Experience Both Algorithm 1 (the iterative algorithm) and Algorithm 2 performed strikingly similar to the corresponding algorithms for the exhaustive service case. We do not report on their performance here. To get some feel for numerical values, we revisit some of the examples given in section 4, and present the results obtained. We present the results for Examples 3 through 6 of section 4. These are the examples with a relatively large number of nodes, and high utilizations. For the gated service system, we can again make use of the "Computational Remark" in section 3. For all cases reported, we ran Algorithm 2 at level 6. Example 3 revisited. Both algorithms obtained the following mean waiting times: W1 = 58.9669, W2 = 46.2956, W3 = 46.2918, W4 = 46.2874, W5 = 46.2822, W6 = 45.4192, W7 = 45.4182, W8 = 45.4171, W9 = 44.5587, W10 = 44.5788. Example 4 revisited. The mean waiting times ranged from 916.6964 at node 1 to 1061.0240 at node 11. Algorithm 2 obtained the exact mean waiting times (correct to 4 decimal places). In fact, it obtained the exact mean waiting times even at level 5. Example 5 revisited. The mean waiting times ranged from 235.8342 at node 24 to 309.7431 at node 1. Algorithm 2 obtained the exact mean waiting times (correct to 4 decimal places). Example 6 revisited. The mean waiting times ranged from 291.2991 at node 48, to 310.0504 at node 2. Algorithm 2 obtained the exact mean waiting times (correct to 4 decimal places) in 861.59 milliseconds. 28

6. Polling systems with a mix of gated and exhaustive service disciplines In this section, we present the iterative algorithm for a polling system in which some of the nodes may be serviced according to the exhaustive service discipline, while other nodes may be served using the gated service discipline. The development of an algorithm similar to Algorithm 2, for this case, is left as a topic for possible future research. The analysis, again, proceeds in a manner similar to the analysis for either the pure exhaustive service polling system, or for the pure gated service polling system. Let G denote the set of nodes using the gated service discipline, and E denote the set of nodes using the exhaustive service discipline. As before, we obtain the mean waiting times at the nodes, after computing the first and second factorial moments of the number of customers present at the node when the server polls it. The mean waiting time at a node, m, is given in terms of the fm and f(2) as: f(m2 1l+pm f(2) Xmb(m2) Wm - 2fm l m m G, and Wm = m +, m e m E. (6.1) km 2Xmfm 22(1-pm) To compute fi and f(2), we redefine the recursive funtion y as follows: Ym(O) = zm, m = 1,...,M, (6.2a) Ym() = P1m( [Xk-kYk(-1)] + [Xk-xkYk(j)]), j > 0, m G, (6.2b) k<m k>m Ym() Tm( = T [Xk-xkYk(J-1)] + [Xk-Xk-k(j)]), j> 0, m E E. (6.2c) k<m k>m The term am(J) is still defined by equation (2.5), with the Ym terms as defined above. Proceeding exactly as before, we obtain: 00 Fl(Y1(O),....,yM(O)) = 1n (rn m(i) (6.3) Using the same definitions as before, equations (2.8), (2.9), and (2.12) still hold. Equation (2.11) is modified as follows. Forj > 0, rnm() = Pm [ Vk(j-1) + Vk(J)], m E G, (6.4a) k<m k>m /m(j) = m [ Z Vk(J-1) + Vk(J)], m E E, (6.4b) k<m k>m (2)0) = Pm [ -1)+ )] + [ v(k( —1) + d)] +Vb( [ k(J)]2, m E G, (6.4c) k<m k>m k<m k>m 29

Xmb(2) m( MI +l+= a5( (j<m vk(-1)+ k()], m e E, (6.4d) k<m ( m)k<m k>m Following a similar analysis as used to obtain fi for either the pure exhaustive service case, or the pure gated service case, we get fl = X -, if node 1 E G, and fl = Xls, if node 1 E E, (6.5) 1-p 1-p The term f(2) has an identical expression as the one for eithereither the pure exhaustive or the pure gated case: 00 f(2) = f + E Sm Yk + Var(Sm) [ w k(j)+ k(i+l)]2 (6.6) m k m j=O k<m k>m where 3Xmb(2) 00 EY~m = m ~ (f -(2 E i). 2 (6.7) m m ~ (^Pm j=l Unlike either the pure exhaustive service case, or the pure gated service case, we must now take care when substituting for the term [Ik<m k(j) +,k>m vkoi+l)]2 For the pure gated service case, we used equation (5.4a) to write this term as fm(+l1)/p2, and we can still do that here. For the pure exhaustive service case, we used equations (2.20) and (2.21) to set this term equal to xVm+(j+l)/pm21, for j > 0. However, this assumes that node m+l also uses the exhaustive service discipline, which may no longer be true for the system we are considering. Nevertheless, we can still simplify this term for m E E, as [Xk:m Vk(J) + Xk>m fk(J+ 1)] = [Nfm(j) + Vm(j+1)/6m]2, where 8m = pJm(l-pm). Therefore, if we let Xm() = 1m() (6.8) XI (Pm(O) = E Xm(j) and m(0) = [Xm(i) + Xm(i+l)/m]2, (6.9) j=1 j=o from equations (6.6) and (6.7) we obtain Lemma 6.1: Lemma 6.1: Xmb(m2) Var(Sm) f(2) = + 2 + 1 s ) + ( (pm(0) + ~ Var(Sm) m(O). Fo Lma61 a eutn2 ( 1 2 w o in m (1P)Pm me G Pm meE From Lemma 6.1 and equations (6.1) and (6.5), we obtain: 30

Theorem 6.2: The mean waiting time at node 1 for the polling system with the mixed service discipline is: l+p + (l+pl)E Pm(O) Xb(m2 ) V1 Var(Sm)pm((0)) 2 1-Pm Pp1 52 2 meG P2 meE J if node 1 e G, (6.10a) s 1-pl (Pm(O) mb() 1-p ( Var(Sm)p(0) ar(S EW -=2 I E + 22 Var(Sm),m(0) 12-p m (_1-pl)p2 2 2s(1-pl) meG p2 meE klb(2) + 1if node 1 E E, (6.10b) 2(1-pl) 6.1 Computational Experience: We report on some of the numerical results for the mixed service discipline system. We use the data presented in examples 2 through 4, from section 3. Example 1 revisited. The same data applies as before except that nodes 1, 3, and 4 use the gated service discipline, while nodes 2 and 5 use the exhaustive service discipline. The mean waiting times are W1 = 139.5932, W2 = 76.1434, W3 = 147.2045, W4 = 152.6066, Ws = 111.6857. Example 2 revisited. Nodes 1, 3, and 6 use the gated service discipline, while the other nodes use the exhaustive service discipline. The mean waiting times are W1 = 340.1163, W2= 216.5708, W3 = 358.8866, W4 = 247.5317, W5 = 290.8190, W6 = 327.9425, W7 = 275.1886. Example 3 revisited. Nodes 1, 2, 3, 6, and 9 use the gated service discipline, while the other nodes use the exhaustive service discipline. The mean waiting times are W1 = 59.3568, W2 = 46.6172, W3 = 46.6183, W4 = 39.7035, W5 = 39.6888, W6 = 45.7251, W7 = 40.5487, W8 = 40.5423, W9 = 44.8510, W0l = 41.4468. 7. Summary and Conclusions We have presented a powerful technique, which we term the Individual Station (IS) technique for obtaining the mean waiting time at one or more nodes in a continuous-time polling system which uses either the pure exhaustive service discipline, the pure gated service discipline, or a system with a mix of both disciplines. The primary advantage of the IS technique is that one can now obtain the mean waiting time for a select subset of the nodes in the polling system. The techniques developed in the past require the simultaneous computation of the mean waiting times at all nodes. We developed two efficient algorithms, one based on an iterative approach (Algorithm 31

1), and the other based on solving a small system of equations of size 6 or less (Algorithm 2), to obtain the mean waiting times. Algorithm 2 requires O(M2) operations to find the mean waiting time at a node. In this paper, Algorithm 2 was developed only for the pure exhaustive, and the pure gated service systems. We presented a number of compuational results for the pure exhaustive service polling system, the pure gated service polling system, and the mixed service polling system. The computational results suggest that the iterative algorithm requires roughly the same number of iterations for a given value of p, for any number of nodes. This also appears to suggest that the iterative algorithm would probably perform better than an O(M3) algorithm in computing the mean waiting times at all nodes, so long as p is less than 0.90. Algorithm 2 does not depend on the value of p for its execution time. This algorithm is exact for the exhaustive service (gated service) system with 7 (6) nodes. However, it performs remarkably accurately in computing the mean waiting times. We experimented with p values very close to 1, and found that the maximum percentage error from the exact mean waiting time was less than 0.0025%. (For all practical purposes, this would characterize the algorithm as an exact algorithm.) These errors could be reduced even further, if needed, with little increase in computational effort, simply by executing the algorithm at a higher level (i.e., for the exhaustive service system, this would involve solving a system of 7 equations, instead of a system of six equations). The IS technique would be especially useful when evaluating several alternative system configurations in the preliminary design stages of, say, a computer network, where the analyst is usually interested in obtaining very accurate measures of performance only for a select few nodes in the network. There are a number of topics which remain to be explored, on applications of the IS technique. It would be of interest to consider applying the IS technique to other polling systems such as, for example, the nondeterministic polling system (Srinivasan 199 la, Boxma and Weststrate 1989). We have recently developed an iterative algorithm (Srinivasan 199 ib) to determine waiting time variances for the polling systems considered in this paper.1 1 The author thanks K.G. Murty and R. Saigal for their helpful comments. 32

References Y. Aminetzah. 1975. An exact approach to the polling system. Ph.D. dissertation, McGill University, Montreal, P.Q., Canada. J.E. Baker and I. Rubin. 1987. Polling with a general-service order table. IEEE Trans. on Communications, v 35, pp. 283-288. O.J. Boxma and J.A. Weststrate. 1989. Waiting times in polling systems with Markovian server routing. Proc. Conference Braunschweig, Springer Verlag. R.T. Carsten, E.E. Newhall, and M.J.M. Posner. 1977. A simplified analysis of scan times in an asymmetrical Newhall loop with exhaustive service. IEEE Transactions on Communications, v 25, pp. 951-957. R.B. Cooper. 1970. Queues served in cyclic order: waiting times. Bell System Tech. Jnl., v 49, pp. 399-413. R.B. Cooper and G. Murray. 1969. Queues served in cyclic order. Bell System Tech. Jnl., v 48, pp. 675-689. R.B. Cooper. 1981. Introduction to queueing theory, Second Edition, Elsevier Science Publishers, New York. M. Eisenberg. 1972. Queues with periodic service and changeover time. Opns. Research, v 20, pp. 440-451. M.J. Ferguson and Y.J. Aminetzah. 1985. Exact results for nonsymmetric token ring systems. IEEE Transactions on Communications, v 33, pp. 223-231. M., Fiedler, and V. Ptak. 1962. On matrices with nonpositive off diagonal elements and positive principal minors. Czechoslovakian Mathematical Journal, v 12, 1962, pp. 382-400. O. Hashida. 1972. Analysis of multiqueue. Rev. Elect. Commns. Lab., v 20, pp. 189-199. P. Humblet. 1978. Source coding for communciation concentrators. ESL-R-798, Electron Syst. Lab., Massachussets Institute of Technology, Cambridge, MA. L. Kleinrock and H. Levy. 1988. The analysis of random polling systems. Opns. Res., v 36, pp. 716-732. A.G. Konheim and B. Meister. 1974. Waiting times and lines in systems with polling. Jnl. Association for Computing Machinery, v 21, pp. 470-490. H. Levy. 1989. Delay computation and dynamic behavior of non-symmetric polling systems. Performance Evaluation, v 10, pp. 35-51. 33

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. I. Rubin and L.F. DeMoraes. 1983. Message delay-analysis for polling and token multiple-access schemes for local communication networks. IEEE Transactions on Selected Areas in Communications, v 1, pp. 935-947. G.B. Swartz. 1980. Polling in a loop system. J. Assoc. for Comp. Machinery, v 27, pp. 42-59. D. Sarkar and W.I. Zangwill. 1989. Expected waiting time for nonsymmetric cyclic queueing systems - exact results and applications. Management Science, v 35, pp. 1463-1474. M.M. Srinivasan. 1991a. Non-deterministic polling systems. Management Science, v 37, no 6, pp. 667-691. M.M. Srinivasan. 1991b. Waiting time variances in cyclic service systems. Technical Report 91 -37, Department of Industrial and Operations Engineering, The University of Michigan, Ann Arbor, MI 48109-2117, December 1991. H. Takagi. 1986. Analysis of polling systems, MIT Press, Cambridge, MA. H. Takagi. 1990. Queueing analysis of polling models: an update. Stochastic Analysis of Computer and Commn. Systems, Elsevier Science Publishers B.V. (North-Holland), Amsterdam. 34