This method was first proposed by Jeffrey P. Buzen in his 1971 PhD dissertation[1] and subsequently published in a refereed journal in 1973.
[3] Performing a naïve computation of the normalizing constant requires enumeration of all states.
For a closed network with N circulating customers and M service facilities, G(N) is the sum of
individual terms, with each term consisting of M factors raised to powers whose sum is N. Buzen's algorithm computes G(N) using only NM multiplications and NM additions.
This dramatic improvement opened the door to applying the Gordon-Newell theorem to models of real world computer systems as well as flexible manufacturing systems and other cases where bottlenecks and queues can form within networks of inter-connected service facilities.
[4] The values of G(1), G(2) ... G(N -1), which can be used to calculate other important quantities of interest, are computed as by-products of the algorithm.
Consider a closed queueing network with M service facilities and N circulating customers.
be the steady state probability that the number of customers at service facility i is equal to ni for i = 1, 2, ... , M .
Buzen's algorithm represents the first efficient procedure for computing G(N).
[2][4] The individual terms that must be added together to compute G(N) all have the following form:
Note that this set of terms can be partitioned into two groups.
, a surprising result emerges: the modified terms in the first group are identical to the terms used to compute the normalizing constant for the same network with one customer removed.
Thus, the sum of the terms in the first group can be written as “XM times G(N -1)”.
As a result, service facility M effectively disappears from all terms in this group (since it reduces in every case to a factor of 1).
This leaves the total number of customers at the remaining M -1 service facilities equal to N. The second group includes all possible arrangements of these N customers.
To express this concept precisely, assume that X1, X2, … XM have been obtained for a given network with M service facilities.
For any n ≤ N and m ≤ M, define g(n,m) as the normalizing constant for a network with n customers, m service facilities (1,2, … m), and values of X1, X2, … Xm that match the first m members of the original sequence X1, X2, … XM .
Given this definition, the sum of the terms in the second group can now be written as g(N, M -1).
In addition, the normalizing constant G(N) in the Gordon-Newell theorem can now be re-written as g(N,M).
Since G(N) is equal to the combined sum of the terms in the first and second groups, G(N) = g(N, M ) = XM g(N -1,M ) + g(N,M -1) This same recurrence relation clearly exists for any intermediate value of n from 1 to N, and for any intermediate value of m from 1 to M .
Buzen’s algorithm is simply the iterative application of this fundamental recurrence relation, along with the following boundary conditions.
g(0,m) = 1 for m = 1, 2, …M g(n,1) = (Xi)n for n = 0, 1, … N The Gordon-Newell theorem enables analysts to determine the stationary probability associated with each individual state of a closed queueing network.
Many of these marginal probabilities can be computed with minimal additional effort.
Clearly, Xi must be raised to the power of k or higher in every state where the number of customers at service center i is greater than or equal to k. Thus Xi k can be factored out from each of these probabilities, leaving a set of modified probabilities whose sum is given by G(N-k)/G(N).
This observation yields the following simple and highly efficient result: P(ni ≥ k) = (Xi)k G(N-k)/G(N) This relationship can then be used to compute the marginal distributions and expected number of customers at each service facility.
The expected number of customers at service facility i is given by
These characterizations of quantities of interest in terms of the G(n) are also due to Buzen.
[2] It will be assumed that the Xm have been computed by solving the relevant equations and are available as an input to our routine.
Note that C[0] remains equal to 1 throughout all subsequent iterations.
In the second loop, each successive value of C(n) for n≥1 is set equal to the corresponding value of g(n,m) as the algorithm proceeds down column m. This is achieved by setting each successive value of C(n) equal to: g(n,m-1) plus Xm times g(n-1,m).