Blanc, J.P.C.

English

Research Papers in Economics

Journal  of  Computational  and  Applied  Mathematics  20  (1987)  119-125 
North-Holland 
119 
On  a numerical  method  for  calculating  state 
probabilities  for  queueing  systems 
with  more  than  one  waiting  line 
J.P.C.  BLANC  * 
Department  of  Mathematics  and  Informatics,  Delft  University  of  Technology,  2600  AJ  De&,  The  Netherlands 
Received  26  May  1986 
Abstract:  Keane,  Hooghiemstra  and  Van  de  Ree  have  proposed  a  new  numerical  method  for  calculating  state 
probabilities  for  queueing  systems  with  more  than  one  waning  line  in  parallel.  The  method  is based  on  power  series 
expansions  of  state  probabilities  as functions  of  the  traffic  intensity  of  a system.  The  coefficients  of  these power  series 
can  be recursively  calculated.  The  coefficients  of  the power series  expansions  of  moments  of  queue  length  distributions 
can  be  derived from  those  of  the state  probabilities  in  a straightforward  manner.  The  above  method  is discussed  for  a 
rather  general  class  of  exponential  queueing  systems.  The  asymptotic  behaviour  of  moments  in heavy  traffic  is used to 
obtain  extrapolations  of  the  coefficients  of  their  power  series  expansions  at  the  origin.  The  calculation  of  moments  is 
strongly  improved  by  means  of  these  extrapolations. 
Keywords:  Queues  in  parallel,  moments,  power  series  expansions,  traffic  intensity. 
1. Introduction 
Queueing  systems  with  more  than  one  waiting  line  are  in  general  hard  to  analyse  if  the joint 
queue  length  distribution  does  not  possess  a  product  form.  During  the  last  decade  a  rather 
general  analytical  method  has  been  developed  for  the  study  of  systems  with  two  waiting  lines. 
This  method  is  based  on  application  of  the  theory  of  Riemann-Hilbert  type  boundary  value 
problems  to  solve  functional  equations  for  bivariate  generating  functions  of joint  queue  length 
distributions,  see  a  [S]  and  [2].  At  present  it  does  not  seem  likely  that  this  method  can  be 
generalised  as useful  tool  for  the  analysis  of  systems  with more  than  two waiting  lines.  Moreover, 
in  many  instances  non-trivial  algorithms  (e.g.  the  numerical  solution  of  integral  equations)  are 
required  to  obtain  numerical  data  from  the  formulas  which  result  by  application  of  this  method. 
Numerical  methods  for  calculating  state  probabilities  without  the  use  of  generating  functions 
have  been  proposed  by  several  authors.  They  include  Conolly  [4],  who  uses  truncation  of  the 
state  space  (i.e.  of  the  waiting  rooms)  and  inversion  of  the  matrix  of  transition  rates,  and 
Gertsbakh  [6], who  uses truncation  of  the  state  space  and  numerical  evaluation  of  the  matrix-ge- 
*  Present  address:  Catholic  University  of  Brabant,  Faculty  of  Economics,  P.O.  Box  90153,  5000  LE  Tilburg,  The 
Netherlands. 
0377-0427/87/$3.50  0  1987,  Elsevier  Science  Publishers  B.V.  (North-Holland) 120  J. P. C. Blanc  /  Numerical  method  for  queueing  systems 
ometric  solution  of  the  truncated  system.  Keane  et al. [7] propose  another  method  for  calculating 
state  probabilities  directly  from  the  equilibrium  equations.  They  consider  power  series  expan- 
sions  of  the  state  probabilities  at  the  origin  as functions  of  the  traffic  intensity.  For  a wide class 
of  queueing  systems  the  coefficients  of  these  power  series  expansions  can  be  calculated  recur- 
sively. Keane  et al. [7] found  that  for  several  models  a (bilinear)  conformal  transformation  has  to 
be  used  to  ensure  convergence  of  the  power  series  over  the  whole  range  of  traffic  intensities  for 
which  a system  is  stable.  The  coefficients  of  the  power  series  expansions  of  the  moments  of  the 
queue  length  distributions  follow  directly  from  those  of  the  state  probabilities.  The  number  of 
coefficients  required  to  reach  a  certain  accuracy  can  be  strongly  reduced  by  extrapolation  of 
these  coefficients  on  the  basis  of  the  singular  behaviour  of  the  moments  near  the  end  of  the 
region  of  stability  of  the  queueing  system. 
The  idea  of  the  method  of  Keane  et  al.  [7]  will  be  outlined  in  Section  2  for  a  class  of 
exponential  queueing  systems;  Section  3 is devoted  to  the  procedure  for  calculating  moments  of 
queue  length  distributions.  In  Section  4  the  procedure  will  be  illustrated  by  means  of  some 
examples,  and  experiences  with  the  method  will be  discussed. 
2.  Calculating  state  probabilities 
Consider  the  class  of  stable  exponential  queueing  systems  with  s  (s  >  1)  waiting  lines  which 
can  be  described  in the  following  way.  Let  Z =  (n,,  . . . , n,)  be  a vector  with  non-negative  entries; 
let  paj(Z)  be  the  arrival  rate  to  queue  j  and  dj(Z)  be  the  departure  rate  from  queue  j, 
j=l  >...,  s,  when  the  system  is in  state  Z. Let  N,  be  the  number  of  customers  in  queue  j  in  the 
stationary  situation,  j  =  1,.  . . , s,  and  define  the  stationary  state  probabilities 
p(p;  E)=Pr{N,=n;,  j=l,...,  s}.  (2.1) 
Let  I(E)  stand  for  the  indicator  function  of  the  event  E,  and  Zj =  ( ejl,  _  . . , ej,)  for  the  vector 
with  e,j  =  1 and  eji =  0 if  i  #j,  j  =  1,.  . . , s.  The  equations  for  the  state  probabilities  (2.1)  read: 
I  p  2 LZj(ii) +  i  dj(n)l(njzO)  P(Pi ‘1 
/=I  j=l  I 
=  C  dj(2+ZJ)p(p;  ii+Zj)  +P  C  LZj(Fi-ej)p(p;  n-zj)l(nj>O)'  (2-4 
j=l  j=l 
The  following  assumptions  are  made: 
(Al)  the  system  is stable  for  0 G p <  1; 
(A2)  the  state  probabilities  p(p;  Fz)  possess  analytic  continuations  as  functions  of  p  into  a 
domain  which  includes  the  disk  1  p  -  i  1 <  :  ; 
(A3)  not  all servers  are idle  when  customers  are  present  in  the  system,  i.e. 
idj(E)l(nj>O)>O  foralln,  Z#O=(O,...,O). 
j=l J. P. C. Blanc  /  Numerical  method for  queueing  systems  121 
Introducing  the  notations,  for  m = 0,  1, 2,.  . . ,  1  E 1 =  C:=lnj, 
Iiil=m  J=l 
A(Pi m>  =  ,_E  P(Pi ‘> C  aj(n>,  (2.3) 
n  m  j=l 
we  obtain  by  summation  of  the  equations  in  (2.2)  over  all  states  Z  with  1  ii  1 =  m,  for 
m=l,2  >..*> 
P&?  m> +  NP;  m> = D(p;  m+l)+pA(p,  m-l),  pA(p;O)=D(p;  1). 
This  implies,  by  induction,  that 
pA(p;  m)=D(p;  m+l),  m=0,1,2  )...  .  (2-4 
From  (2.3)  and  assumption  (A3)  it  follows,  for  m =  1,  2,.  . . , that  if  D(p;  m)  =  0(  pk),  p  JO,  for 
some  k,  then  also  p(p;  Z) =  O(p“),  p JO,  for  all  7i with  I ii  I =  m,  and  hence  A(p;  m)  =  O(pk), 
p 1 O.-Therefore,  relations  (2.4)  imply  that  D(p;  m)  =  O(p”),  p JO,  m =  1,  2,.  . .  (note  that 
p( p; 0) =  O(l),  p  JO).  Hence,  if  state  probabilities  p( p;  n)  satisfy  (2.2)  and  if  assumption  (A3) 
holds,  then 
PC&  n)  =  O(C),  P JO.  (2.5) 
Motivated  by  this  property  Keane  et  al.  [7] propose  the  following  procedure.  Write 
p(p;  Z)  =p”’  5  pku(k;  ii),  (24 
k=O 
and  substitute  these  power  series  expansions  into  (2.2).  Equating  coefficients  of  corresponding 
powers  of  p leads  to:  for  k =  0,  1,.  . . , 
c d,(i)l(nj>  O)u(k;  ii) 
j=l 
=  -J$,o,(i)u(k-l:  Z)I(k>O)+  &,(E+e,)u(k-1;  i+S-,)l(k>O) 
j=l 
+  c  a,(+e,)l(n,>O)u(k;  E-Zj). 
j=l 
(2.7) 
To  obtain  a  relation  for  u( k;  ?I), k =  0,  1,.  . . ,  the  law  of  total  probability  is  used:  substituting 
(2.6)  and  equating  corresponding  powers  of  p  gives 
u(0;  0)  =  1;  u(k;  fi)=  -c...cu(k-  17~1; Z),  k=l,2  ,...  .  (2.8) 
O<Iiil<k 
It  is  readily  verified  that  all  coefficients  u( k;  Z)  can  be  recursively  calculated  from  (2.8)  and 
(2.7).  The  state  probabilities  p(p;  E)  can  then  be  approximated  with  any  degree  of  accuracy  for 
p  smaller  than  the  radius  of  convergence  of  their  power  series  (2.6).  When  these  radii  of 
convergence  are  smaller  than  one,  Keane  et  al.  [7]  propose  to  apply  a  bilinear  conformal 
transformation  of  the  real  interval  [O,l]  onto  itself: 
(j=  l+G  8 
l+GpP’  ‘=  l+G-G8’  G>,  0.  (2.9) 122  J. P. C. Blanc  /  Numerical  method for  queueing  systems 
Studying  the  state  probabilities  as functions  of  6  it  can  be  shown  in  a similar  way as above  that, 
cf.  (2.6) 
(2.10) 
where,  cf.  (24,  (2.9), 
b(0;  Tj) =  1;  W;  Q=  -C**.Cb(k-  InI;  ?I),  k=l,2  )...)  (2.11) 
O<lril<k 
and  for  E #  0, 
(1 +  G)  c  d,(n)Z(nj>O)b(k  n> 
j=l 
+  Cd,(n+e,)[(l+G)b(k-1;  Fi+Zj)Z(k>O)-Gb(k-2;  E+e,)Z(k>l)] 
j=l 
-  ‘&,(Z)-G&)Z(nj>O)  b(k-1;  E)Z(k>O),  k=O,l,...  . 
j=l  j=l  1 
(2.12) 
By  assumption  (A2)  the  constant  G can  be  chosen  in  such  a way that  the  radii  of convergence  of 
the  power  series  (2.10)  are  at  least  one. 
Remark.  Note  that  the  coefficient  of  u(k;  E)  in  (2.7)  and  that  of  b( k;  E)  in  (2.12)  does  not 
vanish  when  E #  8 by  assumption  (A3). 
3.  Calculating  moments 
The  coefficients  of  the  power  series  expansions  of  the  moments  of  queue  length  distributions 
can  be  directly  calculated  from  those  of  the  state  probabilities.  For  i,  j  =  1,.  . . , s,  let,  cf.  (2.9), 
(2.10), 
E{  T.}  =  f  Z+K  E{N,N,}  =  E  hjj(k)8k.  (3.1) 
k=l  k=l 
It  is easily  verified  that  for  i,  j  =  1,.  . . , s,  k  =  1,  2,.  . . , 
f,(k)  =  ~<;;;<;Q(k-  IEli  3,  hjj(k)=  ~..+,njb(k-  IFzl;  Z).  (3.2) 
.  .  O<liil<k 
As  can  be  expected  the  power  series  of  the  moments  converge  very  slowly  for  8  close  to  one 
(heavy  traffic,  assumption  (Al)),  because  the  nth  moments  will have  a pole  of  (at  most)  order  n 
at  8 =  1.  The  latter  implies  for  instance  that  f,(k)  tends  to  some  constant  and  h,J  k)  tends  to J. P. C. Blanc  /  Numerical  method for  queueing  systems  123 
some  linearly  increasing  sequence  as  k  +  00  (i,  j  =  1,.  . . , s).  Therefore,  we  propose  to  use 
&(W  +f,(WS, 
k=l 
(3.3) 
as  approximations  for  E{  Nj}  and  E{  NiNj}  respectively,  when  the  coefficients  of  the  power 
series  expansions  have  been  calculated  up  to  the  Mth  power  of  8.  These  extrapolations  of  the 
power  series  of  moments  strongly  improve  the  rates  of  convergence  to  the  exact  values  of  the 
moments  when  0  is close  to  one.  E.g.  in  a  typical  case  with  p =  0.9  the  relative  error  between  the 
approximate  and  the  exact  value  of  the  average  total  number  of  customers  in  a  system  is 0.03  for 
M  =  60  without  the  use  of  the  extrapolation,  and  0.001  for  A4 =  24  with  the  aid  of  the 
extrapolation. 
4.  Examples 
In  this  section  two  queueing  models  will  be  discussed  for  which  the  equations  to  be  satisfied 
by  the  state  probabilities  fall  within  the  general  framework  of  (2.2).  The  interarrival  times  and 
the  service  times  have  negative  exponential  distributions  in  both  models. 
Example  1.  The  shortest  queue  model:  customers  choose  with  equal  probabilities  one  of  the 
shortest  of  s queues  upon  arrival;  customers  in  queue  j  are  served  by  server  j  with  service  rate  5, 
j=  l,...,  S; Cg=irj  =  1.  For  all  states  Z  in  this  model,  for  j  =  1,.  . . , s, 
dj(E)  =  rj,  ai  =  0  if  3i  n,  <  nj, 
=  l/v,  if  Vi  n, 2  nj,  v=  #{i;  ni=nj}.  (44 
The  special  case  s =  2,  ri  =  r,  =  :,  has  been  studied  by  several  authors;  Conolly  [4]  and 
Gertsbakh  [6]  provided  numerical  data.  In  [l]  moments  of  the  waiting  time  distribution  have 
been  calculated  for  more  general  cases  (S =  2,  3)  based  on  the  present  method.  In  Table  1 values 
of  the  average  total  number  of  customers  in  the  system  and  the  correlation  between  the  number 
of  customers  in  two  queues  are  given  for  the  symmetrical  models  with  two  and  three  servers  for 
several  values  of  the  traffic  intensity  p. 
Table  1 
The  shortest queue model 
P 
0.10 
0.30 
0.50 
0.70 
0.90 
0.98 
s=2,  rI=r2=+ 
E{N,+N,I 
0.2035 
0.6864 
1.426 
2.951 
9.855 
49.97 
CoMN,,  Nz) 
0.09982 
0.3032 
0.5257 
0.7659 
0.9647 
0.9984 
s =  3, rl =  r,  =  r,  = f 
E{N,+N,+N,) 
0.3010 
0.9579 
1.867 
3.589 
10.75 
51.00 
Corr(N,,  N2) 
0.09098 
0.1984 
0.3925 
0.6560 
0.9404 
0.9973 124  J.P.  C.  Blanc  /  Numerical  method  for  queueing  systems 
Table  2 
The  pre-emptive longer  queue model  (s =  2, pL1  = 4~~) 
P  A,  =  4h,  A,  =  A,  4h,  =  A, 
E{Nll  E{N2)  E{NlI  E{NzI  EINlI  EiN21 
0.10  0.05588  0.05547  0.02283  0.08874  0.006823  0.1045 
0.30  0.2211  0.2126  0.1061  0.3378  0.03656  0.4005 
0.50  0.5264  0.4934  0.3283  0.7679  0.1455  0.9195 
0.70  1.230  1.151  1.037  1.724  0.6286  2.073 
0.90  4.621  4.470  5.343  6.314  4.743  7.417 
0.98  24.65  24.46  32.44  33.54  34.67  38.17 
Example  2.  The  pre-emptive  longest  queue  model:  s  independent  customer  streams  form  s 
queues,  one  server  provides  service  to  a  customer  which  is  selected  with  equal  probabilities  from 
one  of  the  longest  queues;  at  each  arrival  instant  the  server  reselects  the  customer  to  be  served 
according  to  this  rule.  In  this  model  we  have  for  each  state  7i, for  j  =  1,.  . . , s, 
ai  =  Aj,  dj(Fi)  =  0,  if  3i  n,  > nj, 
=pj/v,  ifVini<nj,  v=  #{i;  ni=nj}.  (4.2) 
The  arrival  rates  are  normalized  by  the  relation 
i  xj/Pj=l,  (4.3) 
j=l 
in  order  to  satisfy  the  assumption  (Al).  In  the  Tables  2 and  3 values  of  mean  queue  lengths  and 
correlation  between  queue  lengths  are  displayed  for  the  case  of  two  queues,  a  ratio  between 
service  rates  of  4 : 1,  and  ratios  between  arrival  rates  of  4 : 1,  1 : 1,  1 : 4,  respectively.  Models  with 
s queues  and  one  server,  without  a pre-emptive  service  discipline,  such  as  the  longer  queue  model 
studied  by  Cohen  [3], where  a customer  is selected  from  the  longer  of  two  queues  and  then  served 
without  interruption,  and  the  alternating  service  model  studied  by  Cohen  and  Boxma  [2], are  not 
contained  in  the  framework  of  equations  (2.2).  But  by  extending  the  state  space  with  a  variable 
indicating  the  queue  from  which  a  customer  is  being  served  these  models  are  readily  adapted  to 
the  present  procedure  for  calculating  state  probabilities  and  moments. 
In  general,  the  choice  of  G  in  (2.9),  i.e.  the  choice  of  the  conformal  mapping,  is  difficult, 
because  the  state  probabilities  p(  p;  Z),  as  functions  of  p,  may  possess  different  singularities 
Table  3 
The  pre-emptive  longer  queue model  (s =  2, pL1  =  4p2) 
P  h, = 4x,  &=A,  4x,  =  A, 
CON%  N2)  E{IN,-N2ll  CoNN,,  N2)  E{IN,-N2ll  CoMN,,  N2)  E{IN,-N2Il 
0.10  0.06604  0.09978  0.06451  0.1035  0.04175  0.1084 
0.30  0.2629  0.2964  0.3107  0.3310  0.2252  0.3904 
0.50  0.5062  0.4901  0.6055  0.5838  0.4765  0.8089 
0.70  0.7612  0.6858  0.8402  0.8603  0.7131  1.482 
0.90  0.9638  0.8909  0.9793  1.159  0.9307  2.708 
0.98  0.9984  0.9777  0.9991  1.284  0.9957  3.525 J. P.C.  Blanc  /  Numerical  method  for  queueing  systems  125 
inside  the  unit  disk  for  different  states  Z.  G  should  be  at  least  so  large  that  all  singularities  of  all 
state  probabilities  and  moments  (which  are  unknown  in  general)  are  mapped  outside  the  unit 
disk  (when  all  values  of  the  traffic  intensity  p,  0 < p  <  1,  are  considered).  However,  choosing  G 
much  larger  often  leads  to  slower  convergence  of  the  power  series,  because  other  singularities 
(e.g.  on  the  interval  (1,  cc))  are  mapped  closer  to  the  unit  disk  as  G  increases.  In  some  cases  the 
power  series  expansions  of  moments  did  converge  for  a specific  value  of  G,  while  those  of  several 
state  probabilities  did  not.  A  value  which  often  gives  good  results  is  G =  1. 
It  seems  that  limitations  on  application  of  the  method  are  caused  more  by  the  amount  of 
memory  space  required  to  store  the  coefficients  of  the  power  series  than  by  the  amount  of 
computation  time.  When  power  series  are  needed  up  to  the  Mth  power  of  8 coefficients  b( k;  Z) 
have  to  be  calculated  for  all  k  and  E  with  k  +  1  Z 1 G M:  To  avoid  reservation  of  superfluous 
memory  space  points  (k,  Fi)  could  be  mapped  onto  the  integers,  e.g.  in  the  following  way 
/ 
\ 
(k,  n,  ,...,  n,)+  i  k+‘+  i$‘i  . 
j=O 
j+l  I 
Then  (,p+i)  memory  positions  are  required  instead  of  (M  +  I)S+l. 
Remark.  Assumption  (Al)  requires  that  the  stability  condition  of  a  system  is  known  in  order  to 
normalize  the  arrival  rates  a_J Fz) for  all  states  Z  and  j  =  1,.  . . , s.  If  the  stability  condition  is not 
known,  it  can  be  obtained  from  the  coefficients  of  the  power  series  expansion  of  the  average 
number  of  customers  in  the  system  which  will  tend  to  a  geometrical  sequence.  In  a  similar  way 
inspection  of  the  coefficients  of  the  power  series  expansions  of  the  state  probabilities  will  give 
information  about  the  radii  of  convergence  of  these  series.  This  information  can  be  used  to 
choose  an  appropriate  conformal  mapping. 
References 
[l]  J.P.C.  Blanc,  A  note  on  waiting  times  in  systems  with queues  in parallel,  J.  Appl.  Probab.,  to  appear. 
[2]  J.W.  Cohen,  O.J.  Boxma,  Boundary  Value  Problems  in  Queueing  system  Analysis  (North-Holland,  Amsterdam, 
1983). 
[3]  J.W.  Cohen,  A  two  queue,  one  server  model  with  priority  for  the  longer  queue,  University  of  Utrecht,  Dept.  of 
Math.,  preprint  379,  1985. 
[4]  B.W.  Conolly,  The  autostrada  queueing  problem,  J.  AppI.  Probab.  21  (1984)  394-403. 
[5]  G.  Fayolle  and  R.  Iasnogorodski,  Two  coupled  processors:  the  reduction  to  a  Riemann-Hilbert  problem,  Z. 
Wahrscheinlichkeitsth.  47  (1979)  325-351. 
[6]  I.  Gertsbakh,  The  shorter  queue  problem:  a  numerical  study  using  the  matrix-geometric  solution,  Eur.  J.  Oper. 
Rex  15  (1984),  374-381. 
[7]  M.S.  Keane,  G.  Hooghiemstra  and  S.  van  de Ree,  Private  communication. 

A note on waiting times in systems with queues in parallel,

A two queue, one server model with priority for the longer queue,

Boundary Value Problems in Queueing system Analysis (North-Holland,

The autostrada queueing problem,

The shorter queue problem: a numerical study using the matrix-geometric solution,

Two coupled processors: the reduction to a Riemann-Hilbert problem,

On a numerical method for calculating state probabilities for queueuing systems with more than one waiting line.

On a numerical method for calculating state probabilities for queueuing systems with more than one waiting line.

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics