Super-Exponential Solution in Markovian Supermarket Models: Framework and Challenge

Sup er-Exp onen tial Solutio n in Mark o vian Sup ermark et Mo dels: F ramew o rk and Ch allenge Quan-Lin Li Sc ho ol of Economics and Managemen t Sciences Y anshan Univ ersity , Qinh uangdao 06600 4, Ch ina Jan uary 28, 2011 Abstract Marcel F. Neuts op ened a k ey do or in n umerical computation of sto chastic models by means of phase-t ype (PH) distributions and Marko vian ar riv a l pr o cesses (MAPs). T o celebrate his 75th bir thday , this pap er repo rts a more g eneral fra mework of Mar ko- vian sup ermarket mo dels, including a system of differential equations for the fractio n measure and a system of nonlinear eq ua tions for the fixed p o int. T o understa nd t his framework heuristically , this pap er gives a detailed analysis for three imp ortant su- per market e x amples: M/G/1 type, GI/M/1 type and multiple choices, expla ins how to derive the system o f differ e ntial equations by means of density-dep endent jump Marko v process es, and shows that the fixed po int ma y b e simply sup er-ex po nential through solv ing the sys tem o f nonlinear equations. Note that superma rket mo dels are a cla ss of complicated queueing systems and their analys is can no t apply p opular queueing theory , it is necessary in the study of super market models to summarize s uch a more gener al framework which enables us to fo cus on impo rtant research issues. O n this line, this pap er develops matrix- analytical metho ds of Markovian supermar ket mo dels. W e ho p e this will b e able to op en a new av enu e in p erforma nce ev a luation o f sup e rmarket mo dels by mea ns o f matrix-a nalytical metho ds. Keyw ords: Randomized load balancing, sup ermarket model, matrix- analytic method, sup e r-exp onential s olution, densit y-dep endent jump Ma rko v pro cess, Batch Mar ko- vian Arriv al Pro c ess (BMAP), pha s e-type (PH) distribution, fixed p oint. 1 1 In tro duction In the stud y of Marko vian sup ermark et mod els, this pap er pr op oses a more general frame- w ork includin g a system of differen tial equations for the fraction measur e and a system of nonlinear equations for the fixed p oin t, and the b oth systems of equ ations enab le u s to fo cus o n imp ortan t research issues of Marko vian sup ermark et mo d els. A t the same time, this p ap er indicates that it is difficult and c hallenging to analyze the system of differential equations and to solv e the system of nonlinear equations from four key dir ections: Ex- istence of s olution, uniqueness of s olution, stabilit y of s olution and effectiv e algorithms. Since there is a large gap to pr o vide a complete solution to the b oth systems of equations, this p ap er dev otes heuristic un derstanding of how to organize and solv e the b oth systems of equations by means of discuss in g thr ee im p ortan t sup ermark et examples: M/G/1 t yp e, GI/M/1 t yp e and m ultiple c hoices. Sp ecifically , the sup ermark et examp les sho w a key result that the fixed p oin t can b e sup er-exp onent ial for more sup ermark et mo dels. Note that sup ermark et mo dels are a class of complicated queueing sy s tems and th eir analysis can not apply p opu lar queueing theory , while recen t r esearc h ga v e some sim p le and b eau- tiful resu lts f or sp ecial su p ermarket mo dels, e.g., see Mitzenmac her [19], Li and Lu i [11] and Luczak and McDiarmid [14], this motiv ates us in th is p ap er to summarize a more general f ramew ork in order to develo p matrix-analytical metho d s of Marko vian s up erm ar- k et mo dels. W e hop e this is able to op en a n ew a v en ue for p erformance ev aluation of sup er m ark et mo dels by means of matrix-analytical metho d s. Recen tly , a num b er of companies, such as Amazon a nd Go ogle, are offering cloud computing service and cloud m an ufacturing tec h nology . This m otiv ates us in this p ap er to study randomized lo ad balancing f or large-scale net w orks with man y compu tational and manufact uring resources. Rand omized load balancing, wh ere a job is assigned to a serv er f rom a small subset of rand omly c hosen serv ers, is v ery sim p le to implement. It can su rprisingly delive r b etter p erforman ce (for example r educing collisions, w aiting times and bac klogs) in a n u m b er of applica tions includ ing data cen ters, d istributed memory mac h ines, path selection in computer net w orks, and task assignmen t at w eb servers. Su- p ermarket mo d els are extensively us ed to stu dy r andomized load b alancing schemes. In the past ten y ears, sup ermark et mo dels ha v e b een stu died by queueing theory as well as Mark o v p ro cesses. Since the s tu dy of sup ermark et mo d els can not app ly p opular queu eing theory , they ha v e not b een extensive ly stud ied in queu eing committee up to no w. There- 2 fore, this leads to that av ailable queueing results of sup ermarke t mo dels are few up to no w. Some recen t wo rks dealt with the su p ermarke t mo d el w ith Poisson arriv als and ex- p onentia l service times by means of d ensit y-dep endent jum p Mark o v p ro cesses, discu ssed limiting b eha vior of the sup ermark et mo del under a weakly conv ergent s etting wh en the p opulation size go es to infin ite, and indicated th at there exists a doub ly exp onen tial so- lution to th e fixed p oin t thr ough solving the system of nonlinear equations. Readers ma y refer to p opulation p ro cesses by Kurtz [8], and doub ly exp onen tial solution w ith exp onen- tial imp ro v emen t b y Vve densk ay a, Dobrushin and Karp elevic h [27], Mitzenmac her [19], Li and Lui [11] and Luczak and McDiarmid [14 ]. Certain generalization of sup ermark et mo dels has b een explored in, for example, study- ing simple v ariations b y Vvedensk a ya and Su ho v [28 ], Mitzenmac her [20], Azar, Br o der, Karlin and Up fal [1], V¨ oc kin g [26], Mitzenmac h er, Ric ha, and Sitaraman [22] and Li, Lu i and W ang [13]; considerin g non-Po isson arriv als or non-exp onen tial service times by Li, Lui and W ang [12], Li and Lu i [11], Bramson , Lu and Prabhak ar [2 ] an d Li [10]; dis- cussing load information b y Mirc hand aney , T o w sley , and Stanko vic [23], Dahlin [3] and Mitzenmac her [21]; mathematical analysis by Graham [4, 5, 6], Luczak and Norris [16] and Lu czak and McDiarmid [14, 15]; using fast Jac kson n et works by Martin and Suh o v [18], Martin [17] and Su h o v and Vv edensk a y a [25]. The main contributions of the p ap er are t w ofold. The fi rst one is to prop ose a more general framewo rk for Maro vian sup ermark et mo dels. T his framework con tains a system of differentia l equations for th e f raction measur e and a system of nonlinear equations for the fixed p oin t. It is indicated th at there exist m ore d ifficulties and c h allenges for dealing with the system of differential equations and for solving the system of nonlinear equations b ecause of t w o k ey factors: infinite dimension and complicated structure of nonlinear equ ations. Sin ce there is still a large gap u p to b eing able to deal with th e b oth systems of equations systematically , the second con tribution of this pap er is to analyze three imp ortant sup ermark et examples: M/G /1 t yp e, GI/M/1 typ e an d multiple choic es. These examples pr ovide n ecessary u nderstand in g and heuristic methods in order to discuss the b oth systems of equatio ns from practical and m ore general p oin t of view. F or the sup er m ark et examples, this pap er d eriv es the systems of differentia l equations for the fraction measur e by m eans of d ensit y-dep endent jump Mark o v p ro cesses, and illustrates that th e fixed p oin ts ma y b e sup er-exp onent ial through solving the systems of nonlinear equations by means of matrix-analytic metho d s. 3 The remainder of this pap er is organized as follo ws. S ection 2 p r op oses a more general framew ork for Mark ovian sup ermark et mo dels. This f r amew ork conta ins a sys tem of differen tial equations f or the fracti on measure and a system of nonlinear equations f or the fi xed p oin t. In Sections 3 and 4, w e consider a sup ermark et mo del of M/G/1 t yp e b y means of BMAPs and a sup ermarke t mo del of GI/M/1 typ e in terms of batc h P H service p r o cesses, resp ectiv ely . F or the b oth su p ermarket mo dels, we d er ive the sy s tems of d ifferen tial equations satisfied by the fraction measure in terms of densit y-dep end en t jump Mark ov p ro cesses, and obtain the system of nonlinear equations satisfied b y th e fixed p oint whic h is sho wn to b e sup er-exp onentia l. In S ection 5, we analyze t w o su p ermarke t mo dels with multiple choic e n umbers, and giv e su p er-exp onen tial solution to the fixed p oints f or the t w o sup ermark et m o dels. Not e that the su p ermarket examples discuss ed in Sections 3 to 5 can p ro vide a heuristic u nderstand in g for the more general framework of Mark o vian sup ermark et mo del giv en in Section 2. 2 Mark o vian S up ermark et Mo dels In this section, w e prop ose a more general f r amew ork for Marko vian sup ermark et mo d els. This framew ork con tains a system of differentia l equations for th e f raction m easure and a system of n on lin ear equations for the fixed p oin t. Recen t researc h, e.g., see Mitzenmac h er [20] and Li and Lui [11], shows that a Marko - vian su p ermarke t mo d el con tains tw o imp ortan t factors: (1) Contin uous -time Mark o v c hain Q , called sto chastic en vironmen t of the sup ermark et mo del; and (2) Ch oice num b ers, includ ing inpu t c hoice num b ers d 1 , d 2 , . . . , d v and output c hoice n um b ers f 1 , f 2 , . . . , f w . Note that the choic e n um b ers determine d ecomp osed structur e of the sto c hastic en vironmen t Q . W e fi r st analyze sto c hastic environmen t of the Marko vian sup ermark et mo d el. F rom p oint of view of sto c hastic mo dels, we take a more general sto c hastic environmen t whic h is a con tin u ous-time Marko v c hain { X t , t ≥ 0 } with blo ck str u cture. W e assume that the Mark o v c hain { X t , t ≥ 0 } on state space Ω = { ( k , j ) : k ≥ 0 , 1 ≤ j ≤ m k } is irreducible 4 and p ositiv e recurr en t, and that its infinitesimal generator is giv en b y Q =         Q 0 , 0 Q 0 , 1 Q 0 , 2 · · · Q 1 , 0 Q 1 , 1 Q 1 , 2 · · · Q 2 , 0 Q 2 , 1 Q 2 , 2 · · · . . . . . . . . . . . .         , (1) where Q i,j is a matrix of size m i × m j whose ( r , r ′ )th en try is the transition rate of the Mark o v chain from state ( i, r ) to state ( j, r ′ ). It is well- kno wn that Q i,j ≥ 0 for i 6 = j , Q i,i is inv ertible with strictly n egativ e diagonal en tries an d nonnegativ e off-diagonal en tries. F or state ( i, k ) , i is called the level variable and k the phase variable . W e write lev el i as L i = { ( i, k ) : 1 ≤ k ≤ m i } . Since the Mark o v chain is ir reducible, for eac h lev el i there m ust exist at east one left- blo c k state transition: ← lev el i or lev el i ← , and at east one r igh t-block state transition: → leve l i or lev el i → . W e write E left = {← leve l i or lev el i ← : lev el i ∈ Ω } and E right = {→ lev el i or leve l i → : lev el i ∈ Ω } . Note that E left and E right describ e output and input pr o cesses in the sup ermark et mo d el. Based on the tw o blo ck-t ransition sets E left and E right , we write Q = Q left + Q right , (2) and f or i ≥ 0 Q i,i = Q left i + Q right i . Th us we hav e Q left =         Q left 0 Q 1 , 0 Q left 1 Q 2 , 0 Q 2 , 1 Q left 2 . . . . . . . . . . . .         and Q right =         Q right 0 Q 0 , 1 Q 0 , 2 · · · Q right 1 Q 1 , 2 · · · Q right 2 · · · . . .         . 5 Note that Qe = 0 , Q left e = 0 and Q right e = 0, where e is a column v ector of ones w ith a suitable dimension in the con text. W e assu me that the m atrices Q left j for j ≥ 1 and Q right i for i ≥ 0 are all inv ertible, w h ile Q left 0 is p ossibly singular if there is not an output pro cess in lev el 0. W e call Q = Q left + Q right an inpu t-output rate d ecomp osition of the Marko vian sup er m ark et mo del. No w , w e p r o vide a choi ce decomp osition of the Marko vian su p ermarket mod el through decomp osing the t w o matrices Q left and Q right . Note th at th e choice decomp osition is based on the input choic e num b ers d 1 , d 2 , . . . , d v and the outp ut c hoice num b ers f 1 , f 2 , . . . , f w . W e write Q left = Q left ( f 1 ) + Q left ( f 2 ) + · · · + Q left ( f w ) (3) for the outp ut choic e num b ers f 1 , f 2 , . . . , f w , and Q right = Q right ( d 1 ) + Q right ( d 2 ) + · · · + Q right ( d v ) (4) for the in put c hoice num b ers d 1 , d 2 , . . . , d v . T o study the Mark ovia n sup ermark et mo d el, we need to introd u ce t wo v ector n otation. F or a v ector a = ( a 1 , a 2 , a 3 , . . . ), we write a ⊙ d =  a d 1 , a d 2 , a d 3 , . . .  and a ⊙ 1 d =  a 1 d 1 , a 1 d 2 , a 1 d 3 , . . .  . Let S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ) b e the fraction measure of the Marko vian su- p ermarket mo del, where S i ( t ) is a row v ector of size m i for i ≥ 0. Th en S ( t ) ≥ 0 and S 0 ( t ) e = 1. Base d on the inpu t-output r ate decomp osition and th e c hoice decomp osition for the sto c hastic environmen t, w e introdu ce the follo wing sys tem of differenti al equations satisfied by the fraction measure S ( t ) as f ollo ws: S 0 ( t ) ≥ 0 and S 0 ( t ) e = 1 , (5) and d d t S ( t ) = w X l =1 S ⊙ f l ( t ) Q left ( f l ) + v X k =1 S ⊙ d k ( t ) Q right ( d k ) . (6) In the Mark o vian sup ermark et mo del, a row v ector π = ( π 0 , π 1 , π 2 , . . . ) is called a fi xed p oint of th e fraction measure S ( t ) if lim t → + ∞ S ( t ) = π . I n this case, it is easy to see that lim t → + ∞  d d t S ( t )  = 0 . 6 If there exists a fi xed p oin t of the fraction measure, then it follo ws from (5) an d (6) that the fixed p oin t is a nonnegativ e non-zero solution to the follo wing sy s tem of nonlinear equations π 0 ≥ 0 and π 0 e = 1 , (7) and w X l =1 π ⊙ f l Q left ( f l ) + v X k =1 π ⊙ d k Q right ( d k ) = 0 . (8) Remark 1 If d k = 1 for 1 ≤ k ≤ v and f l = 1 for 1 ≤ l ≤ w , then the system of differ ential e qu ations (5) and (6) is given by S 0 ( t ) ≥ 0 and S 0 ( t ) e = 1 , and d d t S ( t ) = S ( t ) Q. Thus we obtain S ( t ) = cS ( 0) exp { Qt } . L et W ( t ) = ( W 0 ( t ) , W 1 ( t ) , W 2 ( t ) , . . . ) = S (0) exp { Qt } , wher e W i ( t ) is a r ow ve ctor of size m i for i ≥ 0 . Then S ( t ) = cW ( t ) , wher e c = 1 /W 0 ( t ) e . At the same time, the system of nonline ar e quations (7) and (8) is given by π 0 ≥ 0 and π 0 e = 1 , and π Q = 0 . L et W = ( w 0 , w 1 , w 2 , . . . ) b e the stationary pr ob ability v e ctor of the Markov chain Q , wher e W i is a r ow ve ctor of size m i for i ≥ 0 . Then π = cW , wher e c = 1 /w 0 e . Note that the stationary pr ob ability ve ctor W of the blo c k -structur e d Markov chain Q is g i ven a detaile d analysis i n Chapter 2 of Li [9] by me ans of the RG-factor izations. If there exist some d k ≥ 2 or/and f l ≥ 2 in th e Mark o v ian s up erm ark et mo del, then the system of differentia l equations (5) and (6) and th e system of nonlinear equations (7) to (8) are tw o deco mp osed p o w er -form generaliza tions of tr ansien t solution and of s ta- tionary p robabilit y of an irredu cible contin uous-time Mark o v c h ain with b lo c k structure 7 (see Chapters 2 and 8 of Li [9 ]). Note that Li [9] can deal w ith transien t solution and stationary probabilit y for an irr educible blo ck-structured Mark o v c hain, where the RG- factorizat ons pla y a key role. How eve r, the R G-factoriza tons can not hold for Mark o vian sup er m ark et mo d els with some d k ≥ 2 or/and f l ≥ 2. Ther efore, there exist more dif- ficulties and challe nges to stu d y the system of d ifferen tial equations (5) and (6) and the system of nonlin ear equations (7) to (8). Sp ecifically , it still ke eps not to b e able to answe r four im p ortan t issues: Existence of solution, uniqueness of solution, stabilit y of solution and effectiv e algorithms. Th is is similar to some researc h on the four imp ortan t issu es of irreducible con tin u ous-time Mark o v c hains with blo c k structure. In the remainder of this pap er, we will stud y thee imp ortant Marko vian sup ermark et examples: M/G/1 t yp e, GI/M/1 t y p e and multiple c hoices. O ur p urp ose is to p r o vide heuristic und erstanding of h o w to set u p and solve the system of differenti al equations (5) and (6 ), and the system of nonlinear equations (7) to (8). 3 A Sup ermark et M o del of M/G/1 Typ e In this section, we consider a sup ermark et mo del with a BMAP and exp on ential service times. Note that the sto c hastic en vironmen t is a Mark o v c hain of M/G/1 t yp e, the sup er m ark et mo del is called to b e of M/G/1 type. F or the sup ermark et m o del of M/G/1 t yp e, we set up the s ystem of differential equations for the fraction measur e by means of densit y-dep endent jump Marko v pro cesses, and derive the system of nonlinear equations satisfied by the fixed p oin t which is sho wn to b e sup er-exp onen tial solution. The su p ermarket mo del of M/G/1 type is describ ed as f ollo ws. Customers arriv e at a q u eueing system of n > 1 serv ers as a BM AP with irreducible matrix descriptor ( nC, nD 1 , nD 2 , nD 3 , . . . ) of size m , where the matrix C is inv ertible and has strictly nega- tiv e diagonal entries and nonnegativ e off-diagonal; D k ≥ 0 is the arriv al rate matrix with batc h size k for k ≥ 1. W e assume that P ∞ k =1 k D k is fin ite and th at C + P ∞ k =1 D k is an irreducible infinitesimal generator with ( C + P ∞ k =1 D k ) e = 0. Let γ b e the stationary probabilit y vecto r of the irreducible Ma rk o v c hain C + P ∞ k =1 D k . Then the sta tionary arriv al r ate of the BMAP is giv en by nλ = n γ P ∞ k =1 k D k e . The service times of eac h cus- tomer are exp onen tiall y distribu ted with service rate µ . Eac h batch of arriving cu stomers c ho ose d ≥ 1 serv ers indep endently and uniformly at random from the n ser vers, and joins for service at the serv er whic h currently p ossesses the fewest num b er of cu stomers. 8 If there is a tie, s erv ers with the few est num b er of cus tomers will b e c h osen rand omly . All customers in ev ery server will b e serve d in the fir st-come-first service (F CFS) manner. W e assume that all the random v ariables defined ab o v e are indep enden t of eac h other and that this system is op erating in the region ρ = λ/µ < 1. Clearly , d is an input c hoice n um b er in this sup ermark et mo del. Figure 1 is d epicted as an illustration for su p ermarket mo dels of M/G/1 type. 2 1 Server 1 Server 2 Server 3 Server n BMAP Each batch of customers Probe servers 1 2 ( , , , ) nC nD nD ! d d 1 d  Figure 1: A sup ermark et m o del of M/G/1 type The sup ermark et mo del with a BMAP and exp onen tial service times is stable if ρ = λ/µ < 1. Th is pro of can b e giv en b y a simple comparison argumen t with the qu eueing system in whic h eac h customer queues at a random serv er (i.e. , where d = 1). When d = 1, eac h serv er acts like a BMAP/M/1 queue wh ic h is stable if ρ = λ/µ < 1, see c hapter 5 in Neuts [24]. Similar to analysis in Winston [30] and W eb er [29], th e comparison argument leads to t w o u seful results: (1) the sh ortest queue is optimal due to the assumptions on a BMAP and exp onential service times in th e sup ermark et m o del; and (2) th e size of the longest q u eue in the sup ermarket m o del is sto c hasticall y domin ated by the size of the longest queue in a set of n indep enden t BMAP/M/1 q u eues. W e defi n e n ( i ) k ( t ) as the n umber of queues with at least k cus tomers, includ ing cus- tomers in serv ice, and w ith the BMAP in ph ase i at time t ≥ 0. Clearly , 0 ≤ n ( i ) k ( t ) ≤ n for k ≥ 0 and 1 ≤ i ≤ m . Let x ( i ) n ( k , t ) = n ( i ) k ( t ) n , whic h is the fraction of queues with at least k customers and the BMAP in ph ase i at 9 time t ≥ 0 for k ≥ 0. W e write X n ( k , t ) =  x (1) n ( k , t ) , x (2) n ( k , t ) , . . . , x ( m ) n ( k , t )  for k ≥ 0, and X n ( t ) = ( X n (0 , t ) , X n (1 , t ) , X n (2 , t ) , . . . ) . The state of the sup ermarket mo del m a y b e describ ed b y the v ecto r X n ( t ) for t ≥ 0. Since the arr iv al pro cess to the queueing system is a BMAP and the service time of eac h customer is exp onential, the sto chastic pro cess { X n ( t ) , t ≥ 0 } is a Mark o v pro cess whose state space is giv en by Ω n = {  g (0) n , g (1) n , g (2) n . . .  : g (0) n is a pr obabilit y v ector, g ( k ) n ≥ g ( k +1) n ≥ 0 for k ≥ 1, and ng ( l ) n is a vec tor of nonn egativ e intege rs for l ≥ 0 } . Let s ( i ) k ( n, t ) = E h x ( i ) k ( n, t ) i , and S k ( n, t ) =  s (1) k ( n, t ) , s (2) k ( n, t ) , . . . , s ( m ) k ( n, t )  for k ≥ 0 , S ( n, t ) = ( S 0 ( n, t ) , S 1 ( n, t ) , S 2 ( n, t ) , . . . ) . As sho wn in Martin and Su ho v [18] and Luczak and McDiarmid [14], the Mark o v pr o cess { X n ( t ) , t ≥ 0 } is asymp totically deterministic as n → ∞ . Th us lim n →∞ E h x ( i ) k ( n, t ) i alw a ys exist by means of the la w of large n umbers for k ≥ 0. Based on this, w e w rite S k ( t ) = lim n →∞ S k ( n, t ) for k ≥ 0, and S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ) . Let X ( t ) = lim n →∞ X n ( t ). Then it is easy to see from th e BMAP and the exp onen tial service times that { X ( t ) , t ≥ 0 } is also a Mark o v pr o cess whose state space is give n by Ω = n g (0) , g (1) , g (2) , . . .  : g (0) is a pr obabilit y v ector , g ( k ) ≥ g ( k +1) ≥ 0 for k ≥ 1 o . If the initial distribution of the Mark o v pro cess { X n ( t ) , t ≥ 0 } approac h es the Dirac delta- measure concen trated at a p oin t g ∈ Ω , th en X ( t ) = lim n →∞ X n ( t ) is concentrat ed on 10 the tra jectory S g = { S ( t ) : t ≥ 0 } . This indicates a law of large n umbers for the time ev olution of the fraction of queues of differen t lengths. F urthermore, the Mark o v p ro cess { X n ( t ) , t ≥ 0 } con v erges weakly to th e fraction vec tor S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ) as n → ∞ , or for a s u fficien tly sm all ε > 0, lim n →∞ P {|| X n ( t ) − S ( t ) || ≥ ε } = 0 , where || a || is th e L ∞ -norm of ve ctor a . In w h at follo ws w e set up a system of differentia l v ector equations satisfied b y the fraction vecto r S ( t ) by means of d ensit y-dep endent jump Marko v pro cesses. W e first p r o vide an example to indicate ho w to derive the different ial v ector equations. Consider the sup ermark et mo del with n servers, and determine the exp ected c hange in the n um b er of queues w ith at least k customers o v er a small time int erv al [0 , dt ). Th e probabilit y ve ctor that an arriving customer joins a qu eue with k − 1 customers in th is time interv al is give n by h S ⊙ d 0 ( n, t ) D k + S ⊙ d 1 ( n, t ) D k − 1 + · · · + S ⊙ d k − 1 ( n, t ) D 1 + S ⊙ d k ( n, t ) C i · n d t, since eac h arriving customer chooses d servers indep end en tly and uniformly at r andom from the n serv ers, and wait s for service at th e server wh ic h currently con tains the fewe st n um b er of customers. Similarly , the p r obabilit y v ector that a cus tomer lea v es a serv er queued b y k customers in this time inte rv al is giv en by [ − µS k ( n, t ) + µS k +1 ( n, t )] · n d t. Therefore, we obtain d E [ n k ( n, t )] = " k − 1 X l =0 S ⊙ d l ( n, t ) D k − l + S ⊙ d k ( n, t ) C # · n d t + [ − µS k ( n, t ) + µS k +1 ( n, t )] · n d t. This leads to d S k ( n, t ) d t = k − 1 X l =0 S ⊙ d l ( n, t ) D k − l + S ⊙ d k ( n, t ) C + − µS k ( n, t ) + µS k +1 ( n, t ) . (9) Since lim n →∞ E h x ( i ) k ( n, t ) i alw a ys exists for k ≥ 0, taking n → ∞ in b oth sides of Equation (9) we can easily obtain d S k ( t ) d t = k − 1 X l =0 S ⊙ d l ( t ) D k − l + S ⊙ d k ( t ) C − µS k ( t ) + µS k +1 ( t ) . (10) 11 Using a similar analysis to that in Equation (10), w e obtain the system of d ifferen tial v ector equ ations for the fr action vec tor S ( t ) = ( S 0 ( t ) , S 1 ( t ) , . . . ) as follo ws: S 0 ( t ) ≥ 0 , S 0 ( t ) e = 1 , (11) d d t S 0 ( t ) = S ⊙ d 0 ( t ) C + µS 1 ( t ) (12) and f or k ≥ 1 d d t S k ( t ) = k − 1 X l =0 S ⊙ d l ( n, t ) D k − l + S ⊙ d k ( n, t ) C − µS k ( n, t ) + µS k +1 ( n, t ) . (13) Let π b e the fi xed p oin t. Then π satisfies the follo wing sys tem of nonlinear equations π 0 ≥ 0 , π 0 e = 1 , (14) π ⊙ d 0 C + µπ 1 = 0 (15) and f or k ≥ 1, k − 1 X l =0 π ⊙ d l D k − l + π ⊙ d k C − µπ k + µ π k +1 = 0 . (16) Let Q right =            C D 1 D 2 D 3 D 4 · · · C D 1 D 2 D 3 · · · C D 1 D 2 · · · C D 1 · · · . . .            and Q left =            0 µI − µI µI − µI µI − µI . . . . . .            . Then the sys tem of differentia l vec tor equations is giv en by S 0 ( t ) ≥ 0 , S 0 ( t ) e = 1 , and d d t S ( t ) = S ⊙ d ( t ) Q right + S ( t ) Q left ; 12 and th e system of nonlinear equations (14) to (16) is given by π 0 ≥ 0 , π 0 e = 1 , and π ⊙ d Q right + π Q left = 0 . Remark 2 F or the sup ermarket mo del with a BMA P and exp onential servic e times, its sto chastic envir onment is a M arkov chain of M/G/1 typ e whose infinitesimal gener ator is given by Q = Q left + Q right . This example cle arly indic ates how to set u p the system of diffe r ential e quations (5) and (6) for the fr action me asur e and the system of nonline ar e quations (7) to (8) f or the fixe d p oint. In the remainder of this section, we pro vide a sup er-exp onen tial solution to the fixed p oint π by m eans of some usefu l relations among the v ectors π k for k ≥ 0. It follo w s from (16) that  π ⊙ d 1 , π ⊙ d 2 , π ⊙ d 3 , . . .          C D 1 D 2 · · · C D 1 · · · C · · · . . .         + ( π 1 , π 2 , π 3 , . . . )         − µI µI − µI µI − µI . . . . . .         = −  π ⊙ d 0 D 1 , π ⊙ d 0 D 2 , π ⊙ d 0 D 3 , . . .  . (17) Let A =         − µI µI − µI µI − µI . . . . . .         . Then ( − A ) − 1 =         1 µ I 1 µ I 1 µ I 1 µ I 1 µ I 1 µ I . . . . . . . . . . . .         . 13 Note that         D 0 D 1 D 2 · · · D 0 D 1 · · · D 0 · · · . . .          − A − 1  =            1 µ ∞ P k =0 D k 1 µ ∞ P k =1 D k 1 µ ∞ P k =2 D k · · · 1 µ ∞ P k =0 D k 1 µ ∞ P k =0 D k 1 µ ∞ P k =1 D k · · · 1 µ ∞ P k =0 D k 1 µ ∞ P k =0 D k 1 µ ∞ P k =0 D k · · · . . . . . . . . .            and  π ⊙ d 0 D 1 , π ⊙ d 0 D 2 , π ⊙ d 0 D 3 , . . .   − A − 1  = 1 µ π ⊙ d 0 ∞ X k =1 D k , 1 µ π ⊙ d 0 ∞ X k =2 D k , 1 µ π ⊙ d 0 ∞ X k =3 D k , . . . ! , it f ollo ws from (17 ) that π 1 = π ⊙ d 0 " 1 µ ∞ X i =1 D i # + ∞ X j =1 π ⊙ d j " 1 µ ∞ X i =0 D i # (18) and f or k ≥ 2, π k = k − 1 X i =0 π ⊙ d i   1 µ ∞ X j = k − i D j   + ∞ X j = k π ⊙ d j " 1 µ ∞ X i =0 D i # . (19) T o omit the terms P ∞ j = k π ⊙ d j  1 µ ∞ P i =0 D i  for k ≥ 1, we assu m e that the system of n onlinear equations (18) and (19) h as a closed-form solution π k = r ( k ) γ ⊙ 1 d , (20) where r ( k ) is an und erdetermined p ositiv e constan t for k ≥ 1. Th en it follo ws f rom (18), (19) and (20) that π 1 = π ⊙ d 0 " 1 µ ∞ X i =1 D i # (21) or r (1) γ ⊙ 1 d = π ⊙ d 0 " 1 µ ∞ X i =1 D i # ; (22) and f or k ≥ 2, π k = k − 1 X i =0 π ⊙ d i   1 µ ∞ X j = k − i D j   or r ( k ) γ ⊙ 1 d = π ⊙ d 0   1 µ ∞ X j = k D j   + k − 1 X i =0 [ r ( i )] d γ   1 µ ∞ X j = k − i D j   . (23) 14 Let θ = 1 /γ ⊙ 1 d e . Then 0 < θ < 1. Let λ k = γ P ∞ i = k D i e an d ρ k = λ k /µ . Then it follo ws from (22) and (23) that r (1) = θ µ π ⊙ d 0 ∞ X i =1 D i e (24) and f or k ≥ 2 r ( k ) = θ µ π ⊙ d 0 ∞ X j = k D j e + θ µ k − 1 X i =1 [ r ( i )] d γ ∞ X j = k − i D j e = θ µ π ⊙ d 0 ∞ X j = k D j e + θ k − 1 X i =1 [ r ( i )] d ρ k − i . (25) It is easy to s ee fr om (24) and (25) that π 0 and r (1) are t w o k ey underd etermined terms for the closed-form solution to the sy s tem of nonlinear equations (22) and (23). Let us first deriv e the vec tor π 0 . It follo ws from (15) and (21) that      π ⊙ d 0 C + µπ 1 = 0 , µπ 1 = π ⊙ d 0 ∞ P i =1 D i . This leads to π ⊙ d 0 C + ∞ X i =1 D i ! = 0 . Th us, it is easy to s ee that π 0 = θ γ ⊙ 1 d , wh ich is a probability vect or with π 0 e = 1. Hence w e h a v e π 1 = − θ d µ γ C = θ d µ γ ∞ X i =1 D i . (26) It follo ws fr om (24) and (25) th at r (1) = θ µ · θ d γ ∞ X i =1 D i e = θ d +1 ρ 1 (27) and f or k ≥ 2 r ( k ) = θ µ π ⊙ d 0 ∞ X j = k D j e + θ µ k − 1 X i =1 [ r ( i )] d γ ∞ X j = k − i D j e = θ d +1 ρ k + θ k − 1 X i =1 [ r ( i )] d ρ k − i . (28) Therefore, we obtain the su p er-exp onentia l solution to the fixed p oint as f ollo ws: π 0 = θ γ ⊙ 1 d 15 and f or k ≥ 1 π k = " θ d +1 ρ k + θ k − 1 X i =1 [ r ( i )] d ρ k − i # γ ⊙ 1 d . 4 A Sup ermark et M o del of GI/M/1 Typ e In this sectio n, we analyze a sup ermark et mo d el with Poisson arriv als and batc h PH service pro cesses. Note that the sto c hastic environmen t is a Mark o v c hain of GI/M/1 typ e, thus the sup erm ark et mo d el is calle d to be of GI/M/1 t yp e. F or the sup ermarke t m o del of GI/M/1 t yp e, w e set u p th e system of differen tial equations f or the fraction m easure b y means of dens ity-dep en den t jum p Mark o v pro cesses, and deriv e the s ystem of nonlinear equations satisfied the fi xed p oin t whic h can b e computed b y an iterativ e algorithm. F urther, it is seen that the sup ermark et mo del of GI/M/1 t yp e is more d iffi cu lt than the case of M/G/1 t yp e. Let us describ e th e sup ermark et mo del of GI/M/1 typ e. Cu stomers arrive at a queu e- ing system of n > 1 serv ers as a Poisso n p ro cess with arr iv al r ate nλ for λ > 0. The service times of eac h batc h of customers are of phase t yp e with irr ed ucible representa tion ( α, T ) of order m and with a batc h s ize distribution { b k , k = 1 , 2 , 3 , . . . } for P ∞ k =1 b k = 1 and b = P ∞ k =1 k b k < + ∞ . Let T 0 = − T e  0. Th en the exp ected service time is give n b y 1 /µ = − αT − 1 e = η T 0 , where η is the stationary p robabilit y vec tor of the Marko v c hain T + T 0 α . E ac h batch of arriving customers c ho ose d ≥ 1 s erv ers in dep end en tly and uniformly at random from the n serv ers, and waits for service at the server w hic h currentl y con tains the fewest n u m b er of customers. If there is a tie, servers with the few est n um b er of cu s tomers will b e chosen randomly . All customers in ev ery server will b e serv ed in FCFS for different batc hes and in rand om service within one b atc h. W e assu me that all random v ariables defined ab o ve are ind ep endent of eac h ot her, and that the system is op erating in the stable region ρ = λ/µ b < 1. Clearly , d is an inp ut choic e num b er in this su p ermarket mo d el. Figure 2 is d epicted as an illustration for su p erm arket mo d els of GI/M/1 type. W e define n ( i ) k ( t ) as the num b er of queues with at least k customers and the service time in ph ase i at time t ≥ 0. Clearly , 0 ≤ n ( i ) k ( t ) ≤ n for k ≥ 1 and 1 ≤ i ≤ m . Let X (0) n ( t ) = n n = 1 , 16 2 1 Server 1 Server 2 Server 3 Server n Poisson arrivals Each customer Probes servers n O d d 1 d  Batch PH service Batch PH service Batch PH service Batch PH service Figure 2: A s up erm ark et mo del of GI/M/1 typ e and k ≥ 1 X ( k ,i ) n ( t ) = n ( i ) k ( t ) n , whic h is the fraction of qu eues with at least k customers and the serv ice time in phase i at time t ≥ 0. W e write X ( k ) n ( t ) =  X ( k, 1) n ( t ) , X ( k , 2) n ( t ) , . . . , X ( k ,m ) n ( t )  , k ≥ 1 , X n ( t ) =  X (0) n ( t ) , X (1) n ( t ) , X (2) n ( t ) , . . .  . The state of the sup ermark et mo d el ma y b e describ ed b y the v ecto r X n ( t ) for t ≥ 0. Since the arriv al pro cess to the queu eing system is P oisson and the service times of eac h server are of ph ase t yp e, { X n ( t ) , t ≥ 0 } is a Mark o v pro cess whose s tate space is giv en by Ω n = {  g (0) n , g (1) n , , g (2) n . . .  : g (0) n = 1 , g ( k − 1) n ≥ g ( k ) n ≥ 0 , and ng ( k ) n is a vec tor of nonnegativ e in tegers for k ≥ 1 } . Let s 0 ( n, t ) = E h X (0) n ( t ) i and k ≥ 1 s ( i ) k ( n, t ) = E h X ( k ,i ) n ( t ) i . Clearly , s 0 ( n, t ) = 1. W e write S k ( n, t ) =  s (1) k ( n, t ) , s (2) k ( n, t ) , . . . , s ( m ) k ( n, t )  , k ≥ 1 . 17 As sh o wn in Martin and S uho v [18] an d Luczak and McDiarmid [14], the Mark o v pro- cess { X n ( t ) , t ≥ 0 } is asymptotically deterministic as n → ∞ . Thus lim n →∞ E h X (0) n ( t ) i and lim n →∞ E h X ( k ,i ) n i alw a ys exist by m eans of the law of large num b ers. Based on this, w e w rite S 0 ( t ) = lim n →∞ s 0 ( n, t ) = 1 , for k ≥ 1 s ( i ) k ( t ) = lim n →∞ s ( i ) k ( n, t ) , S k ( t ) =  s (1) k ( t ) , s (2) k ( t ) , . . . , s ( m ) k ( t )  and S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ) . Let X ( t ) = lim n →∞ X n ( t ). Th en it is easy to see from P oisson arr iv als and batc h PH service times that { X ( t ) , t ≥ 0 } is also a Mark o v pr o cess whose state space is give n by Ω = n g (0) , g (1) , g (2) , . . .  : g (0) = 1 , g ( k − 1) ≥ g ( k ) ≥ 0 o . If the initial distribution of the Marko v pro cess { X n ( t ) , t ≥ 0 } approac hes the Dirac delta-measure concen trated at a p oint g ∈ Ω, then its steady-state distr ibution is con- cen trated in the limit on the tra jectory S g = { S ( t ) : t ≥ 0 } . Th is indicates a law of large n um b ers for the time evol ution of the fraction of qu eues of differen t lengths. F ur- thermore, the Marko v pro cess { X n ( t ) , t ≥ 0 } conv erges wea kly to the fr action vect or S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ), or for a sufficiently small ε > 0, lim n →∞ P {|| X n ( t ) − S ( t ) || ≥ ε } = 0 , where || a || is th e L ∞ -norm of ve ctor a . T o d etermine the fraction ve ctor S ( t ) , we n eed to set u p a system of differenti al v ector equations satisfied b y the fr action measure S ( t ) by means of densit y-dep endent jump Mark o v pro cesses. Consid er the s up erm ark et m o del with n servers, and determine the exp ected c hange in the num b er of queues with at le ast k customers o ver a small time p erio d of length d t . Th e prob ab ility v ector that dur in g this time p erio d, any arriving customer joins a qu eue of size k − 1 is giv en b y n h λS ⊙ d k − 1 ( n, t ) − λS ⊙ d k ( n, t ) i d t. 18 Similarly , the p robabilit y vecto r th at a cus tomer lea ves a s er ver qu eued by k customers is giv en b y n " S k ( n, t ) T + ∞ X l =1 b l S k + l ( n, t ) T 0 α # d t. Therefore, we obtain d E [ n k ( n, t )] = n h λS ⊙ d k − 1 ( n, t ) − λS ⊙ d k ( n, t ) i d t + n " S k ( n, t ) T + ∞ X l =1 b l S k + l ( n, t ) T 0 α # d t. This leads to d S k ( n, t ) d t = λS ⊙ d k − 1 ( n, t ) − λS ⊙ d k ( n, t ) + S k ( n, t ) T + ∞ X l =1 b l S k + l ( n, t ) T 0 α. (29) T aking n → ∞ in b oth sides of Equation (29), we hav e d S k ( t ) d t = λS ⊙ d k − 1 ( t ) − λS ⊙ d k ( t ) + S k ( t ) T + ∞ X l =1 b l S k + l ( n, t ) T 0 α. (30) Using a similar analysis to Equation (30), we obtain a system of differen tial v ector equations f or the fr action vec tor S ( t ) = ( S 0 ( t ) , S 1 ( t ) , S 2 ( t ) , . . . ) as follo ws: S 0 ( t ) = 1 , (31) d d t S 0 ( t ) = − λS d 0 ( t ) + ∞ X l =1 S l ( t ) T 0 ∞ X k = l b k , (32) d d t S 1 ( t ) = λαS d 0 ( t ) − λS ⊙ d 1 ( t ) + S 1 ( t ) T + ∞ X l =1 b l S 1+ l ( t ) T 0 α, (33) and f or k ≥ 2, d d t S k ( t ) = λS ⊙ d k − 1 ( t ) − λS ⊙ d k ( t ) + S k ( t ) T + ∞ X l =1 b l S k + l ( t ) T 0 α. (34) If the r o w v ector π = ( π 0 , π 1 , π 2 , . . . ) is a fixed p oin t of the fraction ve ctor S ( t ), then the fixed p oin t π satisfies th e follo wing system of nonlinear equations π 0 = 1 (35) − λπ d 0 + ∞ X l =1 π l T 0 ∞ X k = l b k = 0 , (36) 19 λαπ d 0 − λπ ⊙ d 1 + π 1 T + ∞ X l =1 b l π 1+ l T 0 α = 0 , (3 7) and f or k ≥ 2, λπ ⊙ d k − 1 − λπ ⊙ d k + π k T + ∞ X l =1 b l π k + l T 0 α = 0 . (38) Let Q right =            − λ λα · · · − λI λI · · · − λI λI · · · − λI λI · · · . . . . . .            and Q left =                  0 T 0 T T 0 ∞ P k =2 b k b 1 T 0 α T T 0 ∞ P k =3 b k b 2 T 0 α b 1 T 0 α T T 0 ∞ P k =4 b k b 3 T 0 α b 2 T 0 α b 1 T 0 α T . . . . . . . . . . . . . . . . . .                  . Then the sys tem of differentia l vec tor equations for the fraction measure is give n by S 0 ( t ) = 1 , and d d t S ( t ) = S ⊙ d ( t ) Q right + S ( t ) Q left ; and th e system of nonlinear equations for the fixed p oint is given by π 0 = 1 , and π ⊙ d Q right + π Q left = 0 . (39) In the remainder of this section, we pr o vide an iterativ e algorithm for computing the fixed p oin t for the s u p ermarket mo d el of GI/M/1 t yp e. Sp ecifically , the iterativ e algorithm indicates that the sup ermark et mo del of GI/M/1 type is more diffi cu lt than the case of M/G/1 type. 20 Let B =         − λI λI · · · − λI λI · · · − λI · · · . . .         and Q service =            T b 1 T 0 α T b 2 T 0 α b 1 T 0 α T b 3 T 0 α b 2 T 0 α b 1 T 0 α T . . . . . . . . . . . . . . .            Then it follo ws fr om (39) that π d 0 ( λα, 0 , 0 , 0 , . . . ) + π ⊙ d L B + π L Q Service = 0 , (40) where π L = ( π 1 , π 2 , π 3 , . . . ). Note that ( − B ) − 1 =         1 λ I 1 λ I 1 λ I · · · 1 λ I 1 λ I · · · 1 λ I · · · . . .         , using π 0 = 1 we obtain π d 0 ( λα, 0 , 0 , 0 , . . . ) ( − B ) − 1 = ( α, α, α, . . . ) and Q service ( − B ) − 1 = 1 λ              T T T T · · ·  T 0 α  b 1 T +  T 0 α  b 1 T +  T 0 α  b 1 T +  T 0 α  b 1 · · ·  T 0 α  b 2  T 0 α  2 P k =1 b k T +  T 0 α  2 P k =1 b k T +  T 0 α  2 P k =1 b k · · ·  T 0 α  b 3  T 0 α  3 P k =2 b k  T 0 α  3 P k =1 b k T +  T 0 α  3 P k =1 b k · · · . . . . . . . . . . . .              , Th us it follo ws fr om (40) that π ⊙ d L = π d 0 ( λα, 0 , 0 , 0 , . . . ) ( − B ) − 1 + π L Q service ( − B ) − 1 = ( α, α, α, . . . ) + π L Q service ( − B ) − 1 . 21 This leads to that for k ≥ 1 λπ d k = λα + k +1 X l =1 π l T + b 1 k +2 X l =2 π l  T 0 α  + b 2 k +3 X l =3 π l  T 0 α  + · · · . (41) T o solve the system of nonlinear equations (41) for k ≥ 1, w e assume th at the fixed p oint has a closed-form solution π k = r ( k ) η , where η is the stationary p robabilit y vecto r of the Mark o v c hain T + T 0 α . It follo ws from (41) th at λr d ( k ) η ⊙ d = λα + k +1 X l =1 r ( l ) η T + b 1 k +2 X l =2 r ( l ) η  T 0 α  + b 2 k +3 X l =3 r ( l ) η  T 0 α  + · · · . T aking θ = η ⊙ d e . Th en θ ∈ (0 , 1) . Noting th at αe = 1 , η T e = − µ and η T 0 = µ , we obtain that f or k ≥ 1 ρθ r d ( k ) = ρ − k +1 X l =1 r ( l ) + b 1 k +2 X l =2 r ( l ) + b 2 k +3 X l =3 r ( l ) + · · · . (42) This giv es ρθ h r d ( k ) − r d ( k + 1) i = r ( k + 2) − ∞ X l =1 b l r ( k + 2 + l ) . (43) Th us it follo ws fr om (43) that ρθ  r d (1) − r d (2) , r d (2) − r d (3) , r d (3) − r d (4) , . . .  = ( r (3) , r (4) , r (5) , . . . ) C, (44) where C =            1 − b 1 1 − b 2 − b 1 1 − b 3 − b 2 − b 1 1 . . . . . . . . . . . . . . .            . Therefore, we hav e C − 1 =              1 b 1 1 b 2 + b 2 1 b 1 1 b 3 + 2 b 2 b 1 + b 3 1 b 2 + b 2 1 b 1 1 b 4 + 2 b 3 b 1 + 3 b 2 b 2 1 + b 4 1 b 3 + 2 b 2 b 1 + b 3 1 b 2 + b 2 1 b 1 1 . . . . . . . . . . . . . . . . . .              , 22 and th e norm || · || ∞ of the matrix C − 1 is giv en by || C − 1 || = s u p i ≥ j ≥ 1 {| ζ i,j |} ≤ ∞ X k =1 b k = 1 , where ζ i,j is the ( i , j )th en try of th e matrix C − 1 for 0 ≤ j ≤ i . Note that the norm || C − 1 || ≤ 1 is usefu l for our follo win g iterativ e algorithm d esigned by the matrix ρθ C − 1 with || ρθ C − 1 || = ρθ < 1. Let X =  r d (1) − r d (2) , r d (2) − r d (3) , r d (3) − r d (4) , . . .  and Y = ( r (3) , r (4) , r (5) , . . . ) . Then X =  r d (1) , r d (2) , Y ⊙ d  −  r d (2) , Y ⊙ d  and it f ollo ws from (44) that Y = X  ρθ C − 1  =  r d (1) , r d (2) , Y ⊙ d   ρθ C − 1  −  r d (2) , Y ⊙ d   ρθ C − 1  . (45) It follo ws fr om (42) that r (1) = ρ − ρθ r d (1) − r (2) (1 − b 1 ) + Y ( b 1 + b 2 , b 2 + b 3 , b 3 + b 4 , . . . ) T (46) and r (2) = ρ − r ( 1) − h b 1 r (2) + ρθ r d (2) i + Y ( b 1 + b 2 − 1 , b 1 + b 2 + b 3 , b 2 + b 3 + b 4 , . . . ) T (47) No w , we use Equ ations (45) to (47) to provide an iterativ e algorithm for computing the fixed p oin t π k = r ( k ) η for k ≥ 1. T o that end, we write Y N = ( r N (3) , r N (4) , r N (5) , . . . ) (48) and R N = ( r N (1) , r N (2) , r N (3) , . . . ) = ( r N (1) , r N (2) , Y N ) . (49) 23 Let r N + 1 (1) = ρ − ρθ r d N (1) − r N (2) (1 − b 1 ) + Y N ( b 1 + b 2 , b 2 + b 3 , b 3 + b 4 , . . . ) T , (50) r N + 1 (2) = ρ − r N (1) − h b 1 r N (2) + ρθ r d N (2) i + Y N ( b 1 + b 2 − 1 , b 1 + b 2 + b 3 , b 2 + b 3 + b 4 , . . . ) T (51) and Y N + 1 =  r d N (1) , r d N (2) , Y ⊙ d N   ρθ C − 1  −  r d N (2) , Y ⊙ d N   ρθ C − 1  . (52) Based on the iterativ e relations given in (50 ) to (52 ), w e provi de an iterativ e algorithm for computing the v ecto r R = ( r (1) , r (2) , r (3) , . . . ). This giv es the fi xed p oint π = (1 , r (1) η, r (2) η, r (3) η , . . . ). An Iterative Algorithm: C omputation of the Fixed P oin t Input: λ, ( α, T ) , { b k } and d . Output: R = ( r (1) , r ( 2) , r (3) , . . . ) and π = (1 , r (1) η , r (2) η , r (3) η , . . . ). Computat ional Step s: Step one: T aking the initial v alue R 0 = 0, that is, r 0 (1) = 0 , r 0 (2) = 0 , Y 0 = 0. Step two: Comp u ting R 1 = ( r 1 (1) , r 1 (2) , Y 1 ) through r 1 (1) = ρ, ← (50) r 1 (2) = ρ, ← (51) Y 1 =  ρ d , ρ d , 0   ρθ C − 1  −  ρ d , 0   ρθ C − 1  , ← (52) Step thr e e: If R N is known, computing R N + 1 = ( r N + 1 (1) , r N + 1 (2) , Y N + 1 ) through r N + 1 (1) = ρ − ρθ r d N (1) − r N (2) (1 − b 1 ) + Y N ( b 1 + b 2 , b 2 + b 3 , b 3 + b 4 , . . . ) T , ← (50) r N + 1 (2) = ρ − r N (1) −  b 1 r N (2) + ρθ r d N (2)  + Y N ( b 1 + b 2 − 1 , b 1 + b 2 + b 3 , b 2 + b 3 + b 4 , . . . ) T , ← (51) Y N + 1 =  r d N (1) , r d N (2) , Y ⊙ d N   ρθ C − 1  −  r d N (2) , Y ⊙ d N   ρθ C − 1  . ← (52) Step four: F or a suffi cien tly small ε > 0, if there exists Step K suc h that || R K +1 − R K || < ε , then our computation is end in this step; otherw ise we go to S tep three for con tin u ous computations. 24 Step five: When our computation is o v er at Step K , computing π = (1 , r K (1) η , r K (2) η , r K (3) η , . . . ) as an app r o ximate fixed p oin t und er an error ε > 0. In what f ollo ws w e analyze t w o numerical examples by means of the ab o ve iterativ e algorithm. In the fir st example, we tak e λ = 1 , d = 2 , α = (1 / 2 , 1 / 2) , T (1) =   − 4 3 2 − 7   , T (2) =   − 5 3 2 − 7   , T (3) =   − 4 4 2 − 7   , T able 1 illustrates how the sup er-exp onent ial solution ( π 1 to π 5 ) dep end s on the matrices T (1), T (2) and T (3), resp ectiv ely . T able 1: The sup er-exp onenti al solution dep ends on the matrix T T (1) T (2) T (3) π 1 (0.2045, 0.1591) (0.1410, 0.1026) (0.312 5, 0.2500) π 2 (0.0137, 0.0107) (0.0043, 0.0031) (0.0500, 0.0400) π 3 (6.193e-05, 4.817e-05) (3.965e-06, 2.884e-06) (0.0013 , 0.0010) π 4 (1.259e-09, 9.793e-10) (3.390e-12, 2.465e-12) (8.446e-07, 6.757e-07) π 5 (5.204e-19, 4.048e-19) (2.478e-24, 1.802e-24) (3.656e-13, 2.925e-13) In the second example, we tak e λ = 1 , d = 5 , α (1) = ( 1 / 3 , 1 / 3 , 1 / 3) , α (2) = (1 / 12 , 7 / 12 , 1 / 3) , T =     − 10 2 4 3 − 7 4 0 2 − 5     , T able 2 sho ws ho w the su p er-exp onen tial s olution ( π 1 to π 4 ) dep end s on the vecto rs α (1) and α (2), resp ectiv ely . 5 Sup ermark et Mo dels with M ultiple Choices In this section, we co nsider tw o sup er m ark et models with m ultiple c hoices: The fir st one is one mobile server with multiple wa iting lines under the service discipline of joint- shortest 25 T able 2: The sup er-exp onenti al solution dep ends on the vect ors α α = ( 1 3 , 1 3 , 1 3 ) α = ( 1 12 , 7 12 , 1 3 ) π 1 (0.0741, 0.1358 , 0.2346) (0.0602, 0.1728, 0.2531) π 2 (5.619e-05, 1.030e-05, 1.779e-04 ) (7.182e-05, 2.063e-04, 3.020e-04) π 3 (1.411e-20, 2.587e-20, 4.469e-20) ( 1.739e-19, 4. 993e-19, 7.311e-19) π 4 (1.410e-98, 2.586e-98, 4.466e-98) ( 1.444e-92, 4. 148e-92, 6.074e-92) queue and serve -longest queue, and the second one is a su p ermarket mo d el with multiple classes of Po isson arriv als, eac h of which has a c hoice n um b er. Our main pur p ose is to organize the system of n onlinear equations for the fixed p oint und er multiple choic e n um b ers, and to b e able to obtain sup er-exp onentia l solution to the fixed p oints for the t w o sup ermark et mo dels. 5.1 One mobile serv er w ith m ultiple waiting lines The su p ermarket mo d el is structured as one mobile serve r with m u ltiple waiting lines, where the Poi sson arr iv als joint a waiting line with the sh ortest queue and the mobile serv er en ters a waiti ng line with the longest queue for his service w oks. Such a system is dep icted in Figure 3 f or an illustration. F or one mobile server with n wa iting lines, customers arrive at this system as a P oisson p ro cess with arr iv al rate nλ , and all customers are ser ved by one mobile serve r with service r ate nµ . Eac h arrivin g customer chooses d ≥ 1 w aiting lines indep end en tly and un iformly at ran d om fr om the n wa iting lines, and wai ts for service at a wa iting lin e which currentl y conta ins the f ew est n umber of customers. I f there is a tie, w aiting lines w ith the few est n um b er of customers will b e c hosen by the arriving cus tomer randomly . The mobile serv er c ho oses f ≥ 1 waiti ng lines ind ep endently and un iformly at random from the n w aiting lin es, and en ters a wa iting line whic h cur ren tly con tains the most num b er of cus tomers. If there is a tie, wa iting lines with the most n um b er of customers will b e c h osen b e the serve r randomly . All customers in ev er y wai ting line will b e served in the F CFS man n er. W e assume th at all rand om v ariables defin ed ab ov e are indep enden t of eac h other, and that the sys tem is op erating in the stable region ρ = λ/µ < 1. Clearly , d and f are inpu t choi ce num b er and output choice num b er in this sup er m ark et mo del, resp ectiv ely . It is clear that the sto chastic en vironmen t of this sup ermark et model is a p ositive recurrent birth-death pro cess with an irreducible infinitesimal generator Q = Q left + Q right , 26 2 1 Server Poisson arrivals Each customer Probes servers n O d d 1 d  Waiting line 1 Waiting line 2 Waiting line 3 Waiting line n 1 2 1 f  f n P Figure 3: A s up erm ark et mo del with in put and outpu t c hoices where Q left =            0 µ − µ µ − µ µ − µ . . . . . .            and Q right =         − λ λ − λ λ − λ λ . . . . . .         . Similar d eriv ation to those giv en in Section 3 or 4, w e obtain that the fi x ed p oin t satisfies the system of n onlinear equations π 0 = 1 and π ⊙ f Q left + π ⊙ d Q right = 0 . (53) Let Q =   Q 0 , 0 U V Q ( L )   , where Q 0 , 0 = − λ, U = ( λ, 0 , 0 , . . . ) , V = ( µ, 0 , 0 , . . . ) T , 27 Q ( L ) arriv al =         − λ λ − λ λ − λ λ . . . . . .         and Q ( L ) service =         − µ µ − µ µ − µ . . . . . .         . It follo ws fr om (53) that π 0 = 1 , (54) − λπ d 0 + µπ f 1 = 0 (55) and π d 0 U + π ⊙ d L Q ( L ) arriv al + π ⊙ f L Q ( L ) service = 0 . (56) It follo ws fr om (54) and (55) th at π 1 = ρ 1 f . Note that h − Q ( L ) service i − 1 =         1 µ 1 µ 1 µ 1 µ 1 µ 1 µ . . . . . . . . . . . .         , it f ollo ws from (56 ) that for k ≥ 2 π f k = π d k − 1 ρ. This leads to π k = ρ k − 1 P i =0 d i f k − 1 − i f k = ρ 1 f k − 1 P i =0  d f  i . (57) Sp ecifically , when d 6 = f , we hav e π k = ρ ( d f ) k − 1 d − f . 28 Remark 3 Equation (57) indic ates differ ent i nfluenc e of the input and output choic e numb ers d and f on the fixe d p oint π . If d > f , then the fixe d p oint π de c r e ases doubly exp onential ly; and if d = f , then π k = ρ k f which i s ge ometric. However, it is very inter esting for the c ase with d < f . In this c ase, lim k →∞ π k = ρ 1 f − d , which il lustr ates that the fr action of waiting lines with infinite customers has a p ositive lower b ound ρ 1 f − d > 0 . This shows that if ρ < 1 and d < f , this sup ermarket mo del is tr ansient. 5.2 A sup ermark et mo del with m ultiple input c hoices No w , w e analyze a s up erm ark et mo d el with m ultiple in put c hoices. There are m types of differen t cu s tomers w ho arriv e at a queueing system of n > 1 servers for receiving their required service. Arriv als of customers of i th t yp e are a P oisson pro cess with arriv al rate nλ i for λ i > 0, and the service times at eac h server are exp onentia l with service rate µ > 0. Note that differen t t yp es of customers hav e the same service time. Eac h arriving customer of i th t yp e c ho oses d i ≥ 1 s erv ers indep enden tly and uniform ly at rand om from the n servers, and waits for service at the serv er w hic h currentl y conta ins the fewest n um b er of cus tomers. If there is a tie, serve rs with the fewe st num b er of customers will b e c hosen r andomly . All customers in eve ry serve r will b e serv ed in the F CFS manner. W e assume that all ran d om v ariables d efined ab o v e are indep end en t of eac h other, and that the system is op erating in the stable region ρ = P m i =1 ρ i < 1, where ρ i = λ i /µ . Clearly , d 1 , d 2 , . . . , d m are multiple inp ut c hoice num b ers in this sup ermark et mo del. Let Q right ( i ) =         − λ i λ i − λ i λ i − λ i λ i . . . . . .         and Q left =         0 µ − µ µ − µ . . . . . .         . Ob viously , the sto c hastic environmen t of this su p ermarket mo del is a p ositiv e r ecurrent birth-death pr o cess w ith an irred ucible infinitesimal generator Q = Q left + P m i =1 Q right ( i ). 29 Similar d eriv ation to those giv en in Section 3 or 4, w e obtain that the fi x ed p oin t satisfies the system of n onlinear equations π 0 = 1 (58) and π Q left + m X i =1 π ⊙ d i Q right ( i ) = 0 . (59) Let π = ( π 0 , π L ) , U i = ( λ i , 0 , 0 , . . . ) , Q ( L ) arriv al ( i ) =         − λ i λ i − λ i λ i − λ i λ i . . . . . .         and Q ( L ) service =         − µ µ − µ µ − µ . . . . . .         . Therefore, the s y s tem of nonlinear equations (58) and (59) is w ritten as π 0 = 1 , (60) − m X i =1 λ i π d i 0 + µ π 1 = 0 , (61 ) m X i =1 π d i 0 ( λ i , 0 , 0 , . . . ) + m X i =1 π ⊙ d i L Q ( L ) arriv al ( i ) + π L Q ( L ) service = 0 . (62) It follo ws fr om (60) and (61) th at π 1 = m X i =1 ρ i = ρ, and f rom (62) that π L = m X i =1 π d i 0 ( λ i , 0 , 0 , . . . ) h − Q ( L ) service i − 1 + m X i =1 π ⊙ d i L Q ( L ) arriv al ( i ) h − Q ( L ) service i − 1 . 30 This leads to that for k ≥ 2 π k = m X i =1 π d i k − 1 ρ i . Let δ 1 = ρ an d δ k = P m i =1 δ d i k − 1 ρ i for k ≥ 2. T hen the fixed p oint has a sup er- exp onenti al s olution π 0 = 1 and f or k ≥ 1 π k = δ k . Ac kno wledgemen ts Q.L. Li was supp orted by the National Science F ound ation of Ch in a under gran t No. 10871 114. References [1] Y. Azar, A.Z. Bro d er, A.R. Karlin and E. Upfal (1999). Balanced allocations. SIAM Journal on Computing 29 , 180–200. [2] M. Bramson, Y. Lu and B. P r abhak ar (2010). Randomized load balancing with general service time distrib u tions. In Pr o c e e dings of the ACM SIGMETRICS international c onfer enc e on Me asur ement and mo deling of c omputer systems , p ages 275–286. [3] M. Dahlin (1999) . Interpreting stale load inf ormation. IEEE T r ansactions on Par al lel and Distribute d Systems 11 , 1033–10 47. [4] C. Graham (2000). Kinetic limits for large comm unication net w orks. In Mo del ling in Applie d Sci-enc es , N. Bellomo and M. Pulvirent i (eds.), Birkh¨ auser, pages 317–370 . [5] C. Graham (2000). C haoticit y on path space for a queueing net w ork w ith selection of the shortest queue among sev eral. J ournal of Applie d P r ob ab ability 37 , 198–201. [6] C. Graham (2004). F un ctional cent ral limit theorems for a large n et work in whic h customers join the sh ortest of sev eral qu eues. Pr ob ability The ory R elate d Fie lds 131 , 97–12 0. 31 [7] M. Harc hol-Balte r and A.B. Downey (1997). Exp loiting pro cess lifetime d istributions for dynamic load balancing. ACM T r ansactions on Computer Systems 15 , 253–2 85. [8] T.G. Kur tz (1981). Appr oximation of P opulation Pr o c esses . SIAM. [9] Q.L. Li (2010 ). Constructive Computation in Sto chastic Mo dels with A pplic ations: The R G-F actorizations . Springer an d T singh ua Pr ess. [10] Q.L. Li (2010). Doubly exp onentia l solution for rand omized load balancing with gen- eral ser v ice times. Submited for p ublication. [11] Q.L. Li and John C.S. Lui (2010). Doubly exp on ential solution for randomized load balancing mo dels with Mark o vian arriv al pro cesses and PH service times. Sub mited for pub licatio n. [12] Q.L. Li, John C.S . Lui and Y. W ang (2010). A matrix-analytic solution for r andom- ized load balancing mo d els with P H service times. In Pr o c e e ding of P ERF ORM 2010 Workshop , Lecture Notes of C omputer Science, Page s 1-20. [13] Q.L. Li, John C.S. Lui and Y. W ang (2010). Su p er-exp onen tial solution for a retrial sup er m ark et Mmo del. Sub mited for publication. [14] M. Luczak and C. McDiarmid (2006). On the maximum qu eue length in the sup er- mark et mo del. The Anna ls of Pr ob ability 34 , 493–52 7. [15] M. Luczak an d C. McDiarmid (2007). Asymptotic distribu tions an d c h aos for the sup er m ark et mo del. Ele c tr onic Journal of Pr ob ability 12 , 75–99. [16] M.J. Lu czak and J.R. Norris (2005). Str on g approximati on for the sup ermark et mo del. The Annals of A pplie d Pr ob ability 15 , 2038–20 61. [17] J.B. Martin (2001). Poin t pr o cesses in fast Jac kson net w orks. The Annals of Applie d Pr ob ability 11 , 650-66 3. [18] J.B. Martin and Y.M Suhov (1999). F ast J ackson net w orks. The Annals of Applie d Pr ob ability 9 , 854–870 . [19] M.D. Mitzenmac h er (1996). Th e p ow er of t w o choic es in rand omized load balancing. PhD thesis, Universit y of California at Berk eley , Department of Computer Science, Berk eley , CA. 32 [20] M. Mitzenmac her (1999). On the analysis of randomized load balancing schemes. The ory of Computing Systems 32 , 361–386 . [21] M. Mitzenmac her (2000). Ho w useful is old information? IEEE T r ansactions on Par al lel and Distribute d Systems 11 , 6–20. [22] M. Mitzenmac her, A. Ric ha, an d R. Sitaraman (20 01). The p o w er of t w o random c hoices: a surv ey of tec hniques and results. In Handb o ok of R andomize d Com puting: V olume 1 , P . Pardalos, S. Ra jasek aran and J. Rolim (eds), pages 255-31 2. [23] R. Mirchandaney , D. T o w s ley , and J.A. Stank o v ic (1989). An alysis of the effects of dela ys on load sh aring. IEE E T r ansactions on Computers 38 , 1513–1525. [24] M.F. Neuts (1989). Structur e d sto chastic matric es of M /G/ 1 typ e and their applic a- tions. Marcel Deck er Inc.: New Y ork. [25] Y.M. Su h o v and N.D. Vv edensk a y a (2002) . F ast Jac kson Net works with Dynamic Routing. P r oblems of Information T r ansmission 38 , 136 { 15 3. [26] B. V¨ oc kin g (1999). How asymmetry helps load balancing. In Pr o c e e dings of the F or- tieth Annual Symp osium on F oundations of Computer Sc i enc e , pages 131–140. [27] N.D. Vv edensk a y a, R.L. Dobrushin and F.I. K arp elevic h (1996). Q ueueing system with selection of the shortest of t wo queues: An asymptoti c approac h . Pr oblems of Information T r ansmissions 32 , 20–34. [28] N.D. Vvedensk a ya and Y.M. Su ho v (1997 ). Dobrushin’s mean-field app ro ximation for a qu eue with dyn amic r outing. Markov Pr o c esses and R elate d Fields 3 , 493–52 6. [29] R. W eb er (1978). On the optimal assignment of customers to parallel servers. Journal of A pplie d Pr ob abiblities 15 , 406–413 . [30] W. Winston (1977) . Optimalit y of the sh ortest line discipline. Journal of Applie d Pr ob abilities 14 , 181–189 . 33