On the design of a finite-capacity queue with phase-type service times and hysteretic control

European Journal of Operational Research 62 (1992) 221-240 North-Holland

221

Theory and Methodology

On the design of a finite-capacity queue with phase-type service times and hysteretic control Marcel F. Neuts * Department of Systems and Industrial Engineering, University of Arizona, Tucson A Z 85721, USA Rao Department of Applied Statistics and Operations Research, Bowling Green State University, Bowling Green, OH 43403, USA B.M.

Received July 1990; revised December 1990 Abstract: We consider the design of a single-server queue with finite buffer capacity. Arrivals follow a Poisson distribution and service times have phase-type distributions. It is possible to switch the service rate between the normal and higher rates. A bi-level hysteretic control policy is considered where two trigger points, say u and 1, are used for changes in service rate. When the number in the system exceeds u, the service rate is increased and it returns to the normal level only when the number in the system drops to l (0 < 1 < u). Algorithms to compute the equilibrium probability vector are presented along with a detailed discussion of various issues arising in their implementation. Numerical results are presented along with discussions of the behavior of the system and of the effect of changes in the various model parameters. The paper concludes with a discussion on how the heuristic understanding of the system behavior gained can be exploited in performing an efficient numerical search for optimal parameter values. Details on extensions of the proposed algorithmic procedure to systems with more general arrival processes are also presented. Keywords: Design; control; queues; phase-type distributions; algorithmic probability

1. Introduction In the design of stochastic service systems it is often inefficient to equip the system to meet the peak demand. One means of improving design efficiency is to operate the server at a lower service rate (possibly with traffic intensity greater than one) and to increase the service rate only when the number of customers in the system builds up. Control of this type is particularly suitable if the arrivals are bursty in nature or if the arrival rate to the system is subject to random fluctuations. Operating the server at a * Research by Marcel Neuts was supported in part by Grants Nrs. DDM-8915235 from the National Science Foundation and AFOSR-88-0076 from the Air Force Office of Scientific Research. 0377-2217/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved

222

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

higher rate only when required may result in substantial reduction in operating costs. These savings must be balanced against additional costs associated with changes in service rates. When the cost of switchover between service rates is negligible, it is appropriate to use a simple control with a single trigger point above which a higher service rate is used. When the switchover cost is significant, hysteretic control [5] may be more appropriate. In that policy, two trigger points, say u and 1, are used. When the number in the system exceeds u, the service rate is increased and returns to the normal slower rate only when the number in the system drops to 1 (0 < 1 < u). Variations and special cases of this model with infinite waiting space have been extensively studied [1-6,13-18]. Gebhard [5] coined the term hysteretic control and for exponential interarrival and service times, derived expressions for the equilibrium probabilities. Yadin and Naor [18] discussed the same model with an infinite number of service rates and gave optimization procedures for such systems. Crabill [2] considered the optimal control of a M / M / 1 queue with an arbitrary number of possible service rates. Considering the general costs dependent on the state of the system and the costs associated with each service rate and zero switching costs, he established that the optimal service rate is a nondecreasing function of the state of the system. Heyman [6] considered an M / G / 1 queue where the server can be turned off to effect some cost savings. He included startup, shutdown and operating costs for the server and a control policy of turning the server on when n customers are present and off when the server is empty. Using a dynamic programming formulation, he derived the optimal value of n as a function of the system parameters and the costs. Tijms [15] studied an M / G / 1 queue with two controllable service rates and hysteretic control with changes in service rates occurring at service completion epochs and developed an algorithm to determine optimal values of u and 1. The general approach relies on the truncation of the system state at some large value N. Federgruen and Tijms [4] developed a computationally tractable procedure for computing the steady-state probabilities for the model in [15] for any given control policy. In controlled queueing systems involving switching costs, determination of an optimal operating policy is in general a difficult theoretical problem [15]. Considering other costs (e.g. cost of lost customers) adds new dimensions to the problem, further complicating its analytical treatment. One of the objectives of this paper is to demonstrate that such problems can be solved directly and efficiently using algorithmic analysis. In this paper we consider an M / G / 1 queue with a finite waiting space and service times having a phase-type distribution [10] and a two-level hysteretic control. We develop algorithms to compute the equilibrium probabilities and use them in deriving expressions for many characteristics of interest. These algorithms can be used interactively to perform an efficient search for optimal values of u, l and k. Such interactive algorithms permit easy modification of the cost structure and frequently lead to the optimal solution with little effort. Extensive numerical results are presented to provide insights into the behavior of the system and to understand the effect of key parameters of the model. We conclude with a discussion of how these algorithms can be extended to systems with more general input processes.

2. Mathematical model

The arrival process of customers is characterized by a Poisson stream of A.There is a finite waiting room of size k and customers arriving when the number in the system is k + 1 (i.e., the waiting room is full) are assumed lost. The normal service time has a phase-type probability distribution with m phases and representation [/3, S]. When the number of customers in the system exceeds u, the service time distribution is changed to one with representation [/31, $1]. The service time distribution returns to normal when the number in the system drops to l (0 < I < u). Let p~ and P2 denote the traffic intensities at the two service rates, defined in the usual manner. In many cases, the change in service time distribution can be viewed as an increase in the service rate by a factor of 0 (0 > 1), and we can write fl~ = fl and S~ = 8S. In some applications, the acceleration of the service may be accompanied by changes in the nature of the service (e.g., elimination of a stage of

M.F. Neuts, B.M. Rao / Finitequeue with hystereticcontrol

223

inspection) which may result in a change in the number of phases in the service time distribution. When S and S 1 have different n u m b e r of phases, we require a more elaborate mathematical notation but the following mathematical analysis remains essentially the same. To avoid such uninteresting complications, we assume that S 1 is also of dimension m in the following analysis. A state of the system is the triple (i, j, s), where i denotes the n u m b e r of customers in the system including the customer in service, j the phase of the service time distribution when the server is busy, and s the state of the server, s is 1 or 2 depending on whether the server is operating with the normal or the modified service time distribution. Let i, 1 < i < u, represent the set of states {(i, j, 1)11 < j < m} and i', l + 1 < i < k + 1, the set of states {(i, j, 2): 1 < j < m}. The system dynamics can then be described by a continuous-parameter Markov chain on the state space [0, 1, 2 . . . . . u, (1 + 1)', (l + 2)', . . . . (k + 1)']. The infinitesimal generator Q of this Markov chain can be represented as shown in Figure 1, where Aoo = - A , Aol = )tO, A1 o = S ° and A~+l,k+ 1 = S 1. The other submatrices are as follows:

Ai,i=S-)tI,

l
A~,i = S 1 - )tI,

l + l <_i <_k,

Ai,i_I=S°O,

2
A'i,i_I=S°O1 ,

l+l
Ai,i+ 1 = ) t I ,

2 <_i < u ,

A'i,i+a=)tI,

l+ l
where S O= - S e , S ° = - S l e , e a column vector of ones a n d I an identity matrix of dimension m. In Figure 2 we display the infinitesimal generator for 1 = 2, u = 4, and k = 6. Let x 0 be the probability of an empty system and xi, 1 < i < u, and Yi, l + 1 < i < k + 1, respectively, m-vectors of equilibrium probabilities of the sets of states i and i'. Define z = [x 0, x 1, x 2. . . . . xu, Yt+ ~, Yt+2,--. ,Yk+ 1}, then z is the unique vector satisfying the equation zQ = O, such that the components of z sum to 1. z can be computed by taking advantage of the block structure of Q. The equation z Q = 0 can be rewritten as the following set of equations: - A x o + x l S ° = O, )tX0fl + Xl(S -- A I ) + x z S ° O = O, Axi_ I + xi( S - A I ) + x i + l S ° f l = O ,

2<_i<_l-1,

Axt_ 1 + xt( S - A I ) + Xt+lS°[3 + yt+lS°[31 = O, AXi_l + X i ( S - A I

)+xi+lS°O=O,

Axu_ 1 + x u ( S - A I )

l+ l
= O,

yl+l(S1 - A I ) +yt+2S°01 = 0,

AYi_l + Y i ( S 1 - A I ) +Yi+IS°01 = 0 , AI u + AS u + Yu+ l( S1

-

l+ 2
AI) +yu+2S°fll=

AYi_I + Y i ( S 1 - A I ) +Yi+IS°O1 = 0 ,

0,

u + 2 <_i<_k,

Ay k +yk+lS1 = O.

(1) (2) (3) (4) (5) (6) (7) (8) (9) (lO) (11)

Post-multiplying (2)-(11) by e and simplifying, we obtain: (12)

AX 0 = Xl SO, A x i e = x i + l SO,

l <_i < l - 1 ,

(13)

-Aoo AOl Alo At]

A~+ l,l

At,t-1

Al2

At,I Au-l,u

Al,l+ 1

2

A'u+l,u

Au,u A~+l,t+l A~+I,I+2 A~÷zS+I A'1 + 2 , / + 3

Au-l,u

Figure 1. Infinitesimal generator

Au,u-1

Au-l,u-1

t A,,+t,u+l

Au,u+l

At

"

k,k - 1

t A.+l,u+2

A~,~ A~,k+1 A'k+I,k A~+l,k+]

¢5 "x.

bO

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

0

0 -A

1

SO

2

1

3'

4

4'

5'

6'

7'

AI

S-hi

S°l]

2

3

225

AI

S-AI

3

S°l]

S - AI

AI

S°l]

S - AI

4

s°l],

3'

AI

S l - AI

AI AI

4'

s°l]~

5' 6'

S 1 - AI

AI

sOft,

S, - AI

AI

sl'l],

Si

7' F i g u r e 2. I n f i n i t e s i m a l g e n e r a t o r f o r 1 = 2, u = 4, a n d k = 6

hXie=Xi+lS°+Yl+lS°

A x , e -- - Y l +

,

-

1,

1S1o ,

=Yi+l SO,

(14) (15)

hYie = Y i + l S ° - h X u e , hYie

l
l+ 1
U + 1 < i < k.

(16) (17)

Substituting (12)-(17) into (1)-(11) and simplifying, we obtain: x i = X o f l R i,

l
x t = [xt_ ~ - x ~ e ( f xi=(xi_

(18)

fl)]R,

1-xuef)R,

l+l
Yt+l = x , e f i R l ,

(19) (20) (21)

yi=Yi_lRl+Yt+l,

l+2
(22)

Yu+l = [Y, + x ~ ] R 1 ,

(23)

y i = y . R ~ -~,

(24)

u+2<_i<_k,

Yk+ I = - - A Y k S ~ 1

(25)

w h e r e the matrices R and R 1 are defined by R = A ( A I - heft - S ) - ' ,

(26)

R 1= &(hi-

(27)

h e f t - S1)-1

F r o m Neuts [10,11] we know that the inverses on the right-hand sides of (26) and (27) exist and are strictly positive. Using (19) and (20), x u can be expressed as follows: xu = XoflRUA - 1 =

It_

l RU -t+ 1A - 1

(28)

w h e r e A = [ I + e f l ( R + R 2 + . . . + R ~-1) + e(fl - f l ) R " - t + 1]. Using (18)-(29), an algorithm for computing z is as follows: 1. Set x 0 to an a p p r o p r i a t e l y chosen initial value. A r e a s o n a b l e choice for x 0 would be the probability of e m p t y n e s s of an M / M / 1 / k q u e u e with traffic intensity -~(P, ' + Pz). 2. C o m p u t e x i, 1 < i < l - 1, using (18). 3. C o m p u t e x u using (28). 4. C o m p u t e x t using (19). 5. C o m p u t e x i, l + 1 < i _< u - 1, using (20).

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

226

6. 7. 8. 9.

Compute Compute Compute Scale the

yi, l + 1 < i < u, using (21) and (22). y~, u + 1 < i < k, using (23) and (24). Yk+l using (25). elements of z to sum to one.

Remarks: 1. Concurrent with the increase in service rate it may be possible to control the input rate. Let A1 be the arrival rate when the server is operating at the higher service rate and define 3' = A1/A as the factor by which the arrival rate is reduced. If the change in service time distribution is treated as an increase in the service rate by a factor of O (0 > 1), then the expression for R 1 as given in (27) can be rewritten as

RI=A(O)[A(O)I-A(~)e~-S]

-1,

which indicates that R1 and hence the values of x/ and Yi depend only on the ratio of 3' and 0 and not on their individual values. Therefore, as far as the steady-state probabilities are concerned, decreasing the arrival rate is equivalent to increasing the service rate by the same factor. That dearly does not hold for waiting time distributions. 2. In the computation of the matrix A in (28), it is necessary to accumulate the sum of the powers of R up to R u-t. Since R is a positive matrix, generally speaking, the values of the elements of the powers of R either increase or decrease with the increasing powers of R depending on whether r/, the spectral radius of R, is greater or less than one. From Neuts [11] it can be seen that r / > 1 when tO1 > 1. Thus, when ~7 and u - l are sufficiently large, the elements of the higher powers of R will increase quite rapidly which may result in overflow and numerical instability. Another factor that limits the range of numerical stability for this matrix-analytic algorithm is indicated by (19)-(20) which require a vector subtraction for each vector x i computed. Here again, when u - l is large, accumulated roundoff errors due to repeated subtractions may lead to numerical instability. For these reasons, when /31 >> 1 a n d / o r u - l is large, the algorithms described in this section should be avoided. For such systems, iterative procedures such as those presented in Section 3 are recommended.

3. Iterative solution

An alternative approach to determine z is to solve zQ -- 0 using an iterative method based on sparse matrix algorithms [8]. Although these algorithms are, in general, computationally more demanding than the matrix-analytic method discussed in Section 2, their numerical stability is almost always guaranteed. Since a heuristic search for the optimal values of u, 1 and k is recommended in this paper, iterative methods have an added advantage in the evaluation of z. During the numerical search, the equilibrium vector z may have to be determined several times. The vector z, computed for a given set of u, 1 and k values, provides a good starting solution for a neighboring set of values of u, I and k, thereby resulting in significant computational economy. The basic iterative principle can be described as follows. Any coefficient matrix T (infinitesimal general Q in the present case) can be split as T = T 1 - T 2 where T 1 is a nonsingular matrix. The equation zT = 0 can then be rewritten as zT 1 = zT 2 which suggests the iterative scheme given by

z ( t + 1)T 1 = z ( t ) T 2 where z(t) is the estimate of z obtained at iteration t. Under fairly general conditions convergence to the correct value of z is guaranteed. Many schemes are described in the literature for the choice of T 1 and T 2 and perhaps the simplest of them is the point Gauss-Seidel scheme, in which T 1 = Tut + Td and

227

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

0

2

1

0

- ;t

a/3

1

SO

S-hi

sot3

2 3

3

4

AI S - AI

AI

S°ft

S- AI

AI

S°t3

S - AI

4 sofia

3'

3'

4'

S t - AI

AI

S°fll

S 1- AI

5p

AI

4'

s°lJ1

5'

Figure 3. Coefficient matrix for

AI S 1 --

AI + R I S ° f l l

IP1

T2 = - T , , where T d, Tut and Tlt a r e the diagonal, strictly upper triangular and strictly lower triangular portions of T. An added advantage of the point Gauss-Seidel scheme is that during each iteration the elements of z can be updated in location, eliminating the need to store two estimates of z simultaneously. At the end of each iteration, the elements of the vector may be scaled to sum to one. While any initial estimate satisfying the conservation law of probability converges to the correct z, a judicious choice will most likely reduce the number of iterations necessary. For the infinitesimal generator under study, one can take advantage of the matrix-geometric structure of the initial and tail portions of z. We propose two iterative schemes which exploit the structure of z in reducing the size of the system of equations to be solved. Firstly, using (24), (9) can be rewritten as: AX u + Ay u

+

Yu+x(S1-

RlS°fll) = 0 .

AI +

(29)

Replacing (9) by (29), it is possible to avoid solving the set of equations given by (10) and (1) iteratively. The coefficient matrix of the resulting system of equations, shown in Figure 3, for the example with u -- 4, l = 2 and k = 6, retains the structure of an infinitesimal generator so that the iterations will be numerically stable. The tail portion of z can then be evaluated using (24) and (25). Finally, the elements of z can be normalized to sum to one. In this procedure, referred to as IP1, the number of equations to be solved is reduced from 1 + (k + u - l + 1)m to 1 + (2u - l + 1)m. An alternate reduction in the number of equations to be solved iteratively can be accomplished by considering the matrix-geometric structure of the initial and tail portions of z, as given by (18) and (24). Using (18) and (24), (4) can be rewritten as follows: - Axte ~

+X/+lS°f3

(30)

+ Yl+ lS°gJ1 ~- O.

Replacing (4) by (30), it is possible to avoid solving the equations derived from (1), (2) and (3) iteratively. The coefficients matrix of the resulting system of equations is shown in Figure 4, for the example with u -- 4, 1 = 2 and k -- 6. As can be seen by the submatrix -Ae/3, this matrix does not retain the structure of an infinitesimal generator. To avoid any possible numerical instability resulting from this, one can take advantage of the matrix-geometric structure of the initial portion of the equilibrium probability vector. For an appropriately chosen value of x0, one can compute the values of x i, 1 < i <_ l, which are accurate to within a multiplicative constant, x t can now be used as part of the initial estimate to start the iterations. However, x t should not be updated during the iterations. After convergence, the 2

2 - Aefl

3

s°gJ

4 3' 4' 5'

4

3 AI S-AI soft

S°fll

3'

4'

5'

S I - AI

AI

S°lll

S 1- AI

AI

sO~]l

S 1 - AI + RIS°I[J1

hi S - AI

AI

Figure 4. Coefficient matrix for

IP2

228

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

exact value of x 0 and the other elements of the vector can be obtained by normalization. In this procedure, referred to as IP2, the size of the system of equations to be solved is reduced from 1 + ( k + u - l + 1)rn to 2(u - l + 1)m. It can easily be seen that in both procedures, updating each element (except x 0 in IP1) during each iteration requires 2m multiplications, 2m - 1 additions and 1 division. Hence, the computational effort during each iteration is directly proportional to the number of equations. Even though IP2 requires fewer computations per iteration, it is not always guaranteed that IP2 requires less computational effort than IP2. Results of our numerical experiments are presented in Section 6.

4. T i m e s p e n t at e a c h service rate

The server may be considered to be in one of three macro states, namely, idle, normal (service) rate or higher (service) rate. Clearly, the sojourn times in the idle state are exponentially distributed with mean 1/A. Let the random variables tn and th, respectively, denote the times spent by the server at the normal and higher rates during each visit to the corresponding macro state, while the system is operating under steady-state conditions. Let the random variable thh be the time between the end of an interval t h and the start of the next interval t h. Properties of t , , t h and thh and a knowledge of how they are affected by the values of u, l, and k provide insights into the system behavior and can be of use in its design and control and also in the heuristic search for the optimum values of u, 1 and k. An interval t n can be initiated either by an arrival when the system is empty or by a service completion when the system is in the set of states l + 1. Similarly, t n c a n end either by the server becoming idle or by the server moving into the set of states u + 1 and switching to higher service rate. Thus, tn can be modeled as the time to absorption in a Markov chain with two absorbing states corresponding to the server becoming idle (a l) and switching to higher service rate (az). The infinitesimal generator for this Markov chain is shown below in Figure 5 for l = 2, u = 4, k = 6. An arbitrary interval t~ can start either in the set of states 1 or in the set of states I. Let 7r~' and ~-~' denote the probabilities of an arbitrary interval t. starting in the sets of states 1 and l, respectively. Let t be an arbitrary point in time when the system is operating under steady-state conditions, rr~' and ~-~' can be written as n 77" 1

xoq(O, dt)

Pn(d t) m

~ , yt+ l j q ( l + 1, j , d t ) j=l

p.(dt) where, q(0, dt) and q ( l + 1, j, dt) are the conditional elementary probabilities of commencement of t n in the interval (t, t + dt), given that at time t the system is in state 0 and (l + 1, j), respectively. P n ( d t ) is the corresponding unconditional elementary probability of commencement of t n in the interval (t, t + dt). al al

1 2

3 4

1

0 s O S-AI S°S

2

AI S - ~I S°13

3

4

a2

AI S-AI S°13

AI S - al

)re

0

a2

Figure 5.

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

229

Clearly, q(0, d t ) = A dt and q(l + 1, j, d t ) = S°l(j) dt, where SOl(j) is the j-th element of the vector SOl. pn(dt) can be expressed as p n ( d t ) = [x0A +yl+lSol] dt. 7r~ and 7r~' can now be written as x0A

x0A + Y t + l S ° ' Yl+ l SO

~7

x0A + Yt+ISOi " The interval t n always starts with the initiation of a service. Hence, given the probabilities that t n starts in a given set of states, the relative probabilities of starting in any given phase are given by the elements of l]. Hence, the random variable t n has a phase-type distribution with m u phases and representation [ a n, Tn]. For l = 2, u = 4, and k = 6, Tn is given by, S°l] Tn =

S - hi

hi

S°l]

S - AI

AI

S°t]

S - hi

'

a n is given by [cry'l] ~'~'l] 0 0] where the two subvectors in the positions corresponding to the c o m m e n c e m e n t of the interval from an empty state and from state l + 1 with higher service rate. Using similar arguments, it is possible to express thh as a phase-type random variable with 1 + m u phases and representations (ahh, Thh). C o m p a r e d to G, thh will have an additional phase for the state 0. Othh will have the vector l] corresponding to the set of states l and zeros elsewhere, because thh always starts in the set of states I. The structures of Thh and Othh are shown below for l = 2, u = 4, and k = 6. SO

S-hi

hi S-hi

S°l]

Thh =

AI

s°l]

S-AI

AI

s°l] 0

Olhh =

ol]

0

s-hi

0].

The interval t h always starts in the set of states (u + 1)', initiated by an arrival to the system when the system is in the set of states u and ends with a service completion when the system is in the set of states (l + 1)'. Unlike t n and thh , which are initiated at the instant of the c o m m e n c e m e n t of a service, I h is initiated at the instant of an arrival. Hence, the relative probabilities of starting in any given phase depend on the service phase at the arrival instant and are given by x * , where xu* is the vector x u normalized by the sum of its elements. Hence, t h c a n be described by a phase-type random variable with m ( k - l) phases and representation [a h, Th]. ah, will have the vector x~* corresponding to the set of states u and zeros elsewhere. The structures of T h and a h are shown below for l = 2, u = 4, and k = 6. S 1 - AI

AI

s°l]~

Sl - ,~ I

hi

s°l],

Sl - hi

hi

s°l]l

S 1-hi

AI

s°l],

S1

L =

tl~h

0

0

x u*

0

0].

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

230

Let M, V and C, with subscripts n, hh, and h, denote the mean, variance and coefficient of variation of the random variables t,, thh and t h. The above results, in principle, provide a means of computing the moments and density function of the random variables t~, thh and t h based on the properties of phase-type random variables [10]. Except for cases when m, l, u and k are very small, it is not practical to compute the density functions. In many cases, even computation of the inverses of the matrices T~, Thh and T h can be a very forbidding task in terms of the storage requirements, CPU time and numerical difficulties. Matrix inversion can however be avoided by recognizing that to compute the first two moments of, say, in, it is necessary to evaluate only the vectors T n l e and T~2e. M~ and V~ can be written as M, = -anTnle

= --atnl~,

Vn = 2anT~Ze - M ff = 2atnoJ - M 2

where 0 = TZ le and co = T~2e. q, and oJ can be evaluated by solving the following two systems of equations sequentially: T,,*=e,

(31)

T.oJ = 0 .

(32)

Higher moments can be computed sequentially by solving additional systems of equations. Solution of the system of equations (31) and (32) can be accomplished efficiently by using block LU decomposition of the matrix T. by taking advantage of its block tridiagonal structure. An LU factorization of T. can be expressed as I L2

I L3

Tn = L U =

II U1

C1

U2

I L,

I L,

I

C2

U3

C3 U,_ 1

C,_ 1 U~

The recurrence relations for U~ and Li, i = 1, 2 . . . . . u, in terms of the submatrices are as follows: U1 =B~,

(33)

L i =AiUi-], U~.=BiLiCi_I,

2 < i < u,

(34)

2
(35)

An LU decomposition for a block tridiagonal matrix exists, provided all the block diagonal matrices in the original matrix are nonsingular [7]. From the properties of phase-type distributions, one can easily verify that this condition is satisfied. The equation TnqJ -- e can be expressed as LUqt=Lv=e

where

v=U~.

The vector qJ can now be solved in two stages by solving two separate systems of equations. In the first stage the intermediate vector v is evaluated by solving the system of equations L v = e. This can easily be accomplished using the following equations: V1 =e,

vi=e-ZiVi_l,

2
where v i, 2 < i < u, is the i-th subvector of v. In the second stage the vector • is evaluated by solving the system of equations Uq, = v.

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

231

This can be accomplished by solving the equations

Ou = V21Vu,

(36)

I~i ~---Ui- 1(/.,'i+ 1 "~- }tl~i+l),

U -- 1 > i >__1.

(37)

where Oi, 2 < i _
5. System characteristics Using the results of Sections 3 and 4, expressions for various system characteristics can be routinely written down as follows: 1. PI, the equilibrium probability of a lost customer, is given by Yk+le. 2. Ls and V~, the mean and variance of the number of customers in the system, including the customer in service, are given by:

u

k+l

L,= Ei(xie) + E i(Yie), i=1

[--~1

k+lll2(Yie)] i2(xie) + i=~l+ " -L2s"

Vs= i 3. Similarly, Lq and

i=1+1

Vq, the mean and variance of the number of customers in the queue, are given by: k+l

Lq= ~ ( i - 1 ) ( x i e )+ i=2 [ u

•

Vq= ~ , ( i - 1 ) 2 ( x i e)+ i=2

(i-1)(yie),

i=1+1 k+l

•

]

(i-1)2(yi e) - L 2.

i=l+l

4. Using Little's theorem, W~, the mean time spent in the system by a customer, is given by: Ls W~= h ( 1 - P t ) " 5. F h, the fraction of time the server operates at the higher service rate, is given by

k+l

Fh= ~7~ yi e. i=1+1

6. Let the costs per unit time of operating the server at the normal and higher service rates be given by c~ an ch, respectively. Let the fixed cost.of increasing [decreasing] the service rate be denoted by c i [ca]. Let c~ be the cost imputed to each lost customer and let cq denote the cost per unit time associated with keeping a customer waiting. Let C 1, C 2, C 3, C 4 denote the expected costs per unit

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

232

time of (i) operating the server, (ii) service rate changes, (iii) waiting customers and (iv) lost customers under equilibrium condition. Expressions for these costs can be obtained as follows:

C 3 = [ L q l c o,

C4=[AY~+le]c,.

6. Numerical results

In this section we present numerical results to provide insight into the behavior of the system and to study the effects of key model parameters on the system performance measures. Considering the entire range of possible values for the parameters l, u, k, Pl and 0 results in a very large p a r a m e t e r space. To limit the p a r a m e t e r space to a manageable size, we fixed the value of k at 50. Values of u and l have been varied in steps of five. Pl is always selected to be greater than or equal to 1 and 0 is always chosen large enough to yield P2 < 1. This restriction on Pl and P2 does not result in any loss of insight into the system behavior because the most interesting behavior patterns are observed only when Pl > 1 and P2 < 1. When Pl and P2 are less than one, the server may enter the higher rate only on rare occasions with the result the system spends most of the time at the normal service rate. On the other hand, if Pl and P2 are greater than one, the system spends most of the time at the higher service rate. In either of these cases, study of the numerical results is not likely to lead to interesting observations. To study the behavior of the system under a wide variety of probability distributions for the service time, we considered seven different phase-type probability distributions with coefficients of variations (CV) ranging from 0.447 to 2.0. These seven distributions, referred to as D 1 - D 7 later in the section, include a subset of the phase-type distributions used by Ramaswami and Latouche [12]. The parameters of these seven distributions described below have been chosen to yield a mean service time of one, so that the traffic intensity can be varied by changing the mean interarrival time. D1 : E6-distribution with a rate of 6 in each phase (CV = 0.44721). D2 : E2-distribution with a rate of 2 in each phase (CV = 0.70711). D3 : Exponential distribution with rate 1 (CV = 1). D4 : Mixture of an E4-distribution with a rate of 6.6778512 in each phase and E2-distribution with a rate of 1.03323700 in each phase, with mixing probabilities of 0.7 and 0.3 (CV = 1.0). D5 : H2-distribution with rates 2.82085228 and 0.50806659 and mixing probabilities 0.6 and 0.4 (CV = 1.5). D6 : Mixture of an E4-distribution with a rate of 8.385265 in each phase and E2-distribution with a rate of 0.607932 in each phase, with mixing probabilities of 0.814075 and 0.185925 (CV = 1.5). D7 : H2-distribution with rates 0.22540333 and 1.77459667 and mixing probabilities 0.11270167 and 0.88729833 (CV = 2.0). The effect of the parameters u and l can be better studied using the Mean [ = ½(u + / ) ] and Spread[ = u - I], denoted by a and b, respectively. As can be seen from Tables 1 and 2, for a fixed value of Mean, varying Spread did not have a significant effect on the mean number in the system and correspondingly the mean time spent by a customer in the system. The variance of the number in the system, however, increased with increasing value of b. For a fixed value of Spread, increasing Mean resulted in an increase in L s as well as Vs. Figures 6 and 7 graphically display the effect of Mean and Spread, respectively, on the equilibrium probabilities of number of customers in the system when the service times follow and E6-distribution. The behavior of t h is governed largely by the value of b. For a fixed value of a, increasing b resulted in an increase in M h and a decrease in Ch. This is not surprising because the dimensions of the matrix T h

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

233

Table 1 Effect of a and b on system performance a

b

xo

Pt

Fh

Ls

V~

Mn

Cn

Mhh

Chh

Mh

Ch

11.381 14.238 17.095 19.953 22.810 19.650 25.317 30.983 36.650 42.317 18.086 24.181 30.276 36.371 42.467 21.727 27.993 34.260 40.527

2.765 3.096 3.395 3.670 3.926 2.086 2.309 2.512 2.701 2.876 2.192 2.391 2.576 2.748 2.911 2.541 2.684 2.846 3.000

20.667 25.667 30.667 35.668 40.668 136.669 174.447 212.226 250.004 287.782 221.667 292.778 363.889 435.000 506.111 377.778 482.222 586.667 691.111

4.274 4.766 5.211 5.621 6.004 1.486 1.669 1.835 1.987 2.129 1.051 1.178 1.298 1.409 1.513 0.962 1.052 1.140 1.224

1.500 1.500 1.500 1.500 1.495 11.334 11.333 11.333 11.333 11.320 21.333 21.333 21.333 21.332 21.319 31.333 31.333 31.332 31.319

2.166 2.186 2.185 2.182 2.150 0.787 0.787 0.787 0.786 0.782 0.573 0.573 0.573 0.573 0.571 0.473 0.473 0.473 0.472

3.062 3.072 3.075 3.075 3.075 21.815 22.499 22.655 22.690 22.698 35.532 40.895 42.287 42.607 42.679 52.181 60.041 62.091 62.564

2.841 2.921 2.950 2.959 2.962 1.044 1.083 1.100 1.107 1.109 0.817 0.791 0.799 0.805 0.808 0.746 0.671 0.665 0.666

3.069 3.074 3.075 3.075 3.075 22.590 22.676 22.695 22.699 22.700 42.119 42.571 42.671 42.693 42.698 62.100 62.566 62.670 62.693

2.932 2.954 2.961 2.962 2.963 1.087 1.104 1.108 1.109 1.110 0.776 0.799 0.806 0.809 0.809 0.648 0.662 0.666 0.668

3.073 3.069 3.054 3.006 2.843 22.692 22.673 22.611 22.403 21.709 42.689 42.666 42.589 42.331 41.472 62.664 62.582 62.309 61.400

3.229 3.207 3.151 3.017 2.723 1.185 1.180 1.166 1.134 1.060 0.865 0.862 0.855 0.838 0.800 0.712 0.708 0.698 0.675

Pa = 1.0, 0 = 1.5, P2 = 0.67, Service time distribution = D2: 20.0 25.0 30.0 35.0 40.0 17.5 22.5 27.5 32.5 37.5 15.0 20.0 25.0 30.0 35.0 17.5 22.5 27.5 32.5 p = 1.25, 0 20.0 25.0 30.0 35.0 40.0 17.5 22.5 27.5 32.5 37.5 15.0 20.0 25.0 30.0 35.0 17.5 22.5 27.5 32.5

0 0.034 0.000 0 0.028 0.000 0 0.023 0.000 0 0.020 0.000 0 0.018 0.000 5 0.038 0.000 5 0.031 0.000 5 0.025 0.000 5 0.022 0.000 5 0.019 0.000 10 0.044 0.000 10 0.034 0.000 10 0.028 0.000 10 0.023 0.000 10 0.020 0.000 15 0.038 0.000 15 0.031 0.000 15 0.025 0.000 15 0.022 0.000 = 1.5, Pz = 0.83, Service 0 0.000 0.000 0 0.000 0.000 0 0.000 0.001 0 0.000 0.002 0 0.000 0.008 5 0.001 0.000 5 0.000 0.000 5 0.000 0.000 5 0.000 0.001 5 0.000 0.004 10 0.002 0.000 10 0.000 0.000 10 0.000 0.000 10 0.000 0.001 10 0.000 0.003 15 0.001 0.000 15 0.000 0.000 15 0.000 0.001 15 0.000 0.002

0.068 0.055 0.047 0.040 0.035 0.077 0,061 0.051 0.043 0.038 0.088 0.068 0.055 0.047 0.040 0.077 0.061 0.051 0.043

10.987 43.067 13.474 63.588 15.965 88.286 18.459 117.148 20.950 150.093 9.776 35.434 12.245 53.804 14.725 76.352 17.211 103.070 19.699 133.922 8.745 31.110 11.162 47.429 13.609 67.908 16.073 92.551 18.545 121.338 10.228 43.976 12.606 62.468 15.024 85.078 17.466 111.816 time distribution = D2: 0.500 20.935 27.314 0.500 25.898 27.148 0.498 30.840 26.177 0.494 35.699 23.960 0.480 40.355 20.126 0.501 18.388 29.535 0.500 23.337 29.907 0.499 28.292 29.344 0.497 33.197 27.703 0.489 37.961 24.528 0.503 15.967 35.132 0.501 20.863 36.385 0.500 25.811 36.306 0.498 30.736 35.068 0.493 35.555 32.297 0.502 18.433 46.186 0.500 23.337 47.018 0.499 28.262 46.249 0.495 33.110 43.780

are independent of a and depend only on the values of b and m. a influences t h indirectly through the vector x u* and will have minimal effect. As might be anticipated, Mhh increased with increasing a or b, and Chh decreased with increasing b and increased with increasing a. W h e n Pl is significantly greater than one, the difference in the values of M n and Mhh was very small. However, Mhh was considerably larger than Mn when Pl is close to one. This can be explained by noting that when Pl is large, the system rarely visits the empty state so that thh is essentially the same as t n. As Pl approaches 1, visits to the empty state become more frequent along with a reduction in the frequency of visits to the higher service rate, leading to significantly larger values for Mhh relative to M,. M n behaved quite predictably by increasing with increasing values of a and b. However, the variability of t, is strongly influenced by Pr When Pt is high, increasing b resulted in a decrease in C,, because the system makes fewer changes in service rate leading to long sojourns at each service rate and hence a

234

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

Table 2 Effect of a and b on system performance

a

b

Xo

Pt

eh

L,

V,

Pl = 1.0, 19= 1.5, P2 = 0.67, Service time distribution = D5: 20.0 0 0.057 0.001 0.227 13.623 93.806 20.0 10 0.057 0.001 0.225 13.884 99.213 20.0 20 0.057 0.001 0.224 1 4 . 3 4 1 110.336 20.0 30 0.057 0.001 0.223 15.052 126.588 22.5 5 0.053 0.001 0.206 14.882 106.412 22.5 15 0.053 0.001 0.206 15.177 114.388 22.5 25 0.053 0.002 0.204 15.696 128.077 22.5 35 0.053 0.002 0.202 16.424 145.973 25.0 0 0.049 0.001 0.191 15.900 115.476 25.0 10 0.049 0.001 0.189 16.122 120.777 25.0 20 0.049 0.002 0.188 16.481 131.432 25.0 30 0.049 0.003 0.185 17.026 146.965 27.5 5 0.045 0.002 0.175 17.142 129.534 27.5 15 0.045 0.002 0.174 17.372 136.987 27.5 25 0.046 0.003 0.171 17.763 149.648 27.5 35 0.046 0.004 0.167 18.288 165.865 30.0 0 0.043 0.002 0.162 18.158 139.890 30.0 10 0.043 0.002 0.161 18.339 144.831 30.0 20 0.043 0.003 0.158 18.601 154.382 30.0 30 0.043 0.004 0.154 18.977 167.916 Pl = 1.25, 19= 1.5, P2 = 0.83, Service time distribution = D5: 20.0 0 0.006 0.003 0.506 21.464 87.090 20.0 10 0.007 0.003 0.506 21.646 94.992 20.0 20 0.009 0.004 0.509 21.774 113.170 20.0 30 0.013 0.005 0.514 22.031 139.508 22.5 5 0.005 0.004 0.499 23.775 89.636 22.5 15 0.006 0.005 0.500 23.810 103.017 22.5 25 0.008 0.006 0.501 23.868 126.207 22.5 35 0.012 0.008 0.504 23.999 156.307 25.0 0 0.003 0.005 0.494 25.830 86.057 25.0 10 0.004 0.006 0.493 25.951 93.994 25.0 20 0.005 0.007 0.492 25.911 112.652 25.0 30 0.007 0.009 0.491 25.482 140.120 27.5 5 0.002 0.007 0.487 28.124 86.683 27.5 15 0.003 0.008 0.485 28.053 100.070 27.5 25 0.004 0.010 0.482 27.885 123.370 27.5 35 0.006 0.014 0.477 27.601 153.848 30.0 0 0.002 0.009 0.481 30.171 81.505 30.0 10 0.002 0.010 0.479 30.228 89.181 30.0 20 0.002 0.012 0.475 30.025 107.471 30.0 30 0.003 0.016 0.467 29.642 134.277

M.

Cn

ghh

Chh

gh

Ch

7.721 11.596 12.027 12.196 12.314 13.324 13.610 13.745 9.626 14.434 14.970 15.179 15.002 16.229 16.577 16.741 11.531 17.273 17.912 18.163

2.409 2.491 2.911 3.327 2.386 2.781 3.173 3.566 2.672 2.664 3.047 3.416 2.576 2.933 3.294 3.645 2.913 2.827 3.148 3.525

21.893 161.406 297.231 433.057 104.087 255.308 406.519 557.729 26.893 198.001 364.595 531.190 125.295 307.288 489.267 671.247 31.893 234.595 431.960 629.324

3.028 1.140 0.906 0.830 1.542 1.031 0.882 0.826 3.349 1.250 0.971 0.866 1.689 1.114 0.933 0.855 3.642 1.353 1.036 0.908

6.419 46.896 85.975 124.256 27.056 66.057 103.296 140.898 6.338 46.219 84.378 120.939 26.541 64.450 101.020 134.225 6.175 44.864 81.182 114.306

2.922 1.079 0.791 0.650 1.398 0.890 0.701 0.594 2.834 1.045 0.766 0.630 1.344 0.856 0.675 0.578 2.695 0.994 0.731 0.610

5.341 24.541 27.704 22.904 20.075 33.598 33.989 26.840 5.770 32.642 42.905 40.429 23.260 45.157 52.442 46.971 6.002 38.761 57.855 62.149

2.267 1.203 1.432 1.892 1.264 1.184 1.451 1.923 2.493 1.110 1.173 1.462 1.286 1.041 1.163 1.469 2.669 1.066 0.998 1.153

5.938 44.930 81.469 115.063 26.510 64.619 101.312 134.984 6.092 46.206 84.462 121.233 27.000 66.050 104.383 141.193 6.170 46.850 85.971 124.343

2.703 0.953 0.699 0.582 1.293 0.820 0.646 0.552 2.840 1.003 0.734 0.603 1.346 0.854 0.671 0.568 2.928 1.036 0.758 0.623

6.074 46.046 84.321 121.552 26.432 64.545 101.679 136.905 5.941 44.938 81.903 117.032 25.624 62.297 97.255 128.854 5.707 42.989 77.649 109.078

3.165 1.145 0.840 0.692 1.477 0.941 0.742 0.632 3.034 1.097 0.807 0.668 1.402 0.895 0.711 0.615 2.849 1.032 0.764 0.645

reduction in variability. However, when Pl is close to one, for a fixed a, increase in b resulted in an increase in Cn. This occurs because, as b is increased, there will be a decrease in the number in the system at the instant the service rate drops to the normal rate. This increases the likelihood of absorption into the empty state a 1, resulting in some relatively small values for t n, causing a marginal increase in variability. In almost all the examples displayed, we notice that the value of x 0 is very close to zero. This is a result of limiting the value of Pl to be greater than or equal to one. On a similar note we can observe that for all the examples chosen the value of Pl is very small. This dearly is the result of the restriction P2 < 1. For systems with P2 > 1, a substantial portion of the arriving stream of customers would be lost to the system. In such systems, the server spends almost all of the time at the higher service rate with the

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

0.1

235

Equilibrium Probability

/~,,'" 008

•

~',,~¢~000"~'~

~

l

/

'

/

Mean = 15

- .....

," / /~ ': ? / f ~: 1t '~'(~ ~>/

"t ~>

1 1

-'~- Mean ` 25 0" Mean = 30

::: 0

5

10

15

20

25

30

35

40

45

50

55

60

Number in system Figure 6. Effect o f M e a n (b = 10, Pt = 1.5, O = 2, k = 50)

result all the x u values will be close to zero. For such systems, the optimal values of the parameters can be determined more easily than the systems considered here. An interesting observation on the system behavior is the relative insensitivity of F h, the fraction of time the system spends at the higher service rate, to the values of u and l. This observation is particularly true when Pl, the probability of a lost customer, is close to zero. Even when the value of Pl is significantly different from zero, F h remained relatively unaffected by the changes in the values of u and l. This observed behavior can be explained by considering a simplified model with exponential service times. Under equilibrium conditions the effective arrival rate must be equal to the effective service rate as given by (1 - P t ) A

= (xo)O + (1 - x o --Fh)lJ, 1 + Fhl.t 2

0.3

Equilibrium

(38)

Probability spread-O ...... spread=lO -spread=20 ...... spread=30

0,25 t 0.2[ 0.15 QI 0.05 0 0

5

10

15

20

25 30 35 40 Number in system

45

50

Figure 7. Effect of Spread ( a = 25, Pt ~ 1.5, 0 = 2, k = 50)

55

60

236

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

Table 3 E f f e c t o f 0 on s y s t e m p e r f o r m a n c e

0

P2

Xo

Pt

Fh

Ls

~

Mh

Ch

0.201 0.113 0.079 0.060 0.049 0.041 0.035 0.031 0.028 0.025

15.268 13.500 12.929 12.654 12.494 12.388 12.314 12.259 12.216 12.182

97.553 71.150 63.929 60.810 59.097 58.020 57.282 56.745 56.337 56.017

152.376 77.192 51.480 38.611 30.889 25.741 22.063 19.306 17.160 15.444

0.836 0.622 0.513 0.447 0.403 0.370 0.345 0.325 0.308 0.294

0.661 0.500 0.400 0.334 0.286 0.250 0.223 0.200 0.182 0.167

26.191 23.231 22.085 21.505 21.157 20.924 20.759 20.634 20.538 20.460

65.980 43.399 36.249 33.439 32.059 31.273 30.776 30.439 30.198 30.019

150.557 77.178 51.509 38.635 30.908 25.757 22.077 19.318 17.171 15.454

0.888 0.677 0.559 0.487 0.438 0.403 0.375 0.353 0.332 0.319

0.826 0.713 0.625 0.556 0.500 0.455 0.417 0.385 0.357 0.333

28.911 25.777 24.394 23.671 23.234 22.942 22.734 22.578 22.456 22.359

64.281 41.714 32.245 28.190 26.173 25.031 24.317 23.839 23.500 23.251

146.804 77.030 51.532 38.662 30.931 25.776 22.094 19.332 17.184 15.466

0.944 0.747 0.621 0.542 0.487 0.447 0.416 0.391 0.370 0.352

Pl = 1.0, u = 30, l = 15, Service t i m e distribution = D I :

1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00

0.910 0.830 0.770 0.710 0.670 0.620 0.590 0.560 0.530 0.500

0.020 0.023 0.024 0.024 0.024 0.025 0.025 0.025 0.025 0.025

0.000 0.000 0.000 0 000 0.000 0.000 0.000 0.000 0.000 0.000

Pl = 1.2, u = 30, 1 = 15, Service t i m e distribution = D I :

1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 Pl

=

1.5, U

1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50

0.920 0.860 0.800 0.750 0.710 0.670 0.630 0.600 0.570 0.550 =

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

30, l = 15, Service t i m e distribution = D I :

0.940 0.880 0.830 0.790 0.750 0.710 0.680 0.650 0.630 0.600

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

where P-1 and /x 2 are the normal and higher service rates. When Pt and x 0 are close to zero, (38) simplifies to yield P2(Pt - 1)

Fh

(Pl --P2) '

(39)

which indicates that F h is independent of the values of u and 1. For systems with values of Pt a n d / o r x 0 significantly different from zero, F h can be expected to display some dependence on the values of u and l as is observed in the numerical results. Even though (38) and (39) cannot be extended to systems with phase-type service times, they can, however, serve to provide an approximate explanation for the observed behavior. Our computations have also confirmed that, for systems with phase-type service times, the values of F h computed by using (39) differed from the values obtained by the methods described in this paper by 0.0001 or less in most cases. The effect of 0, the service rate acceleration factor, is quite predictable (Table 3). Increasing 0 while fixing the other parameters resulted in a decrease in L s as well as Vs. The rate of improvement in L s per unit increase in 0 decreased with increasing 0. Increase in 0 also resulted in a decrease in the mean and variability of T h as well as the fraction of time spent at the higher service rate. Clearly, changes in 0 will not affect the attributes of the random variable t n.

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

237

Table 4 Effect of service time distribution on system performance

Dj u = 30, D1 D2 D3 D4 D5 D6 D7 u = 30, D1 D2 D3 D4 D5 D6 D7 u = 30, D1 D2 D3 D4 D5 D6 D7 u = 30, D1 D2 D3 D4 D5 D6 D7 u ~ 30, D1 D2 D3 D4 D5 D6 D7

xo l = 15, p~ = 0.023 0.028 0.035 0.035 0.050 0.050 0.067 l = 15, Pl = 0.024 0.030 0.038 0.038 0.057 0.057 0.077 l = 15, Pl = 0.000 0.000 0.001 0.001 0.007 0.007 0.016 I = 15, Pl = 0.000 0.000 0.001 0.001 0.006 0.005 0.013 / = 15, Pl = 0.000 0.000 0.000 0.000 0.001 0.001 0.003

Pl

Fh

Ls

1.0, 0 = 1.2, P2 = 0.83: 0.000 0.113 13.500 0.000 0.140 13.935 0.000 0.174 14.592 0.000 0.174 14.523 0.002 0.240 16.046 0.002 0.239 15.956 0.009 0.290 17.090 1.0, 0 = 1.4, P2 = 0.71: 0.000 0.060 12.654 0.000 0.075 12.817 0.000 0.096 13.004 0.000 0.096 13.004 0.000 0.142 13.735 0.000 0.141 13.648 0.002 0.186 14.284 1.25, 0 = 1.5, P2 = 0.83: 0.000 0.333 21.778 0.000 0.334 21.499 0.000 0.335 21.147 0.000 0.335 21.196 0.001 0.342 20.665 0.001 0.341 20.676 0.004 0.348 20.490 1.25, 0 = 1.75, P2 = 0.71: 0.000 0.500 23.229 0.000 0.500 23.337 0.001 0.500 23.493 0.001 0.500 23.533 0.005 0.500 23.810 0.005 0.498 23.804 0.015 0.488 23.875 1.50, 0 = 2.0, P2 = 0.75: 0.000 0.500 23.234 0.000 0.500 23.341 0.000 0.500 23.495 0.000 0.500 23.556 0.001 0.499 23.850 0.001 0.499 23.917 0.008 0.491 24.318

V~

Mn

Cn

Mhh

Chh

Mh

Ch

71.150 79.869 93.937 94.806 128.829 129.962 170.829

35.782 27.993 21.176 21.257 13.324 13.378 9.150

2.664 2.684 2.713 2.710 2.781 2.773 2.845

603.857 482.222 376.000 380.373 255.308 258.578 201.999

1.056 1.052 1.046 1.047 1.031 1.033 1.018

77.192 78.187 79.383 79.995 80.759 81.392 82.484

0.622 0.708 0.806 0.801 0.940 0.932 0.973

60.810 64.909 71.561 72.153 90.746 91.714 122.901

35.782 27.993 21.176 21.257 13.324 13.378 9.150

2.664 2.684 2.713 2.710 2.781 2.773 2.845

603.857 482.222 376.000 380.373 255.308 258.578 201.999

1.056 1.052 1.046 1.047 1.031 1.033 1.018

38.611 39.166 39.993 40.328 42.080 42.464 46.167

0.447 0.521 0.612 0.610 0.771 0.766 0.867

29.439 35.714 46.208 46.213 73.674 74.213 110.639

61.276 60.041 54.072 55.687 33.598 35.455 21.515

0.603 0.671 0.796 0.783 1.184 1.156 1.551

61.798 62.566 63.453 64.022 64.619 65.553 68.621

0.604 0.662 0.728 0.727 0.820 0.823 0.877

30.912 31.349 31.994 32.239 33.630 33.968 36.641

0.447 0.521 0.612 0.610 0.771 0.766 0.872

36.712 47.018 63.769 63.993 103.017 103.555 145.209

61.276 60.041 54.072 55.687 33.598 35.455 21.515

0.603 0.671 0.796 0.783 1.184 1.156 1.551

61.798 62.566 63.453 64.022 64.619 65.553 68.621

0.604 0.662 0.728 0.727 0.820 0.823 0.877

61.800 62.582 63.507 63.951 64.545 65.106 65.500

0.622 0.708 0.806 0.801 0.941 0.932 0.978

26.173 30.652 38.870 38.981 63.872 63.871 97.468

30.930 31.343 31.769 32.082 30.326 31.553 28.051

0.464 0.504 0.558 0.555 0.700 0.681 0.868

30.932 31.362 31.991 32.217 33.580 33.981 36.956

0.464 0.505 0.558 0.556 0.656 0.656 0.731

30.931 31.360 31.981 32.199 33.454 33.791 35.707

0.487 0.565 0.659 0.657 0.819 0.813 0.910

The effect of the service time distribution on the system performance is summarized in Table 4. Increase in variability of the service time distribution resulted in an increase Vs but did not show a significant impact on L s. Given a coefficient of variation, the distributional form (service time distributions D 3 - D 4 and D 5 - D 6 ) did not have significant effect on the overall system performance. This is further confirmed in Figure 8, where the marginal equilibrium probabilities of the system state (i.e., X o e , x i e , 1 < i < u, Yi, I + 1 < i < k + 1) have been displayed. In Figure 7, x-coordinate values 0 and through 35 correspond to states 0, i, 1 < i < 35, and 36 through 71 correspond to states i', 16 < i < 51. In summary, studying the effect of the values of the trigger points in terms of their Mean and Spread values is more revealing than their actual values. As might be expected intuitively, the mean number in the system is controlled largely by the Mean of the trigger points and the variability of number in the system is strongly influenced by Spread. It is also clear that increasing Spread increases the mean time spent by the system at the two service levels during each visit to a state while maintaining the fraction of

238

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control Equilibrium Probability 0031 0.025

•

•-

--I-- D2

f 0,02

D3

f

0.015

D7

0.01

'

/

0.005 / / I , . , r

0

5

10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Slale N u m b e r

Figure 8. Effect of service time distribution (u = 35, l = 15, Pl = 1.5 and P2 = 0.75)

time s p e n t at the higher service rate relatively constant. This indicates that i n c r e a s i n g S p r e a d decreases the n u m b e r of c h a n g e s in service rates. T h e fraction of time s p e n t by the server at the higher service rate is relatively insensitive to the values of a a n d b. T h e n u m b e r of iterations necessary for c o n v e r g e n c e to a n accuracy of 10-1° for the two iterative p r o c e d u r e s IP1 a n d IP2 are s u m m a r i z e d in T a b l e 5. A l t h o u g h n o definitive s t a t e m e n t o n the relative efficiencies of IP1 a n d IP2 can be m a d e b a s e d o n these limited n u m e r i c a l experiments, some g e n e r a l observations c a n be made. W h e n the variability of the service time d i s t r i b u t i o n is low, IP2 p e r f o r m e d b e t t e r t h a n IP1. W i t h i n c r e a s i n g variability of the service time distribution, the relative p e r f o r m a n c e IP1 i m p r o v e d a n d it e v e n t u a l l y o u t p e r f o r m e d IP2. Also, at lower overall traffic intensity IP2 p e r f o r m e d b e t t e r t h a n IP1 while at higher overall traffic intensity IP1 o u t p e r f o r m e d IP2. O n e should, however, c o n s i d e r the c o m p u t a t i o n a l effort p e r i t e r a t i o n t o g e t h e r with the n u m b e r of iterations before choosing a n iterative scheme. V a r i a t i o n s in the expected costs C x, C 2, C 3 a n d C4 as a f u n c t i o n of the m o d e l p a r a m e t e r s d e p e n d largely o n the m a g n i t u d e s of the basic costs c,,, c h, c i, c d, c t a n d cq as d e f i n e d in Section 5. This m a k e s it difficult to draw any g e n e r a l conclusions except the most intuitively obvious. F o r this r e a s o n n o n u m e r i c a l results o n the cost c o m p u t a t i o n s have b e e n p r e s e n t e d .

Table 5 Performance comparison of IP1 and IP2

l

D1 IP1

IP2

Px = 1.5, P2 = 0.75: 10 995 412 15 623 353 20 300 275 25 230 161 Pl = 2 . 5 , P 2 = 0.83: 10 556 717 15 364 600 20 208 444 25 87 241

D2 IP1

IP2

D3 IP1

IP2

D4 IP1

IP2

D5 IP1

IP2

D6 IP1

IP2

D7 IP1

IP2

867 536 259 284

444 370 268 148

524 327 232 258

441 372 275 144

548 343 242 269

434 366 275 152

367 239 293 331

535 419 284 136

383 266 318 355

515 405 278 138

283 301 373 428

510 379 242 110

523 357 208 83

785 633 443 217

361 251 149 81

773 597 393 175

391 275 164 83

764 597 403 193

327 222 128 121

706 506 258 148

358 252 147 126

686 503 309 139

358 256 151 165

569 424 269 120

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

239

7. Concluding remarks The algorithms developed in this p a p e r can be used interactively to perform an efficient search for the optimal combination of trigger points. One can start with a reasonable set of values for l and u and fine tune them until the desired performance measure is optimized. Strategies for the search algorithms can be based on the heuristic understanding of system behavior discussed in Section 6 and any other available information about the system behavior. Such a search procedure permits the analyst to incorporate his or her intuitive understanding of the system behavior into the search process and will likely lead to the optimal solution with little effort. As indicated in Section 1, hysteretic control is most appropriate when the arrival process is bursty in nature, where the arrival rates depend on an external process. Such arrival processes can be very effectively modeled using a Markov-Modulated Poisson Process (MMPP). An MMPP is a doubly stochastic Poisson process in which the arrival rate varies according to the state of an n-state irreducible Markov chain. M M P P is used extensively in many applications particularly in modeling communication systems. An n-state M M P P is parametrized by the n × n matrix q~, which is the infinitesimal generator of the underlying Markov chain and by an n-vector A of rates Ai. At least one of the rates Ai is assumed to be positive. It will be convenient to introduce a diagonal matrix A which has the quantities Ai, 1 < i < n, as diagonal elements. For the system described in Section 2, if the arrival process is replaced by an MMPP, the state description will have to be modified as (i, r, j, s), where r, 1 < r < n, is the state of the Markov chain governing the arrival rates. Arranging the states in lexicographic order, the infinitesimal generator will be of the form shown in Figure 1 with Ate) = q~ - A, A01 A ®/3, A 10 I ® S° and A'k + l,k + I = ~ ® I + I ® S, where I is an identity matrix of order m and ® is the Kronecker product operator. The other submatrices are as follows: =

Ai.i=(~-A)®I+I®S, A~.i=(~-A~)®I+I®S Ai,i_ I =I®S°~,

A~, i I = I ® S I ) ~ I , Ai.i+ ! = A ® I, A'i.i+l=Al®I,

:

l
l+l
2<_i<_u, l+l<_i
The submatrices (except A00, A0~, Ai0) will now be of order nm. Thus, the total number of states is multiplied by a factor of n to n[1 + m ( k + 1)]. The procedure described in this p a p e r can easily be adapted to study the above system. The computations will however be significantly more demanding but can be efficiently organized to take advantage of the special structure. Recently, Lucantoni, Meier-Hellstern and Neuts [9] introduced a new Markovian arrival process which includes the M M P P as a particular case. This arrival process can be used to model many nonrenewal arrival processes and yields mathematical models which are algorithmically very tractable. The discussion presented in this section applies to that arrival process also.

References [1] Cohen, J.W., " O n the optimal switching-level for an M / G / 1 queueing system", Stochastic Processes and Applications 4 (1976) 297-316. [2] Crabill, T.B., "Optimal control of a service facility with variable exponential service times and constant arrival rate", Management Science 18 (1972) 560-566.

240

M.F. Neuts, B.M. Rao / Finite queue with hysteretic control

[3] Crabili, T.B., Gross, D., and Magazine, M.J., "A classified bibliography of research on optimal design and control of queues", Operations Research 25 (1977) 219-232. [4] Federgruen, A., and Tijms, H.C., "Computation of the stationary distribution of the queue size in an M / G / 1 queueing system with variable service rate", Journal of Applied Probability 17 (1980) 515-522. [5] Gebhard, R.F., "A queueing process with bilevel hysteretic service-rate control", Naval Research Logistics Quarterly 14 (1967) 55-67. [6] Heyman, D.P., "Optimal operating policies for M / G / 1 queueing systems", Operations Research 16 (1968) 362-382. [7] lsaacson, E., and Keller, H.B., Analysis of Numerical Methods Wiley, New York, 1966. [8] Kaufman, L., "Matrix methods for queueing problems", SIAM Journal for Scientific and Statistical Computing 4 (1983) 525-552. [9] Lucantoni, D.M., Meier-Hellstern, K.S., and Neuts, M.F., "A single server queue with server vacations and a class of non-renewal arrival processes", Advances in Applied Probability 22 (1990) 676-705. [10] Neuts, M.F., Matrix-Geometric Solutions in Stochastic Models - An Algorithmic Approach, The Johns Hopkins University Press, Baltimore, MD, (1981). [11] Neuts, M.F., "Explicit steady-state solutions to some elementary queueing models", Operations Research 30 (1982) 480-489. [12] Ramaswami, V., and Latouche, G., "An experimental evaluation of the matrix-geometric method for the G I / P H / 1 queue", Stochastic Models 5 (1989) 629-667. [13] Tijms, H.C, "Opmtimal control of the workload in an M / G / 1 queueing system with removable server", Mathematische Operationsforschung und Statistik 7 (1976) 933-944. [14] Tijms, H.C., "On a switch-over policy for controlling the work load in a queueing system with two constant service rates and fixed switch-over costs", Zeitschrift fiir Operations Research 21 (1977) 19-32. [15] Tijms, H.C., "A unified steady-state analysis for controlled Markov drift processes in inventory, queueing and replacement problems", in: P. Bartfai and J. Tomco (eds.), Point Processes and Queueing Problems, North-Holland, Amsterdam, 1981. [16] Tijms, H.C., "An algorithm for average costs denumerable state semi-Markov decision problem with applications to controlled production and queueing systems", in: R. Hartley, L.C. Thomas and D.J. White (eds.), Recent Developments in Markovian Decision Processes, Academic Press, New York, ,1980. [17] Tijms, H.J., and Van der Duyn-Schouten, F.A., "Inventory control with two switch-over levels for a class M / G / 1 queueing systems with variable arrival and service rates", Stochastic Processes and Applications 6 (1978) 213-222. [18] Yadin, M., and Naor, P., "On queueing systems with variable service capacities", Nat,al Research Logistics Quarterly 14 (1967) 43-53.