Buzen's algorithm: Difference between revisions

Revision as of 23:59, 18 April 2009

Latest revision as of 11:22, 27 May 2025

(47 intermediate revisions by 20 users not shown)

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=Buzen, J.P. |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 | access-date = 2006-04-15 | archive-date = 2016-05-13 | archive-url = https://web.archive.org/web/20160513192804/http://www-unix.ecs.umass.edu/~krishna/ece673/buzen.pdf | url-status = dead }}</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>

{{Refimprove|date=July 2007}}

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 <math>\tbinom{N+M-1}{M-1}</math> 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.<ref name=":1">{{cite journal |last1=Denning |first1=Peter J. |date=24 August 2016 |title=Rethinking Randomness: An interview with Jeff Buzen, Part I |url=https://dl.acm.org/doi/10.1145/2986329 |journal=Ubiquity |volume=2016 |issue=August |pages=1:1–1:17 |doi=10.1145/2986329|doi-access=free }}</ref> 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.

In [[queueing theory]], a discipline within the mathematical [[probability theory|theory of probability]], '''Buzen's algorithm''' is an algorithm for calculating the [[normalization constant]] ''G''(''K'') in the [[Gordon–Newell theorem]]. This method was first proposed by [[Jeffrey P. Buzen]] in 1973.<ref name="buzen-1973">{{cite journal

| first = Jeffrey

| last = Buzen

| authorlink = Jeffrey P. Buzen

| year = 1973

| month = September

| title = Computational algorithms for closed queueing networks with exponential servers

| journal = Communications of the ACM

| volume = 16

| issue = 9

| doi = 10.1145/362342.362345

| pages = 527

}} [http://www-unix.ecs.umass.edu/~krishna/ece673/buzen.pdf]</ref> Once ''G'' is computed the probability distributions for the network can be found. In contrast, [[mean value analysis]] is an alternative algorithm that can also be used to derive some performance measures (such as the mean queue lengths) without having to directly compute the normalization constant.

==Problem setup==

The motivation for this algorithm is efficiency: a straightforward Gordon-Newell calculation would require the enumeration of all states that the system can be in, resulting in a combinatorial explosion. Buzen's algorithm is of order ''N''², making the application of the G-N theorem practical, and opening up a large class of queueing systems to accurate modeling.

Consider a closed queueing network with ''M'' service facilities and ''N'' circulating customers. Assume that the service time for a customer at service facility ''i'' is given by an [[exponentially distributed]] random variable with parameter ''μ''_''i'' and that, after completing service at service facility ''i'', a customer will proceed next to service facility ''j'' with probability ''p''_''ij''.<ref name="gn" />

In the queuing literature Buzen's algorithm is sometimes also referred to as the '''convolution algorithm'''.

Let <math>\mathbb P(n_1,n_2,\cdots,n_M) </math> be the steady state probability that the number of customers at service facility ''i'' is equal to ''n_i'' for ''i'' = 1, 2, ... , ''M .'' It follows from the [[Gordon–Newell theorem]] that

==Problem setup==

<math>\mathbb P(n_1,n_2,\cdots,n_M) = \frac{1}{\text{G}(N)}</math><math> \left( X_1 \right)^{n_1}</math><math> \left( X_2 \right)^{n_2}</math> .... <math> \left( X_M \right)^{n_M}</math>

Consider a closed queueing network with ''M'' service facilities and ''N'' circulating customers. Let <math>n_i</math> represent the number of customers present at the ''i''th facility, with <math>\sum_{i=1}^M n_i = N</math>. We assume that the service time for a customer at the ''i''th facility is given by an exponentially distributed random variable with mean <math>1/\mu_i</math> and that after completing service at the ''i''th facility a customer will proceed to the ''j''th facility with probability ''p''_''ij''.

This result is usually written more compactly as

It follows from the [[Gordon–Newell theorem]] that the equilibrium distribution of customers in the network is given by:

:<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>

whereThe thevalues of ''X''<mathsub>X_i''i''</mathsub> are found fromdetermined theby equations:solving

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

and ''G''(''N'') is a normalizing constant suchchosen so that the abovesum probabilitiesof sumall <math>\tbinom{N+M-1}{M-1}</math> values of <math>\mathbb P(n_1,n_2,\cdots,n_M) </math> is equal to 1. <ref> Buzen:'s op.citalgorithm represents the first efficient procedure for computing G(''N'').<ref name="buzen-1973" /><ref name=":1" />

==Algorithm description==

The purpose of the Buzen algorithm is to numerically compute values of ''G''.

The individual terms that must be added together to compute G(''N'') all have the following form:

==Marginal distributions, expected number of customers==

<math> \left( X_1 \right)^{n_1}</math><math> \left( X_2 \right)^{n_2}</math> .... <math> \left( X_M \right)^{n_M}</math>. Note that this set of terms can be partitioned into two groups. The first group comprises all terms for which the exponent of <math> \left( X_M \right)</math> is greater than or equal to 1.  This implies that <math> \left( X_M \right)</math> raised to the power 1 can be factored out of each of these terms.  

The general G-N distribution given above is mainly of theoretical interest. However, expressions for a number of useful performance measures can be derived from it. Buzen showed that the probability that there are exactly ''k'' customers at facility ''i'' is given by:

After factoring out <math> \left( X_M \right)</math>, 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 “''X''_''M'' times G(''N'' -1)”. This insight provides the foundation for the development of the algorithm.<ref name=":1" />  

:<math>P(n_i = k) = \frac{X_i^k}{G(N)}[G(N-k) - X_i G(N-k-1)]</math>

Next consider the second group.  The exponent of ''X''_''M'' for every term in this group is zero.  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.

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

To express this concept precisely, 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, ''m'' 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'' .

:<math>E[n_i] = \sum_{k=1}^N X_i^k \frac{G(N-k)}{G(N)}.</math>

Given this definition, the sum of the terms in the second group can now be written as g(''N'', ''M'' -1).

Note that these expressions also involve G. It is in the evaluation of these expressions that the algorithm is useful.

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'' )”.  

==Derivation==

The algorithm computes not just ''G'' but a two-dimensional function

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

Once the calculation ends the values we are interested in are found by

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

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

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

The definition of ''g'' and a recursive method of calculating it are as follows

G(''N'') = g(''N'', ''M'' ) = ''X_M'' g(''N'' -1,''M'' ) + g(''N'',''M'' -1)

:<math>

\begin{align}

g(n, m) & = \sum_{n_1 + \cdots + n_m = n} \prod_{i=1}^m (X_i)^{n_i} \\

& = \sum_{(n_1 + \cdots + n_m = n) \wedge (n_m=0)}^n \prod_{i=1}^m (X_i)^{n_i} +

\sum_{(n_1 + \cdots + n_m = n) \wedge (n_m>0)}^n \prod_{i=1}^m (X_i)^{n_i} \\

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

\end{align}

</math>

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'' .  

Also

This implies g(''n,m'') = ''X_m'' g(''n'' -1,''m'') + g(''n,m'' -1).  Buzen’s algorithm is simply the iterative application of this fundamental recurrence relation, along with the following boundary conditions.

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

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

and

g(''n'',1)  =  (''X''_i)^''n'' for ''n'' = 0, 1, … ''N''

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

==Marginal distributions, expected number of customers==

This recursive relationship allows for the calculation of all ''G''(''N'') up to any value of ''N'' in [[Big O notation|order]] O(''MN'') time.

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 other important probabilities. For example P(''n_i'' ≥ ''k''), the probability that the total number of customers at service center ''i'' is greater than or equal to ''k'', must be summed over all values of ''n_i'' ≥ ''k'' and, for each such value of ''n_i'', over all possible ways the remaining ''N'' – ''n_i'' customers can be distributed across the other ''M'' -1 service centers in the network.

Many of these marginal probabilities can be computed with minimal additional effort.  This is easy to see for the case of P(''n_i'' ≥ 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 and highly efficient result:

P(''n_i'' ≥ ''k'') = (''X_i'')^''k'' G(''N''-''k'')/G(''N'')

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)}.</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>

These characterizations of quantities of interest in terms of the G(''n'') are also due to Buzen.<ref name="buzen-1973"/>

==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(''gn,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''.

The first loop in the algorithm below initializes the column vector C[n] so that C[0] = 1 and C(n) = 0 for n≥1.   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 ''X_m'' times g(''n-1,m'').  

Note that g(''n,m-1'') is the previous value of C(n), and g(''n-1,m'') is the current value of C(n-1)

<sourcesyntaxhighlight lang="pascal">

C[0] := 1

for n := 1 step 1 until N do

Line 84 ⟶ 88:

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];

</syntaxhighlight>

</source>

At completion, the final values of C[n] correspond to column ''CM'' containsin the matrix g(''n,m'').  Thus they represent the desired values G''G(0),'' G''(1), ... ,'' to G''G(N)''. <ref>Buzen, op.name="buzen-1973" cit.</ref>

==References==

Line 94 ⟶ 98:

{{reflist}}

*[httphttps://www.cs.wustl.edu/~jain/cse567-08/ftp/k_35ca.pdf Jain: The Convolution Algorithm (class handout)]

*[https://cs.gmu.edu/~menasce/cs672/slides/CS672-convolution.pdf Menasce: Convolution Approach to Queueing Algorithms (presentation)]

*[http://cs.gmu.edu/~menasce/cs672/slides/CS672-convolution.pdf Menasce: Convolution Approach to Queueing Algorithms (presentation)]

{{Queueing theory}}

{{probability-stub}}

[[Category:Stochastic processes]]

[[Category:Queueing theory]]

[[Category:Statistical algorithms]]