TI-IE STATIONARY G/G/s QUEUE WITH NON-IDENTICAL SERVERS

We extend a recently developed factorization method to the case of the G/G/s queue with non-identical servers, by presenting three simple proper-lies which lead to a simple numerical calculation method. We compare our results with those determined by classical Markovian (phase) methods in the case of the symmetrical M/G/s queue, and for the mean queueing delay we compare with results given by traffic simulation.


Introduction
In a recent paper [3], we studied the stationary G/G/s queue by means of a new fac- torization method more general than a Wiener-Hopf type of decomposition.In Sec- tion 3, we show how this method may be readily extended to the case of non-identical servers for delayed customers.In Section 4, we then show how the effect of the busy period can be combined with that of the partial occupancies to evaluate the probabili- ty of delay.The calculations are made in both Sections 3 and 4 with the aim of deriv- ing a simple numerical calculation method offered by three simple properties.We apply the results successively to the stationary GI/G/s queue and the M/G/s queue.
In Section 5, we close with numerical comparisons of results obtained by applying classical Markovian methods to the symmetrical M/E2/s and M/H2/s queues, and we compare the calculated values of the mean queueing delay with results given by traf- fic simulation.
We begin, in Section 2, by defining notation and assumptions and by outlining the recently-derived preliminary results for the case of the symmetrical G/G/s queue.164 PIERIE LE GALL Notation, Assumptions and Preliminary Results for the Symmetric G/G/s Queue 2.1 Notation and Assumptions Except for the service time distribution, different for each server, the notations and assumptions will be the same as in Le Gall [3].We consider a queue handled by a multiserver of s non-identical servers.
a) The Arrival Process: We assume a metrically transitive, strictly stationary pro- cess of successive non-negative interarrival times.Let N(t)denote the random num- ber of arrivals in the interval (0,t].We write dN(t)= 1 or 0 depending on whether or not there is an arrival in the elementary interval (t, t + dr).We exclude the possibi- lity of simultaneous arrivals.We can then write E{dN(to).dN(t o + t)} E[dN(to) p(t).dr, where p(t) is the arrival rate at time t + t o if an arbitrary arrival occurred at time t o.We let ezt.p(t).dt-Cl(Z)-O,x(Z),Pt(z) < O, (2)   0 x=l where OO, x(Z corresponds to the xth arrival following the epoch t 0. However, the sta- tionary assumption and the Abelian theorem give that limz__.oZ.Cl(Z -A, where h is the mean arrival rate.In a more general way, we may write for j-1, 2,... E{dN(to)  and for l{e(zj) < 0 and j-1,2,... f e zltl.dtl.../ eZjtj.dtj.fj(7l...7j)--otj(Zl...zj).
b) The Service Times: The successive service times T n are mutually independent and independent of the arrival process.For server j (j-1,...,s) the service times Tn(j) are identically distributed with a distribution function El(t; j); and we let Ti(z;j)-E[eZTn(J)], for Re(z) < 0. We exclude the possibility of batch service and, consequently, F 1(0; j) FI( + 0; j) 0.
We assume Tl(z;j) to be holomorphic at the origin.From Paul Levy's theorem, we deduce that Ti(z;j) exists for Re(z) _ 5, where 5 is a real positive number.
c) The Service Discipline: The servers are supposed to be non-identical with differ- ent service time distributions.But they are indistinguishable for the service discipline which is "first come-first served." d) The Traffic Handled: Loynes [4] demonstrated the existence of the stationary regime.In Section 3, we shall see that the non-identical servers are equivalent (during the busy period in the stationary regime) to different single servers handling the same value y for the traffic intensity per server (at any period), with the necessary and sufficient condition: e) Queueing Delay: Since the term "waiting time" means "sojourn time" in Little's formula, for clarity we prefer to use the term "queueing delay" r for the queue- ing process only and for an arbitrary customer.f) Contour Integrals: In this paper we use (Cauchy) contour integrals along the imaginary axis in the complex plane.If the contour (followed from the bottom to the top) is to the right of the imaginary axis (the contour being closed at infinity to the right), we write f If the contour is to the left of the imaginary axis, we write +0 f Unless it is necessary to specify whether the contour is to the right or to the left -0 of the imaginary axis, we write f. 0 2.2 Preliminary Results for the Symmetrical G/G/s Queue Now, we outline the recent results that were presented in [3] for the case of the sym- metrical G/G/s queue in equilibrium.To avoid very complicated calculations, Le Gall defined the singular points of the function E[e-qr], with Re(q) > 0, relating to the queueing delay v for an arbitrary delayed customer.Secondly, Le Gall established conditions under which this function is holomorphic, these conditions being satisfied by a more general factorization method than the Wiener-Hopf type of decomposition.
The results will very easily enable one to tackle the difficult case of non-identical ser- vers for the evaluation of the queueing delays of delayed customers.
Note the following facts.

The Delayed Customer
In this section, we consider the busy period only and a server's behavior during this period (congestion state).There is generated (for server j) a queueing delay w(j) for j-1,...,s, while the multiserver generates a queueing delay " for an arbitrary delayed customer, in a stationary regime.

The Traffic Per Server
Let T(j) denote the service time of an arbitrary customer in server j (j 1,...,s), in a stationary regime.The termination rate is #j-(1/T(j)).The total termination rate for the multiserver (during a busy period in the steady state) is $ t-Z #j (-), (19)   defining the mean service time T of the multiserver.Let # with In other words, k./ > 0) is defined by the relation (20 The total arrival rate is -[EdN(t)].The traffic intensity handled by the multi- server is A. ' -[EdN(t)]. .Due to the stationary regime and expression (20) for /ti the arrival rate at server j is 1 EdN(t)]. ( The traffic intensity handled by server j is, due to expressions (21) and ( 22), (23 The Stationary G/G/s Queue with Non-Identical Servers 169 This traffic intensity is the same in each server.We may conclude that the following property holds. Property 1 (Behavior of server j for the non-symmetrical G/G/s queue in a sta- tionary regime): Server j behaves as a G/G/1 server, as if an arbitrary arrival is chosen with probability (1/kj) of being handled by server j.It follows that the traf- fic intensity y has the same value in each server during a busy period, in the steady state.
Note the following facts.

1)
The symmetrical G/D/s queue has to be excluded because of a deterministic mechanism for the choice of arrivals.But, when the service times are non- identical, the above property is correct.

2)
For the symmetrical G/G/s queue, we have kj-s, j-1,...,s. (24 We deduce the queueing delay r of an arbitrary delayed customer in the G/G/s queue by applying the expectation operator to expression (9) to evaluate E[e-qr].But, we want to check in the case when kj is valued in set (1, 2,..., s).

Z... Z --1 =1 in expression (18) by
To check this equivalence, we come back to the preceding expressions, related to (17) and (18) and leading to relations (30) through (33) in Le Gall [3], to be applied in the above expression (16) in order to satisfy the factorization method.It follows that the successive residues for (i zi (i s, s -1,..,1)lead to the application of cs_.x to expression (16).The successive residues at the respective poles i-1 z -q-for s,s-1,..., 1, lead to the expression a s_),(q...q), as with the above expression for cs_ .We find again expression (25) or a more general expres- sion corresponding to expression (11).In other words, the singular points are not changed, and the factorization method is still applicable.The new expression for the functional Ee-qr does not change its value, since we know that the solution is unique to determine the G/G/s queue.Now, we shall use the symmetry in expression (17), which leads to the following Expression (16) leads to the expression: l(Zj; j) 1 zj "s-(q"'q)" Finally, considering every variable zj (j > A), we find that (16) reduces to the expres- sion: I I2Ij=A+l(j)l "Os-A(q'''q)' (27) instead of c%_.x(q...q).Property 1 is satisfied and, finally, we may use expression (9).
Property 2 (The distribution of -, the queueing delay of the non-symmetrical G/G/s queue): In a stationary regime, if w(j) is the queueing delay of an arbitrary delayed customer served by server j behaving under Property 1, the queueing delay of an arbitrary delayed customer served by a non-symmetrical (or symmetrical) multi- server is -Min[w(1),..., w(s)], (28) and the expression for the functional Fie-qr] may be deduced from expression (9).From the notation in formula (11) in [3], using Pollaczek's formula in [5], for Re(q) >_ 0, we may write 1 }, with K(() 1 O1( () kj [1(; J) 1], (29 where kj is defined by expression (21).In particular, when ql increases indefinite- ly, we obtain, for the probability of no delay of this server equivalent to server j the expression Qj-Exp -+0 and for the mean busy period size (= mean number of customers served) we deduce the expression 1 nj --, due to some classical relations.In fact, Qj is the probability of initiating a busy period.
a) Case of the Non-Symmetrical GI/G/s Queue: Taking expressions (9), ( 28) and (29) into account, we have that the queueing delay r, of the non-symmetrical GI/G/s queue for an arbitrary delayed customer is given [for Re(q) > 0] by /dZl Jdzs ( -[ -qw(j)) q E[e -qr]-1 (2i)s --i--"" -5-s" Ee Expression (35) below will give a simpler expression for the complementary distribu- tion function of the queueing delay, which will be very convenient for numerical calcu- lations.
b) Case of the Non-Symmetrical G/G/s Queue: For the non-symmetrical G/G/s queue, it is much more difficult to expand the terms between brackets in expression (32), since servers are mutually dependent through the arrival process.Finally, it is simpler to proceed in expressions ( 17) and (18) with the substitutions (15).But the expression in the distribution function of the queueing delay is much more intricate than that of the GI/G/s queue.

The Busy Period
The busy period corresponds to the s servers being busy simultaneously.To evaluate the probability Q0 that the multiserver initiates a busy period, we may note that, among a great number N of successive arrivals, the mean number N. Q0 initiates a busy period of the multiserver and the mean number N.[Qj/kj] corresponds to server j due to Property 1.We deduce the relations: __Qj s 1 Qo j= l-J andln-"kJ' (33) due to expression (31), with n being the mean busy period size of the multiserver, and kj being defined by expression (21).In the case of a symmetrical G/G/s queue, where kjs (see expressions ( 24)), we have nj-n, j-1,...,s. (34)

The Distribution Function
Let W the queueing delay of an arbitrary (delayed or not delayed) customer and let F(t) denote the queueing delay distribution function of a delayed customer.Intro- duce the complementary function G(t)= [1-F(t)] for the multiserver and introduce Gj(t) for server j.
3.4.1 Case of the non-symmetrical G/G/s queue Let 5(t) denote the queueing delay distribution function of an arbitrary (delayed or not delayed) customer of the multiserver.We have (t where P is the probability of delay.As we consider the non-symmetrical GI/G/s queue, expression (28) makes say that since the relation r > t needs simultaneously to have w(j) > t for j 1,..., s.Expres- sion (29) allows us to evaluate Gj().Expressions (a5) and (a6) give for the mo- ments of W: E(Wa) cP.J t a-1.G(t)" dt, 1, 2,...
(37) 0 Gj(t) and this expression (37) can be easily calculated on a computer.
where n is the mean busy period size of the multiserver mean number of customers served during a busy period).Consequently, we may write the following simple expression from expressions (36) and (38) for the complementary distribution func- tion of the delayed customer: G(t)-Exp( for the mean busy period size n of the multiserver.From Pollaczek [5], we may write a.i(0 1 r.E rU" 1 1 F(u; j).du where [-]() denotes the k-fold convolution of the function [-].Finally, (36) gives an intricate expression for the complementary distribution function G(t) of the delayed customer: ii r/ /.
1 Fl(U; j).du which can be evaluated numerically on a computer.need to define the probability of delay P.
To use expression (35) we, now, 4. The Probability of Delay During the busy period (= the congestion state), server's behavior has been defined in a way quite independent of partial occupancy states.For these states, it follows that a busy period appears exactly as a unique congestion state in the lost call model, with n successive service times handled as if there were a unique arrival, with n being the mean value of the busy period size, i.e., of the number of customers served during this busy period.This fact could not be observed with the classical Markovian methods, and it has not been noted in Pollaczek's equation of [6].
With the lost call model in a stationary mode, let Pi Po" h(i) denote the proba- bility that i servers are busy upon the arrival of an arbitrary customer.The probability of loss is Pa P s with 1 =l+l+'"+h(s-l) i-0 With our preceding remark to evaluate the probability of delay P, we have to sub- stitute n. h(s) for h(s). (48) We may now conclude.
Property 3 (The probability of delay)" For a non-symmetrical G/G/s queue in a stationary regime the probability oJ delay is n'Pa (49) P 1 +(n-1).P a where Pa is the probability of loss in the lost call model and n is the mean value of the busy period size as defined by expression (33).
As already seen in [1] and [2], we know that the evaluation of Pa is extremely diffi- cult except in two symmetrical cases: The GI/M/s and M/G/s queues.In particu- lar, for the M/G/s queue we conclude that the delay Erlang formula may be extended for a general service time distribution.In that case, expression (44) gives n [1/(1-r])].From fact 1) after Property 1, we know that the M/D/s queue has to be excluded; however, it has already been noted by C. Palm that the delay Erlang for- mula gives still an excellent approximation. 5. Numerical Comparisons for Some M/G/s Queues Now, we present some numerical comparisons with the results obtained by applying classical Markovian methods to the symmetrical M/E2/s and M/H2/s queues.
In Table 1, we give (for 0.8) our results concerning t(p) and the probability of delay P from expression (49), for s 2, 5, 10, 25, 50 and p 0.5, 0.9, 0.95, 0.99.For the results given by the Markovian methods phase methods), we refer to Table 1 (first part) in Seelen and Tijms [7].These results appeared as approximated (Marko- vian) results in our  1: The symmetrical M/E2/s queue 1) Results a) the conditional queueing delay percentile f(p) for the delayed customer b) probability of delay: P 2) Parameters a) traffic intensity per server: r/= 0.8 b) service time distribution from expression (50): E(T)= 1, C 2 0.5 3) Calculations a) "exact": Section 5.1 b) "Markov': Phase method 5.2 The Symmetrical M/H2/s Queue The service time distribution of any server j is" -blt -b1.b2.(1-r), /2 bl "" 2 b2 r] -/V (bl -b 2 4 2 bl" b2" (1 ). ( 55 As for expression (52), the complementary queueing delay distribution function of any server j is Gl(t The symmetrical M/H2/s queue 1) Results a) the conditional queueing delay percentile t(p) for a delayed customer b) probability of delay: P 2) Parameters a) traffic intensity per server: r/-0.8 b) service time distribution from expression (54): E(T)-1, 3) Calculations a) "exact": Section 5.2 b) "Markov": Phase method Tables 3 (for r/-0.5) and 4 (for -0.9) give some comparisons between values of the mean queueing delay W of an arbitrary (delayed or not delayed) customer, given by calculation and by traffic simulation.Now, the service time distribution of any server j is The Stationary G/G/s Queue wih Non-Identical Servers 177 -blt -b2t Fl(t a I [1 e A-(1 al).Table 3: The symmetrical M/H2/s queue The mean queueing delay W of an arbitrary customer 1) Comparisons between calculations (expression (58)), line C and simulations, line S.
2) Parameters a) traffic intensity per server: t/-0.5 b) service time distribution from expression (57)" E(T) 1, C The mean queueiug delay W of an arbifrary customer 1) Comparisons between calculations (expression (58)), line C and simulations, line S.
[Gl(t)] s. dt, 0 (58 where P is given by expression (49).Tables 3 and 4 consider the cases n-5 and 10.Taking the accuracy of simulations and of calculations of (58) into account, the results given by traffic simulations and by calculations are in good agreement.

Table 1 .
The deviation is not significant.

Table 2 :
2 "e Now, table 2 uses expressions (53) and (56) to give the new values of t(p) and P corresponding to the same values of parameters r/,s and p as in Table1.For the approximated (Markovian) results we refer to Table1(second part) in Seelen and Tijms[7].The deviation is not significant either.