Buzen's algorithm

Consider a closed queueing network with M service facilities and N circulating customers. Write n_i(t) for the number of customers present at the ith facility at time t, such that ${\displaystyle \scriptstyle \sum _{i=1}^{M}n_{i}=N}$. We assume that the service time for a customer at the ith facility is given by an exponentially distributed random variable with parameter μ_i and that after completing service at the ith facility a customer will proceed to the jth facility with probability p_ij.[2]

It follows from the Gordon–Newell theorem that the equilibrium distribution of this model is

${\displaystyle \mathbb {P} (n_{1},n_{2},\cdots ,n_{M})={\frac {1}{{\text{G}}(N)}}\prod _{i=1}^{M}\left(X_{i}\right)^{n_{i}}}$

where the X_i are found by solving

${\displaystyle \mu _{j}X_{j}=\sum _{i=1}^{M}\mu _{i}X_{i}p_{ij}\quad {\text{ for }}j=1,\ldots ,M.}$

and G(N) is a normalizing constant chosen that the above probabilities sum to 1.[1]

Buzen's algorithm is an efficient method to compute G(N).[1]

Write g(N,M) for the normalising constant of a closed queueing network with N circulating customers and M service stations. The algorithm starts by noting solving the above relations for the X_i and then setting starting conditions[1]

${\displaystyle g(0,m)=1{\text{ for }}m=1,2,\cdots ,M}$

${\displaystyle g(n,1)=(X_{1})^{n}{\text{ for }}n=0,1,\cdots ,N.}$

The recurrence relation[1]

${\displaystyle g(n,m)=g(n,m-1)+X_{m}g(n-1,m).}$

is used to compute a grid of values. The sought for value G(N) = g(N,M).[1]

The coefficients g(n,m), computed using Buzen's algorithm, can also be used to compute marginal distributions and expected number of customers at each node.

${\displaystyle \mathbb {P} (n_{i}=k)={\frac {X_{i}^{k}}{G(N)}}[G(N-k)-X_{i}G(N-k-1)]\quad {\text{ for }}k=0,1,\ldots ,N-1,}$

${\displaystyle \mathbb {P} (n_{i}=N)={\frac {X_{i}^{N}}{G(N)}}[G(0)].}$

the expected number of customers at facility i by

${\displaystyle \mathbb {E} (n_{i})=\sum _{k=1}^{N}X_{i}^{k}{\frac {G(N-k)}{G(N)}}.}$

It will be assumed that the X_m have been computed by solving the relevant equations and are available as an input to our routine. Although g is in principle a two dimensional matrix, it can be computed in a column by column fashion starting from the leftmost column. The routine uses a single column vector C to represent the current column of g.

C[0] := 1
for n := 1 step 1 until N do
   C[n] := 0;

for m := 1 step 1 until M do
for n := 1 step 1 until N do
   C[n] := C[n] + X[m]*C[n-1];

At completion, C contains the desired values G(0), G(1) to G(N). [1]

• Jain: The Convolution Algorithm (class handout)
• Menasce: Convolution Approach to Queueing Algorithms (presentation)