Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A queueing analysis of a multichannel, integrated voice and data communications system Haller, Dennis Raymond 1982-04-13

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


831-UBC_1982_A7 H34.pdf [ 4.19MB ]
JSON: 831-1.0065616.json
JSON-LD: 831-1.0065616-ld.json
RDF/XML (Pretty): 831-1.0065616-rdf.xml
RDF/JSON: 831-1.0065616-rdf.json
Turtle: 831-1.0065616-turtle.txt
N-Triples: 831-1.0065616-rdf-ntriples.txt
Original Record: 831-1.0065616-source.json
Full Text

Full Text

A QUEUEING ANALYSIS OF A MULTICHANNEL, INTEGRATED COMMUNICATIONS SYSTEM VOICE AND DATA by DENNIS RAYMOND HALLER B.E., University of Saskatchewan, 1980 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Electrical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June 1982 © Dennis Raymond Haller, 1982 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Electrical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: 28 June 1982 i i Abstract A multichannel radio communications system is modeled as a multiserver queue, with two distinct customer classes representing voice and data messages. Data, the class with the shorter average length, is given non-preemptive priority over voice. The queueing model is analyzed as a continuous-time Markov chain with an infinite state space. The infinite set of steady state balance equations is truncated, then solved numerically using the linear programming (LP) technique of Kotiah. Upper and lower bounds are thereby obtained for the mean waiting times of each customer class, and the probability distribution for the number of messages in the system. Exploitation of the Markov chain's property of irreducibility improves the original algorithm by considerably reducing the computational cost. Simulation is used to help analyze the system and to validate the numerical results. The particular four-channel case of the queueing model is treated in detail; both simulation and numerical results are presented. The LP method produces excellent results when the data traffic intensity is less than about 0.1 . This corresponds typically to the proportion of data messages being less than 90%. For other traffic mixtures, the upper bounds, especially for waiting times, are often too high to be of practical use. However, the lower bounds of the mean waiting times are of greater use; the simulation results show them to be reasonably good approximations to the true steady state values. ... X 111 Table of Contents Abstract ii List of Figures v Acknowledgement vI. INTRODUCTION 1 II. SYSTEM ANALYSIS 3 A . THE APPROACHB. MARKOV MODELINGC. NUMERICAL TECHNIQUES OF MARKOV MODELING 6 1. The Markov Chain Equations 6 2. Finite State Space 7 3. Infinite State Space 8 III. MODEL SELECTION 11 A. BASIC ASSUMPTIONSB. SERVICE DISCIPLINES . . 12 C. ALTERNATING PRIORITY MODELS 14 1. Previous Investigations2. Markov Chain Model — Non-exhaustive Service ..15 3. Markov Chain Model - Exhaustive Service 16 D. NON-PREEMPTIVE PRIORITY MODELS 17 1. Previous Investigations2. Markov Chain Model 18 IV. NUMERICAL SOLUTION TECHNIQUE 20 A. THE STATE SPACE 2B. THE BALANCE EQUATIONS 21 C. THE LINEAR PROGRAMMING SOLUTION TECHNIQUE 24 D. NORMALIZATION CONSTRAINTS 26 E. IRREDUCIBILITY OF THE MARKOV CHAIN 28 1. The Significance of Irreducibi1ity 29 2. Improvements to the LP Algorithm 30 F. IMPLEMENTATION OF THE MODIFIED LP ALGORITHM 31 V. RESULTS 35 A. SYSTEM PARAMETERS 35 B. THE STEADY STATE PROBABILITY DISTRIBUTION p(n) ..38 C. MEAN WAITING TIMES 43 1 . Definition2. The Limiting Cases: a,-M 44 3. Waiting Time Results 46 D. WAITING TIME SIMULATIONS 51 E. COMMENTS ON THE WAITING TIME RESULTS 53 1. Empirical Observations 52. Waiting Times in a Non-preemptive Priority System 54 F. CONVERGENCE OF UPPER AND LOWER BOUNDS 56 G. COMPUTATIONAL COSTS 63 1. LP Solution Cost2. A Comparison with Simulation Costs 65 VI. CONCLUDING REMARKS 67 A. LIMITATIONS OF THE MODEL 6B. SUMMARY AND CONCLUSIONS 9 C. SUGGESTIONS FOR FURTHER STUDY 72 APPENDIX A - A COMPARISON OF NON-PREEMPTIVE PRIORITY AND ALTERNATING PRIORITY SERVICE DISCIPLINES, 74 APPENDIX B - THE STATE SPACE TRANSITION DIAGRAM 77 APPENDIX C - SETTING UP THE BALANCE EQUATION COEFFICIENT MATRIX "A": AN EXAMPLE 79 APPENDIX D - A DERIVATION FOR THE NORMALIZATION CONDITION GIVEN IN [KOTI 77] 84 BIBLIOGRAPHY 86 Glossary of Notation 89 V List of Figures 1. Relationship Between Parameters &, and o, 37 2. The Distribution p(n) for p = 0.7, o,= 0.80, (*,= 0.14) 39 3. The Distribution p(n) for p = 0.5, o,« 0.80, (*,= 0.14) 41 4. The Distribution p(n) for p = 0.5, oy= 0.95, 0.44) 42 5. Waiting Times as a Function of *i : P - 0.7, n = 30 c. 47 6. Waiting Times as a Function of e, : p = 0.5, n = 25 c. 48 7. Waiting Times as a Function of B, : p - 0.3, n = 20 c. ....49 8. Convergence of Bounds for Probabilities p(0), p(4) 58 9. Convergence of Waiting Time Results: p - 0.7, o,= 0.80, (*,,= 0.14) 59 10. Convergence of Waiting Time Results: p = 0.5, o,= 0.80, 0.14) 60 11. Convergence of Waiting Time Results: p = 0.5, o,= 0.95, (p,= 0.44) 1 12. Schematic of a 2-Class Multiserver Queueing System ..74 13. State Space for n Z s 78 14. State Space for n < s15. Matrix "A" For nc=6, s=4 Servers 80 vi Acknowledgement I would like to thank my supervisor, Dr. C.S.K. Leung, for his guidance and constant availability during the course of the research and writing of this thesis. Financial assistance in the form of an NSERC Postgraduate Scholarship, and a University of British Columbia teaching assistantship is also gratefully acknowledged. This work was partially supported by NSERC research grant A-1052. 1 I. INTRODUCTION In this study we analyze a centrally co-ordinated, multichannel radio communications system. Service is provided to both data messages and voice messages, and each type can be queued when all radio channels are busy. One example of such a system may be a modern mobile radio network in which each mobile can communicate with a central office either by voice or through an interactive computer terminal. Requests for transmission of messages are sent via an auxiliary channel to a central controller. If all the message-carrying radio channels are busy, the controller temporarily denies incoming requests and allows a queue to form. In this way, the controller imposes an orderly service discipline on the use of a finite resource — the s available channels. Our objective is to model this communications system to obtain some measure of the average queueing times for each message type and the steady state probability distribution for the number of messages in the system. In a model based on queueing theory the radio channels become servers, while the data and voice messages become customers of types 1 and 2 respectively. The queueing analysis will in fact require the numerical solution of the continuous-time Markov chain state equations. We will specifically obtain numerical results for the case of s = 4 channels. An auxiliary signaling channel is assumed to exist, but it will not be included in our model. The service discipline of the central controller must not only efficiently assign channel use, but also must maintain 2 message queueing times at acceptable levels. Among the many possible service disciplines, the most relevant to the given physical system are: first-come first-served (no priority), cutoff priority, non-preemptive priority, and alternating priority. Only the last two of these will be given serious further consideration. The following chapter describes the analytical approach of this study, and outlines the numerical techniques of Markov modeling. Chapter III reviews previous work done with the alternating priority and non-preemptive priority models, and assesses both with respect to their analytical and numerical tractability. Chapter IV describes in detail the numerical methods use to solve the chosen model, and numerical results for a selected example are presented and discussed in Chapter V. Finally, in Chapter VI the main results from this study are summarized and possible future work is outlined. 3 II. SYSTEM ANALYSIS A. THE APPROACH Analytical solutions are generally unavailable for most multiserver queueing systems with priority service disciplines. Two of the more popular alternatives for systems analysis are simulations and numerical solutions, particularly the solution of Markov chain models. The tradeoff between numerical solution and simulation is in initial analytical effort versus actual program execution time. A numerical solution may require a significant modeling effort, while simulation usually requires only the straightforward translation of system behaviour into a specialized simulation language. However, the length of time a simulation program must run to give acceptable confidence intervals for the measured quantities is usually greater than the execution time of a numerical solution program. As with any computer program, both numerical and simulation programs require independent checks on their validity, it being acceptable to use one as a check on the other. In this investigation we have relied primarily on the numerical approach but have also used simulation to provide additional verification of the results. B. MARKOV MODELING Many probabilistic systems, in which incoming customers demand use of certain system resources, can be analyzed effectively using a Markov chain model. In the following we present a brief summary of Markov chain modeling techniques, and 4 for simplicity, we consider only stochastic systems with time-independent parameters (stationary processes). More detailed treatments can be found in [BHAR 60], [KLEI 75], and [ROSS 72]. The essential feature of a Markov chain whether in discrete-time or continuous-time is the Markov Property. That is, a Markov chain has the property that its present state completely summarizes how its entire history affects its future. A Markov chain is formed of a discrete set of states, each state being a unique configuration of the state variables chosen to characterize the system. There will be a finite number of state variables, say d , and each variable may, depending on the system, take on either a finite or infinite number of discrete allowable values. The d state variables form a d - dimensional state space and, if at least one variable has an infinite range, then there will be an infinite number of states in the state space. For example, in a queueing system with unrestricted queue lengths the Markov chain model will be defined on an infinite state space. The characteristics of the system underlying the Markov chain model will also determine whether a discrete-time or continuous-time representation can be used. In a continuous-time Markov chain, the time between successive state transitions is an exponentially distributed continuous random variable with a known (perhaps state dependent) parameter. An example of a system which can be so modeled is the M|M|s queueing system. The instants at which state transitions may occur in a discrete-time Markov chain are labelled by the integers {0,1,2,... }, and 5 the time between successive state transitions is a geometrically distributed discrete random variable. Discrete-time Markov chain models may be applied to systems which inherently operate in discrete-time, for instance, a digital communications system. If the distribution of time between state transitions is arbitrary (either in continuous or discrete time), then this process is called semi-Markov. The process can still be modeled as a Markov chain if a certain set of time instants can be found which together define an imbedded discrete-time Markov chain. This imbedding technique may be used to model many continuous-time stochastic processes as discrete-time Markov chains. The M|M|s queueing system for example can easily be represented as a continuous-time Markov chain, but a discrete-time model is also possible by imbedding a Markov chain on the set of instants formed of all state transition times. The M|G|1 continuous-time process is not Markov, nor even semi-Markov (with respect to number in the system n(t), t^O), but it may be successfully modeled as a discrete-time Markov chain imbedded on the set of instants consisting of all service completion times. Solving these imbedded discrete-time Markov chains however, gives not the general-time probability distribution of number in the system, but instead a distribution for proportion of time spent in each state of the imbedded chain. The computed probabilities therefore need to be adjusted so that they again correspond to stationary probabilities of the system being in each state. This transformation makes use of the relationship between states of the particular imbedded chain and those in the 6 original general-time system [GROS 74, p.311 ii.]. C. NUMERICAL TECHNIQUES OF MARKOV MODELING 1. The Markov Chain Equations A numerical solution of a Markov chain model aims to calculate the stationary probability that the system is in each state of the state space. From these probabilities it is then possible to compute queue length distributions, moments of waiting time distributions, and other quantities of interest. In a continuous-time Markov chain, an equation is written for each state, balancing the total transition rate (or flow) into the state with the total transition rate out of the state. Denote by tr the row vector of steady state probabilities, each component rri representing the proportion of time the system spends in state i . The set of balance equations is then written (column by column) as: £ Q = 0 (2.1) with Q the matrix of transition rates, and fJ the null (row) vector. Since jr is a probability distribution it is normalized by 1± ir± = 1 . (2.2) In a discrete-time Markov chain a transition probability matrix P is defined such that P is the probability the process will at the next time instant move to state j, given that it is currently in state i. The steady state probability distribution is defined through ir P = n (2.3) 7 and the normalization (2.2). Here the steady state probabilities can be interpreted as the proportion of transitions that take the process into state i . Note that (2.3) may be rewritten in a form structurally similar to (2.1) although the steady state vector ir is not defined in the same way: jr(I-P) = 0 (2.4) If the state space is infinite then (2.1) or (2.4) represents an infinite number of linear equations. 2. Finite State Space Assume first that the state space is finite (N states) so that both (2.1) and (2.4) are (N x N) systems of linear equations. Actually only (N-1) of these equations are independent, so the normalization (2.2) must be used to replace any one of the N equations in (2.1) or (2.4). The linear systems may then be solved exactly to obtain steady state probabilities. Solutions may be obtained either by direct or iterative methods. Direct solutions are usually based on Gaussian elimination and are more efficient when the coefficient matrices are sparse, because then specialized data structures may be used to avoid a computations which involve zeros or ones [DUFF 77]. There are a great variety of iterative methods [STEW 78] the simplest of which is power-iteration [WALL 66], In the method of power-iteration, as applied to the solution of (2.1), the (k+1)'th iterate for the solution vector jr is obtained from the k'th 8 iterate as follows: JLk+1 = JLk (AQ + I) (2.5) where A is a scalar constant. The parameter A controls the rate of convergence and, although guidelines exist for choosing A [STEW 78], its best value must usually be determined experimentally. Iterative methods are in general more efficient than direct methods whenever the coefficient matrices Q or (I - P) are full, or when they are sparse but have a structure that would force a Gaussian elimination algorithm to fill-in many of the zero elements. 3 . Infinite State Space When the state space is infinite, the linear system of equations (2.1) or (2.4) is also infinite and so in most cases cannot be solved exactly. There seem to be only two techniques which have developed to solve truncated versions of (2.1) or (2.4) to provide useful, though approximate, numerical results. The technique of Kotiah [KOTI 77] is best explained by considering the transpose of system (2.1): T T T , Q n = 0 . (2.6) The Markov chain balance equation coefficients now appear row-by-row in Q . This matrix is truncated first to m rows and then to n columns, with n chosen so that no coefficients from the m balance equations are lost. In general n £ m because of the structure of these equations. With the constraint 9 n > 0 (2.7) the problem is now viewed as an (m x n) linear program and standard methods may be used to compute upper and lower bounds on quantities of interest. This technique could presumably also be applied to system (2.4). This method of Kotiah is the one used in the present study and will be discussed in more detail later (see IV - C). The second technique applies to truncations of the discrete-time linear system (2.4). Recall that any system modeled as a continuous-time Markov chain may also, with additional effort, be represented as a discrete-time Markov chain, and hence (2.4) will apply. The "(m x m) northwest corner" truncation method was first developed by Seneta [SENE 67], and was recently extended and much enhanced by Wolf [WOLF 80]. Possible algorithms using Seneta's technique are discussed in [ALLE 77]. The common approach of both Seneta and Wolf is to truncate the transition probability matrix P to m rows and m columns, , .P , and then modify the truncated matrix (m) in some special way. The solution n to the truncated system converges elementwise to JT : lim IT, = ir i=1,2,... ,m (2.8) _ ^ Cm) x i In [WOLF 80] several possible modifications of the truncated matrix , N P are given, one of which includes the technique of Cm; Seneta. For this particular version, a special kind of convergence result can be given [SENE 67,68], in which both upper and lower bounds can be found for the ratios of elements n.. / ir± (i,j =1,2,... ,m). These computations require the 10 inverse of the matrix ( (m)1 ~ (m) p) • For the other modifications to (m)P listed in [WOLF 80], even less is known about the convergence properties of the technique. For both sets of modifications it is not possible to estimate the absolute errors for individual elements: Ifm)"! ~ ffi I i = 1 ,2, . . . ,m . The advantage of the method of Kotiah over that of Seneta and Wolf is its ability to give upper and lower bounds on state probabilities no matter what the truncation m. On the other hand, Wolf's method gives only approximate state probabilities with unknown magnitudes of error. Consequently, before the result IT is judged to be "close enough" to v_ , the (m) convergence behaviour of the problem under study must be determined by obtaining solutions . sn for ever-increasing m . 11 III. MODEL SELECTION A. BASIC ASSUMPTIONS We consider queueing models of the form M|M|s with some priority allocation for two customer classes. Within classes, the customers are served on a first-come first-served basis. The s servers are identical and independent of each other. Customer traffic arrival is modeled as a Poisson process with rate X for customer type i (i=1,2). The arrivals of type 1 i and type 2 customers are assumed to be independent. The total arrival process is then also Poisson, with rate X. = X, + X.2 • The service times of each customer class are exponentially distributed with mean (?± )'1 for customer type i. The utilization factor for each class is p± = k± /(s Kj) and the total utilization is p = p, + p2 In the mathematical model we assume there is no limit on the number of customers allowed in the system. Physical systems to which this analysis applies will of course have finite queue lengths, but these are assumed to be large enough so that the probability of a customer being blocked is negligible . 12 B. SERVICE DISCIPLINES There are chiefly four service disciplines which might be considered for a multichannel, mixed voice and data, communications system. The goal is to select one which will minimize message waiting times, while maximizing use of available resources. We assume that only the distributions of service time are known beforehand, and not the exact service time of each new customer. First-come first-served (FCFS), [SLAT 73], [KOTI 73], is no doubt the simplest discipline for a multiserver system. However, it does have an element of "unfairness" for voice/data systems where the average service times for the two customer types may typically differ by a factor of 20. Data messages, which are inherently quite short, would suffer a delay-to-service-time ratio much greater than for- voice messages. Consequently the quick-response, interactive aspect of data communication would be lost. Cutoff priority, [TAYL 80], [SING 81], makes available a greater share of the system's resources to customers of one type, usually those with shortest service time. High priority customers are always admitted to service if there is a free server. However, low priority customers are allowed service only if the number of busy servers is less than some cutoff level sc . Customers denied service wait in priority-typed queues, with FCFS discipline for each queue. Cutoff priority however, is not work-conserving [KLEI 76, p.113], meaning that it is possible for low priority customers to be denied service 13 even while some servers stand idle. This is a waste of system resources, and seems an unwarranted sacrifice unless the high priority customers have an "emergency" nature. Alternating priority assumes that incoming customers are split into separate queues according to their type. In the case of two queues, the controller alternately attends each queue. A maximum of k, customers are served from queue 1, then the controller switches and a maximum of k2 customers are served from queue 2. If either k, or k2 is unbounded, the service is said to be exhaustive for that particular queue. When there are more than two customer classes, the scheme is generalized to what is called cyclic priority. Because of its flexibility, alternating priority remains a viable option for the radio system under consideration, and so will be further discussed in section III - C. Non-preemptive priority (NPP), more descriptively known as head-of-the-line priority, has a quite simple operating rule. Consider there to be only a single queue, with all high priority customers grouped together at the front, and all low priority customers queued similarly behind. Service is then always from the head of the line. If high priority is given to the customer class with the shortest mean service time, then it is known for single server NPP systems that the overall average customer waiting time is minimized [KLEI 76, p.126]. No such result has yet been proven for multiserver systems but it seems reasonable to assume that this allocation of priority would maintain its advantage for s > 1 . Later, in III - D we investigate non-14 preemptive priority in greater detail. In Appendix A we give a detailed comparison of the NPP discipline and a particular alternating priority system. C. ALTERNATING PRIORITY MODELS 1. Previous Investigations Cyclic priority and especially alternating priority queueing models with only a single server have been analyzed with some degree of success. The first work was that of Avi-Itzhak, et. al. [AVI 65] who analyzed the alternating priority model (two customer classes) with exhaustive service to each queue. In this treatment, the server switches queues without any intermediate delays, but later studies relaxed this restriction [SYKE 70], [EISE 71]. Further generalization of the model was for multiqueue systems [COOP 69,70],[EISE 72], still for the single server case. The more general cyclic service regime in which up to k, customers are served from queue 1 , then up to k2 customers from queue 2 , and so on for v queues, denoted here by C(k,,k2,...,kv) , has only recently been attempted, notably by [NAIR 76] and [SCHL 81] in which v =2 , and [KUEH 79] in which v > 2 but k± =1 for all queues. All of this previous work, however, is applicable only to the single server model, and extensions to multiple servers would seem to be extremely difficult. The analytical intractability of multiserver models arises from the unknown identities of the customers in service at any particular time. For example, in the single server case it is 15 known that the job in service at any time is from the queue the server is attending at that moment. With s > 1 servers, however, the only definite knowledge about the customers in service is that at least one must be from the queue the group of servers is attending at that particular moment. The identities, of the other customers in service are unknown. 2. Markov Chain Model — Non-exhaustive Service Consider now the alternating priority C(k,,k2) multiserver Markovian system with no switching delays. A Markov chain model for such a system has a state space with five independent dimensions, since any system state can be described completely by the following five variables: n, the total number in the system s1f the number of type 1 customers in service q,, the number in queue 1 i, an identifier of the queue currently being served c, a counter for the customers already served from the current queue in the current service cycle. For n fixed, say n = n0 ^ s , then the possible ranges for the other four state variables are: Si,: 0,1,2,... ,s q,: 0,1,2,... , (n0-s) c • 1f2f3f*«* f^ l -1f 2 • Not all combinations of these states can exist, but at least we can say that the number of states with n = n0 will be proportional to (k,+k2) . Observe that (k1+k2) is the number 16 of possible configurations of the two state variables c and i . If it were not for this factor, it will be seen (III - D) that there would be approximately the same number of states, given n = n0 , as in the 2-class non-preemptive priority models. For alternating priority models other than the very simplest, k., = 1 , k2 = 1, the number of states in the Markov chain model can be exceedingly large, and so a numerical solution may become quite costly. 3. Markov Chain Model — Exhaustive Service The exhaustive service variation may be applied to either one or both queues. The resulting numerical model turns out to be "smaller" in the sense that there are not as many possible states for a given system population n = n0 . For example, with exhaustive service to both queues, C( 00, 00 ) , a complete Markov chain representation will require a state space of only four dimensions. The corresponding state variables are n, s1f qlf and i , and for n = n0 > s , the possible ranges of the other variables are: s ^: 0,1,2,... ,s q, : 0,1,2,... , (n0-s) i : 1,2. Again not all combinations of these states exist, but we see that the factor (ki+k2) for the C(k,,k2) model has been replaced by the factor 2 (since i = 1 ,2) for this C(°°,°°) model. Thus the number of states for fixed n = n0 is about the same as if k,=t, k2=1 in the previous alternating priority 17 system C(k1,k2) . The solution of the exhaustive service system would therefore in general require much less computational effort. D. • NON-PREEMPTIVE PRIORITY MODELS 1. Previous Investigations Analytical studies of non-preemptive priority M|M|s queueing models (s > 1) have progressed somewhat further than those of cyclic priority systems. First obtained was the average waiting time for members of each priority class [COBH 54]. Unfortunately this analysis was only possible under the rather restrictive assumption that all customer classes have the same service time distribution. Most subsequent improvements to this work also required uniformity of service times [DAVI 66], [WILL 80]. Those which allow unequal service times must disallow queueing for one customer class altogether [BHAT 76], or abandon the non-preemptive discipline in favor of preemptive priorities [HALF 72], [GOLD 81], or no priorities at all (first-come, first-served). The absence of analytical 0 results for the general case of customer classes with different service time distributions again can be attributed to the unknown distribution of customer types in service at any particular time. When all customers have the same service time distribution, their identities as members of different priority groups disappear after they have been admitted to service. When service times for each class are different, however, then the identities of the customers are maintained even after they have 18 entered service, and the number of customers of each type in service will have some effect on the waiting times of those which have not yet been served. It seems that attempts at analytical solution of both cyclic priority and non-preemptive priority schemes have fallen upon the same problem, namely the heterogeneity of the customers in service at any one time. Cyclic priority models have allowed different service times for different customer types, but have restricted the system to one server. Non-preemptive priority models have restricted all customer classes to a single service time distribution, but have allowed multiple servers. It should be noted that single server NPP systems with different service times for customer classes have been solved [COBH 54]. However, the complementary problem in'cyclic priority, a multiple server model with all service time distributions identical, seems not yet to have been treated. 2. Markov Chain Model A numerical solution of the 2-class M|M|s NPP model requires less computational effort than does the alternating priority model. The Markov chain state space requires only three independent dimensions: n, s,, q, . Compared with the alternating priority model C(k1fk2) , we see that the state variables c and i have been dropped. Consequently the number of states for a given n = n0 > s , is essentially half of what it would be for the simplest alternating priority systems C(1,1) , or. CC00,00). This advantage of simplicity is 19 one reason why the NPP model was selected for closer analysis. 20 IV. NUMERICAL SOLUTION TECHNIQUE A. THE STATE SPACE As mentioned in the previous chapter, the non-preemptive priority model requires a state space of three independent dimensions. For convenience a state will be labelled by four components (s, s2; q, q2) , where s1f s2 are the number of customers of types 1 and 2 in service, and q,, q2 are the number of customers of types 1 and 2 waiting for service. The central controller, the device which selects a waiting customer for service whenever a departure occurs, does not need to be represented by a separate state variable. Its behaviour is totally determined by the given state (s, s2; q, q2) (see Appendix A). It is convenient with two customer classes to think of there being only one queue (pooled storage), with type 1 customers at the front and type 2's at the back. The total number of customers in the system is n = s,+s2+q1+q2 . As this model assumes unbounded queue lengths, the dimensions q, and q2 (or n and q,) are also unbounded. When n < s (s= number of servers) then: 0 < s± < n (i = 1 ,2) , q,= 0, q2= 0 and s,+s2 = n . When n > s , then: 0 < s± < s (i=1,2), q,,q2 > 0 and s,+s2 = s , Qi+Q2 = n - s . One choice of three independent dimensions for this state space is (n,s,,q,) . In what follows we choose rather to use the 21 more convenient state vector (s, s2; qi q2) • B. THE BALANCE EQUATIONS The steady state balance equations are obtained in continuous-time Markov chain models by equating the transition (or flow) rate into a particular state with the transition rate out of the state. This prescription is not always easy to follow for systems of dimensionality greater than one because obtaining the total flow rate into a given state from its neighbours may not be a trivial matter. An easier but more mechanical approach to the construction of the balance equations is given below. Let P(s! s2; qi q2) be the stationary probability that the system is in the state (s, s2; q, q2) . These probabilities for states in a three-dimensional state space are ordered by a single indexing variable i (i=1,2,3,... ) to form the infinite column vector of probabilities x^ . Adapting equation (2.1) then: xT Q = 0T (4.1) — CO — with x and 0 both column vectors. Given that the system is — CO — currently in state i then: Q±± = transition rate out of state i (negative) Q±. = transition rate from state i to state j Note that by construction, all rows of Q sum to zero. The balance equation coefficients appear along the columns of Q . All that is now necessary is to generate the Q ±j coefficients (usually row-by-row) and then from the structure of the columns 22 of Q , deduce the general form of the balance equations. The service discipline for non-preemptive priority may be stated as follows (assume n>s): At the next departure from the system select for service a type 1 customer, unless qi=0 , in which case select a type 2 customer. This service rule governs the allowable state transitions due to departures from the system. Combining these with the transitions caused by system arrivals, a state space transition diagram can be constructed (see Appendix B). Once the service rules for a queueing system have been clearly expressed it is not a difficult matter to generate the Q.. conditional transition rates. Recall that i is the index of the system's current state, written i= (s, s2; qi q2) > with number of customers n = s1+s2+q1+q2 • The object state j will have population (n ± 1). Then: On IT + X.2 +S 1 H 1 + s 2 * 2 ) n > 0 r I (s, + 1 s2 • 9 0 0 ) n < s I p j= (s, s2 m 9 q1+i q2) n > s X-2 I (s, S2 + 1 m f 0 0 ) n < s < (s, s2 • f q, q2 + D n > s S 1 K ! * j= (s, -1 s2 9 9 0 0 ) n < s j= (s, -1 S2 + 1 m t 0 q2-D n > s .j= (s, s2 • • r qi-1 q2) n > s, 'j= (s, s2-1 • 0 0 ) n < s j = (s, + 1 s2-1 • 9 q1-1 q2) n > s, qi* j= (s, s2 • 9 0 q2-D n > s (4.2) These transition rates are a complete list of all allowable transitions indicated by the state space transition diagram in Appendix B. 23 Generalizing the structure of the columns of Q , we can now obtain the steady state balance equations. (x,n2+s,(.,+s2«2) P(S, s2) X, P(s,-1 s2) + X2 P(s, s2-1 ) + (S , + 1 ) * , P(s , + 1 s2) + (S2+1) n2 P(s, s2+1) 0 < n < s-1 (4.3) (X ,+k2 + s, u i+s2>.2) P(s, s2) n = s = X, P(s,-1 s2) + X2 P(s, s2-1) + s,m P(s, s2; 1 0) + (s, +1 ) »», P(s, + 1 s2-1; 0 1) + (s2 +1) u2 P(s,-1 s2 +1; 1 0) + S2K2 Pts, s2; 0 1) (4.4) (X. ,+X2+s , u , + s2u2) P(s, s2; q, q2) n > s+1 X, P(s, s2; q,-1 q2) + X2 P(s, s2; q, q2"D + s , (•, P ( s , s2; q , +1 q2) + 6(q,) (s, + l) 1*1 P(s, + 1 S2-1; 0 q2 + D + (S2 + 1) t.2 P(s,-1 s2+1; q, + 1 q2) + 6(q,) s2„2 P(s, s2; 0 q2 + D (4.5) Here we have used the following conventions: P(s, s2; q, q2) = 0 , when any index is negative, P(s, s2) = P(s, s2; 0 0), 6(q,) = 0 q,= 0 i \0 otherwise . In standard matrix form, the balance equations are written as the transpose of (4.1): QT x = 0 (4.6) The successive rows of Q represent the balance equations for the states i= (s, s2; q, q2) , ordered first with increasing n , then with increasing q, , and finally with increasing s2 . With this ordering of states, the n > s part of matrix T Q possesses only a few different (s+1 x s+1) recurring blocks of coefficients, a fact which simplifies the automatic generation of the balance equations. Other orderings with index q2 replacing q, , or s, replacing s2 would be equally convenient. Refer to Appendix C(1) for further comments and an 24 example. A finite system of balance equations is obtained by including equations of the form (4.5) with those.of (4.3) and (4.4) until all those states for which n = n , a cutoff parameter, have been represented on the left-hand side of (4.5) . The right-hand sides of these equations contain states with n= 0 up to n= n +1 . Consequently, viewed as a system of homogeneous linear equations, the nc-truncated model would have a rectangular coefficient matrix A, with more columns than rows. The dimensions of this matrix are derived in Appendix C(2). C. THE LINEAR PROGRAMMING SOLUTION TECHNIQUE As mentioned in the discussion of Markov modeling, one of the techniques for solving an infinite set of balance equations is the linear programming solution of [KOTI 77]. From (4.6), the infinite linear system of balance equations is written: QT x = 0 (4.6) — oo — with normalization: I X-L = 1 . (4.7) i=1 T The finite matrix A is obtained by selecting from Q the first Nj_ rows and Nj columns. In general Nj i Nj , and the smaller the difference (NJ-NJ- ) , the better will be the results produced (see Appendix C(2)). The minimum value of ^NJ_NI^ for a given NJf is determined both by the structure of the Markov chain model, and by the ordering of the corresponding balance equations. 25 Define x as the finite column vector of the steady state probabilities (xi; i=1,2,... /N^}. This truncation of x^ means that the normalization (4.7) can no longer be applied. However, a weaker but broadly applicable condition, derived from steady state arguments and (4.7), can be used instead (see IV -D): B x = b (4.8) with B a row vector of constant non-negative coefficients, b a scalar positive constant. The truncated problem may now be cast as a linear program (LP): Minimize and maximize each {Z, ; k=1,2,... K] subject to the constraints: A x = 6 B x = b (4.9) x > 0 . The objective functions (Zk} may be chosen as any linear combinations of the state probabilities x^ which are of interest. For example, if mean waiting times are desired then these may be expressed in terms of the xi using Little's Law (see V - C). To obtain both the maximum and minimum of each waiting time the LP (4.9) must altogether be solved four times. Another application is to let the Z^ be components of a marginal probability distribution, Zk = P(k) , k=1,2,. . . ,K in which case each P(k) is a linear combination of the state probabilities x± (see V - B). The LP (4.9) must then be solved 2K times to obtain the maximum and minimum values of the probabilities P(k). These numerical values form the set of upper bounds {PuCO}, and the set of lower bounds (P^(k)}, (k=1,2,... ,K), for the true steady state probability 26 distribution P(k) of the system. That is, Pa (k) < P(k) < Pu(k) , k=1,2,... ,K . (4.10) Note that neither the components P£^) > nor the components Pu(k) form a probability distribution themselves, since in general neither sum to unity. D. NORMALIZATION CONSTRAINTS We now consider the task of choosing appropriate constraints which will ensure that any solution of the linear program (4.9) will be consistent with the normalization (4.7): oo I x. = 1 . (4.7) i=1 1 Equivalently, for (4.7) we write co E p(n) = 1 (4.11) n = 0 where p(n) is the total probability of n customers in the system (n>0). These conditions, collectively called normalization constraints, are represented as a generalization of (4.8): B x = b (4.12) with B the coefficient matrix, of dimension (^ x Nj), b a constant column vector (^ x 1). These NN constraints are assumed to be non-homogeneous equalities involving at most only the first Nj components of x . — 00 Recall that the truncation of the infinite system of balance equations means that the original normalization (4.7) can no longer be exactly represented as a constraint in an LP. In its place, Kotiah [KOTI 77] presents without proof the 27 following: s-1 E (1 - n/s) p(n) = 1 - p . (4.13) n=0 In Appendix D we derive (4.13) under the assumption that a stationary probability distribution for number in the system p(n) exists, with Er p(n) = 1 (4.11). In other words, the system is assumed to reach a steady state. Kotiah suggests that (4.13) and (4.11) are equivalent for a large (unspecified) class of queueing systems. However, we know only that (4.11) implies (4.13); the reverse need not necessarily be true. Consider the multiserver non-preemptive priority queue. It is possible to solve the truncated Markov chain model for p(n) which satisfy (4.13), but which sum to a value much greater that one. Clearly such a solution violates (4.11) because physically all values of p(n) must be non-negative regardless of whether they are included in the truncated Markov chain model (0 ^ n < n +1) or are not (n > nc+l). Thus for the queueing system considered here, equations (4.11) and (4.13) are not equivalent. We observe that none of the four examples in.Kotiah's original paper [KOTI 77] seem to give rise to similar problems. To preclude the possibility of obtaining LP solutions which violate (4.11), an additional constraint must be included: NJ E x± < 1 (4.14) i = 1 or equivalently, n +1 T p(n) < 1 . (4.15) n=0 Modifying Kotiah's original LP solution technique to use both normalization conditions (4.13) and (4.15) instead of just 28 (4.13) alone, then should allow its extension to models which are not so inherently "well-behaved" as those considered in [KOTI 77]. To fit the new constraint (4.14) into the form (4.12) a new variable y , called a "slack variable", is introduced such that: NJ E x + y = 1 i=1 (4.16) y £ 0 This variable is appended to the end of the truncated probability vector x , which now will have (N^+1) components. The column dimension of coefficient matrices A and B must similarly be increased by one. E. IRREDUCIBILITY OF THE MARKOV CHAIN The complete modified linear programming problem may finally be written: Maximize and minimize Zy. , a linear function of x, subject to x = |^0j (4.17) x > 0 where, A is the (Nxx Nj+1) coefficient matrix of the NI balance equations, B is the (%x Nj+1) coefficient matrix of the NN normalization constraints. In the present analysis NJJ = 2, but N and Nj depend on the selected truncation ^ (see Appendix C(2)). The standard LP solution method is the simplex algorithm [LUEN 73], in which a special kind of search is made to determine the optimum of the objective function Zk . However, we will show that for many Markov chain queueing models the simplex algorithm need not be 29 applied to the entire problem (4.17), but only to a small part of it. This results in a substantial savings in computer solution time. 1. The Significance of Irreducibi1ity The notion of irreducibility is a well-known characteristic of some Markov chains. A Markov chain is said to be irreducible if and only if every state can be reached from every other state in a finite number of steps [BHAR 60, p.14]. More mathematically, an irreducible Markov chain is one in which the discrete-time transition probability matrix P cannot be written: P = | P, 0 P2 • Pi 0 Because of the way any continuous-time Markov chain can be represented as a discrete-time chain (see II - B), the previous definition holds with P replaced by Q . Many Markov models of queueing systems are indeed irreducible Markov chains. The present NPP system is no exception, as can easily be seen by the structure of the balance equations (Appendix C(3)). Consider now the N balance equation constraints in the LP (4.17). Since any state i may be reached from any other state j , it follows that if the stationary probability of state i is zero, x.= 0 , then so also must there be zero i probability of the system being in state j : x^= 0 . Continuing this reasoning and using the property of irreducibility, we find x.. = 0 , j = 1,2,... ,H1 . Therefore we conclude that if an irreducible Markov chain is to have a non-null stationary solution, then not one of its states may have a 30 zero steady state probability. This observation, which was indirectly used in [KOTI 77], forms the kernel of our improvements to Kotiah's LP solution method. 2. Improvements to the LP Algorithm Suppose we apply the simplex algorithm to an arbitrary LP problem with m constraints and n variables (n > m): U x = u (4.18) A basic solution x^ is defined as the solution of (4.18) when a certain (n-m) -component subset of {xi} is set to zero. The remaining m nonzero variables are called basic variables, and together they form the basis set. There are at most f n] basic solutions. The simplex algorithm in effect "searches" those basic solutions which are also feasible (x^ > 0) for the one which optimizes the objective function Z [LUEN 73]. Now we can invoke the irreducibility property to show that the number of basic feasible solutions to (4.17) is much less that the general maximum (jj^J • Irreducibility guarantees that if any one of the variables x.^ , 1 < i < N-j- , is set to zero, then all the other variables x^ , 1 < j * i < N-j- , must also be zero. This would violate the normalization constraint (4.13) and so no basic solution to (4.17) can exist if any of the first Nj variables are excluded from the basis set. The basis set must then contain all of the first Nj. variables. The set is completed by choosing (m-N-j. ) further components of xb out of the remaining (n-N^. ) candidates. The number of possible basic solutions has thus been reduced from fn\ to /n-Nj m J \m-Nj 31 From (4.17) we find n = Nj+1 , m = Nj+2 , so that: n - Nj> = fNj+l-Nj^ = (NJ+1-NJ.) (NJ-NJ) /2 (4.19) This is an extraordinary reduction in the number of solutions which must be searched by the simplex algorithm. As a consequence, what would have been a prohibitively expensive numerical analysis has now been made affordable. F. IMPLEMENTATION OF THE MODIFIED LP ALGORITHM A special solution program was written for the LP (4.17) to take full advantage of the computational savings gained by recognizing the Markov chain's irreducibility. The centerpiece of the program is the so-called "revised simplex method" [MURT 81, p.24 ff.], the standard computer algorithm used to solve general LP problems. However, this method is applied only to determine the identities of the last two variables in the optimum basis set; the first Rj. variables in any basis are already known. In practice this means that the steps of the simplex method are applied only to the .(N -N +1) "free" columns of £B] ' "^REE" ^N the sense that they are not known a priori to be a necessary part of all basic solutions. The algorithm works with the LU-factorization of the basis matrix. The basis matrix is the (m x m) square matrix formed from those columns of ^gj which correspond to variables in the current basis set. As each new basic solution is obtained during the search for an optimum, the LU-factorization must be updated. For this purpose we use the scheme of Bartels-Golub [MURT 81]. The solution program requires the user to specify two free 32 variables which will complete the first basis set, ie. yield the first basic feasible solution. We have chosen the variables (Nj+1) (the slack variable y), and (Nj-4) . This choice has a significance which will be explained in the next chapter (see V - E). Besides the property of irreducibility, a perhaps equally important practical aspect of the LP (4.17) is the extreme sparsity of the coefficient matrix . In fact this matrix can be classified as super-sparse because the number of distinct coefficients is very small compared with the total number of nonzero entries. Such a matrix can be very efficiently stored in compacted form [KALA 71]. The LU-factorization and the manipulations of the simplex method will preserve some degree of sparsity (though not super-sparsity). Consequently, standard sparse matrix data structures [GUST 72] can be used throughout the computations. The LU-factorization of the initial basis matrix is based on a standard Gaussian elimination algorithm, modified to take advantage of sparsity. We have used an initial step of "symbolic" factorization (a simplified adaptation of Sherman's NSPIV program [SHER 78a,b]) to determine the changed sparsity sructure of the factored basis matrix. Then with this information the purely numerical steps of factorization are finally performed, using an algorithm given in [GUST 72]. In the factorization step we find again some useful peculiarities in the structure of the LP (4.17). In particular the (N-j-x Nj) submatrix of A , call it A, , (which forms by 33 far the largest part of every basis matrix) exhibits diagonal dominance by columns. Therefore an LU-factorization which proceeds along the diagonal elements (no pivoting) is guaranteed to be numerically stable [REID 74]. The algorithm then need not employ any complicated pivoting strategy nor require the facility for dynamic row and column interchanges. Furthermore, LU-factorization will be stable on any symmetric permutation of the rows and columns of A, , that is, on II A, nT , with n the permutation matrix. The permutation n may then be selected solely to minimize the creation of nonzero elements in locations where there were only zeros. Finding a row and column ordering to reduce this "fill-in" is a somewhat heuristic process, but it is usually a worthwhile search since a good ordering can greatly reduce execution time and memory costs. An efficient permutation was found to be one in which the row and column order of A, is simply reversed: the last row is placed first, the second-last row is placed second, and so on; likewise for the columns. We found experimentally for At of order = 1000 that this row and column permutation results in a very large reduction in execution time, at least by a factor of ten. All of the modifications discussed in this section have been incorporated into the LP solution algorithm. To test the operation of the final program, we applied it to the solution of the M|M|2 FCFS queue. The program successfully duplicated Kotiah's published results [KOTI 77] for the system. In the next chapter we apply these numerical techniques to solve the M|M|4 NPP 2-class queueing model. The effectiveness and 34 efficiency of the modified algorithm will then be assessed. 35 V. RESULTS A. SYSTEM PARAMETERS We now present numerical results for the s = 4 server case of the multiserver 2-class NPP queueing system. To make the correspondence between this queueing model and the communications system introduced in chapter I we identify type 1 and type 2 customers as data and voice messages respectively. The mean service times of these two customer types are ((,,)•' = 5 sec, and (it 2)" 1 = 120 sec. The mean arrival rates X, and X2 are not directly specified. Rather, for each numerical solution we fix the traffic intensity p and the proportion of type 1 arrivals to the system, a,. Recall from (III - A) that X = X,+X2. Then define o±= X±/ X (i = 1 ,2) , (5.1) so that in general a± is the fraction of incoming customers which are of type i (oi+o2=1)• All mean arrival rates are thereby defined since: (from Appendix D) p - X. / O_L + a_2_ "\ . (5.2) S \ V 1 t>2 ) The rate parameters X,, X2, mi fit all of which have dimension (time)"1, may be measured in any convenient unit of time; minutes instead of seconds for example. Because type 2 customers have on average 24 times the service requirement of type 1 customers, the proportions o,, o2, of incoming customers do not directly reflect the division of system resources between the two types. For example, in order to have the system's total service time allotted equally 36 there must be 24 new type 1 arrivals for each new type 2 customer. We formalize this notion by defining p± (i=1,2) as the fraction of system resources ("work") spent serving customers of type i (£i+02=1)« *i= P±/ P (i = 1,2). (5.3) Identify />i as the work expended by the system to serve class i customers (from Appendix D), P±= X.±/ (s „±) (i = 1 ,2). (5.4) Using (5.1 - 5.4) we can show: », = [ 1 + (o2M )/(a,n2) ]" 1 . (5.5) Then, for the above example, = 0.5 when d = 24/25 = 0.96 . Figure 1 shows the relationship between £, and o, for various ratios »-i/i>2 . Each curve begins at ay = ty=0 with slope H2/M1 and ends at c^p^l with slope n^/r<2 . In between, the value of a, at which the slope is unity, o,*, is given by o,* = (M/>2)'/2 . (5.6) 1 + (Y,/^)1'2 For the current system nl/f2 = 24, and o,* = 0.83 . Referring to Figure 1 we see clearly that there is a wide region over which a large change in o, produces only a small change in Only when a« > a,* does the division of work within the system (as measured by £,) become sensitive to changes in a, . However, while the systems analyst is more interested in p1f the user of a communications system is more familiar with o,. In what follows then, both parameters will be specified. 37 Figure 1 - Relationship Between Parameters and o, 38 B. THE STEADY STATE PROBABILITY DISTRIBUTION p(n) We first determine upper and lower bounds for each of the objective functions {Zk} corresponding to the distribution of number in the system: p(n), n=0,1,2,... ,nc + 1. Recall that with the truncation n= n for the infinite set of balance c equations, the LP (4.17) includes states x± for which n= 0 up to n= n£+1 (see IV - B). In terms of the individual state probabilities P(s, s2; qi q2) we define: n E P(s1s2;00) 0 < n < s s,=0 (s2= n-s,) p(n) = (5.7) n-s s E E P(s, s2; q, q2) n > s q,=0 s,=0 (s2= s-s,) (q2= n-s-q,) Figure 2 shows the upper and lower bounds obtained on p(n) for the case />=0.7, o, = 0.B0, £,=0.14, nc=25 . The plot is on a logarithmic scale and for clarity the discrete points are connected by continuous curves, though strictly a histogram representation should be used. The most obvious characteristic of this curve is the upward flare on the tail end of the upper1 bound. This feature seems to be present in some degree no matter what set of parameters (p,c,) are chosen. Neither is it specific to this particular queueing system. A more detailed analysis than that given in [KOTI 77] of the M|M|2 FCFS 2-class system reveals the same behaviour. Such a characteristic arises from the way in which (p^(n)} and (pu(n)} are computed. There are no constraints imposed on the sums 39 Figure 2 - The Distribution p(n) for p = 0.7, o,= 0.80, 0.14) 40 nc +1 nc +1 I P>) , Z pu(n) , n=0 n=0 because each maximum, pu(n0), or each minimum, p^(n0), is the optimum of a unique linear programming problem (see IV - C). Consequently, the upward flare in {pu(n)} is only the result of a general tendency for the relative difference between the bounds (ie. (pu(n)-p^(n))/ P^(n) ) to increase as n increases. In some cases the upper and lower bounds for p(n0) where n0 is near to (nc+1) can be extremely poor indeed. In general, however, if nc is increased with n0 held constant, the bounds for p(n0) will improve. More will be said later about this concept of convergence (see V - F). Figures 3 and 4 are further examples of the upper and lower bounds obtained for p(n). Both of these cases are for />=0.5 and nc=20, while Figure 3 has a,=0.80 ($,=0.14), and Figure 4 has a, = 0.95 (»31=0.44). We observe that for a fixed ci, generally tighter bounds are obtained for less cost (smaller n) when p is decreased; compare Figures 2 and 3. For a fixed n and p, tighter bounds are obtained when o, (or B,) is made smaller; compare Figures 3 and 4. The general character of the bounds obtained for p(n) is not inconsistent with the possibility that p(n) varies geometrically with n for n "large". Many queueing systems show this behaviour, even though the exact expression for p(n) may be more complicated. An example is the M|M|s FCFS queue, for which p(n)a pn when n ^ s. 41 ]0 15 2(T NUMBER IN SYSTEI Figure 3 - The Distribution p(n) for p = 0.5, o,= 0.80, (/»,= 0.14) 42 io-i in-I >-CQo CC-" CD O cr •_ LO-I p= 0.5 a ,=0.950 =0.442 n6=2Q —i 1 1 10 15 20 NUMBER IN SYSTEM 25 30 Figure 4 - The Distribution p(n) for p = 0.5, o,= 0.95, 0.44) ' 43 C. MEAN WAITING TIMES 1. Def in it ion An important performance measure of any queueing system is the time a customer of a given class must wait in the queue before being served. The entire continuous distribution of waiting times is of interest, but the present LP method is best suited to obtain bounds on the average waiting time. Denoting the mean waiting time of class i customers by W±, we have from Little's Law [KLEI 75, p.17]: W. = E(q. ) / \. (i = 1 ,2) , (5.8) where E(qi) is the average number of type i customers in the queue. In terms of the state probabilities P(s, s2; qi q2) then: W. = 1 E 1 \± n=s+1 n-s s I I q± P(s, s2; q, q2) q±=0 s,=0 (5.9) with qi+q2 = n-s, s,+s2 = s, and i=1,2 . There is one very important difference between the expressions (5.9) and (5.7) for Wi and p(n). Each p(n) can be represented exactly as the sum of a finite number of state probabilities, whereas the exact representation of Wi is an infinite series. This same problem occurred in (IV - D) when the implementation of the normalization (4.11) was being considered. Here we apply the same remedy. Truncate the infinite sum to n= nc+1 , and define the objective function Z(Wi) to be the expression given in (5.9) with the infinite upper limit replaced by (n +1). Then, 44 W. > Z(W.) (i=1,2) . (5.10) 1 i The LP method is used to find the maximum and minimum values of the right-hand, side of (5.10). The minimum will serve as a lower bound for W.^ although it will not be as tight as it could be (if nc were greater). On the other hand the maximum may in fact not be an upper bound for Vl± at all, depending on the relative magnitude of n,. We must be careful to choose nc large enough so that Z(Wi) becomes a very close approximation to W (5.9). Then, for practical purposes the LP solutions will provide true upper and lower bounds for W±. The problem of choosing a "large enough" nc will be discussed later in the section on convergence (V - F). 2. The Limiting Cases: o1 -> 0, o 1 In this section we determine the limiting values for both W, and W2 as 0,-^0, and 0,-^-1 . It is easy to determine W^OT-M) and W2(o1-»0) because in these cases the system is very nearly an M|M|s homogeneous type 1 or type 2 system respectively. Then from [GROS 74 p.98]: vMc.-M) = Po (sp. )s s„ , s! W2(oi-»0) = where, Po = S(i2 \1~P2 k r s-1 I (S/>) k=0 k! (Sp)S (1-P) s! -'1 (5.11) (5.12) (5.13) To calculate W^d-^O), imagine a single type 1 customer entering an M|M|s homogeneous type 2 system. No matter how many type 2 customers are in the system, this type 1 customer will, 45 because it Was priority, go directly to the head of the line, and is guaranteed the next free server. If this customer's waiting time is denoted by w1f then w, > 0 iff n > s, and: E(w,|n) = J 0 n < s | (sx2)"1 n £ s . Therefore, E(w,) = E[ E(w,|n)] = E(w1| n>s) Pr(n>s) . We define VMd-^O) = E(w,). From [GROS 74, p.96] we obtain Pr(n>s), and then: W^oi-^O) = pp (s/>2)s . (5.14) Su2 (1 -P2) s! The calculation of W2(o1->1) is somewhat more difficult. Consider a single type 2 customer entering an M|M|s homogeneous type 1 system. This customer does not have priority, and in fact must wait until the queue entirely empties of type 1 customers before it has a chance to occupy the next free server. Let this waiting time be denoted by w2 . This situation is exactly the same as that discussed by Cooper [COOP 72, pp.220-223] for the behaviour of a "polite" customer in an M|M|s homogeneous queue. A polite customer is one which declines service as long as there are customers waiting. Cooper shows, using an analysis of the M|M|s busy period, that E(w2| w2>0) = 1 . (5.15) (,-P,)2 SUy Since w2 > 0 iff n £ s, then: E(w2|n) = J 0 n < s [[(1-p,)2 s,.,]-1 n > s . Therefore, E(w2) = E[ E(w2|n)] = E(w2| n>s) Pr(n>s) . We define W2(c,-»1) = E(w2). 46 Again using Pr(n>s) we have: W2(o1->1) = Po (sp,)s . (5.16) s»i, (1 - p ! ) i s! There are some interesting relationships between the four limiting values of Wi . First observe from (5.2) that for fixed p : /i,(or>1) = p , and />2(oi~>0) = fi • Making these substitutions, we then define R, and R2 as: R, = W2(oi-» 0) = W2(c,-» 1) = 1 (5.17) W, (0l-» 0) W,(o,-» 1 ) ( 1-/)) R2 = W,(o,-» 1 ) = W,Ui^ 1 ) = M2 (5.18) Wj(o,-» 0) W,(o,-» 0) u , ( It can be shown for these two ratios that: R, > 1 for 0 < p < 1 (5.19) R2 < 1 for 0 < p < (1 - u2/u,) The four extreme values of waiting time, Wi (a,—>-0), Wi(a1->1), i=1,2, when plotted versus a, (or £}) form two parallel lines. These lines have slope (R2-1) W^CH-^O), and vertical separation (R,-1) W,(a,-»0). 3. Waiting Time Results Now we present the computed upper and lower bounds on waiting time, plotted as a function of the parameter p,, for the following three cases: Figure 5 p = 0.7 nc = 30 6 0.5 25 7 0.3 20 . The. results from simulation studies of this queueing system are also shown. These will be discussed later in section D. The values of W^^ at the extreme points Pi=0, = are exact calculations using the techniques of section C(2). All other 47 48 Figure 6 - Waiting Times as a Function of 0, : p = 0.5, n = 25 c 49 Figure 7 - Waiting Times as a Function of fi\ : P = 0.3 n = 20 c i: 50 I values of W± have been obtained through numerical solution of the LP (4.17). Numerical solutions were computed for p,= 0.01 and 0!= 0.99 to illustrate the degree to which they approximate the known limiting cases of fi,=0, pi = 1 . As can be seen from the graphs, the only significant failure of this approximation is the upperbound of W2 as fi,—>1 , for p > 0.5 . In fact for p= 0.7 the upper bound of W2 for £,= 0.99 is so far from the known limiting value W2(o1—M) that a numerical solution with a much higher nc would be necessary to obtain a more useful upper bound. To generate the data in Figures 5-7, the parameter nc was selected to ensure that the upper and lower bounds for each of W, and W2 were reasonably close together in the region 0 < < 0.2 . This region is of much physical interest because it corresponds to a very wide range of a, (the proportion of data messages): 0 < o, ^ 0.86 . In addition, nc had to be chosen sufficiently large to ensure that the error due to truncation of the exact expression for W.^ (5.9) could reasonably be neglected (at least in the given region). These waiting time results also show what was apparent by the p(n) graphs, namely that tighter bounds are obtained at less cost (smaller nc) with lower values of p and/or For fixed p and fit it is generally true that better numerical bounds will be obtained only by increasing the number of equations (higher nc) in the LP (4.17). Greater convergence of the bounds is thus directly linked with the computational cost of the numerical solution. These topics will be discussed 51 later in sections F and G . A more extensive LP solution could be a costly option, and may not even be required if some faith can be placed in two very helpful empirical observations. This discussion will, however, be deferred until after the simulation results have been presented. D. WAITING TIME SIMULATIONS Computer simulations of the M|M|4 2-class NPP system were performed for two reasons: (1) To verify the validity of the numerically obtained bounds on the mean waiting times. (2) To compare the relative cost and quality of results obtained with simulations as opposed to numerical solutions using the LP method. The simulation results are indicated in Figures 5-7 by 95% confidence bands around each "measured" mean waiting time. A c% confidence interval is calculated from the statistical variation of the simulation data such that on average the given range will contain the true value of the measured quantity c times out of every 100 such simulations. From among the several methods for statistical analysis of simulation results now available we have chosen the "successive independent batch means" technique [FISH 78]. A single simulation run is made, in which a very large number of customers are allowed to pass through the system. This very long run is divided into a certain number of batches of equal length (ie. each with the same number of customers). From each batch the mean of the quantity of interest is then obtained. By making the batches long enough, successive batch means become 52 statistically independent. The degree of independence is tested using the Ck statistic [FISH 78, p.238], a quantity closely related to the correlation between neighbouring batch means. Then, using standard techniques, 95% confidence levels are obtained on the overall mean of the measured quantity. The initial state of every simulation run is a completely empty system. Therefore to prevent any bias which might be introduced into the final results by the transient start-up behaviour, the results of the first batch are discarded. The simulation program itself was written in the specialized GPSS language [BOBI 76], and compiled by the optimizing compiler GPSS/H. Each simulation run simultaneously produced data for W, and W2. Therefore the data in Figures 5-7 are the results of 12 different simulation runs, each run divided into 20 batches of 4000 customers. It is usually easier, in terms of programming effort, to simulate a queueing system with the aid of a specialized simulation language than it is to set up and solve the Markov chain model of the system. However, whatever advantage is gained in programming effort could be lost when program running times are considered. Generally a large number of events (customer waiting times) must be observed before acceptable confidence intervals for the means can be obtained. As the results indicate, the confidence intervals are not everywhere as narrow as a systems analyst might like. If for certain values of p and 0! better results are desired, then greater costs must be borne; the width of the confidence intervals will on 53 average decrease only as N~1/2f where N is the total number of observations. E. COMMENTS ON THE WAITING TIME RESULTS There are several observations, both empirical and theoretical, which can be made about the character of the waiting time graphs in Figures 5, 6, and 7. 1. Empirical Observations Combining simulation data with the numerically obtained bounds for waiting times, we may conclude with reasonable confidence that the computed lower bounds are very near the true steady state values. This empirical fact is especially useful for W2, in which case the upper bound is often too large to be of practical use. Furthermore, a second empirical observation may be used to virtually eliminate much of the work performed by the LP solution algorithm. We have observed for LP solutions of this particular Markov chain model, that the initial basic feasible solution specified in (IV - F) generates the minimum values for W,, W2, and most of the p(n) objective functions. This has been observed to be generally true for p > 0.5 , and 0, £ 0.1 , a region where acceptable and complete LP solutions become ever more expensive. These two observations seem to allow a convenient and inexpensive numerical approximation to be made for the mean waiting times W±: Solve the linear system (4.17) with the initial basis set given in (IV - F). Calculate W. and assume 54 that these quantities are first of all lower bounds, and in addition, are reasonably close to the unknown true values. This procedure, based as it is solely on empirical evidence, is constrained to apply only to this particular Markov chain queueing model. No theoretical support has been found which would allow its extension to other cases. 2. Waiting Times in a Non-preemptive Priority System The general character of the VJ± vs. results may not at first be expected. As B ^ increases more type 1 customers enter the system while the decrease in the type 2 arrival rate is relatively slight, so why does W, decrease? The answer lies not only in the particular assignment of priority, but also in the general nature of queueing systems with heterogeneous customer classes. Because type 1 customers have non-preemptive priority, their mean waiting time is less than that of the type 2 customers, for any given p and £, . However, the variation of W^^ with 0, shown in Figures 5-7 has more in common with the M|M|s FCFS queue than with a priority queue. In general, waiting times depend on both p and u , increasing with p and decreasing with u . They cannot be expressed solely in terms of p , as is usually the case for the mean number in the system. In an M|M|s FCFS queue with the same two customer classes, the average waiting time would range from that of a homogeneous type 2 M|M|s queue at B,= 0 , down to that of a type 1 queue at p,= 1 . Although p is constant, the effective average service rate continuously increases from s>»2 55 to SMI as 0, increases. The waiting time at 0,= 1 is »ii/c2 = 24 times lower than the value at p, = 0 . This can be seen by comparing W^d-M) and V)2(a,->0) , from equations (5.11, 5.12) in section C(2) . Therefore, even though type 1 customers have priority, they spend more time waiting for type 2 customers to leave service than they do waiting to reach the head of the line. This behaviour is a direct result of assigning priority to those customers with shortest average service times (greatest i>). In this situation, however, the non-preemptive priority service discipline takes on an ever-increasing importance as the traffic intensity p increases. At some point, say p= p* , the type 1 customers will begin on average to spend more time just getting to the head of the line than they actually spend there waiting for a free server. Consequently, the average wait of the priority customers depends more on customers of their own type than on the non-priority customers. Since for fixed p, X., increases with p, , - X, = *, p Sn, (5.20) we may then expect that for p > p* , the waiting times of both types of customers will increase as i1 increases. For p < p* the system is dominated by type 2 customers, those with greatest mean service time. For p > p* , the system is dominated by the priority customers. Although strictly unproven in this context, a good assumption for p* is obtained from section C(2) by setting R2= 1 (5.18): ie. W (c.-M) = W (d-s-0) (i=1,2). (5.21) 56 Both equations (i=1,2) give the same solution, P* = 1 ~ VZ/Ki (n i>e2)• (5.22) For the current model p* = 0.96, which is too high to be of practical significance. Nevertheless it does serve to illustrate how extensively the type 2 customers, with their very long service times, dominate the behaviour of this queueing system. A. CONVERGENCE OF UPPER AND LOWER BOUNDS If the LP solution technique is to have an advantage over simulation as a method of systems analysis, it must be able to provide comparable or "better" results at lower overall cost. A given computed result is seen as "better" if the difference between the upper and lower bounds is smaller than the 95% confidence range obtained from simulation data. In this section we look more closely at the quality of results produced using the LP method, while questions of cost will be discussed"in the next section. The linear programming solution technique guarantees a non-increasing sequence of upper bounds and a non-decreasing sequence of lower bounds as the LP is solved with more and more constraints. Adding more balance equations brings the LP (4.17) closer in form to the infinite linear system (4.6) with the exact solution x^. Consequently the bounds on a particular objective function can be expected to converge steadily toward the true steady state value as nc increases. To this argument we add the proviso that the objective functions must be able to 57 be expressed using only a finite subset of the state probabilities {x±; 0 £ i £ Nj}. Into this class fall the probabilities (p(n); Q£ n < (nc+1)}. A second class of objective functions are those which can only be defined over an infinite subset of the state probabilities. Two examples of interest to us are the mean waiting times W, and W2. As typical examples of convergence behaviour for the first class of objective functions, we have plotted in Figure 8 the computed bounds for p(0) and p(4) as the LP is solved with ever greater n^ . To help interpret this graph, recall that the relation between n^ and m , the number of equations in the LP, is given from (4.17) and (C.2) as: m = Nx + 2 m = [10 + 5(n -3)(n -2)/2. + (n -3)] + 2 . ' (5.23) c c c It is also helpful to know that the overall computational cost of the LP solution is roughly proportional to (see V - G). To show the convergence of bounds for objective functions of the second type, we present Figures 9, 10, and 11. There, the computed bounds for W, and W2 are plotted as a function of nc , for three different sets of system parameters (/>,<»,). Because objective functions of this class have had to be approximated by finite truncations of their infinite expansions, the error induced by such a truncation must be assessed before the results can be interpreted as upper and lower bounds. For example, consider the mean waiting times W . Recall from (V -C(1)) that Z(vO is the objective function obtained by truncating the expression for W. at n= n +1. . Also recall 58 Figure 8 - Convergence of Bounds for Probabilities p(0), p(4) 59 Figure 9 - Convergence of Waiting Time Results: p = 0.7, o,= 0.80, (*,= 0.14) 60 p= 0.5 a =0 .800 P ,=0.142 ~l 1 1 1 1 1 5 10 15 20 ^ 25 30 sim. Figure 10 - Convergence of Waiting Time Results: p = 0.5, o,= 0.80, (*,= 0.14) 61 Figure 11 - Convergence of Waiting Time Results: p = 0.5, 01= 0.95, 0.44) 62 that the computed minimum of Z(Wi) will always be a lower bound for W± . On the other hand, the computed maximum of Z(W±) may or may not be an upper bound for W±f depending on the degree of truncation error. There are thus two opposing forces which affect the convergence of the computed maximum of Z(W±). First, the maximum of Z(Wi) will tend to increase as nc increases because more and more nonzero terms are being added to the expansion for ZfVT) . Secondly, as the LP is solved with more and more constraints there is the tendency for the maximum to decrease toward the true value of W± . Only when the first effect, that of truncation error, has been made negligible is it safe to interpret maxlzfw.^)} as an upper bound for Vv\ . Figures 9-11 show a striking difference in the convergence behaviour of the "bounds" for W, and W2. The maximum of Z(W2) initially seems to be affected by truncation error. Soon a peak is reached, and from there the intrinsic convergence property of the LP method dominates, causing max{Z(W2)}, now an upper bound for W2, to decrease steadily. In contrast, the computed maximum of Z(W,) consistently increases as n • increases (as shown by the precise numerical results), indicating that truncation error is always the dominant factor. Consequently there , is no easily distinguished point beyond which we can say that max{Z(W,)} is an upper bound for W,. Nevertheless, a type of convergence does take place, as the difference betweeen the maximum and minimum values of Z(W,) steadily decreases. For practical purposes this is an entirely acceptable behaviour. One must only be careful to choose n large enough so that the 63 two curves connecting the computed limits of Z(W,) have become nearly horizontal. Figures 9-11 suggest that the tightness of the numerical bounds obtained for a given value of n decreases as either or both p and increase. This was observed also in Figures 2-7. No concrete convergence relationship has been found, however. In consequence, the user of the LP method of analysis will have to assess the convergence characteristics of the system for at least a few different sets of parameters. After this initial experience, the user will be able to select a value of nc which will ensure, for a given set of system parameters, that useful upper and lower bounds will be computed. G. COMPUTATIONAL COSTS This section examines several practical aspects of the LP solution program, chiefly those which affect computational cost. In addition, we compare the relative costs of simulation and solution of the LP. 1. LP Solution Cost A very important practical characteristic of any LP program using the simplex algorithm is the average number of iterations, J, executed before an optimum basic solution is achieved. For a general (m x n) LP, experience shows that on average at least J = m iterations are required to optimize each objective function [LUEN 73, p.55]. In our more specialized algorithm however, in which all but two of the basis variables are known 64 a priori, many fewer iterations are required. On average (over parameters p, p,) we have found that J ranges in a roughly linear fashion from J = 4 at nc= 6 (m=45), to J = 6.5 for nc= 30 (m=l929). Recall that m is the number of constraints in the LP, which for this particular queueing model is given by (5.23). It is apparent from these statistics that an enormous computational savings has been realized through the use of the modified LP technique. A complete LP solution must maximize and minimize a number of different objective functions for a given set of parameters {p, p,, nc ) : W1f W2, and p(n), 0 < n < (nc+l). The total number of iterations required for this task can be reduced considerably by employing a "parallel search" strategy in which at each iteration the new basic solution is tested for optimality against all remaining objective functions. An economy is achieved whenever a single basic solution optimizes more than one objective function. We have found that with the complete set of (n^+4) objective functions listed above, the average number of iterations required per objective function is approximately 1.5 . This value is roughly independent of the parameters (/>, p1r ), although for very close to 0.0 or 1.0 some decrease is noticeable. To summarize then, J has been reduced from J = m for a general LP algorithm, to J = 6 for the specialized LP, and further to J = 1.5 when a parallel search is employed. Computational cost primarily depends on the total number of iterations necessary to optimize all the objective functions. 65 The only additional expense is the initial set-up of the basis matrix and its LU-factorization. This is not usually a significant factor unless only a small number of objective functions are being optimized. With the complete set of (nc+4) objective functions we have found that the total execution cost (CPU time plus memory charges) of an LP solution increases roughly as m2 . Since m is itself proportional to for nc>> s (see 5.23), then cost is proportional to n," . The actual solution program was written in FORTRAN and compiled with the optimizing compiler FORTRAN-H. Single precision arithmetic was generally used, except when e: was near 0.0 or 1.0 , or when the computed upper and lower bounds were to .be very close together. Use of double precision increases execution time by about 50%, and of course incurs extra memory costs as well. As a particular example of program performance, consider the run with p = 0.5, B : = 0.44, n£= 20 (m=794), in single precision. It required 31 iterations to maximize and minimize the 24 objective functions. Total execution time was 19 seconds on an AMDAHL-470 V/8 computer. 2. A Comparison with Simulation Costs Included in Figures 9-11 are simulation results for each mean waiting time. Rigid cost comparisons of simulation and numerical solution are of limited generality because of the large number of options available to users of either technique. Additionally, there may be a large difference in effort initially required to develop the respective computer programs. 66 However, for completeness we present the following comparison. The simulation data for each of Figures 9-11 were gathered from runs of 20 blocks of 4000 customers each. No other statistics besides waiting times were collected. To place the LP solution on an equal footing, we consider only the two objective functions W,, W2 . A single precision solution with computational cost equal to one of the above simulations would occur for n somewhere between 25 and 30 . While numerical c solutions for these two values of nc are separated by a factor of two in cost, the convergence of the bounds is sufficiently slow, so that the following question can be answered: For simulations and solutions of the LP (4.17) of equal execution cost, which method provides the better results for the mean waiting times? The answer is wholly unambiguous only for Figure 10 (p = 0.5, o,= 0.8), in which the LP solution produces better bounds. For the other two cases the answer depends on whether W, or W2 is considered. We conclude that the LP solution seems to provide better results than a simulation of equal cost over the (approximate) range of parameters: p & 0.5, c, £ 0.9, S 0.2 . 67 VI. CONCLUDING REMARKS A. LIMITATIONS OF THE MODEL Simplifying assumptions are necessary during the course of any modeling process. They allow the operation of the physical system under study to be made amenable to mathematical or computational methods. Upon completion of the modeling process the systems analyst must, however, re-assess the assumptions of the model in the context of the original physical problem. The applicability of the present analysis of a communications system is constrained by some of the assumptions traditionally made for such systems. In addition, a numerical solution of the Markov chain equations further reduces the potential generality of the analysis. Perhaps the most obvious limitation of the model is that it provides only a steady state solution. In reality, the transient behaviour and other time-dependent solutions for periodic input traffic may be as important as the equilibrium condition. The choice of a continuous-time Markov model is relevant to physical systems which operate in continuous-time, or to those with relatively short packet lengths. Even in steady state, there is the question of how well the assumed exponential distribution for interarrival times (Poisson arrivals). represents the actual message arrival process. However, this assumption is usually an accurate representation of reality when a communications system serves many equal and independent message sources. The assumption of exponentially 68 distributed service times may in practice not be so well justified. If a more general representation of service times is required, one of the family of Erlang distributions may be used. A continuous-time Markov chain model is still possible although it would have greater complexity and require a more costly numerical solution. A second set of constraints arises from the demands of the LP solution method itself. For example, the number of channels (s) may be so large that the increased number of states renders an economical LP solution impossible. In general, this comment may be made of the entire Markov chain modeling process. If a particular system is of such complexity that a Markov chain model would require a large number of state variables, then it might be well to consider in advance the use of alternate modeling techniques, particularly simulation. One assumption which is always violated by a physical system is that which allows infinite queue lengths. A real system will have a finite number, say L, of storage positions available for waiting, customers. When all L positions are filled, arriving messages are blocked from entering the queue. To make the blocking probability negligible, L is chosen large enough so that the system functions almost as if it had infinite waiting room. Unless L' is quite small, an exact numerical solution of the finite Markov chain would not be economical. Otherwise, the LP solution requires a truncation of the set of balance equations. Because these balance equations are ordered with increasing n (total number in the system), it is apparent 69 with a truncation at n = nc, that L should be interpreted as the maximum total number in queue, regardless of customer type. This favors a system of pooled storage, in which any customer forced to queue may be assigned any one of the L storage positions. The method of truncation, and the structure of the Markov chain state space itself, is not as well suited to a system with separate queues for each customer type. B. SUMMARY AND CONCLUSIONS In this thesis we have set up a Markov chain queueing model of a multichannel, combined voice and data communications system, with non-preemptive priority for the data messages. No analytical results are available for queueing systems of this type. The most closely related 2-class M|M|s NPP models for which some analytical results exist either require all service time distributions to be identical, or require s = 1. A linear programming technique was therefore used to obtain numerical bounds on quantities of interest, most notably, the average waiting times of each customer class. Kotiah's LP solution method could not be applied in its original form because of the complexity and structure of the Markov chain model. First of all, the additional normalization constraint (4.15) was added: nc + 1 E p(n) < 1 . (4.15) n=0 This serves to generalize the LP method to Markov chains for which the first constraint (4.13), 70 s-1 I (1 - n/s) p(n) = 1 - p , (4.13) n=0 is alone insufficient to ensure acceptable solutions. A second improvement is the computational savings which are realized when the Markov chain's property of irreducibility is exploited. Both of these changes are widely applicable to LP solutions of Markov chain models other than the one studied here. This thesis contains several results which are not directly related to our use of Kotiah's LP technique. A proof is given for equation (4.13) which shows it to be a valid steady state condition for a wide class of queueing systems. Although it was indicated in [KOTI 77] that (4.13) could be derived through direct summation of the balance equations, the proof given in Appendix D is a more general and intuitively pleasing result. In addition, the present 2-class M|M|s NPP queue serves as a counter-example to the proposition that (4.13) is equivalent to (4.11). A second important analytical result of this thesis is the clarification of mean waiting time behaviour in a heterogeneous 2-class NPP queueing system. We refer especially to the variation of waiting times with the parameter (at constant p), and the relations between the limiting values Wi(o1->0), W±(a1-»1), (i=1,2). The particular change in the waiting time behaviour predicted for p = p* further illustrates the behaviour of queueing systems of this type. The primary concern in using the LP method to determine bounds on mean waiting times is ensuring that the computed results are affected only negligibly by truncation error. 71 Herein lies one of the limitations of a numerical analysis: if an objective function has an infinite expansion it cannot be exactly computed, the degree of error in a truncated approximation usually unknown except through empirical studies. This is the reason why higher order moments of the waiting time distributions were not computed. Although Little's Law can be generalized for these quantities [KLEI 75, p.240], truncation error would be much more severe. For quantities like p(n), however, which have finite representations, true upper and lower bounds can always be computed. Increasing the number of balance equations (m) in the LP reduces the truncation error for the mean waiting time objective functions. For objective functions which can be expressed exactly in terms of a finite number of state probabilities, the upper and lower bounds converge towards each other as the size of the LP increases: each component variable becomes more and more constrained as additional equations are included. A second class of objective functions are those, such as the mean waiting times, for which an exact representation requires an infinite number of the state probabilities. Adding more balance equations to the LP means the expressions for these objective functions span more nonzero terms. Nevertheless, the results show that at least the difference between the upper and lower bounds decreases as the size of the LP increases. These benefits of increasing the size of the LP are gained at the expense of additional computational effort: computational costs were found to increase roughly as m2. 72 Compared with simulation,' the LP solutions were found to provide better results at less cost over the approximate range of parameters p < 0.5, a, <, 0.9 . The use of simulation could be recommended for the analysis of systems operating outside of this region. Whichever method is employed, the use of extrapolation and other heuristic techniques may render an exhaustive formal analysis unnecessary. Such heuristics may only suggest themselves, however, after great familiarity has been acquired with the particular Markov chain model. C. SUGGESTIONS FOR FURTHER STUDY There are two main areas which deserve further study. These are: (1) the study' of solution techniques for infinite Markov chain models, (2) the comparison of multiclass, multiserver systems with various service disciplines. As discussed in chapter II, the two main numerical techniques for truncating and solving infinite Markov chains are those of Seneta-Wolf, and Kotiah. Although we have employed the LP method of Kotiah, it would be interesting to compare our results with those which could be obtained using the alternate approach of Seneta and Wolf. It would be particularly useful to compare the convergence behaviour of the two methods. Further effort could also be directed towards a search for convenient and useful system approximations. One such approximation might be to somehow combine states whose individual probabilities are not explicitly required, but whose aggregrate probability is 73 important. For example, combine sets of states with fixed q,, q2, but with different server configurations, into the single state (*; q, q2). Aside from further work in the methodology of Markov chain solutions, the techniques themselves might profitably be applied to the analysis and comparison of other multiclass, multiserver queueing systems. For example, a comparison could be made among the various service disciplines: cutoff priority, preemptive priority, non-preemptive priority, and general alternating priority. A less comprehensive comparison might include only an analysis of mean waiting times. Limiting arguments similar to those employed in this thesis could be used to estimate the dependence of W± on the system's mixture of customer classes (fi±). Alternatively, a single system might be studied in more detail, perhaps by making the model correspond more realistically to a true communications system, or by determining the minimum number of servers to ensure a given level of system performance. 74 APPENDIX A - A COMPARISON OF NON-PREEMPTIVE PRIORITY AND  ALTERNATING PRIORITY SERVICE DISCIPLINES In this appendix we discuss and compare two specific service disciplines, non-preemptive priority for type 1 customers (denoted NPP), and the alternating priority case C(°°,1). Recall that C(°°,1) means that queue 1 is served exhaustively, then queue 2 has at most one customer admitted for service. The two service disciplines, general alternating priority C(k,,k2) and NPP are in general quite different; it is only the particular case C(°°,1) that is nearly indistinguishable from the NPP system. In fact, when the rules for queue switching are properly chosen, the C(°°,1) system can be made to function exactly the same as the NPP discipline. In the following we assume for sake of convenience, that n>s , because all work-conserving service disciplines should function identically whenever n<s . queue 1 , [ x arrivals central s servers controller Figure 12 - Schematic of a 2-Class Multiserver Queueing System 75 To explain the necessary queue switching rules we make use of the diagram in Figure 12 . Each service discipline may be interpreted as a set of rules by which the central controller of Figure 12 selects waiting customers for service. The controller is assumed to be able to sense the idle/busy state of each server, as well as the empty/not-empty state of each queue. With this information the controller can direct waiting customers to free servers, switching queues if necessary. Queue transitions, represented by changes in the state variable i , and all other operations of the controller are assumed to be instantaneous. In the NPP discipline the central controller normally attends queue 1 (i=l). Whenever a server becomes idle the controller will, if queue 1 is not empty, select a type 1 customer for service. If queue 1 is empty and queue 2 is not, then a type 2 customer is immediately served (i=2), and the controller returns to monitor queue 1 (i=l). In one version of the C(°°,1) system which .is not identical to the NPP discipline, the controller spends slightly less time on average observing queue 1. This occurs because the controller switches from queue 1 (i=1) to queue 2 (i=2) as soon as it senses that queue 1 is empty and that queue 2 is not. In general at this instant there need not be any free servers. If a customer is in queue 2 the controller reserves it to be served as soon as a server becomes idle, even if a type 1 customer arrives in queue 1 in the meantime. When one type 2 customer is finally sent to the free server, the controller returns to 76 monitor queue 1 (i=l). In this same situation the NPP controller would have first been observing queue 1 until a server became available and so would have selected for service the type 1 customer just arrived in queue 1. To make the C(°°,1) discipline equivalent to non-preemptive priority for type 1 customers, the queue switching of the controller must be made server dependent rather than queue dependent. That is, the controller should switch from queue 1 to queue 2 independently of exactly when queue 1 becomes empty, and only when a server becomes available for the next potential type 2 customer. The switching rules then turn out to be exactly equivalent to those described above for the NPP system. Consequently, this version of the Ci™,}) discipline would select for service the same customers in the same order as would the NPP service discipline. In fact, it is not even necessary in a Markov chain model of these simple systems that the operation of the controller be explicitly represented by a state variable (i), as long as the controller is understood to operate as described above. The reason is that once the state variables n, s,, q, have been specified, the position of the controller is completely determined, and therefore i is not an independent state variable. 77 APPENDIX B - THE STATE SPACE TRANSITION DIAGRAM As a visual aid to the development of the steady state balance equations it is useful to construct the state space transition diagram. The utility of such a diagram, depends crucially on the choice of state variables used to represent the state space. After some experimentation we have chosen to split the state space of the NPP model into two pieces: (1) n < s , 2-dimensional, with variables (s1f s2) (2) n 2: s , 3-dimensional, with variables (s, , q1f q2). For typical states in each segment, the conditional transition rates (Q^ ) to neighbouring states are indicated. There can be at most four such transitions from each state, representing arrivals and departures of each of the two customer types. Forbidden transitions are those prohibited by the service discipline, or those naturally impossible (for example, no type 1 departures when s^O). Figures 13 and 14 below are one set of state space transition diagrams for the multiserver 2-class NPP queueing system. The two diagrams have only the (s+1) states for n=s (<3i=Q2 = 0) in common, and transitions between the two diagrams along this boundary are noted appropriately. 78 Figure 14 - State Space for n ^ s 79 APPENDIX C ~ SETTING UP THE BALANCE EQUATION COEFFICIENT MATRIX "A": AN EXAMPLE This appendix presents further details concerning the structure of the balance equations. We also show with an example how the infinite system of equations is ordered and truncated to form matrix A (see (4.9) in IV - C). The three numbered sections of this appendix are in the order they are referenced in the main text. Part (1) Recall that we found the state space of the 2-class multiserver non-preemptive priority system to be three-dimensional (see III - D). However, to set up and solve a system of balance equations the states must be ordered sequentially in one dimension. States can be indexed in order of increasing population n , but within a group of states with the same n , an ordering should be chosen which: (1) facilitates automatic generation of the balance equations, (2) minimizes the effort required to factorize A into triangular components: A = LU (e.g. by minimizing the fill-in of zero elements). One ordering found satisfactory was to index the states first with increasing n , then with increasing q, , and then with increasing s2 . The result for the case of s = 4 servers, using the balance equations (4.3 - 4.5) of IV - B is seen in Figure 15. This diagram displays the upper left-hand corner of the coefficient matrix QT for the infinite set of balance equations. The shaded squares represent nonzero coefficients. Figure 15 - Matrix "A" For T\.=6, S=4 Servers 81 Indicated along the top and side are the groups of states for which n = 0,1,2,... . The balance equations, for states with n > s = 4 , display a very regular matrix structure. This is easily seen when the coefficients are blocked into square submatrices of dimension s+1 = 5 . When the actual balance equation coefficients are written it is seen that the entire matrix of coefficients for n > s could be represented using only five different (s+1 x s+1) submatrices. The number of submatrices (five) is independent of s . With such a predictable structure, the truncated matrix A can be easily generated up to any preselected cutoff n = nc . Part (2) The figure shows for the particular case n = 6 how the infinite matrix of coefficents would be truncated to N-j. equations and Nj variables. The heavy black line encloses all of the balance equations for n =0,1,2,.. ,nc=6 . Before deriving and we first must determine the number of states which exist in the state space (s, s2; qi q2) , for any fixed population n = n0 ^ 0 . For each n < s there will be (n+1) states since q,=0t q2=0, and s,+s2 = n . For n > s , then qi+q2 = n-s , and therefore there are (n-s+1) configurations for (q, q2) . Since s,+s2 = s for n £ s , the number of configurations for (s, s2) is (s+1) . There are then (s+1)(n0-s+1) states for any given n = n0 ^ s . The number of variables Nj is equal to the number of states for which 0 ^ n < n +1 . Therefore, with n £ s , c c 82 nc +1 N = 1 + 2 + ... + s + E (s+1)(n-s+1) J n=s or, N = s(s+1)/2 + (s+1)(n-S+2)(n_-s+3)/2 (C.1) J C l_ Then Nj =60 for the example nc = 6 . The number of equations Nx might at first appear to be only those for which 0 ^ n < n . But on examination of the c balance equations in the nc+1 =7 group, three more equations are found which do not involve any more than the N^ variables already included (see Figure 15). These three additional equations can be included in the total number N^ of balance equation constraints. For this example, then N =40+3 =43 . Maximizing the number of constraints for a fixed number of variables in this way, we are able to obtain the best possible bounds [KOTI 77]. In general, for a truncation at nc , there will be (n,-s+l) additional equations, and the truncated matrix A will have dimensions (Nj_ x Nj) ' Nj given by (C.1), and n N = 1 + 2 + ... + s + E° (s+1)(n-s+1) + (nc-s+l) n=s or, N = s(s+l)/2 + (s+1 )(nc-s+1 )(nc-s+2)/2 .+ (n,-s+1). (C.2) Part (3) Figure 15 also illustrates the irreducibility [BHAR 60, p.14] of the continuous-time Markov chain model. In any given row, the state on the matrix diagonal can be reached from any of the states with nonzero coefficients (shaded squares) on the left or right. This is the basic structure of balance equations for many queueing systems; any state can be reached from any other state. Hence, the Markov chain is rreduc ible. 84 APPENDIX D - A DERIVATION FOR THE NORMALIZATION CONDITION GIVEN IN [KOTI 77] Consider a queueing system with s identical servers and v customer types. The average customer arrival rates are Xi = a±\ (i=1,2,...,v) , where Z±<*± = 1 > anc5 is the total average arrival rate. The average service times are y± , (i=1,2...,v). The following reasoning is guided by that of Kleinrock [KLEI 75, p.19]. Assume that a stationary probability distribution for the number of customers in the system p(n) , exists, with oo I p(n) = 1 . (D.1) n=0 Assume that all customers are eventually served, 0 < p < 1 . In addition, the system must be work-conserving [KLEI 76, p.113], in particular, no customer can have a service time shorter or longer than that which is obtained when its service time distribution is sampled. Beyond this proviso, there are no other restrictions to be placed on either the service time or interarrival time distributions. The only additional restriction on service discipline is that individual customers have an equal chance of being served by any of the s identical servers. First we show that the utilization factor p , may be interpreted as: p = Pr(any specific server is busy at any random time). Or, that p = 1 - z , where z = Pr(any specific server is free at any random time). During an arbitrarily long time interval r , the expected number of arrivals to one server is (X./S)T . Also during T , the busy time of the server is 85 T(1-Z) , and so the average number of customers served during T will be T(1-Z) / E(y) . Here, E(y) is the expected service time of the customers being served: E(y) = E[ E(y | customer type in service)] = E^)^ . In interval T , the number of arrivals roughly matches the number served: U/s) T = T(1-Z) / E oiyi . As T -> °° , there is equality, and: U/s) E. o y = 1 - z i i i E± xiyi/s = 1 ~ z • Defining the utilization factor as *  = z±  k±y±/ s (D- 2) then z = 1 - p . (D.3Now we develop an' alternate expression for z using (D.1) and conditional probability. z = Pr(any specific server is free at any random time) = E[ Pr(any specific server is free ... | number in system)] = 1 p(0) + (1- 1/s) p(D + (1- 2/s) p(2) + ... ... + (1- (S-1)/S) p(s-1) +0 + 0 + 0 + ... s-1 z = E (1 - n/s) p(n) (D.4) n=0 Therefore, combining (D.3) and (D.4), the desired result is obtained: s-1 1 - ft = E (1 - n/s) p(n) , 0<p<1 . (D.5) n=0 This equation appears in the text as (4.13) in IV - D. 86 BIBLIOGRAPHY [ALLE 77] B. Allen, R.S. Anderssen, E. Seneta, "Computation of Stationary Measures for Infinite Markov Chains" TIMS Studies in the Management Sciences M.F. Neuts, ed., 7 1977 13-23. [AVI 65] B. Avi-Itzhak, W.L. Maxwell, L.W. Miller "Queueing with Alternating Priorities" Opns. Res., J_3 1965 306-318. [BART 81] P. Bartfai, J. Tomko, Point Processes and  Queuing Problems, Colloquia Mathematica Societatis Janos Bolyai, 24 (1978) North-Holland, New York, 1981. [BHAR 60] A.T. Bharucha-Reid, Elements of the Theory of Markov Processes and Their Applications, 1st ed. McGraw Hill, New York, 1960. [BOBI 76] P.A. Bobillier, B.C. Kahan, A.R. Probst Simulation with GPSS and GPSSV Prentice-Hall, Englewood Cliffs, N.J., 1976. [BUNC 76] J.R. Bunch, D.J. Rose, Sparse Matrix Computations Academic Press, New York, 1976. [COBH 54] A. Cobham, "Priority Assignment in Waiting Line Problems", Opns. Res., 2 1954 70-76. [COOP 69] R.B. Cooper, G. Murray "Queues Served in Cyclic Order" Bell Syst. Tech. J.,- 48(3) 1969 675-689. [COOP 70] R.B. Cooper, "Queues Served in Cyclic Order: Waiting Times" Bell Syst. Tech. J., 49(3) 1970 399-413. [COOP 72] R.B. Cooper, Introduction to Queueing Theory, Macmillan, New York, 1972. [DAVI 66] R.H. Davis, "Waiting-Time Distribution of a Multi-server Priority Queueing System" Opns. Res., J_4 1966 1 33-136. [DUFF 77] I.S. Duff, "A Survey of Sparse Matrix Research" Proc. IEEE, 65(4) 1977 500-535. [EISE 71] M. Eisenberg, "Two Queues with Changeover Times" Opns. Res., 19(2) 1971 386-401. 87 [EISE 72] [EVAN 74] [FISH 78] [GOLD 81] [GROS 74] [GUST 72] [HALF 72] [KALA 71] [KLEI 75] [KLEI 76] [KOTI 73] [KOTI 77] [KUEH 79] M. Eisenberg, "Queues with Periodic Service and Changeover Time" Opns. Res., 20(2) 1972 440-451. D.J. Evans, ed., Software for Numerical  Mathematics, Academic Press, New York, 1974. G. S. Fishman, Principles of Discrete Event  Simulation, Wiley, New York, 1978. H. M. Goldberg, "Computation of State Probabilities for M|M|s Priority Queues with Customer Classes having Different Service Rates" INFOR, J_9(1) 1981 48-58. D. Gross, CM. Harris Fundamentals of Queueing Theory Wiley, New York, 1974. F.G. Gustavson, "Some Basic Techniques for Solving Sparse Systems of Linear Equations" in [ROSE 71], 41-52. S. Halfin, "A Priority Queueing Model for a Mixture of Two Types of Customers" SIAM J. Appl. Math., 23(3) 1972 369-379. J.E. Kalan, "Aspects of Large-Scale, In-Core Linear Programming" Proc. Annual ACM Conference, Chicago, 1971 304-313. L. Kleinrock, Queueing Systems, Volume I:  Theory, Wiley, New York, 1975. L. Kleinrock, Queueing Systems, Volume 11:  Computer Applications, Wiley, New York, 1976. T.C.T. Kotiah, N.B. Slater, "On Two-Server Poisson Queues with Two Types of Customers" Opns. Res., 2J_ 1973 597-603. T.C.T. Kotiah, "On a Linear Programming Technique for the Steady-State Behaviour of some Queueing Systems", Opns. Res., 25(2) 1977 289-303. P.J. Kuehn, "Multiqueue Systems with Nonexhaustive Cyclic Service" Bell Syst. Tech. J., 58(3) 1979 671-698. [LUEN 73] D.G. Luenberger, Introduction to Linear and  Nonlinear Programming Addison-Wesley, Reading,Mass., 1973. 88 [MURT 81 ] [NAIR 76] [REID 74] [ROSE 71] [ROSS 72] [SCHL 81] [SENE 67] [SENE 68] [SHER 78a] [SHER 78b] [SING 81] [SLAT 73] B.A. Murtagh, Advanced Linear Programming; < Computation and Practice McGraw-Hill, New York, 1981. S.S. Nair, "Alternating Priority Queues with Non-Zero Switch Rule" Computers and Opns. Res., 3 1976 337-346. J.K. Reid, "Direct Methods for Sparse Matrices" in [EVAN 74], 29-47. D. J. Rose, R.A. Willoughby Sparse Matrices and Their Applications Plenum Press, New York, 1971. S.M. Ross, Introduction to Probability Models Academic Press, New York, 1972. W. Schlee-KSssler, "A Single Server Queue with Two Types of Customers and Alternating Service in Pieces of ^ and " in [BART 81], 343-357. E. Seneta, "Finite Approximations to Infinite Non-Negative Matrices" Proc. Camb. Phil. Soc, 63 1967 983-992. E. Seneta, "Finite Approximations to Infinite Non-Negative Matrices II: Refinements and Applicat ions" Proc. Camb. Phil. Soc, 64 1968 465-470. A.H. Sherman, "Algorithms for Sparse Gaussian Elimination with Partial Pivoting" ACM Trans, on Mathematical Software 4(4) 1978 330-338. A.H. Sherman, "Algorithm 533: NSPIV, A Fortran Subroutine for Sparse Gaussian Elimination with Partial Pivoting" ACM Transactions on Mathematical Software 4(4) 1978 391-398. R. Singh, S. Subba Rao, S.C. Gupta "Analysis of a Mobile Radio Communication System with Two Types of Customers and Priority" Int. Conf. on Comm., Denver, 1981 IEEE (ICC-1981), v.3 57.4.1-57.4.5 . N.B. Slater, T.C.T. Kotiah "The Steady-State of a Multiserver Mixed Queue" Adv. in Appl. Prob., 5 1973 614-631. 89 [STEW 78] W.J. Stewart, "A Comparison of Numerical Techniques in Markov Modeling" Communications of the ACM, 2±(2) 1978 144-152. [SYKE 70] J.S. Sykes, "Simplified Analysis of an Alternating-Priority Queue with Setup Times" Opns. Res., _[8 1970 1182-1 1 92. [TAYL 80] I.D.S. Taylor, J.G.C. Templeton "Waiting Time in a Multi-server Cutoff-Priority Queue, and its Application to an Urban Ambulance Service" Opns. Res., 28(5) 1980 1168-1188. [WALL 66] V.L. Wallace, R.S. Rosenberg "Markovian Models and Numerical Analysis of Computer System Behaviour" AFIPS Spring Joint Comp. Conf, 28 1966 141-148 [WILL 80] T.M. Williams, "Nonpreemptive Multi-server Priority Queues" J. Operational Res. Soc. 31(12) 1980 1105-1107. [WOLF 80] D. Wolf, "Approximation of the Invariant Prob ability Measure of an Infinite Stochastic Matrix" Adv. in Appl. Prob., 12 1980 710-726. 90 Glossary of Notation A coefficient matrix of the Nx balance equations A, first (N-j-x Nx) submatrix of A b, b right-hand side constant (vector) of the normalization constra int(s) B coefficient matrix of the ^ normalization contraints c counter in the state space model for cyclic priority c% general confidence interval for simulation results C(k,,k2,..,kv) cyclic priority system with v queues d number of state variables in a Markov chain E['] expectation operator FCFS first-come, first-served service discipline i,j,k integer indices I identity matrix ^1 (m x m) identity matrix J average number of simplex iterations required per objective function ^ maximum number served during phase i of cyclic priority K total number of objective functions for an LP problem L a fixed maximum queue length LP linear program LU a factorization of a matrix as a product of lower and upper triangular matrices m number of rows in the LP M|G | 1 single-server queue with a Poisson arrival process and a general service time distribution M|M|s s-server queue with a Poisson arrival process and exponentially distributed service times 91 n (1) number of columns in the LP, (2) index for number of customers in a queueing system n0 fixed value of n, number in the system nc truncation parameter for the set of balance equations N (1) total number of states in a finite state space, (2) total number of measured events in a simulation number of balance equations (rows) in A N number of state probability variables associated with columns of A N number of normalization constraints (rows) in B N NPP non-preemptive priority p(n) steady state probability for n customers in a queueing system p£(n), pu(n) lower and upper bounds on p(n) for each n P discrete-time Markov chain transition probability matrix P,, P2 submatrices of P , .P (m x m) truncation of P P(k) general discrete probability distribution P^(k), Pu(k) lower and upper bounds on P(k) for each k P(s, s2) same as Pfs, s2; 0 0) P(s, s2; q, q2) probability of state (s, s2; q2) q± state variable for number of type i customers in the queue Q continuous-time Markov chain transition rate matrix R, , R2 ratios of limiting values of average waiting times s total available number of servers sc cutoff level for cutoff priority queues s± state variable for number of type i customers in service t time (superscript) matrix transpose 92 v (1) number of queues in cyclic priority, (2) number of customer classes in a general system w± random variable for waiting time of type i customers Wi mean waiting time for type i customers x± an indexed member of the set {Pts, s2; qi q2)} x finite column vector of variables in the LP, (state probabilities plus the slack variable) x^ infinite column vector of the probabilities xi x^ a basic solution for a linear program y (1) slack variable, (2) random variable for service time y_^ random variable for service time of a type i customer y average service time of a type i customer z probability that any given server is free at any random time Z, Z^. one of the objective functions for the LP Z(Wi) truncation of VJ± used as an objective function 0 null vector a± fraction of arriving customers which are of type i Oi* value of in a 2-class system at which the p, vs. Oi curve has slope=1 fi± fraction of system resources allocated to customers of type i 6(') discrete delta function A convergence parameter in a matrix iteration technique Poisson arrival rate of type i customers X total Poisson arrival rate of customers n± exponential service rate of type i customers TT_ infinite row vector of Markov chain probabilities 93 m-component solution of a truncated Markov chain model permutation matrix utilization parameter for type i customers total utilization (traffic intensity) parameter a critical value of p for NPP multiserver systems arbitrary time interval 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items