Buzen's algorithm: Difference between revisions

Revision as of 19:58, 22 May 2023

Revision as of 00:14, 8 June 2023

Line 1:

In [[queueing theory]], a discipline within the mathematical [[probability theory|theory of probability]], '''Buzen's algorithm''' (or '''convolution algorithm''') is an algorithm for calculating the [[normalization constant]] G(''N'') in the [[Gordon–Newell theorem]]. This method was first proposed by [[Jeffrey P. Buzen]] in his 1971 PhD dissertation<ref name=":0">{{Cite book |last=Defense Technical Information Center |url=http://archive.org/details/DTIC_AD0731575 |title=DTIC AD0731575: Queueing Network Models of Multiprogramming |date=1971-08-01 |language=english}}</ref> and subsequently published in a refereed journal in 1973.<ref name="buzen-1973">{{Cite journal | last1 = Buzen | first1 = J. P. | author-link = Jeffrey P. Buzen| title = Computational algorithms for closed queueing networks with exponential servers | doi = 10.1145/362342.362345 | url = http://www-unix.ecs.umass.edu/~krishna/ece673/buzen.pdf| journal = Communications of the ACM | volume = 16 | issue = 9 | pages = 527–531 | year = 1973 | s2cid = 10702 }}</ref> Computing G(''N'') is required to compute the stationary [[probability distribution]] of a closed queueing network.<ref name="gn">{{Cite journal | last1 = Gordon | first1 = W. J. | last2 = Newell | first2 = G. F. | author-link2 = Gordon F. Newell| doi = 10.1287/opre.15.2.254 | jstor = 168557| title = Closed Queuing Systems with Exponential Servers | journal = [[Operations Research (journal)|Operations Research]]| volume = 15 | issue = 2 | pages = 254 | year = 1967 }}</ref>

Performing a naïve computation of the normalisingnormalizing constant requires enumeration of all states. For a closed network with ''N'' circulating jobs and ''M'' service centers, G(''N'') is the sum of <math>\tbinom{N+M-1}{M-1}</math> individual terms, with each term comprised of ''M'' factors raised to powers whose total sum is ''N''. Buzen's algorithm computes G(''N'') using a total of ''NM'' multiplications and ''NM'' additions. This is a significant improvement andthat allowsopened forthe computationsdoor to beapplying performedthe withGordon-Newell muchtheorem to models of largerrealistic networkssize.<ref name="buzen-1973" /> TheIn addition, the values of G(1), G(2) ... G(''N''-1), which can be used to express other important quantities of interest,<ref name=":0" /> are also computed as by-products of the algorithm.

==Problem setup==

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 ''i''th facility at time ''t'', such that <math>\scriptstyle \sum_{i=1}^M n_i = N</math>. We assumeAssume that the service time for a customer at the ''i''th facility is given by an [[exponentially distributed]] random variable with parameter ''μ''_''i'' and that, after completing service at the ''i''th facility, a customer will proceed next to the ''j''th facility with probability ''p''_''ij''.<ref name="gn" />

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

::<math>\mathbb P(n_1,n_2,\cdots,n_M) = \frac{1}{\text{G}(N)}\prod_{i=1}^M \left( X_i \right)^{n_i}</math>

where the ''X''_''i'' are found by solving

::<math>\mu_j X_j = \sum_{i=1}^M \mu_i X_i p_{ij}\quad\text{ for }j=1,\ldots,M.</math>

where

and ''G''(''N'') is a normalizing constant chosen that the above probabilities sum to 1.<ref name="buzen-1973" />

<math>\mathbb P(n_1,n_2,\cdots,n_M) </math> is the steady state probability that the number of customers at the ''i''th service facility is equal to ''n_i'' for i = 1, 2, ... , ''M''

::and where the ''X''_''i'' are found by solving<math>\mu_j X_j = \sum_{i=1}^M \mu_i X_i p_{ij}\quad\text{ for }j=1,\ldots,M.</math>

and ''G''(''N'') is a normalizing constant chosen so that the sum of all the above probabilities sumis equal to 1.<ref name="buzen-1973" />

Buzen's algorithm is an efficient method to compute G(''N'').<ref name="buzen-1973" />

Line 17 ⟶ 23:

==Algorithm description==

The algorithm is based on partitioning the set of terms that add up to G(''N'') into two groups.  The first group is comprised of all terms for which the exponent of ''X''_''M'' is greater than or equal to 1.  This implies that ''X''_''M'' can be factored out of each of these terms.  

Write g(''N'',''M'') for the normalizing 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<ref name="buzen-1973" />

::<math>g(0, m) = 1 \text{ for }m=1,2,\cdots,M</math>

After factoring out ''X''_''M'' , a surprising result emerges: the sum of the modified terms in the first group are exactly equal to 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 “''X''_''M'' times G(''N'' -1)”.

::<math>g(n, 1) = (X_1)^n \text{ for }n=0,1,\cdots,N.</math>

The recurrence relation<ref name="buzen-1973" />

Now consider the second group.  The exponent of ''X''_''M'' for every term in this group is zero.  In effect, the ''M'' ^th service facility 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''.

::<math>g(n, m) = g(n,m-1)+X_m g(n-1,m).</math>

is used to compute a grid of values. The sought for value G(''N'')&nbsp;=&nbsp;g(''N'',''M'').<ref name="buzen-1973" />

To express this result mathematically, assume that X₁, X₂, … ''X''_''M'' 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, service facilities 1,2, … m, and values of  X₁, X₂, … ''X''_''m''  that match the first m members of the original sequence X₁, X₂, … ''X''_''M''  

Given this definition, the normalizing constant G(''N'') in the Gordon-Newell theorem can now be re-written as g(''N'',''M'').

It also follows immediately that “X''_M'' times G(''N'' -1)”, the sum of the terms in the first group, can be re-written as “X''_M'' times g(''N'' -1,''M'' )”.   More importantly, the sum of the terms in the second group can now be written as g(''N'', ''M'' -1).

Since the combined sum of the terms in the first and second groups is equal to G(''N''),

G(''N'') = g(''N'', ''M'' ) = X_M 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'' .  This implies

g(n,m) = X_m g(n-1,m) + g(n,m-1).  

The thought process that led to the discovery of this recurrence relation is discussed in the final sections of a 2016 interview. Buzen’s algorithm is simply the iterative application of this this fundamental recurrence relation, along with the following boundary conditions.

g(0,m) = 1 for m = 1, 2, …''M''

g(n,1)  =  (X_i)ⁿ for n = 0, 1, … ''N''

==Marginal distributions, expected number of customers==

The Gordon-Newell theorem enables analysts to determine the stationary probability associated with each individual state of a closed queueing network.  These individual probabilities must then be added together to evaluate marginal probabilities such as P(n_j≥k), the probability that the total number of customers at service center i is greater than or equal to k (summed over all values of n_j≥k and all possible ways the remaining N – n_i customers can be distributed across the other M-1 service centers in the network).

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

::<math>\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,</math>

Many of these marginal probabilities can be computed with minimal additional effort.  This is easy to see for the case of P(n_j≥k).   Clearly, X_i 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 (X_i)^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 result:

::<math>\mathbb P(n_i = N) = \frac{X_i^N}{G(N)}[G(0)].</math>

the expected number of customers at facility ''i'' by

::P(n<mathsub>\mathbb E(n_ij</sub>≥k) = \sum_{k=1}^N X_i^(X_i)^k \frac{G(N-k)}{/G(N)}.</math>

This relationship can then be used to compute the [[marginal distribution]]s and [[expected value|expected]] number of customers at each service facility.<math>\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,</math><math>\mathbb P(n_i = N) = \frac{X_i^N}{G(N)}[G(0)].</math>

The expected number of customers at service facility ''i'' is given by<math>\mathbb E(n_i) = \sum_{k=1}^N X_i^k \frac{G(N-k)}{G(N)}.</math>

==Implementation==

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''(n,m) is in principle a two dimensional matrix, it can be computed in a column by column fashion starting from the top of the leftmost column and running down each column to the bottom before proceeding to the next column on the right. The routine uses a single column vector ''C'' to represent the current column of ''g''.

<syntaxhighlight lang="pascal">

Line 45 ⟶ 74:

</syntaxhighlight>

At completion, ''C'' contains the desired values ''G(0), G(1), ... ,'' to ''G(N)''. <ref name="buzen-1973" />

==References==