CLOSURE METHODS FOR T H E SINGLE-SERVER RETRIAL QUEUE by D A V I D W A R R E N S A B O B. Sc. (Honours), The University of Alberta, 1972 Ph. D., The University of British Columbia, 1977 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Mathematics) Wc accept this thesis as conforming to the required standard: THE UNIVERSITY OF BRITISH COLUMBIA June, 1987 ©David Warren Sabo, 1987 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the requirements f o r an advanced degree a t the U n i v e r s i t y of B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head o f my department o r by h i s or her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department of The U n i v e r s i t y o f B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date w v j ^ e DE-6 (3/81) ABSTRACT This work focusses on the development and evaluation of so-called "closure methods" for solving the equations governing the time-dependent behaviour of single-server retrial queues. These methods involve assuming that particular known algebraic relationships between various characteristics of the correspondin steady-state queue also apply approximately when the queue is not at steady-state. The objective is to replace a problem requiring the solution of dozens or hundreds of simultaneous linear differential equations with a system of a few differential equations that has a solution that approximates those queue characteristics of immediate interest. The viability of such closure methods is assessed by examining the results of a series of test calculations. The methods described in this thesis apply to a retrial queue in which inter-arrival times for new customers, inter-retrial times, and service times are all assumed to be exponentially distributed. The steady-state solution for such a queue is described in some detail. A survey of the literature indicates that the description of this steady-state retrial queue has become quite sophisticated, whereas only very tentative steps have been taken in the study of the time-dependent behaviour of such queues. On the other hand, the time-dependent behaviour of the simple M/M/s queues have been studied to a much greater extent. The apparent value of closure methods in computing approximations to various basic time-dependent M/M/s queue characteristics motivated this examination of the extension of such methods to the single-server retrial queue. i i i . After discussing the basic approach to be used in devising and testing prospective closure methods for the single-server retrial queue, a variety of such methods is presented, with each being tested in considerable detail. It is found that three of the methods devised give results of comparable or better accuracy than those closure methods for the simple M / M / s queues which motivated this study. A l l recommended closure methods developed here involve systems of either two or three differential equations and permit the calculation of good approximations to four of the characteristics of greatest interest for non-stationary queues: the probability that the server is idle, the mean queue length, the variance of the queue length, and the conditional mean number of customers in the system given that the server is idle. Each of the methods presented is tested for queues with constant mean arrival, retrial and service rates, as well as for queues in which arrival and retrial rates vary sinusoidally with time. i v . TABLE OF CONTENTS Abstract ii List of Tables vi List of Figures vii Acknowledgements x CHAPTER 1 Review of the Retrial Queue, Time Dependence and the Notion of Closure Methods 1 1.1 The Single Server Retrial Queue 2 Description of the Retrial Queue 2 Equations Governing Queue Behaviour 4 The Steady State Solution 10 Evaluation of Steady State Moments 13 Related Queueing Models 15 1.2 Retrial Queues to Date 18 Cohen 18 Keilson, Cozzolini and Young 19 Alexsandrov 21 Falin 21 Hashida and Kodaira 22 Choo and Conolly 22 Kulkarni 24 Conclusions 25 1.3 Time-Dependence of Queues 27 Why Bother With Explicit Time Dependence? 27 Koopman's Work 28 Solution by Power Series Expansion 29 Determination of Time-Dependent Probability Generating Functions 32 Transform Methods 34 Previous Work on Time-Dependent Queues 35 Numerical Approaches 37 Truncation of Systems of Differential Equations . . . . 41 1.4 Closure Approximations for M/M/s Queues 46 Details of the Rothkopf and Oren Closure Method . . . . 46 Numerical Results 49 Period-Mean Values and Period-Mean Errors 54 Extensions and Other Approaches 57 V . CHAPTER 2 Closure Techniques for the Time-Dependent Single-Server Retrial Queue 66 2.1 The Basic Equations 67 Retrial Queue Characteristics of Immediate Interest . . 67 Derivation of the Basic Equations 68 Additional Criteria 70 Notation 75 2.2 One-Equation Closure Methods 76 Method IA ( n0(t) ) 76 Method IB ( m(t) ) 80 2.3 Two-Equation Closure Methods: II and m 87 Methods 2A and 2B ( II0(t) and m(t) ) 87 Numerical Performance of Methods 2A and 2B 89 Method 2C (n 0 (t) and m(t) ) 90 Constant Parameter Calculations Using Method 2C . . . 91 Oscillatory Parameter Calculations Using Method 2C .. 93 Computation of Period-Averages of ^ (t) and m(t) . . 100 2.4 Two-Equation Closure Methods: II , m and v 114 Method 2D ( I I 0 , m, v) 114 Numerical Performance of Method 2D 116 Method 2E ( I I 0 , m, v) 118 Constant Parameter Calculations Using Method 2E . . .119 Oscillatory Parameter Calculations Using Method 2E . . .122 Computation of Period-Averages of m(t) and v(t) . . . .123 Estimation of m 0(t) Using Method 2E 127 2.5 Closure Methods With Three Or More Equations 145 Three Equation Methods 145 Closure Methods Involving Four or More Equations . . . .149 2.6 Summary and Conclusions 153 Appendix A Details of Calculation of Exact Queue Properties. . . .157 Bibliography 159 v i . LIST OF TABLES TABLE 1 Period-Mean Values of P,,(t), m(t), and v(t) and Their Relative Errors for the M/M/l Queue When the Rothkopf and Oren Closure Method is Used . . . . 59 TABLE 2 Formulas for Moment Derivatives for the Single Server Retrial Queue 71 TABLE 3a Exact Period-Mean Values of n0(t) and m(t) and Mean Absolute and Relative Mean Absolute Errors When These Quantities are Computed Using Closure Method 2C for Various Values of 9Q 102 TABLE 3b Exact Period-Mean Values of H0(t) and m(t) and Mean Absolute and Relative Mean Absolute Errors When These Quantities are Computed Using Closure Method 2C for Various Values of c or v 103 0 0 TABLE 4a Period Averages of n0(t) and m(t) for Various New Customer Arrival Rates, Based on Method 2C 104 TABLE 4b Period Averages of n0(t) and m(t) for Various Retrial Rates for the Single-Server Retrial Queue, Based on Method 2C 105 TABLE 5a Period Averages and Mean Absolute Errors in n0(t), m(t), and v(t) for Various New Customer Arrival Rates Using Method 2E 129 TABLE 5b Period Averages and Mean Absolute Errors in II0(t), m(t), and v(t) for Various Retrial Rates Using Method 2E .131 TABLE 6 Period Averages and Relative Mean Absolute Errors in n0(t), m(t) and v(t) Using Method 2E With Various Amounts of Mixing, n , in the Closure Formula for m0(t) 133 TABLE 7 Summary of Serviceable Closure Methods for Computing the Time-Dependent Behaviour of the Single-Server Retrial Queue 155 v i i . LIST OF FIGURES FIGURE 1.1 Schematic Representation of a Single-Server Retrial Queue 3 FIGURE 1.2 Possible State Transitions and Differential Transition Rates for the Single-Server Retrial Queue 8 FIGURE 1.3 The Runge-Kutta-Felberg-45 Algorithm 39 FIGURE 1.4 Relative Error in Po(t) Computed Using Rothkopf and Oren's Closure Method for the Simple M/M/l Queue With Constant Arrival and Service Rates 60 FIGURE 1.5 Relative Error in m(t) Computed Using Rothkopf and Oren's Closure Method for the Simple M/M/l Queue With Constant Arrival and Service Rates 61 FIGURE 1.6 Relative Error in v(t) Computed Using Rothkopf and Oren's Closure Method for the Simple M/M/l Queue With Constant Arrival and Service Rates 62 FIGURE 1.7 Comparison of Exact p~(t) and pQ(t) Computed Using Rothkopf and Oren*s Closure Method for the M/M/l Queue in the Oscillatory Parameter Case . . . . . 63 FIGURE 1.8 Comparison of Exact m(t) and m(t) Computed Using Rothkopf and Oren's Closure Method for the M/M/l Queue in the Oscillatory Parameter Case 64 FIGURE 1.9 Comparison of Exact v(t) and v(t) Computed Using Rothkopf and Oren's Closure Method for the M/M/l Queue in the Oscillatory Parameter Case 65 FIGURE 2.1 n0(t) Computed Using Method IA in the Constant Parameter Case 83 FIGURE 2.2 Comparison of n0(t) Computed Approximately Using Method IA With Its Exact Values in the Oscillatory Parameter Case 84 FIGURE 2.3 Mean Number of Customer, m(t), in the Retrial Queue Computed Using Method IB in the Constant Parameter Case 85 v i i i . FIGURE 2.4 Mean Number of Customers, m(t). in the Retrial Queue Computed Using Method IB in the Oscillatory Parameter Case 86 FIGURE 2.5 Relative Errors in II0(t) and m(t) When Computed Using a Two-Equation Closure Approximation with Basic Methods 'A' and 'B' in the Constant Parameter Case . . 1 0 6 FIGURE 2.6 Relative Errors in n0(t) and m(t) When Computed Using a Two-Equation Closure Approximation with Basic Methods 'A' and 'B' in the Oscillatory Parameter Case 107 FIGURE 2.7 Relative Errors in no(t) Using the Two-Equation Closure Method 'C in the Constant Parameter Case With Various Arrival and Retrial Rates 108 FIGURE 2.8 Relative Errors in m(t) Using the Two-Equation Closure Method 'C in the Constant Parameter Case With Various Arrival and Retrial Rates 109 FIGURE 2.9 II (t) Computed Using the Two-Equation Closure Method 'C in the Oscillatory Parameter Case for Various Values of P 0 110 FIGURE 2.10 m(t) Computed Using the Two-Equation Closure Method 'C in the Oscillatory Parameter Case for Various Values of P Q I l l FIGURE 2.11 n0(t) Computed Using the Two-Equation Closure Method 'C in the Oscillatory Parameter Case for Various Values of c 112 0 FIGURE 2.12 m(t) Computed Using the Two-Equation Closure Method 'C in the Oscillatory Parameter Case for Various Values of c 113 FIGURE 2.13 Relative Errors in II0(t), m(t) and v(t) Using Closure Method 2D for Various Values of p in the Constant Parameter Retrial Queue . 1 3 7 FIGURE 2.14 Relative Errors in Jl0 (t), m(t) and v(t) Using Closure Method 2D for Various Values of cQ in the Constant Parameter Retrial Queue 138 FIGURE 2.15 Comparison of Exact and Approximate Quantities for the Retrial Queue Using Method 2D 139 i x . FIGURE 2.16 Comparison of Exact and Closure Method 2D Approximation Values of n0(t) and m(t) for Oscillatory Parameter Retrial Queues 140 FIGURE 2.17 Relative Errors in v(t) Computed Using the Closure Method 2E in the Constant Parameter Case 141 FIGURE 2.18 Performance of the Closure Method 2E When the Retrial Rate is Much Smaller Than the New Customer Arrival Rate 142 FIGURE 2.19 Comparison of Exact and Method 2E Approximation v(t) for Various Values of n in Equation (2.4.9) in the Oscillatory Parameter Case 143 FIGURE 2.20 Comparison of Exact m0(t) with m0(t) Given by Equation (2.4.9) for Selected Values of n in an Oscillatory Parameter Retrial Queue 144 FIGURE 2.21 Relative Errors and Comparison of Exact Quantities With the Results of the Closure Method 3A for a Constant Parameter Retrial Queue 150 FIGURE 2.22 Comparison of Exact Quantities With the Results of the Closure Method 3A for an Oscillatory Parameter Retrial Queue .151 FIGURE 2.23 Comparison of Exact Quantities With the Results of the Closure Method 3B for an Oscillatory Parameter Retrial Queue 152 X . A C K N O W L E D G E M E N T S I would like to take this opportunity to express my gratitude to Dr. P. E . Greenwood for her guidance and many helpful suggestions regarding the work which lead to the completion of this thesis. I would also like to acknowledge the assistance of Dr. G. Hooghiemstra, who introduced me to the retrial queue. Dr. R. A. Restrepo and Dr. M. Sion have had significant influence on the progress of my studies for this degree from the beginning. My studies for this degree have involved nearly five years as a part-time student in the program, and one year before that as an unclassified student. I would like to thank my full-time employer, the British Columbia Institute of Technology, for allowing me to rearrange my work schedule to permit me to pursue these studies, and for providing financial support for some aspects of the program. Most of all, I wish to acknowledge the tolerance displayed by my wife, Marlene, and by my children, Marlis and Robert, of a part-time program which all too often dominated evenings, weekends, and summers that could have been used very well in activities of greater interest to them. 1. CHAPTER 1 Review of the Retrial Queue, Time Dependence, and the Notion of Closure Methods 'Consider your verdict,' the King said to the jury. 'Not yet, not yet!', the Rabbit hastily interrupted. 'There's a great deal to come before that!' (Lewis Carroll, Alice in Wonderland) 2 . 1.1 The Single Server Retrial Queue / have seen something under the sun: The race is not to the swift or the battle to the strong, nor does food come to the wise or wealth to the brilliant or favour to the learned; but time and chance happen to them all. (Eccleaiastes 9:11. NIV) Description of the Retrial Queue This thesis deals with an approach to computing the time dependent behaviour of busy telephones and similar systems. The main results apply to the simplest such system, denoted technically as the M / M / l / 1 retrial queue. The situation to be dealt with can be described schematically as in Figure 1.1. Here, the box labelled S represents the 'server', which at any given time can be idle ( £ = 0) or busy (jr = 1). If the server is idle when a customer approaches, the customer receives immediate service. On the other hand, if the customer finds the server busy, he must abandon his attempt to engage the server, and goes into 'orbit'. Blocked (or orbiting or milling) customers are presumed to continue to approach the server repeatedly until they find the server idle, receive service, and exit the system, though it is possible to formulate more sophisticated models in which an orbiting customer eventually tires and abandons further attempts at service [Cohen, 1957]. The characteristics of this system which we will be interested in calculating are quantities like the number or mean number of customers in the system and its variance as a function of time, the fraction of the time that the server is idle, the length of time a customer can expect spend in the system to receive service, and so on. 3 . . , Server 1 ^ j ~ , > • — New customers arrive at an average rate of A per unit time. Blocked customers (arriving when the server is engaged) must go into 'orbit' or 'mill'. Each will attempt to engage the server an average of v times per unit time until successful in receiving service. Served customers leave the system. The server is assumed to be able to accommodate an average of u customers per unit time. F I G U R E 1.1 Schematic Representation of a Single-Server Retrial Queue. 4 . What makes the consideration of this problem nontrivial is that in actual systems of this type, the exact • times at which new customers enter the system • times at which orbiting customers will 'retry' for service • length of time the server will spend with each customer are random variables. This means that rather than being able to state precisely when a particular customer (the first or the tenth, etc.) will first approach the server, or precisely when the server will be approached by a new or orbiting customer, or precisely how long the server will spend with a particular customer, we will only be able to speak of the probability of such events occurring in time intervals of some length or of the expectation of these quantities. The degree of complexity this indeterminacy brings to the analysis of the problem depends on the nature of the probability distributions of the random variables involved. Equations Governing Queue Behaviour We assume that the arrival times of new and retrying customers are statistically independent and that the probability of more than one arrival of each type occurring in some time interval approaches zero as the length of the time interval approaches zero. While it is possible to envision unusual circumstances under which, say, telephone calls placed to a particular receiver (server) do not 5 . conform to these assumptions1, they do reflect a sufficiently wide range of experience to justify adopting. These assumptions are consistent with adopting the exponential probability distribution for interarrival times for new as well as retrying customers. A similar, though weaker, case can be made for service times having the exponential distribution.2 An immediate consequence of assuming that arrival, retrial, and service rates are exponentially distributed is that the state of the system is completely specified by the pair of numbers (£ ,n), with E, = 0 or 1 denoting the server as idle or busy respectively, and n indicating the number of customers in orbit.3 Let A be the mean number of new customers entering the system per unit time, v be the mean number of retrials a customer makes per unit time and u be the mean number of customers the server can handle per unit time. For the xThe telephone calls of a Wall Street stock broker on the afternoon of October 24, 1929 comes to mind as a potential example of a of a retrial queue in which both assumptions are violated. 2In his first paper on the subject, Cohen (1957) refers to experimental results which 'fully justify' the adoption of the exponential distribution for interarrival times as well as service times. Because the exponential distribution does reflect the behaviour of a considerable variety of systems which can be modelled as queues, and because of its mathematical simplicity, it is usual to adopt this distribution for interarrival times and service times as a first approximation at least. Since the initial goal here is to develop computational techniques for a model already appearing in the literature, nothing further will be said in attempting to justify the model (though the issue is not entirely irrelevant). 3No other continuous probability distributions share the 'lack of memory' or Markov property of the exponential distribution (see, for example, Ross (1983), pages 24 -25). If interarrival times or service times were assumed to follow some other such distribution, then in addition to the values of £ and n, a complete description of the state of the system might also require the time at which the most recent customer arrived or the time at which the currently served customer entered service. 6 . moment, we can think of X , u, and v, as being constants, though the end result of the derivation below does not change materially if they are sufficiently smooth functions of time. Using the properties of the negative exponential distribution, we can state the following for some small time interval At: Pr(a new customer arrives in the interval At) = A At + o(At) , (la) Pr(service of a customer is completed in the interval At| the system is not empty) = u At + o(At), (lb) Pr(a customer retries for service in the interval At| there are n customers in orbit) = vnAt + o(At), (ic) Pr(more than one new customer arrives in At) = o(At) , (Id) Pr(more than one customer completes service in At) = o(At), (le) Pr(more than one customer retries for service in A t) = o (At) , (If) Pr(more than one event of any type in At) = o(At) . (lg) The notation o( A t) indicates a quantity that becomes negligible compared to At as At approaches zero (in an abuse of notation, we can say that o(A t) is any 'quantity' satisfying lim o( At)/ At = 0 . ) At-*0 It is now a simple matter to derive a system of differential equations for the probabilities {p (t), E, = 0 , 1 ; n = 0 , 1 , 2 , . . . } where p (t) = Preserver is idle and there are n customers in orbit at time t), On n = 0, 1, 2, ... and p (t) = Pr(server is busy and there arc n customers in orbit at time t), I n n = 0, 1, 2, ... To do this, we simply use the so-called law of total probability, which can be written here as oo Jl p r ( t + At)= z z p. f t ; n . r r A t ; , 5 = o,i n = 0,1,2,... where 11^ ^^(At) is the conditional probability that the system will be in state ( 5, n) at time t + At given that it was in state (CJ , k) at time t.4 If we are concerned only with retaining terms to first order in At in equations (2), then the conditions given in equations (1) drastically reduce the number of nonzero terms in the summation above. The solid arrows in Figure 1.2 represent the n (AtJ which arc first order in A t or larger as At becomes small. Using that diagram, for example, we have for any n = 0, 1, 2, 4Here, I I ^ f ^n ('ht) depends only on At because all the interarrival and service times are assumed exponentially distributed. If that were not true, the II , _ would also depend on the value of t. Uik, E,n ( l - u - A ) d t ( l - u - A ) d t ( l - u - A ) d t F I G U R E 1.2 Possible State Transitions and Differential Transition Rates for the Single-Server Retrial Queue The six pairs ( £,n) above represent possible states of the retrial queue, with £ = 0 or 1 corresponding to the server being idle or busy respectively, and 'n' being the number of customers in orbit. The arrows represent possible system state transitions for this single-server retrial queue with the corresponding differential transition rates as indicated. P„ (t + At) = Pn f t ; - n „ C A t ; + p f t ; - n . C A t ; On On On,On In ln,On P0n(t)* Pr(the system is in state (0,n) at time t and no arrivals of either type occur in the interval At) + p (t) • Pr(the system is in state (l,n) at time t and there In is one service completion in A t but no arrivals of either type) p ( t ) { l - AAt + o ( ( \ t ) } { l - nvAt + o(At)} On + P l n ( t ) { ] i k t + o( A t ) } { l - \ k t + o ( A t ) } { l - n v A t + o ( A t ) } ? 0 n ( t ) ~ X p 0 n ( t ) ^ ~ n V p Q n ( t ) h t + M p l n ( t ) h t + ofAt; . Subtracting PQn(t) from either side, and dividing this last equation through by At gives At ~ ' v / r O n ' ' H ^ 2 n ' At In the limit as At -y 0, this reduces to • p n (t) = -(X + nv)pn (t) + up f t ; , dt *On On In 10 . In a similar fashion we find that (t) + Xp„ (t) - — — p ( t ) = -(X + \l)p (t) + XP ~rQ dt In In 1 , n-1 un (3b) + (n+l)Vp„ ( t ) , n=0,l,2,... , 0 ,n + l with the understanding that P^n(t) is identically zero whenever n < 0. The characteristics of the system which are of practical interest are all calculable from the p_ (t) ( £ = 0,1; n = 0,1,2, ...) in a straightforward manner. In However, these p (t) can be obtained only by solving the combined infinite sets of differential equations in (3a) and (3b). The Steady State Solution A n exact solution of the combined system of equations above has been obtained only when X, p, and v are constants, A < y , and in the limit that t becomes infinite. In this case, the p r (t) become independent of the value of t, and the values of the constants p = p_ ( ° ° ) are unaffected by the initial state of the system. This so-called steady-state solution can be obtained in one of several ways [Gross & Harris, 1974, pages 44-51]5. While we will be primarily interested in 5 That such a steady-state solution of equations (3a) and (3b) exists is inferred from the construction of it which follows. This is a common approach (see, for example, Gross and Harris (1974), pages 38 - 40 and Cooper (1981), section 2.3). It is recognized that the existence of such steady state solutions involves the concept of ergodicity, for which some rather general results are examining techniques for approximating the time dependent behaviour of the retrial queue, the major technique studied makes extensive use of this steady state solution. Our intent is to use the characteristics of the calculable steady state solution to develop approximation techniques which will allow bypassing the full systems of equations in (3a) and (3b). Perhaps the simplest approach to obtaining the steady state solution to this single server retrial queue problem is to assume the solution exists, replace the derivatives in (3a) and (3b) by zeros, and solve the resultant difference equations recursively. With the derivatives replaced by zeros, the systems (3a) and (3b) can be rearranged to give p = ~^~(c + n)p , (4a) In c On c (Q + 1) c . ,„,, - P , „ " „ , , f P , + , P n J > (4b> 0,n + l (n + l)p "In n + 1 * 1,n-1 "On n = 0 , 1 , 2 , where A , A , - . p = and c = . (5) y v known (see Shiryayev (1984), section 1.12 and Chapter V). To quote Gross and Harris (1974, page 40), "It is generally not necessary to appeal to these [ergodic] theorems to determine the conditions under which queueing processes are ergodic, since they usually fall out from other considerations." We can apply the remarks of Cooper (1981, page 19 and following) for the M/M/s queue, to indicate here that if limiting values of Porft) and Pj„(t) exist as t->-°° , then equations (3a) and (3b) imply that their derivatives p'0 (t) and Pj n(0 must also have limiting values as t ->- 0 0 . But the existence of limiting values for the probabilities is possible only if the limiting values of their derivatives are zeros. If on replacing the derivatives in equations (3a) and (3b) by zeros, we obtain a system of difference equations with no solution or with no acceptable solution (as when A > y , when the probabilities diverge with increasing n to »), we will presumably have encountered an example of "other conditions" indicating that the queue is not ergodic. 12 . Alternating between (4a) and (4b) for increasing values of n, we get successively P10 = ppoo poi = c p p o o p n = p 2 ( c + 1 ) p o o 2 P c(c + 1) P02 2 P00 p3(c + 1)(c + 2) P12 2 P00 p3c(c + 1) (c + 2-J P03 2 3 P00 and so on. A reasonably clear pattern in the formulas for these steady state probabilities has already emerged, leading us to hypothesize that for general n, n [c + n - l\ - P c(c+l)(c+2) •••(c+n-1) _ [ n Pan ~ n~> P00 ~ \ n / P P00' ( 6 a ) n = 1 , 2 , 3, and p = Pn + 1(c + l ) ( c + 2 ) - - - ( c + n) = (c + n\ n+l pIn n! F00 \ n J H y00 n = 0, 1,2, ... . Substitution of (6) into (4a) and (4b) confirm that the formulas (6) do indeed satisfy the system of difference equations. Finally, since the sum of all of the steady state probabilities must be 1, a bit of routine algebra allows determination 13 . of p0Q as P Q 0 = (1 - P)C + 1 (?) which together with the formulas in (6) completes the specification of the steady state probabilities. We note in passing that p , the probability that at steady state there are a n total of n customers in the system, either in service or in orbit, is given by ( c + n \ n Pn = P0n + P T - = I ~ ! " J P " ' 1 " p ; C + J ( 8 ) These probabilities form the so-called negative binomial distribution with parameters c+1 and 1 - p. Evaluation of Steady-State Moments In what follows, we shall make considerable use of formulas for steady state moments for the system. The evaluation of such moments is greatly facilitated by first,obtaining corresponding probability generating functions oo , , , c + l P (z) = Z p z n = ( 1 ~ P ) (9a) 0 n=0 ° n (1 - Oz)C P l ( Z ) - z P l y = -n=0 (1 - pz) 14. and c+1 P(z) = E p z n = P ( Z ) + Z P ( Z ) = ( 1 ~ P ; n = 0 " 0 1 ( 1 _ p z ; < ^ C9c; The steady state quantities of particular importance are = P r C s e r v e r idle) = E = P (1) = 1 - p , n = 0 f 20; II = P r e s e r v e r busy; = E p = 1 - 1 1 = p, 1 „ _ n i n 0 n = 0 (11) mQ = E n n = 0 n dz 0 P„(z) c p , z = l (12) E (n + l)p = n = 0 — P.(z) dz 1 z = l P(cp + 1) 1 - P (13) m = mean number of customers in the system mo + mi p f c + 1 ) 1 - P (14) y 2 L n p n = 0 On dz 2 0 z = l p c f p c + 1) + m = — r i 5 ; 2 - p ( Z P , r z ; } dz 2 J + m. z = l p { f c p + 1) (cp + 2) - (1 - p) } (1 - P)2 (16) 15 . C = expectation of the square of the number of customers in the system 2 2 (c + 1)(cp + p + p) 2 (1 - P) (17) T. n p n = 0 On dz z = l + 3£, + 2m 0 0 pc {(pc + 1)(pc + 2) - p(l - p)} (18) 2 a - P; We also note that the variance, v, of the number of customers in the system is given by v = £ - m P(c + 1) (1 - P) (19) Related Oueueing Models Before leaving this rather condensed description of the basic system we will be examining in greater detail shortly, it is worth recognizing that the problem is in a sense intermediate between two of the very small number of queueing models for which not only the steady state problem is easily solved, but analytic solutions have been found for all positive values of t. 16 . If we regard the orbiting customers as forming a sort of balky waiting line, the single server retrial queue is seen to be a generalization of the well-known M/M/ l /oo single server queue with service in random order (SIRO).6 The big difference between the two is that whereas in the ordinary M/M/l/°° queue, the server will never be idle as long as there are any customers in the system, in the retrial queue, it is quite possible for the server to be idle for a time even though quite a large number of customers are orbiting unserved. We would expect that as the mean retrial rate, V , becomes larger and larger, the retrial queue will look more and more like the M / M / l / oo/SIRO queue. We note that the system of equations governing the time-dependent behaviour of this queue have been solved formally only in the case that the arrival and service rates are constants. Alternatively we can regard the orbiting customers as similar to customers who are not waiting in the system at all, giving the impression that the retrial queue is a bit like an ordinary single server queue with no waiting room at all (ie. the ordinary M / M / l / l queue). This queueing model is perhaps too simple to provide much useful insight into the characteristics or behaviour of the single server retrial queue. 6Queues are usually described with a notation developed by Kendall (1953) and others, making use of a collection of symbols separated by slashes: A/B/C/D/E. Here 'A' indicates the interarrival time distribution, 'B' indicates the service time distribution, 'C indicates the number of parallel service channels, 'D' indicates the capacity of the system including customers in service, and 'E' indicates the waiting line discipline. An extensive list of the symbols commonly used can be found in most basic books on queueing theory (for example, Gross and Harris (1974), page 9 and pages 464-471). The symbol 'M' is used to indicate the exponential distribution, so the notation M/M/ l /oo denotes a queue with exponentially distributed interarrival and service times, a single server, and infinite capacity. To indicate that customers in the system are served in random order, we could expand the notation used to read M / M / l / "VSIRO. 1 7 . The similarity between the single server retrial queue and the ordinary M / M / l / °°/SIRO queue is what prompted the investigations described in Chapter 2 below regarding the adaptation to the retrial queue of so-called closure methods which have been applied successfully to computing the approximate time dependent behaviour of the simpler M / M / l / «> queues. It should be noted also that the retrial queue is distinct from so-called 'feedback queues' in which customers may return to the queue to await service again after they have completed service. For example, while a particular customer might have been released from service, the service provided may have been incomplete or unsatisfactory in some way. In the retrial queue, customers return to the queue after an unsuccessful attempt to engage the server, but once served, they are viewed mathematically as leaving the system forever. 18, 1.2 Retrial Queues to Pate Miller's Law: 'You can't tell how deep a puddle is until you step into it.' (The Official Rules. Paul Dickson, ed.) Cohen The literature on the subject seems to agree that Cohen [1957] was the first to consider the retrial queue problem in detail. In this lengthy paper, he deals both with the underlying assumptions of the formalism as they apply to telephony as well as with developing the mathematical formalism in great detail. Cohen actually treats the multiple-server retrial queue in the paper. He uses an approach similar to that outlined in the section 1.1 of this thesis leading to equations (3a) and (3b) to derive the underlying system of differential equations which can be solved for the probabilities of finding the system in any of its possible states. . However, to cope with the considerably greater complexity of the resultant equations whose solutions are the steady state probabilities, he first transforms the equations for the probabilities into a system of a smaller number of equations for probability generating functions. Despite this stratagem, obtaining the formal solution of the problem is no mean feat. The expression Cohen gets for the mean number of 'congested calls' at steady state (what we would refer to here as the mean number of customers in orbit) is a fraction with both the numerator and denominator containing double index summations over products of generalized Laguerre polynomials and integrals 19 . involving components of probability generating functions for the system. Cohen then goes on to demonstrate how solutions for a variety of special subcases of the problem can be obtained from his general solution. His results are so detailed and comprehensive that no attempt will be made to itemize or summarize them here. Keilson. Cozzolini and Young The next appearance of the retrial queue in the literature seems to be in Keilson, Cozzolini and Young [1968]. These authors consider both the single-server and the two-server retrial queue with exponentially distributed interarrival times and general service time distributions. They consider the time-dependent behaviour of the queues but obtain formal and numerical results for both types of system only at steady state. Note that when service times have a general (non-exponential) distribution, the single-server retrial queue is characterized by something like a pair (n, £), where n is a non-negative integer giving the number of customers in orbit, E = 0 indicates the server is idle, and £ > 0 indicates that the server has been engaged by the current customer for £ time units already. Instead of the two discrete values of E (0 and 1 above) required to characterize this system when service times are exponentially distributed, E is now needed to be a non-negative real quantity. For the steady-state single server retrial queue, Keilson, Cozzolini and Young derive formulas for the probability that the server is idle (P = 1 - A T), the mean number of customers waiting in orbit (m = [ i- \ 2 T2 + ( — )AT]/[I - AT]) 2 U the mean waiting time in orbit of a customer ( "f = m/ A), and the mean number of attempts at service per customer for each successful call ( = 1 + v"r )• All of these quantities are determined by the mean service time, T, which is simply I/u when service times are exponentially distributed. For the multiple-server case, Keilson, Cozzolini and Young assume exponentially distributed service times as well as exponential arrival and retrial rates. This permits them to describe the system by a pair of numbers, (n, £ ), where n is the number of customers in orbit and g is the number of servers engaged. Thus, g can only assume values from a finite set of non-negative integers. They obtain formal results for the steady-state of a two-server retrial system using an approach borrowed from the study of current flows in networks of resistors in electrical engineering. Their analysis gives formulas for quantities P , Q , and R , which are the probabilities that there are m customers in m m m orbit and zero, one, or two servers busy respectively. The characteristics of the retrial queue mentioned in the previous paragraph for the single-server case can all be computed for this two-server case as sums of multiples of the P , Q , m m and R^ . While use of a computer is necessary because these sums cannot be rearranged into a simple closed algebraic form, the necessary formulas arc not difficult to program. 21 . Alexsandrov Chronologically, the next treatment of the retrial queue appears in Alexsandrov [1974], who derives a number of results for the M/G/ l system considered earlier by Keilson et. al. [1968]. This short paper concentrates on obtaining expressions for the first two moments of the system at steady state. He also derives expressions for the expected time of sojourn of a customer in the system, the effective arrival rate of retrying customers, and the expected number of retrials each customer will make before succeeding in obtaining service. Alexsandrov's paper provides a number of references to earlier papers on the subject appearing in the Russian literature. Falin The Russian torch, as far as the retrial queue is concerned, is then taken up by G. I. Falin, culminating in a lengthy paper (Falin, [1979]) dealing with, in his words, the 'non-stationary regime of operation of the system', and concentrating on characterizing a 'system busy period'. Falin defines the busy period for the system as the time interval from when a customer enters service in an otherwise empty system until the system becomes totally free of customers either in service or in orbit. While for non-retrial queues, the system busy period defined in this way would correspond to the server busy period, for the retrial queue there may be many intervals during the system busy period when the server is idle. Expressions are obtained for the 22 . length of the system busy period and the number of service completions in such a busy period. It is noteworthy that this difficult paper makes reference only to the Russian literature. Hashida and Kodaira Meanwhile, Hashida and Kodaira [1975] published a treatment of the retrial queue in which new customer interarrival times were assumed exponentially distributed, but inter-retrial times and service times had general distributions. They were interested in the problem of the inefficient use of circuit switching networks because of line time wasted by customers during repeated unsuccessful attempts to obtain service. Their approach is similar to that of Keilson et. al. [1968], and they produce some numerical results for the probability that incoming calls will find the server busy. Chop and Conolly Choo and Conolly [1979] present some results for the system time of customers, service idle time, and system busy period for the steady state M/G/ l retrial queue. It is in this paper that the term 'orbit' is first used in reference to the situation in which those customers unsuccessful in obtaining service on their first attempt find themselves. Choo and Conolly attempt to extend Alexsandrov's work by obtaining formulas for the customer's system time (waiting time in orbit plus time in eventual service), the server idle time and the system busy period. Their derivations assume service times are exponentially distributed 2 3 . as are new customer arrival times and retrial times, though some of the results are valid when service times have a general distribution. Choo and Conolly derive an infinite system of equations for probability distribution functions of the first passage times between various states of the queue. This is motivated by the observation that the system busy period is the first passage time from the state ( £, n) = (1,0) to the state (£ , n) = (0, 0). Unfortunately, it does not appear possible to find an exact time-dependent solution to this system of equations. The authors present computational algor-ithms for calculating the first and higher moments for these first passage times, from which moments of the system busy period can be obtained. For example, they find that the mean system busy period is given by B = * [—1 - 1], s * Poo where p00 is the probability that the system is empty (see equation (7) in section 1.1 above). Choo and Connolly apparently have greater success in obtaining simple algebraic formulas for the mean system time of a customer (the time from when the customer first attempts to obtain service until service is completed for him) and its variance. Unfortunately, Kulkarni [1982] and Falin1 point out an apparent error in the approach of Choo and Conolly to the analysis of a customer's waiting xPrivate Communication with B. W. Conolly. In a short letter to the editor of the Journal of Applied Probability, appearing immediately after Kulkarni [1982], Conolly acknowledges the 'serious difference of opinion' expressed by Kulkarni [1982] as an error in the waiting time analysis appearing in Choo and Conolly [1979] and laments the apparent lack of collaboration between people in Britain, the United States, and the USSR, interested in the retrial problem . 24 . time in the system. Kulkarni presents a corrected analysis of the problem and finds it impossible to obtain a simple closed form expression for the waiting time. Instead, his result is in the form of a numerically well-behaved convergent infinite series whenever P < 1, with terms defined by a set of recursion relations. Choo [1978] considered the retrial queue in his Ph. D. thesis as well. Kulkarni The most recently appearing paper on the retrial queue seems to be Kulkarni [1983]. Kulkarni develops a simple formula, AR = ^SRS, where A is the mean arrival rate for new customers, A^ is the mean service rate, R is the average number of unsuccessful retrials made by a customer, and R^ is the average number of unsuccessful retrials made during one service period. This formula is reminiscent of the so-called Little's Formula (L =AW)2, though its intuitive 2The intuitive basis of Kulkarni's formula comes from recognizing that each unsuccessful retrial belongs to some customer and takes place during some service period (according to the model, an unsuccessful retrial cannot take place when the server is idle). Little's formula as stated, with L being the mean total number of customers in the system at steady state, including customers in service, A the mean arrival rate, and W the mean time a customer spends in the system including time in service, applies with great generality to a wide variety of steady state queues and is one of the most celebrated and discussed results in queueing theory. It too has a simple intuitive base (the mean number of customers in the system must be the mean number of customers that arrive during the time a single customer spends in the queue and in service). Kulkarni's proof of his formula is modelled on a proof developed by Stidham [1972, 1974] for Little's formula, and Kulkarni claims the same sort of generality for his formula relative to retrial systems as is the case for Little's formula for non-retrial queues. 25 . plausibility belies the intricacy of the analysis required in a rigorous proof. This is another feature shared by the two formulas. Kulkarni applies his formula to the steady state analysis of a retrial queue involving two types of customers having exponentially distributed arrival, retrial and service times, but with differing mean values. He obtains expressions for the expected waiting time for each type of customer in the system and the expected number of each type of customer in the system. One of the interesting special cases to which Kulkarni's results can be applied is a situation in which one customer decides to conduct retrials at a different average rate than all other customers in the system. Kulkarni presents a formula for the expected waiting time (before service) of such a customer. One would expect that a detailed understanding of the characteristics of such a result would have some very useful applications.8 Conclusions Though a number of the works cited in this section contain derivations of equations governing the time-dependent behaviour of single-server or multiple-server retrial queues (for example, Keilson, Cozzolini and Young [1968]), the resultant equations invariably do not have simple closed form solutions. Because 3This result of Kulkarni's would seem to be a useful starting point for developing strategies for a customer who finds himself in orbit, and would like to minimize his expected waiting time in orbit. It would provide some way of estimating the value of 'keen-ness' in such a customer. 2 6 . of this, major results are presented only for the much more tractable steady-state case for such queues. This twenty-five or more years of work on the steady-state retrial queue has led to some very sophisticated results. However, it is not difficult to think of actual systems corresponding to retrial queues to which these steady-state results do not apply because arrival, retrial, and/or service rates vary in time to the extent that a steady-state regime cannot be established. Because the main objective in this thesis is to examine some viable approaches to determining the time-dependent behaviour of a single-server retrial queue, it is unnecessary to discuss any further the results described briefly in this section for the steady-state retrial queue. However, we note that these steady-state results could bear a more detailed examination, since there appear to be some discrepancies outstanding between the forms of results obtained by different authors. One example is the relatively simple expression Keilson, Cozzolini and Young [1968] obtain for the mean waiting time in orbit of a customer, which contrasts markedly with, for example, the complexity of the corresponding result obtained by Kulkarni [1982]. 27 . 1 . 3 T ime-Dependence of Queues Anderson's Law: 'I have yet to see any problem, however complicated, which, when you look at it in the right way, did not become still more complicated.' (The Official Rules, Paul Dickson, ed.) Why Bother With Explicit Time Dependence? The main part of this thesis is concerned with formulating and testing a particular set of approximation methods for studying the time-dependence of the single server retrial queue. In this section, the motivation for these studies is developed by discussing the importance of considering time-dependent behaviour of non-stationary queues, by noting the inapplicability of various traditional approaches to solving the relevant differential equations or systems of differential equations, and finally, by noting the computational problems encountered using standard numerical approaches to the problem. For some years it has been recognized that there were many applications of qucueing theory for which the steady state solutions were inappropriate in principle if not in practice. The problem was not so much a need to take account of the transient time-varying regime as systems evolve from some initial state to a steady state, but rather that many systems have intrinsic non-stationarity because arrival rates and/or service rates vary with time. Further, this variation in time may not (in fact, likely will not) be as simple low order polynomials in t, but more likely as a periodic function involving trigonometric functions. 28. Koopman's Work Koopman [1972] pointed out that the queue formed by aircraft wishing to take-off or land at an airport had a definite 24-hour periodicity. Within each 24-hour period, demand for access to runways varied from zero during the early morning hours to a level during rush hours which could for some airports be well in excess of the allowable runway capacity. To model such a system as a steady state queue is clearly inappropriate. To deal with such a significant time variation in the characteristics of the queue, Koopman first noted that the smallest time interval over which a change in the state of the queue could occur was strictly positive. Since an aircraft takes of the order of one or two minutes to enter the runway and leave it, either during take-off or landing, changes in the state of the queue can only occur on that sort of time scale, and so the whole problem is discrete. Secondly, he noted that for even the largest airports, the number of aircraft that can be waiting on the ground to take-off, or the number which can be waiting in the vicinity to land is necessarily quite small. Thus one can realistically assume that there is a specific finite number (maybe a few dozen, and certainly not more than a few hundred), with the probability of the number of aircraft in the system ever exceeding that number being zero in practice. As a result, the infinite system of linear differential equations analogous to (1.1.3a) and (1.1.3b) above which needs to be solved to describe the time-dependent behaviour of the system can be truncated to a finite system. The solution of the equations using numerical methods then becomes feasible. Though Koopman was concerned primarily with the repercussions his computational results had for airport design and operation, and though he considered a number of other issues1 in the application of queueing models to systems like this, the main impact of his paper seems to have been to establish the necessity of taking time dependence of queues into account explicitly in modelling many real systems. His quantitative results showed that quantities like the expected queue length and its standard deviation, the probabilities of finding the queue empty or finding it full, the cumulative number of customers turned away from a full queue, and the expected waiting time for a customer entering the queue, varied appreciably with variations in arrival and service rates. Solution by Power Series Expansion There are several straightforward approaches to determining analytic solutions for systems of linear differential equations like (1.1.3a) and (1.1.3b). Power series 1Koopman also considered the issue of choosing an appropriate probability distribution for the service times. He carried out a number of calculations using first the exponential distribution (corresponding to minimum information regarding the time of the next service completion) and then assuming constant service times (corresponding to maximum information about the time of the next service completion). He found that the quantitative results were relatively insensitive to which of these service time distributions was used. Considering these alternative models to represent extremes, he concluded that it was an acceptable approximation to adopt the exponential distribution to take advantage of simplifications in the formalism, in place of a more realistic but likely much more complicated distribution. Koopman's devotion to realism is further demonstrated by his use of actual airport records to evaluate his computed results. 30 . solutions are sometimes useful. In this case, we can write 2 P ( t ) = Z a ( n ) t j , P ? (t) = Z b ( n ) t j , (1) j=0 3 l n j=0 J n = 0 , 1 , 2 , . . . Substituting these two expansions into equations (1.1.3a) and (1.1.3b) and requiring that the coefficients of a given power of t on either side of each equation be identical, we obtain , , V*(n) - (X + n v ) a ( n ) <n> = 1 1 (3a) j + 1 and X b ( n - 1 } + X a ( n ) + (n + l ) V a ( n + 1 ) - ( X + V ) t ( n ) b [ n ? = 2 2 2 1 ( 3 b ) 3 + 1 j + 1 When n = 0, the first term in the numerator of (3b) must be dropped. These formulas express coefficients with lower index 'j+1' completely in terms of 2Carter and Cooper [1972] employ a power series approach in a study of the stationary waiting time distributions for various queues. They factor a factorial out of the coefficients in the series. In our case, that approach would make formulas (1) look like (n) „ d ( n ) c p ( t ) = Z 2 tJ , p 7 (t) = Z —J- t ° . (2) on . . . , In . - . , 3=0 j ! 3=0 3 I The coefficients c j D and d j " are then given by equations (2a) and (2b) respectively, but with the denominators 'j+1' dropped. One suspects that this approach will lead to difficulties in calculating factorial-sized quantities explicitly when large values of j are handled. 31 . coefficients with lower indices 'j', and so form a feasible basis for computing these coefficients recursively. The above formulas apply when A, u , and v arc all constants, and the situation would certainly become more complicated if any or all of these quantities were time-dependent. One can readily imagine the additional complexity involved in separating powers of t if, for example, arrival rates were sinusoidal. Of course, the real problem with this approach is that even formally convergent power series expansions form hazardous bases for computation except when the value of the independent variable is very small. In this case, the numerical sum of these series will always be a number between zero and one (since the sum represents a probability), and yet, for values of t much greater than one, the small sum may be constructed out of terms having magnitudes much greater than the final sum, and so serious round-off error can occur in summing the series as a result of the finite precision with which all computing machinery stores numbers. For example, when A = 0.5, u = 0.7, and v = 0.3 and all constant, the probability PQO(20) is very near its steady state value of 0.035. Yet in evaluating this quantity using the power series expansions above, one finds that more than 40 60 twenty terms in the sum including those involving t through t are all greater in magnitude than 10^ 2 , and alternate signs in an irregular fashion. This means that something like the fifteen lowest order significant figures of the computed series sum will be garbage. Since the highest precision (in the computerese sense) normally used in scientific calculations is around 15 or 16 significant 32 . figures, the inappropriateness of this approach to computing these probabilities is clear.3 Determination of Time-Dependent Probability Generating Functions An alternative approach to solving equations (1.1.3a) and (1.1.3b) exactly is to transform each set of equations into a single partial differential equation for the corresponding probability generating function. All probabilities and moments of the system can be computed straightforwardly from the probability generating functions. For the retrial queue, we define the two generating functions, Pn(z,t) = Z p (t)z , P(z,t) = Z p 7 (t)z . (4) —0 „ On —1 „ ln n=0 n=0 Multiplying equations (1.1.3a) and (1.1.3b) each as written by z", summing the resultant equations over n, and rearranging terms appropriately, we obtain (5) ~ - Pi(Z/t) = -(\+]x)P1(z,t) + \zp1(z,t) + A F 0 r z , t ; + -^p_o(Z,t). 3The basis of this rough quantitative estimate of the effect of numerical roundoff error on the accuracy of the computed sum of these power series is somewhat simplistic. More sophisticated analyses of the problem are available (see, for example, Lyashenko and Nikulin [1986]). Such analyses appear to indicate that the true cumulative round-off error in computed sums can be significant even when there are not such large differences in the magnitudes of the individual summands. 33 . So, the doubly infinite system of ordinary differential equations has been replaced by a system of two partial differential equations. These equations can be solved directly (see Zachmanoglou [1976], chapter 10), but that approach seems to result in a system of two integral equations which are solvable only by iteration, and the complexity of the expressions arising out of the iterative approach may make attainment of adequate accuracy infeasible. 4 4 Following the method of characteristics as outlined by Zachmanoglou [1976], we obtain: •vt v Q ( z , t ) - V + X -ze • V ( t - t ) v 0 ( z o ( t ) , t ) ze •v(t-t) v i ( z o ( t ) , t } s d t ' ttft-v (z,t) = where 2X - + u v 0 < z 1 ( t ) ,t) j j - j - - (X + y - AzTj v 1 ( z i ( t ) ,tJ^d~i Zq (t) = ze V ( t t } , z (t) = z are the equations of the characteristics of equations (5) passing through the point 34 . Transform Methods Alternatively, one could attempt to solve equations (5) using Laplace transforms, pQ(z,s) = 0 -st — P_Q(z,t)e dt, P1(z,s) 0 J -st P_ (z , t) e dt, (6) which, on substitution into (5), gives the two equations sPQ(z,s) - PQ(z,0) = -XPQ(z,s) + \XP1(z,s) - V Z ^ - P Q (z,s) (7) and ;Pi(z,s) - P (z,0)-= - (X+\i~^z) P 1( z ,s) + XPQ(z,s) + V T ^ P 0 ( Z , S ) if we assume A., u, and v, are independent of time. In principle, we could solve the first equation for P (^z,s), substitute this into the second equation, and obtain a single first order partial differential equation for Po (z,s). At best the resultant equation is of daunting complexity. This (z,t). Here, P_0(z,t) = zvQ(z,t) and P_±(z,t) = vi(z,t) - vQ(z,t) It would appear possible to evaluate the necessary integrals analytically if an iterative method was employed to solve these equations, even if some or all of A, u , and v are non-constant, though the amount of work required may increase rapidly with successive iterations. An as yet unresolved problem here is that in employing the successive approximation approach described in Zachmanoglou [1976], we seem to get an expression for P (^z,t) involving negative powers of z, which contradicts equation (4). 35 . approach is not likely to be feasible when any or all of X , u , or v , are functions of time, since then equations (7) would become even more complicated, and potentially involve higher order derivatives of the Laplace transforms. Previous Work on Time-Dependent Queues Even though comparatively little attention has been given to the time-dependent behaviour of queues in the literature, the number of papers which have appeared on the subject is quite large. No attempt will be made here to give a comprehensive review of the entire history of the subject. However, it is worthwhile to note the types of problems which have been solved successfully, and the sorts of approaches used. Among the early works on the subject was that of Morse [1958], who was able to obtain the time-dependent probabilities for the M/M/ l queue with constant service and arrival rates using the probability generating function approach. These analytic solutions were very complicated, involving sums over terms containing Bessel functions. Later, Saaty [1961] succeeded in obtaining the Laplace transforms of the time-dependent probabilities for the M/M/s queue with constant arrival and service rates, but was able to invert the transforms only for the two-server case. In both works, the time-variation of the probabilities was a transient effect as the systems 'decayed' to their steady state. Similar approaches would not have produced corresponding solutions had the arrival and/or service rates themselves been functions of time. 36 . With the failure of frontal attack for more complicated queue models, an alternative approach was developed by Newell [1968a, b, c], with extensions more recently by Keller [1982] and Massey [1985]. Newell was interested in modelling the behaviour of a traffic system going through a rush hour during which the 'arrival rate' significantly exceeded the capacity of the system. Since the length of a queue experiencing more arrivals than departures on average will grow without limit, a steady state cannot be established. Newell recognized that in such a saturated queue, to a good approximation one could consider arrivals and departures to be a continuous process because the time intervals between successive such events are very small. This approach led to a kind of diffusion equation for the time-dependent distribution function for queue length. That equation could not be solved generally, but its form facilitated the drawing of various conclusions about the behaviour of the system. Keller [1982] employed a formal asymptotic analysis amounting to adopting a time scale over which changes in arrival and service rates are very small. Keller found that his analysis reproduced many of Newell's results. Massey has apparently put Keller's work on a more solid theoretical foundation, particularly in proving that several of Keller's results were in fact asymptotic to the exact solutions. These analyses are quite complicated, and neither Keller nor Massey attempt to compare the results of their approach with exact numerical solutions. All three authors deal with some variant of the M/M/l queue, though Keller states that his method of analysis applies to other cases as well. To extend this approximation approach to the retrial system would be quite a major undertaking. 37 . Kelton and Law [1985] have restructured the M/M/s problem in terms of the probabilities: P (n,i) = Pr{X = i k customers were already present at t = 0}, k n t h where X is the number of customers in the system at the instant the n 'new' n customer arrives, including the new arrival. Thus we are looking at the system at a monotonely increasing sequence of random, discrete times, rather than in a continuous time framework. Using this approach, Kelton and Law are able to compute the expected delay in the queue for each arriving customer exactly. This approach is worth examining for the retrial queue, but that will not be done here because our immediate goal is to produce a viable approach to describing the continuous time behaviour of the system. Numerical Approaches As far as being able to predict the continuous-time behaviour of actual single server retrial queues, there seems to be little hope of finding computationally feasible exact analytic solutions. Turning attention to strictly numerical methods, we find that while they introduce difficulties of their own, these can be minimized at a price. None of the factors of the retrial queueing problem which dash hopes of an exact analytic solution are of any great consequence as far as implementation of numerical methods are concerned. If we adopt a specific 38 . numerical approach to solve the systems (1.1.3a,b), the algorithm will simply require us to be able to supply numerical values of A , u , and v , for any specified value of t. Whether these parameters are constants throughout or are generated by complicated formulas involving transcendental functions or worse, is effectively irrelevant to the structure of the numerical algorithm.5 The available algorithms for the numerical solution of ordinary differential equations easily number in the hundreds if not thousands. For the work discussed following this section, a modification of the Runge-Kutta-Felberg-45 or RKF45 method [Felberg, 1970] based on Algorithm 6.3 of Burden, Faries, and Reynolds [1978] was used exclusively (see Figure 1.3). This method generates both a so-called fourth order estimate and a fifth-order estimate of the solution efficiently at the end of each interval in t. There are good reasons for considering the difference between these two solution estimates as an acceptable estimate of the numerical error in the fourth-order estimate (see Burden, Faries, and Reynolds [1978], section 6.5). Since the fifth order estimate is the one retained, this difference between the two estimates is generally a very pessimistic estimate of the actual numerical error in the adopted solution. The main reason for choosing this method was that it not only produces an estimate of the numerical error being incurred, but it can be implemented so that 5Most numerical methods for integrating systems of differential equations arc based on some combination of polynomial extrapolation and interpolation. Changing the functional form of various coefficients in the equations would not require modifying the structure of the algorithm, though it would likely affect the detailed performance of the algorithm. As these coefficients, and so the solution, become less like a smooth low-order polynomial, the polynomial extrapolation/interpolation approach may become less effective or efficient. F I G U R E 1.3 T he Runge-Kutta-FehIberg-45 Algorithm The following is the implementation of a Runge-Kutta method as developed by Fehlberg (see Burton, Fairies, and Reynolds [1978], section 6.3) as used exclusively to solve systems of differential equations for the calculations presented in this thesis. The algorithm produces a fourth order and a fifth order estimate of the solution, requiring a total of only six derivative evaluations per internal step. The difference between the two estimates is taken as an upper limit on the error in the fifth order estimate, and internal step lengths are adjusted automatically to meet a preset tolerance on this error limit. I N I T I A L I Z A T I O N : The algorithm computes an approximation to y(b) = (yi(b), y 2 (b), y (b)} with n an integer greater than or equal to 1, given a function subprogram which"computes y'(t) = {y.'(t), y »(t), y '(t)} = {f (t), f (t), .... f (t)} -1- 2 n 1 2 n = f(t) and the initial values y(a) = (y,(a), y (a), y (a)}. •L 2 n It is also necessary to specify e = the desired maximum absolute difference ( e > 0) between the generated fourth and f ifth order estimates of y i (t) , i = 1, 2, .... n, over the range a ^ b . 1 /4 S T E P 1: Estimate the initial internal step length: h = (e) Set t Q = a. Set (w ) = y.(a), i = 1, 2, n . O i l S T E P 2: Set t = t Q + h . If t > b, reset t = b, h = b - t . o S T E P 3a: Set (k± ) ± = h-f (tQ,w o), i = 1, 2, n . S T E P 3b: Set (k ) - h*f (t + 4i, w + ^ k , ), i = 1, 2, .... n . 2 i i o 4 o 4 x S T E P 3c: Set (k ). = h«f (t + |"h, w + ~ k + ), i = 1, 2, n. 3 l i o ° o 3 2 1 32 Z S T E P 3d: S e t ( k 4 ) . = h . f . ( t 0 + i | . h , w o + f f f f ^ - f f f f k 2 + f f f ^ ) . i = 1, 2, n 40 . Figure 1.3 continued S T E P 3e: Set ( k - ) . - h « f . ( t + h, w + - 8k_ + 34f?-k - -JLilk ), 5 ' i i v 0 o 216~ 2 513 3 4104 4 i = 1, 2, n . S T E P 3 f : Set (kg ). - h- f . ( t o + ^, W Q - j | k , + 2 k 2 - | | | | k 3 + i M ^ 4104 4 i i - lc. ), i = 1, 2, n . S T E P 3 g : S e t w = w + + ^ + f f ^ k • I t . 4 o 216 1 2565 2 4104 4 5 5 S T E P 3 h : Set w = w + r ^ S + r l f f l * , + - + - | k c 5 o 13 5T. 12825 3 56430 4 50 5 55 6 S T E P 3 i : C o m p u t e r = max{ |(w g ) - ( w 4 ) j , i = 1, 2 n } /h S T E P 4: U p d a t e i n t e r n a l step l eng th : C o m p u t e 6 = 0.84 ( e / r j 1 ' 4 - 1 4 I f r is e f f e c t i v e l y zero (say r < 10 ), go to step 6. I f 6 < 0.1, r ep lace h by O. lh . I f 6 > 4.0, rep lace h by 4h. I f 0.1 < 6 < 4.0, rep lace h by 6h. S T E P 5: I f r > e , r e t u r n to step 2 (the f i f t h o rde r es t imate w has been re jec ted) . 5 S T E P 6: We h a v e r ^: e , t he re fo re the c u r r e n t f i f t h o rde r es t imate at t is accep ted . I f t = b, go to step 7. C o p y w 5 i n t o w . Set t Q e q u a l to t. R e t u r n to step 2. S T E P 7: T h e n u m b e r s ( w 5 )^ , i = 1 ,2 , n, es t imate y^ (b), i = 1 ,2 , n , w i t h l o c a l e r ro r es t ima ted to be less t h a n e ( the s u m o f the va lues o f r c o m p u t e d f o r each i n t e r n a l step g ives a rough es t imate o f the g l oba l e r r o r o v e r the i n t e r v a l [a,b]). 41 . internal step sizes are adjusted automatically to keep this numerical error below a pre-set threshold. This is a very important feature since we wish to be able to control numerical errors entering computed results so that discrepancies between results obtained from the exact equations (1.1.3a,b) and our approximations can be reliably attributed to the nature of the approximation rather than to carelessness in handling the equations.6 Questions of efficiency of implementation and speed of computation are important, but secondary to the objectives here, so no studies have been done to determine how the RKF45 algorithm performs compared to other types of algorithms. Truncation of Systems of Differential Equations No numerical algorithm can cope with an infinite system of differential equations, and so the system (1.1.3a,b) had to be truncated. In practice, this truncation is done fairly arbitrarily. As calculations reported in this thesis were carried out, the sum of the computed probabilities was subtracted from one at each point where the solution was computed to give an estimate of the residual probability 'lost' due to truncation. If that residual probability became greater than I.OxlO-6, the calculation was repeated with the system of equations 6In the calculations reported in this thesis, this error threshold was taken to be 0.00001 throughout on the probabilities pg (t). Except in rare cases when the expected number of customers in the system can become very large, this ensures that the computed first and second moments at least should be accurate to three or four significant digits. Of course, it must be remembered that this error estimate threshold is an upper limit, and the actual errors incurred may be very much smaller than this. 4 2 . truncated at a larger dimension. It would not be difficult to implement the algorithm to permit the dimension of the truncated system to be adjusted from step to step automatically in order to keep the residual probability below a pre-set value, but this was not done. Also, at each step, the residual probability was added into the highest index probability computed (or pro rata into the two highest index 'n' probabilities for the retrial queue) to reduce system truncation errors in computed moments, and the buildup of these errors from step to step as the complete specified range of time was covered. Truncation of the system to a finite size is justified on the basis that the equations dropped are effectively all identities, 0 = 0. In order for this to be true, we must retain equations for all states of the system which have a probability greater than some criterion of negligibility. A problem with the retrial system is that for a given utilization of the queue, the number of equations is roughly twice or more the number which would have to be included for the simple M/M/l-type queue, since the probabilities are indexed both by the number of customers in the system and by whether the server is idle or busy. Also, we can have a relatively lightly utilized retrial system ( P = A/u significantly smaller than 1) but if the retrial rate, v , is also small, there can be significant probabilities of finding larger numbers of customers in the system and so the dimension at which (1.1.3a,b) can be truncated has to be correspondingly greater. Steady-state probabilities illustrate this problem. The p for the M/M/l queue n and the p„ and p, for the M/M/l/1 retrial queue normally exhibit a maximum On In value for some value of the index n and then decrease with increasing n, eventually to zero. A relative indication of the necessary dimensions of truncated systems of the differential equations can be obtained by noting the number of elements of these sequences of steady-state probabilities with value at or above a specific threshold e. For example, take the "negligibility" threshold e as 1.0 x 10 6 When P = 0.5 and c = 1, we find that p < e for n > 19, whereas p < e for n > 26 and n On p. < e for n > 20. So, while only the first 18 steady-state probabilities for the i n M/M/ l queue are non-negligible by this criterion, there are 25 + 19 = 44 non-negligible steady-state probabilities for the corresponding retrial queue. Using p = 0.5 and c = 10, the number of non-negligible steady-state probabilities for the retrial queue is 86 by this criterion (and is still 18 for the M/M/l queue, since the p depend only on p). n Alternatively, as c gets nearer to zero (indicating retrial rates are much greater than new customer arrival rates), the number of non-negligible P0n's for the retrial queue becomes similar to the number of non-negligible p 's for the n simple M/M/ l queue, and the number of non-negligible p, 's approaches zero. 1 n For example, when p = 0.5 and c = 0.1, there are 18 non-negligible p and 10 On non-negligible p ^ for a total of 28 non-negligible steady-state probabilities for this retrial queue (compared to the total of 18 for the corresponding simple M/M/l queue). Of course, for time dependent calculations, the truncated systems may have to have higher dimensions than indicated here to ensure that at times before a steady-state is approached, no non-negligible probabilities are excluded from the calculation. However, this rough analysis gives an indication of the differences in dimension of the system of differential equations which must be solved in computing the time-dependent behaviour of the two queues starting with systems (1.1.3a) and (1.1.3b). This problem is of significance in computing moments of the distribution, since errors in high index probabilities (including the 'error' of neglecting them) are magnified. Thus, for retrial queues, we will find that it is necessary to truncate the system of differential equations to be solved numerically at dimensions two or three or more times the size as would be necessary for a simple M/M/l queue. In the calculations reported in this thesis, necessary dimensions for truncations of the systems (1.1.3a,b) were typically in the range of 50-100 though a few calculations required handling systems of 130 equations to attain the desired accuracy.7 Details are given in Appendix A. 7Again, the actual computed results were considered more important than how fast or efficiently they were computed. Programs were written in FORTRAN and calculations were carried out on an IBM mainframe using the WATFIV compiler for convenience. With the numerical error threshold set at 0.00001 over a step of length At = 0.5, of the order of 0.1 to 1 second of cpu time was required for integration of systems of equations with dimensions in the 50 - 100 range. On an IBM PC/XT without a numerical coprocessor and using the Microsoft FORTRAN compiler, execution times for calculations like these would require of the order of one to ten minutes per step of length A t = 0.5. These numbers are of order of magnitude only since no attempt was made to write FORTRAN code optimized for execution speed, and since the algorithm will break 45 . It is noted that there are multi-step methods for the numerical solutions of systems of differential equations which also permit dynamic step length/error control (see for example, Algorithm 6.5 in Burden, Faries and Reynolds [1978], based on a fourth order predictor-corrector approach). These algorithms require a single-step routine like the Runge-Kutta methods to start them off, and they would be worthwhile evaluating if a search was being made for the most efficient algorithms. Similarly there are higher order single-step methods (see, for example, Shanks [1966], who gives several eighth-order Runge-Kutta methods). These higher order formulas in principle give greater accuracy for a given step length, or, more importantly, they allow longer step lengths for given error tolerances. Higher order formulas are mainly used to avoid unacceptable round-off error due to the finite precision with which a digital computer stores numbers,8 which grows as the step length gets shorter. an interval of length At = 0.5 up into the required number of smaller subintervals to meet the error tolerance. The number of such subintervals required may vary widely with the value of t for certain functional forms of A(t), u(t), and v(t). 8As distinguished from what we have been calling the numerical error, which is due to the approximation of the function of interest by a low order polynomial over small intervals. 46 . 1.4 Closure Approximations for M / M / s Queues Snafu Law #1: Given any problem containing n equations, there will be n+1 unknowns. (1001 Logical Laws.... Peers & Bennett) Details of the Rothkopf and Qren Closure Method The rationale behind the use of so-called closure techniques to the approximation of solutions of infinite systems of differential equations is well described in Rothkopf and Oren [1979]. The essential feature in applying such a technique in a queueing context or in any other context is to make one or more assumptions about the form of the functional relationship between various quantities with the result that an infinite system of equations can be reduced to a small number of equations determining a few quantities of direct interest. We illustrate this approach as applied by Rothkopf and Oren to the M/M/l queue (with contributions by earlier authors as cited in their paper). Defining p^ (t) as the probability that there will be n customers in the M/M/l/°° system at time t, including any customer in service at that time, we can proceed much as in section 1.1 above to obtain the equations " f t - p o ( t ) = - X p o ( t ) + » p i ( t ) (1) (t) + up . ( t ) " p ( t ) = - ( X + U)p (t) + Xp . . r r „ . i dt n n n-1 n+1 n = 1 , 2 , 3 , . . . Here, ^ is the average arrival rate and u the average service rate. In order to compute any characteristic of the system, such as the mean or the variance of 47 . the number of customers at time t, it is necessary in principle to solve the infinite system (1). The difficulties to be encountered have been described in detail in the preceding section of this thesis. What Rothkopf and Oren propose is the following. If each of the equations in (1) is multiplied by n, then all equations summed over n, we get first oo oo oo oo E 4 r ^ - p (t)}= -(X+u) E n-p (t) + A E n-p (t) + y E n-p (t) n d t n n n n n - 1 n n+l n=0 n=0 n=0 n=0 With some rearrangement, this can be written as a m (t) = A - y {l - pn ft; } (2a) dt " r 0 oo where m(t) = E np (t) is the mean number of customers in the system at time t. n n = 0 2 In a similar way, the equations can be multiplied by n , summed, and using the 2 fact that the variance of the number of customers in the system = v(t) = E(n ) -2 m , so that cl d 2 dm v(t) = E(n ) - 2m at dt dt we get v(t) = A + y - up. ft; {2m(t) + 1} . (2b) dt 0 Equations (2a) and (2b) contain three unknown functions: PQiX)-. m(t) and v(t). There is no obvious way to eliminate one of the three exactly. Rothkopf and Oren then proposed an approximation to eliminate PQ(t) in terms 48 . of m(t) and v(t), so that equations (2a) and (2b) would become a system of two differential equations in these two unknown functions only. This involves introducing a so-called "closure assumption", in the form of an approximate relation (or "closure formula") between the three quantities. We will use the terms closure assumption and closure formula interchangeably. Rothkopf and Oren noted that the steady-state probabilities for the M/M/l/<» system were a geometric distribution, which was a special case of the negative binomial or Pascal distribution. They noted further that for the Pascal distribution, one could write (3) They then adopt this functional relationship to eliminate P Q(t) from equations (2a) and (2b).1 The choice of this functional relationship is prompted and justified in part by the resultant equations correctly reproducing the steady-state values of m and v when the derivatives in (2a,b) are replaced by zeros. However, as Rothkopf and Oren note, "the key to the success of a closure method is its closure assumption, while the test of any such assumption is the empirical usefulness of the resulting equation." The success of the approximation will depend on how close to a negative binomial distribution the p (t) are over all finite values of t. Intuitively, we n 1A modification of this formula is given for use in the case that v=m — this usually only occurs at the beginning of a calculation when the system is assumed to be empty so that m(0) = v(0) = 0. 49 . would expect an approach like this to perform most poorly when at least some of the probabilities, p (t), are changing rapidly with time, since then it is plausible n that the deviation from steady-state relationships between these quantities is likely to be greatest. Numerical Results The largest part of Rothkopf and Oren [1979] consists of the presentation and discussion of various computed results. Their assessment of this closure approximation as giving an acceptably accurate practical computational method is generally favorable. However, the results they present are mainly graphical and so cannot be subjected to detailed quantitative analysis. A set of standard calculations were performed here, to which corresponding results for retrial queues will later be compared. For one set of calculations, the mean arrival and service rates were assumed constant, so that over time, the closure solutions obtained should describe the transient decay of the system to steady state. At the other extreme, arrival rates were assumed to have the functional form \(t) = X {1.1 + sin(^p-) } . (4) where A^ is a constant and P represents the period of oscillation of the arrival rate (taken to be 10 time units throughout this thesis). We will refer to the first type of calculations as constant parameter calculations, and to the second type as 50 . oscillatory parameter calculations. The number '1.1' was used in (4) to avoid A(t) becoming zero. This seemed to cause numerical algorithm problems. For both types of calculations, the service rates were held constant in time. The RKF45 algorithm was used to solve all systems of differential equations here to ensure that the discrepancies found between the quantities computed from the effectively exact solutions of the system (1) and those obtained by solving (2a,b) with (3) are due to errors in the closure approximation, and not due in any significant way to the numerical approximation of the solutions of the equations. Rothkopf and Oren used the simple forward Euler method for the work reported in their paper. While the adequacy of this simplest of all algorithms may not be strenuously questioned for very well-behaved systems of only two differential equations at the level of accuracy required, it is not necessary to adopt an algorithm so prone to numerical error, given the power of commonly accessible computing machinery these days. For the constant parameter case, graphs depicting the relative errors2 in P^CO, m(t) and v(t) respectively as calculated by the closure approximation are given in Figures 1.4, 1.5, and 1.6. Calculations were done for p = 0.1, 0.3, 0.5, 0.7, 0.9, and 0.95 (where p = A / u is a measure of the mean utilization of the server) to compare the accuracy of the closure approximation for queues approaching steady-state with a range of characteristics. Only the results for the three largest throughout this thesis, errors are always computed as 'approximate - exact', and relative errors are that quantity divided by the exact value of the quantity being considered. 51 . values of p are graphed to avoid clutter. In all calculations, it was assumed that the system was empty at t=0 (that is, po(0)=l, m(0)=0, and v(0)=0) and in all cases, the closure approximation solution converged to the steady-state solution p (°°) = 1 - p , m(°°) = P , v(<*>) = 2 — ' (5) P ( 1 - 0 ) monotonically or with small oscillations about a monotonically converging path. As can be seen from Figures 1.4 - 1.6, the relative errors in each of the three quantities computed with this closure approximation start out at zero, increase in magnitude to a maximum, and then decrease slowly and in a somewhat oscillatory fashion to zero. On the scales adopted in the figures, the relative error in v(t) when P = 0.1 shows up only as a barely noticeable blip near t = 1.5. Relative errors increase in magnitude with p . For PQ(t) and m(t), the maximum relative errors when p = 0.95 would appear to be under 4%, whereas for v(t), this maximum is nearer to 15%. When p is significantly less than 0.9 or 0.95, these maximum relative errors are much smaller. These constant parameter calculations are not difficult. Even when p = 0.95, the RKF45 algorithm was able to produce a solution to a truncated system of 80 equations as in (1) with error estimates in the range 10" to 10" using step lengths between 0.5 and 1.0. Nevertheless, some fairly significant changes in values are occurring (for example, when p = 0.95, p (1.5) = 0.251, a drop of 75% 0 in value in 1.5 time units, m(10) = 5.62 against m(0) = 0, and v(25) = 53.4 against v(0)=0, so that over this range, the derivative of v(t) averages more than 2.0). It 52 . is encouraging to see that the much simpler closure approximation does reproduce these changes rather accurately. It is more difficult to get an overall impression of the "pattern" of accuracy of the closure approximation for the oscillatory parameter calculations. Figures 1.7a, 1.8a, and 1.9a show the relative errors in P0(t), m(t) and v(t) respectively for p Q = X^/y values of 0.5, 0.7, and 0.9. Now the relative errors as a function of t are oscillatory with rather complicated shapes. It is somewhat disquieting to note that the relative errors for larger values of Po attain rather large magnitudes for certain ranges of values of t. For = 0.9 and in the range t = 0 to t = 25, the relative error in PQ(t) approaches 0.9 three times (always near the mid-period where X(t) is greatest and so PQ(t) is smallest), and the error seems to be exhibiting an increasing sequence of maxima like 0.5, 0.7, .... at the beginning of each period (where the value of A(t) is midway between its greatest and its smallest values). Similar types of behaviour of the relative errors in m(t) and v(t) are also observed, but with one important qualification. The extreme values in the relative error plots are not greatest for the largest values of p Q over the range of values of p^ represented in Figures 1.8a and 1.9a. For smaller values of p^ like 0.1 and 0.3, the relative error extremes for m(t) and v(t) are smaller than they are for P Q = 0.5. For the values of pQ for which calculations were done, p^ = 0.5 leads to the greatest extremes in the relative error in m(t) and p^ = 0.7 leads to the greatest extremes in the relative error in v(t). 53 . The reason for this shape of the relative error plots can be seen from Figures 1.7b-1.7d for p o(t), Figures 1.8b-1.8d for m(t), and Figures 1.9b-1.9d for v(t). In these figures, exact and approximate values are plotted for each of p Q = 0.5, 0.7, and 0.9. It is seen that the correspondence between the approximate and exact curves is not all that bad. Though the gaps between the two are readily observable in each instance, in most of these plots the approximate curve tracks the exact curve rather well. These quantities are all oscillatory, reflecting the oscillatory mean arrival rate. Near the minimum values in each curve, where the graphs reverse direction, the gaps between the exact and approximate curves widen. As a result, absolute errors tend to be largest where the values of these quantities are smallest, leading to larger relative errors in those regions. What is happening is that the approximation is poorest where the queue is not only changing rapidly, but where the rates of change of the queue are changing rapidly (which, as remarked earlier, is precisely where we would expect these closure techniques to break down since they are based on steady state relationships). This happens to coincide at times with regions where the quantities being computed attain their minimal absolute values. The instantaneous relative error is thus not a very good indication of the overall accuracy of the approximation, since its graph will be dominated visually by small regions of relatively large errors in small quantities. 54 . Period-Mean Values and Period-Mean Errors Perhaps a better approach is to get some estimate of the 'average' relative error over a single period of oscillation. Using the known values of P Q(t), m(t) and v(t) in steps of At = 0.5 forming the basis for the plots in Figures 1.7 - 1.9, we can apply Simpson's Rule to estimate the integral of each of these three quantities and their absolute errors over 'periods' of length 10 time units, and so obtain their mean value over that time interval. 3 Dividing the mean absolute error over such a period by the mean value of the quantity itself then gives a relative mean absolute error, which will be more representative (it is hoped) of the overall discrepancy between the exact and approximate solutions than are the point-by-point relative errors. This has been done in Table 1 for each of = 0.1, 0.3, 0.5, 0.7, and 0.9, and in each case, for the two ranges t = 0 to t = 10 and t = 10 to t = 20. It is worth taking note of several features of the information in Table 1. First, the relative mean absolute errors do seem to reflect the degree of agreement between the graphs of approximate and exact values of these quantities. Except for m(t), these relative errors increase with the value of p as happens in the constant parameter case. 3In effect, we would be using a 21-point Simpson's formula on an interval of length 10. When the function sin(2-rrt/10) is integrated on the interval [0,5] using points spaced at intervals of 0.5 in width, the error in the integral obtained is 0.000174 on a value of 3.1830989 or about 0.0055%. While the quantities we arc trying to integrate above are not simple sine functions, their shape is not so different that we would expect the Simpson formula to perform drastically poorer. It would appear that this approach will produce integrals of adequate accuracy for our purposes here. 5 5 . A good example of the detailed correspondence between the Figures 1.7 - 1.9 and Table 1 occurs for the variance when aQ = 0.9. In Figure 1.9d, it is clear that the two curves are moving apart steadily as t increases. From Table 1, we see that the relative mean absolute error over the interval [0,10] is 0.1953, whereas over the interval [10,20] it is 0.2470. One does not get the same sense from the corresponding plot of the point-by-point relative error in Figure 1.9a, which consists of a sequence of peaks at about t = 9.5, 19.5, ... . In fact, from Figure 1.9a, one probably gets the opposite impression since the peak at t = 19.5 has the value 0.457 against a value of 0.493 at t = 9.5, and we tend to perceive peak height as most significant in graphs of errors. A closer look indicates that while the peak at t = 9.5 is slightly higher than the peak at t = 19.5, the former peak is also sharper than the latter, and so the average error over an interval containing the first peak could well be smaller than the average error over an interval containing the second peak. In order that Figure 1.9a appear to agree with Table 1, we have to focus on both the heights and the widths of the peaks in Figure 1.9a. In work dealing with the oscillating parameter case for the retrial queue below, we will rely heavily on this relative mean absolute error to indicate the goodness of a particular closure approximation. What is most striking about the relative mean absolute errors listed in Table 1 is that they are so small in general. The mean queue length, m(t), seems particularly well determined by this closure approximation. Errors in v(t) and 56 . p Q(t) only become disquietingly large when is 0.7 or larger. Because of the factor '1.1' in equations (4), this corresponds to a time average ratio of interarrival to service times of 0.77 or greater. Also included in Table 1 are values of the relative errors found in computing mean values of P^it), m(t), and v(t) over 10-unit time periods. These relative errors cannot be greater than the relative mean absolute errors discussed above, and in fact, will be equal to those only when the approximation produces uniformly positive or uniformly negative errors over the period in question. When actual approximation errors change sign, we find that errors in the computed mean value of these queue characteristics can become very small. In the case of p Q(t), the relative mean absolute errors over the time intervals [0,10] and [10,20] are 23% and 38% respectively when p^ = 0.9, whereas the relative error in the mean value of PQ(t) over each of these intervals are only 5% and 7% respectively. The difference between these two types of error for P^ CO are even greater when p = 0.5 . 0 On the other hand, for v(t) with QQ = 0.5, 0.7, and 0.9, both types of relative error are identical because the Rothkopf and Oren closure method consistently overestimates v(t). Relative errors in computed period-mean values are useful if the primary interest is to compute such mean values of queue properties. However, the relative mean absolute error is a more sensitive single number measure of the overall accuracy of an approximation method here, since it takes into account all errors with no cancellation between positive and negative errors in different subintervals of time. 57 . The influence of the initial state of the system shows up in the difference between the averages computed in Table 1 applying to the time interval [0,10] compared to the time interval [10,20]. It is expected that after a sufficiently long time has elapsed, all of the properties of the queue will attain a sort of periodic steady state. Koopman [1972] refers to this periodic solution, with the statement that it is a consequence of classical theorems that only one solution of the basic system of differential equations will have the periodicity of the mean arrival rate and all other solutions will approach this periodic one exponentially with increase of time. He notes that this periodic solution is analogous to the steady state solution approached exponentially with increasing time in the constant parameter case [see Odoni and Roth, 1983]. For = 0.1, this periodicity is easily observed to six figures or better in the numerical results b y the time t = 15. For other values of p ^ , the graphs demonstrate a tendency towards establishing such a periodicity, but it has not occurred to similar precision over the time ranges for which solutions were generated. Extensions and Other Approaches Finally, we note that Rothkopf and Oren [1979] also present a closure method for use with simple multiserver queues. In the case of s > 1 servers, they require expressions for p (^t), p (t), p (t) in terms of m(t) and v(t) in order to 1 2 s -1 'close' their pair of equations. They compute these probabilities using the negative binomial distribution, which is no longer quite correct. This means that the closure equations no longer produce the correct steady state solution in the 58 . constant parameter case. Calculations indicate that the error due to this problem is small, and Rothkopf and Oren include a table of corrections to be applied (which depend on the number of servers and the value of A/u ). The errors they observed in the corrected results obtained using their approximation were apparently of an acceptable magnitude. Clark [1981] has also considered closure approaches to the multiple server problem, using what he calls a 'surrogate distribution'. He found that his approach generally was more accurate than the method of Rothkopf and Oren, though it required the solution of a system of five differential equations as opposed to their two. Since this thesis will not deal with the multiple server retrial queue, it is not necessary here to examine these works in greater detail. TABLE 1: Period-Mean Values of Pa(t), m(t) and v(t) and Their Relative Errors for the M/M/l Queue When the Rothkopf and Oren Closure Method is Used. The first column in each triple gives exact mean values over the time intervals indicated. The second column in each triple gives the relative error in the period-mean computed using values of the quantity given by the Rothkopf and Oren closure method. The third column in each triple gives the relative mean absolute error in the Rothkopf and Oren approximation in each case. time -fj i n t e r v a l m(t) v (t) 0.1 0.3 0.5 0.7 0.9 [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] .8927 .8900 .6790 .6700 .4677 .4501 .2722 .2357 .1304 .07323 -7.246X10 1.231x10 -3.194x10 -1.526X10 -7 - 6 .003103 .0004911 .03590 .01114 -.05368 -.06970 1.640x10 1.494x10" .001772 .001755 .02301 .02684 . 1 2 4 6 . 1 4 9 6 .2253 .3817 .1271 .1304 .5965 .6162 1.7427 1.8373 4.1337 4.7874 7.9358 11.3146 3.934x10 4.858x10" -.002217 -.002126 -.01787 -.01126 -.04107 .03815 -.03627 -.01401 .004351 .004146 .008117 .007849 .03886 .03665 .05356 .05179 .03627 .01536 .1058 .1549 1.1086 1.1518 4.9362 5.3220 14.4653 18.9005 28.5147 55.2801 7.231x10 6.350x10" .02674 .02564 .1547 .1607 .2638 .3799 .1953 .2470 .002532 .002411 .03398 .03292 .1547 .1607 .2638 .3799 .1953 .2470 On ^5 60 . 0 5 10 15 20 25 FIGURE 1.4 Relative Error in pQ(t) Computed Using Rothkopf and Oren's Closure Method for the Simple M / M / l Queue with Constant Arrival and Service Rates. The relative error is computed here as (approximate value - exact value)/(exact value), plotted here for various values of p. The relative error when P = 0.3 is discernable on the scale above only for small values of t, and that when p = 0.1 is not distinguishable from the t-axis on this plot. The service rate u is fixed at 4.0/unit . time in each calcultion and the system is assumed empty at t = 0. Details of calculations to obtain 'exact values' are given in Appendix A. 61 . F I G U R E 1.5 Relative E r r o r in m(t) Computed Using Rothkopf and Oreii 's Closure Method for the Simple M / M / l Queue with Constant Arrival and Service Rates. T h e re l a t i ve e r ro r is compu ted here as (app rox ima te v a l u e - exac t va lue ) / (exac t va lue) , p lo t ted here f o r va r i ous va lues o f the u t i l i z a t i o n p . T h e re l a t i ve e r ro r w h e n P =0.1 shows up o n l y as a s l igh t t h i c k e n i n g o f the t -ax is f o r s m a l l va lues o f t on th is plot . T h e se rv i ce rate u is f i x e d at 4 .0 /un i t t ime i n each c a l c u l t i o n a n d the sys tem is assumed empty at t - 0. Because m(0) = 0, the re l a t i ve e r ro r is not d e f i n e d at t = 0 (but at t = 0, the abso lu te e r ro r i n m is zero by d e f i n i t i o n here). De ta i l s o f c a l c u l a t i o n s to o b t a i n 'exac t va lues ' are g i ven i n A p p e n d i x A . 62 . 0.15 FIGURE 1.6 Relative Error in v(t) Computed Using Rothkopf and Orcn's Closure Method for the Simple M / M / l Queue with Constant Arrival and Service Rates. The relative error is computed here as (approximate value - exact value)/(exact value), plotted here for various values of the utilization p. The relative error when p = 0.1 is-barely discernable along the t-axis for small values of t on this plot. The service rate u is fixed at 4.0/unit time in each calcultion and the system is assumed empty at t = 0. Because v(0) = 0, the relative error is not defined at t = 0 (but at t = 0, the absolute error in v is zero by definition here). Details of calculations to obtain 'exact values' are given in Appendix A. 63 . 0 5 10 15 20 25 0 5 10 15 20 25 (c) p Q = 0 . 7 (d) p 0 = 0 . 9 FIGURE 1.7 Comparison of Exact pQ(t) and pQ(t) Computed Using Rothkopf & Oren's Closure Method of the M/M/ l Queue in the Oscillatory Parameter Case. Figure 1.7a gives the relative error in the closure approximation p0(t) for pQ = 0.5, 0.7, and 0.9. Figures 1.7b,c,d plot the exact P O ( 0 along with the closure approximation (t) for pQ = 0.5, 0.7, and 0.9, respectively. In each case, the plot of the approximate p0(t) tends to swing wider at the maxima or minima relative to the plot of the exact pG(t). Here the period of oscillation is 10 time units and the system is assumed to be empty at t = 0 (i.e. Po(0) = 1.0). F I G U R E 1.8 Comparison of Exact in(t) and m(t) Computed Using Rothkopf & Oren's Closure Method of the M/M/l Queue in the Oscillatory Parameter Case. Figure 1.8a gives the relative error in the closure approximation m(t) for p = 0.5, 0.7, and 0.9. Figures 1.8b,c,d plot the exact m(t) along with the closure approximation m(t) for p Q = 0.5, 0.7, and 0.9, respectively. At least in part because m(t) generally increases with p Q, the magnitude of relative errors generally decreases with p Q as shown in Figure 1.8a. Here the period of oscillation is 10 time units and the system is assumed to be empty at t = 0 (i.e. m(0) = 0.0). 65 . (c) p 0 = 0.7 (d) p 0 = 0.9 FIGURE 1.9 C o m p a r i s o n of E x a c t v(t) a n d v(t) C o m p u t e d U s i n g R o t h k o p f & O r e n ' s C l o s u r e M e t h o d of the M / M / l Queue in the O s c i l l a t o r y P a r a m e t e r C a s e . Figure 1.9a gives the relative error in the closure approximation v(t) for pQ= 0.5, 0.7, and 0.9. Figures 1.9b,c,d plot the exact v(t) along with the closure approximation v(t) for p Q = 0.5, 0.7, and 0.9, respectively. The approximate v(t) produces the upper curves throughout (note that the relative errors over the time range shown in Figure 1.9a are positive). Here the period of oscillation is 10 time units and the system is assumed to be empty at t = 0 (i.e. v(0) = 0.0). 66 . CHAPTER 2 Closure Techniques for the Time-Dependent Single-Server Retrial Queue 'Write that down', the King said to the jury, and the jury eagerly wrote down all three dates on their slates, then added them up, and reduced the answer to shillings and pence. (Alice in Wonderland. Lewis Carroll) 67 . 2.1 The Basic Equations Our little party of travellers awakened next morning refreshed and full of hope (The Wizard of Oz. L. Frank Baum) Retrial Queue Characteristics of Immediate Interest We wish to do now for the M/M/l/1 retrial queue what was done for the simple M/M/ l queue in section 1.4. The retrial queue is more complicated and so the system of differential equations (1.1.3a,b) contains more detail than does the corresponding system (1.4.1) for the M/M/l queue. As a result, the variety of potential closure techniques is much greater for the retrial queue. The characteristics of the retrial queue of greatest likely interest are: oo (1)11 (t) = Pr(server idle at time t) = Z p (t) (2) m(t) = mean number of customers in the system at time t, including any customer i n s e r v i c e oo CO E n = 0 n • p On (t) + E n = 0 (n+1)p In It) (3) v(t) = variance from the mean number of customers in the queue at time t 6 8 . •• I n 2 p n (t) + Y. (n + l ) 2 P i (t) - { m ( t ) } 2 n On . ^ In n=0 n=0 c o ( t ) + ^ i ( t ) ~ ^ m ( t ) }2 (4)Y0(t) = conditional mean number of customers in the system at time t given that the server is idle at that time * 0 ( t ) The quantities m (t), m7(t), £ (t), 5 (t), and so on are defined in section 1.1. 0 0 1 The four quantities above are all determined by second or lower order moments of the probability distribution for the state of the retrial queue. Derivation of the Basic Equations The first step in developing these closure techniques is to obtain equations for the derivatives of the quantities of interest. This is easily done by summing appropriate multiples of equations (1.1.3a) or (1.1.3b). For example, we have p' (t) = - (X + nv)pn (t) + VP-, ( t ) , n = 0 , 1 , 2 , . . . . (1.1.3) On On In Summing both sides of this equation over all values of n, n = 0, 1, 2, we get 69 Z p' (t) = n = 0 •X Z p (t) - V Z n - p o ( t ) + y Z p 7 (t) On On ._ „ I n n = 0 n = 0 n = 0 d d t n = 0 P 0 n ( t ) \(t) m Q ( t ) or n ( t ; = -xno(t) - vmQ(t) + u n 1 c t ) (i) Similarly, we can first multiply equation (1.1.3a) by n on both sides before summing, to get ^ n-p (t) n = 0 -X Z n - P n (t) - v Z n P n ( t ) + y Z n - P i ( t ) On . "On n " ln n=0 n=0 n=0 We notice that Z n-p^ ( t ) = Z {(n + 1) - I } p 7 ( t ) = Z (n + l ) p . ( t ) - Z p 7 ( t ) In n In . i n n I n n = 0 n=0 n=0 n=0 m 1 ( t ) - J[JL(t) , and, since n is independent of t on the left hand side, n = 0 Thus, we obtain Ln=0 Z n . p 0 n ( t ) = m 0 ( t ) m'Q(t) = -Xm (t) - V ? f t j + umi(t) - y l L j f t ; (2) We can proceed in the same way to derive equations for rn^'(t), m'(t), and so 70 . on. The derivative of Y Q(t) = m^ftj/njt) can be obtained by applying the basic formula for differentiating a quotient and using (1) and (2). Formulas for the derivatives of moments through second order are listed in Table 2. Additional Criteria The next step is to take sets of one or more equations from Table 2 and eliminate enough of the quantities they contain using relationships between those quantities suggested by steady-state formulas so that the resultant system can be solved uniquely. These relationships ("closure formulas") can often be constructed in a number of different ways, and the assessment of their goodness or badness will ultimately depend on the outcome of test calculations. However, it seems useful to impose two ground rules on the process to eliminate approximations which are very likely to fail the computational tests. (1) Whenever possible, exact relationships are to be used to eliminate quantities from equations in preference to presumed approximate relationships based on inspection of the steady state solution. Thus, to eliminate ^j(t) in terms of either 11^(0 or m 0(t) from equation (1), we would always use n (t) = 1 - II (t), since this relation must always hold for 1 0 valid probability distributions. At steady state, it is true that 1 1 ^ = p = m Q / c -However, we would not consider eliminating n (t) from equation (1) by writing T A B L E 2. Formulas for Moment Derivatives for the Single Server Retrial Queue. (i) n = -An + yn - vm 0 0 1 o (2) = -Amo - VC Q + yn^ - y l ^ (3) = AmQ + VC Q + A - ym^ (4) m* = X - u ( l - IT ) (5) = - X ^ Q - v n Q + y C 1 " 2y m ; L + y H 1 (6) ^ = X C Q + v n Q - \ir,i + 2ym + X (7) C' = 2(X - y)m + X + yII + 2ymQ (8) v' = X + y - yll 0(2m + 1) + 2 ymQ 72 . n^(t) = m0(t)/c because this would allow the resulting closure equations to have a solution violating the fundamental relation n o ( t ) + J [ i ( t ) = l f ( 3 ) whereas using n^(t) = 1 - TT^ Ct) would eliminate that possibility. Other inviolable relations here are m ( t ) + m Ct). = m(.t) f (4) and ^o(t) + ^ i ( t ) = ^ ( t ) • ( 5 ) (2)We will insist that proposed closure methods correctly produce steady state values as t becomes large in constant parameter calculations. This amounts to requiring that when derivatives in the closure equations are replaced by zeros, and A, p, and v are constants, the resulting system of algebraic equations give uniquely the correct steady state expressions for the quantities whose derivatives were originally present. For example, in section 2.3 below, we propose a closure technique for computing 11^ (0 and m(t) based on the two-equation system n ^ f t ; = -XRQ(t) + y l ^ f t ; - vm(t) + vm (t) (6) m' ct) = x - u n 7 ct) 73 . The steady state solution is obtained correctly by replacing the derivatives by zeros here to get 0 = - X I l ^ + u l l - Vm + Vm^ , o = x - u n ^ The second equation gives immediately n 7 = — = P , i v or, using 1 - ^ 0 = n , we get 1 - n^ = p , so that the solution can be rewritten as II = 1 - p , which is the correct expression for the constant parameter steady state case. Then from the first equation, we have = {-xn o + u n 2 > + m1 V = — ^ - { - X f l - p; + up} + m1 X p + m V 1 n , p e e p + i) p r e + i) = c p + = , 1 - p I - p as required. The substitution of the steady state expression for in the last line above is justified by the fact that we will eliminate m (^t) from the system (6) using steady state relationships to obtain a viable closure method for computing ^(t) and m(t). 7 4 . An example of a system of equations which does not satisfy this criterion is the three-equation system ft'(t) = - A l l f t ; + ]iH (t) - Vm(t) + vm (t) , m' (t) = A - u { l - I l p f t ; } , ( ? ) v' (t) = A + U - uRQ (t) {2m(t) + 1} + 2umQ(t) which might be used as a basis for a closure method for computing nQ(t), m(t) and v(t). When the derivatives on the left hand side are replaced by zeros, we obtain a system of equations which contains no reference to the variance v. The resultant system can be solved for the correct steady state expressions for say nQ 5 mQ, and m once 1 1 ^ = 1 - J[Q and m^ = m - m^ are eliminated (no approximations are used here). But, since the steady state expression for the variance cannot be determined in this fashion, we would reject this system of equations as a basis for a closure method for computing ^(t), m(t) and v(t). In fact, if one proceeds to. eliminate say nQ(t), m (^t), and mQ(X) from equations (7) in terms of 11^ Ct) and m(t), the transient time-dependent solution of the system for the constant parameter case is not too bad for ^(t) and m(t). The solution for v(t) is not at all good, however, and as t increases, v(t) does not apparently converge to a value anywhere near its actual steady state value.1 1As mentioned earlier, Rothkopf and Oren [1979] did propose a closure method for multiserver queues which did not converge in the constant parameter case to the correct steady state values. However, they were able to devise a table of approximate numerical corrections, which, when applied to the solutions produced by their method, did result in correct steady state values. These 7 5 . In a small number of test calculations done for five or six methods which violated either or both of these ground rules, the numerical results were invariably unrealistic, with moments negative at times, probabilities significantly exceeding unity at times, and clear convergence to incorrect steady state values for increasing values of t in constant parameter calculations. The first rule above amounts to avoiding use of inexact relationships when viable exact relationships are available. The second rule amounts to requiring that the methods proposed be capable of handling the simplest aspect of the simplest type of time-dependent queue behaviour if we are to consider it seriously for describing time-dependent queue behaviours in more complicated situations. Notation The closure methods examined in this thesis are described in detail in the following four sections. We categorize them according to the number of equations involved. For example, methods IA and IB are alternative approaches based on a single equation. Similarly methods 2A, 2B, 2C, would be alternative approaches based on systems of two equations. numerical corrections were then carried over to variable parameter calculations. In their case, Rothkopf and Oren were able to trace the source of the errors to specific assumptions in the underlying probability distribution on which they based their approach. Here, the error in v(t) at steady state is being approached is an artifact of the numerical algorithm employed to solve the system (7), and so general correction factors do not exist. 7 6 . 2.2 One-Equation Closure Methods Committees of twenty deliberate plenty, Committees of ten act now and then, But most jobs are done by committees of one. (old rhyme quoted by Laurence Peter in Why Things Go Wrong) It is possible to devise relatively simple single-equation closure methods for calculation of either n^ (t) or m(t). Though neither method produces highly accurate results, it is interesting to note how close their solutions are to the exact solution when so little detail is included in the computation. Method IA (IT (t)) To obtain a single equation for approximating no(t), we begin with the first equation in Table 2, n f t ; = -\nQ(t) + .yn j f t ; - vm (t) . f i ; n^ (t) is replaced here by the exact equivalent, 1 - H (t). At steady state, mQ = cp = c f J - II ; , (2) as seen from equation (1.1.12). Thus, as an approximation, we can write equation (1) as i T f t ; = -xn f t ; + y U - n Q f t ; } - vc { i - n Q f t ; } , or, simplifying, as n^ft; = - y n o f t ; + y - A O) 77 . There is a certain degree of plausibility to this formula. When the server is idle, n^ (t) = 1, then n^ 'Ct) = - ^, indicating that the probability of the server being idle will decrease with increasing new customer arrival rate. When the server is busy, ^ (t) = 0, then ^ '(t) = u - X , indicating that II^ (t) will tend to increase if the service completion rate is greater than the new customer arrival rate. The term - u«II0(t) in equation (3) simply indicates that one cannot expect the difference in service and arrival rates to lead to high rates of change in H-0(t) if the system is largely empty ( H (t) near 1). One significant deficiency in this formula is that it is independent of the retrial rate, v(t). In a way, this is what you would expect of a formula based on steady-state relationships. When the number of customers in the system is no longer changing, then over any given time interval, the number of customers served (whether they are new customers served the first time they request service, or old customers who have been attempting to obtain service over a very long time) must equal the number of new customers arriving in that time interval. It is difficult to speculate about the effect of this deficiency other than to say that for very low or for very high retrial rates, the value of will be primarily determined by X and u when the server is not saturated. Under these conditions, when v is very small, the arrival of an orbiting customer is a rare event and most of the customers seen by the server are first-trial customers. The system would look much like an M/M/l queue with no waiting room. 7 8 . Similarly, when v is very large compared to A , orbiting customers would tend to be served before the next new customer arrives, again giving something like the M/M/l queue with no waiting room and a random delay of onset of service after the customer arrives. On this basis, we would expect the lack of dependence of (3) on v(t) to be most critical at intermediate values of v(t), when compared with A and JJ . It would be quite straightforward to test this suggestion by a series of calculations, but that has not been done. When A and y are constants, equation (3) can be solved easily in closed form to get nQ(t) = (1 - p) + {Jl0(0) - (1 - p)}e~Vt . (4) In particular, if 11(0) = 1 (that is, if the server is initially idle), then this becomes n o f t ; = (1 - p) + pe'^ , ( 5 ) Figure 2.1a shows plots of nQ(t) against time for 11^ (t) computed exactly and then using equation (5), for the constant parameter case with A = 2 and y = 4. It is seen that equation (5) gives a solution here that decays to the steady state value of much too fast "with time. Figure 2.1b gives a plot of the relative error in RQ(t) computed using (5). While this plot indicates a maximum relative error of only 18% at t = 1, it should be noted that for these values of A and y , equation (5) does not represent a solution significantly more accurate numerically that simply defining n 0 « » = i n f t ; = i - p, t > o 79 . It will be recalled from section 1.4 that maximum relative errors in IT (t) in 0 the constant parameter case using Rothkopf and Oren's closure method and p = 0.5 was significantly below 1%, so equation (5) cannot seriously be considered as an acceptably accurate approximation. When A or u vary with time, equation (3) must be solved numerically. Figure 2.2 is a plot of H^t) computed exactly and II0(t) computed using (3) with u(t) = 4 throughout and A (t) given by equation (1.4.4) with A ^ = 2.0 and P = 10.0. The retrial rate is also taken as oscillatory now, using the formula The 't-1' was included in the oscillatory part to give the effect of retrial rate changes lagging hew customer arrival rate changes.1 Equation (3) does lead to a solution which reproduces the periodic behaviour of 1One might attribute this kind of variation in retrial rates to unsuccessful customers (of which there will be more when the new customer arrival rate is higher) initially increasing their rate of reapplication for service. Then, eventually the ones who still don't succeed in obtaining service will decrease their rate of retrial out of frustration or exhaustion. Of course, the retrial rate, v(t), is applied to all customers in orbit — those who have just entered the system and those who have been orbiting for some time, and so it represents a kind of average collective 'enthusiasm' of all customers in orbit at a particular time. This is not the place to discuss in greater detail the psychology of unsatisfied customers, but these remarks underline again the fact that the usefulness in actual application of all of the computational techniques discussed or used in this thesis may well be overshadowed by non-mathematical considerations. The main goal of the work in this thesis is to find viable ways of circumventing the solution of large systems of first order linear differential equations. This is only a sub-goal of the much larger goal of being able to predict the behaviour of real-life single server retrial queues. 0 1 . 1 + (6) 80 . the exact solution rather well, but the approximate solution overshoots badly at the relative extrema of ^^ (t). This overshoot causes the approximate ^(t) t 0 become slightly negative in this case in the vicinity of the relative minima of H0 (t). The mean relative error over the time interval [0,10] is 0.336 and over the time interval [10,20] is 0.347, which is also not very good. Method IB f m(t) ) A second single-variable closure approximation can be obtained from the fourth equation in Table 2, m'(t) = A - - T[Q(t)} , (7) to determine m(t). We note that at steady state, i - JIO = n2 = p and n .(c + i) P(c + 1) _ i m = --1 - P \ Solving (8) for n^= 1 - II , in terms of m, we get (8) i - n = • , <9> 0 (c + 1 + m) which, on substitution into (7) gives the equation m'(t) = A - U_m_f_t2 (10) c + 1 + m(t) As with equation (5) above for I I (t), equation (10) leads to a transient time-0 81 . dependent solution in the constant parameter case which converges much faster than the exact solution to the steady state value of m(t). This can be seen in Figure 2.3a. Figure 2.3b shows that the maximum relative error is again much greater than for Rothkopf and Oren's method applied to the M/M/l queue. The results of a calculation based on (10) when the arrival rate is oscillatory are shown in Figure 2.4. The tracking of the exact graph of m(t) by the solution of (10) is remarkable, though again the closure approximation overshoots on the bends. The average relative error on the time interval [0,10] is 0.269 and on the time interval [10,20] is 0.196. These compare with values of 0.0389 and 0.0367 respectively for the Rothkopf and Oren method applied to the simple M/M/l queue under similar loads. For what it's worth, equation (10) is separable when X , y, and v are constants, and can be rewritten as r + 1 + m X(c + 1) + (X - \x)m dm = dt (11) Integration of both sides gives y f c + i) 7 T~ 1o? X - y (X - u) X(c + 1) + (X - u)m X(c + 1) m(0) = t (12) For example, when ^ = 2, y = 4, v =2, and m(0) = 0, this reduces to •m m -y- - 2 • log (1 - —) = t (13) 82 . We could regard this as a closure method for computing approximate times at which m(t) attains particular allowed values in the constant parameter case. Since it is not a particularly difficult nonlinear equation to solve numerically, it can also be used as an alternative to equation (10) in computing values of m(t) for specific numerical values of t. In fact, formula (12) is more efficient to use when m(t) is desired for only a few isolated values of t, since the integration of equation (10) requires determining m(t) for a sequence of values of t between the initial time and the time of interest. Neither of these closure approximations would be suitable for serious study of the time-dependent behaviour of a single-server retrial queue. However, it is noteworthy that the solutions of the approximating single equations mimic the exact behaviour of the corresponding characteristics of the queue so well, even when nothing resembling a steady state is occurring. 83 . FIGURE 2.1 II0(t) Computed Using Method IA in the Constant Parameter Case. Figure 2.1a compares exact values of JI0(t) for a constant parameter calculation (A = 2.0, u = 4.0) with Ilf/t) computed using the approximation Method IA, equation (2.2.5). Both curves converge quickly to the steady state value II0= 0.5. Figure 2.1.b is a plot of the relative error in the approximate njt) shown in Figure 2.1a. For these calculations, the system was assumed initially empty, i.e., n</o) = i.o. 20 1 0.3 - I e x a c t a p p r o x i m a t i o n I I T T I I | 1 1 I I 1 I 1 I I 1 1 I I T 10 (a) 15 1 1 T I I I i i i i i i i i i i i i i i T " T r i " 10 ( b ) 15 20 84 . 1 t F I G U R E 2.2 Comparison of Il^t) Computed Approximately Using Method IA With Its Exact Values in the Oscillatory Parameter Case. The graph of the approximation is the one which swings wider at all local maxima and minima. The calculation was done with Ao= 2.0, v 0 = 2.0, u = 4.0, and a period of oscillation of 10 time units. Note the small but disconcerting excursions of the approximation curve below the horizontal axis. 85 . F I G U R E 2.3 Mean Number of Customers, m(t), in the Retrial Queue Computed Using Method IB in the Constant Parameter Case. The lower curve in Figure 2.3a is the exact m(t) when X = 2.0, u = 4.0, and v = 2.0. The upper curve is m(t) computed using Method IB, equation (2.2.10). Figure 2.3b shows the corres-ponding relative error in this approximation m(t). The closure approximation estimate of m(t) approaches the steady state value of 2.0 much too fast. In these calculations, it was assumed that the system was initially empty, i.e. m(0) = 0. 20 e 0 - III | i i i i i i i i i | 1 1 1 I 10 ( a ) 15 86 . 7 FIGURE 2.4 Mean Number of Customers, m(t), in the Retrial Queue Computed Using Method IB in the Oscillatory Parameter Case. For the curves plotted here, A 0 = 2.0, u(t) = 4.0, V 0 = 2.0, and the period of oscillation is 10 time units. Initially, the system is empty, i.e., m(0) = 0.0 . 87 . 2.3 Two-Equation Closure Methods: 11^ and m 7 should like to buy an egg please,' she said timidly. 'How do you sell them?' 'Fivepence farthing for one — two pence for two,' the sheep replied. 'Then two are cheaper than one?' Alice said in a surprised tone, taking out her purse. 'Only you must eat them both, if you buy two,' said the Sheep. [Through the Looking Glass, Lewis Carroll] Methods 2A and 2B (IIQ (Q and m(tl) The first and fourth equations in Table 2 can be used as the basis of a much more accurate closure approximation than those examined for retrial queues in the previous section. We have n ^ r t ; = - A n ( t ; + u l ^ r t ; - vm(t) + v m ^ f t ; , (1) m'(t) = A - u n j ( t ) In order to turn these equations into a system involving just Il^t) and m(t), it is necessary to eliminate K (^t) and mj(t). As before, we have the exact relation IT (0 = 1 - II (t). However, for 1 0 m (^t) it is necessary to resort to some sort of approximation. At steady state, we have P f c p + 1) , m = , (2) whereas P f c + 1) m = n 0 = i - p , = p) . o) 88 . We would like to express at steady state in terms of both m and II , since then account of the values of both quantities will be taken in computing ^(t). Just by inspection, we note that m - P ( c + 1 } . + 1 _ m C l I I * 1 . . . - • - i n - — . (4) 1 - p c + 1 c + 1 This formula reflects the type of relationship expected between changes in and corresponding changes in m and in . As m increases, we expect both of its components, m Q and m^ ,, to increase since in this situation there would be larger numbers of customers present when the server is busy as well as when the server is idle. Also, m^ , is the component of m corresponding to the server being busy. As n increases, the proportion of the time the server is busy is larger and so the component m^ is expected to reflect a larger fraction of m. A second obvious approach is simply to replace p by n everywhere in (2) (or 1 - p by 11^ ) to obtain the closure formula i ^ r d i + i) m i = - ( 5 ) which effectively expresses m^ in terms of only. Again, this formula has qualitative plausibility since it expresses the fact that, for reasons cited just above, m will tend to increase directly (though not necessarily in strict direct proportion) with increases in and inversely with increases in II . We will have a problem here if II (t) attains the value zero. However, in that 89 . case all P 0 j ( t ) , n = 0, 1, 2, must be zero, and so m 0(t) is zero and we have that m^(t) = m(t). It is unlikely that n^(t) will vanish except when the system is started with the server assumed busy (then 11^ (0) = 0). However, it is likely that this proposed closure formula will be at its worst when Il^Ct) gets closer to zero. Adopting (4) as the relationship between m^t), n_j(t) and m(t) at all times in equations (1) leads to method ' 2 A \ say, and (5) leads to method '2B'. Numerical Performance of Methods 2A and 2B A set of test calculations were done using (1) with method 2A and with method 2B in the constant parameter case with A = 2, y =4, and v = 2. Relative errors are plotted against time in Figures 2.5a and 2.5b for each of these two options. The maximum relative error in H ^ t ) in this example is around 4% for both methods. The maximum relative error in m(t) is around 10% for both methods. It appears that the 2B method is slightly better in this constant parameter case than is the 2A method, though the difference is very slight. A corresponding set of calculations was done using these two closure methods in the oscillatory parameter case with X ^ = 2 and \>Q = 2, leading to the relative error plots in Figures 2.6a and 2.6b. The oscillation period here, as usual, is ten time units with maximum arrival rates occurring at t = 2.5, 12.5, and minimum arrival rates occurring at t = 7.5, 17.5, The 2A-method gives H^(t) with relative errors oscillating between roughly +0.1 and -0.1, which is comparable 90 . to the relative errors in P0(t) computed using Rothkopf and Oren's approach (sec Figure 1.7a). The relative errors in nQ(t) computed using the 2B-method here oscillated between + 0.2 and -0.2 roughly. On the other hand, relative errors in m(t) as computed using the 2A-method vary between 0 and 0.6 (except for an initial region where the errors in m(t) arc negative1), whereas the 2B-method results in relative errors in m(t) lying in the much narrower range of -0.2 to +0.1 . Keep in mind that magnitudes of relative errors can be somewhat misleading, particularly where the quantity under study has rather small numerical magnitudes. Method 2C ( II ft) and mft^ 0 What is more interesting about the two sets of Figures, 2.5 and 2.6, is that where the 2A-method results in a positive and maximal error, the 2B-method roughly gives a negative and minimal error. This strongly suggests that an intermediate method '2C should prove even more accurate than either of these two. For any one of the exact solutions we have for particular test cases, we could construct a combination of (4) and (5) to minimize in some sense the errors in the solution produced by the hybrid method '2C\ However, here we will simply test the potential value of such a hybrid closure method by quite arbitrarily throughout, errors have always been computed as 'approximate - exact', and relative errors are this quantity divided by the exact value. 91 . adopting an equal mix of the two approximations above, and for method 2C simply writing 2 d f t ; + J n (t){cTi (t) + 1} m (t) = - <m(t) - + - - > . (6) ' c + i n (t) Constant Parameter Calculations Using Method 2C Straightforward application of (6) in the constant parameter calculations done to produce Figures 2.5a,b indicates the maximum relative error in H (t) has been reduced to something like 0.013 and in m(t) to something like 0.04. This represents a decrease in the maximum relative error by factors of three roughly. A more extended series of calculations was carried out using this method 2C. With v fixed at 2.0 and JJ at 4.0, calculations were run for A = 0.4, 1.2, 2.0, 2.8, 3.6, and 3.8 (corresponding to p = 0.1, 0.3, 0.5, 0.7, 0.9, and 0.95) in the constant parameter case. Then, with y fixed at 4.0 and A fixed at 2.0, constant parameter calculations were run for V = 10, 4, 2, 1, 0.4, and 0.2 (corresponding to c = 0.2, 0.5, 1, 2, 5, and 10 respectively). Plots of the relative errors in I I o (t) for each of these values of A appear in Figure 2.7a and for each of these values of v in Figure 2.7b. Similar plots of the relative error in m(t) are given in Figures 2.8a and 2.8b. 2The value of varying the way simple closure formulas are 'mixed' is considered at greater length in section 2.4 . 92 . There is nothing very complicated about these plots. Generally, the maximum relative errors increase in magnitude for increasing A and for decreasing v for both n (0 and m(t). Also, the usually observed decrease in relative error as t gets large takes longer to reach the vicinity of zero for those calculations resulting in larger peaks in the relative error, mainly because this decrease has farther to go and not because it is remarkably slower. The other main thing to note is that over this fairly wide range of parameter values, the relative errors do not become unreasonably large. In the worst cases shown, they peak in the vicinity of 0.10, whereas in the application of Rothkopf and Oren's method to the simple M/M/l queue, these maximal relative errors were more in the vicinity of 0.04, as seen in Figures 1.4 and 1.5 . A simple direct comparison of Figures 2.7 and 2.8 with Figures 1.4 and 1.5 is not easy, since the relative errors found here for the retrial queue calculations are quite sensitive to the retrial rate -- much more so for m(t) than for n (t). The current method 0 appears viable over a wide range of operational conditions for the retrial queue. It becomes quite difficult computationally to obtain exact solutions when A or v are very large. For example, for the calculation with v = 10, X = 2,and J J = 4, the exact calculation of the P n^(t) from which H C^O and m(t) were computed in half-integer steps over the time range t = 0 to t = 25.5 required just over 10,000 evaluations of the derivative vector and about 25 seconds of computer processor time. On the other hand, the closure method required under 2000 derivative evaluations and less than 0.4 seconds of computer processor time to obtain the approximate values of n (t) and m(t) at these points. As mentioned 9 3 . before, no rigorous attempt was made to optimize any of the programs developed for these calculations and so the actual amounts of processor time required arc not significant. But, what is significant is that at a cost of incurring relatively small errors, this closure approximation in the case cited reduced the computing time required by a factor of 60. For most queue parameter value combinations, such a dramatic reduction in computing time would not be observed, yet the reductions are still significant.3 Oscillatory Parameter Calculations Using Method 2C A corresponding series of calculations was carried out using this method in the oscillatory parameter case. Calculations were done with y(t) = 4, v 0 = 2, and XQ = 0.4, 1.2, 2.0, 2.8 and 3.6 (corresponding top o = 0.1, 0.3, 0.5, 0.7, and 0.9, respectively). Another set was done where u(t) = 4, A Q = 2, and varied over the six values 10, 4, 2, 1, 0.4, and 0.2. No calculations were done with A ^ = 3.8 in the oscillatory parameter case, since the term '1.1' in equation (1.4.4) results in an average arrival rate of approximately 4.2 over a period of oscillation. This exceeds the average service rate of 4.0 over each period, and so the system is divergent. In fact, with A 0 3.6, the time-average arrival rate is 3.96, corresponding to a time-averaged utilization of P Q = 0.99, which is a system very near saturation. Since our 3In general, exact calculations become very costly for both large values of A (near y) and large values of v . Large values of A require system truncation at large dimensions for acceptable accuracy. Large values of v lead the RKF45 algorithm to use very short internal step lengths, for reasons which are not clear, and so calculations require large numbers of derivative evaluations per time increment. 9 4 . calculations were limited to the first 2-1/2 periods, and we start with an empty system, we do not have to deal with the large mean queue lengths that would occur by the time this system reaches the equilibrium periodic state. In Figure 2.9a are plotted exact values of II (t) against time, giving an indication of the nature of the variation of this quantity with time in the oscillatory parameter case. It is seen that the influence of the initial state here is more pronounced as increases. Also, the regions of maximal value of n^ (t) become increasingly sharper peaks as p ^ increases, whereas the troughs become wider. Remember that large values of II (t) correspond to less customers in the system. The wide smooth troughs (compared to the sharp peaks) are likely due to the more uniform utilization of the server because of the increased effective retrial rate, nv(t), when the number of customers present is larger. It would be interesting to repeat these calculations using an oscillatory service rate (either in-step with the arrival rate or with a slight delay relative to the oscillating arrival rate). In Figure 2.9b the relative error in ^0(t) as computed using this hybrid closure formula is plotted against t for p ^ = 0.1, 0.3, 0.5, 0.7, and 0.9 . The higher values of P lead to wider swings in this plot of relative error. On the scale used, relative errors in nfl(t) when QQ = 0.1 show up as a periodic thickening of the t-axis only. A plot of the exact and approximate IT^ Ct) against time is given in Figure 2.9c. The closure approach does not produce as wide swings from minimum to maximum values of II (t) around t = 4, 14, 24, ... units. 95 . The plots in Figure 2.9c give the impression that the closure approach ^(t) ' s decaying to the periodic steady state more quickly than the exactly computed n Q (t). This pattern of error was also noted for the one-variable closure methods examined in the preceding section. Relative mean absolute errors are listed in Table 3a for these calculations. They too increase sharply for increasing as well as for increasing time at each value of 0Q. However, this increase in the relative mean absolute errors is again due at least partly to the decreasing size with time of the quantity to which the mean absolute errors are being compared, since, the mean absolute errors appear to decrease with time. Figures 2.10a-d deal with the computation of m(t) in the oscillatory parameter case. Figure 2.10a shows m(t) against time for = 0.1, 0.3, 0.5, 0.7, and 0.9. For p 0 = 0.1, 0.3, and 0.5, the periodic steady state seems to be established very quickly, and there is little significant difference between successive periods of variation from the beginning. A more quantitative indication of this is given by the mean values of m(t) over the first two ten-unit time intervals, listed in Tabic 3a. For p = 0.7 and 0.9, the local maxima and minima of m(t) are still increasing very markedly after two such periods, and by comparison with steady state values of m computed using p = 1.1 p^ and c = 1.1, we would expect the periodic steady state will not be established until perhaps twenty or more additional periods of ten time units have elapsed. By that time, m(t) will be oscillating 96 . around a mean value of something in excess of 200 when PQ= O.9.4 In Figure 2.10b is compared m(t) computed exactly when using this closure method with \ Q = 3.6. The approximate m(t) tracks the exact m(t) very well, giving an estimate of the exact m(t) which is slightly too large in this case over the time interval considered. At least over the first 2-1/2 periods, there is no sign of the accuracy of the approximation deteriorating. The largest gaps occur in the vicinity of the local maxima of m(t), which occur at t = 5, 15, 25, ... . When t=5, the approximation error is 1.442 (with m(5) = 15.167, this gives a relative error of 0.095), when t = 15, the approximation error is 1.515 (corresponding to a relative error of 0.065 since m(15) = 23.168), and when t = 25, the approximation error is 1.308 (corresponding to a relative error of 0.045, since m(25) = 29.145). It is seen from Table 3a that p Q = 0.9 leads to the greatest mean absolute errors in m(t) over all of the test calculations represented in Table 3a, though not the greatest relative mean absolute errors because the mean values of m(t) over one period are much greater when p Q = 0.9 than for any of the other values of p for which calculations were done here. 4We can compare the period-averaged values of n o and m given in Table 3a with the corresponding steady-state values of JI 0 and m computed using p = 1.1PO and c = 1.1 . For Tl0 these latter are 0.89, 0.67, 0.45, 0.23 and 0.01, respectively, for Po = 0.1, 0.3, 0.5, 0.7 and 0.9 . For m, the corresponding steady-state values are 0.2596, 1.034, 2.567, 7.030, and 207.9 . The period-averaged values of IT0(t) are similar to the steady state values obtained when p and c are period-averaged in those cases when p0 is small enough for our calculations to give an apparent periodic steady state. For m(t), the similarity between the early time period-averages and steady-state values computed using mean values of P and c is less obvious from the data in Table 3a. 9 7 . Figures 2.10c and 2.10d illustrate variations of relative errors in values of m(t) computed using this closure method. Plots of these relative errors for the cases pQ = 0.7 and 0.9 were done separately in Figure 2.10d to allow adoption of a larger scale and to make it easier to trace the rather complicated curves than would have been possible if all five relative error curves in these two figures had been plotted together. It is noteworthy that the local maxima in these relative error plots quite clearly decrease in value as t increases (though the magnitude of the local minima do increase in some cases shown — the ^ 0 = 1.2 and \ Q = 2.0 cases being examples). It would be interesting, though potentially very expensive computationally, to push these comparisons between exact and approximate for, say, A Q = 2.8, to much greater values of t to determine if these oscillatory relative errors tend to a mean of zero over a complete period of oscillation, with the extremes of the relative error over such a period tending to zero as well. In Figure 2.10d, there are two swings of the relative error in m(t) from relative maximum to relative minimum. The difference between these two successive pairs of extrema are 0.982 and 0.969 respectively, perhaps indicating a trend. One other feature of Figures 2.10c and 2.10d is the complexity of these relative error curves. They are by no means primarily simple sine curves. This complexity seems to repeat itself each period though some smoothing out of irregularities is noticeable on comparing the portions of the curves in the second full period to their parts in the first full period. This complexity seems to be a 98 . feature of relative error plots. As mentioned earlier, a second series of oscillatory parameter calculations was carried out, keeping A^ at 2.0, jj(t) constant at 4.0, and considering the values v Q = 10, 4, 2, 1, 0.4, and 0.2 . These calculations were done to assess the accuracy of this closure method when the mean retrial rate changed relative to new customer arrival rates and service rates. While the steady state value of II ^ does not depend upon the retrial rate for reasons already given, ^(t) does show some variation with v (t) in the oscillatory parameter case, as illustrated in Figure 2.11a. The server is more likely to be idle if orbiting customers attempt only rarely to obtain service than if they do so relatively frequently. The degree to which II0(t) is affected by varying is not as great as when A Q is varied. Il^ Ct) is much more sensitive to the rate at which new customers enter the system than to the rate at which unsuccessful customers reapply for service. Figure 2.11b shows the variation of the relative error in the approximate ^(t) computed using the current closure method for the sequence of values of considered. Figure 2.11c compares the approximate no(t) with its exact counterpart when = 10. The worst correspondence between the two occurs in the vicinity of the local minima of R (t) which tends to magnify local maxima of relative errors. The closure method gives local minima in H (t) which are somewhat shallower, broader, and shifted compared to those of the exact II (t). 99 . Corresponding results for m(t) computed for this sequence of values of VQ arc depicted in Figures 2.12a-c. Other than the observation that relative errors in the approximate m(t) for the smaller values of show very weak oscillatory behaviour as a function of time, these figures require little discussion. Relative mean absolute errors are listed in Table 3b for this series of oscillatory parameter calculations corresponding to various retrial rates. The interpretation of these numbers corresponds to that applying to Table 3a, and the discussion will not be repeated.5 The relative mean absolute errors given in Tables 3a and 3b for oscillatory parameter calculations done using this closure method can be compared with their counterparts for p (t) and m(t) in Table 1 where the Rothkopf and Oren method 0 has been applied to the simple M/M/l queue. Generally the values in Table 1 are smaller than the corresponding ones in Tables 3a and 3b, in some cases by large factors. The relative mean absolute errors for m(t) when = 0.1 are smaller in Table 1 by a factor of near 50 than those in Table 3a. However, in this instance, the relative mean absolute errors in m(t) from Table 3a are still under 2% when p^ = 0.1, and so this large difference in accuracy is really the difference between an exceedingly accurate approximation and an acceptably accurate approximation. Where relative mean absolute errors become larger, this retrial queue closure method gives results which do not exceed the numbers in Table 1 by more than a 5For the calculations represented in Table 3b, the period-averaged value of p is 0.55. It would be very interesting to determine if for sufficiently large values of t, the period-averaged value of (t) approached nQ = 1 - p = 0.45. There is some support for expecting such a trend from the data in Table 3b, though the data given apply only to the first two full periods of oscillation. 100 . factor of two, and in most cases, the correspondence between the relative mean absolute errors in the two sets of calculations is much better. Computation of Period-Averages of ^ (Q a n ( j m(t) Before leaving discussion of this two-variable closure method, we note the comparison of exact and approximate period-averaged values of nQ(t) and m(t) given in Tables 4a and 4b. Exact and approximate values of nQ(t) and m(t) were computed for integer and half-integer values of t. The period-averages were obtained from these equally time-spaced values using Simpson's rule. In a footnote in section 1.4 above, we have explained why we believe that this approach leads to no significant error in the numerical integration step, so that discrepancies between period-averages of exact quantities and period-averages of their approximations computed using this closure method are effectively totally attributable to the closure approximation. It is seen from Tables 4a and 4b that these period averages can be computed very accurately from the closure method results. The relative error in the mean of IT (t) is far below 1% for most of the test values of A and v and rises 0 0 0 to 2.5% in the first ten-unit time period when = 0.2 and 4.3% in the first such period when p = 0.9. These cases are near the computationally most difficult end of the range of values A and v may have. The relative errors in the period average of m(t) are quite uniform over the set of computations for various values of v , being in the 8% - 10% range for the first ten-unit time period, and in the 6% - 10% range for the second such period. With v fixed at 101 . 2.0, the relative errors in the period averaged m(t) vary with from a low of just under 1% when p Q = 0.1, to a 6% - 9% range when QQ = 0.5, 0.7, and 0.9 . In both Tables 4a and 4b, the relative errors in the period averages of m(t) generally are smaller in the second ten-unit time period than in the first. This occurs in part because these averages are increasing in successive periods as the system evolves from its initial empty state to its eventual periodic steady state. The absolute errors are also given in these tables and are seen to be relatively unchanged in the two periods for which numbers were obtained and in most cases slightly larger in the second period compared to the first. The same trend in relative errors is not observed for the period averages of n (t), which generally decrease slowly through successive periods. 102 . T A B L E 3a Exact Period-Mean Values of IT0(t) and m(t) and Mean Absolute and Relative Mean Absolute Errors When These Quantities are Computed Using Closure Method 2C for Various Values of p . o The closure approximation values of II0(t) and m(t) on which these numbers arc based were calculated using Method 2C, equations (2.3.1) and (2.3.6). The first column gives the value of p 0. The second gives the time interval to which the computed mean values apply. The next two columns give the exact mean values of IT 0(t) and m(t) respectively over these time intervals, the next two give the mean absolute errors in the closure approximation estimate of these quantities over the indicated time interval (computed as integrated absolute errors divided by the interval length), and the last two columns give the corresponding relative mean absolute errors (computed as the mean absolute error divided by the mean value of the quantity in question). Exact Mean Mean Absolute Relative Mean Value Error Absolute Error p o ( V (m) < V (ra) ( V (m) 0.1 [0,10] [10,20] .8930 .8900 .1509 .1571 .0005177 .0004485 .002962 .002699 .000579 .000504 .01962 .01718 0.3 [0,10] [10,20] .6840 .6701 .8979 .9789 .007470 .008760 .03718 .03943 .01092 .01307 .04141 .04028 0.5 [0,10] [10,20] .4969 .45595 2.7446 3.4260 .02868 .02946 .2226 .2257 .05771 .06461 .08112 .06588 0.7 [0,10] [10,20] .3560 .2799 5.9320 9.0347 .05283 .05231 .5610 .5591 .1484 .1869 .09457 .06189 0.9 [0,10] [10,20] .2638 .1717 10.1175 18.0434 .06480 .05733 .8653 1.0238 .2457 .3340 .08552 .05674 103 . T A B L E 3b Exact Period-Mean Values of n Q(t) and m(t) and Mean Absolute and Relative Mean Absolute Errors When These Quantities are Computed Using Closure Method 2C for Various Values of c„ or V . o o The closure approximation values of II0(t) and m(t) on which these numbers arc based were calculated using Method 2C, equations (2.3.1) and (2.3.6). The first column gives the value of v Q, the mean retrial rate.and the second column gives the value of c Q = A Q / V Q . The third column gives the time interval to which the computed mean values apply. The next two columns give the exact mean values of II0(t) and m(t) respectively over these time intervals, the next two give the mean absolute errors in the closure approximation estimate of these quantities over the indicated time interval (computed as integrated absolute errors divided by the interval length), and the last two columns give the corresponding relative mean absolute errors (computed as the mean absolute error divided by the mean value of the quantity in question). Exact Mean Mean Absolute Relative Mean Value Error Absolute Error V 0 c 0 (m) (m) ( V (m) 10.0 0.2 [0,10] .4734 1.9893 -.002197 .1648 .08629 .1427 [10,20] .4505 2.1655 -.0004056 .1346 .09539 .1285 4.0 0.5 [0,10] .4821 2.3078 -.002023 .1852 .06878 .1075 [10,20] .4519 2.6483 -.0008732 .1459 .07753 .09209 2.0 1.0 [0,10] .4969 2.7446 .00006752 .2224 .05771 .08112 [10,20] .4559 3.4260 -.001340 .1992 .06461 .06588 1.0 2.0 [0,10] .5230 3.3794 .005831 .2816 .05080 .08332 [10,20] .4676 4.7900 -.001536 .3062 .05807 .06392 0.4 5.0 [0,10] .5714 4.3977 .008987 .4167 .04722 .09476 [10,20] .5032 7.5169 .0002891 .5898 .05478 .07846 0.2 10.0 [0,10] .6088 5.1257 .01527 .5568 .04565 .1086 [10,20] .5434 9.8624 .005137 .9539 .05165 .09672 104 . T A B L E 4 a Period Averages of II0(t) and m(t) for Various New Customer Arrival Rates, Based on Method 2C. This tabic gives the exact period-averages of (t) and m(t) (third column), the corresponding averages of the same quantities calculated using t he two-equation closure approximation 'C* (fourth column), and the actual errors (fifth column) and errors relative to the exact values of these approximate values (sixth column) for various values of p 0. For these calculations, v was fixed at 2.0, and u(t) was constant in time with value 4.0 . The s e c o n d column indicates the time interval over which the averages were computed, and represent complete periods of oscillation in the new customer arrival rate and retrial rate. ( t ) : Exact Approximate Absolute Error Relative Err< p 0 1 Mean Mean in Mean in Mean 0. [0,10] .8930 .8930 .00001828 .00002047 [10,20] . 8900 .8900 .00000025 -.000000281 0. 3 [0 , 1 0 ] . 6840 . 6842 .0001832 .0002678 [10 , 2 0 ] . 6701 . 6713 .001154 .001722 0. 5 [0,10] .4969 . 4970 .00006752 .0001359 [10 , 2 0 ] . 45595 .45461 - .001340 -.002938 0. 7 [0 , 1 0 ] . 3560 .3596 .003644 .01024 [10,20] .2799 .2757 .004167 -.01489 0. 9 [0 , 1 0 ] . 2638 .2752 .01144 .043355 [10,20] . 1717 .1708 .0008418 -.004904 m ( t) : 0. 1 [0 , 1 0 ] . 1509 . 1521 .001120 .007419 [10,20] . 1571 . 1585 .001488 .009475 0. 3 [0,10] . 8979 .9284 . 03057 . 03405 [10,20] .9789 1 . 0101 .03118 .03185 0. 5 [0,10] 2 .7446 2 .9671 .2224 .08104 [10,20] 3 .4260 3 . 6252 . 1992 .05813 0. 7 [0 , 1 0 ] 5 .9320 6 .4930 .5572 . 09457 [10 , 2 0 ] 9 .0347 9 .59196 . 5610 .06168 0. 9 [0,10] 10 . 1175 10 .9828 .8653 .08552 [10,20] 18 .0434 19 .0672 1.0238 .05674 105 . T A B L E 4b Period Averages of II0(t) and m(t) for Various Retrial Rates for the Single-Server Retrial Queue, Based on Method 2C. This tabic gives the exact period-averages of II0(t) and m(t) (fourth column), the corresponding averages of the same quantities calculated using the two-equation closure approximation 'C (fifth column), and the actual errors (sixth column) and errors relative to the exact values of these approximate values (seventh column) for various values of v 0 , the mean retrial rate for orbiting customers. For these calculations, X 0 was fixed at 2.0, and u(t) was constant in time with value 4.0 . The third column indicates the time interval over which the averages were computed, which comprise complete periods of oscillation in the new customer arrival rate and retrial rate. The second column gives c 0 = X 0 / v 0 . n o U ) _ : 10. o 4 . o 2 . 0 1.0 0.4 0.2 m(t) : V 0 10. 0 4 . 0 2 . 0 1.0 0.4 0.2 Exact Approximate Absolute Error Relative Error 0 0.2 0.5 1.0 2 . 0 5.0 10.0 0 0.2 0.5 1.0 2 . 0 5.0 10. 0 {0,101 [10,20] r o , i o i [10,20] [0,101 [10,20J [0,101 [10,20] [0,101 [10,20] [0,101 [10,20] [0,101 [10,20] [0,101 [10,20] [0,101 [10,20] [0,101 [10,20] r o , i o i [10,20] [0,101 [10,20] Mean Mean in Mean in Mean . 4734 .4505 .4711 .4501 -.002197 -.0004056 -.004640 -.0009004 . 4821 .4519 . 4801 .4510 -.002023 -.0008732 -.004195 -.001932 . 4969 .4559 . 4970 .4546 .0006753 -.001340 .0001359 -.002915 . 5230 .4676 . 5262 .4661 . 005831 -.001536 .01121 -.003284 .5714 .5032 . 5804 . 5035 .008987 .0002891 .01573 .0005746 . 6088 . 5434 . 6241 . 5485 .01527 .005137 .02508 .009453 1. 9893 2 . 1541 . 1648 . 08285 2 . 1655 2 ..3001 . 1346 .06218 2 . 3078 2 . 4930 . 1852 . 08025 2 . 6483 2 .7942 . 1459 .05508 2 . 7446 2 .9671 . 2224 .08104 3 . 4260 3 . 6252 . 1992 .05813 3 . 3794 3 . 6609 . 2816 .08332 4 . 7900 5 . 0962 . 3062 .06392 4 . 3977 4 .8145 .4167 .09476 7. 5169 8 . 1067 . 5898 . 07846 5. 1257 5 . 6826 . 5568 . 1086 9. 8624 10 .8164 .9539 .09672 106 . -0.12 I i i i i i i i i i | i i i | i i i i i i i i i | i i i i i i i i i 0 5 10 15 20 (a) Method 'A' 0.1 -0.01 — i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 0 5 10 15 (b) Method 'B' FIGURE 2.5 Relative Errors in no(t) and m(t) When Computed Using a Two-Equation Closure Approximation With Basic Methods 'A' and 'B' in the Constant Parameter Case. Figure 2.5a results from a calculation using equation (2.3.4) to eliminate m± (t) (the so-called method 'A'), whereas Figure 2.5b results from a calculation using equation (2.3.5) (the so-called method 'B'). What is noteworthy here is not so much the magnitudes of the relative errors as the fact that the two figures are very nearly reflections of each other in the t-axis. This observation suggests the hybrid method 'C described in section 2.3 . Calculations here correspond to X = 2.0, y = 4.0, and v = 2.0 . 107 . t o "0.4 | i i 11 i i i 11 i i 11 i i i 11 i i i 11 i i 11 i i i 11 | i i , i i i 11 i | 0 5 10 15 20 25 (b) Method 'B' FIGURE 2.6 Relative Errors in n o(t) and m(t) When Computed Using a Two-Equation Closure Approximation With Basic Methods 'A' and 'B' in the Oscillatory Parameter Case. These two figures are the oscillatory parameter counterparts of Figures 2.5a and 2.5b, corresponding to A Q = 2.0, v 0= 2.0 and u(t) = 4.0 . Because these relative error curves oscillate, it is not so easy to see at a glance that those in Figure 2.6a go through their maximal and minimal values in the regions where those in Figure 2.6b go through their minimal and maximal values respectively. However, comparison of the t-coordinates of these local extrcma indicates that this correspondence in their locations is virtually exact. The shapes of the corresponding curves here are not as nearly mirror images as in Figures 2.5a and 2.5b. However, the overall range in values of relative errors observed for each quantity is not much different between Figure 2.6a and Figure 2.6b if one ignores the first five or so units of time. It is also noteworthy that the maximum magnitudes of relative errors in this oscillatory parameter case are two to six times those observed for the corresponding constant parameter cases. 10 8 . -0.02 11' 11111 1111111111111111111111[ 1111111111 | 0 5 10 15 20 25 (a) X ( o r p) v a r y i n g 0.028 -U.CD5 | i i i i i i i i i 11 i i I 11 i i 11 i i i 11 i i I 11 I i 11 i i i 11 | i i i i i i i i i | 0 5 10 15 20 25 (b) V ( o r c) v a r y i n g F I G U R E 2.7 Relative Errors in n o ( t ) Using the Two-Equation Closure Method ' C in the Constant Parameter Case With Various Arrival and Retrial Rates. Calculations represented in Figure 2.7a were done with u (t) = 4.0 and v(t) = 2.0 in all cases, with X (t) taking on the values 0.4, 1.2, 2.0, 2.8, 3.6, and 3.8 in sequence. Calculations represented in Figure 2.7b were done with u(t) = 4.0 and X(t) = 2.0 in all cases, with v (t) taking on the values 10.0, 4.0, 2.0, 1.0, 0.4, and 0.2 in sequence. The p = 0.5 calculation in Figure 2.7a is the same as the c = 1 calculation in Figure 2.7b, corresponding to the calculations shown in Figure 2.5 for n o(t). Note that the maximum relative error here is smaller by a factor of three to four for Hd(t) in that particular case. Appendix A gives details of calculations done to obtain exact values required in determining these relative errors. 109 . (a) X ( o r p) v a r y i n g 0.1 (b) v ( o r c) v a r y i n g FIGURE 2.8 Relative Errors in m(t) Using the Two-Equation Closure Method 'C in the Constant Parameter Case With Various Arrival and Retrial Rates. Calculations represented in Figure 2.8a were done with u (t) = 4.0 and v(t) = 2.0 in all cases, with X (t) taking on the values 0.4, 1.2, 2.0, 2.8, 3.6, and 3.8 in sequence. Calculations represented in Figure 2.8b were done with u (t) = 4.0 and *A(t) = 2.0 in all cases, with V(t) taking on the values 10.0, 4.0, 2.0, 1.0, 0.4, and 0.2 in sequence. The p = 0.5 calculation in Figure 2.8a is the same as the c = 1 calculation in Figure 2.8b, corresponding to the calculations shown in Figure 2.5 for m(t). Note that the maximum relative error here is smaller by a factor of 2.5 to 3 for m(t) in that particular case. Appendix A gives details of calculations done to obtain exact values required in determining these relative errors. 10 T -15 -Tr-io T 25 (a) X ( o r p ) v a r y i n g O O -as 10 ( b) A ( o r 15 V ' " I «• W 25 v a r y i n g 1 at a 07 a« -05 -,04 OJ • 02 • ai e x a c t . FIGURE 2.9 n o(() Computed Using the Two-Equation Closure Method ' C in the Oscillatory Parameter Case for Various Values of p Q . Figure 2.9a simply displays the exact n0(t) for various values of p0 . Figure 2.9-b displays the corresponding relative error in n0(t) calculated using closure method 2C. Figure 2.9c compares the exact and approximate U0 (t) when p 0 = 0.9. Note that on the scale employed in Figure 2.9b, the relative error curve for p0 = 0.1 is indistinguishable from the t-axis. Appendix A gives details of calculations done to obtain exact values used in preparing these plots. I l l . FIGURE 2.10 m(t) Computed Using the Two-Equation Closure Method ' C in the Oscillatory Parameter Case for Various Values of p . 0 Figure 2.10a simply displays the exact m(t) for various values of pQ. Figure 2.10b compares the closure approximation m(t) with its exact counterpart when p^ = 0.9. The plots of relative errors in the closure approximation m(t) are displayed in Figures 2.10c and 2.10d. Appendix A gives details of calculations done to obtain exact values required in determining these relative errors. 112 . (a) V ( o r c ) v a r y i n g 0 o I 09 0.3 0.7 0.6 4-> 04 OJ 02 01 0 5 10 15 20 25 (c) V =10 ( C =.2) 0 O F I G U R E 2 . 1 1 II 0 (t) Computed Using the Two-Equation Closure Method ' C in the Oscillatory Parameter Case for Various Values of c Q . F i g u r e 2.11a shows the exact II 0(t) for var ious values of c 0 . A t the local m i n i m a , c o = 0.2 corresponds to the lowest curve and c 0 = 10 to the highest. F i g u r e 2.1 l b shows cor respond ing re lat ive errors in the closure approx imat ion n o ( t ) , with the widest extremes cor respond ing to c 0 = 0.2, and the narrowest to c 0 = 10 . In F igure 2.11c is compared the exact and approximate n Q(t) in the worst case c o = 0 . 2 . (b) V ( o r c ) v a r y i n g O o 113 . 0 5 ( 0 15 20 25 ( c ) C =10 (V =10) O 0 FIGURE 2.12 m(t) Computed Using the Two-Equation Closure Method ' C in the Oscillatory Parameter Case for Various Values of cQ. Figure 2.12a shows the exact m(t) for various values of c Q. Throughout, c Q = 0.2 corresponds to the lowest curve and c Q = 10 to the highest. Figure 2.12b shows corresponding relative errors in the closure approximation m(t), with the curve having the lowest local minima corresponding to c Q = 0.2, and the one with the highest local minima corresponding to c = 10 . In Figure 2.12c is compared the exact and approximate m(t) when c Q = fo.O . 114 . 2.4 Two-Equation Closure Methods: . m. and v Then he got up on top with a tip of his hat. 'I call this game FUN-IN-A-BOX', said the cat. 'In this box are two things I will show to you now. You will like these two things', said the cat with a bow. [The Cat in the Hat. Dr. Seuss] Method 2D ( II m. v) Because the steady-state probabilities for the single-server retrial queue form the negative binomial distribution, it is possible to write down an analogue of the Rothkopf and Oren method using the following equations from Table 2: m' (t) = A - \i{l - II (t) }, v' (t) = A + U - \lT[Q(t) {2m( t) + 1} + 2]Sn>o(t) . As before, m(t) is the mean number of customers in the system at time t, m0(t) is the contribution to m(t) by server idle terms, v(t) is the variance of the number of customers in the system, and n^Ct) is the probability that there arc no customers in the system at time t. At steady state, p r e + i) (i) n o = 1 - p, ' mQ = cp, m ' 1 - P and (2) P(c + 1) v = _(1 - P) The object here is to eliminate 11^ (0 and m0(t) in equations (1) to obtain a system of two equations which can be used to compute values of m(t) and v(t), at least approximately. 115 . Following Rothkopf and Oren, we note that 2 m / ( v-m) p 0 0 - u - v"1 - l-s-) Thus, l/(c+l) ( ^ m /{(v-m)(c+1)} \ = a - P) = P q o = l ^ r ) . (4) When v - m = 0, we follow Rothkopf and Oren in writing n -m/(c+l) n Q = e (5) In the simple M/M/l case, the last term in the equation above for v'(t) does not appear. Examination of the formulas in (2) easily yield a number of relationships for m at steady state. Among these are mQ - cp = c (1 - JlQ) , (6) obtained by direct replacement of p by 1 - II^ mo = m]lo - n i ' ( 7 ) obtained by replacing cp by mQ in the numerator of the formula for m, solving for m^ , and replacing- p by or 1 - p by II ; and finally, mo = vJ[o2 ~ n i ' ( 8 ) obtained as is (7) but starting with the formula for v. 116 . Formula (7) has the most intuitive appeal, since it relates m 0(0 in a direct way to m(t) and ^(t). We would expect m0(t) to increase with both m(t) and with n^CO in general. Unfortunately, substitution of (7) into the right-hand side of the equation for v'(t) causes it to become identical to the right-hand side of the equation for m'(t) in (1), which is unrealistic. The subtracted terms in (7) and (8) are also cause for concern, and computational tests quickly indicate that using (8), v'(t) is prone to become negative when v is near zero (as is the usual case at the beginning of a calculation when the system is in a known definite state). This leads to v(t) becoming negative, and so (m/v) is negative, and the right hand side of (3) is not real. As a result, we will initially consider using only formula (6) in conjunction with (5) and (1) to comprise a closure method for computing approximations to m(t) and v(t). We will call this approach Method 2D. Numerical Performance of Method 2D A series of calculations was carried out using Method 2D. The results were not very encouraging. In constant parameter calculations, relative error maxima were several times larger than they were for corresponding M/M/l queue calculations using the Rothkopf and Oren method. This is best seen by comparing Figures 2.13a, 2.13b, and 2.13c directly with Figures 1.4, 1.5, and 1.6. Figure 2.13 shows the relative errors in IT m, and v, as functions of time for a sequence of new customer arrival rates. 117 . For the calculations represented in Figure 2.14, the new customer arrival rate is held constant and the retrial rate is varied. When the retrial rate is comparable or smaller than the new customer arrival rate, relative accuracy over the time range shown in the figures is poor. Of particular interest are Figures 2.13c and 2.14c, dealing with the variance, since no other closure technique has been discovered for computing this quantity for the retrial queue so far. Figure 2.15a shows the exact v(t) and v(t) obtained using this proposed closure method when p = 0.5 and c = 10 (actually A = 2.0, u = 4.0, and v = 0.2). This figure gives an alternative impression of the approximation errors arising. The situation is much the same for oscillatory parameter calculations using this approach. Figures 2.15b, 2.15c, and 2.15d contain plots of the exact v(t) and the approximate v(t) for the most difficult cases studied (p^ = 0.7, c^ = 1.0, in Figure 2.15b; P = 0.9, c = 1.0, in Figure 2.15c; and p = 0.5, c = 10.0 in O 0 0 0 Figure 2.15d). The approximation errors in n^t) and m(t) computed using this method are not as relatively large as for v(t). However, they are significantly greater than when the best of the methods developed in the preceding section arc used. Figures 2.16a and 2.16b show the approximate n^(t) along with its exact values in the oscillatory parameter case when A = 3.8, u = 4.0, v = 2.0, and 0 0 when = 2.0, u = 4.0, and = 0.2, respectively. The corresponding comparisons of approximate m(t) with exact m(t) in these two worst-of-the-cascs-J i g . studied situations are displayed in Figures 2.16c and 2.16d. Figures 2.16a and 2.16c go with Figure 2.15c, and Figures 2.16b and 2.16d go with Figure 2.15d. In Table 5 are given period averages of n0(t), m(t), and v(t) as well as relative mean absolute errors. These numbers indicate further the breakdown of the accuracy of this approach as the problematic large or large c^ situation is approached. Error measures in Table 5 should be compared with corresponding quantities in Tables 4 and 1. Method 2E ( II . m. v) 1—0 — There are two approximations in the method just described. The first is to use (4) to express njt) in terms of m(t) and v(t). This rather complicated formula is unlikely to be a significant source of error since its analogue in the Rothkopf and Oren method for simple M/M/l queues worked extremely well. Also, relative errors in H (0 computed using this method are much smaller generally than for m(t) and v(t). The second approximation in the computations just described was in using equation (6) to eliminate m 0(0 from the formula for v'(t). Equation (6) was obtained simply by substituting 1 - 11^ (0 for p in the steady-state formula for m ,^ and so is analogous to formula (2.3.5) in the previous section. For the closure method developed in section 2.3, there was a second expression, equation (2.3.4), which expressed the component rn^(t) in terms of m(t). The analogue of equation (2.3.4) in the present context is equation (7), which, as already 119 . explained, cannot be used by itself to eliminate nip(t) in equations (1). However, in section 2.3, it was found that a superior approximation for m (^t) was obtained by combining formulas (2.3.4) and (2.3.5). And, in fact, the difficulty in using equation (7) here does not arise if a combination of (6) and (7) is used: mQ(t) = r\c(i - J[Q) + a - n ; ( m n o - n J } , (9) 0 < T) < 1 This formula, (9), is the analogue of equation (2.3.6) in the preceding section. There, the value n = 0.5 was adopted somewhat arbitrarily. Here we will look a little more closely at the effect of varying the mixing factor n because of the large errors which can arise in the approximate v(t) when n = 1. Of course, one could stretch the reader's indulgence further by starting with a mixture of formulas (6), (7), and (8), but the author's advancing age precludes an adequate study of such a scheme. In addition, formula (8) would not appear to have a strong intuitive appeal, and as well, contains v(t), which is a quantity we are not able to estimate very accurately at present. Constant Parameter Calculations Using Method 2E The calculations already described in this section correspond to i] =1. The highest p or p and the highest c or c f l (corresponding to smallest v or ) calculations and the intermediate p (or ) = 0.5, c (or c^ ) = 1.0 calculations were repeated with n = 0.75, 0.5, 0.25, 0.1, and 0.05. Since the 120 . formula (7) is exact when ^(t) = 1.0, a scheme was also tested in which the mixing factor n varied with time by computing n(t) = 1 - ^ (t). Results by Method 2E in constant parameter calculations were very encouraging. Figure 2.17 illustrates the relative errors in v(t) for the various mixing factor values. In the three situations illustrated, 1 = 0.25 seems to give the lowest overall errors of all of the values of n tested. When A = 2 and v = 2, the smallest of the local maxima of absolute relative errors observed is 0.06, occurring in the calculation with n = 0.25. This is less that one-quarter of the maximum relative error occurring when n = 1 is used. Actually, for this particular calculation, the relative error starts off with a value of 0.102 when t = 0.5, but its value drops rapidly toward zero. For low values of t, v(t) also has a smaller value and so larger relative errors are not as significant. When A = 3.8 and v = 2, the n = 1 calculation gives a maximum absolute relative error in v(t) of approximately 0.48, whereas with n = 0.25, this maximum is likely to be around 0.121, again smaller by a factor of near four. When A = 2.0 and v = 0.2, the maximum relative error in v(t) decreases monotonically with the values of n used, and has a value of 0.52 when n = 0.05, the smallest value of H for which calculations were done. This compares with a xThe absolute relative error in v(t) when n = 0.25 occurs at a time greater than t = 25.5, the largest value of t for which numbers were obtained in these calculations. However, as seen in Figure 2.17b, the relative error curve is very nearly horizontal by t = 25.5, and which point the relative error is -0.1138. It is reasonable to infer that the actual turning point of the curve will not be far below -0.12.-121 , maximum of 0.81 when 0= 1 . Thus, as v becomes small relative to A , use of the formula (9) with n < 1 does lead to more accurate estimates of v(t), though the improvement is not as great as for the other two examples illustrated in Figure 2.17. Only in the intermediate case ( A = 2.0, V = 2.0) is the variable n. approach competitive with fixing n at some constant value. In the other two cases studied, it is actually no better than a poor intermediate between the worst case H = 1 calculation and the calculation using the absolute relative error minimizing constant value of n. . Jumping ahead a bit, we note here that this conclusion regarding the accuracy of the variable n approach has also been born out by comparing relative mean absolute errors given in Table 6 for corresponding oscillatory parameter calculations. Figure 2.18 gives a more graphic impression of the limitations of formula (9) when v is small relative to A in the constant parameter case. For the calculations represented in Figures 2.17a and 2.17b, A = v and A - 2 respectively, whereas in Figure 2.17c, A = 10v. If, for the first two cases, curves of one approximate quantity are plotted on the same grid as the exact quantity for different values of 1, the gap between the n = 1 approximation and the exact curve is nearly evenly filled by curves corresponding to a sequence of values of n between 1 and the best value of n for the given queue parameter values. Figure 2.18d illustrates this by plotting v(t) when A = 3.8 and v = 2.0, and the approximations to v(t) computed with n.. = 1, 0.75, 0.5, and 0.25. 122 . Figures 2.18a, 2.18b, and 2.18c give corresponding plots for n (t), m(t), and v(t), respectively, when A = 2.0 and v = 0.2. In each case, it is seen that the approximate curves for the various values of 1 form a tight band well separated from the exact curve. This means that the comparative lack of improvement observed in using formula (9) when V [s quite small relative to A is due to the intrinsic limitations of formula (9) and not the result of an inadequate search for the best value of n . For constant parameter calculations with A not large compared to v , the formula (9) gives an closure approximation of comparable accuracy to those of the previous section when a value of n somewhere between 0.5 and 0.25 is used. When the retrial rate, V , is low relative to the new customer arrival rate A , we are still without a closure method for computing v(t) with truly acceptable accuracy, though formula (9) represents best (and significant) improvement over (6) if a value of n near 0.1 is used. Oscillatory Parameter Calculations Using Method 2E The oscillatory parameter case is, as usual, more complicated to analyze here. Figure 2.19a shows the exact v(t) plotted on the same grid as the approximate v(t) curves obtained using n = 1. 0.75, 0.5 and 0.25. When A Q = VQ = 2.0, this figure indicates that a highly accurate calculation of v(t) should result if a value of H approximately midway between 0.5 and 0.25 is used. Figure 2.19b shows a similar comparative plot when A ^ = 3.6 and = 2.0. Again, it appears that the best value of H is somewhere between 0.5 and 0.25. 12 3 . Figure 2.19c gives a similar comparative plot for v(t) in the oscillatory parameter case when A = 2.0 and v = 0.2. Here, the n = 0.75 curve has been 0 0 left out and n = 0.1 and 0.05 curves included. Unlike the corresponding constant parameter calculations when the new customer arrival rate is much larger than the retrial rate, here the use of formula (9) with a value of n in the neighborhood of 0.1 results in remarkable improvement in accuracy in the computation of v(t). The improvement on using (9) is not as dramatic for n (^t) and m(t) in this oscillatory parameter case, though still significant. Without belabouring the discussion, we note that a closure method based on (9) appears to give acceptably accurate approximations for v(t) in this oscillatory parameter case using similar values of n. to those found best for corresponding constant parameter calculations. This tends to allay suspicion that formula (9) is simply a numerical trick, and to strengthen the presumption that it actually approximates the relationship between ^ (^t) and II (^t) and m(t) in some generality. Computation of Period-Averages of m(t) and v(t) In Table 6 are found period-averages and mean absolute relative errors of Ilj(t), m(t) and v(t) for various values of n. in the three cases of oscillatory parameter calculations done. These entries can be compared with the entries for n. = 1 in each case to get a numerical measure of the improvements in accuracy formula (9) represents over formula (6). 124 . The entries for v(t) in Table 6 are of greatest interest, and the results are quite remarkable. When the parameters have intermediate values, \ Q = = 2.0 with u(t) = 4.0 as usual, the smallest of the relative errors observed in the mean value of v(t) and in the mean absolute errors occurred for n = 0.35.2 When n = 0.35 was used, the relative error in the computed mean value of v(t) over the time intervals [0,10] and [10,20] were just under 5% and 6% respectively, compared to just over 33% and 18% respectively when n = 1 was used. This compares well with the values both near 16% found for the corresponding calculation using Rothkopf and Oren's method on the M/M/ l queue (see the P Q = 0.5 entries in Table 1). Thus the accuracy of the computed mean values of v(t) over these two periods of oscillation improves by factors of more than six and more than three respectively. Similarly, the relative mean absolute error in v(t) over these two time intervals drops from 40% and 28% to 8% and 8%, respectively. These compare with 15% and 16% as the relative mean absolute errors in v(t) observed using Rothkopf and Oren's method on the simple M/M/l queue in the corresponding pQ = 0.5 case (see Table 1). These are the first instances observed in which the relative accuracy of a closure technique for the retrial queue is better than that of the Rothkopf and Oren method for the corresponding M/M/l queue, and the increase in accuracy here is significant. 2The "mean values" of TT^ Ct), m(t), and v(t) or their absolute errors referred to here are always computed over 10-unit time intervals, the period of oscillation of A(t) and v(t) used throughout in oscillatory parameter calculations. For discussion purposes throughout this thesis, such mean values are always computed over two "standard" time intervals, [0,10] and [10,20]. 125 . The corresponding improvements in accuracy shown when A Q = 3.6 and = 2.0 are even more dramatic. When 1 = 1 , the relative errors in the computed mean values of v(t) over the first two complete periods of oscillation are 78% and 55% respectively, whereas when P = 0.35, these relative errors drop to 8% and 0.5% respectively, an improvement by factors in the neighborhood of 10 or better. The relative mean absolute errors in v(t) in this case drop from 78% and 55% when n = 1 (the fact that both types of relative errors discussed here are equal indicates that the approximate v(t) is always greater than the exact v(t) ) to 8% and 2% when n = 0.35 . These numbers compare with 20% and 25% respectively for the Rothkopf and Oren method applied to the simple M/M/l queue. Because the relative error in the computed period-mean of v(t) is positive when n = 1 and negative when n = 0.1, and is apparently continuous,3 one could in principle find values of n which would make this period-mean of v(t) exact in any pre-specified period. Thus, it is noteworthy that errors in the mean value of v(t) in the two disjoint periods for which calculations were done become small for the same or nearly the same values of n_. But more significant is the fact that the mean absolute error in v(t) decreases so dramatically here, since in computing this quantity, there is not the potential for large positive error components in one 3The claim of "apparent continuity" is based on the numbers in Table 6 for various values of n between 1 and 0.1 . This looks like the sort of assertion that could be proven rigorously without a great difficulty by someone more familiar with the theory of ordinary differential equations, n is effectively a simple constant numerical parameter in equations (1), and the claim would be that the solution-of the system (1) under these circumstances is continuous in n • 126 . region to cancel large negative error components in another region as can occur in computing errors in the period-mean itself. The numbers in Table 6 for the calculation of v(t) when A^ = 2.0 and v 0 = 0.2 are not as good as for the two cases just described, though use of equation (9) with an appropriately chosen value of n still does give significant improvements in the accuracy of computed approximate mean values compared to calculations when n = 1. When n = 1, the relative error in computed mean values of v(t) are 63% and 40% respectively for the two standard periods, when n = 0.05, these are -4% and -18%, respectively. One would expect that at some value of 0 between 0.1 and 0.05, both relative errors would be found roughly between 10% and -10%, which is an improvement by factors of 4 - 6 compared to the r) = 1 case. The same is true for the relative mean absolute errors in this case. In all three of the cases covered in Table 6 for v(t), the variable n method gives, if not the poorest relative accuracy, then close to the poorest as far as the computation of these period-mean values is concerned. The dependence of the relative error in the computed period-mean values of II0(t) and m(t), and their relative mean absolute errors on n is not as great as for v(t). Because these relative errors are generally much smaller than for v(t), we simply note here that inspection of Table 6 indicates that while minimal relative errors for ^(t) and m(t) seem to occur when n is closer to 1 rather than 0.35 when retrial rates are not much smaller than new customer arrival rates, their values are generally well under 10% over most intermediate values of 127 . H and compare well with the numbers observed for PQ(t) and m(t) of the M/M/l queue with Rothkopf and Oren's method. Estimation of mJt) Using Method 2E It appears that in these cases where the retrial rate is not much smaller than the new customer arrival rate, equation (9) is not only compensating for deficiencies in equation (6), but also to some extent for deficiencies in equation (4), since with optimal values of n, we are observing significantly smaller relative errors here than for the corresponding results for v(t) for the M/M/l queue with . Rothkopf and Oren's original method. Evidence for this claim is found in the fourth part of Table 6 dealing with mj[t) and in Figure 2.20 . If it was simply a matter of equation (9) giving more accurate estimates of m 0(0. then we would expect that the value of H leading to the most accurate v(t) should also give much improved estimates of m 0(t). Not only is this not so according to Table 6, but the relative errors found for rn. (t) in virtually all cases are extremely large, indicating that equation (9) generally would be a completely unacceptable way to compute m 0(0 even though it is an integral part of a method which produces quite acceptable estimates for n^(t), m(t), and v(t). The numbers in Table 6 for m (t), and even more so, the plots in Figure 2.20, 0 are disturbing because they cast at least some doubt on the entire rationale of closure approaches as described by Rothkopf and Oren in their paper, and in this 12 8 . thesis. Of the several closure methods examined so far, this is the first which required two distinct closure formulas, and so it is the first instance where accuracy in some of the quantities computed may be enhanced by taking advantage of the interplay between closure formulas used. What is also remarkable after seeing the plots in Figure 2.20 is that the n = 1 calculations give as accurate values for n (t), m(t), and v(t) as they do, since the estimates of ni (t) when n = 1 are almost as bad as they can be without being totally absurd numerically. However, it is not valid to conclude that closure formulas are nothing but fancy or complicated "fudge-formulas" , since here, ^ ( t ) is also estimated using such an imposed formula, and the resulting estimates reflect exact values very well. A more viable approach for computing m (t) is described in the next section. Also, the method described in the preceding section produces an estimate of m^(t) as well as m(t), and the difference of these will give an estimate of m 0 ( t ) . No attempt has been made to determine the accuracy of such estimates of m 0 ( t ) however. If one was going to do that, it would probably be worthwhile to first explore the improvements possible by varying the mixing factor in equation (2.3.6). Since such an approach to computing m 0 ( 0 would involve taking the difference of two approximate values, it would be desirable to maximize the accuracy of both approximate values first. 129 . T A B L E 5a Period Averages and Mean Absolute Errors in II (t), m(t) and v(t) for Various New Customer Arrival Rates Using Method 2E. The entries in this table correspond to oscillatory parameter calculations using Method 2E with v 0 = 2.0 and u(t) = 4.0, a constant, in all cases. The first column below gives the value of p0 = A 0 /u , the second gives the time interval over which the mean values were computed, the third column gives the exact mean, the fourth column gives the closure method mean value, the fifth column gives the error in these period-mean values, and the sixth column gives the relative error in the closure method mean values. The last two columns give, respectively, the mean absolute error over the indicated time intervals, and the corresponding relative mean absolute error. 0.1 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 8 9 3 0 . 8 9 0 0 . 8 9 3 5 . 8 9 0 0 . 0 0 4 9 8 7 . 0 0 0 0 0 1 . 0 0 0 5 5 8 5 . 0 0 0 0 0 1 1 . 0 0 2 4 5 2 . 0 0 1 8 4 2 . 0 0 2 7 4 6 . 0 0 2 0 6 9 0.3 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 6 8 4 0 . 6 7 0 1 . 6 8 5 7 . 6 7 0 0 . 0 0 1 7 4 3 . 0 0 0 1 2 0 9 . 0 0 2 5 4 8 . 0 0 0 1 8 0 4 . 0 2 1 5 9 . 0 1 8 3 4 . 0 3 1 5 7 . 0 2 7 3 6 0.5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 6 9 . 4 5 6 0 . 4 9 1 5 . 4 5 2 3 . 0 0 5 4 1 4 . 0 0 3 6 0 8 . 0 0 1 0 8 9 . 0 0 7 9 1 2 . 0 5 2 8 6 . 0 4 8 2 4 . 1 0 6 4 . 1 0 5 8 0.7 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 3 5 6 0 . 2 7 9 9 . 3 4 5 3 . 2 6 3 1 . 0 1 0 7 1 . 0 1 6 8 0 . 0 3 0 0 8 . 0 6 0 0 0 . 0 8 5 3 5 . 0 7 8 8 0 . 2 3 9 8 . 2 8 1 5 0.9 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 2 6 3 8 . 1 7 1 7 . 2 5 3 0 . 1 4 5 2 . 0 1 0 8 4 . 0 2 6 4 9 . 0 4 1 0 8 . 1 5 4 3 . 0 9 8 2 2 . 0 7 5 2 1 . 3 7 2 3 . 4 3 8 1 m(t) : 0.1 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 1 5 0 9 . 1 5 7 1 .1545 .1604 . 0 0 3 5 3 6 . 0 0 3 3 7 1 . 0 2 3 4 2 . 0 2 1 4 6 . 0 1 1 3 1 . 0 1 1 2 1 . 0 7 4 9 1 . 0 7 1 3 5 0.3 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 8 9 7 9 . 9 7 8 9 . 9 2 3 7 . 9 9 3 4 . 0 2 5 8 0 . 0 1 4 4 2 . 0 2 8 7 4 . 0 1 4 7 3 . 1 2 6 8 . 1 2 1 7 . 1 4 1 2 . 1 2 4 3 0.5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 2 . 7 4 4 6 3 . 4 2 6 0 2 . 8 2 0 0 3 . 2 6 1 1 . 0 7 5 3 8 . 1 6 5 0 . 0 2 7 4 6 . 0 4 8 1 5 . 3 8 3 8 . 3 0 7 8 . 1 3 9 8 . 0 8 9 8 3 0.7 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 5 . 9 3 2 0 9 . 0 3 4 7 6 . 1 0 6 2 8 . 4 9 4 6 .1741 .5401 . 0 2 9 3 5 - . 0 5 9 7 8 . 6 2 8 0 . 5 8 0 6 . 1 0 5 9 . 0 6 4 2 7 0.9 [ 0 , 1 0 ] 1 0 . 1 1 7 5 1 0 . 4 7 6 3 . 3 5 8 8 . 0 3 5 4 6 . 7 7 0 7 . 0 7 6 1 6 [ 1 0 , 2 0 ] 1 8 . 0 4 3 4 1 7 . 4 7 9 9 - . 5 6 3 5 - . 0 3 1 2 3 . 6 1 2 7 . 0 3 3 9 6 130 . Table 5a , continued.... v ( t ) : 0 . 1 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 1 7 8 7 . 1 8 6 2 . 1 8 3 4 . 1 9 0 8 . 0 0 4 6 4 7 . 0 0 4 5 3 8 . 0 2 6 0 0 . 0 2 4 3 7 . 0 1 3 9 9 . 0 1 4 1 2 . 0 7 8 2 6 . 0 7 5 8 3 0.3 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 1 . 5 7 1 3 1 . 7 4 2 0 1 . 7 4 5 0 1 . 9 0 0 2 . 1 7 3 7 . 1 5 8 2 . 1 1 0 6 . 0 9 0 8 3 . 3 3 9 6 . 3 3 9 8 . 2 1 6 1 . 1 9 5 1 0.5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 6 . 4 2 5 6 8 . 8 2 3 6 8 . 5 7 5 7 1 0 . 4 1 2 4 2 . 3 9 8 0 1 . 5 8 8 8 . 3 7 3 2 . 1 8 0 1 2 . 5 9 8 1 2 . 4 4 9 0 . 4 0 4 3 . 2 7 7 6 0.7 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 1 5 . 3 7 1 0 2 8 . 5 5 6 7 2 4 . 4 5 3 2 3 8 . 4 2 9 9 9 . 0 8 2 2 9 . 8 7 3 2 . 5 9 0 9 . 3 4 5 7 9 . 2 3 8 3 1 0 . 6 9 6 2 . 6 0 1 0 . 3 7 4 6 0.9 [ 0 , 1 0 ] [ 1 0 , 2 0 ] 2 5 . 5 9 5 0 5 7 . 8 2 3 3 4 5 . 5 8 2 1 8 9 . 3 6 7 6 1 9 . 9 8 7 2 3 1 . 5 4 4 3 . 7 8 0 9 . 5 4 5 5 1 9 . 9 8 7 2 3 1 . 5 4 4 3 . 7 8 0 9 . 5 4 5 5 2 31 . T A B L E 5b Period Averages and Mean Absolute Errors in n 0(t), m(t) and v(t) for Various Retrial Rates Using Method 2E. The entries in this table correspond to oscillatory parameter calculations using Method 2E with XQ = 2.0 and u(t) = 4.0, a constant, in all cases. The first column below gives the value of v0 and the second column gives the value of CQ = ^o/^o- The remaining seven columns are arranged as are the last seven columns of table 5a. 1 0 . 0 0 . 2 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 7 3 4 . 4 5 0 5 . 4 7 4 5 . 4 5 0 6 . 0 0 1 1 1 5 . 0 0 0 0 7 7 5 . 0 0 2 3 5 4 . 0 0 0 1 7 2 0 . 0 1 0 7 8 . 0 1 0 6 6 . 0 2 2 7 6 . 0 2 2 3 2 4 . 0 0 . 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] .4821 . 4 5 1 9 . 4 8 1 5 .4511 . 0 0 0 5 5 7 8 . 0 0 0 7 6 8 8 . 0 0 1 1 5 7 . 001701 . 0 2 7 4 2 . 0 2 4 8 5 . 0 5 6 8 7 . 0 5 4 9 8 2 . 0 1 .0 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 6 9 . 4 5 6 0 . 4 9 1 5 . 4 5 2 3 . 0 0 5 4 1 4 . 0 0 3 6 0 8 . 0 1 0 8 9 . 0 0 7 9 1 2 . 0 5 2 8 6 . 0 4 8 2 4 ,1064 . 1 0 5 8 1 .0 2 . 0 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 5 2 3 0 . 4 6 7 6 . 5 0 9 0 . 4 5 6 0 . 0 1 4 0 2 . 0 1 1 5 9 . 0 2 6 8 0 . 0 2 4 7 8 . 09094 . 0 8 2 0 0 . 1 7 3 9 .1754 0 . 4 5 . 0 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 5 7 1 4 . 5 0 3 2 . 5 5 5 2 . 4 7 4 2 . 01614 . 0 2 9 0 2 . 0 2 8 2 5 . 0 5 7 6 8 .1482 . 1290 . 2 5 9 4 . 2 5 6 4 0 . 2 1 0 . 0 [ 0 , 1 0 ] C 1 0 . 2 0 ] . 6 0 8 8 . 5 4 3 4 . 6 1 4 6 . 5 0 9 6 . 0 0 5 7 9 8 . 0 3 3 8 6 . 0 0 9 5 2 3 .06231 .1864 . 1 6 0 8 . 3 0 6 2 . 2 9 5 9 m ( t ) 1 0 . 0 0 . 2 [ 0 , 1 0 ] 1 .9893 1 .9789 - . 0 1 0 4 0 - . 0 0 5 2 2 7 .06651 . 0 3 3 4 3 [ 1 0 , 2 0 ] 2 . 1 6 5 5 2 . 4 1 7 0 - . 0 1 8 4 7 - . 0 0 8 5 3 1 . 0 7 2 1 2 . 03331 4 . 0 0 . 5 [ 0 , 1 0 ] 2 . 3 0 7 8 2 . 3 2 8 0 [ 1 0 , 2 0 ] 2 . 6 4 8 3 2 . 5 9 6 1 . 0 2 0 1 7 . 0 5 2 2 3 . 0 0 8 7 3 8 . 01972 .1663 .1484 . 0 7 2 0 5 . 0 5 6 0 5 2 . 0 1.0 [ 0 , 1 0 ] 2 . 7 7 4 6 2 . 8 2 0 0 [ 1 0 , 2 0 ] 3 . 4 2 6 0 3 . 2 6 1 1 . 0 7 5 3 8 . 1650 . 0 2 7 4 6 . 0 4 8 1 5 . 3 8 3 8 . 3 0 7 8 . 1 3 9 8 . 0 8 9 8 3 1 . 0 2 . 0 [ 0 , 1 0 ] 3 . 3 7 9 4 3 . 6 0 6 7 [ 1 0 , 2 0 ] 4 . 7 9 0 0 4 . 4 1 8 4 . 2 2 7 3 . 3 7 1 7 . 0 6 7 2 6 . 0 7 7 5 9 . 7 2 8 5 . 5 7 2 5 . 2 1 5 6 . 1 1 9 5 0 . 4 5 . 0 [ 0 . 1 0 ] 4 . 3 9 7 7 5 . 2 1 7 2 [ 1 0 , 2 0 ] 7 . 5 1 6 9 7 . 2 4 9 0 . 8 1 9 5 - . 2 6 7 9 . 1 8 6 3 . 0 3 5 6 3 1 . 2 9 7 8 .9721 .2951 . 1293 0 . 2 1 0 . 0 [ 0 , 1 0 ] 5 . 1 2 5 7 6 . 8 2 0 7 [ 1 0 , 2 0 ] 9 . 8 6 2 4 1 0 . 7 4 0 6 1 .6949 . 8 7 8 2 . 3 3 0 7 . 0 8 9 0 5 1 .7413 1 . 5 2 1 9 . 3 3 9 7 . 1 5 4 3 132 . Table 5b, continued.... V ( t ) : 10.0 0.2 [0,10] [10,20] 4.0 0.5 [0,10] [10,20] 2.0 1.0 [0,10] [10,20] 1.0 2.0 [0,10] [10,20] 0.4 5.0 [0,10] [10,20] 0.2 10.0 [0,10] [10,20] 5.3642 6.4752 6.0539 7.2249 5.8539 7.4099 7.1297 8.5744 6.4256 8.5757 8.8236 10.4124 7.0719 10.1098 11.3991 13.2489 7.8444 12.2899 15.2618 19.0427 8.3430 13.5835 17.6819 24.7199 1.1110 .2071 1.1710 .1934 1.5565 .2659 1.4447 .2026 2.3980 .3732 1.5888 .1801 3.0379 .4296 1.8499 .1623 4.4455 .5667 3.7809 .2477 5.2405 .6281 7.0380 .3980 1.1223 .2092 1.1834 .1955 1.7028 .2909 1.6789 .2355 2.5981 .4043 2.4490 .2776 3.8331 .5420 3.7199 .3263 5.3012 .6758 6.1736 .4045 5.7930 .6944 8.6585 .4897 133 . T A B L E 6 Period Averages and Relative Mean Absolute Errors in n o (t), m(t) and v(t) Using Method 2E With Various Amounts of Mixing, n , in the Closure Formula for m^(t). The three triads of numbers in the rightmost nine columns below are the approximate mean value of the quantity indicated (as calculated from closure method values), the relative error in this approximate mean value, and the relative mean absolute error over the intervals indicated, respectively, for three different oscillatory parameter retrial queues with the parameter values listed at the bottom of each section of this table. The number, n , tabulated in the first column, is the mixing factor, n , in equation (2.4.9) . (a) II „ (t) : 1.0 [ 0 . 1 0 ] [ 1 0 , 2 0 ] . 4 9 1 5 . 4 5 2 3 • . 0 1 0 8 9 - . 0 0 7 9 1 2 . 1064 . 1 0 5 8 . 2 5 2 9 . 1 4 5 2 - . 0 4 1 0 8 - . 1 5 4 3 . 3 7 2 3 .4381 . 6 1 4 6 . 5 0 9 6 . 0 0 9 5 2 3 - . 0 6 2 3 0 . 3 0 6 2 . 2 9 5 9 0 . 7 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] .4941 . 4 5 3 0 - . 0 0 5 7 6 6 - . 0 0 6 4 8 6 . 1 0 8 7 . 1 0 8 2 . 2 4 6 8 . 1 4 2 5 - . 0 6 4 5 5 - . 1 6 9 8 .3901 . 4 7 0 4 . 6 1 4 2 . 5 0 9 0 .008771 - . 0 6 3 3 9 .3074 . 2 9 6 6 0 . 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 6 6 . 4 5 3 6 - . 0 0 0 6 4 6 6 - . 0 0 5 0 4 1 . 1 1 5 9 . 1 2 2 4 . 2 5 2 8 . 1 3 1 8 - . 1 0 2 6 - . 1 9 5 6 . 4 1 1 9 . 5 0 1 5 . 6 1 2 8 . 5 3 0 3 . 0 0 6 5 4 5 - . 0 6 5 8 1 . 3 1 0 7 . 3 0 3 2 0 . 4 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 7 2 . 4 5 3 9 . 0 0 0 4 5 7 3 - . 0 0 4 3 9 3 . 1 1 9 8 . 1 3 1 5 . 2 3 0 9 .1351 - . 1 2 4 6 - . 2 1 2 6 . 4 1 9 9 . 5 1 2 2 0 . 3 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 7 3 .4541 . 0006131 - . 0 0 4 0 3 5 . 1 2 1 9 . 1 3 6 6 . 2 2 7 4 . 1 3 3 2 - . 1 3 7 8 - . 2 2 3 9 . 4 2 4 3 . 5 1 7 2 0 . 2 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 6 8 . 4 5 4 5 - . 0 0 0 2 6 0 8 - . 0 0 3 2 7 0 . 1 2 9 9 . 1 4 7 4 . 2 3 4 9 . 1 2 1 0 - . 1 7 0 5 - . 2 5 5 2 .4351 . 5 2 4 8 . 6 0 7 4 .5301 - . 0 0 2 4 4 9 - . 0 7 1 6 9 . 3 1 8 2 . 3 2 1 9 0 .1 [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 4 9 3 4 . 4 5 4 7 - . 0 0 7 1 9 4 - . 0 0 2 7 8 3 .1571 .1694 .1981 . 1 1 1 3 - . 2 4 8 9 - . 3 5 1 8 .4631 . 5 2 1 2 . 5 9 5 4 . 4 9 9 8 - . 0 2 2 1 2 - . 0 8 0 3 7 . 3 3 1 0 .3461 0 . 0 5 [ 0 , 1 0 ] [ 1 0 , 2 0 ] .4901 . 4 5 4 2 - . 0 1 3 7 2 - . 0 0 3 8 4 5 . 1 7 3 2 . 1 8 2 3 . 1 8 6 2 . 0 9 9 2 4 - . 2 9 4 2 - . 4 2 1 9 . 4 7 3 9 . 5 2 8 3 . 5 8 6 0 . 4 9 4 7 - . 0 3 7 4 7 - . 0 8 9 6 1 . 3 4 0 9 . 3 5 6 3 v b l e [ 0 , 1 0 ] [ 1 0 , 2 0 ] . 5 0 4 5 . 4 5 6 5 . 0 1 5 1 0 . 0 0 1 2 7 5 .1185 . 1 5 1 5 . 2 5 7 6 . 1 5 6 4 - . 0 2 3 3 0 - . 0 8 8 9 9 .3951 . 5 0 6 4 . 6 1 7 6 . 5 1 6 6 . 0 1 4 4 5 - . 0 4 9 3 9 . 3 1 2 7 .3181 E x a c t [ 0 , 1 0 ] [ 1 0 , 2 0 ] p o c 0 X 0 V 0 u . 4 9 6 9 . 4 5 6 0 1 .0 1 .0 2 . 0 2 . 0 4 . 0 . 2 6 3 8 . 1 7 1 7 0 . 9 1.0 3 . 6 2 . 0 4 . 0 . 6 0 8 8 . 5 4 3 4 1 .0 1 0 . 0 2 . 0 0 . 2 4 . 0 134 . Table 6, continued... (b) m(t) 1.0 0.75 0.5 0.4 0.35 0.25 0.1 0.05 vbte E x a c t [0.10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] 2.8200 3.2611 2.7710 3.2515 2.6836 3.2367 2.6294 3.2220 2.5959 3.2092 2.5106 3.1613 2.3004 2.9319 2.1819 2.7225 2.7856 3.6729 2.7446 3.4260 1.0 1.0 2.0 2.0 4.0 .02746 .1398 .04815 .08983 .009625 .1284 .05095 .09186 .02224 .05526 .04199 .05955 .05420 .06327 .08527 .07728 .1618 .1442 .2050 .2053 .01492 .07206 .1280 .1065 .1336 .1164 .1400 .1221 .1582 .1357 .2166 .1745 .2545 .2120 .1123 .1377 10.4763 17.4799 10.3148 17.2285 10.0988 16.8022 9.9897 16.5391 9.9286 16.3753 9.7895 15.9511 9.5172 14.8405 9.3953 14.1551 10.3706 17.9813 10.1175 18.0434 0.9 1.0 3.6 2.0 4.0 .03546 .03123 .01950 .04516 .07617 .03396 .07767 .04516 .001849 .08200 .06879 .06879 .01263 .08337 .01867 .09245 .03242 .1160 .05933 .1775 .07138 .2155 .08533 .08337 .08741 .09245 .09287 .1160 .1068 .1775 .1144 .2155 .02501 .07113 .003442 .03371 6.8207 10.7406 6.7776 10.6638 6.6984 10.5263 6.5200 10.2291 6.2550 9.7145 6.0817 9.2632 6.6968 10.8931 5.1257 9.8624 1.0 10.0 2.0 0.2 4.0 .3307 .3397 .08905 .1543 .3223 .3343 .08126 .1523 .3068 .3272 .06732 .1499 .2720 .03718 .2203 .01500 .1865 .06076 .3065 .1045 .3198 .1485 .3257 .1466 .3335 .1508 .3231 .1633 135 . Table 6, continued.... (c) v(t) [0.10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] [0,10] [10,20] E x a c t [0,10] [10,20] 8.5757 10.4124 8.1823 10.1664 7.4808 9.7751 7.0277 9.5223 6.7415 9.3387 6.0015 8.7741 4.2022 6.6994 3.2710 5.1419 9.1767 14.1814 6.4256 8.8236 1.0 1.0 2.0 2.0 4.0 .3346 .1801 .2734 .1522 .1642 .1078 .09371 .07918 .04918 .05838 .06600 .005612 .3460 .2407 .4909 .4173 .4282 .6072 .4043 .2776 .3032 .2016 .1645 .1098 .09371 .07918 .07762 .08096 .1333 .1186 .3489 .2566 .4929 .4173 .4282 .6072 45.5821 89.3676 40.2226 79.9766 33.1423 67.8332 29.6465 61.6341 27.7192 58.0919 23.4414 49.7807 15.6918 32.2777 12.6375 23.9212 45.9507 98.9332 25.5950 57.8233 0.9 1.0 3.6 2.0 4.0 .7809 .5455 .5715 .3831 .2949 .1731 .1583 .06590 .08299 .004645 .08414 .1390 .3869 .4418 .5062 .5863 .7953 .7110 .7809 .5455 .5715 .3831 .2949 .1731 .1583 .06590 .08299 .02388 .08515 .1390 .3874 .4418 .5066 .5863 .7953 .7110 13.5835 24.7199 13.3437 24.1721 12.8761 23.1664 11.6088 20.9121 9.4407 17.3672 8.0350 14.5742 13.7861 27.1265 8.3430 17.6819 1.0 10.0 2.0 0.2 4.0 .6281 .3980 .5994 .3671 .5433 .3102 .6944 .4897 .6511 .4423 .5657 .3538 .3914 .1827 .1316 .01779 .03691 .1758 .6524 .5341 .3914 .1899 .1365 .03091 .1083 .1758 .6524 .5341 1 36 . Table 6, continued.... (d) raQ(t) n 1.0 [0.10] [10.20] .5085 .5477 .1948 .2874 .5402 .4738 .7470 .8548 .5761 .6442 .6444 .6675 .3854 .4904 .8619 .9002 .8619 .9002 0.75 [0,10] [10,20] .4995 .5477 .2090 .2873 .4247 .3644 .7951 .9213 .5488 .6165 .5860 .6165 1.1089 1.5051 .6026 .4606 .6241 .6938 0.5 [0,10] [10.20] .4738 .5377 .2498 .3003 .3636 .3261 .8309 .9113 .5712 .6208 .6183 .6208 1.8010 2.5505 .3546 .4976 .5056 .4976 0.4 [0,10] [10,20] .4905 .5488 .2233 .2860 .3417 .3069 .4480 .5420 .7458 .7744 .7867 .7744 2.1212 2.9248 .2398 .4050 .4745 .4289 0.35 [0.10] [10,20] .4649 .5351 .2638 .3037 .3727 .3412 .3355 .3738 .8096 .8444 .8519 .8444 2.2252 3.1747 .2026 .3542 .4676 .4201 0.25 [0,10] [10,20] .3893 .4880 .3835 .3650 .4898 .4408 .6690 .7144 .6858 .6985 .7409 .7152 2.3869 3.4673 .1446 .3302 .4721 .4243 0.1 [0,10] [10,20] .2493 .3666 .6053 .5230 .7050 .6119 .2710 .3687 .8462 .8466 .9016 .8472 2.5424 3.5388 .08890 .2801 .5079 .4514 0.05 [0,10] [10,20] .1588 .2608 .7485 .6607 .8311 .7123 .1148 .1307 .9349 .9456 .9903 .9456 2.4843 3.4292 .1097 .3024 .5433 .4753 v b l e [0,10] [10,20] .6483 .7903 .02662 .02836 .3014 .2723 1.3016 1.4561 .2613 .3939 .3581 .4106 2.2512 3.0418 .1932 .3812 .4792 .4379 Exact [0,10] .6315 [10,20] .7685 1.7622 2.4026 2.7905 4.9159 1.0 1.0 2.0 2.0 4.0 0.9 1.0 3.6 2.0 4.0 1.0 10.0 2.0 0.2 4.0 (c) FIGURE 2.13 Relative Errors in n0(t), m(t) and v(t) Using Closure Method 2D for Various Values of P in the Constant Parameter Retrial Queue. The graphs were prepared from calculations based on equations (2.4.1) with (2.4.6). The six curves in each Figure correspond from lowest to highest to P =0.1, 0.3, 0.5, 0.7, 0.9, and 0.95, respectively. In Figures 2.13a and 2.13b, the relative error curves corresponding to P = 0.9 and P = 0.95 are effectively indistinguishable at the resolution employed. Appendix A gives details of calculations done to obtain exact values required in determining these relative errors. 138 . ( c ) F I G U R E 2.14 Relative Errors in II0(t), m(t) and v(t) Using Closure Method 2D for Various Values of c in the Constant Parameter Retrial Queue. The graphs were prepared from calculations based on equations (2.4.1) with (2.4.6). The six curves in each Figure correspond from lowest to highest to c = 0.2, 0.5, 1.0, 2.0, 5.0, and 10.0, respectively. As always in this thesis, the system is assumed to be empty at t = 0 (i.e., IIo(0) = 1.0, m(0) = 0.0 and v(0) = 0.0). Appendix A gives details of calculations done to obtain exact values required in determining these relative errors. 139 . F I G U R E 2.15 C o m p a r i s o n of E x a c t and A p p r o x i m a t e Q u a n t i t i e s f o r the Retrial Queue U s i n g M e t h o d 2 D . The calculations corresponding to these Figures were based on equations (2.4.1) and (2.4.6). Figure 2.15a applies to the constant parameter case, whereas the other three Figures on this page apply to oscillatory parameter cases. The emphasis is on v(t) here, and the accuracy of the closure approximation is clearly not great in any of the cases displayed. Appendix A gives details of calculations done to obtain these exact quantities. 140 . 0 5 10 15 20 25 0 5 10 15 20 25 (a) A =3.6 , V =2.0 (b) V =0 .2 , A =2.0 O O 0 o F I G U R E 2.16 Comparison of Exact and Closure Method 2D Approximation Values of ^ ( t ) and m ( 0 for Oscillatory Parameter Retrial Queues. For all of the calculations represented here, u(t) = 4.0 . The discrepancies between the exact and approximate H0(t) and m(t) are not as great as those between the exact and approximate v(t) as displayed in Figures 2.15c and 2.15d. However, particularly when v = 0.2 and A = 2.0, agreement between exact and approximation is still less than acceptable. Appendix A gives details of calculations done to obtain these exact quantities. 141 . t ( c) A = 2 . 0 , V = 0 .2 F I G U R E 2.17 Relat ive E r r o r s in v(t) C o m p u t e d U s i n g the Closure M e t h o d 2 E in the C o n s t a n t Parameter Case . T h e calculations represented here were done using equation (2.4.9) along with (2.4.1). Particularly when v is not smal l compared to A , the relative errofi are very sensitive to the value of n used i n (2.4.9). In Figure 2.17c, the relative error curves for n = 0.1 and n = 0-05 peak lower yet than those s h o w n , at values of approximately 0.54 and 0.52, respectively. Note that the vertical scale in Figure 2.17c is clipped. In all calculations represented here, u = 4.0 . 142 . 0 5 10 15 20 25 0 5 10 15 20 25 (c) V = 0 . 2 , A = 2 . 0 (d) v = 2 . 0 , A = 3 . 8 FIGURE 2.18 Performance of the Closure Method 2E When the Retrial Rate is Much Smaller Than the New Customer Arrival Rate. Figures 2.18a, b, c, illustrated the limited improvement represented by equation (2.4.9) when the retrial rate is significantly lower than the new customer arrival rate. The notation 'all n ' means the sequence of values n = 1, 0.75, 0.5, 0.25, 0.1, and 0.05, as well as the time varying n (t) = 1 - II (t). Figure 2.1 8d provides a contrast to the first three Figures here with a similar plot when A and v have more similar values. 143 . ( c ) A = 2 . 0 , V = 0 . 2 O o F I G U R E 2.19 Comparison of Exact and Method 2E Approximation v(t) for Various Values of n in Equation (2.4.9) in the Oscillatory Parameter Case. For all of the calculations represented here, u(t) = 4.0 . Appendix A gives details of calculations done to obtain these exact quantities. 144 . e x a c t ( a ) X = V O 0 2 . 0 ( b ) X = 3 . 6 , V = 2 . 0 o o 7 • 4-> 4 O e x a c t /Mr| v b l e n = ' 3 7 \ / A W ' \ | / n v b l e . n = i 20 ( c ) X = 2 . 0 , V = 0 . 2 0 O F I G U R E 2.20 Comparison of Exact m 0(t) with m^t) Given by Equation (2.4.9) for Selected Values of n in an Oscillatory Parameter Retrial Queue. Thereason equations (2.4.1) with (2.4.6) give the poor results observed earlier seems clear from comparing the n = 1 curves here with the exact m0(t). What is surprising is the accuracy of the results achieved when equation (2.4.9) is used with (2.4.1) and intermediate values of n • It is clear that equation (2.4.9) docs not generally constitute an acceptably accurate way to compute mQ(t) for its own sake. 145 . 2.5 Closure Methods With Three or More Equations Grandma Sonderquist's Conclusion: A chicken doesn't stop scratching just because the worms are scarce. [1001 Logical Laws..., Peers and Bennett) Three Equation Methods Between the two two-equation closure methods described and analyzed in the preceding two sections, we are able to compute approximate values for all four the quantities listed at the beginning of section 2.1 . Little has been said about the computation of Y^t) so far, but since both of these methods produce estimates of ^^(t) and m 0 ( 0 , it is easy to compute y JS). We have used only some of the equations listed in Table 2, and that Table can be extended indefinitely. A reasonably concerted search has been made for additional viable closure methods, perhaps based on systems of more than two differential equations. One such method is described below, but in general the results of this search have not been particularly encouraging. The main problem encountered in widening the scope of equations from Tabic 2 which were considered is in finding viable closure formulas. The quantities which must be eliminated from many of the equations not considered until now usually include moments or components of moments of order two or greater. These quantities have much more complicated steady state formulas, and they tend to vary over much wider ranges of values as a function of time. 146 . We end up trying to approximate higher order moments by complicated expressions involving lower order moments which are not only difficult to test comprehensively, but also difficult to justify intuitively. It is not surprising that such approximations are usually not simply poor, but lead to equations with solutions that violate some very basic requirements (such as the occurrence of negative values for moments of queue length, or probabilities with values exceeding 1). However, a three-equation closure method was found for computing (X), m(t), and Y0(t). We start with the first, second, and fourth equations in Table 2, modifying them trivially to get the system n ' ( t ) = - xn - yn - vn y o o 1 o' o m' ( t ) = X - \lTL1 , (1) ffl'n„ - m n' 0 0 0 0 where m 0 ( t ) = -XII o Y 0 - VC 0 + U J n 2 - . (2) Since we have the exact relations m 0 lli = 1 - nQ , m j = m - mQ , and T Q = ^ the only closure formula we need is to eliminate ?0(t) from equation (2). 147 . Now, at steady state, p c f p c + 1) ^0 1 - P and so with a little effort, we can come up with a list of potential closure formulas, including the two simplest, rQ = cmi , (3) since m =^ p(cp + 1)/(1 - p), giving Method 3A, and m (m + 1) c „ - S S - , (4, 0 n o since m 0 = cp, and ftQ = 1 - p , giving Method 3B. A number of more complicated flights of fancy were included briefly in test calculations, but all resulted in catastrophic numerical breakdown of the solution of the resulting system of differential equations within one or two steps of the starting point . What formula (3) has going for it is simplicity. In constant parameter calculations with A = 2.0, u = 4.0, and v = 2.0, results were almost misleadingly good. Figure 2.21 includes a plot of relative errors as well as plots comparing the exact and approximate quantities for this set of constant parameters. Note from Figure 2.21d that the approximate m 0(0 takes a bit of a negative dive for a short interval initially, and that the relative error in m (t) is 148 . rather large initially.1 Figure 2.22 shows what happens when the corresponding oscillatory parameter calculation is done using equation (3). There is no real need to discuss magnitudes of relative errors here, though, for what it's worth, the first clipped peak in Figure 2.22d actually rises to a value of approximately 150, which is around 60 times the true value of Y^C1) a t that point. Reservations regarding use of formula (3) above in a closure method for any sort of calculations are difficult to avoid in view of Figure 2.22. This is perhaps a good example of the risks in using a closure formula based only on simple algebraic coincidence between steady state formulas, since it is difficult to think of reasons for z; to vary in simple direct proportion to m^ in general. On the other hand, the formula (4) leads to a very accurate closure method here. This could well be due to the fact that not only is (4) true at steady state in the constant parameter case, but one could perhaps argue that r, will vary in 2 a manner similar to the square of . At any rate, for the constant parameter retrial queue with A = 2.0, u = 4.0, and v = 2.0, relative errors in II , m, and m ,^ never exceed a magnitude of 0.03, and when plotted on the same grid at the level of resolution employed in figures in this thesis, graphs of the approximate and exact quantities are virtually indistinguishable. Figure 2.23 compares plots of exact and approximate values of II , m, and m ,^ using (4) for the oscillatory f igure 2.21a has been clipped somewhat. At t = 0.5, the first positive value of t for which an approximate value of m (^t) was computed, the relative error in m^t) is -1.7. While this is a rather large value, the magnitude of the relative error does decrease rapidly to a range of values comparable with the relative errors in the other two quantities being computed. 149 . parameter case with A ^ = 2.0, u = 4.0, and = 2.0. Clearly, when formula (4) is combined with equations (1) and (2), a closure method results which is of comparable accuracy to those methods discussed in previous sections.2 Closure Methods Involving Four or More Equations Two attempts were made to find viable four-equation closure methods. One attempt involved the equations for IT^(t), m'(t), C^( r )> a n o " from Table 2 (a round-about way to get at II m, and v). This required a closure formula for m^ and , the server-idle component of the third order moment of the system. A second line of work involved the equations for (t), m 0 ' ( t ) , m'(t), and v'(t), which required closure formulas for m and r . A l l attempts here 1 0 produced ludicrous results. A five equation system using the equations for II '(t) , m '(t), m'(t), 0 0 C 0 ( t ) , and £ j ( t ) was tested briefly, but abandoned on the basis of early results. This system required only a single closure formula for , but none of the four such potential expressions tested showed any promise at all. 2 Actual ly , one can use equations (1) and (2) as written to solve f o r l l 0 (t), m(t) and Yn(t), or simply use the original equations for n^t) , m'(t) and m0'(t) in Table 2 with formula (4) to solve for II (t), m(t) and m (t), from which yQ(t) can be computed easily. 0 150 . The calculations represented here correspond to A = 2.0, y = 4.0 and v = 2.0. Appendix A gives details of calculations done to obtain the exact quantities used in the preparation of these plots. 151 . The calculations represented here correspond to X 0 = 2.0, u(t) = 4.0, ^0 = 2.0 and a period of oscillation of 10.0 time units. Appendix A gives details of calculations done to obtain these exact quantities. 152 . The calculations represented here correspond to X0 = 2.0, u(t) = 4.0, v0 = 2.0 and a period of oscillation of 10.0 time units. Appendix A gives details of calculations done to obtain these exact quantities. 153 . 2.6 Summary and Conclusions In each of sections 2.2 through 2.5, we have described and analyzed one or more distinct closure methods for the single server retrial queue along with computations illustrating the accuracy of the approximate quantities produced and discussed why we might expect a certain degree of accuracy in each case. Of the methods studied, three were found to give results of the same degree of accuracy as the closure method due to Rothkopf and Oren [1979]. These are summarized in Table 7. Only one of these three serviceable methods produces estimates of the variance of the queue length, v(t), whereas all three give estimates of the probability that the system is empty, ^ ( t ) , and the mean number of customers in the system, m(t). It is also possible to compute estimates of m 0 ( t ) , the server-idle contribution to the mean number of customers in the system or Y 0 ( 0 , the conditional expected number of customers in the system given the server is idle, using any of the three best methods described, but the accuracy of these estimates has been examined only for Method 3B. The quantities ^^(t) and m(t) can be obtained using any of the three closure methods found viable. Experience so far does not indicate that any one of these three methods is markedly superior to the others if I I (t) and m(t) are the only quantities of interest. Use of Method 2E in this case may be attractive, since that method has been studied most comprehensively. 154 . The primary value of these closure techniques is in the reduction in computational labour they allow compared to determining the same quantities 'exactly' by solving an appropriately truncated subsystem of equations (1.1.3). Instead of solving a system of dozens or perhaps hundreds of simultaneous ordinary linear differential equations, we can get good approximations for these quantities by solving a system of two or three such equations only. The reduction in dimension of the systems to be solved when closure methods are used brings their solution within the capability of currently available microcomputer systems. This work also demonstrates that Rothkopf and Oren's closure approach can be usefully modified to apply to more complicated queuing regimes, and no doubt to other probability models involving systems of ordinary linear differential equations. We do note however, that the work in this thesis applies only to the single server retrial queue, with all inter-arrival, inter-retrial, and service times assumed exponentially distributed. The work of Rothkopf and Oren would seem to imply that extension of these techniques to multi-server retrial queues with arrival and service times exponentially distributed would be easier to do than to allow alternative and inevitably more complicated probability distributions. One would not be surprised if the development of closure methods turned out to be very difficult if not impossible for systems which could not be described by a system of ordinary linear differential equations. In Chapter 1, mention was made of a number of alternative approaches proposed for non-exponential queue systems. 155 . T A B L E 7. S u m m a r y of Serv i ceab le C losu re Me thods f o r C o m p u t i n g the T i m e -Dependent B e h a v i o u r o f the S i n g l e - S e r v e r R e t r i a l Queue . M E T H O D 2 C : II (t), m ( t) ( T w o E q u a t i o n s ) o n^ct; = -n r t ; + yn f t ; - v m f t ; + vm f t ; m'(t) = x - yn f t ; (2.3.1) w i t h Jl1(t) = i - n Q f t ; c n (t)+i TI (t) {cJi (t) +i] m f t ; = n m f t ; + f i - n ; f 2 . j . 6 ; c+1 ] [ o ( t ) ( C a l c u l a t i o n s r e p o r t e d i n t h i s t h e s i s u s e d t h e a r b i t r a r y v a l u e n = 1/2, a n d n o a t t e m p t w a s m a d e t o d e t e r m i n e i f t h i s w a s a n o p t i m a l v a l u e t o u s e f o r n , n o r w a s a n y s e a r c h m a d e t o f i n d v a l u e s o f n l e a d i n g t o g r e a t e r a c c u r a c y . ) M E T H O D 2E: n^f t ; , m(t), v(t) ( T w o E q u a t i o n s ) m' f t ; = x - y { i - n f t ; } (2.4.1) v'(t) = X + y - y n o f t ; { 2 m f t ; + 1} + 2umo(t) wi th m2/{(v-m) (c + 1) } if v m (2.4.4) Jlo(t) -m/(c+1) . -e ' if v = m (2.4.5) mQ(t) = n c { i - n o f t ; } + f i - n ; {mf t; n f t; - n f t ; (2.4.9) n f t ; = i - n r t ; (when X - V, choose T) between 0.25 and 0.5, when X >> V, choose n = 0.1) 156 . T A B L E 7, continued., M E T H O D 3B : J[ (t), m (t) , m(}(t) o r Y0(t) ( T h r e e E q u a t i o n s ) n 0' (t) = -XU0(t) - \illi(t) - vJ[Q(t)yo(t) m* (t) = X - y l l (t) (2.5.1) m' (t) = - \ H n ( t ) y n ( t ) - vr, (t) + yzn (t) - yl l (t) (2.5 o 0 0 0 1 l 2) O I m'Q(t) HQ(t) - mo(t)Tl'o(t) y ' n ( t ) ~ 2 ° {mo(t)V wi th n1(t) = i - j[Q(t) m1(t) = m(t) - mQ (t) y0(t) = n,o(t)/no(t) m (t) {m (t) + 1 } r,o(t) = —2 °- (2.5.4) JlQ(t) 157 . APPENDIX A DETAILS OF CALCULATION OF EXACT QUEUE PROPERTIES As described in Section 1.3, exact values of time dependent probabilities required to compute exact properties of the various queues considered were determined by solving truncations of the systems (1.1.3a,b) for the M/M/l/1 retrial queue and (1.4.1) for the simple M/M/l queue using the RKF45 algorithm. The dimension of the actual system solved numerically was set by specifying a maximum value of n (the number of customers in the system for the simple M/M/l queue and the number of customers in orbit for the M/M/l/1 retrial queue). The RKF45 algorithm then extrapolated a numerical value for each of the probabilities corresponding to all values of n smaller than the specified maximum value, beginning with a fixed initial value at time zero. In all cases, solution estimates were obtained in steps of At = 0.5 starting with t = 0. At the end of each such step, all of the computed probabilities were summed, and the residual obtained by subtracting this sum from unity was a measure of the error caused by solving only the finite subset of the infinite system of differential equations describing the queue. The maximum value of n to be considered was generally chosen so that this residual probability never became greater than 10-6, though in some calculations used, this residual attained larger values over brief time intervals. To partially compensate for this truncation error, the residual was added into the maximum index probability for the M/M/l queue or pro rata into the maximum index probabilities for the M/M/l/1 retrial queue. The actual maximum values of n used in calculation of exact probabilities for test calculations reported in this thesis and the maximum values attained by the probability residuals in each case are tabulated below. It should be noted that the quoted residuals are the maximum rather than typical values that occurred over the entire time range considered. Simple M/M/l Queue: P Constant Parameter Oscillatory Parameter 0.1 25 30 30 35 65 80 <10-i2 <10-i2 25 30 30 50 70 <10-i2 0.3 0.5 0.7 0.9 0.95 4.6 x 10-10 8.7 x 10-7 5.3 x 10-8 1.4 x 10-9 1.0 x 10-i2 5.6 x 10-7 2.6 x 10-7 2.6 x 10-5 158 . M / M / l / 1 Retrial Queue fc = 1): Oscillatory p Constant Parameter Parameter 0.1 20 <10-i2 20 6.0 x IO-" 0.3 20 2.4 x 10-n 20 8.1 x 10-8 0.5 30 3.0 * 10-9 29 3.9 x 10-5 0.7 35 6.1 x 10-6 50 2.2 x 10-5 0.9 50 5.0 x 10-5 65 5.4 x 10-4 0.95 60 1.1 x 10-5 M / M / l / 1 Retrial Queue f p = 0.5): Oscillatory c Constant Parameter Parameter 0.2 20 4.1 x 10-7 20 2.5 x 10-4 0.5 20 8.4 x 10-7 20 4.4 x 10-4 1.0 30 3.0 x 10-9 29 3.9 x 10-5 2.0 25 4.7 x 10-7 25 4.1 x 10-4 5.0 30 4.1 x 10-7 30 5.0 x 10-4 10.0 35 2.3 x 10-7 35 4.0 x 10-4 159 . BIBLIOGRAPHY Alexsandrov, A. M., [1974], "A Queueing System With Repeated Orders", Engineering Cybernetics Rev., 1_2, pages 1 - 4 . Burden, Faires, and Reynolds, [1978], Numerical Analysis, Prindle, Weber and Schmidt, Boston. Carter, G. M., and Cooper, R. B., [1972], "Queues with Service in Random Order", Operations Research, 20, pages 389 - 405 . Choo, Q. H., [1978], The Interaction of Theory and Simulation in Queueing Analysis, Ph. D. Thesis, Chelsea College, University of London. Choo, Q. H., and Conolly, B., [1979], "New Results in the Theory of Repeated Orders Queueing Systems", J. Appl. Prob., J_6, pages 631 - 640 . Clark, G. M , [1981], "Use of Polya Distributions in Approximate Solutions to Nonstationary M/M/s Queues", Communications of the ACM, 24, pages 206 - 217 . Cohen, J. W., [1957], "Basic Problems of Telephone Traffic Theory and the Influ ence of Repeated Calls", Phillips Telecomm. Rev. 18. pages 49 - 100 . Cooper, Robert B., [1981], Introduction to Queueing Theory, Second Edition, North Holland, New York. Falin, G. I., [1979], "A Single-Line System With Secondary Orders", Engineering Cybernetics Rev., J_7, pages 76 - 83 . Felberg, E., [1970], "Klassiche Runge-Kutta Formeln Vierter und Neidrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Warmeleitungs Probleme.", Computing, 6, pages 61 - 71 . Gross, Donald, and Harris, Carl M., [1974], Fundamentals of Queueing Theory, John Wiley & Sons, New York. Hashida, O., and Kodaira, K., [1975], "Repeated Calls Resulting from Data Terminal Busy and Delayed Delivery Service", Electronics and Communications in Japan, 58A. pages 1 - 9 . Keilson, J., Cozzolini, J., and Young, H., [1968], "A Service System With Unfilled Requests Repeated", Operations Research, 6, pages 1126 - 1137 . Keller, J. B., [1982], "Time-Dependent Queues", SIAM Review, 24, pages 401 - 412. Kelton, W. D., and Law, A. M., [1985], "The Transient Behaviour of the M/M/s Queue, With Implications for Steady-State Simulation", Operations Research, 33, pages 378 - 396 . 160 . Kendall, D. G., [1953], "Stochastic Processes Occurring in the Theory of Queues and Their Analysis by the Method of Imbedded Markov Chains", Ann. Math. Statist., 24, pages 338 - 354 . Koopman, B. O., [1972], "Air-Terminal Queues Under Time-Dependent Conditions", Operations Research, 20, pages 1089 - 1114 . Kulkarni, V. G., [1982], Letter to the Editor, J. Appl. Prob., 19, pages 901 - 904 . Kulkarni, V. G., [1983], "On Queueing Systems With Retrials", J. Appl. Prob., 20, pages 380 - 389 . Lyashenko, N. N., and Nikulin, M S., [1986], Computer and Some Statistical Inference, Technical Report #32, Department of Statistics, University of British Columbia. Massey, W. A., [1985], "Asymptotic Analysis of the Time Dependent M/M/l Queue", Mathematics of Operations Research, 1_0, pages 305 - 327 . Morse, P. M., [1958], Queues, Inventories, and Maintenance, John Wiley & Sons, New York. Newell, G. F., [1968a], "Queues With Time-Dependent Arrival Rates I - The Transition Through Saturation", J. Appl. Prob., 5, pages 436 - 451 . Newell, G. F., [1968b], "Queues With Time-Dependent Arrival Rates II -- The Maximum Queue and The Return to Equilibrium", J. Appl. Prob., 5, pages 579 - 590 . Newell, G. F., [1968c], "Queues With Time-Dependent Arrival Rates III -- A Mild Rush Hour", J. Appl. Prob., 5, pages 591 - 606 . Odoni, A. R., and Roth, E., [1983], "An Empirical Investigation of the Transient Behaviour of Stationary Queueing Systems", Operations Research, 34, pages 432 - 455 . Ross, Sheldon M., [1983], Stochastic Processes, John Wiley & Sons, New York. Rothkopf, M. H, and Oren, S. S., [1979], "A Closure Approximation for the Nonstationary M/M/s Queue", Management Science, 25, pages 522 - 534 . Saaty, T. L., [1961], Elements of Queueing Theory with Applications, McGraw-Hill, New York. Shanks, E. B., [1966], "Solutions of Differential Equations by Evaluations of Functions", Math. Comp., 20, pages 21 - 38 . Shiryayev, A, N., [1984], Probability, Springer-Verlag, New York. 161 . Stidham, S., [1972], "L = MV: A Discounted Analogue and a New Proof", Operations Research, 20, pages 1115 - 1126 . Stidham, S., [1974], "A Last Word on L = \W\ Operations Research, 22, pages 417 - 421 . Zachmanoglou, E. C. and Thoe, D. W., [1976], Introduction to Partial Differential Equations With Applications, Williams & Wilkins, Baltimore.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Closure methods for the single-server retrial queue
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Closure methods for the single-server retrial queue Sabo, David Warren 1987
pdf
Page Metadata
Item Metadata
Title | Closure methods for the single-server retrial queue |
Creator |
Sabo, David Warren |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | This work focuses on the development and evaluation of so-called "closure methods" for solving the equations governing the time-dependent behaviour of single-server retrial queues. These methods involve assuming that particular known algebraic relationships between various characteristics of the corresponding steady-state queue also apply approximately when the queue is not at steady-state. The objective is to replace a problem requiring the solution of dozens or hundreds of simultaneous linear differential equations with a system of a few differential equations that has a solution that approximates those queue characteristics of immediate interest. The viability of such closure methods is assessed by examining the results of a series of test calculations. The methods described in this thesis apply to a retrial queue in which inter-arrival times for new customers, inter-retrial times, and service times are all assumed to be exponentially distributed. The steady-state solution for such a queue is described in some detail. A survey of the literature indicates that the description of this steady-state retrial queue has become quite sophisticated, whereas only very tentative steps have been taken in the study of the time-dependent behaviour of such queues. On the other hand, the time-dependent behaviour of the simple M/M/s queues have been studied to a much greater extent. The apparent value of closure methods in computing approximations to various basic time-dependent M/M/s queue characteristics motivated this examination of the extension of such methods to the single-server retrial queue. After discussing the basic approach to be used in devising and testing prospective closure methods for the single-server retrial queue, a variety of such methods is presented, with each being tested in considerable detail. It is found that three of the methods devised give results of comparable or better accuracy than those closure methods for the simple M/M/s queues which motivated this study. All recommended closure methods developed here involve systems of either two or three differential equations and permit the calculation of good approximations to four of the characteristics of greatest interest for non-stationary queues: the probability that the server is idle, the mean queue length, the variance of the queue length, and the conditional mean number of customers in the system given that the server is idle. Each of the methods presented is tested for queues with constant mean arrival, retrial and service rates, as well as for queues in which arrival and retrial rates vary sinusoidally with time. |
Subject |
Queuing theory |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080170 |
URI | http://hdl.handle.net/2429/26528 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A6_7 S22.pdf [ 7.06MB ]
- Metadata
- JSON: 831-1.0080170.json
- JSON-LD: 831-1.0080170-ld.json
- RDF/XML (Pretty): 831-1.0080170-rdf.xml
- RDF/JSON: 831-1.0080170-rdf.json
- Turtle: 831-1.0080170-turtle.txt
- N-Triples: 831-1.0080170-rdf-ntriples.txt
- Original Record: 831-1.0080170-source.json
- Full Text
- 831-1.0080170-fulltext.txt
- Citation
- 831-1.0080170.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080170/manifest